EDP Sciences
Free Access
Issue
A&A
Volume 604, August 2017
Article Number A70
Number of page(s) 22
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201630290
Published online 08 August 2017

© ESO, 2017

1. Introduction

Understanding the cycle of matter and energy within galaxies is a necessary step for our knowledge of the structure formation in the universe. Because of the wide diversity of spatial and temporal scales, which govern these cycles, it is not possible to directly simulate a galaxy (e.g. Tasker & Bryan 2006; Dubois & Teyssier 2008; Bournaud et al. 2010; Kim et al. 2011; Dobbs et al. 2011; Tasker 2011; Hopkins et al. 2011; Renaud et al. 2013), with a well-resolved interstellar medium (ISM) although the spatial resolution is continuously improving. To address this question, an alternative approach has been developed which consists of simulating a small portion of a galactic disc leading to a better spatial resolution (Korpi et al. 1999; Slyz et al. 2005; de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; Hill et al. 2012; Kim et al. 2011, 2013; Gent et al. 2013; Hennebelle & Iffrig 2014; Gatto et al. 2015) although at the expense of solving the large galactic scales. Clearly these two approaches are complementary and must be used in parallel.

Generally speaking the most important motivation for performing this type of calculations is to study the impact of various feedback mechanisms such as supernovae on the star formation rate (SFR), the structure of the galactic disc and the outflows that are launched. Due to technical difficulties, in particular the very small time steps induced by the high velocities and temperatures that feedback generates, these studies are usually limited to relatively modest numerical spatial resolution and as a consequence the detailed flow properties have been more rarely considered (de Avillez & Breitschwerdt 2007; Padoan et al. 2016), being mainly addressed through smaller scale simulations in which the energy input has to be prescribed (Vázquez-Semadeni et al. 2006; Kritsuk et al. 2007; Banerjee et al. 2009; Audit & Hennebelle 2010; Inoue & Inutsuka 2012).

In this paper we continue from our previous study (Hennebelle & Iffrig 2014) by performing a series of high-resolution simulations varying the magnetic intensities from purely hydrodynamical to strongly magnetized flows. In order to better constrain the influence that feedback exerts onto the galactic disc evolution, we also reduce its strength and compare with the fiducial feedback values. Our aim is to characterize not only the SFR and disc structure but also the properties of the turbulence which develops. Indeed, the flows in such an environment are complex because not only they are self-gravitating and magnetized but also because energy is injected at intermediate scales through the stellar feedback processes, and last but not least, because the galactic discs are strongly stratified making them highly non isotropic. Characterizing the turbulence and more generally the physical and statistical properties of these flows is a necessary step to understand the star formation process and more generally the energy cycles that are taking place in the galaxies. Therefore we calculate the various bidimensional power spectra at various altitudes finding for some of them a strong dependence with it. We also investigate the alignment between the velocity and the magnetic field. Indeed, not only this latter may be playing a role in the theory of MHD turbulence (Boldyrev 2006) but also the strong influence it may have on the molecular cloud evolution and collapse has been recently stressed by Inoue & Inutsuka (2009) and Körtgen & Banerjee (2015). Finally, we extract the clumps and study their properties in particular how they depend on the magnetization and the feedback strength.

In the second section of the paper we describe the numerical setup that we use. In the third section we give a general description and analysis of the mean quantities such as the star formation rate as well as the thickness of the galactic plane. In the fourth section, we study the cold clumps extracted from the simulations and discuss their physical properties. The fifth section is dedicated to the physical properties of the turbulence which develops in our simulations and we study the power spectra, the alignment between the velocity and magnetic field as well as the structure properties. The sixth section concludes the paper.

2. Numerical setup

2.1. Numerical code and resolution

We ran our simulations with the RAMSES code (Teyssier 2002; Fromang et al. 2006), an adaptive mesh code using a Godunov scheme and a constrained transport method to solve the MHD equations, therefore ensuring the nullity of the magnetic field divergence. We used a 5123 or 10243 grid with no adaptive mesh refinement in order to get consistent power spectra: the interplay between an adaptive mesh and numerical diffusivity may introduce non-trivial biases into the spectra on a large range of scales because of variable spatial resolution. This uniform grid has a cell size of approximately 1 or 2 pc.

Because of the small time steps induced by high velocities and temperatures generated by supernova explosions, these simulations are quite demanding. They required a total of approximately 15 million CPU hours on a BlueGene supercomputer.

2.2. Physical processes

Our simulations include various physical processes known to be important in molecular clouds. We solve the ideal magnetohydrodynamics (MHD) equations with self-gravity and take into account the cooling and heating processes relevant to the ISM. We also added an analytical gravity profile accounting for the distribution of stars and dark matter. The corresponding gravitational potential is (Kuijken & Gilmore 1989): (1)with a1 = 1.42 × 10-3 kpc Myr-2, a2 = 5.49 × 10-4 Myr-2 and z0 = 180 pc, as used by Joung & Mac Low (2006).

The equations we solved are with ρ, v, P, B, Φ, and E respectively being the density, velocity, pressure, magnetic field, total gravitational potential, and total (kinetic plus thermal plus magnetic) energy. The loss function , includes UV heating (kept constant and uniform through the simulations) and a cooling function with the same low-temperature part as in Audit & Hennebelle (2005) and the high-temperature part based on Sutherland & Dopita (1993), resulting in a function similar to the one used in Joung & Mac Low (2006). The gravitational potential Φ has two terms as stated before: the one due to stars and dark matter φext, and the one due to the gas itself φ, hence Φ = φ + φext.

2.3. Initial and boundary conditions

We initialized our simulations with a stratified disc: we used a Gaussian density profile: (7)where n0 = 1.5 cm-3 and z0 = 150 pc. This leads to a total column density, Σ, through the disc that is equal to where ρ0 = mpn0 and mp is the mean mass per particle. We get Σ = 4 × 10-3 g cm-2 = 19.1M pc-2. For convenience, we defined the scale height of the gas .

The temperature was set to an usual warm neutral medium (WNM) temperature, around 8000 K. In order to prevent this disc from collapsing, an initial turbulent velocity field was generated with a RMS dispersion of 5 km s-1 and a Kolmogorov (Kolmogorov 1941) power spectrum with random phase. We added an initial horizontal magnetic field: (8)with B0 ≃ 0,3,6or12 μG for different runs (see Table 1). We note that our setup was not designed to provide a detailed equilibrium along the z-axis since the latter must result from a self-consistent feedback and star formation cycle. The boundary conditions are periodic in the x and y-directions and outflowing vanishing gradient along the z-axis, that is to say that gas can leave the box but cannot enter.

2.4. Supernova feedback

The feedback scheme used for this set of simulations is similar to the fiducial run (C1) in Hennebelle & Iffrig (2014). We used Lagrangian sink particles (Krumholz et al. 2004; Federrath et al. 2010a; Bleuler & Teyssier 2014) to track star-forming regions. The specific implementation of Bleuler & Teyssier (2014) uses a reconstruction of the clumps on top of which sinks are placed if the clump is gravitationally unstable and collapsing. When such a sink particle accreted more than 120 M of gas, we assumed a massive star has been formed, and we triggered a supernova within a 10 pc radius around the sink particle, the exact position being determined randomly. This was done by injecting a few 1043 g cm/s (Iffrig & Hennebelle 2015a) of radial momentum (later denoted pSN, see also Table 1) in a sphere (a few computational cells in radius) around the chosen location. The thermal energy was however saturated to a corresponding temperature of 105 K to prevent very high sound speeds outside the galactic plane, as such speeds enforce a very small time step, and such a run would be computationally too expensive (the number of required time steps is an order of magnitude higher). In the same way, we also limit the velocity produced by the supernovae to a maximum value of 200 km s-1. We have tested the influence of this threshold by increasing its value and find that it does not have a strong influence. These approximations give reasonable results for the gas in the disc, but the low-Mach outflowing gas is not treated correctly essentially because the hot phase produced by supernova explosions, which is largely responsible of gas expulsion, is absent from these simulations.

