Large-scale turbulent driving regulates star formation in high-redshift gas-rich galaxies II: Influence of the magnetic field and the turbulent compressive fraction

The observed star formation rate (SFR) in galaxies is well below what it should be if gravitational collapse alone were at play. It has recently been shown that one candidate that might regulate star formation, the feedback from massive stars, is suitable only if the mean column density at the kiloparsec scale is lower than $\approx 20 M_\odot\cdot\mathrm{pc}^{-2}$. On the other hand, intense large-scale turbulent driving might slow down star formation in high-density environments to values that are compatible with observations. In this work, we explore the effect of the nature and strength of the turbulent driving, as well as the effect of the magnetic field. We performed a large series of feedback-regulated numerical simulations of the interstellar medium in which bidimensional large-scale turbulent driving was also applied. We determined the driving intensity needed to reproduce the Schmidt-Kennicutt relation for several gas column densities, magnetization, and driving compressibility. We confirm that in the absence of turbulent forcing and even with a substantial magnetic field, the SFR is too high, particularly at a high column density, compared to the Schmidt-Kennicutt relation. We find that the SFR outcome strongly depends on the initial magnetic field and on the compressibility of the turbulent driving. As a consequence, a higher magnetic field in high column density environment may lower the energy necessary to sustain a turbulence that is sufficiently intense to regulate star formation. Stellar feedback does not seem to be sufficient to regulate star formation in gas-rich galaxies where large-scale turbulent driving may be needed. The sources of this large-scale turbulence as well as its characteristics, such as its intensity, compressibility, and anisotropy, need to be understood and quantified.


