The SILCC project - V. The impact of magnetic fields on the chemistry and the formation of molecular clouds
Magnetic fields are ubiquitously observed in the interstellar medium (ISM) of present-day star-forming galaxies with dynamically relevant energy densities. Using three-dimensional magneto-hydrodynamic (MHD) simulations of the supernova (SN) driven ISM in the flux-freezing approximation (ideal MHD) we investigate the impact of the magnetic field on the chemical and dynamical evolution of the gas, fragmentation and the formation of molecular clouds. We follow the chemistry with a network of six species (H, H, H, C, CO, free electrons) including local shielding effects. We find that magnetic fields thicken the disc by a factor of a few to a scale height of , delay the formation of dense (and molecular) gas by and result in differently shaped gas structures. The magnetised gas fragments into fewer clumps, which are initially at subcritical mass-to-flux ratios, , and accrete gas preferentially parallel to the magnetic field lines until supercritial mass-to-flux ratios of up to order 10 are reached. The accretion rates onto molecular clouds scale with . The median of the inter-cloud velocity dispersion is and lower than the internal velocity dispersion in the clouds (). However, individual cloud-cloud collisions occur at speeds of a few .
keywords:ISM: clouds – ISM: evolution – ISM: kinematics and dynamics – ISM: magnetic fields – ISM: molecules – (magnetohydrodynamics) MHD
Magnetic fields are ubiquitously observed in the interstellar medium (ISM, Beck 2009; Crutcher 2012; Haverkorn 2015) and might play a vital role for the evolution of galaxies (e.g. Naab & Ostriker, 2017). With mean field strengths in the spiral structures of galaxies of a few the magnetic energy density is typically lower than the kinetic but overall comparable to the mean thermal energy density (e.g. Boulares & Cox, 1990; Cox, 2005). The high conductivity and the partial ionisation in the ISM result in an efficient coupling of the magnetic field with the gas and the assumption of ideal MHD is valid for a large fraction of the ISM. This allows for a dynamical impact of the magnetic field on the motions of the gas from scales of the entire galaxy (Beck, 2001; Hanasz et al., 2009; Kotarba et al., 2009; Pakmor & Springel, 2013; Rieder & Teyssier, 2017), representative parts of the ISM (de Avillez & Breitschwerdt, 2005; Koyama & Ostriker, 2009a, b; Hill et al., 2012; Kim & Ostriker, 2015b; Walch et al., 2015; Girichidis et al., 2016; Pardi et al., 2017), molecular clouds (Shu et al., 1987; Padoan & Nordlund, 1999; Heitsch & Burkert, 2002; Mac Low & Klessen, 2004; Banerjee et al., 2009) down to star-forming regions (Hennebelle & Fromang, 2008; Hennebelle & Teyssier, 2008; Hennebelle & Ciardi, 2009; Peters et al., 2011; Seifried et al., 2011b, a, 2012; Klassen et al., 2017; Körtgen & Banerjee, 2015).
Turbulent or turbulence-like motions, that are trans- or supersonic (Goldreich & Kwan, 1974; Zuckerman & Evans, 1974; Elmegreen & Scalo, 2004; Scalo & Elmegreen, 2004) stir the ISM. Supernovae (SNe, e.g. Gatto et al. 2015; Walch & Naab 2015), stellar feedback processes like radiation (Walch et al., 2013; Kim et al., 2016; Peters et al., 2017; Haid et al., 2018) and winds (Gatto et al., 2017), gravitational instability, and differential rotation in the galactic disc form a multitude of complex dynamical interactions that are expected to shape the ISM on different scales and with different efficiencies (Klessen & Glover, 2016). SNe and gravitational contraction are regarded as the most energetic dynamical drivers.
Theoretical considerations describe the interstellar medium in present-day star-forming disc galaxies as a multi-component fluid consisting of two (cool and warm) stable (Field, 1965; Field et al., 1969; Wolfire et al., 1995, 2003) and one hot meta-stable phases (Cox & Smith, 1974; McKee & Ostriker, 1977; Ferrière, 2001). However, due to the strong dynamics a significant fraction of the gas is in the unstable regime. In addition, there is a constant exchange between the phases. The hot gas with temperatures of the order of originates from supernovae and fills most of the volume. The warm and cold component coexist in approximate pressure equilibrium (e.g. Cox, 2005). The hot medium is expected to be dominated by thermal and kinetic rather than magnetic energy. Contrary, the cold regions can be permeated by strong magnetic fields that can retard or theoretically completely prevent gravitational contraction and the formation of dense clouds and stars (Crutcher et al., 2010; Crutcher, 2012; Pardi et al., 2017). The transition from magnetically sub- to supercritical regions, which is needed for gravitational forces to overcome pressure balance is still a matter of active discussion (e.g. Körtgen & Banerjee, 2015; Valdivia et al., 2016). In the case of strong fields the contraction perpendicular to the field lines is prohibited due to the magnetic pressure. However, condensations can form along the field lines (Field, 1965; Heitsch & Hartmann, 2014; Walch et al., 2015; Girichidis et al., 2016; Iffrig & Hennebelle, 2017).
This paper is part of the SILCC project111https://hera.ph1.uni-koeln.de/~silcc/ (SImulating the Life Cycle of molecular Clouds) that investigates the processes in the interstellar medium using hydrodynamical simulations. The first two studies (Walch et al., 2015; Girichidis et al., 2016) focus on different modes of SN driving and varying SN rates. Both find that SNe partially need to explode in low-density environments to create a realistic chemical composition and launch outflows from the disc. Gatto et al. (2017) uses active star clusters that self-consistently form in dense regions and accrete gas. The study also includes stellar winds in addition to SN feedback. Peters et al. (2017) further added radiation from star clusters using the same formalism. All of these studies have a maximum grid resolution of , which is too coarse to study molecular clouds. The SILCC-zoom study (Seifried et al., 2017) specifically investigates the condensation of gas and the formation of clouds with a maximum resolution of in the non-magnetised ISM and without the formation of star or star cluster particles. The high resolution accounts for a more detailed chemical evolution but restricts the study to only a few clouds. With this study we specifically aim at probing the impact of magnetic fields on the chemical evolution along with the formation of molecular clouds while keeping the dynamical connection to the large-scale ISM flows (see dicussion in Seifried et al., 2017). We increase the resolution to , which allows us to more accurately follow the evolution in dense regions and at the same time resolve many clouds with higher resolution to have a statistically relevant sample of molecular clouds.
Following the numerics and the simulation setup (Sections 2 and 3) we first give an overview of the morphological evolution of the simulations (Section 4), followed by an analysis of the degree of magnetisation of the gas (Section 5). We then focus on the formation of shielded and molecular gas (Section 6) and the fragmentation into clouds and their properties (Section 7). Discussions and conclusions are presented in Sections 8 and 9.
2 Numerical Methods
A detailed description of the numerical setup, the equations to solve as well as the numerical methods used is given in Walch et al. (2015) and Girichidis et al. (2016). Here we only present a brief summary. We use the hydrodynamical code Flash in version 4 (Fryxell et al. 2000; Dubey et al. 2008, http://flash.uchicago.edu/site/) which is parallelised using the Message Passing Interface (MPI). The Eulerian grid divides the computational domain in blocks of cells which can be adaptively refined (Adaptive Mesh Refinement, AMR).
The magneto-hydrodynamic (MHD) equations are solved using the HLLR5 finite-volume scheme for ideal MHD (Bouchut et al., 2007, 2010; Waagan, 2009; Waagan et al., 2011). The directionally split solver ensures positivity of the density and pressure by construction and is suitable for high-Mach-number flows.
Gravitational effects are included in two ways. We solve the Poisson equation for the effects of self-gravity using the tree-based method by Wünsch et al. (2018). In addition we include an external potential that accounts for the acceleration due to the stellar component of the disc. We use tabulated values based on an isothermal sheet (Spitzer, 1942) with a stellar surface density of and a scale height of . The acclerations are comparable to the Milky-Way values by Kuijken & Gilmore (1989) for low heights above the midplane.
In order to both follow the chemical state of the gas and compute radiative cooling accurately, we use a chemical network that includes ionised (H), atomic (H) and molecular hydrogen (H) as well as singly ionized carbon (C) and carbon monoxide (CO). This allows us to include non-equilibrium abundances as in Glover & Mac Low (2007) and Micic et al. (2012) using the carbon chemistry by Nelson & Langer (1997). We assume a temporally constant interstellar radiation field of (Habing, 1968; Draine, 1978), which is locally attenuated based on how strongly shielded the computational cell is. The column densities and the derived optical depth in every cell are computed using the Flash implementation (Wünsch et al., 2018) of the TreeCol algorithm (Clark et al., 2012). For every cell we integrate the column density out to a radial distance of . The column density dependent attenuation factor follows Glover & Clark (2012).
For radiative cooling we follow the atomic and molecular cooling functions of Glover et al. (2010) and Glover & Clark (2012). High-temperature cooling above is based on the rates described in Gnat & Ferland (2012). Heating includes a constant cosmic ray ionisation and corresponding heating rate (Goldsmith & Langer, 1978), X-ray heating (Wolfire et al., 1995), and photoelectric heating that is coupled to the optical depth (Bakes & Tielens, 1994; Bergin et al., 2004; Wolfire et al., 2003). We assume a constant dust-to-gas mass ratio of 0.01 using dust opacities of Mathis et al. (1983) and Ossenkopf & Henning (1994).
3 Simulation setup and parameters
Shown are the name of the simulation, the initial magnetic field strength in the midplane, the effective numerical resolution, the corresponding spatial resolution of the smallest cell as well as the threshold values for refinement and derefinement.
We set up a stratified box with a size of . The boundary conditions in and are periodic. For the boundary in we use diode boundary conditions, which allow gas to leave but not enter the box. We emphasise that recent numerical studies of the SN-driven ISM with different mean densities and SN positioning in periodic boxes discuss in detail the problem of an over-pressurised ISM (thermal runaway, Gatto et al. 2015; Li et al. 2015), which is aggravated by a closed box. The gas density follows a Gaussian distribution in with a scale height of . The temperature is set such that the gas is in pressure equilibrium. In the densest regions of the disc at the composition of the gas is purely atomic with a minimum temperature of . At the lowest density () at large we set and assume that all the gas is ionised. The gas surface density is , which yields a total mass of . Initially the gas is at rest.
We distinguish between a non-magnetic setup and two magnetic setups with different initial field strengths. The magnetic field is oriented along the direction an with initial field strengths in the midplane of and . The magnitude of scales with the square root of the density, . We do not introduce a random component of the magnetic field.
Feedback from stars is included as individual SNe exploding at a constant rate. We refrain from using sink particles and an accretion based star formation and SN rate (e.g. Schmidt, 1959) because of the uncertainties of how to treat the magnetic field at the interface between sink particles and the gas. By using the same SN rate and positioning, it is also easier to directly compare the different magnetisations. Based on the total gas surface density in the box we use the Kennicutt-Schmidt (Kennicutt, 1998) relation to derive a star formation rate, which we convert into a SN rate using the stellar initial mass function of Chabrier (2003). For our surface density of this yields SNe per Myr. We apply the clustered SN driving from Walch et al. (2015) and Girichidis et al. (2016), which we summarize briefly. The total SN rate is split into a fraction of 20% type Ia SNe, which explode with random and positions and a Gaussian distribution in with a scale height of . The remaining 80% are SNe of type II with a scale height of . We split the latter component into of individual SNe (representing run-away OB stars) and of SNe, which are associated with clusters (Cowie et al., 1979; Gies & Bolton, 1986; Stone, 1991; Hoogerwerf et al., 2001; de Wit et al., 2005; Tetzlaff et al., 2011) and are expected to be an important agent for driving a hot phase and superbubbles (e.g. Mac Low & McCray, 1988; Gatto et al., 2015; Li et al., 2015). The clusters have an associated life time of and contain SNe, randomly drawn from a power law distribution , where is the number of SNe. All SNe of a cluster explode at the same position. We ensure a constant SN rate in the box, which requires the temporal spacing between SNe in individual clusters to vary. The random positions of all SNe are precomputed beforehand and stored in a table, so that all simulations have exactly the same SN positions. If we resolve the Sedov-Taylor radius with at lest four cells at the highest refinement level (), we inject of thermal energy per SN into a dynamically computed volume with a total mass of . If we do not meet this resolution criterion, we heat the gas in the minimum volume () to and inject momentum according to Blondin et al. (1998), see also Gatto et al. (2015). The resolution in all simulations is high enough to resolve the Sedov-Taylor radius with at least 4 cells for more than of all SNe, which ensures that the efficiency of the explosions is not impaired by numerical over-cooling or errors in the momentum injection, e.g. for overlapping SN remnants with cancelling momenta (see also the resolution requirements discussed by Kim & Ostriker, 2015a).
The base resolution is cells corresponding to a cell size of . We allow for one and two additional levels of refinement with cell sizes of (name tag 2pc) and (1pc), respectively. The local adaptive refinement is based on two quantities: The density field following Lohner (1987), which is the standard in Flash, as well as on the fraction of molecular hydrogen. We refine a block if the molecular mass fraction of a cell exceeds and derefine if the values fall below . For comparison we perform two additional runs with a ten times higher refinement/derefinement threshold, . We perform eight different simulations listed in Table 1. The first column shows the name, which encodes the initial field strength (B0, B3, B6), the resolution (2pc and 1pc) as well as the additional higher refinement threshold (ht), where needed. All simulations use the same SN rate as well as the same clustering and SN positioning. We run all simulations for .
4 Morphological evolution and global dynamics
The morphological evolution differs early on when comparing the simulations of varying initial magnetic field strength. Changes of the maximum resolution and refinement threshold result in significantly smaller effects than the differences in magnetisation. We therefore concentrate here on the differences in and only discuss resolution and refinement details in the appendix. We show face-on column density plots in Fig. 1 for simulations B0-1pc (left), B3-1pc (centre) and B6-1pc (right) at simulation times (, , , and from top to bottom) and the corresponding edge-on view in Fig. 2. In the non-magnetic runs, the SNe can quickly form over-dense regions, which cool and condense to filamentary structures and clouds near the midplane under the effect of the external potential and self-gravity. After the simulation reaches column density contrasts of more than six orders of magnitude ranging from in voids to in numerous small dense clouds forming at the intersections of elongated filaments. After dense individual clouds have formed in the non-magnetic run with a sharp transition from the dense cloud to the low-density voids. Over time the small individual clouds in the non-magnetic run merge to eventually result in approximately 10 clouds after with masses of each. In the first the two magnetised runs mostly form diffuse overdensities. The maximum column densities of are found in only a few clouds in B3-1pc and are barely reached in B6-1pc. In addition, the magnetised clouds are embedded in larger coherent structures with column densities of . The field provides direct magnetic support in addition to thermal and turbulent pressure, which prevents or delays the condensation of gas into thin filaments and sheets. In addition, the magnetic pressure attenuates the impact of the SNe leading to lower maximum densities in the swept-up gas, which in turn results in less efficient cooling and more thermal pressure support. The delayed formation of dense gas has also been reported by Heitsch et al. (2001) and more recently by Hennebelle & Iffrig (2014) and Iffrig & Hennebelle (2017) as a later accumulation of gas in sink particles and a lower star formation rate. The smaller number of clouds and the more diffuse envelope in the magnetic runs is consistent with similar simulations by Pardi et al. (2017) in a periodic box. After the magnetised simulations end up with a similar number and mass distribution of clouds as the non-magnetic counterpart due to continuous merging of clumps.
A common behaviour in all simulations is that early overdensities created by the first SNe continuously grow in mass and size. The SN feedback influences the formation process but does not provide an efficient destruction mechanism for the clouds. This means that we do not simulate the full cycle of a molecular cloud but focus on the formation process. To some extent the destruction can be enhanced by placing SNe always in the densest regions, which mimics the destruction of clouds from inside out. However, previous models like the peak and mixed driving runs in Walch et al. (2015) and Girichidis et al. (2016) lead to an unrealistic composition of the ISM (peak driving) and in general do not solve the problem of the indestructible clouds, see also Section 8.
4.2 Vertical structure
Fig. 2 shows column densities edge-on and reveals the different evolution in the vertical structure. Without magnetic pressure support the overdensities in the intersections of swept-up SN shells can cool, lose their thermal pressure support and contract towards the midplane. After only the gas has formed a thin sheet with few vertical fluctuations. Including magnetic fields, the gas can delay the thermal compression by SNe as well as the gravitational attraction until the first dense clouds form at for B3-1pc and at for B6-1pc. The gas thus remains in a diffuse state for a longer time and provides a larger effective cross section for the SN blast waves, which results in a thicker disc with elongated structures of diffuse gas along the vertical direction. The dense clouds form later and partially in the extended elongated structures. As a consequence the dense gas has a larger spread in the vertical distribution. Once compact clouds have formed out of the diffuse gas and the effective cross section is reduced, the gravitational attraction dominates over the outward pointing thermal and magnetic support of the environment and the clouds move towards the midplane. This causes the vertical heights to decrease towards the end of the simulation.
Before quantifying the results of the simulations in terms of different density and chemical composition, we would like to illustrate the relation between them. Fig. 3 presents the volume weighted density distribution for simulation B3-1pc at . The lines show the total gas as well as the individual hydrogen species, H, H, and H. All chemical distributions widely overlap over several orders of magnitude (see also discussion in Walch et al., 2015). In the following discussion we therefore investigate various quantities with weighting by mass, volume and chemical species separately.
Fig. 4 depicts the chemical mass fractions, () over time (top) as well as the vertical extent of the gas (bottom). Initially, all the gas is in atomic form, which starts condensing to form molecular hydrogen at different times depending on magnetisation, see Section 6 for details. The sum of atomic and molecular gas makes up for more than 98 per cent of the total hydrogen mass. The lower panel shows the heights of 90% enclosed mass for all the gas (solid lines) and 90% of the molecular gas (dotted lines). The visual impression of the column density plots is well reflected in the time evolution plots. The vertical scales for the gas reach a few s of parsecs in the non-magnetic runs. The magnetic field helps to lift the gas to heights of for the weak magnetic field and up to for the strong magnetic field. In the case of H the vertical heights show the same trend but at a lower absolute height with , , and up to for B0-1pc, B3-1pc, and B6-1pc, respectively. We note that the scale height of the gas is always well below . Observations of the vertical mass distribution in the Milky Way based on CO (corresponding more to our molecular gas rather than the total gas) reveal an exponential scale height of (Langer et al., 2014), which would be more consistent with our magnetic runs. In other galaxies the scale height can range from for the molecular gas (Caldú-Primo et al., 2013; Yim et al., 2014). This range is large and besides different physical effects from within the disc, global disc dynamics like warps or accretion onto the galaxy are likely to play a role for the vertical distribution of gas in individual systems. Ignoring the global effects, the observations seem to be more consistent with our magnetic runs. In the Milky Way giant molecular clouds (GMCs) tend to populate regions much closer to the midplane compared to low-mass clouds (Stark & Lee, 2005). This agrees well with our evolution of the cloud population in the magnetic runs. The smaller clouds form at larger altitudes and merge to GMCs close to the midplane, which causes the scale height to decrease along with the formation of massive clouds via accretion and merging.
The applied SN driving is expected to result in fountain flows that exceed the boundary at (Hill et al., 2012; Walch et al., 2015; Girichidis et al., 2016). We notice outflowing gas in all simulations, which reduces the total amount of gas in the box. Fig. 5 shows the total mass in the simulation box over time for the main simulations in this study (B0-1pc, B3-1pc, and B6-1pc; solid lines) as well as simulation B3-2pc (long-dashed line). In addition we plot the total mass of the same volume (up to the height of above the midplane) in the clustered non-magnetic simulation S10-KS-clus and the clustered magnetic run S10-KS-clus-mag3 of Girichidis et al. (2016, G16) with a resolution of . In the non-magnetic run the mass loss is over the total simulation time in very good agreement with the previous simulation at a four times coarser resolution. The larger vertical extent of the gas in simulations B3-1pc and B6-1pc result in losses of and , respectively. The comparison between B3-1pc, B3-2pc and S10-KS-clus-mag3 (all red lines) indicate that the mass loss through the boundary at is resolution dependent. With higher resolution we can resolve the Sedov-Taylor phase of more SNe, which allows for a more efficient stirring of the gas. In addition, higher resolution allows for a more efficient small-scale dynamo and a resulting stronger magnetic support, which helps lifting the gas above . We note that with the outflows from the box, we also lose magnetic flux through the boundary. However, we are primarily losing low-density gas through the boundary, which is weakly magnetised. As a result the accompanying loss of magnetic flux with the outflow is unlikely to affect the dynamics in the region around the midplane.
4.3 Global dynamics
Before analysing the global dynamics in the box it is worth noting that the gas is initially at rest. The missing shear and the absent dynamical coupling to galactic dynamics in our setup results in little large scale motions besides the ones driven by the SNe and a vanishing centre-of-mass velocity.
We investigate the geometric weighted mean of the gas velocity
where we use weights for the velocity of cell , . We distinguish between weighting by the cell mass (), the cell volume (), as well as by chemical mass fractions , , for ionised, atomic and molecular hydrogen. We use the geometric mean in order to avoid a disproportionately strong impact of a few cells with very large values. In addition we correct for outflow motions by investigating only the velocities in and , , and computing an isotropic 3-dimensional velocity . Fig 6 presents the time evolution of the gas velocities (lines) for mass and volume weighting in the top panel and chemical weights in the bottom. The corrected values excluding the outflow are shown in the grey area, which is bounded by the smallest and largest average of the three main simulations. The noticeable differences in the vertical structure and the amount of outflowing gas between the different magnetisations are not reflected in the globally averaged velocities. The temporal variations are larger than the differences between the different runs. The hot gas can increase the average speeds in the volume weighted case to values above . Weighting by mass, atomic and molecular hydrogen shows a very similar evolution with a steady increase from a few to about . The outflow motions have a noticeable effect on the global gas velocities. However, the denser gas (weighting by mass, H, and H) is not strongly affected, which causes the grey area to overlap with the lines. For the low-density gas (weighting by volume and H) the averages for the global and the corrected velocities vary by a factor of two.
In order to relate the velocities to the sonic and Alfvénic speeds we also compute the Mach number,
with the same weights as for the velocity, . Here, is the adiabatic sound speed of cell , where we take the chemical composition into account and compute the cell-averaged thermal pressure according to the local abundances. The corresponding Alfvénic Mach number reads
with the Alfvén speed that depends on the magnetic field strength and the density in cell . The weighting factors are again the same. Here, we directly use the velocities without subtracting a locally averaged mean velocity. The numbers thus include the effects of cloud motions and outflow velocities. The dynamics in the dense gas is investigated in Section 7 in more detail. In Fig. 7 we plot the geometric mean of the sonic Mach number (left) and the Alfvénic Mach number (right) as a function of time. The upper panels show the mass and volume weighted quantities in the entire box. The lower panels depict the Mach numbers weighted by chemical mass fractions. The sonic Mach number indicates transsonic and supersonic motions in all cases. The volume weighted numbers indicate transsonic motions with small variations over time. The mass weighted sonic Mach number increases over time and flattens after at a value of . Due to the time delay in the formation of clouds and dense gas the mass and H weighted Mach numbers differ with magnetisation. The increase over time is relatively similar with a delay of the order of for the magnetic runs. After the curves for molecular gas flatten at a value of . For atomic hydrogen the numbers slightly increase over time to reach with little variations. The ionised gas is slightly supersonic () with very little variations over time and between the non-magnetic and magnetic runs. The similarity between atomic and ionised hydrogen indicates that the faster velocity dispersion in H compensates the hotter temperatures and the resulting sound speed. We note very similar numbers of for weighting by mass and molecular hydrogen as well as for the volume and the H weighted mean.
We measure super-Alfvénic motions with steadily increasing with time. For all weightings the differences between the initial magnetisation is not significant. The mass weighted values quickly reach unity and increase to at the end of the simulation. The volume weighted numbers with overall lower magnetic fields and faster motions in the low-density cells start at values of order 5 and reach after . Contrary to the sonic Mach number the relation between mass and volume weighted numbers and the chemistry weighted counterparts is more complicated. Atomic hydrogen has the lowest Mach number ranging from . Molecular hydrogen corresponds to the densest gas with the strongest fields, which would lower . However, the gravitational contraction and the corresponding acceleration as well as the inter-cloud motions outweigh the stronger field and the mass weighted Alfvénic Mach number increases to . The non-negligible fraction of partially ionised gas with stronger average fields compared to most of the volume results in Mach numbers for the ionised hydrogen of , between the values for molecular and atomic hydrogen. Concerning the volume and H weighted Alfvénic Mach numbers we may be limited by our small box compared to global disc models, which we discuss in Section 8.
As the simulations have multiple collapsing regions, we would like to separate motions related to the individual clouds and centres of collapse. In order to do so, we identify clumps and investigate their physical quantities in detail in Section 7. Here, we illustrate the local dynamics by means of the line-of-sight (los) velocity,
and the corresponding velocity dispersion
We use a uniform grid with the maximum effective resolution of the simulation (, i.e. a cell size of ) and sum up the mass weighted () velocity () of cell along the and direction (). The view along looks very similar despite the fact that the magnetic field is mainly oriented along the direction. Fig. 8 shows the los velocity, the dispersion and the correlation with the column density. The top plot shows simulation B0-1pc, the bottom one B3-1pc both at . The left column depicts the column density projected along (top row) and (bottom row). The central panels show , the right ones . The view along reveals that the largest line-of-sight velocities () are found in regions of low column density, which is consistent with the low-density gas being driven out of the midplane to form a fountain or outflow. Dense regions show smaller values of . For the view along we find overall lower numbers. This is not surprising as we average over multiple independent regions rather than along the stratification of the gas. The velocity dispersion shows a strong inverse correlation with the column density, in particular for the view along . In low-density voids the values increase to . In the regions of molecular clouds can be as low as .
Fig. 9 shows the correlation of the line-of-sight velocity dispersion with the column density of the gas for a projection along all three principal directions (top) as well as the projection along split into the three hydrogen components (bottom) for simulation B3-1pc at . Over time more gas moves to higher column densities but the overall shape of the distribution does not change. The differences between the simulations are minor. The views along and look similar with a broad distribution and a weak anti-correlation of with . A wide spread is not surprising because we integrate over several individually contracting and independently moving regions. For the line-of-sight along we note a clear anti-correlation which we can approximate with a power-law,
independent of magnetisation and time, except for the highest column densities at late times, where gravitational collapse of the clouds keeps the velocity dispersion above , see also Section 7. The scaling indicates a constant kinetic energy per unit area (), which is likely to be connected to our SN driving scheme. We inject thermal energy at a constant rate and in an area-filling manner because the SNe have uniformly distributed random positions in and . The fraction of clustered SNe is not dominating the positioning in the -direction. As a result, the energy injection rate is reflected in the vertical motions.
At the lowest column densities, which correspond to the hot low-density gas, we find , which is in good agreement with the averaged volume-weighted box velocities (see Fig. 6). At high column densities we find in good agreement with the velocity dispersion measured in individually identified clouds, see Section 7.
The strong correlation in for the total gas is not observed, see for example Heyer et al. (2009). However, there are several aspects that weaken this conflict. Firstly, our SN driving scheme is clearly idealised and simplified despite the fact that the spatial distribution and clustering properties are motivated by observations. Secondly, our box is not exposed to the larger galactic environment. Missing components like shear, galactic infall or warp of the disc reduce the overall dynamics within the box, which provides better conditions for a nicely correlated energy input by SNe. The projections along and , which correspond to views much closer to the perspectives in our Galaxy and include integration effects through several independently moving regions, already indicate a weaker anti-correlation. Finally, observations rarely cover all the gas in a continuous fashion, but are rather biased towards a special component with possibly much weaker or vanishing correlation. In fact, split into individual chemical species the correlation is not visible any more. Ionised hydrogen occupies a relatively narrow low- regime with a high line-of-sight velocity dispersion. Atomic and molecular hydrogen cover a significantly broader range in with intermediate (H) and low (H) , respectively.
5.1 Global magnetic properties
Fig. 10 depicts the time evolution of the globally averaged geometric mean of the magnetic field strength in the simulation box,
The weights are again cell mass, cell volume and chemical mass fractions. The averaged mass weighted field starts at a few and slightly increases to values of . The differences in the initial field do not result in a systematic difference throughout the simulation. The H weighted field behaves very similarly with slightly higher values. The field in the atomic hydrogen follows the curves of the mass weighted field and starts decreasing after to . H and volume weighted fields are the weakest decreasing from to and from to The combination of expanding volume filling SN remnants and the ideal MHD approximation are intuitive explanations for the decreasing values over time. However, as for the Alfvénic Mach number, the small box and the missing large scale galactic dynamics might significantly alter the volume filling mean magnetic field, see Section 8. It is interesting to note that the morphological differences between the two magnetic runs are not reflected in the evolution of the magnetic field strength.
Initially, the gas in the midplane with a temperature of is embedded in a magnetic field of and for B3-1pc and B6-1pc, which corresponds to and . Due to cooling most of the dense gas reaches after a few Myr, except for the hot SN remnants. The SNe create volume-filling high- regions in the box. Over time more volume is dominated by hot gas with , only the dense regions remain at low . The compression of the gas into dense regions and the accompanying amplification of the magnetic field results in a similar evolution for the gas with low . However, the Alfvén speed does not increase as much as the gas velocities, which results in very small regions of the simulation box with . A more quantitative view is presented in Fig. 11 showing the total fraction of the gas with (solid lines) and (lashed lines). After of the gas is cold enough to have . This fraction slightly declines to at before it drops perceptibly to at the end of the simulation. The evolution is attributed to the effects of SNe, which do not just create low-density regions with little mass but also heat gas at intermediate densities which contains a noticeable fraction of the total mass and is dominated by thermal pressure. The regions with show a different evolution. The total fraction of gas with decreases quickly to 20 percent of the total mass after leaving a negligible fraction of the mass with sub-Alfvénic motions () at the end of the simulation. However, we should not over-interpret this small fraction. As most of the gas is collected in a small region of the box with limited resolution, we are likely to underestimate the effects of a small-scale dynamo. In addition, the collapse motions around the molecular clouds are poorly resolved.
5.2 Scaling properties with the density
Before discussing the scaling of the magnetic field with the density as well as the angle between the magnetic field vector and the velocity, we would like to stress an important property of the ideal MHD approximation. Fig. 12 illustrates two extreme cases of gas flow with respect to the magnetic field orientation. The top part of the figure presents the initial magnetic field configuration (blue lines), an initial gas cloud (black circle) and a region of fixed size that we investigate at different times (grey circle). Gas motions are indicated by the arrows, where we distinguish between flows along and perpendicular to the magnetic field lines. The centre part depicts the situation for a perpendicular flow. The field lines are compressed with the gas cloud. The observed region thus shows a larger density and magnetic field compared to the initial configuration (, ) in that region. Compression along the field lines (bottom panel) results in a denser cloud with unaltered magnetic field configuration (, ). For very strong magnetic fields the magnetic pressure will channel the flow and restrict it to the parallel configuration. Weak magnetic fields can easily be compressed by the gas flow. It is thus instructive to investigate the scaling of the field strength together with the energy ratios at different densities. Concerning the formation of dense clouds and cores, the distribution of the angle between gas velocity and magnetic fields is a quantitative measure of the accretion mode.
Two dimensional histograms of the gas distribution as a function of density are presented in Fig. 13 for B3-1pc at . The distributions for B6-1pc are very similar. From top to bottom we plot , plasma and . Colour coded is the gas mass. The black lines are the median values along the ordinate. The sonic Mach number has a local minimum at , which marks the low-density (high-temperature) end of the thermally unstable region. At lower densities the gas moves at larger velocities. For the median of the Mach number increases from transsonic motions to , which is mainly due to cooling of the gas and the resulting decrease in the thermal sound speed. Plasma decreases from at to unity at and below unity for densities up to . At the highest densities the thermal pressure increases linearly with the density (isothermal behaviour of the gas). However, the magnetic pressure follows a flatter scaling with (see discussion below), which causes plasma to increase and exceed unity above . The Alfvénic Mach number (bottom panel) shows a similar overall behaviour. At low densities is of order 100. The region of the thermal instability marks the minimum with noticeable fractions of the gas located at sub-Alfvénic speeds and a median of order unity. At high densities the strong increase in kinetic energy density outweighs the increase in magnetic pressure, which causes the ratio to increase again to reach median values of order 10.
We expect the scaling of the magnetic field strength to depend on the dynamical and thermal properties in different density regimes. In the case of very weak fields, turbulent motions and motions driven by thermal pressure (e.g. an expanding SN shell) dominate the evolution of the magnetic field. For compressive motions perpendicular to the field line the field strength in the MHD approximation scales directly with the density with . For flows along the field line the strength is not altered, . Isotropic turbulence statistically results in a scaling of for super-Alfvénic systems (Mestel, 1966). For sub-Alfvénic motions the classical scaling law by Chandrasekhar & Fermi (1953, CF53) results in or (e.g. Mouschovias & Ciolek, 1999).
We analyse the strength of the magnetic field (Fig. 14, top) as well as the fraction of the gas that is magnetically dominated (bottom) as a function of density. The field strength is shown as the median values with errorbars indicating the 25 and 75 percentile of the distribution. The grey points show the low-density range of the observational estimates by Crutcher (2012), which are in good agreement with our simulated values. In the bottom panel we distinguish between gas with and gas with . We also mark the maximum density we can resolve based on the Truelove criterion (Truelove et al., 1997), which states that the Jeans length has to be resolved with at least 4 cells222We compute the maximum density using the isothermal sound speed, a mean molecular weight of 2.35 and a fiducial temperature of 50K. The cell size is taken at the highest level of refinement, i.e. . in order to avoid artificial fragmentation. We compare simulations B3-1pc and B6-1pc at . At low densities we note a steep scaling which can be considered approximately consistent with the statistical average for weak fields. In this density regime only a small fraction of the gas is magnetically dominated. At a density of the scaling becomes flatter with , which coincides with an increase in -dominated gas. Above the relative fraction of gas with drops again, but with motions mainly parallel to the field line, so the field strength does not have to increase, see also next section. The simulations thus do not reproduce the classical scaling law of Chandrasekhar & Fermi (1953), where . However, their model assumes sub-Alfvénic turbulence, which is not the case in our simulations, where we have supersonic and super-Alfvénic motions (see Fig. 7) and resulting MHD shocks, see also de Avillez & Breitschwerdt (2005). Recent ISM simulations by Padoan et al. (2016) show a similar scaling behaviour compared to our runs. They find a very flat scaling exponent of below and above , respectively. Fitting a power law to the data in the density range yields scaling exponents between for simulation B3-1pc and slightly smaller numbers of for B6-1pc, all slopes increasing over time from . A common feature in all of our simulations is a locally flatter scaling around the density of the thermal instability, . A fit over one order of magnitude around reveals scaling exponents that are smaller with ranges of for B3-1pc and for B6-1pc, again measured from . At ISM densities of around the median field strength is of order a few in our simulation, which is in good agreement with the observations by Crutcher (2012). Observationally, the field strength increases to values of up to a few hundred at molecular cloud densities of , consistent with recent other numerical studies (Gent et al., 2013; Kim & Ostriker, 2015b; Körtgen & Banerjee, 2015; Padoan et al., 2016; Iffrig & Hennebelle, 2017; Evirgen et al., 2017) but not resolved in our setups.
We should add a word of caution with regard to the densest regions and the credibility of the strongest magnetic fields because of magnetic dissipation and numerical reconnection. Close to the resolution limit, i.e. in entities that are resolved with cells (see comments in Hennebelle et al., 2011), the magnetic fields are prone to numerical diffusion and (turbulent) reconnection (Lazarian & Vishniac, 1999; Kowal et al., 2009, 2012), which causes flux to be lost. Lazarian et al. (2012) conclude that the reconnection model is consistent with observations and that there is little evidence for subcritical self-gravitating clouds. Although there is gas with low Alfvénic Mach numbers in our simulations, most of the gas is dominated by the kinetic energy density, both infall as well as turbulent motions. This means that the centres of the clouds at the resolution limit are likely to be affected. In which way is difficult to estimate as the loss of magnetic flux will be (partially) balanced by a magnetic dynamo, which is underresolved inside our clouds.
Our measured range in magnetic field strength is consistent with observations in the ISM. Crutcher (2012) reports values ranging from for densities up to . Stronger fields of up to a are only found in gas that is denser than what we simulate. At a density of the median field is , which is consistent with recent observations of Taurus (Chapman et al., 2011), Serpens South (Sugitani et al., 2011), clouds in the Goult Belt (Li et al., 2013), the Lupus I molecular cloud (Franco & Alves, 2015) or the Polaris Flare (Panopoulou et al., 2016).
5.3 Alignment of the magnetic field with the velocity
In order to investigate the direction of the flow with respect to the magnetic field lines, we compute the normalised distribution function of the angle between velocity and magnetic field vector,
As in ideal MHD the sign of the orientation of the field lines does not matter – a flow parallel to the field is numerically and physically the same as an anti-parallel flow – we use this symmetry and plot . In Fig. 15 we show the distribution of angles for the entire simulation box for simulation B6-1pc at . Colour coded are the mass (top panel) and the volume (bottom panel). The solid line marks the median of the distribution along the distribution. The dashed lines are the 25 and 75 percentile. The median varies between 0.4 and 0.6 and becomes flatter over time (not shown). For simulation B3-1pc the variations are similar but less pronounced. Gas at the lowest and highest densities does not move with any preference with respect to the magnetic field orientation. In the range there is a tendency for a flow parallel to the field lines. The denser gas () moves with a slight preference perpendicular to the field lines. This is an effect that coincides mostly with the overall thickening of the disc, the launching of outflows at early times, and gravitational contraction of the disc at late times, where gas is moved along together with the field lines, which are mainly oriented along the -direction. The evolution changes slightly over time with flatter distributions at later times, but the general trend remains. This behaviour is very similar to the results in Iffrig & Hennebelle (2017), see their Fig. 11 and 12. We emphasise that there is a difference between the globally averaged distribution of angles and the distribution in dense clouds, where the overall effect of the thickening of the disc is not visible. There, most of the mass in the accretion flow is funneled along the field lines, which we investigate in more detail in Section 7.
6 Formation of shielded gas and molecular hydrogen
6.1 Dense and shielded gas
This section is dedicated to the details of the formation of dense gas, the correlation with shielding from the interstellar radiation field and the formation of molecular gas. The visual impression of the differences in the column density (cf. Fig. 1 and 2) manifests in a significant delay in the formation of dense gas in the magnetised simulations. In Fig. 16 we present the volume weighted (top) and mass weighted (bottom) probability distribution function (PDF) of the total gas density at (solid lines) and (dashed lines). The volume weighted PDF shows typical features like the broad lognormal-like maxima at and , which is a natural result of turbulent gas dynamics (e.g. Vazquez-Semadeni, 1994; Passot & Vázquez-Semadeni, 1998). In addition there is a high-density excess, which is attributed to the gravitational attraction and local collapse (e.g. Klessen, 2000; Slyz et al., 2005; Kainulainen et al., 2009; Schneider et al., 2012; Girichidis et al., 2014). The earlier snapshot indicates that the distribution in the low-density gas shows similarities in all simulations. The non-magnetic run converts more gas into dense regions with indicated by the broad peak in the mass weighted PDF with an order of magnitude larger masses compared to B6-1pc. After of evolution the PDFs look very similar for all three simulations, so statistically the contraction to high densities seems to only experience a time delay because of magnetic fields, but no significant difference in the PDF after dense gas has started to form.
The formation of dense gas in clouds is accompanied by the formation of regions that are shielded from the interstellar radiation field. In the following we consider gas with an attenuation of as shielded, where is computed in every grid cell by averaging over multiple lines of sight by means of the TreeCol algorithm. This is approximately the threshold for the formation of molecular hydrogen (Röllig et al., 2007; Glover et al., 2010). Fig. 17 shows the column density of the regions with at different times for simulations B0-1pc (left), B3-1pc (centre), and B6-1pc (right). The maps reflect the impact of the magnetic field on the clumping properties. We can estimate the characteristic size of the regions by computing the two-dimensional power spectrum of the column density map in the -plane with , see Fig. 18 for the spectra at and . The corresponding spectra for the -plane look very similar. In the magnetic runs the gas accumulates in larger coherent structures, which quantitatively becomes apparent in the power. At the strongly magnetised run is dominated by the power on the largest scales. For B3-1pc the spectrum is relatively flat. In the non-magnetic case the dense clumps show up as a broad peak at around . At the spectra look very similar with most of the power at spatial scales of . We interpret this position as the characteristic size for molecular structures and as our fiducial cloud radius for the subsequent analysis.
6.2 Correlation between shielded and molecular gas
The total fraction of molecular gas in the simulation box as a function of time is plotted in Fig. 19. The stronger the magnetic field is the longer it takes to form molecular hydrogen. Towards the end of the simulation the majority of the gas in the box is in the form of H. The non-magnetic run reaches a fraction of . The magnetic runs show a delay in the formation of molecular gas, which vanishes after to approach a total fraction of .
Fig. 20 (top) shows the ratio of molecular gas over total shielded gas. The bottom panel shows the fraction of shielded atomic () and molecular gas (). We only plot the ratio if the mass fraction of molecular gas exceeds 10% of the total mass in order to avoid numerical noise in the case of negligible H gas. At early times the ratio (top) is small and increases to reach unity at the end of the simulations. This illustrates that first gas with condenses out of the ISM before molecular gas forms within those regions. At late times the ratio approaches unity in all cases, which means that there is the same amount of molecular and gas with . The bottom panel shows the relative fractions of shielded gas in molecular and atomic gas. The lines for H are close to unity, which emphasises that basically all molecular gas forms in shielded regions. The atomic gas has a more complicated evolution. The non-magnetic run forms of H with at . After that molecular gas forms out of the shielded gas and the fraction declines to give only at the end of the simulation. Weakly magnetised gas (B3-1pc) results in approximately half of the atomic hydrogen being shielded for about before dense condensations leave the majority of the atomic gas in a diffuse state that is not shielded. In the case of B6-1pc the fraction of atomic gas with is initially lower () and increases to a peak fraction of at , which coincides with the onset of H formation. At the end of the simulations the differences are small for all runs.
7 Fragmentation and Cloud Statistics
7.1 Number of fragments
We identify clouds by finding the local minima in the gravitational potential. In order to avoid effects due to initial transient fluctuations, we exclude the first from this analysis. In addition to this necessary criterion, we also require that at the position of the local minima the gas is shielded () and molecular () in order to be identified as a molecular cloud. We have chosen the gravitational potential to identify clouds because it is less prone to temporal fluctuations compared to the density field (Smith et al., 2009). There is a variety of other algorithms to identify clouds and connected structures, among which Clumpfind (Williams et al., 1994) is probably the most popular one, however, with weaknesses (Pineda et al., 2009; Berry, 2015). We apply Clumpfind to our data sets and compare the differences in Section B. However, as it has more free parameters and shows strong dependencies of the results on them, we refrain from using Clumpfind for the main analysis.
Fig. 21 shows the total number of minima in the potential (top, dashed lines) as well as the number with the additional restrictions (top, solid lines). The non-magnetic run quickly develops numerous () local minima in the potential, which merge during the evolution. The magnetised runs show fewer clouds with only and for B3-1pc and B6-1pc. In simulation B6-1pc the number of minima remains constant for about . Due to dissolution and merging of the individual condensations the number decreases during the simulation and drops to about with very little difference between the runs. The number of molecular gravitational minima strongly varies between the runs in the first half of the simulated time. This is expected because the magnetic runs need a longer time before H forms, so the secondary selection criteria are not fulfilled. After the difference between the runs vanishes, see also Fig. 2. We also note that after the formation of even a small fraction of molecular gas the additional constraints are not very restrictive, i.e. the majority of all gravitational minima are molecular and have after .
The differences in the dynamical evolution become apparent if we show the number of molecular gravitational minima as a function global H fraction in the box (see Fig. 21, bottom). The two magnetic runs show a very similar behaviour concerning the total number as well as the evolution with increasing fraction of molecular gas. We note a peak in the distribution at around with clouds. In the non-magnetic run, in which the gas is less filamentary and forms more spherical clouds, the peak is at almost 500 molecular clouds.
7.2 Fiducial cloud radius
In order to analyse the formation of molecular clouds from the very beginning we only use the gravitational minima without the other restrictions of optical thickness and molecular fractions unless explicitly indicated. We investigate a spherical volume of radius around the minima of the gravitational potential as the fiducial cloud size. The motivation for the size is the peak in the spectra of the column density (cf. Fig. 18). However, as the peaks are very broad and we would like to include not only the densest parts of the clouds, we increase the investigated volume and use the spectral peak at as the radius instead of the diameter of the analysis volume. In addition we also probe radii of and , which are presented and discussed in the appendix (Section C). We would like to stress that clouds can encompass other local minima of the gravitational potential besides the one in their centre. Despite variations in cloud sizes, probing radii from seems to be a reasonable choice consistent with observations (Solomon et al., 1987; Bolatto et al., 2008; Miville-Deschênes et al., 2017). Simulations using a very similar setup to ours (Ibáñez-Mejía et al., 2016, 2017) identify cloud radii with more sophisticated methods and find similar sizes (). In addition, the temporal changes of the cloud radii in Ibáñez-Mejía et al. (2016) are less than a factor 2 for an evolution of more than 2 free-fall times after they switched on self-gravity, see dotted lines in their Fig. 7. Also simulations of a global disc (e.g. Tasker et al., 2015) find a distribution of cloud radii with a prominent peak at .
7.3 Cloud masses and dynamics
Fig. 22 shows the mass function of MCs for simulation B3-1pc at different times. The evolution for the other simulations are similar. At the distribution is narrow ranging from and widens to after . Clouds merge and massive clouds accrete gas yielding final masses of up to . Overall the number of clouds is relatively small in particular at the end of the simulation, so a comparison of features of the distribution is inappropriate. Also, the entire gas mass in our simulations is , which is lower than the mass of observed massive GMCs. However, a few global numbers seem to agree with observational data. Miville-Deschênes et al. (2017) find median cloud masses of for all clouds with more massive clouds being located closer to the Galactic centre. Observations of clouds in the Goult Belt show as well median masses of the order of a few (Loren, 1989; Cambrésy, 1999; Pineda et al., 2010; Lee et al., 2015).
Fig. 23 shows the arithmetic mean and median mass of all clouds (top), the mass weighted internal and inter-cloud velocity dispersion (centre) as well as the mass-to-flux ratio (bottom) as a function of time. In all three panels the solid lines are the median of the distribution. The shaded area is bounded by the 25 and 75 percentile. The mean mass in the top panel (dashed lines) is dominated by the most massive clouds and increases from to . The formation of low-mass clouds influences the median mass, which increases by a factor of in the first half of the simulation time. The merging – inparticular of massive clouds in the second half of the simulation time – causes the media to decrease approximately to the initial values at the end. For the velocity dispersion we distinguish between the mass weighted internal dispersion and the inter-cloud velocity dispersion. The former one is computed as
Here, , are the mass weighted velocity dispersions in each direction,
where we sum over all the cells within a radius of around the local gravitational minimum with being the cell mass, the cloud mass, and is the centre-of-mass velocity of the cloud. The inter-cloud dispersion is computed analogously,
Here we weigh by the individual cloud masses and measure the dispersion with respect to the centre-of-mass velocity of all clouds, . The internal velocity dispersions are similar for the magnetic and non-magnetic runs with values of the order of at and at . These averages over the identified clumps with a radius of are in good agreement with the simple line-of-sight estimates (cf. Fig 8). At the velocity dispersion increases dramatically indicating the dominating impact of gravitational collapse, see also the results of high-resolution cloud MC1 in Seifried et al. (2017). The inter-cloud velocity dispersions are overall smaller than the internal values. We find the biggest difference between and in the non-magnetic run, which has the largest internal and the lowest intra-cloud velocity dispersion over the majority of the simulated time. This is in agreement with the morphological evolution in the simulation. The fast formation of dense cold clouds allows for locally more efficient gravitational contraction and stronger accretion flows. However, the clouds are all close to the midplane and the external potential does not accelerate them towards the midplane.
The fact that the inter-cloud velocity dispersion is smaller than the internal one is not surprising considering that we start with the gas being at rest and that the simulations do not include the effects of differential rotation of the galactic disc. However, including the galactic environment still seems to be consistent with our results. Simulations of an entire disc by Dobbs et al. (2015) find cloud merging velocities of , which in agreement of what we find without galactic shear. In addition, those values are computed for directly merging clouds and are likely to be higher than the average inter-cloud velocities. Observations reveal slightly higher velocities () for interacting clouds (Furukawa et al., 2009; Ohama et al., 2010; Torii et al., 2011, 2015; Fukui et al., 2014, 2015, 2016; Saigo et al., 2017). However, in those observations there is also a focus on cloud-cloud collisions, which possibly marks the higher end of inter-cloud velocities. In fact, computing the maximum of all pairwise relative velocities reveals that the clouds are approaching each other with velocities in the range of , in good agreement with the observational values.
where is the gravitational constant. Objects with a ratio below the critical value are supported by magnetic pressure and are expected not to collapse. Above the critical number gravitational forces dominate. For the two magnetic runs the clouds start as subcritical objects with lower values for stronger magnetic fields as expected. After the median mass-to-flux ratio of the clouds exceeds the critical value and grows to values of order at . From first-order theoretical considerations of the ideal MHD approximation, the mass-to-flux ratio can only change by mass flowing along the magnetic field lines. Mass flow perpendicular to the field lines compresses the field and increases the magnetic flux. However, magnetic dissipation and turbulent reconnection can provide an efficient loss of magnetic flux in dense turbulent clouds, which can increase the mass-to-flux ratio without significant accretion. Observationally, there seem to be more supercritial than subcritial clouds (Crutcher, 2012). This is in agreement with the cloud properties in the second half of the simulation. In the first the majority of our clouds is subcritical. This apparent disagreement might stem from the different analysis focus in our simulations compared to the observations. Diffuse structures that clearly appear as gravitational minima in the simulation might still be too diffuse to be classified as proper clouds and thus be excluded as observational targets.
7.4 Mass accretion
Fig. 24 presents the mean value for the accretion rate as a function of time (top) and as a function of cloud mass (bottom) starting from to avoid spurious fluctuations at the beginning of the simulation. The accretion rates of individual clouds show strong fluctuations over time with positive and negative values. Averaged over the mean accretion rate of the cloud ensemble,
slowly increases over time for all simulations with values from , which again agrees well with the high-resolution simulations by Seifried et al. (2017) as well with observations of star-forming regions (Fukui et al., 2009; Kawamura et al., 2009; Peretto et al., 2013). The non-magnetic clouds show noticeably higher accretion rates, which is consistent with the larger mean cloud mass as well as the higher internal velocity dispersion. Plotted as a function of cloud mass (bottom) reveals a scaling of the accretion rate with the cloud mass of approximately , approximately consistent in all simulations. At this point we need to emphasise that the effective accretion rate computed as stated above includes the gain of gas by merging of the clouds and thus cannot be understood as a smooth gas accretion model like, e.g. Bondi accretion (Bondi & Hoyle, 1944; Bondi, 1952), where . Consequently, we cannot apply this scaling to the time evolution of the mass of an individual cloud but have to understand it as the scaling of the ensemble, i.e. statistically a more massive clouds accretes more gas with a relation following , but not smoothly over time.
We use again the angle between the velocity and the magnetic field vector to investigate the accretion flow onto the clouds. A uniform distribution of angles would result in a uniform distribution of . Fig. 25 shows the volume weighted distribution of at for the two magnetic runs, divided into clouds that are subcritical with as well as the supercritical clouds with . Temporal variations as well as variations for different cloud radii are not significant. The distributions show a very similar structure with a peak at indicating that the gas flow is stronger along the field lines than perpendicular to them. The variations across the distribution in are not very large, but nonetheless reveal that the effective integrated accretion from is higher than the one in the range of . Combined with the fact that the accretion rates strongly vary with also temporally negative rates, this could support the scenario, in which the clouds need some time to effectively accrete gas along this parallel channel before becoming supercritical. However, accretion perpendicular to the field lines in combination with magnetic flux loss can as well increase the mass-to-flux ratio. Fig. 26 shows the mass weighted magnetic field strength (top) and the mean mass-to-flux ratio (bottom) of the clouds as a function of their mean density for simulation B6-1pc. Colour coded is the simulation time. The upper panel only includes the data points for a cloud analysis radius of and includes the theoretical threshold line for the critical mass-to-flux ratio. We note that the clouds evolve towards higher masses without significantly increasing the mean magnetic field strength. This supports the accretion scenario, in which the gas flow mainly follows the magnetic field lines. The bottom panel shows the evolutionary tracks for all cloud radii. The critical ratio for is indicated by a dashed line. The differences in the cloud analysis radii mainly manifest in different mean densities. The transition from subcritical to supercritical clouds occurs at approximately the same time ().
7.5 Magnetic flux
An important hint on the dominant accretion scenario is given by the scaling of as a function of , which is equivalent to the same scaling with mass for our fixed cloud radii, . Accretion perfectly along the field lines results in a linear scaling, , because the flux does not change. A flatter scaling with requires a stronger increase in the flux than in the mass. Significant flux losses would dramatically increase . For the scaling is approximately , which is significantly stepper than in previous simulations results, , (Banerjee et al., 2009; Inoue & Inutsuka, 2012; Iffrig & Hennebelle, 2017) and still steeper than recent numerical simulations by Hennebelle (2017), who report a slightly sub-linear scaling. We can relate the effective flux to the mass as and the time derivative
We approximate the time evolution of the cloud mass and the accretion rate as
Assuming a constant we find for the temporal change in the effective magnetic flux
which gives an e-folding time of for . This means that the effective flux in our clouds drops by approximately to over the simulated time.
8 Discussion and Caveats
8.1 Outflows and mass loss
Due to the diode boundary conditions in (gas can leave the box but not enter) we include the effects of outflows but fountain flow activity is thus not captured in this setup. However, our previous studies (Girichidis et al., 2016) as well as other studies (e.g. Hill et al., 2012; Kim & Ostriker, 2018) indicate that the fountain flow cycle occurs on time scales larger than . With our comparably short evolution time of we are not able to capture this effect. In our previous studies the formation of a steady outflow takes for clustered SN driving, which is the most realistic of our employed SN driving modes. Gatto et al. (2017) also report outflows after approximately . This is consistent with the onset of mass loss that we measure in the setups discussed here. Due to the short simulation time we do not simulate the larger box as in previous studies and refrain from investigating outflow properties like temperature, outflow velocities and chemical composition.
8.2 Continuous increase of cloud masses
Clouds in our simulations mainly grow, but are not destroyed by the SN feedback. This is partially due to the positioning of the SNe (60% clustered SNe and 40% individual SNe) at random positions. Thus there is no correlation between the local gas density and the SNe. Therefore, there are no (or very few) explosions that occur directly in the centre of molecular clouds. However, previous studies with a very similar setup (Walch et al., 2015; Girichidis et al., 2016) show that the clustered SN driving results in the most realistic ISM conditions. Studies by Li et al. (2015) elucidate that the fraction and distribution of runaway stars justifies a SN positioning that is not directly related to the parental density peaks of the corresponding molecular cloud core.
Apart from the positioning problem of SNe our study is missing other physical effects like direct radiation from massive stars and stellar winds. Both effects have been included into similar setups using sink particles (Hennebelle & Iffrig, 2014; Gatto et al., 2017; Peters et al., 2017). All simulations including sink particles are able to (partially) disperse molecular clouds. However, the interaction between the gas and sink particles should cover at least a few cells in radius on the highest refinement level, which corresponds to a region of with our current resolution. In particular the lack of a reliable prescription for the magnetic interaction between the Eulerian grid and the Lagrangian sink particles dissuades us from including sink particles in this study.
Numerical experiments on scales of a full galaxy also find similar effects concerning the evolution of GMCs. Tasker & Tan (2009), Tasker (2011), Renaud et al. (2013), Bournaud et al. (2014) and Tasker et al. (2015) simulate full galactic discs with varying ISM and feedback models. They all conclude that feedback is not able to disperse massive condensations.
Detailed simulations of isolated clouds including internal feedback processes show that the destruction or dispersal of molecular clouds strongly depends on the cloud and feedback details. Banerjee et al. (2007) show that jets and outflows from individual sources are not able to disrupt clouds because the energy injection is primarily in compressive modes that are damped too efficiently. Outflows from multiple sources stir the gas more efficiently but still do not provide enough energy transfer to disrupt the clouds (Wang et al., 2010; Nakamura & Li, 2014; Peters et al., 2014). Stellar winds and ionising radiation are more efficient in removing gas from molecular clouds. However, ionising radiation and the resulting thermal pressure are able to change the internal structure but are unlikely to be able to disrupt clouds more massive than (Krumholz & Matzner, 2009; Murray et al., 2010). Numerical simulations by Walch et al. (2012) and Dale et al. (2012) confirm this limited dynamical effect. Simulations including stellar winds allow to unbind of the mass (Dale et al., 2013, 2014) but without completely disrupting them on time scales of . Both radiation and stellar winds are able to perforate the cloud’s envelope, see e.g. the one-dimensional study by Rahner et al. (2017). As a result, SNe placed in the centre of molecular clouds have a limited effect in disrupting them because of the high porosity and efficient radiative cooling (Rogers & Pittard, 2013; Walch & Naab, 2015). The missing central stellar feedback in our simulations might thus change the structures in and around the clouds (see also Gatto et al., 2017; Peters et al., 2017; Butler et al., 2017), but would likely allow the massive clouds to survive.
8.3 Missing galactic shear
We do not use shearing boundary conditions in our setup. For the evolution over long time scales shear might be an important effect for both GMC evolution as well as the global magnetic field properties. However, estimates based the Rossby number for this setup in Girichidis et al. (2016) suggests that the simulated setup is not dominated by shear effects. Concerning the magnetic dynamo properties, simulations by Piontek & Ostriker (2005), Gent et al. (2013), and Kim & Ostriker (2015b) in a similar simulation box including shear indicate that the time scales for an effective dynamo is much longer than the simulated .
8.4 Influence of the galactic environment
More important than the dynamo effect due to large scale galactic dynamics might be the missing galactic large-scale structure that anchors the field at the boundaries of our box. The efficient SN driving can easily drive magnetic flux out of the box, which would be more difficult in a full galactic context. Although this effect is unlikely to change our results for the dense regions and the cloud analysis, it might help the magnetic field to persist in the low-density environment. As a result the field strengths and the Alfvénic Mach numbers in the low-density gas (volume weighted curves of in Fig. 10 and in Fig. 7) might be lower and upper limits, see also observations by Heiles & Troland (2005) who find magnetism and turbulence to be in approximate equipartition. Observations summarized by Haverkorn (2015) also report strong fields above the disc in the Milky Way, which would correspond to significantly lower Alfvénic Mach numbers assuming that the overall dynamics in our simulations is comparable to the observations.
8.5 Long-term effect of magnetic fields
We primarily note differences between the magnetic and non-magnetic runs in the first half of the simulation time. Towards the end the differences significantly reduce. That suggests that the magnetic fields might not have a long-term effect on the evolution of the ISM. However, the initial differences are significant. In particular the formation of dense gas and the onset of the first gravitationally collapsing clouds are significantly delayed if magnetic fields are included. As we are not following the star formation based on dense or collapsing gas but apply a constant SN rate, we are not including possibly reduced star formation and SN rates in our model. Whether (temporally) reduced feedback will enhance or further reduce the impact of magnetic fields will depend on the details of the feedback like the positioning of the SN and cannot be answered with the current model.
We perform three-dimensional simulations of the SN-driven magnetised interstellar medium with focus of the formation of molecular clouds in the presence of magnetic fields. The simulated volume covers a stratified gas distribution along with a magnetic field strength of , , and in the midplane (). The total volume covers with a gas surface density of . We employ a chemical network including ionised, atomic and molecular hydrogen as well as CO, C and free electrons. The effects of (self-)shielding and the attenuation of the interstellar radiation field are taken into account via the TreeCol algorithm. We evolve the system for a total time of with two different effective resolutions ( and ). The results can be summarised as follows.
In the presence of magnetic fields the morphological evolution of the ISM differs from non-magnetic environments. The magnetic pressure thickens the galactic disc leading to much larger scale heights ( for the total gas, for molecular gas) compared to the non-magnetic environments, where 90% of the total gas is confined to and 90% of the molecular gas is located at a height of only from the disc midplane. As a result the non-magnetic simulation forms small clouds with a spherical shape close to the midplane. The thick magnetised disc allows the formation of structures at larger heights above the midplane. In addition the magnetised structures form as more massive and elongated entities.
Both the direct support by magnetic pressure as well as the resulting lower central densities due to the extended disc structure retard the formation of dense structures and molecular clouds. In the non-magnetic simulation molecular gas starts forming after only and reaches a total mass fraction of 0.5 after . In contrast the simulation with an initial magnetic field strength of needs before significant amounts of H form. Once the formation of molecular hydrogen sets in the formation rates are comparable. At the end of the simulation we do not notice a difference in total H fraction between the non-magnetic and magnetic runs.
Most of the gas is moving at supersonic and super-Alfvénic velocities. The average Mach numbers in the ionised hydrogen are around with very little temporal variation and negligible difference between the magnetic and non-magnetic simulations. The Mach numbers in atomic hydrogen increase over time from to values of , again with little dependence on the magnetisation. The highest Mach numbers () are reached in the molecular gas. The final values are very similar for all simulations but due to the different formation times of molecular gas, the intermediate evolution differs with temporally lower values for the magnetic runs. The globally averaged Alfvénic Mach numbers indicate super-Alfvénic motions for both magnetic simulations. Over time increases to values of for the molecular gas and for H. The lowest Alfvénic Mach numbers develop in atomic hydrogen with .
The mean magnetic field increases from the initial values of and to in the dense gas. The low-density environment loses magnetic intensity to minimum values of . The low-density gas with weak magnetisation shows a scaling of the field strength with density of order with as expected from dimensional arguments for a negligible magnetic field impact. At densities of order the mean ISM density () the field becomes stronger such that of the mass is in regions with and a fraction of with . Above this density the scaling of the field becomes flatter with , which indicates that the magnetic field is able to channel the flow along the field lines suppressing the compressional effect.
The non-magnetic simulations form significantly more molecular clouds (peak of 500 at a global molecular gas fraction of ) compared to the magnetic simulations with a maximum of clouds at approximately the same total molecular gas fraction. Over time the clouds merge and accrete material, which allows them to grow from an average mass of to . The mean accretion rate scales approximately with the mass of the cloud as , independent of the magnetisation of the gas. The magnetised initial clouds form with a subcritical mass-to-flux ratio ( in units of the critical value), i.e. are supported by against gravitational collapse. Accretion flows that are stronger along the field lines than perpendicular to them allows the clouds to transition from sub to supercritical clouds at a time of . At the end of the simulation the clouds reach median mass-to-flux ratios of more than times the critical value. This transition is a potentially important process to delay the collapse of clouds and the formation of stars and might alter the star formation efficiency.
The internal velocity dispersion of the clouds ranges from and constantly increases over time. The inter-cloud velocity velocity dispersion is slightly lower with values from and also increases. The non-magnetic run shows the largest internal and smallest inter-cloud values. For the simulation with strong magnetic fields both quantities are the same within the temporal scatter.
We thank the referee for helpful comments and suggestions that improved the paper. PG, DS, SW, TN, SCOG, and RSK acknowledge support from the DFG Priority Program 1573 Physics of the Interstellar Medium. PG acknowledges funding from the European Research Council under ERC-CoG grant CRAGSMAN-646955. DS and SW acknowledges the support of the Bonn-Cologne Graduate School, which is funded through the Excellence Initiative and the German Science Foundation (DFG) for funding through the Collaborative Research Center 956 The conditions and impact of star formation, project C5. SW further acknowledges the support from the European Research Council through the ERC Starting Grant RADFEEDBACK (project number 679852) under FP8. TN acknowledges support from the DFG cluster of excellence Origin and Structure of the Universe. RW acknowledges support by the Albert Einstein Centre for Gravitation and Astrophysics via the Czech Science Foundation grant 14-37086G and by the institutional project RVO:67985815 of the Academy of Sciences of the Czech Republic. RSK and SCOG thank the DFG for funding via the SFB 881 The Milky Way System (subprojects B1, B2, and B8). RSK furthermore acknowledges support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant STARLIGHT (project number 339177). The authors thank the Max Planck Computing and Data Facility (MPCDF) for computing time and data storage. The software used in this work was developed in part by the DOE NNSA ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago. Parts of the analysis are carried out using the YT analysis package (Turk et al. 2011, yt-project.org).
Appendix A Numerical Resolution & Refinement Threshold
We compare two different maximum refinement levels, corresponding to a minimum cell size of and . In all simulations we refine based on the density as well as on the mass fraction of H. For the simulations B0-1pc, B3-1pc and B6-1pc we refine the grid at a threshold H fraction of and derefine at a value of . We perform two control runs (B0-1pc-ht, B3-1pc-ht), where we use the values of for refinement and for derefinement. Fig. 27 shows the total H content in the box over time for all simulations. In all cases the differences due to the magnetic field are significantly larger compared to the refinement criterion or the maximum refinement level. However, we note that the interior structure of the molecular clouds is not resolved. The weak differences for the two spatial resolutions is thus not comparable to the much larger resolutions in Seifried et al. (2017) in a very similar setup because their study resolves structures down to .
Fig 28 shows the number of fragments (number of local gravitational minima) as in Fig. 21 for the (solid lines) and (dashed lines) resolutions. In all cases the higher resolution results in approximately twice as many local gravitational minima (upper panel). Including the restrictions of optical thickness () and molecular composition () does not fundamentally change the behaviour (bottom panel). At early times the magnetic high resolution runs form up to an order of magnitude more molecular clouds but this difference reduces at later stages of the simulations.
Appendix B Defining Clouds
We have used local minima of the gravitational potential to identify clouds. However, there are also other methods for defining clouds. One popular method is Clumpfind (Williams et al., 1994), which investigates contours of the density to define connected clouds. Alternatively, friends-of-friends methods or dendrograms have been used to find groups of connected regions. One difficulty with other methods are the additional free parameters. In this section we discuss the differences between our method and the classical Clumpfind algorithm.
For Clumpfind, we vary the minimum number of cells that are needed to identify clumps () as well as the minimum average density of the clump (). In all cases we set the density ratio between the contours to 10. Fig. 29 shows the column density of the gas together with the identified clumps. The top plots compare different required to be identified as clump. The bottom one depicts the differences based on minimum average densities, , with a fixed . The left-hand column in both parts shows the clumps based on the local minima of the gravitational potential and is identical in both parts. The three right-hand columns in the top plot show the clumps for , , and . The number of identified clouds has a large range from . The bottom part shows the number of clouds with and the additional constraint of , , and from left to right. The number of clumps is much less dependent on the minimum average density. For and there is reasonable agreement between the clumps identified by the gravitational potential and Clumpfind for the most-massive clumps.
Fig. 30 shows the time evolution of the total number of clouds (top), the median mass (centre) and the derived size (bottom) based on spherical symmetry, . The lines indicate the Clumpfind results; the shaded area shows the numbers for the gravitationally found clumps and is bound by the smallest and largest median of the three main simulations. The left-hand panels are for different , the right ones for different with a fixed . The number of clumps is a strong function of , which indicates that most of the identified structures are small density fluctuations. None of the parameters yields a comparable evolution to found for the gravitationally detected clumps. The median masses even vary by three orders of magnitude when changing from to . Here, is comparable to our preferred method within a factor of a few. The radius also reflects the strong dependence on , which is smaller than our fiducial radius of . To investigate the effects of we set , which corresponds to a radius of at the highest level of refinement. The variations for different minimum densities are overall smaller. The total number of clouds shows a similar qualitative time evolution as for our main method. The masses are smaller except for the highest minimum density, which only extracts the densest gas as clouds. On one hand this suggests that the properties the found clumps are not a strong function of . On the other hand this reflects that all numbers are again set by