We recall that the exact way supernova feedback is implemented has been found to have a drastic influence on the results (Hennebelle & Iffrig 2014; Gatto et al. 2015). In particular if the supernovae are not sufficiently correlated with the dense gas, they do not exert any substantial influence on the star forming gas and this leads to very high SFR. Since treating precisely the dense gas and the supernova correlation implies resolving not only the star formation well but also following the detailed star trajectories, this constitutes a very difficult task, that renders prescriptions like the one we are using unavoidable. More generally, other sources of feedback such as HII radiation and stellar winds should be considered as well (Walch et al. 2012; Dale et al. 2013, 2014; Geen et al. 2015, 2016).

2.5. Runs performed

Table 1

Summary of the runs.

The different runs are summarized in Table 1. To understand the influence of the magnetic field we considered four values of the initial magnetic intensities, while to study the impact of feedback we performed runs with two values of pSN, the momentum injected by supernovae. The first value, 4 × 1043 g cm-1, corresponds to the typical momentum thought to be injected in the ISM by supernovae. We will refer to this as “standard”. We will refer to the second value of 1043 g cm-1 as “weak”. By combining low resolution (labeled “L”) and high resolution runs, we can check for numerical convergence. All simulations, except B0, have been run up to 80 Myr. Apart from the weak feedback cases, our simulation are similar to the ones performed by Kim & Ostriker (2015b) except that we did not treat the shear and we have higher numerical resolution.

3. Global properties

In this section we give a general description of the simulations and quantitative estimate of the SFR and disc scale height.

3.1. Qualitative description

Figure 1 shows column density maps for the various runs. Apart from the highest resolution, the general features are similar to the ones of the lower resolution runs presented in Hennebelle & Iffrig (2014) and also seen in other works (e.g. Kim et al. 2011; Walch et al. 2015; Kim & Ostriker 2015a,b). In particular, a stratified multi-phase disc developed. As expected, its density profile varies with magnetic field and feedback strength. By comparing the snapshots at 40 and 80 Myr, it is clear that the disc profile significantly evolves with time and becomes progressively thinner, a point that is closely analysed below.

thumbnail Fig. 1

Column density maps for the various runs. From left to right: runs B0, B1, B2 and B4 (increasing initial magnetic field). From top to bottom: fiducial pSN = 4 × 1043 g cm/s at 40 and 80 Myr, reduced pSN = 1043 g cm/s at 80 Myr. First row: edge-on view. Second row: face-on view.

Open with DEXTER

thumbnail Fig. 2

Mass in sinks (first row) and estimated star formation rate (second row) as a function of time. The dotted line for the sink mass corresponds to the initial gas mass in the simulation box. The dotted line for the star formation rate corresponds to the typical observational value (e.g. Kennicutt et al. 2007) given our column density. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER

The most visible effect of the variation of the initial magnetic field intensity is the change of the disc thickness at time 40 Myr. While this is expected because qualitatively a stronger magnetic field entails more magnetic support, and therefore the disc should be thicker, we will see later that quantitatively this may not be so simple. Another effect of the magnetic field is also that it reduces the star formation rate and delays the beginning of the star-forming phase. This explains the thicker disc in the B0L case, where star formation starts earlier, which in turn implies more feedback from supernovae. One can also see that the medium is rather filamentary in the magnetized run, whereas it looks clumpy in the hydrodynamical case (Hennebelle 2013).

Comparing the middle and bottom panels of Fig. 1, the effect of reduced momentum input is relatively clear: there is less support against gravity, hence a thinner disc. The horizontal structure is also affected, resulting in a dense structure surrounded by diffuse medium, whereas a stronger supernova feedback disperses the high-density clumps more easily, resulting in lower column density contrasts. Moreover we see that a very dense region has developed in runs B0W (at x ≃ 250 pc and y ≃ −200 pc) and B1W and is probably developing in runs B2W and B4W. This shows that global collapse is going on because feedback is too weak to prevent it.

thumbnail Fig. 3

Magnetic field in the equatorial plane as a function of time for the three runs B1L, B2L and B4L. Solid line is the x-component while dotted and dot-dashed are the y and z-components respectively. As can be seen the x-component decreases with time which indicates that the initial magnetic field gets expelled from the galactic plane.

Open with DEXTER

3.2. Star formation history

The star formation rate is estimated by following the mass accreted by the sink particles as a function of time. Figure 2 displays the total mass and the SFR (obtained by taking the time derivative of the sink mass with respect to time). For reference the corresponding estimate is given for a typical Milky-Way type galaxy. Typical SFR in the simulations are on the order of a few 10-2M yr-1 kpc-2. Given the total column density of about 19 M pc-2, these SFR appears to be comparable to observed values (e.g. Kennicutt 1998; Kennicutt et al. 2007; Bigiel et al. 2008).

As expected, Fig. 2 shows that the SFR is reduced in the presence of a stronger initial magnetic field by a factor of up to four in the standard feedback case and ten in the weak feedback one. Another effect of stronger fields is that star formation is delayed by roughly 30 Myr between runs B1L and B4L. While these effects are substantial, one may speculate that they could be even stronger and that star formation could, in principle be almost entirely quenched. This is because for a 10 μG magnetic field, the Alfvén speed within the WNM is about 20 km s-1 and therefore between two and three times larger than the sound speed that is equal to about 7 km s-1. Therefore the magnetic Jeans mass is 8–20 times the thermal Jeans mass. To investigate this, Fig. 3 displays the mean value of the magnetic intensity in the equatorial plane for the runs B1L, B2L and B4L. As can be seen for the runs B1L and B2L the magnetic intensity at 30 Myr is above their initial values (respectively about 2.5 and 5 μG), this is due to the initial contraction of the galactic plane and the magnetic flux freezing, that tends to amplify the magnetic field. Clearly the Bx component decreases with time and comes close to the initial field intensity for run B1L and B2L. For run B4L, it falls below (by almost a factor of two) the initial intensity. This indicates that some magnetic flux is expelled from the galactic plane and does not remain frozen into the gas. Let us stress that the galactic disc is clearly contracting after 40 Myr and nearly stationary between times 60 and 80 Myr (see bottom panel of Fig. 4), therefore the flux decrease is not a consequence of global expansion of the disc. It is likely a consequence of the turbulence, which triggers reconnection and magnetic diffusivity a process largely observed in other contexts (Lazarian & Vishniac 1999; Joos et al. 2013). On the contrary, in the three cases the By and Bz components are amplified likely because of turbulence, to values on the order of 1–2 μG. To conclude, magnetic field does reduce the SFR but its impact is mitigated by turbulent magnetic flux leakage. We stress that since shear is not included in our study, we may neglect or underestimate the generation of magnetic field.

