© ESO, 2014

1. Introduction

Star formation is a multi-scale and multi-physics problem, which is only partially understood. In particular, what physical process is responsible for the relatively low star formation of the Milky-way (e.g., Zuckerman & Evans 1974; Dobbs et al. 2013) remains a subject of controversy.

Historically, three main physical processes have been emphasized, namely magnetic field (e.g., Shu et al. 1987), turbulence (e.g., Mac Low & Klessen 2004), and stellar feedback (e.g., Mac Low 2013; Agertz et al. 2013), which include supernovae explosions, ionising radiation, heating by stellar radiation, stellar outflows, and stellar winds. While magnetic field has been measured to have substantial intensities (Crutcher 2012), it may be nevertheless too weak to reduce the star formation rate (SFR) by orders of magnitude. The effects of turbulence and feedback are not straighforward to discriminate. In particular, turbulence decays in a few crossing times (e.g., Mac Low & Klessen 2004) and must be fed at large scales either through galactic large scale spiral waves or by the various sources of stellar feedback. The feedback may therefore play a dual role in limiting the amount of mass that is eventually accreted by stars, while at the same time triggering large scale turbulence. In any case, previous studies, which have been simulating a whole galaxy found a clear impact of the feedback onto the SFR (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), which is able to reduce star formation substantially, may be up to observed values. Although there is a general agreement that the simulations without feedback present an SFR that is too high, the amount by which it is reduced when feedback is introduced depends on which feedback is introduced and how. For example, Tasker & Bryan (2006), who perfomed simulations with supernovae feedback found that the SFR is reduced by a factor of about 2. Tasker (2011) included the UV radiation feedback found that the SFR is typically reduced by a factor of 1.5−2. Finally, Hopkins et al. (2011) have introduced the radiative feedback assuming that the radiation of stars can efficiently communicate its momentum to the gas. They found that the SFR can be reduced by a factor on the order of 10 to 30.

In spite of these studies, the exact roles played by feedback, for both triggering the turbulence and for limiting the star formation are only partially understood. This is also the case for the role played by magnetic field on the various steps of the star formation process. First of all, given the limited resolution of large scale studies (typically a few pc), the exact way feedback is applied remains partially arbitrary. In particular, the first pioneering studies, which have first focussed on kpc scales regulated by supernovae explosions (de Avillez & Breitschwerdt 2004, 2005; Joung & Mac Low 2006) did not include self-gravity and, therefore, could not associate supernovae with star formation events accurately, though they unanbigously show the relevance of these studies. Even when self-gravity is included, the exact influence of the choices made to determine their locations has not been explored clearly. Second of all, the influence that the magnetic field has on the SFR at the kpc scale is less explored. The only simulations that include both self-gravity and magnetic field, which have been performed to date are presented in Wang & Abel (2009), Pakmor & Springel (2013), Beck et al. (2013). These authors have concluded that the magnetic field reduces the SFR by a factor of a few. Moreover these studies modelled a whole galaxy implying that the spatial resolution is limited in describing the ISM structure.

In parallel to the numerical studies (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), a few analytical models have been developed and compared with observations and simulations (Ostriker et al. 2010; Kim et al. 2011; Faucher-Giguère et al. 2013). For example, in their model, Ostriker et al. (2010) consider vertical mechanical equilibrium between gravity and pressure (mainly thermal and kinetic) as well as thermal equilibrium in the ISM (i.e. equilibrium between heating from stars and cooling). It is then further assumed that the thermal and turbulent supports are proportional. This leads them to predict the SFR as a function of the column density through the galactic planes, and the model compares well with a sample of observations (Leroy et al. 2008) and simulations (Kim et al. 2011). In their models, Faucher-Giguère et al. (2013) consider a galactic disc, which is also in equilibrium along the vertical direction but assumes that the galactic disc has a Toomre parameter (Toomre 1964) of about Q ≃ 1, that is to say a disc which is in marginal equilibrium and self-regulates. They then perform an energy budget between the energy dissipated by units of time through turbulent cascade and the energy injected through stellar feedback. In both cases, important assumptions are made regarding how energy and momentum are injected within the dense gas. While it is clearly unavoidable to make such assumptions, given the difficulty of the problem, it is nevertheless important to understand how exactly are momentum and energy injected, more precisely how they distribute between diffuse and dense gas. More generally, what are the uncertainties induced by our incomplete understanding of the correlation between massive stars at the origin of most of the feedback, and the surrounding dense material ? It is the purpose of the present paper to address these issues.

In this paper, we adopt a similar setup to the one adopted by 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), that is to say a kpc simulation of a galactic disc in which turbulence is driven by supernovae remnants. Previous studies have been finding that they can reproduce many of the interstellar medium feature, such as a multi-phase ISM, approximate energy equipartition among thermal, magnetic and kinetic energies as well as galactic outflows and the formation of molecular clouds, which therefore demonstates the interest of performing this type of simulations. Indeed, this range of scales is a good compromise between the need for enough numerical resolution to describe sufficiently well the ISM physics and the amount of computational power available on present computers. While the results obtained previously are encouraging, only a few of these works have treated self-gravity, and none of these works have considered magnetic field and self-gravity, which is mandatory for a proper physical description.

We include both self-gravity and magnetic field, and we explore various schemes for the supernovae feedback going from a random distribution to a distribution in which supernovae are correlated both spatially and temporally with star formation. This spatial and temporal correlation turns out to be drastically important when self-gravity is self-consistently treated. The primary reason is that, a dense core becomes largely decoupled from the surrounding medium when it undergoes gravitationnal collapse and therefore little influenced by the supernovae, which may be exploding nearby. This occurs only if a supernova explodes within the collapsing region and while collapse is occuring feedback has a significant impact and can reduce the mass that is eventually accreted.

The second section of the paper presents the numerical setup and, in particular, discusses the various schemes we have been developing to implement supernovae feedback. In the third section, we describe the disc structure in the various models, while in the fourth section we investigate the properties of the multi-phase ISM in two of the simulations. In Sect. 5, we discuss the SFR and the mass distribution of the sink particles formed in the different simulations. Section 6 concludes the paper.

2. Numerical setup

Fig. 1 Column density along the z-axis (left panels) and along the y-axis (right panels) for MHD run C1 (upper panels) and hyrodynamical run C2 (lower panels). The left panels illustrate the complex multi-phase structure of the galactic plane, while the right panels show the stratification induced by the gravitational field of the galaxy. In the hydrodynamical run, the interstellar medium is more fragmented than in the MHD one.

Fig. 2 Column density along the z-axis (left panels) and along the y-axis (right panels) for the run with no feedback, NF1, (upper panels) and the run with scheme D (lower panels). Comparison with Fig. 1 clearly shows the drastic impact of the feedback. In particular, the run NF1 shows thin and long filaments in which self-gravitating fragments develop, indicating that the gas is primarily organised by self-gravity and undergoes run away collapse. The galactic disc is considerably thinner than when supernovae are included. On the other hand, run D shows a broader galactic disc than in runs C1 and C2, illustrating the importance of the spatial correlation between supernovae explosions and dense gas.

