© The Authors 2023

1 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 sub-parsec 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. 2011, 2013; Hennebelle & 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 models 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 (Bournaud et al. 2010; Renaud et al. 2012; Goldbaum et al. 2015, 2016; Krumholz 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, 2010; Daddi 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 Sect. 2. Next, we present our numerical setup in Sect. 3, and we explore the effects of the turbulence and the magnetic field in Sects. 4 and 5, respectively. We then discuss the main goal of this study, that is, the reproduction of the SK relation in Sect. 6. The caveats of this study are discussed in Sect. 7. Finally, we discuss the results and conclude in Sect. 8.

2 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⁻² and a velocity dispersion σ ≃ 10–100 km s⁻¹. Since the mass contained in a surface L² is ΣL², whereas the specific kinetic energy is σ²/2 and the crossing time is L/σ, we obtain $$\begin{array}{l}{P}_{\text{diss}}\simeq \text{Σ}{L}^2\frac{{\sigma }^3}{2L},\hfill \\ \simeq 3\times {10}^{40}\text{erg\hspace{0.17em}}{\text{s}}^{-1}\left(\frac{\text{Σ}}{100M_{\odot }{\text{pc}}^{-2}}\right){\left(\frac{\sigma }{100\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^3\left(\frac{L}{1\text{\hspace{0.17em}kpc}}\right).\hfill \end{array}$$(1)

The largest amount of energy that can be extracted per unit of time from galactic rotation is typically equal to the rotation energy $V_{\text{rot}}^2/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 $$\begin{array}{l}{P}_{\text{rot}}\simeq \text{Σ}{L}^2\frac{V_{\text{rot}}^3}{4\pi R_{\text{gal}}},\hfill \\ \simeq 2\times {10}^{40}\text{erg\hspace{0.17em}}{\text{s}}^{-1}\hfill \\ \left(\frac{\text{Σ}+{\text{Σ}}_{\star }}{100M_{\odot }{\text{pc}}^{-2}}\right){\left(\frac{V_{\text{rot}}}{200\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^3{\left(\frac{L}{1\text{\hspace{0.17em}kpc}}\right)}^2{\left(\frac{R_{\text{gal}}}{2\text{\hspace{0.17em}kpc}}\right)}^{-1}\hfill \end{array}$$(2)

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), $$Q=\frac{\sigma \kappa }{\pi G\text{Σ}},$$(3)

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 υ_inf not only brings mass into the galaxy, but also delivers energy at a rate ${\dot{M}}_{\text{accr}}{\upsilon }_{\inf }^2/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, $${\dot{M}}_{\text{accr}}\simeq L^2{\text{Σ}}_{\text{SFR}}.$$(4)

We estimate the surface density of the SFR by assuming that the SK relation (Kennicutt & Evans 2012) is verified, $${\text{Σ}}_{\text{SFR}}\simeq 0.1M_{\odot }{\text{yr}}^{-1}{\text{kpc}}^{-2}{\left(\frac{\text{Σ}}{100M_{\odot }{\text{pc}}^{-2}}\right)}^{1.4}.$$(5)

On the other-hand, we compute the infall velocity from the gravitational potential, $$\frac{{\upsilon }_{\inf }^2}{2}\simeq G\left(\text{Σ}+{\text{Σ}}_{\star }+{\text{Σ}}_{\text{dm}}\right)R_{\text{gal}},$$(6)

where Σ_dm is the column density of the dark matter. We thus find that $$\begin{array}{l}{P}_{\text{acc}}\simeq {\text{Σ}}_{\text{SFR}}{L}^{2}G\left(\text{Σ}+{\text{Σ}}_{\star }+{\text{Σ}}_{\text{dm}}\right)R_{\text{gal}}\hfill \\ \simeq 5\times {10}^{37}\text{erg\hspace{0.17em}}{\text{s}}^{-1}\hfill \\ \left(\frac{\text{Σ}+{\text{Σ}}_{\star }+{\text{Σ}}_{\text{dm}}}{100M_{\odot }{\text{pc}}^{-2}}\right){\left(\frac{\text{Σ}}{100M_{\odot }{\text{pc}}^{-2}}\right)}^{1.4}{\left(\frac{L}{1\text{\hspace{0.17em}kpc}}\right)}^2\left(\frac{R_{\text{gal}}}{2\text{\hspace{0.17em}kpc}}\right).\hfill \end{array}$$(7)

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 (Σ_SFRL²)/120 M_⊙, as given by the Salpeter initial mass function, multiplied by 10⁵¹ erg.

By again assuming that the SFR is in line with the SK relation (Eq. (5)), we can estimate $$\begin{array}{l}{P}_{\text{sn}}\simeq \frac{{\text{Σ}}_{\text{SER}}{L}^2}{120M_{\odot }}{10}^{51}\text{erg}\hfill \\ \simeq 2\times {10}^{40}\text{erg\hspace{0.17em}}{\text{s}}^{-1}{\left(\frac{\text{Σ}}{100M_{\odot }{\text{pc}}^{-2}}\right)}^{1.4}{\left(\frac{L}{1\text{\hspace{0.17em}kpc}}\right)}^2.\hfill \end{array}$$(8)

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).

3 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³ 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 4pc. The impact of the resolution on the results is discussed in Appendix C.

3.1 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, $$n\left(z\right)=n_0\exp \left(-\frac{1}{2}{\left(\frac{z}{z_0}\right)}^2\right),$$(9)

with n₀ a free density parameter and z₀ = 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 $${\text{Σ}}_{\text{gas},0}=\sqrt{2\pi }{m}_p{n}_0{z}_0,$$(10)

where m_p = 1.4 × 1.66 × 10⁻²⁴ 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\sqrt{n_0/1.5}\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$ was also added. It presents a Kolmogorov power spectrum and was generated using a random phase (Kolmogorov 1941). The dependence on n₀ 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, $$B_x\left(z\right)=B_0\exp \left(-\frac{1}{2}{\left(\frac{z}{z_0}\right)}^2\right),$$(11)

where B₀ is a parameter of our simulation with a value of a few microgauss (see Table 1).

3.2 Numerical models

The gas is subject to an external gravitational potential, corresponding to old stars and dark matter, of the form $$g_{\text{ext}}\left(z\right)=-\frac{a_{1}z}{\sqrt{z^2+z_0^2}}-a_{2}z,$$(12)

with a₁ = 1.42 × 10⁻³ kpc Myr⁻², a₂ = 5.49 × 10⁻⁴ Myr⁻², and z₀ = 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³ cm⁻³. The simulation with an initial column density higher than 50 M_⊙ pc⁻² used a higher threshold of 5 × 10³ cm⁻³. 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.

3.3 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 $${\tau }_{*}\left(M\right)={\tau }_0\exp \left(-a{\left(\log \left(\frac{M}{M_0}\right)\right)}^b\right),$$(13)

with τ₀ = 3.265 Myr, M₀ = 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⁻¹. 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⁵¹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⁴³ g cm s⁻¹ 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⁶ K and of the velocity to 300 km s⁻¹ 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^{\prime }}_0$ can then be written as $${G^{\prime }}_0=\frac{{\text{Σ}}_{\text{SFR}}}{{\text{Σ}}_{\text{SFR,}\odot }}=\frac{{\text{Σ}}_{\text{SFR}}}{2.5\times {10}^{-9}{M}_{\odot }{\text{pc}}^{-2}{\text{yr}}^{-1}}.$$(14)

In our model, ${G^{\prime }}_0$ has a minimum value of 1 (as a background contribution) and follows Eq. (14) when the SFR increases.

3.4 Injection of turbulence

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. 2006, 2009; Federrath et al. 2010).

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 $$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \upsilon \right)=0,$$(15) $$\rho \left(\frac{\partial \upsilon }{\partial t}+\left(\upsilon \cdot \nabla \right)\upsilon \right)=-\nabla P+\frac{\left(\nabla \times B\right)\times B}{4\pi }+\rho g+\rho f,$$(16) $$\rho \left(\frac{\partial e}{\partial t}+\left(\upsilon \cdot \nabla \right)e\right)=-P\left(\nabla \cdot \upsilon \right)+\rho f\cdot \upsilon -\rho ℒ,$$(17) $$\frac{\partial B}{\partial t}-\nabla \times \left(\upsilon \times B\right)=0.$$(18)

In these equations, ρ is the gas density, υ 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 ℒ is the energy loss function. The evolution of the Fourier modes of the turbulence driving the acceleration field $\widehat{f}$ (k, t) follows a stochastic differential equation (Schmidt et al. 2006, 2009), $$\text{d}\widehat{f}\left(k,t\right)=-\widehat{f}\left(k,t\right)\frac{\text{d}t}{T}+F_0\left(k\right)P_{\zeta }\left(k\right)\text{d}{W}_t.$$(19)

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₀ 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), $$F_0\left(k\right)=\{\begin{array}{l}1-{\left(\frac{k}{2\pi }-2\right)}^2\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}}1\text{\hspace{0.17em}<\hspace{0.17em}\hspace{0.17em}}\frac{\left|k\right|}{2\pi }<\text{\hspace{0.17em}\hspace{0.17em}}3\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}{k}_z=0\hfill \\ 0\text{\hspace{0.17em}if\hspace{0.17em}not.}\hfill \end{array}.$$(20)