In the weak feedback runs (bottom panels of Fig. 2), the SFR is higher than for the standard feedback case by a factor of between three and ten (depending of time and runs). It is interesting to note that contrary to the standard feedback case, where a plateau is reached, the SFR keeps increasing with time indicating that no stationary behaviour has been reached. This is likely because the feedback is too weak and therefore the galactic disc is undergoing a catastrophic collapse as it is indeed suggested by Fig. 1.

We conclude that the scheme employed here is able to reproduce the observed SFR, which is a fundamental property for galaxies. While the total amount of momentum delivered in the flow is roughly correct, we caution that changing the correlation between supernovae and dense gas may alter this conclusion. This, however, suggests that the present setup can be employed to study further the flows and the dense structures arising in galactic discs.

thumbnail Fig. 4

Vertical density profiles. Top: disc profile for the different runs with strong feedback, at 40 Myr. Bottom: evolution with time in the B1 and B1L runs. The dashed lines are the lower-resolution profiles.

Open with DEXTER

3.3. Vertical structure

3.3.1. Density and pressure profiles

In order to study the vertical structure of the simulated galactic disc, we compute vertical profiles by averaging over horizontal slices. The resulting density profiles are shown in the top panel of Fig. 4. The corresponding full widths at half maximum are respectively 160, 110, 110, and 240 pc for runs B0L, B1L, B2L, and B4L at 40 Myr. These values are comparable to the observed values for the Milky Way: 120 pc for the molecular gas and 230 pc for the atomic gas (Ferrière 2001), although there is no distinction between those components in the simulations. In an attempt to better quantify the scale height in our simulations, we have calculated the mean height per density bin and displayed the corresponding result in Fig. 5. Using the definition of (Ferrière 2001; see her Eq. (1)), which is slightly different from the one we use, and the values Hm = 71−81 pc that she quotes, this corresponds for the molecular gas to a mean height of about 40−45 pc and 100 pc for the HI. These numbers are in reasonable agreement for the gas of densities between 10 and 100 cm-3 and less dense than 1 cm-3 respectively. They are about two times larger for the gas denser than 100 and 1 cm-3 respectively. It is therefore hard to conclude whether the disc is too narrow by a factor of two or compatible with the observations as this would imply a good description of the molecular gas in the simulation (as well as of the CO emission), which in particular likely requires a better spatial resolution. We stress that our simulations do not include the cosmic rays, which likely contribute to support the galactic disc against gravity. In particular Girichidis et al. (2016) have recently performed simulations which suggest that this cosmic rays could indeed have a significant contribution, mainly because they dissipate less easily that the hot gas produced by supernova explosions. Definite conclusions are, however, prevented given the difficulties of accurately measure this scale height.

thumbnail Fig. 5

Mean height as a function of gas density for runs B0L, B1L, B2L and B4L at time 60 Myr.

Open with DEXTER

The bottom panel of Fig. 4 shows the density profiles at various time steps for run B1. As can be seen it follows a complex evolution. First of all between 20 and 40 Myr, the disc starts expanding. This is a consequence of the increasing star formation between these two time-steps. Then the disc contracts and finally reaches an equilibrium, as shown by the comparison between times 60 and 80 Myr. The comparison between runs B1 and B1L shows that numerical convergence seems to have been reached, at least for this particular aspect.

To get a better understanding of the vertical equilibrium, we show the various pressures, namely thermal, kinetic and magnetic across the galactic disc and for the four simulations B0, B1, B2L and B4L about 25 Myr after star formation has begun. As found in other studies (e.g. Kim et al. 2013; Hennebelle & Iffrig 2014), we find that the thermal pressure is negligible while the kinetic one tends to be dominant. We find that the magnetic pressure is typically comparable to the kinetic one except in the most magnetized case for which it dominates. While the values at z = 0 are relatively similar, the profiles present significant differences. In particular the stronger the magnetic intensity, the smoother its distribution. For example while for the B1 run the magnetic pressure falls to very low values above z = 300 pc, at this height it is still one tenth of its value at the origin for run B4L. This clearly shows the magnetic flux expulsion discussed in the previous section (Fig. 3). In particular these profiles do not correlate well with the density ones while flux freezing would imply strict proportionality.

3.3.2. Analytical modelling

Although the various profiles display a great complexity, it is worth building an analytical model to get a better quantitative insight. For this purpose we assumed that the disc is stationary. While this is obviously not the case at early times, Figs. 24 show that this is a very reasonable assumption at time 60 and 80 Myr for the standard feedback case. We can derive an estimate for the disc thickness using the approach of Kim & Ostriker (2015b, Sect. 2). By averaging the momentum conservation Eq. (3)in time and horizontal directions (see e.g. Boulares & Cox 1990), and integrating it between the mid-plane and some altitude zmax, where we assume the density to be negligible, we can write: (9)where σz is the vertical component of the turbulent (thermal support being neglected) velocity dispersion, ΔΠmag is the difference of magnetic support between zmax and the mid-plane, is the weight of the gas disc in its own gravitational potential, and is the weight of the gas disc in the external potential. Expressions for these quantities will be derived in the following paragraphs.

In order to estimate σz, we assumed that a stationary energy equilibrium is established between on one hand the supernova energy injection and on the other hand the turbulent energy dissipation, we write (10)where SN is the number of supernova explosions per unit of time and per unit area that we assumed to be equal to ΣSFR/ 120 M that is to say the mass of stars produced per unit of time divided by the mass necessary to get a massive star (assuming a canonical Salpeter IMF), ESN is the energy of the supernova taken to be equal to 1051 erg and ϵSN is the fraction of this energy that is effectively communicated to the gas. As a canonical value we adopted 5%. Then finally, ηturb is a coefficient of the order of one, that will be adjusted. It reflects our ignorance of the exact way that turbulence is decaying, in particular since the galactic disc is stratified, the choice of the crossing time, while proportional to R/σ is not straightforward. It also depends on the exact efficiency of the energy injection due to the supernovae since ϵSN is difficult to estimate accurately.

thumbnail Fig. 6

Vertical kinetic, thermal and magnetic pressure profiles 25 Myr after the beginning of star formation. From top to bottom: runs B0L, B1L, B2L and B4L.

Open with DEXTER

We further assume energy equipartition between all three directions of the velocity. Thus, we can write (11)where is the total velocity dispersion in the mid-plane. While isotropy may not be entirely achieved given the strong stratification, in the context of protostellar formation Lee & Hennebelle (2016a,b) have investigated the differences between a model which assumed isotropy and a model which considered a smaller velocity dispersion along the z-axis. They found that these differences remain limited.

Given our initial magnetic field distribution, we can assume that the field intensity vanishes at zmax. Then, (12)where Bz = 0 is the magnetic field at the mid-plane and Bz = 0,z is its vertical component. The two terms respectively account for the magnetic pressure and the magnetic tension. This latter term describes the compression induced by the curvature of the field lines. However in our case, the field is quite turbulent and it is not clear that the net effect of the Bz component is captured by this simple relation. As it is small anyway (Fig. 3), we did not consider it further. Since the magnetic field evolves in a complex way, we simply used the values displayed in Fig. 3 to perform our estimate.

To compute the gravitational energy, let us define the weight of the gas disc with respect to a given gravitational potential φ: (13)We then write and , where φext is the potential given by Eq. (1). Using Eq. (6), we get (14)To be able to extract the scale height of the gas disc analytically, we approximated the external potential φext by φext(z) ≈ 2πGρextz2, where (15)describes the density of stars and dark matter in the mid-plane, and we supposed a Gaussian density profile. Then, (16)The equilibrium equation then becomes (17)This equation can easily be solved using a simple bisection method, once ηturb, nSN and B0 are specified.