2.1. Physical processes and initial conditions

The physical processes and initial conditions are similar to what has been described by previous authors (e.g., de Avillez & Breitschwerdt 2004, 2005; Joung & Mac Low 2006; Hill et al. 2012; Kim et al. 2013). We consider a 1 kpc computational box in which a gravitational field is applied along the z-axis. Its value is identical to the choice made by Joung & Mac Low (2006) taken from Kuijken & Gilmore (1989) and is given by $\begin{array}{ccc}\mathit{g}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{1}}\mathit{z}}{\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}\mathrm{-}{\mathit{a}}_{\mathrm{2}}\mathit{z,}& & \end{array}$(1)where a₁ = 1.42 × 10^-3 kpc Myr^-2, a₂ = 5.49 × 10^-4 Myr^-2, and z₀ = 0.18 kpc. This gravitational field represents the contribution of the stars and dark matter in our Galaxy.

We solve the ideal magneto-hydrodynamical (MHD) equations in the presence of self and external gravity and include cooling and heating processes relevant for the ISM. The equations are given by $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \mathit{\rho }\left[\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{v}\mathrm{·}{\nabla }\mathrm{)}{v}\right]& \mathrm{=}& \mathrm{-}{\nabla }\mathit{P}\mathrm{+}\frac{\mathrm{(}{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}}{\mathrm{4}\mathit{\pi }}\\ & & \\ \mathit{\rho }\left[\frac{\mathit{\partial }{e}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{v}\mathrm{·}{\nabla }\mathrm{)}{e}\right]& \mathrm{=}& \\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& \\ \mathrm{\Delta }\mathit{\phi }& \mathrm{=}& \end{array}$where all notation have their usual meaning. The heating and cooling terms, which appear in the loss function ℒ, are identical to the one considered in Audit & Hennebelle (2005), which include UV heating due to photoelectric effect on grains, Lyman-α, oxygen, ionised carbon cooling and cooling due to the recombination onto grains (see e.g., Wolfire et al. 2003). At this stage we use a constant UV field, which is not correlated to the SFR (e.g., Tasker 2011; Kim et al. 2013). At temperature larger than 10⁴ K, we use the fit provided in Joung & Mac Low (2006) for the cooling function inferred by Sutherland & Dopita (1993). Coriolis and centrifugal forces are not included at this stage.

At the beginning of the simulation, the density distribution is given by$\begin{array}{ccc}\mathit{n}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{0}}\exp \left(\mathrm{-}{\left(\frac{\mathit{z}}{{\mathit{z}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\right)\mathit{,}& & \end{array}$(7)with n₀ = 1.5 cm^-3 and z₀ = 150 pc. The temperature is initially equal to about 8000 K, which corresponds to the temperature of the warm neutral gas (WNM). A turbulent velocity field is generated using random phase and a Kolmogorov powerspectrum. Its rms amplitude is equal to 5 km s^-1. Finally, the magnetic field is initially aligned along the x-axis and is proportional to the density field. Its value in the equatorial plane is about 2.5 μG for our fiducial run (later named run C1), the value of 0.5 μG is also explored.

2.2. Supernovae prescriptions

To take the supernovae feedback into account, we first select a position as described below. We then increase the thermal energy of all the cells located at a distance smaller than three grid cells from the supernova centre in such a way that the thermal energy is uniform in this sphere and equal to 10⁵¹ erg. In most, but not all, of our runs we have also introduced a kinetic feedback. This is achieved by adding to the corresponding cells a radial homologous velocity field (proportional to the distance from the supernova centre). The total kinetic energy is equal to 5% of the thermal energy, which corresponds to the typical momentum that is injected at the end of the Sedov phase, that is to say, during the phase for which the supernova expansion remains nearly adiabatic. Indeed, when the shell surrounding the supernova bubble becomes radiative, most of the energy is radiated away, and the expansion proceeds at constant momentum (e.g., Chevalier 1977). Fortunately enough, this momentum is largely independent of the density field (Blondin et al. 1998) even when it is highly irregular (Iffrig & Hennebelle, in prep.). In principle note that, the momentum should be self-consistently generated during the Sedov expansion. However, the lack of resolution does not guarantee that this phase is well treated when the supernovae explodes in a dense region. In any case, as described below, we have also performed the case without kinetic feedback for comparison. These two runs probably constitute upper and lower limits.

The spatial location of the supernovae is another important aspect. Previous authors (e.g., de Avillez & Breitschwerdt 2005, 2007; Joung & Mac Low 2006) have been distributing them randomly or in correlation with density. These authors also tested the case where the supernovae are clustered (see e.g., de Avillez & Breitschwerdt 2004; Joung & Mac Low 2006,for a description). In this work, four spatial and temporal supernovae distributions have been tested, hereafter scheme A, B, C, and D. We recall that an important difference with these classical studies is that self-gravity is treated.

Scheme A is very similar to the scheme described, for example, in Joung & Mac Low (2006). The supernovae are distributed randomly in the x and y-directions. To mimic the observed supernovae distribution in the Milky Way, their z-coordinate follows a Gaussian distribution of thickness equal to 150 pc. Their rate is equal to the observed galactic rate and is equal to 1/50 per year. One difference is that we use a fixed radius for the supernovae remnant rather than using a radius which enclosed a fixed mass. This implies that the timesteps in the simulation can be quite low, since the temperature is higher when a supernova explodes in a diffuse medium. Another difference is that we do not redistribute the mass within the supernova radius as Joung & Mac Low (2006), who impose a uniform gas density inside the supernova bubbles. In principle, it could be worth testing all these choices but here we focus on a different issue.

Scheme B consists looking at the density maximum in the simulations, and choosing its position for the supernova centre. The rate is also equal to the galactic one. A supernovae is introduced only if the density peak is larger than 10 cm^-3 in the simulation. However, denser gas develops rapidly in the simulation and except for the very first, the supernovae are generally associated to gas of densities 10² − 10³ cm^-3, which is present at all time. This scheme has the advantage to have a good spatial correlation with star formation events. However, it does not have any temporal correlation, since the supernovae rate is fixed and equal to the galactic one.

Schemes C and D are different and take advantage of the sink particles used in our simulations. Each time a sink particle has accreted 120 solar masses, we place a supernova in its neighborhood. This prescription is motivated by the typical abundance of stars more massive than the eight solar masses needed to give raise to a supernova explosion. However, we do not place the supernova centre directly at the sink position for various reasons. First of all, as described later, the sink particles have a radius of four computational cells, which corresponds to 16 pc with our current resolution. This number of cells for the sink radius is typical to what is usually assumed (e.g., Krumholz et al. 2004). By definition, it sets the limit of the resolution in the simulation and inside the sink, the gas distribution is not well described. Second of all, it takes at least 4 Myr for the most massive stars to explode. During this time both the star and the cloud have evolved. For example in 10 Myr, a star, which moves at a velocity of 1 km s^-1, will have to cover a distance of about 10 pc. Finally, because of the other sources of feedback, massive stars may have pushed the dense gas away before supernovae explode. To test the importance of the spatial correlation between supernovae and sink particles, we have implemented two prescriptions. First (scheme C), we randomly place the supernovae within the sink particles radius. Second (Scheme D), we place them within a shell of inner radius equal to the sink radius and outer radius equal to two times this value. As seen later, these two schemes lead to similar but not identical results.

Another important issue is the time delay that should be taken into account, since supernovae typically explode between about 4 and 40 Myr after the formation of the massive star. Although we note that some delay is introduced with our scheme since at least 120 solar masses of gas have to be accreted before supernovae take place, it is in most of the time shorter (10⁵ − 10⁶ yr) than a few Myr. However, as emphasized in other studies (e.g., Matzner 2002; Dale et al. 2012, 2013; Agertz et al. 2013; Kim et al. 2013), other sources of feedback, namely ionising radiation, winds, and jets, start influencing the surrounding clouds much earlier. These sources should be taken into account as well. Since we feel it is important to go step by step, we postpone studying these effects and concentrate for now on supernovae only.

2.3. Numerical code and resolution

To carry out our numerical simulations, we employed RAMSES (Teyssier 2002; Fromang et al. 2006), an adaptive mesh refinement (AMR) code that uses Godunov schemes to solve the MHD equations and the constrained transport method to ensure that divB is maintained at zero within machine accuracy. For most runs, we do not use the AMR capacity and keep the resolution fix using 256 grid points in each direction. However, we also perform a run with one more AMR level introduced when the density reaches a value of 10 cm^-3. The computational box size is equal to 1 kpc, and the spatial resolution is 4 pc. This choice is dictated by the very large number of time steps (≃50 000−100 000), which are required for these calculations. This is because supernovae feedback introduced very high velocities of the order of a few 100 km s^-1 and temperatures as high as 10⁸ K in a few cells. Moreover, we integrate far enough to make sure that some quasi-stationary regime has been reached. This is assessed by verifying that the mean profile of various quantities, such as densities, does not vary significantly with time. It should be stressed, however, that no strict stationarity can be reached since self-gravity is treated and accretion is occuring onto the sink particles.

We use periodic boundary conditions in the x and y directions and and outflow condition at the z boundaries. In particular, this implies that the gas ejected from the galactic disc can escape the computational box.

In this work, sink particles (implemented in the public version of RAMSES) are being used to follow the dense regions, which have collapsed under the influence of self-gravity. They closely follow the implementation of Krumholz et al. (2004). The sinks are introduced when a density threshold of 10³ cm^-3 is reached. Their radius is equal to four grid cells. A new sink can be created only if it is not located closer than ten grid cells from another sink. When the sinks get too close, i.e. closer than one grid cell, they get merge using a friend of friend procedure. Finally, the sinks accrete gas from surrounding cells, if they are located at less than a sink radius and if the density is larger than 10³ cm^-3. The sink particles interact with the surrounding medium through the gravitational field. The contribution of the sink to the gravitational potential is included using a particle-mesh approach, which implies that the mass within the sink is projected onto the grid and added to the gas density when the Poisson equation is solved.

2.4. Comparison between various setups

It is worth emphazising the differences between the various setups, which have been used so far. Keeping in mind that given the complexity of the problem under investigation, i) it is hard to include every relevant process; and ii) it is important to perform studies, which make different choices and approximations to disentangle the effects of the different physical processes.

In the setup used by de Avillez & Breitschwerdt (2005); Joung & Mac Low (2006) and Hill et al. (2012), hydro or MHD equations are solved using a 500 pc or a 1 kpc size box in the equatorial plane but a much larger scale height (typically 5 to 10 times these values). This insures a good description of the galactic wind, and the disc-halo connection though large scale modes may be filtered out by the elongated box. Gent et al. (2013) proceed somewhat similarly but also include the galactic shear in their study. In these studies, the supernovae explode randomly or are correlated with density peaks. None of these studies include self-gravity.

Kim et al. (2011, 2013) do include self-gravity and, like us, consider a cubic computational box. These choices are clearly putting more emphasis on the disc itself and, therefore, on the star formation, which takes place, than on the galactic outflows and the halo. There is an important difference regarding the forcing. They estimate the number of massive stars that should form in a given region rather than following the accretion onto sink particles, and they consider only the mechanical feedback from supernovae.

As seen below, the different prescriptions lead to quite different results. In particular, the correlation between star forming gas and the feedback is a necessary condition to prevent very efficient star formation. This is however mostly true if self-gravity is included.

2.5. Runs performed

In the present paper, we perform various runs to test the influence that the magnetic field has onto the galactic disc evolution and to study the influence of the various prescriptions for the supernovae feedback. In our fiducial run (later referred as C1), the magnetic intensity has an intensity in the equatorial plane of about 2.5 μG, and scheme C is used for the supernovae with both thermal and kinetic feedback.

To study the influence of the magnetic field, we perform an hydrodynamical run (later referred as C2) and a run with a lower magnetisation initially equal to about 0.5 μG in the midplane (run C3). Apart from the strength of the field, these runs are identical to run C1.

To study the influence of the feedback scheme, we perform a series of runs identical to run C1 apart for the feedback scheme. First, we run two cases without any feedback, one purely hydrodynamical (NF2) and one with 2.5 μG initially (NF1). Second, we consider two cases with scheme A and B, respectively (simply labelled runs A and B), and both thermal and kinetic feedback. Third, we perform a simulation with supernovae scheme C but with thermal feedback only (run C4). Fourth, we carry out a calculation with scheme D (run D).

Finally, to investigate the important issue of numerical resolution, we also present a run identical to run C1 but use another level of refinement, which leads to an effective resolution of about 2 pc. The refinement is performed when the cell density reaches a threshold of 10 cm^-3 leading to a total number cells, which is comparable to the number of cells in the non-refined runs at this level. The results are presented in the appendix.

Table 1 summarizes the runs performed in the paper and provide consistent labels.

We show four time steps for each run. The first one is at about 25−30 Myr, the second at 40 Myr and the third at about 55−60 Myr. For the fourth, we select a time step of about 90 Myr for runs C1 and C2 and about 65−70 Myr for runs A, D, and NF1. The first time step, have been chosen at the beginning of the star forming phase (see Fig. 13) and the second at about 10 Myr later because the SFR is typically close to its maximum. For the third time step, the SFR is nearly constant with time, and, thus, the profiles correspond to the quasi-stationnary regime. This is confirmed by the last time step which show no significant evolution with respect to the third one although for runs A, D and NF1 the evolution is faster because accretion is higher and the profiles keep evolving rapidly at later times.

For the sake of conciseness, we will select the runs which we think are most relevant to emphasize the impact of the physics and of the schemes. When investigating the global disc structure (Sect. 3), we concentrate mainly on runs NF1, A, C1, C2 and D. When we investigate the multiphase ISM (Sect. 4), we restrict to runs C1 and C2. The other runs (NF2, C3 and C4) are used to quantify the influence of the supernovae scheme on the SFR (Sect. 5).

3. Global structure

In this part, we characterize the global structure of the galactic discs.

3.1. Qualitative description

Figure 1 shows the column density field along the z-axis (left panels) and y-axis (right panels) for MHD run C1 (upper panels) and hydrodynamical run C2 (lower panels). In both cases, the disc is clearly visible, although its structure is quite irregular and varies significantly from place to place. The disc is slightly thicker in the magnetized run than in the hydrodynamical one, which is a clear consequence of the magnetic support. The column density distributions appear different, with and without magnetic field. In particular, small scale fluctuations are more pronounced in the hydrodynamical run. This trend is also similar to what has been found at smaller scales (e.g., Hennebelle et al. 2008) and has been interpreted by Hennebelle (2013) as a consequence of the Kelvin-Helmholtz instability, being stabilized by the magnetic field (Ryu et al. 2000). As discussed later in the paper, this has consequences on the mass of self-gravitating objects, that form.

Figure 2 shows column densities for the run without supernovae feedback, run NF1 (upper panels), and with our feedback scheme D as described above. The disc structure is very different in both cases. When no feedback is applied, very long filaments develop across the computational box. They converge towards a region, where most of the mass accumulates. There are much less dense regions compared to run C1, and the disc is obviously thinner. This behaviour is very similar to what is reported in Hopkins et al. (2011) (see for example their Fig. 2). On the other hand, scheme D, leads to a disc whose structure is much more irregular than the structure of the disc obtained with scheme C, although the number of supernovae and their energy are identical in the two cases. This constitutes a clear confirmation to the works of previous authors (e.g., Bournaud et al. 2010; Dobbs et al. 2011; Hopkins et al. 2011) that feedback is playing a crucial role for the structure of galactic discs and (see next section) for regulating star formation. It is also clear that the correlation between the gas and the supernovae does influence significantly the galactic disc structure.

3.2. Disc density profile

One fundamental aspect for the galactic structure is the density profile and the typical thickness of the gas distribution. In the case of the Milky Way, it has been measured (e.g., Ferrière 2001) that the different phases have different distributions. They tend to be roughly Gaussian, but their thicknesses vary. The molecular gas has a full width at half maximum (FWHM) of about 120 pc, it is about 230 pc while for the atomic gas. Both the molecular and the atomic gas have a mean density of about 0.5 cm^-3 for a total of 1 cm^-3.

Figure 3 shows the density distribution along the z-axis for five different models. The run C1 presents a density of about 3−4 cm^-3 and a FWHM of about 50 pc. This is also roughly the case for the hydrodynamical run C2. Comparison with run C1b shown in the Appendix reveals that numerical resolution may be partly at the origin of this discrepency. In the run based on scheme D, the maximum density varies with time because of stronger accretion but is around 1−2 cm^-3 after 40 Myr. The FWHM is larger than in the previous cases and equal to about 100 pc. When no feedback is included, the disc tends to be thinner. For example, at times 55.8 and 64.8 Myr, the thickness is about 30 pc. A similar distribution is found at time 57.1 and 65.4 Myr with scheme A, which, as seen later, has an accretion behaviour that is very similar to the run without feedback.

3.3. Pressure support

Since the disc thickness is a direct consequence of the various supports, we present the profile along z-axis of the three relevant pressures, namely thermal, kinetic, and magnetic ones, noted P_th, P_kin and P_mag, respectively. While P_th and P_mag have their standard definitions, P_kin is taken as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{kin}}\mathrm{=}\frac{\sum {\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathit{\rho }\mathrm{d}\mathit{V}}{\sum \mathrm{d}\mathit{V}}\mathrm{·}& & \end{array}$(8)Figures 4−6 show P_th, P_kin, P_mag respectively, as a function of altitude.

Fig. 3 Mean density profile along the z-axis for five different models (see label) at four different time steps. The disc profile is much thinner for runs NF1 and A than for runs C1, C2, and D.

Fig. 4 Thermal pressure profile along the z-axis for five different models (see label) at four different timesteps. The largest values are obtained for run D and the smallest for run NF1, which has no feedback.

Fig. 5 Kinetic pressure profile along the z-axis for five different models (see label) at four different timesteps. Run D presents the largest values while run NF1 presents values significantly smaller than the other runs. Moreover, kinetic pressure quickly drops at higher altitude for run NF1.

Fig. 6 Magnetic pressure profile along the z-axis for four different models (see label) at four different time steps. All runs show similar values in the midplane. For runs A, C1, and D, the profile tends to become shallower with time illustrating the transport of the magnetic flux towards higher altitude.

In all models, the thermal pressure ranges between a few 10^-13 and a few 10^-12 erg cm^-3, which corresponds to the typical pressure in the ISM (e.g., Ferrière 2001). Its variation with the altitude, z, closely follows the density variation with the notable exception of runs NF1 and A (i.e. no or randomly distributed supernovae). For the latter, the thermal pressure is roughly constant up to an altitude at which the density is about 0.1 cm^-3. This is clearly due to an efficient production of warm and hot gas by supernovae explosions. Note that run D presents more variability than run C1. This is likely a consequence of the supernovae being less spatially correlated to the sink particles. When a supernova explodes in a dense regions, it tends to mimic what happens for most supernovae in run C1. When it explodes in a diffuse regions, it tends to mimic run A, where most supernovae explode in the WNM, which has the largest volume filling factor.

The kinetic pressure is typically a few times larger than the thermal one and, depending on the model, reaches values of the order of ≃1 − 3 × 10^-12 erg cm^-3. These values are also very similar to what is reported in related studies (e.g., Joung et al. 2009; Kim et al. 2011). Interestingly, the scale height of P_kin is larger than the density scale height by a factor on the order of 1.5−2 for run C1 and up to 3−4 for run C2. In the case of scheme D, it even slightly increases with altitude. As expected, in the absence of supernovae feedback, P_kin drops to small values rapidly.

The magnetic pressure is comparable to the thermal pressure and reaches values of the order of a few 10^-12 erg cm^-3. It therefore contributes to support the galactic disc against gravity. Interestingly, while the magnetic intensity tends to decrease in the midplane as time goes on, it tends to increase with time at high altitude. This is a consequence of the generation of magnetic field through turbulence but also a consequence of the transport of the field lines by galactic outflows.

Fig. 7 Probability density function for MHD run C1 (upper panel) and hydrodynamical run C2 (lower panel). The small drop at n ≃ 2 − 3 cm-3 corresponds to the thermally unstable regime, which persists in spite of the strong turbulence. The two peaks correspond to the WNM and CNM phases.

Fig. 8 Mean temperature as a function of density for run C1 (upper panel) and C2 (lower panel). As expected the temperature lies mainly in three ranges of temperature, namely 106, 104 and 102 K, corresponding to the three phases of the ISM, the HIM, the WNM, and the CNM.

4. Multiphase ISM

In this section, we study the density, temperature and magnetic field distribution of the gas in the simulations. For conciseness, we restrict our attention to runs C1 and C2, and only stresses the most interesting differences with other runs. The comparison between runs C1 and C2 emphasizes the role that the magnetic field has to determine the characteristic of the ISM.

4.1. Density and temperature distributions

As is the case in other works (e.g., de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; Hill et al. 2012), the interstellar medium, which is initially uniform in density and temperature, quickly breaks up into a multiphase medium in which the density varies from less than 10^-3 cm^-3 to almost 10³ cm^-3, while the temperature can be as high as 10⁷ K and as low as a few tens of Kelvin.

Figure 7 shows the mass contribution of the various gas densities in the computational box. It is dominated by the dense gas (n> 10 cm^-3), but a non-negligible fraction lays at lower densities, which corresponds to the WNM regime and thermally unstable gas (0.1 <n< 10 cm^-3). Interestingly, there is a small dip in the thermally unstable regime (n ≃ 2 − 3 cm^-3). This indicates that in spite of the relatively high level of turbulence, typical of galactic discs (see below), the 2-phase structure (Wolfire et al. 2003) that would be obtained in a static medium, is not erased by the dynamical processes and able to persist. At early time, there is a little more dense gas in the hydrodynamical case. This is due to the magnetic support that reduces the amount of self-gravitating gas. For run D (not displayed here for conciseness), the transition between WNM and cold neutral medium (CNM) is less pronounced, which is a consequence of the stronger turbulence in this run. These behaviours are reminiscent of what has been found in colliding flow calculations, which attempted to model the ISM at scales on the order of 10−50 pc (e.g., Audit & Hennebelle 2005, 2010; Vázquez-Semadeni et al. 2006; Hennebelle et al. 2008; Heitsch et al. 2008; Banerjee et al. 2009; Inoue & Inutsuka 2012) and also in similar supernovae regulated galaxy simulations (e.g., Dib et al. 2006; Hill et al. 2012; Kim et al. 2013). In particular, in these simulations, it has been found that the ISM quickly breaks up into a multi-phase, clumpy medium where the two-phase behaviour (i.e. present an excess of gas in thermodynamical states close to the two stable branches of thermal equilibrium) is maintained, even though the medium is largely turbulent.

Figure 8 displays the temperature distribution as a function of density. It clearly shows the existence of three main domains corresponding to the hot ionised medium (T ≃ 10⁶ K), the warm neutral medium (T ≃ 10⁴ K), and the cold neutral medium (T< 100 K). This is in good agreement with the classical three phase model of the interstellar medium, as described earlier, for example, in McKee & Ostriker (1977). As noted by previous authors (e.g., Gazol et al. 2001), there is gas in the thermally unstable regions (T ≃ 10³K), whose existence is permitted by the turbulent motions. The two runs present similar distributions, although the transition between the warm and the cold phase (at density 1−10 cm^-3) is a little more shallow for the MHD run (C1) than for the hydrodynamical one (run C2). This is because the magnetic field contributes to the total pressure and can therefore stabilize the pieces of gas that are thermally unstable. Again in run D, the three regimes are less clearly separated though the global temperature range in the simulation is similar.

Fig. 9 Density (left column) and temperature (right column) fields in the equatorial plane for the MHD run C1 (upper pannel) and the hydrodynamical run C2 (lower panels). The figures illustrate the multi-phase nature of the interstellar medium. Most of the volume is filled by warm neutral gas with temperature of about 8000 K. A tiny fraction is occupied by the hot phase at temperature larger than 106 K.

The spatial distribution is illustrated by Fig. 9, which displays a cut through the equatorial plane of the density and temperature fields both for MHD run C1 and hydrodynamical run C2. As can be seen, most of the volume is found to be occupied by the warm neutral gas with temperatures of the order of 10⁴ K and densities of the order of 1 cm^-3. The hot gas, produced within supernovae explosions, occupies only a small fraction of the volume. Interestingly, the structures in the hydrodynamical and MHD runs have a slightly different shape. As already discussed in the previous section, the dense clouds in the hydrodynamical simulations tend to be more fragmented and on average slightly smaller (see also Fig. 1).

4.2. Velocity dispersion

In the supernovae regulated numerical simulations (e.g., 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), the kinetic energy, which decays through the turbulent cascade, is replenished by the supernovae explosions. The velocity dispersion in the computational box is the result of a balance between injection and dissipation, as emphasized, for example, in Mac Low & Klessen (2004), who present orders of magnitude suggesting that supernovae explosions can explain the velocity dispersion observed in the Milky Way to be of the order of 6 km s^-1.

Figure 10 displays the rms z-component of the velocity field weighted by density $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{z}}\mathrm{=}\sqrt{\frac{\sum {\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathit{\rho }\mathrm{d}\mathit{z}}{\sum \mathit{\rho }\mathrm{d}\mathit{z}}}\mathrm{·}& & \end{array}$(9)Close to the equatorial plane, σ_z is of the order of 4−5 km s^-1 for run C1 and 5−6 km ^-1 for run C2. At a higher altitude, the velocity dispersion increases to values of about 8−10 km s^-1 at 100 pc for run C1 and 8−15 km s^-1 for run C2. This is essentially due to the gas density getting lower at higher altitude. These values are again similar to what has been previously reported in similar studies (see e.g., Fig. 4 of Kim et al. 2013) and can be understood by relatively simple considerations. Following Mac Low & Klessen (2004), we can simply estimate the amount of mechanical energy, which dissipates in the turbulent cascade as Ė_diss = Mσ²/τ, where M is the total mass of the system, τ is the crossing time, and σ the total velocity dispersion. Assuming that τ = h/σ, where h is the disc scale height, we get Ė_diss = Mσ³/h, where Ė_diss is the energy dissipated per units of time. The amount of energy, which is injected by the supernovae into the system, is simply Ė_inj = ϵṄ_sn × 10⁵¹ erg, where Ṅ_sn is the density of supernovae per units of time and ϵ is the efficiency at which turbulence is triggered. Equating these two rates, we get $\begin{array}{ccc}\mathit{\sigma }\mathrm{=}{\left(\frac{\mathit{ϵ}\mathit{Ṅ}\mathrm{sn}\mathit{h}\mathrm{\times }{\mathrm{10}}^{\mathrm{51}}\hspace{0.17em}\mathrm{erg}}{\mathit{M}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}& & \end{array}$(10)To estimate the value of σ, we take values typical for the Milky way. These values are also representative of our simulation parameters. We take a mass of 10¹⁰ M_⊙, a frequency of supernovae Ṅ_sn = 1/50 yr^-1, a height h = 100 pc, and an efficiency ϵ = 0.1, we get $\begin{array}{ccc}\mathit{\sigma }\mathrm{\simeq }\mathrm{8}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}{\left(\frac{\mathit{ϵ}}{\mathrm{0.1}}\frac{\mathit{Ṅ}\mathrm{sn}}{\mathrm{1}\mathit{/}\mathrm{50}{\mathrm{yr}}^{-1}}\frac{\mathit{h}}{\mathrm{100}\mathrm{pc}}\frac{{\mathrm{10}}^{\mathrm{10}}{\mathit{M}}_{\mathrm{\odot }}}{\mathit{M}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{.}& & \end{array}$(11)This value is thus in good agreement with the velocity dispersion inferred from the simulations and from the observations. Due to the weak dependence in all parameters (to the power 1/3), it is relatively unsurprising to find that the velocity dispersion generally does not undergo large variations. We note that the fluctuations in the hydrodynamical run C2 (lower panel) appear to be quite large with respect to the MHD run C1 (upper panel). This is likely a consequence of the higher SFR (see below) found in the hydronamical case. This results in a stronger feedback.

Fig. 10 Root mean square value of vz as a function of z (see text) for MHD run C1 and hydrodynamical run C2. Typical values are about 4−5 km s-1 at the midplane.

Fig. 11 Root mean square velocity as a function of density for MHD run C1 and hydrodynamical run C2.

Fig. 12 Mean magnetic intensity as a function of z and as a function of density. Note in particular that typical values of about 5 μG are being obtained in the midplane. Between 1 and 103 cm-3, the magnetic intensity weakly varies with the density.

Finally, since stars mainly form in the dense gas, it is important to understand the star formation process to know more accurately how velocity dispersion depends on the gas density. For that purpose Fig. 11 displays the rms velocity field (taking the three components into account) weighted by density as a function of density in the whole computational box. As expected, the velocity dispersion is weaker in the dense gas than in the diffuse one by a factor of about 2. Typical velocity dispersion in the dense gas is on the order of 4−5 km s^-1.

4.3. Distribution of magnetic intensity

Figure 12 shows the magnetic intensity as a function of altitude (upper panel) and density (lower panel). At the mean density of the galactic disc, i.e. n ≃ 2 − 3 cm^-3, the magnetic intensity is about 3−5 μG, which is also the mean value up to an altitude of about 100 pc. This is coherent with the values of about 5 μG reported by Heiles & Troland (2005). At higher densities, the magnetic intensity increases and reaches values of about 10 μG at densities of about 10² cm^-3. Note that this corresponds to a rather shallow variation of the magnetic intensity with density as observed in the diffuse gas (e.g., Troland & Heiles 1986). On the one hand, the exact reason of this weak correlation is most likely due to the Lorentz force, that resists contraction perpendicular to the field lines (Hennebelle & Pérault 2000; Passot & Vázquez-Semadeni 2003). On the other hand, it is also partly due to the turbulent diffusivity, which has also been observed to play an important role in numerical simulations (e.g., Lazarian & Vishniac 1999; Santos-Lima et al. 2010; Hennebelle et al. 2011). For densities below n ≃ 10^-2 cm^-3, a steep drop is observed with density. This is due to the fast expansions produced by supernovae explosions, which tend to dilute the magnetic intensity very significantly.

5. Star formation rate, sink mass function, and outflows

We now investigate the characteristics of the star formation in the simulations. This is achieved through the sink particles described in Sect. 2.3. We first quantify the total mass of the sink particles, which represents the SFR in the simulations. We then study the mass distribution, which is the sink mass function of some of our models. Finally, we study the outflows, which are launched at high altitude and eventually escape the computational box.

5.1. Star formation rates

Fig. 13 Total mass of sink particles as a function of time for the various models. Upper panel shows the influence of the supernova feedback scheme, while middle and lower panels show the influence of the magnetic field. The lower panel shows the SFR (i.e. time derivative of mass) corresponding to middle panel. In the upper panel, run B corresponds to the solid curve, which starts at 120 Myr. The runs with feedback present SFR that are typically 10 to 30 times smaller. Magnetic field reduces the SFR by a factor on the order of 2.

Figure 13 shows the total mass of the sink particles as a function of time in the simulations. Upper panel shows the influence of the sink particle prescription, while lower panel shows the influence of the magnetisation. Before commenting on the difference between the various models, we first discuss the main trends and numbers. For all models (except run B), accretion onto sink particles starts between 20 and 30 Myr. In about 10−20 Myr, the total accreted mass varies between a few 10⁵ to ≃10⁷ M_⊙. As discussed below, these differences are due to the various feedback prescriptions and magnetisations. At later times, all models tend to reach a phase of stationary accretion at a rate, which ranges from about 10^-2 to ≃10^-1 M_⊙ yr^-1. It is worth comparing these values with the typical 3 M_⊙ yr^-1 at which Milky way is forming stars. To do so, one must first correct for the volume of our computational box. Since most stars form in the Milky way inside the solar circle whose radius is about 8 kpc and since the size of the computational domain is equal to 1 kpc, a geometrical factor of π × 8² ≃ 200 should be taken into account. However, it should also be accounted for that the efficiency of the mass eventually accreted into the stars is only a fraction of the mass that is accreted onto the cores. This value is not known with great accuracy but has been estimated to be of the order of 1/3 (e.g., Alves et al. 2007). We note that the sink particles used in this study are at this stage much larger than dense molecular cores. Therefore, it could be that the efficiency should be even lower than this value. Combining these two numbers, we find that for a galaxy like the Milky way, our models would predict a SFR ranging from about 1 to 20 M_⊙ yr^-1. Given the large uncertainties, the first value appears in reasonable agreement with the galactic one.

5.1.1. Influence of feedback prescriptions

All models displayed in the top panel of Fig. 13 have initial conditions identical to run C1, which have an initial magnetic field whose initial intensity in the midplane is about 2.5 μG.

First of all, the large difference between the solid line (run NF1) and the dashed line (run C1) confirms the drastic influence of the feedback on the SFR, which is reduced by a factor of 20-30. This constitutes a strong hint that feedback can be largely responsible to solve the long standing issue of the so-called Zuckermann-Evans catastrophe (Zuckerman & Evans 1974). If all the molecular gas of the Milky way was collapsing in a free-fall time, about 300 M_⊙ yr^-1 of stars would form in the Galaxy.

Second of all, when the supernovae are not correlated to density (run A), the feedback is not only unable to reduce the SFR but also this latter is even slightly higher. Since most of the volume is occupied by warm gas, most of the supernovae therefore explodes in low density regions. Their net effect is to further compress the dense gas. When the supernovae are correlated with the density peak (run B), it takes a long time before sinks can form because the dense gas is efficiently dispersed. However, once sinks start forming, they are unable to reduce the accretion rate since the supernovae are not correlated locally in space and in time with accretion but simply with the densest cell in the simulation. Therefore SFRs comparable to the run without feedback are obtained.

Third, runs C4 and D show that SFR larger by a factor of 2−3 are obtained when the feedback is either purely thermal or less tightly correlated to the sink particles. Given that these two aspects are largely uncertain, this illustrates the limit of this modeling and suggests that the typical accuracy of these models is at best on the order of a factor 2−3. Note that another severe source of uncertainties comes from the time at which supernovae are introduced. In particular, if a delay of tens of Myr is introduced, SFR comparable to the ones of run NF1 are obtained. This suggests that to get more accurate models, it is necessary to have a better description of the small scales and in particular of the formation and evolution of massive stars up to the point where they explode. Ideally, this would require running a set of specific small scale simulations to quantify the impact of the feedback more accurately.

5.1.2. Influence of magnetisation

All models displayed in middle and bottom panel of Fig. 13 are performed with either no feedback (runs NF1 and NF2) or with the same feedback scheme (scheme C). Different levels of magnetisation are compared.

In the hydrodynamical run C2, stars start forming a few Myr before run C1. The SFR is initially significantly reduced compared to run C2. At later time, they become, however, comparable. These effects are a consequence of the magnetic support, which contribute to resist the gravitational contraction but also to the probability density function (PDF) that is narrower in the presence of a magnetic field (Molina et al. 2012). A similar effect is obtained for the two runs without feedback (runs NF1 and NF2), for which it is seen that the SFR is a little higher in the hydrodynamical case than in the MHD one.

Interestingly, even when the magnetic field is rather weak (0.5 μG initially), it still has a visible impact and reduces the SFR by a factor of about 50% during the first 20 Myr after star formation has started. This is because the magnetic field is quickly amplified to larger values.

This effect is quantitatively comparable to what has been inferred at smaller scales by various teams, who investigate star formation in substancially magnetized, though supercritical clouds. For example, Price & Bate (2008) simulate the collapse of a self-gravitating clump while Dib et al. (2010) and citetpadoan+2011 perform self-gravitating, MHD calculations within periodic boxes. They all infer that magnetic field reduces the SFR by half. The exact reason of this lower value has not been analysed in great detail so far, but it is likely a consequence of the magnetic support, which tends to resist gravity and the somehow narrower density PDF, which tends to reduce the SFR (Hennebelle & Chabrier 2013).

Fig. 14 Sink mass spectra for four models at four time steps. In the case with no feedback, there are typically few massive sink particles, while a broader distribution develops when feedback is included.

Fig. 15 Mean flux of mass along the z-axis for five different models (see label) at four different time steps. The largest values are obtained for scheme D.

5.2. Mass function of sink particles

Figure 14 shows the sink mass function for various models and at four time steps from which one can verify that the trends discussed below are not due to a time selection. For runs C1 and C2, a large number of sinks form (about 400 and 700 respectively for run C1 and run C2). Their masses span about three orders of magnitude. Given the limited numerical resolution of the present study, many features of the distribution must be taken with great care. In particular, the peak at about 10³ M_⊙ would certainly shift to smaller values in more resolved runs (e.g., Hennebelle & Audit 2007). There is a possible trend for a power-law developing at large masses (in the range ≃3 × 10³ − 10⁴ M_⊙) with an exponent compatible with ≃−1. However, the limited resolution precludes a firm conclusion. We note that the number of fragments has also been found to be reduced by a factor of about two in massive collapsing magnetised clumps. This is a clear consequence of the cold gas being more coherent and less fragmented. As noted previously, the reason is that the magnetic field makes the flow more coherent, since it tends to connect fluid particles linked by the magnetic field lines (e.g., Hennebelle 2013).

The sink mass function obtained when no feedback (run NF1) is included is quite different. There are much less sink particles (about 70), and most of them have a mass larger than 10⁴ M_⊙. Indeed, the most massive sink particle has in this case a mass equal to a few 10⁶ M_⊙. This behaviour is again a consequence of the absence of feedback. The gas tends to concentrate in a few locations under the influence of gravity. The sink mass function obtained for run D is inbetween the one obtained for run C1 and run NF1. This illustrates again that the stellar feedback is less efficient in supporting the gas against gravitational collapse in run D.

5.3. Galactic outflows

The existence of galactic outflows is now well established (see Veilleux et al. 2005, for a recent review) in many galaxies. The typical scale height at which these outflows are observed is of the order of several tens of kpc, which is much larger than the scale of the present simulation. The box length is also equal to only 1 kpc, and we do not have a proper halo structure, which influences the flow launching (e.g., Dubois & Teyssier 2008). It is nevertheless worth quantifying them, because supernovae are believed to be largely responsible of their launching.

Figure 15 shows the mass flux as a function of the altitude z for four time steps and four models. As can be seen, it varies significantly with time and altitude from typically a few 10^-2 to a few 10^-1 M_⊙ yr^-1. Taking into account that the surface of the box is 1 kpc², this would lead to a flux of about π × 8² ≃ 200 times larger for a galaxy similar to ours, which is a few solar mass to a few tens or solar mass per year. These values are typical of what is measured for galactic outflows (Veilleux et al. 2005).

Another interesting trend is that the mass flux broadly correlates with the SFR (see Fig. 13). For run C1 (top panel), the peak value of the mass flux is about three times smaller than what is obtained for run D (third panel). Comparing the mean value at the computational box edges (z = 500 pc), the ratio between the fluxes of run C1 and D leads to somewhat larger values of about 5−10. This is likely a consequence of the dual role of supernovae explosions, which are responsible for the regulation of star formation through energy and momentum injection in the dense gas but also for the launching of the galactic outflows through injection onto the diffuse gas. Since the two processes are linked, it is expected that larger SFR lead to stronger outflows, as they imply more feedback. Indeed the SFR ratio for run D and C1 is about 3 (from upper panel of Fig. 13), which is comparable with the value of 3 quoted above, but a little too small to explain the second value obtained at the box edges. This may indicate that another effect must be considered. We believe that more energy and momentum tend to be injected in the diffuse gas in run D than in run C1 since the supernovae explode further from the sink particles than in run C1. Since the outflows are primarily made by diffuse quickly expanding material, it seems reasonable that the efficiency with which they are produced is higher in run C1 than in run D.

6. Conclusion

We have performed a series of numerical simulations describing a galactic disc regulated by supernovae feedback at kpc scale. Our simulations include both magnetic field and self-gravity. In particular, we have explored the influence of various schemes to prescribe the supernovae feedback. Our simulations reproduce many features already found by other authors, such as multi-phase density, temperature distributions or velocity dispersion, typically of the order of 5 km s^-1 in the galactic plane. Our results are as follows. When the supernovae are randomly distributed, they drive the interstellar turbulence but are unable to resist self-gravity efficiently, and the SFR is as high (even slightly higher) as when no feedback is included. When supernovae are correlated to the density peaks, they efficiently limit star formation by preventing the gas to become too dense. However, as time goes on, dense gas eventually develops. When sink particles are being introduced, then the SFR is as high as its value without feedback.

When supernovae are spatially and temporally correlated to star formation events, the SFR is significantly reduced by a factor of the order of ten or more. However, we find that the exact implementation of the supernovae does influence the galactic disc structure and the SFR significantly. In particular, if the supernoave are distributed in a shell of about 16 pc around the sink particles, the accretion rate is higher by a factor of about three than if they are randomly placed within a sphere of radius equal to 16 pc. In a similar way, if the feedback is purely thermal, the SFR is about twice than if it had 5% of kinetic feedback. This implies that detailed knowledge of how the feedback operates on small scales is mandatory to understand its impact with sufficient precision. In particular, the correlation between the massive stars and the dense star forming gas should be determined using small scale simuations.

The magnetic field has a significant impact. It delays and reduces star formation by a half. It also tends to reduce the number of star formation regions (e.g., sink particles) by half, therefore producing slightly bigger star forming regions. Finally, it should be kept in mind that the magnetic field has an important impact on the fragmentation of massive cores that it tends to reduce significantly (Commerçon et al. 2011; Myers et al. 2013). This implies that more massive stars form when a magnetic field is strong. Since feedback is a non-linear function of the stellar masses and since feedback drastically influences the galactic structure and evolution, it is likely that the impact of the magnetic field on galaxy evolution is probably even larger than what is estimated here.

Acknowledgments

We thank the anonymous referee for a careful reading of the manuscript which has significantly improved the paper. This work was granted access to HPC resources of CINES under the allocation x2014047023 made by GENCI (Grand Equipement National de Calcul Intensif). P.H. acknowledge the finantial support of the Agence National 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: The issue of numerical convergence

To investigate the issue of numerical convergence, we present the various profiles for run C1b, which is identical to run C1 but has an effective resolution that this twice as large. As can be seen from Fig. A.1 and Figs. 3−10, which displayed the results for run C1, the profiles present some moderate differences for the two cases implying that numerical convergence is not fully reached. In particular, the density profile is slightly less peaked for run C1b than for run C1. Similarly, the rms velocity is about 5 km s^-1 for run C1b, while it is equal to about 4−4.5 km s^-1 for run C1. The pressures, however, present very comparable values and profiles between the two runs. Altogether, this shows that the quantities, which characterize the disc structure, are reasonably described at the resolution used in the paper.

Fig. A.1 Mean density, z-velocity dispersion, kinetic, thermal, and magnetic profile along the z-axis for run C1b at four different time steps. These results should be compared with the corresponding quantities for run C1 displayed in Figs. 3−10.

All Tables

All Figures

Fig. 1 Column density along the z-axis (left panels) and along the y-axis (right panels) for MHD run C1 (upper panels) and hyrodynamical run C2 (lower panels). The left panels illustrate the complex multi-phase structure of the galactic plane, while the right panels show the stratification induced by the gravitational field of the galaxy. In the hydrodynamical run, the interstellar medium is more fragmented than in the MHD one.

In the text

Fig. 2 Column density along the z-axis (left panels) and along the y-axis (right panels) for the run with no feedback, NF1, (upper panels) and the run with scheme D (lower panels). Comparison with Fig. 1 clearly shows the drastic impact of the feedback. In particular, the run NF1 shows thin and long filaments in which self-gravitating fragments develop, indicating that the gas is primarily organised by self-gravity and undergoes run away collapse. The galactic disc is considerably thinner than when supernovae are included. On the other hand, run D shows a broader galactic disc than in runs C1 and C2, illustrating the importance of the spatial correlation between supernovae explosions and dense gas.

In the text

	Fig. 3 Mean density profile along the z-axis for five different models (see label) at four different time steps. The disc profile is much thinner for runs NF1 and A than for runs C1, C2, and D.
In the text

	Fig. 4 Thermal pressure profile along the z-axis for five different models (see label) at four different timesteps. The largest values are obtained for run D and the smallest for run NF1, which has no feedback.
In the text

	Fig. 5 Kinetic pressure profile along the z-axis for five different models (see label) at four different timesteps. Run D presents the largest values while run NF1 presents values significantly smaller than the other runs. Moreover, kinetic pressure quickly drops at higher altitude for run NF1.
In the text

	Fig. 6 Magnetic pressure profile along the z-axis for four different models (see label) at four different time steps. All runs show similar values in the midplane. For runs A, C1, and D, the profile tends to become shallower with time illustrating the transport of the magnetic flux towards higher altitude.
In the text

	Fig. 7 Probability density function for MHD run C1 (upper panel) and hydrodynamical run C2 (lower panel). The small drop at n ≃ 2 − 3 cm-3 corresponds to the thermally unstable regime, which persists in spite of the strong turbulence. The two peaks correspond to the WNM and CNM phases.
In the text

	Fig. 8 Mean temperature as a function of density for run C1 (upper panel) and C2 (lower panel). As expected the temperature lies mainly in three ranges of temperature, namely 106, 104 and 102 K, corresponding to the three phases of the ISM, the HIM, the WNM, and the CNM.
In the text

Fig. 9 Density (left column) and temperature (right column) fields in the equatorial plane for the MHD run C1 (upper pannel) and the hydrodynamical run C2 (lower panels). The figures illustrate the multi-phase nature of the interstellar medium. Most of the volume is filled by warm neutral gas with temperature of about 8000 K. A tiny fraction is occupied by the hot phase at temperature larger than 106 K.

In the text

	Fig. 10 Root mean square value of vz as a function of z (see text) for MHD run C1 and hydrodynamical run C2. Typical values are about 4−5 km s-1 at the midplane.
In the text

	Fig. 11 Root mean square velocity as a function of density for MHD run C1 and hydrodynamical run C2.
In the text

	Fig. 12 Mean magnetic intensity as a function of z and as a function of density. Note in particular that typical values of about 5 μG are being obtained in the midplane. Between 1 and 103 cm-3, the magnetic intensity weakly varies with the density.
In the text

Fig. 13 Total mass of sink particles as a function of time for the various models. Upper panel shows the influence of the supernova feedback scheme, while middle and lower panels show the influence of the magnetic field. The lower panel shows the SFR (i.e. time derivative of mass) corresponding to middle panel. In the upper panel, run B corresponds to the solid curve, which starts at 120 Myr. The runs with feedback present SFR that are typically 10 to 30 times smaller. Magnetic field reduces the SFR by a factor on the order of 2.

In the text

	Fig. 14 Sink mass spectra for four models at four time steps. In the case with no feedback, there are typically few massive sink particles, while a broader distribution develops when feedback is included.
In the text

	Fig. 15 Mean flux of mass along the z-axis for five different models (see label) at four different time steps. The largest values are obtained for scheme D.
In the text

	Fig. A.1 Mean density, z-velocity dispersion, kinetic, thermal, and magnetic profile along the z-axis for run C1b at four different time steps. These results should be compared with the corresponding quantities for run C1 displayed in Figs. 3−10.
In the text