In other words, the driving modes (divided by 2π) were $$\left(\begin{array}{c}2\\ 2\\ 0\end{array}\right),\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}2\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}1\\ 2\\ 0\end{array}\right),\left(\begin{array}{c}2\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 2\\ 0\end{array}\right),$$

and the respective weights were approximately

0.31, 0.65, 0.65, 0.94, 0.94, 1, 1.

The projection operator P_ζ is the weighted sum of the components of the Helmholtz decomposition of compressive versus solenoidal modes, projected in 2D, $$P_{\zeta }\left(k\right)=\left(\begin{array}{ccc}\zeta +\left(1-2\zeta \right)\left(\frac{k_x^2}{k^2}\right)& \left(1-2\zeta \right)\left(\frac{k_x{k}_y}{k^2}\right)& 0\\ \left(1-2\zeta \right)\left(\frac{k_x{k}_y}{k^2}\right)& \zeta +\left(1-2\zeta \right)\left(\frac{k_y^2}{k^2}\right)& 0\\ 0& 0& 0\end{array}\right)$$(21)

In other words, $$P_{\zeta }\left(k\right)=\zeta P^{\perp }\left(k\right)+\left(1-\zeta \right)P^{\parallel }\left(k\right),$$(22)

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, $$f\left(x,t\right)=\text{att}\left(z\right)g\left(\zeta \right)f_{\text{rms}}{\displaystyle \int \widehat{f}\left(k,t\right)e^{ik\cdot x}\text{d}{t}^{3}k.}$$(23)

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, $$\text{att}\left(z\right)=\exp \left(\frac{-z^2}{2z_t^2}\right),$$(24)