To compare the estimates of H given by the above model and the simulations, it is useful to compute the half width at half maximum (HWHM) of the density distribution. Again assuming a Gaussian density profile, the HWHM is: (18)The results are summarized and compared to the simulations in Fig. 7, which shows the measured HWHM of the density profiles at 40, 60, and 80 Myr for the weak and standard feedback runs (respectively red and cyan diamonds) as well as the HWHM obtained by solving Eq. (17) taking the values of SN and B0 from Figs. 2 and 3. We note the values of the magnetic field of the x-axis correspond to the mean magnetic field in the equatorial plane at time 80 Myr for the standard feedback run (the same value is used for the three snapshots and the two models for more clarity). The value of ηturb is adjusted to match the hydrodynamical case (run B0L) at time 60 Myr and has been estimated to ηturb ≃ 0.3 (labelled model D1 in Fig. 7). We note that the results vary rather sensitively with it. The corresponding models show reasonable agreement except for the highest magnetization runs (B4) and at time 80 Myr in the hydrodynamical and weakly magnetized cases (B0 and B1). In an attempt to obtain a better fit at high magnetization, we have considered a model in which ηturb varies with B. Indeed in Iffrig & Hennebelle (2015a) we found that the momentum injected by a supernova in the dense gas, does increase with B by a factor on the order of ≃ 1.5. This is likely due to the correlation that magnetic field induces into the flow. When a magnetized fluid particle moves it entrains more fluid in its wake. The model labelled D2 (dashed line in Fig. 7) assumed where B5 is the magnetic intensity in units of 5 μG, that is to say the most magnetized run (B4L and B4W) has ηturb ≃ 1. As can be seen a better agreement is reached for the runs B4 at time 60 and 80 Myr (for run B2 the agreement is slightly better at 60 Myr and worse at time 80 Myr). The disagreement at time 40 Myr persists because the SFR is zero at 40 Myr for the most magnetized runs (B4). At this epoch the disc is probably contracting because of the flux expulsion and out of equilibrium.

Finally, we see that the weak feedback cases do not agree well with the model at time 80 Myr. This is likely because as clearly seen in Fig. 1 (bottom panel), the disc is undergoing a global infall. This is also consistent with the continuous increase of the SFR seen in the bottom panel of Fig. 2 (while on the contrary the standard feedback models display a nearly constant SFR). Therefore the simulations performed with weak feedback are out of equilibrium and the model is not expected to be valid.

thumbnail Fig. 7

Comparison between the half width at half maximum (HWHM) of the density profiles measured in the simulations (diamonds) and the analytical models at time 40, 60 and 80 Myr. Red lines represent the standard feedback case and cyan ones the weak feedback case. Solid lines stand for constant dissipation while dashed ones represent magnetically dependent dissipation (see text).

Open with DEXTER

We conclude that a simple hydrostatic model is reasonably successful in reproducing the disc thickness, though such a model is not very accurate. It suggests that the supernovae driving is not very efficient to inject energy in the system as the corresponding efficiency required is ηturb × ϵSN ≃ 0.015. We stress that this number reflects not only the energy injected by the supernovae into the turbulent gas but also the decay time of the turbulence. Our results are compatible with a strong field increasing this efficiency by a factor of the order of between two and three, though given the fact that the overall agreement between the model and the simulation is qualitative, this needs to be confirmed.

4. Turbulence properties

We now turn to the study of the detailed properties of the MHD supernova-driven turbulence which takes place in the stratified galactic disc.

4.1. Density distribution

We first consider the density distribution as it is a fundamental parameter for the ISM physics and the star formation process. Previous studies found that the density distribution associated to supersonic isothermal turbulence is log-normal (Vázquez-Semadeni 1994; Kritsuk et al. 2007; Federrath et al. 2008, 2010b), with a power-law tail due to gravity (Kritsuk et al. 2011b). Larger-scale simulations (e.g. Hill et al. 2012; Kim et al. 2013; Hennebelle & Iffrig 2014) find such a log-normal distribution only for the gas with density typically above 10 cm-3, which corresponds to the cold phase (the densest gas showing the power-law tail is not visible due to the lack of resolution).

The top panel of Fig. 8 displays the density PDF at time 40 Myr for the runs B0, B1 and B4. It shows a good agreement with a log-normal distribution for densities above 10 cm-3. Given the presence of sink particles, densities above 1000 cm-3 are not present in the simulations. The most diffuse gas is a consequence of the supernova feedback, and the dip at around 1 cm-3 is due to the thermal instability.

Since column densities are observationally inferred, though usually for individual clouds (e.g. Kainulainen et al. 2009; Schneider et al. 2015) rather than for a fraction of the galactic plane, we also present them here for future references and comparisons. The bottom panels show the column density distributions obtained by integrating face-on and edge-on through the galactic plane. The shape is also broadly lognormal. As can be seen, the most magnetized runs present the highest values, which is a consequence of the dense structures being larger.

4.2. Power spectra

In order to understand the details of the supernova-driven turbulence in the simulations, we take advantage of the uniform resolution to compute power spectra. Since the simulations include a stratified structure, we calculated both the full three-dimensional spectra and two-dimensional spectra on slices of constant altitude.

For the velocity, we also computed a Helmholtz decomposition between compressive (vanishing curl) and solenoidal (vanishing divergence) modes by projecting the Fourier transform of the velocity parallel and perpendicular to the wave vector: v = vcomp + vsol where and . For the two-dimensional power spectra, these quantities are computed on the whole 3D Fourier cube, and then an inverse Fourier transform is performed along the vertical axis. Then the power spectrum of the corresponding field is computed at a given altitude.

While power spectra of turbulent flows in the context of the ISM have been widely explored (see references below), it is important to study them in our simulations i) for a consistency check; ii) because there are specific aspects such as the stratification that have not been widely explored before.

4.2.1. Three-dimensional power spectra

Contrary to most studies of turbulence where the energy is injected at large scales, for instance by an external force at wave vectors 1 ≤ k< 3 (e.g. Kritsuk et al. 2007; Federrath & Klessen 2012; Federrath 2013), the turbulence in our simulation is supernova-driven, and energy is likely injected at scales on the order of 50 pc although since supernovae are correlated, they may also inject energy at larger scales. This work (Iffrig & Hennebelle 2015b, and this paper) shows similarities to the work by Padoan et al. (2016), but differs from it in the following respects: i) their 250 pc box does not include vertical stratification and is periodic; ii) in their study the supernovae are not correlated to star-forming regions; iii) their resolution is not uniform (1283 base grid). This will introduce some discrepancies that will be discussed. Given our resolution, the wave vectors k above 64kmin (where kmin = 2π/L, L = 1 kpc being the size of the simulation box) are affected by the numerical dissipation. Besides, the large-scale stratification will affect the three-dimensional power spectra at low k. Therefore, the power-law scalings are accurate only for (roughly) 10 <k/kmin< 64. The three-dimensional power spectra for various variables are shown in Fig. 9. Since the most magnetized runs have been performed at a resolution of 5123, we show here the runs B0L-B4L. For reference the high resolution simulation power spectra are given in Appendix A.

thumbnail Fig. 8