Introduction
Stars are complex objects that connect small-and large-scale physics. Key characteristics of galaxies, such as their structure and shape, the amount of gas, their chemistry, and the velocity field, depend on the efficiency and on the rate of star formation. However, understanding and predicting the star formation rate (SFR) remains both an observational and theoretical challenge. On the theoretical side, the difficulty is the very broad span of scales involved that cannot be captured in just one simulation (Dubois & Teyssier 2008;Hopkins et al. 2011), and the high degeneracy between the multiple processes at play.
A considerable amount of work has been done to understand the interstellar medium (ISM) at the kiloparsec scale, which is an obligatory step on the path towards determining how the various physical processes are combined and influence the structure of the ISM. Since modelling a full galaxy with parsec or subparsec resolution is very challenging, many studies have considered computational boxes with a size of about 1 kpc. The advantage is that this allows describing the disk stratification and the large molecular cloud complexes while still resolving the clouds sufficiently well.
The complexity has been progressively increased in these simulations. The first series of calculations (e.g. Rosen & Bregman 1995;de Avillez & Breitschwerdt 2005;Joung & Mac Low 2006;Hill et al. 2012) solved the equations of Magnetohydrodynamics with an interstellar cooling curve in the presence of an external gravitational potential representing the influence of stars and dark matter. In these calculations, supernova (SN) explosions were placed randomly at a rate aiming to reproduce the expected galactic SN rate. These models produced appealing turbulent multi-phase ISM structures with a velocity dispersion close to the observations. In a second series of calculations (e.g. Kim et al. 2011Kim et al. , 2013Hennebelle & Iffrig 2014;Walch et al. 2015;Iffrig & Hennebelle 2017;Gatto et al. 2017), self-gravity was introduced, together with sink particles that were to represent stellar clusters. This important step allowed addressing the question of the SFR and also the issue of the galactic thickness, which requires the correct balance between gravity and turbulence. In these calculations, the stellar feedback essentially due to SNe was spatially correlated to the sink particles, but no time delay of a few to some tens of million years was introduced. That is to say, the feedback was immediately delivered when enough mass was accreted to form a massive star. Whereas these mod-Article number, page 1 of 21 arXiv:2305.18012v2 [astro-ph.GA] 6 Jun 2023 A&A proofs: manuscript no. brucy2023 els could successively reproduce the star formation rate typical of Milky Way-type galaxies, it was noted that the scheme used for the SN explosions had a strong influence on the results. In particular, if the SN exploded far away from its parent sink particle or simply was forced to explode in the diffuse gas, it would become very inefficient.
In the last generation of models (e.g. Padoan et al. 2016;Kim & Ostriker 2017;Colling et al. 2018;Kannan et al. 2020;Ostriker & Kim 2022;Rathjen et al. 2023), the massive stars were tracked more self-consistently and the SNe explode after a few millions to tens of million years, depending on the mass of their progenitor. In combination with other processes such as stellar wind, HII regions, or galactic shear, it has been concluded by various teams that the SFR was compatible with Milky Way values.
A fundamental assumption made in these feedback-regulated models is that the ISM regulation is relatively local. More precisely, the energy injected in the medium at a given point is entirely due to stellar feedback, and it typically comes from less than 1 kpc. However, several recent studies have shown that injection of turbulence from large galactic scales has to be taken into account in order to explain the observed velocity dispersion and SFR, particularly in gas-rich galaxies Renaud et al. 2012;Goldbaum et al. 2015Goldbaum et al. , 2016Krumholz et al. 2018;Meidt et al. 2020;Nusser & Silk 2022). There is also evidence that large-scale driving is needed to reproduce the statistics of the dense gas, in particular, the large-scale structures in the Large Magellanic Cloud (e.g. Colman et al. 2022). Several observational studies have indeed derived turbulence-injection scales of more than 1 kpc (Dib et al. 2021;Chepurnov et al. 2015) and indicated large-scale driving (Szotkowski et al. 2019;Besserglik & Goldman 2021). Possible sources of turbulence include the orbital energy (Wada et al. 2002;Agertz et al. 2009) and the mass accretion onto the galaxies (Klessen & Hennebelle 2010;Forbes et al. 2023). The former in particular requires a mechanism such as an instability to convert this source of free energy into turbulent energy, which might be the gravitational instability, for instance (Wada et al. 2002;Agertz et al. 2009;Goldbaum et al. 2015;Fensch et al. 2023), or the magneto-rotational instability (Piontek & Ostriker 2007).
An important challenge of numerical simulations for the SFR in particular is reproducing the Schmidt-Kennicutt relation (hereafter SK relation;Kennicutt 1998;Kennicutt & Evans 2012;Kennicutt & De Los Reyes 2021) that links the SFR to the column density of gas. In a recent paper, Brucy et al. (2020) (hereafter Paper I) tested the effect of the injection of turbulence by adding a large-scale turbulent driving similar to the one used by Schmidt et al. (2009). Importantly, Brucy et al. (2020) have shown that whereas the SK relation can be reproduced if the added turbulence is strong enough, the SFR appears to be too high, particularly for gas-rich galaxies, when stellar feedback alone is accounted for.
The goal of the present paper is to complement the study presented in Paper I by studying the effects of the magnetic field strength and the compressive fraction of the turbulent driving on the SFR. Both aspects are known to significantly influence the gas distribution (e.g. Molina et al. 2012) and likely influence the SFR values. We ran simulations of a local region of a galactic disk within a cubic box of 1 kpc. We used a numerical setup very similar to the one used by Colling et al. (2018) and Paper I. We simulated regions of galaxies with a wide range of gas column densities that are representative for Milky Way-like galaxies up to gas-rich galaxies at redshift z = 1-3 (Genzel et al. 2008(Genzel et al. , 2010Daddi et al. 2010).
The paper is organised as follows: We start by discussing the orders of magnitude of the turbulent energy that can be dissipated by various processes in the ISM in section 2. Next, we present our numerical setup in section 3, and we explore the effects of the turbulence and the magnetic field in sections 4 and 5, respectively. We then discuss the main goal of this study, that is, the reproduction of the SK relation in section 6. The caveats of this study are discussed in section 7. Finally, we discuss the results and conclude in section 8.

Energy dissipation and turbulent driving in the ISM, orders of magnitude
As discussed above, we further investigate whether large-scale turbulent driving can help explain the SK relation. It is thus of primary importance to discuss the possible sources of turbulent driving in galaxies. Below, we discuss three mechanisms, namely gravitational instability (Wada et al. 2002;Agertz et al. 2009), accretion-driven turbulence (Klessen & Hennebelle 2010;Forbes et al. 2023), and SN driving (e.g. Kim & Ostriker 2017;Gatto et al. 2017). We estimate the amount of energy P that can be extracted per unit of time through each of the mechanisms, as well as their efficiency ϵ, which is the fraction of that energy that is eventually converted into turbulent kinetic energy. We start by estimating the amount of turbulent energy that is dissipated per unit of time in a region of size L ≃ 1 kpc within a galactic disk, with a column density Σ = 10 − 100 M ⊙ pc −2 and a velocity dispersion σ ≃ 10 − 100 km s −1 . Since the mass contained in a surface L 2 is ΣL 2 , whereas the specific kinetic energy is σ 2 /2 and the crossing time is L/σ, we obtain The largest amount of energy that can be extracted per unit of time from galactic rotation is typically equal to the rotation energy V 2 rot /2, divided by a rotation time that is simply given 2πR gal /V rot , where R gal is the galactic radius. This leads to ≃ 2 × 10 40 erg s −1 where Σ ⋆ is the stellar column density. The efficiency of turbulent driving ϵ rot through gravitational instability clearly depends on the disk stability, which is quantified by the Toomre parameter (Toomre 1964), where κ is the epicyclic frequency. Since a galactic disk is composed of gas and stars, it is likely that the efficiency ϵ rot depends on the gas column density Σ and the gas fraction Σ/(Σ + Σ ⋆ ) (Fensch & Bournaud 2021). Galaxies with a high gas fraction are more prone to gravitational instabilities, and we can therefore expect the efficiency ϵ rot in these galaxies to be high, probably close to 1. However, Milky Way-type galaxies can make several rotations without shrinking significantly due to angular momentum loss. For these galaxies, we expect the efficiency ϵ rot to be well below 1. Accretion at a given mass rateṀ accr and a infall velocity v inf not only brings mass into the galaxy, but also delivers energy at a rateṀ accr v 2 inf /2. To estimate the mass accretion rateṀ accr , we assume that the total amount of gas is roughly constant, that is, the accretion of gas balances the SFR, We estimate the surface density of the SFR by assuming that the SK relation (Kennicutt & Evans 2012) is verified, On the other-hand, we compute the infall velocity from the gravitational potential, where Σ dm is the column density of the dark matter. We thus find that It is complicated to determine the efficiency of turbulent driving by accretion, ϵ acc . Based on colliding-flow simulations, Klessen & Hennebelle (2010) estimated it to be proportional to the density contrast between the accreted material and the medium in which the turbulence is driven. Forbes et al. (2023) have estimated it in the context of a galactic disk and found that it varies with galactic radii, with a typical efficiency of about ϵ acc ≃15-20%. The stationary hypothesis is rather conservative, and it is possible that the accretion rate is much higher, in which case the resulting turbulent driving will be stronger as well.
The amount of energy injected per unit of time by SN remnants is given by the number of SNe produced per unit of time multiplied by the amount of energy a single SN delivers to the ISM. This is roughly equal to (Σ SFR L 2 )/120 M ⊙ , as given by the Salpeter initial mass function, multiplied by 10 51 erg.
By again assuming that the SFR is in line with the SK relation (Eq. (5)), we can estimate The efficiency of turbulent driving by SN remnants, ϵ sn , has been estimated from kiloparsec turbulent box calculations (e.g. Iffrig & Hennebelle 2017) to be about a few percent.
To summarise, the three envisaged mechanisms could all provide a substantial amount of energy to drive turbulence and compensate for turbulent dissipation. They present different scaling with Σ, however, more precisely, P rot ∝ Σ, P acc ∝ Σ 2.4 and P sn ∝ Σ 1.4 . We stress that given the orders of magnitude and the dependence on Σ, it is likely the case that for high values of Σ, energy rotation of both gas and stars is likely the dominant mechanism that triggers turbulence. This is also the source of energy envisioned by Krumholz et al. (2018) and Nusser & Silk (2022).

Numerical setup
We used the RAMSES code (Teyssier 2002) to solve the equations of Magnetohydrodynamics with a Godunov solver (Fromang et al. 2006) on a uniform cubic grid of 256 3 cells with periodic boundaries in the x and y directions parallel to the disk and open boundaries in the vertical direction. The box represents a cubic region of the galactic disk of size L = 1 kpc, so that the resolution is about 4 pc. The impact of the resolution on the results is discussed in Appendix C.

Initial conditions
We used initial conditions similar to those described in previous works (Colling et al. 2018;Brucy et al. 2020). The interstellar gas was initially distributed as a Gaussian along the z-axis, with n 0 a free density parameter and z 0 = 150 pc the typical scale height. The column density of gas (hydrogen and helium), integrated along the z-axis (perpendicular to the disk) is then where m p = 1.4×1.66·10 −24 g is the mean mass of a gas particle. We set the initial temperature at 8000 K to match the typical value of the temperature of the warm neutral medium phase of the ISM. An initial synthetic turbulent velocity field with a 3D dispersion of 5 √ n 0 /1.5 km · s −1 was also added. It presents a Kolmogorov power spectrum and was generated using a random phase (Kolmogorov 1941). The dependence on n 0 ensures that the initial ratio of the kinetic and gravitational energy due to the gas stays the same. Finally, we added a Gaussian magnetic field, oriented along the x-axis, where B 0 is a parameter of our simulation with a value of a few microgauss (see Table 1).

Numerical models
The gas is subject to an external gravitational potential, corresponding to old stars and dark matter, of the form with a 1 = 1.42 × 10 −3 kpc/Myr 2 , a 2 = 5.49 × 10 −4 Myr −2 , and z 0 = 0.18 kpc (Kuijken & Gilmore 1989;Joung & Mac Low 2006). This gravitational force adds up to the self-gravity of the gas. Sink particles (Bleuler & Teyssier 2014) were used to follow the dense gas and model star formation. Sink creation was triggered when the gas number density passed a threshold of 10 3 cm −3 . The simulation with an initial column density higher than 50 M ⊙ · pc −2 used a higher threshold of 5 × 10 3 cm −3 . Since the SFR in these simulations is high, this allowed us to reduce the number of sinks without changing the SFR much (Colling et al. 2018). Gas with a density higher than the threshold within a radius of four cells (16 pc) was then progressively accreted into the sink, with no more than 10 % of the gas being accreted in one time step.
The gas cooling was as described in Audit & Hennebelle (2005), which takes into account all the standard processes taking place in the atomic gas. It is essentially identical to the cooling described in Koyama & Inutsuka (2000), for instance.

Stellar feedback
The simulations included models for the formation and expansion of HII regions, SN explosions, and the far-ultraviolet (FUV) feedback. They are the same as in Paper I. The HII and SN feedbacks were attached to the sinks and were described in length in Colling et al. (2018) and references therein. Each time a sink had accreted a mass of 120 M ⊙ , a massive star particle with a mass randomly determined from the Salpeter initial mass function (Salpeter 1955) between 8 and 120 M ⊙ was created and attached to the sink. The lifetime τ * of this star was computed using the model with τ 0 = 3.265 Myr, M 0 = 148.16 M ⊙ , a = 0.238, and b = 2.205 (Woosley et al. 2002). When this massive star reached the end of its lifetime τ * , it exploded in a random location within a sphere with a radius τ * × 1 km · s −1 . The gas that was located inside a sphere of 12 pc radius around the location of the SN was heated up to inject a thermal energy of 10 51 erg. Since the cooling radius was almost never resolved, the gas would cool down immediately and the SN would have no effect. To avoid this, the explosion also injected 4 × 10 43 g · cm · s −1 momentum into the same region. When SNe explode in a low-density environment, a very high temperature and very high velocities can be generated, which translate into very short time steps due to the Courant condition. A limitation of the temperature to 10 6 K and of the velocity to 300 km · s −1 generated by a SN explosion was therefore implemented. We discuss the impact of this limitation on our results in Appendix B.
We also included self-consistent feedback from HII regions. Energy and momentum were injected according to the flux of ionising photons emitted by the massive star (Vacca et al. 1996). The evolution of the HII regions themselves was computed via the so-called M1 radiative transfer method (Rosdahl et al. 2013).
The FUV heating was uniform. However, it was not kept constant at the solar neighbourhood value because young O-B stars contribute significantly to the FUV emission. As a first approximation, we considered the UV heating to be proportional to the SFR (Ostriker et al. 2010). The mean FUV density relative to the solar neighbourhood value G ′ 0 can then be written as In our model, G ′ 0 has a minimum value of 1 (as a background contribution) and follows equation 14 when the SFR increases. Bournaud et al. (2010), Krumholz & Burkhart (2016), and Krumholz et al. (2018) showed that for galaxies with a high column density or high SFR, large-scale gravitational instabilities are significant sources of turbulent energy and might even dominate stellar feedback. The ISM in a kiloparsec region of the galaxy is far from an isolated environment. It interacts with the rest of the galaxy. Adding turbulent driving to simulations of restricted regions of a galaxy enables us to take a part of these complex interactions into account while keeping the computational cost of the simulation reasonable. The turbulence added in the kiloparsec box is clearly not fully self-consistent. We numerically investigated the effect of turbulent driving on star formation. We used a model for turbulent driving that we adapted from the generalisation of Ornstein-Uhlenbeck that was developed and used by several authors (Eswaran & Pope 1988;Schmidt et al. 2006Schmidt et al. , 2009Federrath et al. 2010).

Injection of turbulence
The turbulent forcing was implemented by adding an external force density f that accelerated the fluid on large scales. Our main hypothesis is that this force is generated by large-scale processes at the galactic scale within the galactic disk. As a consequence, the driving is bidimensional (the force has no vertical component and no vertical mode), and it was applied only at low altitude. The driving procedure is slightly different from the procedure used in Paper I and is therefore described below at length.
The fluid equations that we solved were In these equations, ρ is the gas density, v is the gas speed, P is the thermal pressure, e is the internal energy, B is the magnetic field, g is the gravitational potential, and L is the energy loss function. The evolution of the Fourier modes of the turbulence driving the acceleration fieldf (k, t) follows a stochastic differential equation (Schmidt et al. 2006(Schmidt et al. , 2009, In this equation, dt is the time step for integration, and T is the autocorrelation timescale. In our simulations, we chose T = 5 Myr and dt/T = 1/100. The value of 5 Myr was chosen to be of the same order of magnitude as the crossing time of the over-densities created by the turbulent driving. dW t is a small vector that was randomly chosen following the Wiener process, as described in Schmidt et al. (2009). The main modification to the Ornstein-Uhlenbeck process is the weighting function of the driving modes F 0 and the projection operator P ζ . We selected only large 2D modes, and more weight was given to modes with a wavelength of 500 pc (half the box size), In other words, the driving modes (divided by 2π) were and the respective weights were approximately 0.31, 0.65, 0.65, 0.94, 0.94, 1, 1. Table 1. List of simulations. Each group of simulations corresponds to an experiment discussed in the article. The simulations of the group noturb were discussed in Paper I. The parameter n 0 is the initial midplane density and sets the initial column density Σ 0,gas (see Eqs. (9) and (10)). The strength of the driving is given by parameter f rms (see Eq. (23)). A value of f rms = 10 4 in code units corresponds to an RMS acceleration due to the driving of 1.46 km · s −1 · Myr −1 . The parameter χ = 1 − ζ, introduced in Eqs. (21) and (22), is the compressibility fraction of the turbulent driving. Finally, B 0 is the initial value of the magnetic field, which was initialised as described in Eq. (11). The projection operator P ζ is the weighted sum of the components of the Helmholtz decomposition of compressive versus solenoidal modes, projected in 2D, In other words, with P ⊥ and P ∥ the projection operators perpendicular and parallel to k, respectively (Federrath et al. 2010). The weight ζ is the solenoidal fraction of the driving. In the following, we also refer to the compressive fraction, χ = 1 − ζ.
The forcing field f (x, t) was then computed from the Fourier transform, The parameter f rms controls the total power that is injected by the turbulent force into the simulation. The g(ζ) factor is an empirical correction so that the resulting time-averaged root mean square (RMS) of the power of the Fourier modes is equal to f rms , regardless of the solenoidal fraction ζ (see the discussion in Appendix A). The attenuation function att ensures that the driving occurs only on the disk, where z t = 75 pc.
Article number, page 5 of 21 A&A proofs: manuscript no. brucy2023 There are two main differences with the driving method used in Paper I. First, the modes aligned with the x -and y -axis were not used in Paper I, leading to global motions than can be seen through the diagonal features in the lower right panel of Fig. 1 of that paper. Second, the attenuation function (24) was not used in Paper I. It is introduced here so that the driving occurs only within the scale height of the disk.

Computation of the SFR and the column density
The total mass of gas accreted by the sinks was considered to form stars and is denoted M ⋆ (t). As described in section 3.2, only a fraction of the dense gas is accreted onto the stars so that the efficiency of the conversion of dense gas into stars is strictly lower than 1. The gas in the box is continually consumed by sinks or expelled from the box through the vertical boundaries. As a consequence, the gas available for star formation decreases with time. Since we used the simulations to test the star formation relation that links the SFR with the column density of gas, and to avoid too huge variations of the gas reservoir, we restricted the analysis to a short enough time span. In detail, we started the analysis at t 3% when 3% of the gas was consumed or lost, and we stopped it when this depletion fraction reached 40% (t 40% ). The surface density of SFR was then computed via The uncertainties on this value were computed by taking the standard deviation of all values of the SFR between t 3% and t 40% , that is, The mean column density in the considered time span varies between 0.97Σ 0 and 0.6Σ 0 , where Σ 0 is the initial column density. For each simulation, the associated column density is then Σ = 0.8Σ 0 ± 0.2Σ 0 . The determination of Σ SFR and of Σ is different from the determinations used in Paper I, where the SFR was averaged over a fixed period of time and only the initial column density was considered. This means that the depletion of the gas during the time span covered by the simulation was not taken into account. More precisely, f rms is the RMS acceleration of the force density in code units. Top: Mass accreted in sinks. The SFR computation started when 3% were consumed or lost (t 3% , left triangle) and ended when 40% of gas was consumed or lost (t 40% , right triangle). Bottom: Surface density of the SFR as a function of the strength of the driving f rms . The plotted value is the accretion rate between t 3% and t 40% , and the error bars reflect the standard deviation of all the accretion rates we obtained by choosing starting and end points between these two values. The dashed line is the value of the SFR given by the SK relation (Kennicutt & Evans 2012).

List of simulations
We ran several simulations that we list in Table 1. We separated them into groups, each of which corresponds to an experiment. There are two families of experiments. The first family of experiments (strength, comp, and mag) aimed to quantify the effect of one given parameter on the SFR, and the properties of the ISM in general. For these experiments, we set all parameters except for the tested parameter to fiducial values. The tested parameters were the strength of the turbulent driving, its compressibility, and the initial magnetic field.
The goal of the second family (noturb, noturb_hB, turb, turb_hB, turb_Bvar, and turb3D) was to test whether it is possible to derive an SFR that is compatible with the SK relation in simulation boxes of 1 kpc 3 even in a dense environment under various conditions. Each of the groups contained three to seven simulations with increasing column density. We computed the SFR for each of the simulations and compared it to what is expected from the SK relation. The main parameters we varied are the strength of the turbulent driving and the intensity of the ini-tial magnetic field. For the latter, several values were considered for each column density. For the simulation with turbulent driving, we chose the driving strength by trial and error to find an acceptable SFR. The goals here were to show that it is possible to obtain an SFR that was similar to the observations when turbulent driving was added, but also to analyse the resulting velocity dispersion and gas distribution.

Effects of turbulence
In the following experiments, we consider dense regions of galaxies with an initial column density of 38 and 77 M ⊙ · pc −2 . They are referred to as medium and high column density in what follows. We kept all parameters constant except for one parameter.
We started by investigating the properties of turbulence, in particular, the strength of the driving and the fraction of compressible modes. Whereas we expect that higher compressibility leads to a higher SFR, the impact of higher Mach numbers is less straightforward. One of the difficulties that arises when trying to understand the effect of the turbulent driving on the SFR is that turbulence can both trigger and quench star formation (Mac Low & Klessen 2004). In classical models for star formation (Padoan & Nordlund 2011;Hennebelle & Chabrier 2011;, the SFR is principally set by the width of the density probabilty distribution function (PDF), the density threshold above which the gas collapses under its own gravity, and the time needed for the collapse. When the Mach number is increased, the width of the PDF as well as the density threshold also increase.

Influence of the strength of the turbulent driving
We used our numerical setup to probe the effect of the strength of turbulence (simulation set strength in Table 1) by varying the RMS acceleration of the turbulence between 1 × 10 4 to 1 × 10 5 in code units 1 for medium column densities and from 5 × 10 4 to 8 × 10 5 for a high column density.
These values were selected because they lead to SFR values that are not too far to the value inferred from the SK relation. Figure 1 portrays column density images for the medium column density simulation with several values of the turbulent forcing strength, projected along the x-axis (top panel) and z-axis (bottom panel). When the turbulent forcing is weak, the gas and the recently formed stars are organised in a few dense regions. That is to say, the gas is organised as few giant molecular complexes. As turbulent driving is increased, the density field is less strongly concentrated and the stars tend to be distributed more uniformly. We also note that their number is clearly reduced.  Figure 2 shows that with stronger turbulent driving, accretion onto the sinks is indeed slower, and the SFR decreases. Fig. 1 shows that for the highest forcing values, the anisotropy that results from the stratification is even more pronounced. This is because not only the forcing alone applies in the xy-plane, but also because the star formation and therefore the stellar feedback are reduced.
Increasing the strength of the driving naturally leads to an increase in the velocity dispersion (Fig. 3). Because the driving is 2D, projected on the plane of the disk, only the velocity disper- sion along the disk plane, defined as σ disk = σ 2 x + σ 2 y / √ 2, is increased. We would expect the vertical speed dispersion, which is assumed to be mainly due to the stellar feedback, to be reduced with stronger driving, but no such trend is clearly observed. The velocity dispersion scales linearly with the driving strength, but saturates for very strong driving. Interestingly, the velocity field is generally significantly non-isotropic, with a ratio of the velocity dispersion parallel to the disk σ disk and the vertical velocity dispersion σ z that can reach a factor of 5.
The outcome of this increased velocity dispersion is a reduction in the SFR (Figs. 2 and 3). An increase in the velocity dispersion by a factor of a few leads to a decrease in the SFR by nearly a factor 10, which is a very significant drop. This illustrates the fact that the externally driven turbulence can substantially modify the ISM evolution in this regime.

Influence of the compressive fraction on the turbulence
The dual effect of turbulence, which both triggers and quenches star formation, is particularly clear when the compressive fraction of the driving is varied. The velocity field can be divided through the Helmotz decomposition into compressible modes v c and solenoidal modes v s . Compressible modes verify ∇×v c = 0 and trigger compression and dilation of the density field. By contrast, solenoidal modes verify ∇ · v s = 0 and thus do not affect the density field. When the compressive fraction χ of the turbulent driving is changed (see Eq. (22)), the amount of energy in the different velocities modes is changed accordingly. However, it is important to note that the correspondence between the compressive fraction of the turbulent driving χ and the resulting ratio of the compressive and solenoidal modes in the velocity field is not perfect. Rotation motions induced by solenoidal driving eventually generate compression as well (Vazquez-Semadeni et al. 1996), especially as the turbulence cascades to lower scales (Federrath 2013). With solenoidal driving, a smaller fraction of the kinetic energy is globally in compressive modes that can produce density fluctuations. This is well illustrated in the study of the PDF in purely hydrodynamic simulations (Federrath et al. 2008). The density PDF is broader and explores higher densities when the compressive fraction χ is closer to one (compressive driving) than when it is closer to zero (solenoidal driving). To quantify the effect of the compressibility of the turbulent driving on the SFR, we ran the comp experiment. We again selected galaxies with an initial column density of 38 and 77 M ⊙ · pc −2 . With every other parameter fixed, we modified the compressive fraction χ from 0 to 1.
The column density maps of snapshots at 58 Myr of these simulations are shown in the lower panel of Fig. 1. From left to right, the compressibility is increased from 0 to 1. In the edgeon view (top), the purely solenoidal case produces horizontal filaments, while more clumpy structures are produced when the compressibility increases. The face-on view shows that dense structures tend to be smaller when the compressive fraction is smaller. These larger and denser structures are actively forming stars. As a result, and this is depicted in Fig. 4, increasing the compressibility increases the SFR by a factor of about 10. However, even a purely compressive run produces fewer stars than a run without turbulent driving (compared with Fig. 9). Therefore, the effect of turbulent driving at this strength is always to quench star formation, regardless of the compressibility.
Interestingly, even though the driving was calibrated to produce an identical velocity dispersion regardless of the compressibility (see Appendix A), a dependence of the velocity dispersion on the compressibility (see Fig. 5) is nevertheless observed. The variations in the velocity dispersion remain well below the variations induced by the variations of the driving strength, however.
Overall, these results clearly indicate that the nature and thus the sources of the turbulence are important to identify, not only regarding the driving intensity of the turbulence, but also regarding the compressibility of the ISM turbulence. Both play a significant role in setting the SFR.

Influence of the magnetic field
The magnetic field plays an important role on the gas dynamics and in the structure of the ISM (Hennebelle & Inutsuka 2019). In the following, we study the impact of stronger initial magnetic fields on the SFR in our simulations when all other parameters remain the same. We considered the simulations from the group mag (cf Table 1). We varied the initial midplane magnetic field B 0 from 0 µG to 15 µG. As a reminder, the initial magnetic field is uniform in the galactic plane and follows a Gaussian distribution in the vertical direction (see Eq. (11)). The initial column density for these runs was also 33.7 (medium) or 77.4 M ⊙ · pc −2 (high).  Figure 6 shows that an increased magnetic field has a quasilinear effect on the SFR. For the two cases (medium and high column density), as B 0 increased from 0 to 7 µG, the SFR has dropped by a factor of nearly 10. Since high magnetic intensities lead to a substantial reduction in SFR, this implies that a strong magnetic field could reduce the turbulent energy that would otherwise be required to sustain a low SFR. This point is further investigated in section 6.3.
The effect of the magnetic field on the SFR is a consequence of magnetic support, which locally resists the impact of selfgravity, but also of the global magnetic support that leads to a large disk scale height (Iffrig & Hennebelle 2017) and therefore to reduced densities. To further verify this point, we investigate the effect of the magnetic intensity on the density profile in Fig. 7. For low magnetic fields, the SFR is so high that the gas in the midplane of the disk is accreted in stars or is swept up by the SNe, and this depletes the gas near the equatorial plane, as  (Kennicutt & Evans 2012). The red zone corresponds to an SFR that is three times higher and lower and aims to represent the spread of the SK relation. The orange histogram is the mass (multiplied by 100 for readability) of stars that blow up into SNe in a time span of 20 Myr. The SFR computation for Fig. 10 starts when 3% of the gas is consumed or lost (t 3% , left triangle) and ends when 40% of the gas is consumed or lost (t 40% , right triangle). Kennicutt+2012 (Normal/irregular) Kim+2020 (R models) TURB (driven turbulence) NOTURB (self regulated) Fig. 10. Star formation relation for the groups noturb and turb. The SFR and the associated uncertainties are computed as explained in section 3.5. We added data from observations of normal and irregular galaxies, and the blue line represents the SK relation (Kennicutt & Evans 2012). We also compare our results with the results from the R models of Kim et al. (2020).
we show in the column density map and in the upper left panel, which features the averaged density profile.
The scale height, denoted with plus marks, increases when we increase the initial magnetic from 3 to 15µG. However, the density profile is not smooth, and the effect of the magnetic field is better understood by comparing the profile at t ≈ 58 Myr (in color) with the initial profile (dashed black line). The gas is lifted, and a significant amount of gas is lost through outflows (upper right panel). A stronger magnetic field triggers stronger outflows. This effect is also seen in simulations reported by Girichidis et al. (2018). Magnetically driven outflows can be distinguished from the SN outflows that are observed in the low magnetic field cases because they occur earlier in the simulation. It is important to note that increased outflows also contribute to the diminution of the SFR. Less gas remains available for the formation of gas.
The bottom panel of Fig. 7 shows that the value of the initial magnetic field does not seem to affect the distribution of the gas in the horizontal direction, where turbulent motion generated by the driving dominates the dynamics of the gas. For the same reason, it has very little effect on the velocity dispersion, as we show in Fig. 8. The regulation effect of the magnetic field is thus mainly due to how it changes the vertical structure.
Because the driving is not computed self-consistently, hypothetical stabilising effects of the magnetic field analogous to what was reported by Heitsch et al. (2007) and Zamora-Avilés et al. (2018) at smaller scales would not be captured by our model. However, the turbulent driving aims to model turbulence that is generated by instabilities at the galactic scale. Even if the magnetic field intensity is locally high, the averaged magnetic field at the galactic scale is probably not strong enough to have a significant influence on the development of the instability.

Reproducing the Schmidt-Kennicutt relation
6.1. Self-regulation by stellar feedback?
In Paper I we found that stellar feedback alone cannot regulate the SFR to match the observed values, especially when the column density is high. In this section, we focus on the simulation of the noturb group, without injection of turbulence from the large scale.
Self-regulation works well for galaxies with column densities similar to the mean value of the Milky Way (≈ 10 M ⊙ ·pc −2 ). In our simulation with a starting column density of approximately 12 M ⊙ · pc −2 (n 0 = 1), star formation is moderate and well in line with the observations. However, for higher column densities, Fig. 9 reveals that the SFR, that is, the slope of the blue solid line, is far above the observationally determined values represented by the slope of the thick red line. Moreover, the star formation versus column density relation that can be derived is much steeper (see Fig. 10), with typicallyΣ ∝ Σ 2.5 . Interestingly, a similar result was found by Kim et al. (2020). However, this comparision should be taken carefully because the simulations from Kim et al. (2020) have a different box size and shape, a different integration time, a different feedback implementation, and a higher external potential for high column density. Figure 9 shows that high SFR generates a peak in SN mass. After this peak, star formation is reduced for a while. However, the accumulation of gas at the large scale is not destroyed, and the very dense gas reservoir is replenished over a period of a few tens of million years. As a consequence, star formation can continue at a fast rate. This adds up to the fact that was already highlighted in previous studies (Brucy et al. 2020;Nusser & Silk 2022), that the energy injection via stellar feedback is lower than what would be needed to maintain the SFR at a rate compatible with the SK relation, as well as the energy injected via gravoturbulence to maintain the Toomre parameter Q ∼ 1. Therefore, this constitutes strong evidence that stellar feedback cannot regulate star formation sufficiently in high column density environments. is computed in the whole simulation box and averaged between 4 and 60 Myr. In these simulations, the driving strength was calibrated so that the SFR was close to the SK relation (see section 6.2 and 6.3, as well as Figs. 10 and 11). The mean value is the total energy injected during the period divided by the time span, while the error bar is the standard deviation from the power computed at each time step. The solid lines refer to the estimates of the possible source of turbulent driving, and the notation is the same as in section 2. The efficiency coefficients (ϵ sn ∼ 0.05, ϵ acc ∼ 0.2, ϵ rot ∈ [0.1, 1]) are reasonable estimates, as discussed in section 2 as well. To draw the blue line, we assumed that the velocity dispersion was proportional to the column density.

Impact of turbulent driving on the Schmidt-Kennicutt relation
Paper I has shown that it was possible to obtain a star formation relation that was compatible with the SK relation when the turbulent driving was scaled with the initial column density. In this work, we reproduced this result with a slightly different method for turbulent injection that we presented in section 3.4 and with the improved SFR computation presented in section 3.5. For the group of simulations turb, we kept the same magnetic field for all the simulations, and we studied the increasing column densities. The compressive fraction was set at 0.25. This choice can be discussed: it is lower than the natural mix, which would be 0.5 (Federrath et al. 2008), but the velocity field in the galactic scale simulation of Jin et al. (2017) indicates that the solenoidal mode may dominate at large scales. For each column density, we ran several simulations with different driving strengths (as in section 4.1) and selected the simulation in which the SFR was more similar to the SK relation. Figure 10 shows the result, where the differences between the SFR from simulations with and without driving are rather clear.
This proves that the turbulent driving from the large scale is a valid candidate for closing the gap between the observations and the simulations with stellar feedback alone. The question clearly is whether the added driving is realistic in terms of in-jected energy and resulting speed dispersion, and also regarding the consequence on the structure of the ISM.
One possible diagnostic is to compute the energy injected within the simulation, as has been done in Paper I. Another diagnostic is to consider the velocity dispersion. Both are discussed below.

Impact of the magnetic field on the Schmidt-Kennicutt relation
The actual value of the magnetic field in galaxies (Beck 2015;Han 2017) remains controversial because it is very difficult to measure. However, based on energy equipartition, it is reasonable to think that gas-rich galaxies are subject to stronger magnetic fields. As shown in section 5, a stronger magnetic field results in a lower SFR. It is thus interesting to investigate how the magnetic field, in conjunction with turbulence, may contribute to shaping the SK relation. To investigate this, we ran three series of simulations called noturb_hB, turb_hB, and turb_Bvar. The description of all simulations can be found in Table 1. We recall that noturb_hB and turb_hB are the same as noturb and turb, respectively, but with doubled magnetic intensities. For the group turb_Bvar, the initial magnetic field, parametrised by its midplane value B 0 , scales as the square root of the column density. In this way, the ratio of the kinetic and magnetic energy at the beginning of the simulation remains the same. Figure 11 displays the SFR for the runs . In spite of the initial field, which is twice as intense, the trend remains similar to the SFR of the noturb runs. In particular, the SFR dependence on the column density is clearly too steep.
To investigate the impact that magnetic fields have on the SK relation in conjunction with turbulence, we proceeded as with the turb group. For each value of the initial column density and magnetisation, we made an educated guess of what would be a good strength of the turbulent RMS acceleration to obtain an SFR similar to the SK (equal to within a factor of 2). We then launched the simulation, computed the SFR, and corrected the RMS acceleration if necessary. The simulations we finally selected are described in Table 1 . They produce a mean SFR that is more similar to the SK value. The goal is to determine for an initial given magnetic field intensity and column density how much turbulent energy is needed to reproduce SK.
First, Fig. 11 shows that with our choices of RMS acceleration, we reproduce the SK relation for all groups, as expected. The SFR is globally slightly lower for the turb_Bvar group, meaning that we slightly overestimated the turbulent energy required for this group.
Second, Fig. 12 shows how much power is injected into the simulation by the driven turbulence for all groups as a function of the column density. The injected power was computed at each step by computing the difference between the kinetic energy before and after applying the turbulent driving and dividing by the time step. This reveals that the required power is significantly reduced in the presence of substantial magnetic fields. Compared to the series of runs turb, the power with which it is necessary to drive the turbulence is more than a factor of ten lower when the magnetic field is doubled, and the slope is roughly the same (group turb_hB). When the magnetic field scales as the square root of the column density (group turb_Bvar), the power-law slope of the power needed to reproduce the SK relation versus the column density is lowered from almost 5 to 4. This means that magnetic fields may significantly contribute to the origin of the SK relation. The values of the magnetic field we considered are likely reasonable, but it is entirely possible that stronger intensities have to be considered, which would further reduce the required turbulent driving. Figure 12 also shows the various estimates for the energy injection made in section 2. In gas-rich galaxies, rotation could in principle provide enough energy up to Σ ≃ 100 M ⊙ pc −2 . For higher values, however, this does not seem to be the case, meaning that a source of energy that could explain the SK relation for column densities higher than 100 M ⊙ pc −2 is currently lacking.

Velocity dispersion
The velocity dispersion is also a clear and simple diagnostic to assess the realism of simulations. Figure 13 displays the massweighted velocity dispersion computed at the scale of the full kiloparsec box as a function of column density for the various runs we performed.
Several trends are worth noticing. First, in the absence of turbulent driving (runs noturb and noturb_hB), the velocity dispersion remains broadly isotropic while in the presence of driving, particularly at high column densities, the velocity dispersion in the disk plane is several times higher than the velocity along the z-axis. Interestingly, in the presence of driving, the velocity dispersion in the disk plane is higher than without driving, whereas the reverse is true for the z-velocity. This is because in the latter case, self-gravity is strong and triggers isotropic collapse in various places. In the presence of strong magnetic fields (runs turb_hB and turb_Bvar), the velocity dispersion necessary to reproduce the SK relations is lower by a factor of about 2. Since the energy dissipation rate is expected to be proportional to σ 3 , this agrees well with the factor of about ten induced for the driving power presented in Fig. 12.
It is important to stress that the estimated injected power of Fig. 12 and the velocity dispersion of Fig. 13 are not predictions, but an estimate of the strength of the turbulence needed to quench star formation down to the SK relation. The way in which the turbulence is injected and the parameters of the model, such as the strength of the magnetic field, may change these requirement estimates (see sections 6.5 and 6.3, respectively).
Comparisons with observations are difficult because postprocessing steps beyond the scope of the work are needed to ensure that the same quantities are compared in terms of scale and tracked gas. Data from the ShiZELS galaxies (Swinbank et al. 2012), where properties also measured at scales close to 1 kpc in the Hα band, feature regions with a surface density of SFR Σ sfr around 0.1 M ⊙ ·yr −1 kpc −2 and measured velocity dispersion that can reach up to 200 km · s −1 . Similar results were found by Law et al. (2009) and in the WiggleZ survey (Wisnioski et al. 2011). The values from the turb group, with a initial magnetic field at 3.8 µG, are quite extreme, indicating that for these column densities, the magnetic field is probably stronger.

Three-dimensional driving
The high anisotropy of the velocity field found in the driven simulations is a consequence of the 2D turbulent driving. The driving force is parallel to the disk. Furthermore, it does not depend on the vertical coordinate because the wave vectors of the Fourier decomposition of the force are also aligned with the disk. Both of these assumptions are quite strong and can be relaxed.
Even though large-scale motions are expected to be aligned with the galactic plane, turbulent eddies can spread in the vertical direction if their size is below the scale height of the disk. In Article number, page 13 of 21  this last part, we try to modify our driving model to take this into account. The driving force follows the same equations as in section 3.4 with the following modifications. We changed the power spectrum F 0 to add vertical modes with a wavelength of about 32 pc (the disk scale height is initially 150 pc), The projection operator P ζ was the normal 3D projection operator defined in Schmidt et al. (2009). However, we applied a factor of 0.3 to the vertical motion of the resulting force to avoid motions that were too strong in the vertical direction that may destroy the disk.
We were able to find the value of the turbulent driving for the 2D driving that leads to SFR compatible with the SK relation (Fig 14, left panel). Compared with the most similar set of simulation, turb_hB, the energy needed to achieve SK is slightly increased, as we show in Fig. 12. However, Fig. 13 shows that the velocity dispersion in the disk is lower, as is the anisotropy. The vertical velocity dispersion does not increase significantly, probably because of the strong attenuation of the vertical component of the force. With this experiment, a less anisotropic velocity dispersion can be produced by releasing the assumption on 2D turbulence. When more accurate constraints on the turbulent driving are available, the presence and strength of a vertical component are an important item to monitor.

Shear and Coriolis effect
The fiducial model includes neither shear nor the Coriolis effect because the driven turbulence module and the shearing box used in aprevious study ( Colling et al. (2018)) require incompatible boundary conditions. Our 2D and mainly solenoidal turbulent driving reproduces some features of the shear and the Coriolis effect, but in a less self-consistent way. In its current state, it does not allow us to distinguish between the effect of large-scale unordered turbulence and ordered rotation.
The effect of shear and the Coriolis effect were studied by Colling et al. (2018) and Kim et al. (2020). Both studies showed that a shear value comparable to it did not reduce the SFR enough to match the SK relation and just slightly improved the discrepancy observed in group turb in Fig. 10. Colling et al. (2018) showed that for an initial column density of 19.4 M ⊙ · pc −2 , the shear associated with an angular speed of 28 km · s −1 · kpc −1 , which is typical of the solar neighbourhood, only reduced the SFR by a factor 1.2. This is almost negligible given the high uncertainties of the type of computations. On the other hand, a higher angular speed of 56 km · s −1 · kpc −1 , which may be found in denser environments at a galactocentric radius around 4 kpc, leads to a reduction of the SFR of a factor 3.6. The turbulent driving we chose for the same initial column density of 19.4 M ⊙ · pc −2 reduced the SFR by a factor 4.5 and led to an SFR below the SK relation. This means that a region with this moderate column density undergoing a strong shear would not need additional driving. The question that arises is whether the shear may reduce the much higher energy needed to reproduce the SK law (see Fig. 12) for the cases with a high column density. While we cannot answer this question for the shear with our current setup, it is possible to answer this for the Coriolis effect alone, as we did in Appendix D. The conclusion is that for a high column density of 77 M ⊙ · pc −2 , the Coriolis effect has a limited impact on the SFR and cannot replace or complete the large-scale turbulent driving to match the SK relation.

Nature of the turbulence
Following this line of thought, an important point to keep in mind is that the turbulence injected in groups turb, turb_Bvar, turb_hB, and turb3D is calibrated so that the resulting SFR stays close to the observed value at a given column density and is not generated self-consistently. We showed in section 4 that the exact nature of the turbulence matters greatly for the SFR. The choice of the Ornstein-Uhlenbeck process makes it hard to model a very coherent driving. These drivings may produce velocity fields that are even more solenoidal that the field yielded by a purely solenoidal Ornstein-Uhlenbeck driving (e.g. with the magneto-rotational instability in Gong et al. 2020, Fig. 2). To proceed with more precise studies of gas-rich galaxies, a better modelling of the driving source is required.

Sources of feedback
Several sources of stellar feedback are not included in the model, such as stellar winds and cosmic rays. Both have been studied in other works (Gatto et al. (2017) and Rathjen et al. (2021), respectively). Stellar feedback effects do not add up linearly for the SFR: including at least one kind of early feedback (photoionisation, stellar winds, or infrared irradiation) may change the result dramatically, but combining them does not change the SFR by more than a factor of three.

Integration time
Another possible caveat is that unlike other similar studies (e.g. Kim & Ostriker (2017)), the time span over which the simulations were run was only 150 -200 Myr. This ensured that the total mass of gas within the box did not vary by more than 40% because of outflows and star formation. Previous studies have often solved the issue of outflowing gas by using very elongated boxes that kept a large fraction of the outflowing gas within the computational domain, maybe artificially.

Numerical parameters
Finally, as always with numerical simulations, the influence of the choice of the numerical parameters may be a source of concern. This includes the influence of the sink creation criterion, which was discussed in Colling et al. (2018). They showed in their figure 11b that the SFR is unchanged even when they modified the threshold for sink accretion by a factor 16, but the addition of a test to verify the boundedness of the gas was able to increase the SFR by almost a factor of two. We discuss the details of the SN recipe in Appendix B. The limitation we applied for numerical reasons to the SN subgrid model altered the properties of the low-density gas, but had a limited impact on the dense gas and the SFR. We present the convergence study with resolution in Appendix C and show that an increase in the resolution only modifies the SFR by about 10%.

Discussion and conclusions
We extended the study of Paper I of what sets the SFR at the kiloparsec scale. Using a very similar numerical setup of a kiloparsec cube region of a galaxy, we explored the parameter space, and more particularly, the effect of the compressibility of largescale turbulence and the strength of the magnetic field. In Paper I we showed that stellar feedback alone was not sufficient to quench star formation to match the observed SFRs, and that large-scale turbulent driving may be the main effect that is missing, explaining the gap. To refine the study, both the mechanism of turbulence injection (section 3.4) and the computation of the SFR (section 3.5) were improved to allow a better comparison with the observed SFR. Compared to the reworked method for turbulence injection with the method used in Paper I, quenching via large-scale turbulence is more efficient in low column density regions and slightly less efficient in high column density regions.
The main results of this study are the following: 1. In high column density environments, stellar feedback is unable to destroy a large accumulation of dense gas and quench star formation efficiently enough to yield a rate that is compatible with observations (section 6.1). 2. If it is strong enough, 2D large-scale turbulent driving can destroy these giant clouds and can efficiently quench star formation down to levels that are compatible with observations (section 6.2). 3. Increasing the 2D turbulent driving force linearly increases the velocity dispersion parallel to the disk up to a maximum limit at which it saturates. Below this saturation limit, stronger turbulent driving translates into a lower SFR (section 4.1). In the case of a low magnetic field, the velocity dispersion associated with the turbulent-driving strength that is required to match the SK relation is high compared with observations (section 6.4). 4. The compressibility of the turbulent driving force matters as more compressive driving schemes are ten times less efficient than solenoidal ones in quenching star formation (section 4.2). 5. An increased magnetic field can also have a dramatic effect on the SFR, where an increase of roughly 10 µG reduces the SFR by a factor of ten (section 5). It cannot explain the observed SFR with stellar feedback alone, but it reduces the amount of energy needed from large-scale driving (section 6.3) and the associated velocity dispersion (section 6.4).
With this work, we provided an overview of how the largescale driving and the mean value of the magnetic field can influence the SFR. The next step towards comprehending what regulates the SFR at the kiloparsec scale in general and the role of the turbulent driving in particular is to obtain more precise constraints on the energy that can be tapped from the large scale, as well as the nature of the driving.

Software
We made use of the following software and analysis tools: GNU/Linux, Ramses (Teyssier 2002), Python, Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), Pymses (Guillet et al. 2013, Astrophysix, Astropy (Astropy Collaboration et al. 2013. Many thanks to the authors for making them publicly available. Special thanks to all the contributors of ramses, and in particular to Andrew McLeod for the introduction of the turbulence driving module that has been heavily used in this work.

Data availability
The data underlying this article are available in the Galactica Database at http://www.galactica-simulations.eu, and can be accessed with the unique identifier ISMFEED. Additional data and the source code used to run simulations and perform analysis will be shared on reasonable request to the corresponding author.
The time-averaged RMS depends on the solenoidal fraction (see Fig. A.1), and the role of the normalization factor g ζ of Eq. (23) is to correct for this. Because we used 2D driving instead of 3D driving, we cannot use the normalisation factor given by Schmidt et al. (2009) in their Equation (9). Instead, we empirically computed the relevant factor by running the turbulent generation code without any normalisation and compared the measured RMS in the Fourier space with the input f rms (see Fig. A.1). A fit of the obtained curve then gives the g ζ factor we used in the code,

Appendix B: Limitation of the supernova feedback
When SNe explode in a low-density environment, they can generate very high velocities and temperatures. This leads to a numerical issue because it implies very short time-steps, to a point that the simulation can hardly progress. In order to avoid this in our simulation, the energy from SNe was mainly injected as a momentum, and both the gas temperature and the gas velocity cannot exceed T sat = 10 6 K and V sat = 300 km · s −1 . As a consequence, we did not model the very hot phase of the gas, which, despite representing only a small fraction of the total mass of the ISM, fills the majority of the volume. However, the hot gas cools down very slowly and thus cannot collapse to form stars at the timescale we studied. Furthermore, with our 4 pc resolution, the cooling radius is not resolved for most of the SNe, especially for the most massive simulations.
Limiting the velocity means that overall, the momentum injected is lower than prescribed by the SN    in Fig. B.1. The figure depicts the ratio of the total momentum effectively injected by SNe over the total expected momentum, which is the number of SNe multiplied by the reference value of 4 × 10 43 g · cm · s −1 . On average, 30 to 60% of the reference momentum is injected, with a higher ratio in a dense environment. This may seem very low, but it is worth noting that SNe exploding in a low-density environment have a quasi-negligible im- pact on dense gas, as was demonstrated by Iffrig & Hennebelle (2015).
In order to verify that the limitation of the velocity of the gas does not change the outcome of the SFR significantly, we ran two simulations analogous to noturb with n 0 = 1.5 and 6 without the limitation. When the SNe ignited, the time step of the simulation without the limited velocity was ten times lower than the fiducial time step. Therefore, the simulations without a limited velocity evolved for a shorter time. Since our main focus is the computation SFR, we compared the amount of gas accreted by the sink particles in both cases in Fig. B.2. It shows that for the period we considered, the limitation of the velocity due to SN feedback leads to an overestimation of the SFR by a factor 1.6. This is comparable to the uncertainty due to other factors, such as specific choices of the sink creation recipe, implementation of other stellar feedback, or resolution. We were unable to run the simulation longer, but similar simulations run by Kim & Ostriker (2017) and Ostriker & Kim (2022) with an elongated box in the vertical direction indicates that the SFR is bound to oscillate, with bursts of similar or slightly lower amplitude as the first event of star formation depicted in Fig. B.2.
The column density map (Fig. B.3) reveals that our fiducial simulation fails to capture an important feature of the gas structure: superbubbles of hot gas. However, the dense gas does not seem to be strongly affected by the change in the SN feedback recipe: it is displaced, but its internal structure and its quantity are roughly the same. The latter point can be verified through the mass-weighted probability distribution function of the gas density shown in Fig. B. This figure shows that the amount of very low density gas increases when the SN speed is not limited, but the amount of dense gas stay the same.

Appendix C: Convergence with the resolution
All the simulations we discussed were run with a uniform grid resolution of about 4 pc, with 256 3 cells. This allowed the wide parameter study presented here, with 63 simulations discussed out of a total of more than 200. It is important, however, to determine the influence of the resolution and whether convergence is reached. To this end, we selected two simulations (group turb with n 0 = 6 and group noturb with n 0 = 6) and reran them with a doubled resolution (512 3 cells of size ∼ 2 pc) and half resolution (128 3 cells of size ∼ 8 pc). The sink creation threshold n sink was adapted so that n sink ∝ 1 |∆x| 2 , where ∆x is the size of a cell (Kim & Ostriker 2017). Figure C.1 shows the mass accreted by sinks in these simulations. In the case without turbulent driving (top), increasing the resolution results in a increased SFR. The case with driving is less clear, especially because star formation is more stochastic in this configuration. In both cases, the variations in the amount of mass accreted by the sinks and the resulting SFR (Fig. C.2) due to the resolution are about 10%. This is small compared with the variations measured when changing the physical parameters (sections 4 and 5). It is within the same range as the other source of uncertainties, such as the choice of other numerical parameters (sink creation threshold and seed for the turbulence) or temporal variations. This convergence result is in line with resolution studies that were carried out for a similar setup at a lower column density ( (Colling et al. 2018, Fig. 9) and (Kim & Ostriker 2017, Fig. 14)). Figure C.3 shows that similarly, the variation in the massweighted velocity dispersion are small between our fiducial resolution of 4 pc and a doubled resolution of 2 pc.  Appendix D: Impact of the Coriolis effect on the SFR As stated in section 7, our model does not include shear, mainly for numerical reasons. However, it is rather easy to implement the Coriolis effect in our simulations and assess how much it impacts the SFR. To make this experiment, we added a source term ρf cor in the momentum equation (16). We have where e z is a unit vector in the vertical direction, and Ω is the angular speed, a constant and uniform parameter in our setup. In practice, we updated the velocity field after each hydro and gravity step using an implicit Crank Nicholson scheme, v x,new = v x,old 1 − (Ωdt) 2 + 2v y,old Ωdt The implementation was similar to that used in Colling et al. (2018) and was dutifully tested. We chose an initial column density of Σ 0 = 77.4 M ⊙ · pc −2 and an initial magnetic intensity B 0 = 3.8 µG. We considered values of the angular speed of 27.5 and 110 km·s −1 ·kpc −1 . When we assume a constant rotation velocity of 220 km · s −1 , these values would correspond to the angular speed at 8 and 2 kpc, respectively. As a first step, we ran the simulations without additional large-scale turbulent driving 2 to determine whether the Coriolis effect constituted an alternative mechanism to the largescale driving. Fig. D.1 (top) shows that the Coriolis effect slightly stabilises the disk, with a reduction of the SFR of nearly a factor 1.5 with Ω = 110 km · s −1 · kpc −1 compared to the simulation without the Coriolis effect. However, the resulting SFR is still much higher than the observed values. This means that even with a very strong Coriolis effect, an additional source of driving is still needed. This result is slightly surprising. Because the angular speed Ω is constant and uniform in our simulation domain, the implied Toomre parameter is given by If it is computed at the box scale (1 kpc), the velocity dispersion σ = σ 2 x + σ 2 y + σ 2 z / √ 3 is about10 km · s −1 at 20 Myr for the two runs. For the chosen values of the angular speed Ω = 27.5 and 110 km · s −1 · kpc −1 , we obtain Q = 0.5 and 2.1, respectively. We would expect the disk to be stabilised against gravitational instabilities for the simulation with Ω = 110 km · s −1 · kpc −1 , which has a value of the Toomre parameter Q above one, the stability limit. However, because of the density fluctuations and because the velocity dispersion is lower when it is computed at a smaller scale, it is possible to reach lower values of Q locally. Fig. D.2 shows the value of Q when it is computed in squares of 100 pc × 100 pc at t = 20Myr. For all the tested values of the angular speed Ω, at least one region has a Toomre parameter well below 1. It is also important to have in mind that the value of the velocity dispersion σ used in Eq. (D.5) also encompasses compressive motions, which tend to enhance star formation.
2 Self-consistent SN driving is still present.
Because of the slight stabilisation effect observed, we may also wonder whether the additional turbulent energy required to reproduce the SK relation discussed in section 6 is lowered when the Coriolis force is added. To test this, we ran the same simulations, but with a turbulent driving f rms of 10 5 (code units) with a compressibility χ = 0.25. We compared them with the run with the same parameters in the group strength, which has no Coriolis effect and yields an SFR that is ten times higher that the SFR predicted by the SK relation. Fig. D.1 (bottom) shows that the Coriolis effect has almost no impact on the SFR when it is combined with the large-scale turbulent driving. Despite a slightly different star formation history, the slope of the curve is indeed roughly similar, regardless of the angular speed, without any systematic trend. As a consequence, we find that the additional turbulent energy required to reproduce the SK relation is not affected by the presence or absence of the Coriolis effect. A more complete study of the exact role of the Coriolis force and why in this case it does not seem to play a major role in the stability of the disk is beyond the scope of this article, but can be the topic of further research work. and with turbulent driving (bottom) at t = 20 Myr and for different angular speeds Ω (in km · s −1 · kpc −1 ).