where z_t = 75 pc.

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.

3.5 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 Sect. 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 $${\text{Σ}}_{\text{SFR}}\left(t_{3\%},t_{40\%}\right)=\frac{M_{\text{⋆}}\left(t_{40\%}\right)-M_{\text{⋆}}\left(t_{3\%}\right)}{L^2\left(t_{40\%}-t_{3\%}\right)}.$$(25)

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, Σ_sfr(t_i, t_j) for t_3% ≤ t_i < t_j ≤ t_40%.

The mean column density in the considered time span varies between 0.97Σ₀ and 0.6Σ₀, where Σ₀ is the initial column density. For each simulation, the associated column density is then Σ = 0.8Σ₀ ± 0.2Σ₀. 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.

3.6 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³ 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 initial 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.

4 Effects of turbulence

In the following experiments, we consider dense regions of galaxies with an initial column density of 38 and 77 M_⊙ pc⁻². 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; Federrath & Klessen 2012), 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.

4.1 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⁴ to 1 × 10⁵ in code units¹ for medium column densities and from 5 × 10⁴ to 8 × 10⁵ 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.

Figure 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 dispersion along the disk plane, defined as ${\sigma }_{\text{disk}}=\sqrt{{\sigma }_x^2+{\sigma }_y^2}/\sqrt{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 σ_dįsk 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.

Fig. 1 Column densities for the group STRENGTH (top) and COMP (bottom) for a medium column density. For each panel, the top subpanel is viewed edge-on along the x-axis and the bottom subpanel is viewed face-on. White dots represent sink particles, with a radius proportional to their mass. Top: RMS acceleration of the turbulent driving is increased from left to right. Bottom: compressibilty of the turbulent driving is increased from left to right. χ = 0 and χ = 1 refer to purely solenoidal and compressive driving, respectively.

4.2 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 υ_c and solenoidal modes υ_s. Compressible modes verify ∇ × υ_c = 0 and trigger compression and dilation of the density field. By contrast, solenoidal modes verify ∇ · υ_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⁻². 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 edge-on 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.

Fig. 2 Effect of the intensity of the turbulent driving on the SFR (group STRENGTH). frms is a measure of the driving strength. More precisely, frms 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 (t3%, left triangle) and ended when 40% of gas was consumed or lost (t40%, right triangle). Bottom: surface density of the SFR as a function of the strength of the driving frms. The plotted value is the accretion rate between t3% and t40%, 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).