Density (upper panel) and column density (lower panels) distributions for the runs with strong feedback. The dotted parabolas are log-normal fits between 1 and 100 cm-3 and the dashed lines are tentative power-law fits above 100 cm-3, with slopes −1.1, −1.5, −1.6, and −1.5 for runs B0L, B1L, B2L, and B4L respectively.

Open with DEXTER

thumbnail Fig. 9

Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER

The velocity power spectra show reasonable agreement with the well-known models (Kolmogorov 1941; Iroshnikov 1964; Kraichnan 1965; Sridhar & Goldreich 1994; Lee et al. 2010; Grappin & Müller 2010; Mason et al. 2012; Beresnyak 2011), while not allowing discrimination between the two slopes k− 7/2 and k− 11/3. In the context of supersonic turbulence power spectra on the order of k-3.9 have been inferred (e.g. Kritsuk et al. 2007), which has been interpreted as a consequence of compressibility. In the present case, the somewhat shallower power spectrum may be interpreted as a consequence of the fact that most of the gas is warm and has Mach number on the order of one, meaning that the effect of compressibility could be less pronounced than in high Mach number flows.

The slopes of the density power spectra are very flat and compatible in the inertial range with E(k) ∝ k0, which is usually interpreted as a consequence of the steep variations due to supersonic shocks and thermal instability since the Dirac function has a power spectrum with an index equal to 0. Interestingly, the slope of the hydrodynamical runs (B0L) is slightly shallower than the MHD ones. The power spectra of log ρ present an index and a general behaviour close to the velocity power spectra (Schmidt et al. 2009; Audit & Hennebelle 2010).

The density-weighted velocity ρ1/3v, supposed to be the compressible equivalent of the velocity for the power spectra (see Kritsuk et al. 2007), has a slope that is even shallower than the velocity field power spectrum. At small k its value approaches k-2. At intermediate scales it becomes stiffer and approaches a slope on the order of k− 11/3 between k = 10 and 100, particularly for the magnetized runs, although the limited range does not allow a solid conclusion. This behaviour is a possible signature of the injection of turbulence at scale of about 100 pc. Another alternative explanation is that this signs the transition from 3D turbulence to a more 2D one since the scale height of the disc is about one tenth of the total box length. In this latter case, an energy power spectrum k−8/3 would be expected. Let us stress that the enstrophy cascade that would lead to k2k-3 requires conservation of the vorticity which is not ensured in the magnetized and non-barotropic flows like the ISM (Hennebelle & Audit 2007; Padoan et al. 2016).

thumbnail Fig. 10

Two-dimensional velocity power spectra with the Helmholtz decomposition. From left to right: runs B0L, B1L, B2L, and B4L. From top to bottom: altitude 0, 50, and 100 pc. The spectra are multiplied by k8/3 for comparison with the Kolmogorov scaling law.

Open with DEXTER

Finally we note that in the magnetized simulations, the magnetic field power spectra show the same behaviour as the velocity field ones, a behaviour found in turbulent periodic boxes (Kritsuk et al. 2011a).

4.2.2. Two-dimensional power spectra

The vertical stratification of our simulations induces strong contrasts as a function of altitude. Therefore, we computed two-dimensional power spectra on several horizontal planes to get a more detailed view of the consequences of this stratification. For most quantities like the density and the magnetic field, the power spectra do not look too different at different altitudes apart from their amplitude. For the sake of conciseness, we focus on the velocity power spectra, shown in Fig. 10. The two-dimensional power spectra of other quantities are shown in Appendix B.

Star formation simulations with turbulent forcing (Federrath & Klessen 2012) have shown that the star formation rate can be reduced by an order of magnitude between a purely compressive and purely solenoidal stirring. Energy equipartition between those two components would show up as a ratio Psol/Pcomp ≈ 2, since there are two solenoidal modes for one compressive.

In the galactic mid-plane, Fig. 10 shows that the compressive and solenoidal power spectra are comparable, which means that the compressive modes are a factor of approximately two higher than expected if energy equipartition was achieved. According to the results of Federrath & Klessen (2012), this means that star formation is neither in the most favorable regime (which would be purely compressive), nor in the least favorable (purely solenoidal). There is a possible trend that in the galactic plane at z = 0 the compressible mode power spectra are slightly steeper than the solenoidal one although the lack of statistics makes the noise level significant. This is broadly compatible with the work of Padoan et al. (2016) in which Psolk-3.31 and Pcompk-3.98. For the compressive modes, this corresponds to a Burgers’ (pressureless) fluid.

At higher altitudes, the solenoidal modes dominate, and the ratio of solenoidal to compressive power is stronger for higher magnetic field. This is expected since the magnetic field helps creating solenoidal motions and tends to impede the gas compression.

These results are different from the work of Padoan et al. (2016) since they find that the solenoidal modes dominate while we find that this depends on altitude but in the equatorial plane the compressible modes are dominating. One major difference with this work is the stratification induced by the galactic gravitational field. Another important one is the supernova scheme: they use a random supernova distribution, similar to scheme A of Hennebelle & Iffrig (2014), which injects supernovae mostly in the diffuse gas. The solenoidal power spectrum is slightly shallower than the Kolmogorov power spectrum, which is comparable to our results. It must be kept in mind that the exact correlation between supernova explosions and the dense gas is not well constrained. There is a possibility that the scheme used here is not accurate enough and that the supernovae should be less correlated with the dense gas. However in this case, the SFR is far too high as stressed in Hennebelle & Iffrig (2014). Generally speaking, the exact way feedback operates remains to be clarified. One very important constraint is that the SFR should be compatible with the observational rates.

thumbnail Fig. 11

Relative orientation of the velocity and the magnetic field as a function of time and magnetic field intensity for the strong feedback runs. The fields spontaneously align with time. The histograms are mass-weighted to give more importance to the gas in the galactic disc.

Open with DEXTER

thumbnail Fig. 12

Relative orientation of the velocity and the magnetic field as a function of time and magnetic field intensity for the weak feedback runs. The fields spontaneously align with time. The histograms are mass-weighted to give more importance to the gas in the galactic disc.

Open with DEXTER

4.3. Alignment of velocity and magnetic fields

Star formation simulations based on colliding flows (Inoue & Inutsuka 2009; 2012; Körtgen & Banerjee 2015) show the importance of the alignment between velocity and magnetic field: if the inflow velocity is not along the magnetic field lines, star formation is reduced efficiently. This is expected since the transverse component of the magnetic field is amplified by the converging velocity field. However, it is known that in magnetized flows the velocity and magnetic fields tend to align (Boldyrev 2006; Matthaeus et al. 2008; Banerjee et al. 2009). This effect is due in part to the Lorentz force which vanishes along the magnetic field and therefore is expected to be smaller when the magnetic and velocity fields are parallel. It is also a consequence of how the velocity and magnetic field get transported and generated by the flow.

To what extent colliding flow calculations with an inclined magnetic field are representative of the ISM remains to be clarified. To understand the exact role magnetic field is playing in the ISM, knowing its orientation with respect to the velocity field is crucial. For this purpose, we study the angle between velocity and magnetic fields, defined as (19)A uniform distribution of relative orientations would lead to a uniform distribution of this cosine.

The results are shown in Figs. 11 and 12 where we see that the mass weighted angle distribution clearly shows an excess in the aligned configuration. The amplitude of the effect increases with the magnetic intensity and the feedback strength. For the B1 run, the distribution is two times higher for cos(v,B) = 1 than cos(v,B) = 0, this value is approximately three to four for run B2 and 6–10 for run B4. For runs B1W, B2W and B4W, the effect is even more pronounced which demonstrates that feedback is playing an important role there. Stronger feedback tends to dealign the velocity and magnetic fields which is expected since in a super-Alfvénic shock, the transverse component of the magnetic field is amplified and the velocity tends to be perpendicular to it. Interestingly, we see in runs B1, B2 and B4 that the distribution at time 80 Myr shows less alignment that at time 40 Myr, which is also consistent with feedback reducing the alignment. This effect was reported in Passot et al. (1995). To get a more accurate understanding of the magnetic and velocity alignment we have also investigated the dependence of the alignment with gas density in Appendix C.

Altogether these results show that the magnetic field and the velocity field are well correlated and that it is necessary when estimating the magnetic field influence to consider this effect. In particular, even in the case of the less magnetized runs (B1) about half of the converging flows are expected to present an angle that is below 45°. This number is even higher for the more magnetized runs.

5. Structure formation

We now turn to the dense clouds, which form under the influence of gravity and turbulence. We are primarily interested here by obtaining a statistical description and to determine to what extent their properties vary with feedback and magnetization. In particular the existence of scaling relations, the so-called Larson relations, is well established (Larson 1981; Solomon et al. 1987; Falgarone et al. 2009; Heyer et al. 2009; Roman-Duval et al. 2010; Hennebelle & Falgarone 2012; Miville-Deschênes et al. 2017, etc.): where L, σ and M are respectively the size, velocity dispersion and mass of the clumps. Typical values of these parameters are σ0 = 1.1 km s-1, α = 0.5 (Falgarone et al. 2009), M0 = (228 ± 18)M and β = 2.36 ± 0.04 (Roman-Duval et al. 2010). We note that these observations are performed in the CO lines, which is typically tracing molecular gas of densities of the order of few 102 cm-3. Heyer et al. (2009), Miville-Deschênes et al. (2017) infer a relation that entails the cloud column density, namely: (22)where ΣR is expressed in M per pc. They show that this distribution presents less dispersion that the σR one and even less than the σM ones, suggesting that gravity could be playing an important role here. We note that the effect is clear for massive structures but less apparent for low mass ones.

The structures are obtained with a density threshold of 50 cm-3 using a simple friends-of-friends algorithm. The reason for this threshold is that at a density of 103 cm-3 the sink particles are being introduced. Therefore, to get a significant dynamical range, we adopt a value that is well below this threshold. This means that we may not be tracing exactly the same gas. However, observations of the atomic gas have also been performed in external galaxies as the LMC (Kim et al. 2007) and revealed similar behaviour (though different numbers are obviously inferred). In principle structures should be identified in the same way that observers proceed. However this would imply several steps and in particular the calculation of the CO molecules abundances (e.g. Duarte-Cabral et al. 2015). This latter point is particularly difficult because the CO abundances predicted by PDR codes for intermediate density gas (column densities smaller than a few 1020 cm-2) are underestimated by almost one order of magnitude (see Fig. 11 of Levrier et al. 2012). These issues would require a dedicated study and are clearly beyond the scope of the paper.

We define the velocity dispersion, σ, the cloud radius, the virial α parameter and the mass-to-flux over critical mass to flux ratio (Mouschovias & Spitzer 1976) as

(23)To compute Φ, the magnetic flux, we first compute the cloud center of mass, then we compute the flux accross the three planes parallel to xy, xz and yz and passing through the center of mass. We then take the largest of these three fluxes. To compute the radius, R, we use the eigenvalues, λi, of the inertia matrix, Iij defined by (24)While this choice is reasonable, it is not unique and we have tried different definitions such as using the largest eigenvalues of the mean spherical radius R = (V/ (4π/ 3))1/3 and the resulting distributions do not vary very significantly.

thumbnail Fig. 13

Clump scaling relations at 60 Myr. From left to right: runs B0, B1, B2, and B4. First row: mass-size relation. Second row: size-velocity dispersion relation. Top panels: strong feedback. Bottom panels: weak feedback. The solid red lines show the power-laws stated by Eq. (21).

Open with DEXTER

thumbnail Fig. 14

Virial α parameter and mass-to-flux over critical mass-to-flux ratio of clumps at 60 Myr. From left to right: runs B0, B1, B2, and B4. First row: mass-α relation. Second row: mass-μ relation. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER

thumbnail Fig. 15

Images of the four most massive clouds identified in the simulation B1 at time 60 Myr. Column density and mean velocity field in the z-plane.

Open with DEXTER

thumbnail Fig. 16

Clump mass spectra at 60 Myr. The lines represent the best power-law fit for M> 100 M. From left to right: runs B0, B1, B2, and B4. First row: strong feedback. Second row: weak feedback.

Open with DEXTER

The mass-size and size-velocity dispersion relations are shown in Fig. 13. The power-laws stated by Eq. (21) are show in the figure. As can be seen the mass-size relation of the structures follow a similar power-law to the observed one with β of the order of 2.3. The total mass is below the one inferred from CO survey but as explained above it is likely a consequence of the density threshold being too low. To verify this, we have extracted the clumps using a threshold of 200 cm-3. In this case the cloud masses is as expected about four times larger and present the same power-law behaviour. Interestingly, the number of small clumps is much higher in the hydrodynamical run B0 than in the MHD ones and decreases with magnetic intensity, while the power-law behaviour does not change. This effect, which has already been observed in smaller-scale simulations (Hennebelle 2013) is likely a consequence of magnetic tension, which makes the flow more coherent. The velocity dispersion is also displayed in Fig. 13. The values present a significant dispersion. The largest velocity dispersion of the clouds in simulations B0, B1 and B2 are comparable with the largest velocity dispersion inferred from observations (Falgarone et al. 2009). Both distributions present a large spread and some clouds have a velocity dispersion significantly below the mean value. For the sake of completeness, the velocity dispersion, σ as a function of RΣ is shown in Appendix D. As can be seen the velocity dispersion is significantly lower in the most magnetized case, simulation B4. This is a consequence of stronger field which makes it difficult to bend the field lines but also likely of the reduced star formation rate (as shown in Fig. 2). Interestingly the simulations with weak feedback B0W-B4W present very similar properties to the standard feedback case. This is indeed expected since, as discussed before, the amount of momentum delivered in the ISM are comparable because the SFR are higher in the weak feedback simulations.