5 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₀ 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⁻² (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₀ 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 Sect. 6.3.

The effect of the magnetic field on the SFR is a consequence of magnetic support, which locally resists the impact of self-gravity, 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 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.

Fig. 3 Relation between the strength of the driving, the velocity dispersion, and the SFR. Top panel: mass-weighted velocity dispersion measured in the STRENGTH simulations at 58 Myr at the kiloparsec scale. The velocity dispersion parallel to the disk plane σdisk is defined by ${\sigma }_{\text{disk}}=\sqrt{{\sigma }_x^2+{\sigma }_y^2}/\sqrt{2}$. Bottom panel: SFR as a function of the velocity dispersion parallel to the disk. The dashed line is the value of the SFR given by the SK relation.

Fig. 4 Effect of the compressibility of the turbulent driving on the SFR (group COMP). Surface density of the SFR as a function of the compressibility χ of the driving. The dashed line is the value of the SFR given by the SK relation.

Fig. 5 Mass-weighted velocity dispersion measured in the COMP simulations at 58 Myr at the kiloparsec scale.

Fig. 6 Effect of the intensity of the initial magnetic field on the SFR (group MAG). B0 is the initial midplane intensity of the magnetic field. The legend is the same as in Fig. 2. The dashed line is the value of the SFR given by the SK relation.

Fig. 7 Effect of the magnetic field on the structure of the disk and the outflows (group MAG) for a high column density. Top left panel: averaged density profile (slightly smoothed for readability) at t ≈ 58 Myr. For runs in which the midplane density was close to the maximum density, the position at which the density is half of the midplane density is marked with a plus. The dashed black line corresponds to the initial profile. The colours of the lines refer to the same colour bar as in the plot of the top right panel. Top right panel: mass loss via outflows through the open up and down boundaries. Lower panel: edge-on (top) and face-on (bottom) column densities.

Fig. 8 Velocity dispersion measured in the MAG simulations at 58 Myr.

6 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⁻²). In our simulation with a starting column density of approximately 12 M_⊙ pc⁻² (n₀ = 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 gravo-turbulence 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.

6.2 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 Sect. 3.4 and with the improved SFR computation presented in Sect. 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 Sect. 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 injected 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.

6.3 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 Sect. 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₀, 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−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 Sect. 2. In gas-rich galaxies, rotation could in principle provide enough energy up to Σ ≃ 100 M_⊙ pc⁻². 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⁻² is currently lacking.

Fig. 9 Stellar mass for the group NOTURB. The red line is the stellar mass creation when the SFR in the kiloparsec box is equal to 2.5 × 10−3(Σ/10 M⊙ pc−2)14 M⊙ × yr−1 (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 (t3%, left triangle) and ends when 40% of the gas is consumed or lost (t40%, right triangle).

Fig. 10 Star formation relation for the groups NOTURB and TURB. The SFR and the associated uncertainties are computed as explained in Sect. 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).

Fig. 11 Star formation relation for the simulations with a high magnetic field.

Fig. 12 Power injected by the turbulent driving in the simulations compared with estimations. Markers represent the time-averaged power injected by the turbulent driving force in order to reproduce the SK relation. The power 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 Sect. 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 Sect. 2. The efficiency coefficients (ϵsn ~ 0.05, ϵacc ~ 0.2, ϵrot ϵ [0.1, 1]) are reasonable estimates, as discussed in Sect. 2 as well. To draw the blue line, we assumed that the velocity dispersion was proportional to the column density.

6.4 Velocity dispersion

The velocity dispersion is also a clear and simple diagnostic to assess the realism of simulations. Figure 13 displays the mass-weighted 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 σ³, 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 Sects. 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⁻¹kpc⁻² and measured velocity dispersion that can reach up to 200 km s⁻¹. 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.

6.5 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 this last part, we try to modify our driving model to take this into account. The driving force follows the same equations as in Sect. 3.4 with the following modifications.

We changed the power spectrum F₀ to add vertical modes with a wavelength of about 32 pc (the disk scale height is initially 150 pc), $$F_0\left(k\right)=\{\begin{array}{ll}{\left(\frac{k}{2\pi }\right)}^{-5/3}\hfill & \text{if\hspace{0.17em}}2\le \frac{k_x}{2\pi },\frac{k_y}{2\pi }\le 3\text{\hspace{0.17em}and\hspace{0.17em}}30\le \frac{k_z}{2\pi }\le 32,\hfill \\ 0\text{\hspace{0.17em}if\hspace{0.17em}not.}\hfill & \hfill \end{array}$$(26)

The projection operator Ρ_ζ 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.

Fig. 13 Mass-weighted velocity dispersion in the full computational volume of 1 kpc3 computed at t = 58 Myr. Left panel: velocity parallel to the galactic plane. Right panel: vertical velocity.

7 Caveats

7.1 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⁻², the shear associated with an angular speed of 28 km s⁻¹ kpc⁻¹, 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⁻¹ kpc⁻¹, 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⁻² 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⁻², the Coriolis effect has a limited impact on the SFR and cannot replace or complete the large-scale turbulent driving to match the SKrelation.

Fig. 14 SFR surface density for the TURB3D group compared to the TURB_HB. The legend is the same as in Fig. 10.

7.2 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 TURB_3D 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 Sect. 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.

7.3 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.

7.4 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.

7.5 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 Fig. 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%.

8 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 large-scale 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 (Sect. 3.4) and the computation of the SFR (Sect. 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 (Sect. 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 (Sect. 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 (Sect. 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 (Sect. 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 (Sect. 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 (Sect. 5). It cannot explain the observed SFR with stellar feedback alone, but it reduces the amount of energy needed from large-scale driving (Sect. 6.3) and the associated velocity dispersion (Sect. 6.4).

With this work, we provided an overview of how the large-scale 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.

Several analytical works have investigated the possible source of this turbulence, gravitational instabilities (Nusser & Silk 2022), mass transfer within the galactic disk (Krumholz et al. 2018; Meidt et al. 2020), or accretion onto the disk from the circumgalactic medium. Observations from local (e.g. Sun et al. 2022) or high-redshift galaxies, as well as isolated and cosmological simulations, are other ways to better know how energy is transferred from large to small scales.

Acknowledgements

The authors thank the referee, Enrique Vazquez-Semadeni, for his useful comments which helped to significantly improve the article. The authors acknowledge Interstellar Institute’s program “With Two Eyes” and the Paris-Saclay University’s Institut Pascal for hosting discussions that nourished the development of the ideas behind this work. N.B. thanks Jérémy Fensch, Chang Goo Kim and Munan Gong for interesting and fruitful discussions. This work was granted access to HPC resources of CINES and CCRT under the allocation x2020047023 made by GENCI (Grand Equipement National de Calcul Intensif) and the special allocation 2021SA10spe00010. N.B., P.H. and T.C. acknowledge financial support from the European Research Council (ERC) via the ERC Synergy Grant “ECOGAL: understanding our Galactic ecosystem – From the disk of the Milky Way to the formation sites of stars and planets” (grant 855130). 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), RAMSES (Guillet et al. 2013), ASTROPHYSIX, ASTROPHY (Astropy Collaboration 2013, 2018). 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 (http://www.galactica-simulations.eu/db/STAR_FORM/ISMFEED). Additional data and the source code used to run simulations and perform analysis will be shared on reasonable request to the corresponding author.

Appendix A Calibration of the turbulent driving

In section 4.2 we studied the effect of changing the solenoidal fraction ζ of the driving on the SFR. In order to have a fair comparison, it is important to ensure that the strength of the driving is the same, regardless of ζ. A way to check this is to consider the RMS power of the Fourier modes, $$\text{RMS}\left(\text{t}\right)=\sqrt{{\displaystyle \int {\left|\widehat{f}\left(k,t\right)\right|}^2{\text{d}}^{3}k}}.$$(A.1)

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, $$g_{\zeta }=0.93{\left(1-\zeta \right)}^2-0.87\left(1-\zeta \right)+\mathrm{0.76.}$$(A.2)

Fig. A.1 Ratio of the time-averaged RMS (as in Eq. (A.1)) over the frms parameter for unnormalised simulations. The fit of the curve gives the normalisation factor gζ.

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⁶ K and V_sat = 300 km s⁻¹. 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.

Fig. B.1 Averaged ratio of the injected momentum over the expected reference momentum of 4 × 1043 g · cm · s−1 for all the SNe that blew up during the simulation. For simulations with higher initial column densities, the injected momentum approaches the reference value.

Fig. B.2 Comparison of the SFR with (blue line) and without (orange line) limiting the velocity of the SN ejecta for n0 = 1.5 (top) and n0 = 6 (bottom).

Limiting the velocity means that overall, the momentum injected is lower than prescribed by the SN model, as we show 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⁴³ g · cm · s⁻¹. 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 impact on dense gas, as was demonstrated by Iffrig & Hennebelle (2015).

Fig. B.3 Comparison of the edge-on (top) and face-on (bottom) column density map with and without (left and right of each panel, respectively) limiting the velocity of the SN ejecta for n0 = 1.5 (left panel) and n0 = 6 (right panel).

Fig. B.4 Mass-weighted probability distribution function of the gas density for the same snapshots as in Fig. B.3.

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₀ = 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.4. 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³ 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₀ = 6 and group NOTURB with n₀ = 6) and reran them with a doubled resolution (512³ cells of size ~ 2 pc) and half resolution (128³ cells of size ~ 8 pc). The sink creation threshold n_sink was adapted so that $n_{\text{sink}}\propto \frac{1}{{\left|\text{Δ}x\right|}^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 mass-weighted velocity dispersion are small between our fiducial resolution of 4 pc and a doubled resolution of 2 pc.

Fig. C.1 Comparison of the total mass accreted by sinks as a function of time for different resolutions for a simulation with turbulence driving (top) and a simulation without it (bottom).

Fig. C.2 Comparison of the computed SFR for different resolutions for a simulation with turbulence driving (top) and a simulation without it (bottom). The dashed line is the value of the SFR given by the SK relation.

Fig. C.3 Comparison of the mass-weighted velocity dispersion for different resolutions for a simulation with turbulence driving at t = 72 Myr.

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 $$f_{cor}=-2\text{Ω}{e}_z\times \upsilon ,$$(D.1)

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, $${\upsilon }_{\text{x,new}}=\frac{{\upsilon }_{\text{x,old}}\left(1-{\left(\text{Ω}dt\right)}^2\right)+2{\upsilon }_{\text{y,old}}\text{Ω}dt}{{\left(\text{Ω}dt\right)}^2},$$(D.2) $${\upsilon }_{\text{y,new}}=\frac{{\upsilon }_{\text{y,old}}\left(1-{\left(\text{Ω}dt\right)}^2\right)-2{\upsilon }_{\text{x,old}}\text{Ω}dt}{{\left(\text{Ω}dt\right)}^2},$$(D.3) $${\upsilon }_{\text{z,new}}={\upsilon }_{\text{z,old}}.$$(D.4)

The implementation was similar to that used in Colling et al. (2018) and was dutifully tested.

We chose an initial column density of Σ₀ = 77.4 M_⊙ pc⁻² and an initial magnetic intensity B₀ = 3.8 µG. We considered values of the angular speed of 27.5 and 110 km s⁻¹ · kpc⁻¹. When we assume a constant rotation velocity of 220 km s⁻¹, 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² to determine whether the Coriolis effect constituted an alternative mechanism to the large-scale 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⁻¹ · kpc⁻¹ 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 $$Q=\frac{2\text{Ω}\sigma }{\pi G\text{Σ}}.$$(D.5)

If it is computed at the box scale (1 kpc), the velocity dispersion $\sigma =\sqrt{{\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2}/\sqrt{3}$ is about10 km s⁻¹ at 20 Myr for the two runs. For the chosen values of the angular speed Ω = 27.5 and 110 km s⁻¹ · kpc⁻¹, 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⁻¹ · kpc⁻¹, 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 = 20 Myr. 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.

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⁵ (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.

Fig. D.1 Comparison of the accreted mass in sink particles for different angular speeds Ω (in km s−1 · kpc−1) for the simulations without turbulent driving (top) and with turbulent driving frms of 105 in code units (bottom). The time is rescaled so that 0 Myr corresponds to the time at which 105 M⊙ were accreted into sinks.

Fig. D.2 Value of the Toomre parameter Q (Eq. (D.5)) computed in a square of 100 pc × 100 pc for the simulations without turbulent driving (top) and with turbulent driving (bottom) at t = 20 Myr and for different angular speeds Ω (in km s−1 · kpc−1).

References

All Tables

All Figures

Fig. 1 Column densities for the group STRENGTH (top) and COMP (bottom) for a medium column density. For each panel, the top subpanel is viewed edge-on along the x-axis and the bottom subpanel is viewed face-on. White dots represent sink particles, with a radius proportional to their mass. Top: RMS acceleration of the turbulent driving is increased from left to right. Bottom: compressibilty of the turbulent driving is increased from left to right. χ = 0 and χ = 1 refer to purely solenoidal and compressive driving, respectively.

In the text

Fig. 2 Effect of the intensity of the turbulent driving on the SFR (group STRENGTH). frms is a measure of the driving strength. More precisely, frms 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 (t3%, left triangle) and ended when 40% of gas was consumed or lost (t40%, right triangle). Bottom: surface density of the SFR as a function of the strength of the driving frms. The plotted value is the accretion rate between t3% and t40%, 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).

In the text

Fig. 3 Relation between the strength of the driving, the velocity dispersion, and the SFR. Top panel: mass-weighted velocity dispersion measured in the STRENGTH simulations at 58 Myr at the kiloparsec scale. The velocity dispersion parallel to the disk plane σdisk is defined by ${\sigma }_{\text{disk}}=\sqrt{{\sigma }_x^2+{\sigma }_y^2}/\sqrt{2}$. Bottom panel: SFR as a function of the velocity dispersion parallel to the disk. The dashed line is the value of the SFR given by the SK relation.

In the text

	Fig. 4 Effect of the compressibility of the turbulent driving on the SFR (group COMP). Surface density of the SFR as a function of the compressibility χ of the driving. The dashed line is the value of the SFR given by the SK relation.
In the text

	Fig. 5 Mass-weighted velocity dispersion measured in the COMP simulations at 58 Myr at the kiloparsec scale.
In the text

	Fig. 6 Effect of the intensity of the initial magnetic field on the SFR (group MAG). B0 is the initial midplane intensity of the magnetic field. The legend is the same as in Fig. 2. The dashed line is the value of the SFR given by the SK relation.
In the text

Fig. 7 Effect of the magnetic field on the structure of the disk and the outflows (group MAG) for a high column density. Top left panel: averaged density profile (slightly smoothed for readability) at t ≈ 58 Myr. For runs in which the midplane density was close to the maximum density, the position at which the density is half of the midplane density is marked with a plus. The dashed black line corresponds to the initial profile. The colours of the lines refer to the same colour bar as in the plot of the top right panel. Top right panel: mass loss via outflows through the open up and down boundaries. Lower panel: edge-on (top) and face-on (bottom) column densities.

In the text

	Fig. 8 Velocity dispersion measured in the MAG simulations at 58 Myr.
In the text

Fig. 9 Stellar mass for the group NOTURB. The red line is the stellar mass creation when the SFR in the kiloparsec box is equal to 2.5 × 10−3(Σ/10 M⊙ pc−2)14 M⊙ × yr−1 (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 (t3%, left triangle) and ends when 40% of the gas is consumed or lost (t40%, right triangle).

In the text

	Fig. 10 Star formation relation for the groups NOTURB and TURB. The SFR and the associated uncertainties are computed as explained in Sect. 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).
In the text

	Fig. 11 Star formation relation for the simulations with a high magnetic field.
In the text

Fig. 12 Power injected by the turbulent driving in the simulations compared with estimations. Markers represent the time-averaged power injected by the turbulent driving force in order to reproduce the SK relation. The power 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 Sect. 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 Sect. 2. The efficiency coefficients (ϵsn ~ 0.05, ϵacc ~ 0.2, ϵrot ϵ [0.1, 1]) are reasonable estimates, as discussed in Sect. 2 as well. To draw the blue line, we assumed that the velocity dispersion was proportional to the column density.

In the text

	Fig. 13 Mass-weighted velocity dispersion in the full computational volume of 1 kpc3 computed at t = 58 Myr. Left panel: velocity parallel to the galactic plane. Right panel: vertical velocity.
In the text

	Fig. 14 SFR surface density for the TURB3D group compared to the TURB_HB. The legend is the same as in Fig. 10.
In the text

	Fig. A.1 Ratio of the time-averaged RMS (as in Eq. (A.1)) over the frms parameter for unnormalised simulations. The fit of the curve gives the normalisation factor gζ.
In the text

	Fig. B.1 Averaged ratio of the injected momentum over the expected reference momentum of 4 × 1043 g · cm · s−1 for all the SNe that blew up during the simulation. For simulations with higher initial column densities, the injected momentum approaches the reference value.
In the text

	Fig. B.2 Comparison of the SFR with (blue line) and without (orange line) limiting the velocity of the SN ejecta for n0 = 1.5 (top) and n0 = 6 (bottom).
In the text

	Fig. B.3 Comparison of the edge-on (top) and face-on (bottom) column density map with and without (left and right of each panel, respectively) limiting the velocity of the SN ejecta for n0 = 1.5 (left panel) and n0 = 6 (right panel).
In the text

	Fig. B.4 Mass-weighted probability distribution function of the gas density for the same snapshots as in Fig. B.3.
In the text

	Fig. C.1 Comparison of the total mass accreted by sinks as a function of time for different resolutions for a simulation with turbulence driving (top) and a simulation without it (bottom).
In the text

	Fig. C.2 Comparison of the computed SFR for different resolutions for a simulation with turbulence driving (top) and a simulation without it (bottom). The dashed line is the value of the SFR given by the SK relation.
In the text

	Fig. C.3 Comparison of the mass-weighted velocity dispersion for different resolutions for a simulation with turbulence driving at t = 72 Myr.
In the text

	Fig. D.1 Comparison of the accreted mass in sink particles for different angular speeds Ω (in km s−1 · kpc−1) for the simulations without turbulent driving (top) and with turbulent driving frms of 105 in code units (bottom). The time is rescaled so that 0 Myr corresponds to the time at which 105 M⊙ were accreted into sinks.
In the text

	Fig. D.2 Value of the Toomre parameter Q (Eq. (D.5)) computed in a square of 100 pc × 100 pc for the simulations without turbulent driving (top) and with turbulent driving (bottom) at t = 20 Myr and for different angular speeds Ω (in km s−1 · kpc−1).
In the text