Figure 14 displays the virial α parameter, which allows to quantify the importance of gravity in the cloud. The run B0 (hydrodynamical and standard feedback) presents a broad distribution of nearly three decades for 100 M clouds for which α goes from 100 to 0.1, though typical values are on the order of 10. The distribution is narrower for larger masses and for clouds of mass larger than 104M, it typically ranges between approximately one and ten (we note that exact value depends on the chosen definition of the radius, R). The run B1 shows similar trends except that the spread is significantly reduced for clouds of small masses. The same is true for runs B2L and B4L except that since they have been performed at a lower resolution, 100 M clouds are absent. The behaviour for the weak feedback runs is also similar with a trend for slightly lower values. Altogether these results suggest that the turbulence within molecular clouds is not primarily due to gravity for most low-mass clouds simply because the values of α can be larger than ten and are on average larger than the values for the most massive clouds. The relatively weak dispersion of α for the most massive ones on the other hand suggests that they are primarily driven by a combination of self-gravity and feedback. It is likely that for these objects, α tends to be self-regulated. First of all, we may expect a selection effect. While the most massive and unstable clouds are born out of an ensemble of structures, that on average are dominated by turbulence individually (as shown in Fig. 14), they would not have become strongly self-gravitating if α was too large. Second, there is likely an evolutionary effect induced by gravity which tends to produce collapse motions, that also have α on the order of a few although this may preferentially occur at scales not well described in the present simulations. Finally, when feedback becomes strong enough to unbind the clouds, that is to say when the diverging motions have α again equal to a few, the cloud is disrupted. These points are illustrated by Fig. 15 that displays the column density and the integrated velocity for four most massive clouds in simulations B1. Only the most massive one (top left panel) shows some clear sign of global infall and even there, it is obvious that there are plenty of non-infalling and disordered motions. For the three other clouds shown the most obvious trends are the diverging motions which are due to supernova feedback. This implies that in the present simulations, feedback processes, which are both spatially and temporally correlated with star formation, start to destroy the clouds before a global collapse takes place. Although this behaviour likely depends on the details of the feedback processes which are not sufficiently accurately described in this work, it must be reminded that the star formation efficiency is observed to be rather low in molecular clouds (Lada et al. 2010), which is barely compatible with global infall being dominant on large scales. At smaller scales (not well described in this work), the situation may be different (Peretto et al. 2007; Ballesteros-Paredes et al. 2011).

The mass-to-flux over critical mass-to-flux ratio, μ, is also displayed in Fig. 14. It increases from values of about 0.3 for 100 M clouds to about 8–10 for cloud masses of 104M and presents a rough scaling μM1/2. A similar relation has been obtained by Banerjee et al. (2009) and Inoue & Inutsuka (2012) where a relation μM0.4 has been inferred. As magnetic field is playing a significant role, it is important to understand the origin of this relation keeping in mind that getting the normalization factor (that is to say the value of μ at some specific value) is not straightforward, because, as discussed above, magnetic flux is getting expelled from the galactic disc. Let us consider that the gas has reached some statistical equilibrium constituted by a mixture of warm and cold dense gas, that it has a mean density ρISM and that it is threaded by a mean magnetic field, BISM. As a structure becomes assembled out of a radius R, its mass, M, is typically ρISMR3, while the flux is R2BISM. Therefore the mass-to-flux ratio, μ, is thus R. Since we get a mass size relation MR2.3, we get (25)which is in good agreement with the observed behaviour. This relation is important as it leads to a prediction of the field intensity in ISM structures. As we see it, this is essentially due to a simple geometrical effect, larger structures having a larger volume over surface ratio than smaller ones.

Finally Fig. 16 shows the mass spectrum of the clumps for all simulations. The shape observed in smaller scale simulations is recovered (e.g. Hennebelle & Audit 2007; Heitsch et al. 2008; Banerjee et al. 2009; Inoue & Inutsuka 2012; Padoan et al. 2016; Valdivia et al. 2016) with a plateau at small masses and a power-law at high masses. While the former is a consequence of numerical dissipation, the latter likely reflects the properties of turbulence as discussed in Hennebelle & Chabrier (2008). This good agreement between simulations performed at scales of 50 pc and the present ones which resolve the galactic disc is consistent with the idea that a large scale turbulent cascade is taking place and that the limited range of structure distribution, a clear consequence, of limited resolution, can be extrapolated to the regime of smaller structures.

Figure 16 confirms that stronger fields tend to diminish the number of small scale structures (see for example runs B0W, B1W and B2W which all have a resolution of 10243). Interestingly we also see that the power-law becomes shallower when the magnetic intensity increases going from about dN/ dMM-2 in the hydrodynamical simulations (run B0) to dN/ dMM-1.5 (run B4). This is also consistent with the idea that the fluid particle being partially linked by the field lines, they tend to form bigger clumps. Observationally a slope of about 1.7 has been inferred from CO survey (e.g. Kramer et al. 1998; Heithausen et al. 1998). Since run B1 presents an exponent close to this value, this is consistent, as the large scale magnetic field in this run is on the order of 3 μG and is therefore close to the mean galactic field.

6. Conclusion

We have performed a series of high resolution tridimensional numerical simulations with a resolution up to 10243, aiming to describe self-consistently the vertical structure of a galactic disc and a self-regulated star-forming ISM through supernova feedback. We considered four magnetizations and two feedback injections, one using canonical momentum injected by the supernovae and one four times below this value.

The measured SFR are comparable to the observational values, particularly with the standard feedback and magnetization. It is roughly four times larger when the weak feedback scheme is used. The hydrodynamical runs present SFR two times larger than the intermediate magnetization and the run with the strongest field presents SFR two or three times lower than in the intermediate field case. We found that while significant, the impact of the magnetic field tends to be limited by two effects. First of all magnetic flux tends to be expelled from the galactic plane probably because of the turbulent motions arising there. Second of all the magnetic and velocity fields are preferentially aligned reducing the effect of the Lorentz force. Comparison between an analytical model and the measured scale height, shows that indeed, except for the most magnetized runs, the magnetic field does not increase the disc scale height significantly. This allows us to also estimate the efficiency of the energy injection by the supernovae onto the gas within the galactic disc and we find it to be on the order of a few percents.

We computed tridimensional power spectra of various flow quantities such as density, magnetic field and velocity, finding classical behaviour although the slopes are closer to the canonical 11/3 than the 3.9 inferred for supersonic turbulence. As the simulations are strongly stratified, we also computed bidimensional power spectra in a series of horizontal planes at various heights. In particular, we performed a Helmholtz decomposition and found that in the equatorial plane, even for the strongly magnetized runs, the compressible modes tend to dominate the solenoidal ones. At higher heights the former becomes negligible. We stress that the dominance of the compressible modes in the galactic plane is possibly biased by our particular choice of supernovae driving.

Finally, we extracted the dense clouds and computed their physical properties, finding them to be reminiscent of the observed clouds though we do not exclude that their internal velocities may be too low, which may indicate that either feedback is not strong enough, either there is further energy injection from the large galactic scales. The mass-to-flux ratio is found to be M0.4−0.5 and a simple explanation has been proposed. The virial parameter, α, has been estimated and the shape of massα distribution is also similar to observations. At masses M ≃ 102−3Mα presents a large spread and is typically equal to 10–100. At masses M ≃ 104−5Mα is of the order of a few and presents a narrow distribution.

Acknowledgments

We thank the anonymous referee for a thorough and constructive report, which has improved the manuscript. This work was granted access to HPC resources of CINES and IDRIS under the allocation x2016047023 made by GENCI (Grand Equipement National de Calcul Intensif). P.H. acknowledge the finantial support of the Agence Nationale pour la Recherche through the COSMIS project. This research has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013 Grant Agreement No. 306483).

References

Appendix A: Three-dimensional power spectra of high resolution runs

Here we provide for comparison and reference the power spectra of the high resolution runs corresponding to the standard feedback (Fig. A.1) and to the weak feedback (Fig. A.2). As can be seen these spectra are very similar to the ones presented in Fig. 9, which shows that reasonable convergence has been reached but also that the weak feedback does not affect too much the flow properties.

thumbnail Fig. A.1

Runs B0 and B1 (high resolution with standard feedback). Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER

thumbnail Fig. A.2

Runs with weak feedback (B0W, B1W, B2W and B4LW). Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER

Appendix B: Two-dimensional power spectra of density and magnetic field

thumbnail Fig. B.1

Two-dimensional velocity power spectra of ρ, B and log ρ. From left to right: runs B0L, B1L, B2L, and B4L. From top to bottom: altitude 0, 50, and 100 pc. The spectra are multiplied by k8/3 for comparison with the Kolmogorov scaling law.

Open with DEXTER

As the simulations presented here have a strong stratification, we show for the sake of completeness a series of two-dimensional power spectra obtained at three altitudes. As can be seen the index of the power spectra are broadly compatible to the three-dimensional ones presented in Fig. 9 with some noticeable difference. In particular the index of the magnetic field power spectra varies with altitude.

Appendix C: Alignment of velocity and magnetic fields for density range

Figure C.1 displays the relative orientation of the velocity and magnetic fields as a function of time and for five density bins allowing a more detailed analysis than Fig. 12. As can be seen in the very diffuse gas (n< 0.1 cm-3), there is no alignment. Clearly this is because this gas is produced by supernova explosions. For denser gas, the relative orientation distribution is nearly the same for the three bins 0.1–1, 1–10 and 10–100 cm-3. As expected the alignment is stronger when the field intensity is higher. Interestingly, there is a clear trend for the alignment being less pronounced for n> 100 cm-3. This is consistent with the contraction occurring mainly along field lines at low and intermediate densities and becoming less focused at higher densities, either because gravity leads to global contraction or because the Alfvénic Mach numbers tends to be higher at higher densities.

thumbnail Fig. C.1

Relative orientation of the velocity and the magnetic fields as a function of time for five density bins.

Open with DEXTER

Appendix D: Velocity as a function of RΣ

Figure D.1 shows the velocity dispersion as a function of RΣ for comparisons with Miville-Deschênes et al. (2017). The numbers inferred are pretty similar. However, the correlation is not visibly better than the σR one displayed in Fig. 13.

thumbnail Fig. D.1

Clump scaling relations at 60 Myr. From left to right: runs B0, B1, B2, and B4, size times column density-velocity dispersion relation. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER

All Tables

Table 1

Summary of the runs.

All Figures

thumbnail Fig. 1

Column density maps for the various runs. From left to right: runs B0, B1, B2 and B4 (increasing initial magnetic field). From top to bottom: fiducial pSN = 4 × 1043 g cm/s at 40 and 80 Myr, reduced pSN = 1043 g cm/s at 80 Myr. First row: edge-on view. Second row: face-on view.

Open with DEXTER
In the text
thumbnail Fig. 2

Mass in sinks (first row) and estimated star formation rate (second row) as a function of time. The dotted line for the sink mass corresponds to the initial gas mass in the simulation box. The dotted line for the star formation rate corresponds to the typical observational value (e.g. Kennicutt et al. 2007) given our column density. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER
In the text
thumbnail Fig. 3

Magnetic field in the equatorial plane as a function of time for the three runs B1L, B2L and B4L. Solid line is the x-component while dotted and dot-dashed are the y and z-components respectively. As can be seen the x-component decreases with time which indicates that the initial magnetic field gets expelled from the galactic plane.

Open with DEXTER
In the text
thumbnail Fig. 4

Vertical density profiles. Top: disc profile for the different runs with strong feedback, at 40 Myr. Bottom: evolution with time in the B1 and B1L runs. The dashed lines are the lower-resolution profiles.

Open with DEXTER
In the text
thumbnail Fig. 5

Mean height as a function of gas density for runs B0L, B1L, B2L and B4L at time 60 Myr.

Open with DEXTER
In the text
thumbnail Fig. 6

Vertical kinetic, thermal and magnetic pressure profiles 25 Myr after the beginning of star formation. From top to bottom: runs B0L, B1L, B2L and B4L.

Open with DEXTER
In the text
thumbnail Fig. 7

Comparison between the half width at half maximum (HWHM) of the density profiles measured in the simulations (diamonds) and the analytical models at time 40, 60 and 80 Myr. Red lines represent the standard feedback case and cyan ones the weak feedback case. Solid lines stand for constant dissipation while dashed ones represent magnetically dependent dissipation (see text).

Open with DEXTER
In the text
thumbnail Fig. 8

Density (upper panel) and column density (lower panels) distributions for the runs with strong feedback. The dotted parabolas are log-normal fits between 1 and 100 cm-3 and the dashed lines are tentative power-law fits above 100 cm-3, with slopes −1.1, −1.5, −1.6, and −1.5 for runs B0L, B1L, B2L, and B4L respectively.

Open with DEXTER
In the text
thumbnail Fig. 9

Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER
In the text
thumbnail Fig. 10

Two-dimensional velocity power spectra with the Helmholtz decomposition. From left to right: runs B0L, B1L, B2L, and B4L. From top to bottom: altitude 0, 50, and 100 pc. The spectra are multiplied by k8/3 for comparison with the Kolmogorov scaling law.

Open with DEXTER
In the text
thumbnail Fig. 11

Relative orientation of the velocity and the magnetic field as a function of time and magnetic field intensity for the strong feedback runs. The fields spontaneously align with time. The histograms are mass-weighted to give more importance to the gas in the galactic disc.

Open with DEXTER
In the text
thumbnail Fig. 12

Relative orientation of the velocity and the magnetic field as a function of time and magnetic field intensity for the weak feedback runs. The fields spontaneously align with time. The histograms are mass-weighted to give more importance to the gas in the galactic disc.

Open with DEXTER
In the text
thumbnail Fig. 13

Clump scaling relations at 60 Myr. From left to right: runs B0, B1, B2, and B4. First row: mass-size relation. Second row: size-velocity dispersion relation. Top panels: strong feedback. Bottom panels: weak feedback. The solid red lines show the power-laws stated by Eq. (21).

Open with DEXTER
In the text
thumbnail Fig. 14

Virial α parameter and mass-to-flux over critical mass-to-flux ratio of clumps at 60 Myr. From left to right: runs B0, B1, B2, and B4. First row: mass-α relation. Second row: mass-μ relation. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER
In the text
thumbnail Fig. 15

Images of the four most massive clouds identified in the simulation B1 at time 60 Myr. Column density and mean velocity field in the z-plane.

Open with DEXTER
In the text
thumbnail Fig. 16

Clump mass spectra at 60 Myr. The lines represent the best power-law fit for M> 100 M. From left to right: runs B0, B1, B2, and B4. First row: strong feedback. Second row: weak feedback.

Open with DEXTER
In the text
thumbnail Fig. A.1

Runs B0 and B1 (high resolution with standard feedback). Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER
In the text
thumbnail Fig. A.2

Runs with weak feedback (B0W, B1W, B2W and B4LW). Three-dimensional power spectra. From top to bottom: density, velocity, density-weighted velocity, magnetic field. The spectra are multiplied by k2 (such that the Kolmogorov scaling corresponds to a slope of −5/3).

Open with DEXTER
In the text
thumbnail Fig. B.1

Two-dimensional velocity power spectra of ρ, B and log ρ. From left to right: runs B0L, B1L, B2L, and B4L. From top to bottom: altitude 0, 50, and 100 pc. The spectra are multiplied by k8/3 for comparison with the Kolmogorov scaling law.

Open with DEXTER
In the text
thumbnail Fig. C.1

Relative orientation of the velocity and the magnetic fields as a function of time for five density bins.

Open with DEXTER
In the text
thumbnail Fig. D.1

Clump scaling relations at 60 Myr. From left to right: runs B0, B1, B2, and B4, size times column density-velocity dispersion relation. Top panels: strong feedback. Bottom panels: weak feedback.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.