A new star formation recipe for magneto-hydrodynamics simulations of galaxy formation

Star formation has been observed to occur at globally low yet locally varying efficiencies. As such, accurate capture of star formation in numerical simulations requires mechanisms that can replicate both its smaller-scale variations and larger-scale properties. Magnetic fields are thought to play an essential role within the turbulent interstellar medium (ISM) and affect molecular cloud collapse. However, it remains to be fully explored how a magnetised model of star formation might influence galaxy evolution. We present a new model for a sub-grid star formation recipe that depends on the magnetic field. We run isolated disk galaxy simulations to assess its impact on the regulation of star formation using the code RAMSES. Building upon existing numerical methods, our model derives the star formation efficiency from local properties of the sub-grid magnetised ISM turbulence, assuming a constant Alfv\'en speed at sub-parsec scales. Compared to its non-magnetised counterpart, our star formation model suppresses the initial starburst by a factor of two, while regulating star formation later on to a nearly constant rate of $\sim 1~M_{\odot}~\mathrm{yr}^{-1}$. Differences also arise in the local Schmidt law with a shallower power law index for the magnetised star formation model. Our results encourage further examination into the notion that magnetic fields are likely to play a non-trivial role in our understanding of star and galaxy formation.


INTRODUCTION
Star formation in disk galaxies results from a complex interaction between various mechanisms, such as gravitational collapse, cooling, and turbulence.However, a comprehensive model of both the sum and its parts has yet to be wholly established.Star-forming regions in spiral, gas-rich galaxies generally convert available gas to stars at a surprisingly low instantaneous efficiency (Persic & Salucci 1992;Fukugita et al. 1998;Balogh et al. 2001), with numerical simulations and observational constraints suggesting a parent dark matter halo mass dependence that peaks at ∼20% efficiency for Milky Way sized objects (Behroozi et al. 2013a,b).These regions can also be characterised by the ratio of available gas mass to star formation rate, otherwise known as the depletion time scale  d (Leroy et al. 2008;Bauermeister et al. 2010).
While observed star formation rates can vary dramatically across disk galaxies, empirical correlations emerge at scales larger than 100 pc.The Schmidt (1959) law describes the relation between observed star formation rate surface density Σ SFR and gas surface density Σ gas as a simple power law, where  is an empirically determined normalisation constant.Local Schmidt laws, as measured radially within individual galaxies, demonstrate power law correlations with widely ranging slopes and normalisations (Wong & Blitz 2002).However Kennicutt (1998) first ★ E-mail: egirma@princeton.edufound that when averaged over the entire galaxy, a global Kennicutt-Schmidt relation arises and appears consistent across galaxies, following a power law with slope  ≈ 1.4.
It remains unclear what physical processes, or interplay of processes, result in these observed relations, though presently debated models invoke the role of gravitational instabilities (e.g.Quirk 1972;Larson 1988;Elmegreen 1994;Kennicutt 1989Kennicutt , 1998)), cloud-cloud collisions (e.g.Wyse 1986;Wyse & Silk 1989;Silk 1997;Tan 2000), turbulence and its influence on the gas density probability distribution function (PDF) (e.g.Scalo et al. 1998;Passot & Vázquez-Semadeni 1998;Ostriker et al. 1999;Klessen 2000;Wada & Norman 2001;Ballesteros-Paredes & Mac Low 2002;Li et al. 2003;Kravtsov 2003;Mac Low et al. 2005), and gas thermodynamics (e.g.Struck-Marcell 1991;Struck & Smith 1999).Due to the difficulty of encoding all these physics into cosmological simulations, as well as the inability to resolve individual molecular clouds or multiphase ISM, galaxy formations frequently assume a sub-grid star formation recipe based on the Schmidt law.In absence of an explicit Schmidt-like star formation prescription, however, Kennicutt-Schmidt relations have successfully emerged in more recent stratified disk models (Li et al. 2006;Bieri et al. 2023).These models notably suggest that a dynamic interplay between gravitational instabilities, supernovae (SNe) feedback, and gas cooling somehow drives star formation rates towards what is locally observed.
Schmidt-like sub-grid star formation recipes can assume a uniform and constant  ff , the value of which can be adjusted to reproduce the observed Kennicutt-Schmidt relation (see e.g.Li et al. 2006;Wada & Norman 2007).The star formation rate can then be modelled as where  gas is the local gas density and  ff is the local free-fall time (Stinson et al. 2006).However, extreme variation of  ff has been observed in individual molecular clouds of radii  MC ≲ 100 pc, with values ranging from 1% to 30% (Murray 2011).A physical model for star formation must necessarily explain and capture this smaller-scale variability, while accurately reproducing the observed larger-scale relations and the overall low efficiency in galaxies.
Attempts to identify what external and internal factors regulate star formation have long been underway.Several authors have explored supersonic turbulence in the interstellar medium (ISM) as a promising candidate, as it both counteracts gravitational collapse via turbulent pressure and induces gravitational collapse by sporadically creating dense regions (Chandrasekhar 1951;Bonazzola et al. 1987;Krumholz & McKee 2005;Hennebelle & Chabrier 2008;Federrath & Klessen 2012).The presence of turbulence in molecular clouds can be inferred from the empirical proportionality between velocity dispersion and cloud size (Larson 1981) as well as from measurements of velocity and density power spectra (Heyer & Brunt 2004).
Yet without any driving forces the turbulent energy decays far too quickly relative to the longest estimates of molecular cloud lifetimes (Blitz & Shu 1980;Stone et al. 1998;Mac Low et al. 1998).Potential mechanisms for turbulence production in molecular clouds exist on both large scales (≳  MC ) and small scales (≲  MC ).Models in which larger scale processes like SN feedback and spiral shock forcing dominate turbulence driving appear most consistent with observations (Brunt et al. 2009).Turbulence in the ISM may also be driven at smaller scales by stellar outflows (Bally et al. 1996;Knee & Sandell 2000), stellar winds (Norman & Silk 1980;Lada & Gautier 1982;Franco & Cox 1983), and compact photoionizing HII regions (Matzner 2002).
One factor that remains significant, and is known to operate on a variety of scales as both a regulator and a trigger of star formation is the magnetic field.Magnetic fields are capable of supporting molecular clouds up to a critical mass value where Φ is the magnetic flux,  is the magnetic field strength, and  is the molecular cloud radius (Mouschovias & Spitzer 1976;Shu et al. 1987).For clouds of mass  <  cr (known as "subcritical"), ambipolar diffusion can trigger fragmentation and star formation within their interior by adequately redistributing magnetic flux (Mestel & Spitzer 1956;Mouschovias 1976bMouschovias ,a, 1987)).As clouds with sufficiently low ion-neutral collisions settle into a gravitationally unstable neutral core and ionized envelope, magnetic fields regulate timescales of collapse by virtue of the ambipolar diffusion timescale, where   is the ion density,   is the neutral core density, and   is the neutral core mass (Mouschovias 1977).
However, this process requires strong initial magnetic fields to explain rapid star formation rates observed over ≲ 1 Myr (Hartmann et al. 2001).In addition, most recent observations suggest that the median molecular clouds are actually "supercritical", where magnetic fields alone cannot counteract collapse (Crutcher 2012).For these clouds, magnetised supersonic turbulence appears to decrease star formation by factors of 2-3 when compared with non-magnetised turbulent flows (Price & Bate 2009;Dib et al. 2010;Padoan & Nordlund 2011;Federrath & Klessen 2012;Padoan et al. 2012).Reasons for this seem to be at least twofold.Simulations have demonstrated that the presence of magnetic fields narrows the PDF imposed on gas density by turbulent flows, making it more difficult to reach higher densities (Cho & Lazarian 2003;Kowal et al. 2007;Burkhart et al. 2009;Molina et al. 2012;Mocz et al. 2017).Furthermore, magnetic fields change the criteria for collapse itself, as they provide additional support against self-gravity.
Past studies have added magneto-hydrodynamics (MHD) to turbulence-regulated star formation theoretical models (Krumholz & McKee 2005;Padoan & Nordlund 2011), which appear to agree well with numerical simulations and observations of molecular clouds (Federrath & Klessen 2012).Computational experiments on the evolution of turbulence in isothermal ISM, considering the observed magnetic field distribution  ∝  1/2 gas (Crutcher 1999;Crutcher et al. 2003), have additionally revealed more precise relationships between gas logarithmic density variance, sonic Mach number, and magnetic field strength (Molina et al. 2012).However, the possible effects of magnetic fields and magnetised turbulent ISM models on galaxy evolution and properties have yet to be fully explored.MHD simulations of galaxy formation have thus far demonstrated that magnetic fields amplify rapidly within the turbulent rotating disk over the first ∼ 500 Myr, before reaching a point of saturation and self-regulation.Their force then contributes to lower late-time star formation rates, relative to those in hydrodynamical disk galaxy simulations (Wang & Abel 2009;Dubois & Teyssier 2010;Pakmor & Springel 2013;Steinwandel et al. 2020).It remains an open question exactly what mechanisms drive this magnetic field evolution from its theorised weak primordial state to the ≳ G field strengths observed today (Beck 2016;Han 2017), though numerical methods have been employed to probe a range of different possible dynamos: namely the galactic -Ω dynamo (Wang & Abel 2009;Dubois & Teyssier 2010), turbulent feedback driven dynamo (Pakmor & Springel 2013;Vazza et al. 2014;Rieder & Teyssier 2016), magneto-rotational instability driven dynamo (Gressel et al. 2013;Machida et al. 2013), and cosmic-ray driven dynamo (Lesch & Hanasz 2003;Hanasz et al. 2009).One persistent challenge in galaxy simulations to successfully model both this magnetised dynamo amplification and its effects on local galactic processes, as well as global galaxy evolution, is due to the dramatic spatial range these physics involve.In the case of star formation, certain studies opt to exclude it in their simulations in order to more rigorously model magnetic fields and the entire galaxy (e.g.Wang & Abel 2009;Steinwandel et al. 2022).Others instead base their star formation routines on purely hydrodynamical methods, that then do not directly take into account local magnetic field properties (e.g.Cen & Ostriker 1992;Springel & Hernquist 2003).Katz et al. (2021) examine primordial magnetic fields (PMFs) in a cosmomlogical context through a series of multi-frequency radiation-MHD simulations, and while their findings offer constraints on how PMFs may impact galaxy formation and evolution, the large spatial scales of their simulations prevent any detailed study of star-forming regions.
To date, very few computational experiments exist that incorporate MHD-based star formation models into isolated disk galaxy simulations or simulations within a cosmological context.Martin-Alvarez et al. (2020) have extended thermo-turbulent star formation models by Trebitsch et al. (2017), Mitchell et al. (2018) and Rosdahl et al. (2017) into a magneto-thermo-turbulent (MTT) model, which they implement into a cosmological disk galaxy simulation.Their recipe identifies regions of gas collapse through defining an MTT Jeans length  J, MTT dependent in part on the ratio of ther-mal to magnetic pressure, such that when the length of the gas cell surpasses  J, MTT a certain mass fraction of gas is converted into stellar particles.This fraction is determined by a local  ff , computed via the multi-scale model of Padoan & Nordlund (2011).However the method of computing star formation efficiency remains dependent on solely hydro-dynamical parameters.Additionally, this work along with others (e.g.Vazza et al. 2014;Rieder & Teyssier 2017) suffer from unrealistically slow magnetic field dynamos, due to the effects of limited spatial resolution on turbulent dissipation.There is still great need for a magnetised star formation recipe that can fully integrate the impact of magnetic fields on star formation efficiency, ideally with respect to new and improved models of generally unresolved small-scale turbulent dynamo.
In this work, we develop a new MHD-based star formation model that operates on sub-grid scales, deriving the local star formation efficiency  ff from the properties of the sub-grid magnetised ISM turbulence.We then use this star formation model in MHD simulations of isolated disk galaxies using the code RAMSES (Teyssier 2002).Our aim is to study the regulating effects of magnetic fields on star formation, identify how our magnetised model might impact the properties of galaxies on larger scales, and explain certain features observed in star-forming galaxies.
In Section 2, we derive our MHD sub-grid scale (SGS) star formation model from the pure hydrodynamics case previously developed by Kretschmer & Teyssier (2020).We present in Section 3 the numerical methods used in this paper, detailing additional SGS recipes for turbulence and stellar feedback.In Section 4, we describe the resulting galaxy disk morphology, magnetic field topology, and star formation history of our simulations.These properties and their implications are discussed in Section 5.

SUB-GRID SCALE STAR FORMATION MODEL
Our star formation model draws upon that of Kretschmer & Teyssier (2020), which requires the knowledge of the unresolved turbulent kinetic energy in each cell.This can be obtained using the so-called Large Eddy Simulations (LES) framework (see e.g.Schmidt 2014;Semenov et al. 2016, for relevant equations).In this framework a new hydro variable ( T ) is introduced to represent the sub-grid turbulent kinetic energy, together with a new equation describing its evolution with proper source and sink terms according to mixing length theory (see Kretschmer & Teyssier 2020, for details).We can then recover the one-dimensional velocity dispersion of turbulence from this new variable at the scale of the mixing length, chosen here equal to the size of the cell Δ: This allows us to prolong the properties of turbulence down to unresolved scales ℓ ≪ Δ using the classical scaling for supersonic Burgers turbulence with We consider the fast magnetosonic length ℓ f , defined as the scale at which the 1D velocity dispersion  1D (ℓ f ) is equal to the fast magnetosonic speed, where  A = / √︁ 4 is the Alfvén speed, and  s is the sound speed.This is in essence the MHD counterpart of the sonic length, the length scale at which turbulence becomes transonic.A critical assumption in this definition is that  A is constant at sub-grid scales.This is justified by simulations of supersonic and super-Alfvenic turbulence in molecular clouds, in which  A approaches a mean value that is independent of density and corresponds to a power-law relation  ∝  −0.4 (Padoan & Nordlund 1999).More recent simulations of supersonic magnetised turbulence performed by Padoan & Nordlund (2011) have shown that  A is nearly constant within molecular clouds for densities  ≳ 2 0 , where  0 is the mean density.Importantly, our assumption that  A is constant within our sub-grid model has no bearing on its properties at larger scales.It remains an open question as to whether  A is constant throughout the galaxy.
Introducing the parameter  =  2 s / 2 A , one can further relate the sound speed to the fast magnetosonic speed, as well as the sonic Mach number of the turbulence M s ≡  1D / s to the fast magnetosonic Mach number, For brevity, we introduce from here on the magnetic correction factor  ≡ (1 + )/, such that  f =  s √  and M f = M s / √ .Following for example Federrath & Klessen (2012), we assume that the unresolved gas density distribution is described by a log normal PDF.Note that the use of a log normal PDF is justified by the multiplicative process through which density is distributed: shocks in the isothermal supersonic turbulent ISM are amplified via collisions and interactions with other shocks (Vázquez-Semadeni 1994;Passot & Vázquez-Semadeni 1998;Kritsuk et al. 2007;Federrath et al. 2010).This corresponds to an additive process in log space, which by the central limit theorem produces a Gaussian distribution for logarithmic density.
Defining  = ln(/ ρ) as the logarithmic density, where  is the unresolved gas density and ρ is the mean gas density of the cell, this PDF takes the form where s as the mean logarithmic density and   is the standard deviation.Molina et al. (2012) derive a dependence between   and M f as where the parameter  depends on the nature (compressive or solenoidal) of the turbulence forcing in the medium.It is evident here that the stronger the magnetic field, the narrower the variance of the logarithmic density PDF  2  will be, leading to fewer higher gas density regions.This is further illustrated in simulations of MHD turbulence performed by Federrath & Klessen (2012), where they found that as magnetic field strength increased, density contrast decreased and fewer fragments were formed overall.We recognise that while our model assumes  2  is isotropic, this does not necessarily hold for clouds with a trans or sub-Alfvénic magnetic field (where M A ≡  1D / A ≲ 1).We hope that future work may capture such anisotropic regimes by incorporating density variance models that take into account preferential directions of density fluctuations (e.The left panel displays the shape of these contours in the non-magnetised sub-grid star formation recipe by Kretschmer & Teyssier (2020), while the right panel displays their shape as per our magnetised model in the case of  = 0.5 (  = 3).At the low M s limit, star formation efficiency depends on the entire sonic virial parameter of the cell, occuring at  vir,  ∼ 1 for the non-magnetised model and  vir,  ∼   vir,  ∼ 1 for the magnetised model.For M s ≳ 1 we see efficiency increasing in both models as  vir,  decreases, though the efficiency is generally dampened in the magnetised model as seen in the shift leftwards of efficiency contours.
We now assume that our resolution element is described by an ensemble of small uniform clouds of diameter equal to the fast magnetosonic length ℓ f and whose constant density follows the lognormal PDF introduced just above.Following Krumholz & McKee (2005), we identify star-forming regions as the population of these small clouds for which the virial parameter is less than one.Using the classical spherical transonic cloud model of Bertoldi & McKee (1992), for a cloud mass  = 4/3 3 and a cloud radius  = ℓ f /2, we have We can simplify this expression using the fast magnetosonic Mach number as We can also express this MHD-based virial parameter as a function of the old sonic length virial parameter   from Kretschmer & Teyssier (2020) as   =   /.Given that these clouds will collapse and form stars only if the virial parameter is smaller than one, this defines a critical collapse density at   = 1: Δ 2  This critical density can be expressed in terms of the mean density ρ and the virial parameter of the cell  vir,  as where Therefore the critical logarithmic density can be calculated as Note that this model breaks down at low Mach number (Federrath & Klessen 2012).Following Kretschmer & Teyssier (2020), we modify Equation 14 such that for M f ≤ 1, the lognormal PDF becomes a delta function and the stability criterion must be applied to the cell as a whole.Just as in Kretschmer & Teyssier (2020), this requires re-defining the critical logarithmic density as With the critical logarithmic density now properly defined in the general case, we can compute the local star formation efficiency using the multi-free fall time approach of Federrath & Klessen (2012), Figure 1 displays the differences in star formation efficiency between the original star formation model of Kretschmer & Teyssier (2020) and our MHD version for a typical magnetised case with  = 0.5.While the efficiency contours have similar shapes, the magnetic field effectively reduces the star formation efficiency in the parameter space defined by the cell sonic virial parameter and the sub-grid turbulence sonic Mach number.Note that where  vir, is the virial parameter of the whole cell in the nonmagnetised model.The sonic virial parameter emerges in the nonmagnetised limit,  → ∞ or  → 1.The presence of magnetic fields in a potentially collapsing molecular core will increase the value of the virial parameter.Hence smaller effective sonic virial parameters are necessary for the magnetised model to produce the same efficiencies of star formation.In Figure 1, this corresponds to a critical value of  vir, ∼ 0.1 in the magnetised case with  = 0.5, versus  vir, ∼ 1 in the non-magnetised case.
We see in both models the dramatic break in efficiency at low sonic Mach numbers, as in this regime both models allow for star formation via a single threshold virial parameter.This threshold differs by a factor of 1/, which corresponds to the relation  vir, =  vir,  / in the low M s limit.At larger sonic Mach numbers, the virial parameters  vir, and  vir,  approach equality but the critical logarithmic density increases in the magnetised model relative to the non-magnetised model: As such, efficiency is suppressed when magnetic fields are present, as lower effective sonic virial parameters are required in the magnetised model to reach the same efficiencies as before.This is illustrated by the slight translation leftwards in the efficiency contours between the two plots.

NUMERICAL METHODS
To compare our new star formation model to previous models, we have performed a suite of isolated disk galaxy simulations with the octree-based Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002).In these simulations, dark matter and stars are modelled as a collisionless fluid, while the gas component follows the ideal MHD equations (see below).These different fluids are coupled through self-gravity.
Our initial conditions are taken from the high-resolution AGORA simulation project (see Kim et al. 2016).They correspond to a typical Milky-Way-sized, gas-rich disk galaxy embedded in its parent dark halo.The simulation is initialised with 5 different components: the dark matter halo, the gaseous halo, the stellar disk, the stellar bulge, and the gaseous disk.We use the DICE code to generate our initial particle and gas cell distributions using a Metropolis Hasting Monte Carlo Markov Chain algorithm (see details in Perret 2016).
The dark matter halo is distributed according to a Navarro et al.The gaseous halo is filled with a very low uniform density of  H = 10 −6 cm −3 and a uniform temperature of 10 6 K, zero initial velocity, and zero metallicity.
The disk component contains a total mass of  d = 4.297 × 10 10  ⊙ with scale length  d = 3.432 kpc, scale height  d = 0.1 d , gas mass fraction  gas =  d,gas / d = 0.2, which corresponds to a disk gas-to-star ratio of 1/4.
Both the stellar and gas component of the disk follow an exponential density profile where index  stands here either for gas or stars,  is the cylindrical radius, and  is the vertical height from the disk plane.The central disk densities are computed using Both gas and stellar disks are truncated at 3 d and 3 d .The stellar bulge is modelled with a Hernquist (1990) density profile and stellar mass  b, * = 4.297 × 10 9  ⊙ .Our initial gas disk is turbulent.We add random velocity fluctuations on top of the overall rotating flow.These fluctuations follow a Burgers  −2 kinetic energy spectrum normalized to  1D = 20 km/s.We found this to be a crucial addition compared to previous work (see Kim et al. 2016).Indeed, a disk in perfect hydrostatic equilibrium would collapse vertically as it cools.The resulting razor-thin gas slab will trigger an unrealistically strong starburst.
Our new star formation model is quite sensitive to the magnitude of the sub-grid turbulence.We initialise the sub-grid turbulence to the stationary value obtained by balancing source terms and sink terms in the LES sub-grid kinetic energy equation (Kretschmer & Teyssier 2020).These precautions explain why our star formation history is so smooth and slowly declines in time without a spurious prominent starburst at startup (see the Results section).
We are particularly interested in the evolution of the magnetic field and its impact on star formation.The outcome of the simulation will depend on the initial field strength and topology of the field.In this paper, we restrict ourselves to a very weak initial field, so that the final field strength and topology are solely a consequence of our galactic dynamo.We will study the impact of the initial (or fossil) magnetic field in a follow-up paper.The smallest value one can realistically consider for the initial field is the one produced by the Biermann battery during the linear regime of cosmological density perturbations (Naoz & Narayan 2013).We initialise our magnetic field with a constant magnitude of 10 −20 G and parallel to the  direction.

Details on the adaptive grid
Simulations begin from a uniform grid covering our entire computational box of size (320 kpc) 3 , with a coarse grid of resolution 128 3 corresponding to our coarsest refinement level ℓ min = 7 and a coarsest cell size of 2.5 kpc.Our finest refinement level was set to ℓ max = 14 so that our finest resolution element size is ∼ 20 pc.We employ a quasi-Lagrangian strategy: we trigger a new cell refinement if the baryonic mass (gas plus stars) within the cell exceeds 10 3  ⊙ , or if the cell contains more than 8 collisionless particles.Boundary conditions are set to isolated for the gravity solver and out-flowing for the MHD solver.

Details on the MHD and gravity solvers
We give her more details on our MHD solver, as it entails the core of our study.We solve the following ideal MHD equations, assuming that the gas is ionised enough to justify neglecting all non-ideal MHD effects.Written in conservative form, these equations are where  is the mass density, v is the momentum, B is the magnetic field, g is acceleration due to gravity,  = 1 2 v 2 +  + 1 2 B 2 is the total fluid energy,  tot =  gas + 1 2 B 2 is the total pressure.In the energy equation, Γ is the heating function and Λ is the cooling function,  gas is the pressure of the gas, and  is the specific internal energy.We consider an ideal gas equation of state with adiabatic exponent  = 5/3.
Our numerical scheme is a standard second-order Godunov scheme with the MUSCL-Hancock scheme coupled to the Constrained Transport method for the induction equation (Teyssier et al. 2006).The CT scheme enforces the solenoidal constraint ∇ • B = 0 to machine precision accuracy.As a result, our numerical scheme conserves mass, linear momentum and total fluid energy, as well as magnetic flux.
As explained in Fromang et al. (2006), we use the HLLD Riemann solver to compute intercell fluxes of mass, momentum and energy.We use a 2D version of the HLLD Riemann solver to compute the electric field at the edge of the faces where the magnetic field is defined (Fromang et al. 2006).Second-order reconstruction of primitive variables is performed using the MinMod slope limiter (see e.g.Roe 1986).When new cells are refined, we use straight injection of conservative variables from the parent cell.Last but not least, we use for self-gravity the Multigrid Poisson solver available in RAMSES (Guillet & Teyssier 2011).

Details on the sub-grid physics
In this paper, we compare two simulations with two different star formation models.The first employs, as a reference, the non-magnetised sub-grid star formation recipe developed originally in Kretschmer & Teyssier (2020).The second uses our new sub-grid star formation model, in which the effects of magnetic fields on the local efficiency are considered (see previous section).Both simulations incorporate the LES model for turbulent kinetic energy, providing us with the sub-grid turbulent kinetic energy that enters our sub-grid star formation models.
We also use the same turbulent kinetic energy information for our implementation of a sub-grid magnetic dynamo, as described briefly below.Beyond these important sub-grid models, we have used the traditional galaxy formation ingredients provided by RAMSES code, such as equilibrium H and He gas cooling (e.g.Katz et al. 1996) plus metal cooling at both low and high temperature (Sutherland & Dopita 1993), heating by a standard Haardt & Madau (1996) UV background and a self-shielding density of  H = 0.01 cm −3 (Aubert & Teyssier 2010).

Sub-grid Turbulent Dynamo
As in Liu et al. (2022), we model the impact of unresolved turbulence on the magnetic field evolution using a mean-field approach (see e.g.Schmidt & Federrath 2011, and references therein).Assuming homogeneity and isotropy of the turbulence, one can derive a tensor  relating the mean magnetic field B to the turbulent electromotive force caused solely by unresolved turbulent fluctuations in the velocity and magnetic fields, Ē: A simple model for  has been proposed by Liu et al. (2022) that depends on the local sub-grid turbulent velocity dispersion where  mag = B2 /8 is the mean-field magnetic energy density,  is a turbulent dynamo quenching parameter, and  T = 3 2  2 1D is the sub-grid turbulent kinetic energy density.
For these simulations, we adopt the fiducial quenching parameter  = 10 −3 (Liu et al. 2022).This means that our sub-grid turbulent dynamo will vanish entirely as soon as the magnetic energy reaches 0.1% of the local turbulent energy.For typical gas disk conditions with  H ≃ 1 cm −3 and  1D ≃ 1 − 10 km/s, this corresponds to a mean field of roughly 0.02 − 0.2 G.
We also limit the sub-grid turbulent dynamo model to only cooling and star-forming gas, by setting  = 0 for regions with number densities less than the self-shielding threshold   = 10 −2 cm −3 .In this manner, strong magnetic fields in the circumgalactic medium can only arise from galactic outflows and not from a local turbulent dynamo, unless it is resolved by the code, which is unlikely given our limited resolution.Previous works have shown that when the resolution is high enough, dynamo amplification could emerge from resolved motions without relying on a sub-grid turbulent dynamo model (see Rieder & Teyssier 2016;Vazza et al. 2018;Steinwandel et al. 2022).

Sub-grid Stellar feedback
Our stellar feedback model follows that of Kretschmer & Teyssier (2020), which quantifies for each stellar particle (with typical mass 10 5  ⊙ ) the scalar momentum due to sub-grid supernovae explosions, and injects this momentum into the surrounding gas.Assuming a uniform distribution of individual supernova (SN) explosion times for a given stellar particle between  start = 3 Myr and  end = 20 Myr, the rate of supernovae can be computed as where  SN = 0.2 is the mass fraction of the stellar particle producing supernovae,  ini, * is the initial stellar particle mass, and  SN = 10  ⊙ is the typical supernova progenitor mass.Hence the expected number of supernovae within a timestep  is For each timestep, we can draw the actual number  SN of supernovae for a stellar particle of mass  , * from a Poisson distribution with average ⟨ SN ⟩.The minimum scale at which momentum must be injected is determined by the cooling radius  cool , which is modelled as where  is the gas metallicity and  H the number density of hydrogen.To account for the cases in which  cool is either resolved or unresolved by at least one grid cell with length Δ min , we calculate the total resulting momentum as where This momentum is then injected isotropically into the surrounding gas cells using the numerical techniques described in Kretschmer & Teyssier (2020).

SIMULATION RESULTS
In this section, we describe in detail our simulation results.Our final snapshot has reached a time of 1.4 Gyr, for which the galaxy is forming stars steadily at a rate of 1 solar mass per year.A galactic fountain is clearly visible and plays an important role in regulating the gas available for forming stars in the disk.After an initial phase of fast amplification due to our sub-grid mean field dynamo model, the magnetic field quickly saturates slightly below equipartition with the gas turbulent energy and develops a strong toroidal structure, although in our particular case, it appears to be still dominated by small scale fluctuations.

Disk morphology
Figure 2 shows edge-on and face-on images of the gas surface density and the mass-weighted temperature of our simulated galaxy at  = 1.4 Gyr.The gas surface density reaches its maximum value around 10 3  ⊙ pc −2 in the nuclear region, with multiple cold and dense star-forming gas clumps visible along the spiral arms.We note that the apparent spiral structure of the gas clumps themselves, while differing from observed complex morphologies of such objects, is in fact expected.The limited resolution of our simulations (∼ 20 pc scale) prevents us from resolving the clouds' internal structure, hence these clouds cannot fragment into the necessary substructures.However, the Coriolis force associated with galactic rotation still moulds clouds into a spiral shape.We expect that, were our simulations able to resolve scales below typical gas cloud radii, the gas clouds of Fig. 2 would fragment further into various dense clumps.Such fragmentation is visible in sub-parsec resolution simulations of spiral galaxies, with smaller gas clumps roughly distributed along spiral arms like the "beads of a string" (see e.g.Renaud et al. 2013).
The multiphase nature of the interstellar gas (ISM) appears evident in Figure 2 with low density, hot supernova-driven bubbles ( ≃ 10 7 K) coexisting alongside intermediate density warm gas ( ≃ 10 3 − 10 4 K) and high density cold clumps ( ≃ 10 − 100 K).In the edge-on view, we can clearly see the galactic fountain as thin streams of gas are ejected outwards from and fall back towards the galactic plane.These processes heat the low density circumgalactic medium (CGM), which as seen in Figure 2 ranges in temperature from 10 4 − 10 7 K. Due to our simulation's limited resolution outside of the galactic disk, the multiphase structure of the CGM -which has been observed to comprise of not only hot virialised gas, but also atomic hydrogen (with temperatures  ∼ 10 4 K) and molecular hydrogen ( ∼ 50 K) in circumgalactic gas (for review, see Tumlinson et al. 2017) -cannot be fully reproduced.Given that our work is primarily focused on star formation within the disk, we leave more robust treatments of and refinement within the CGM for future computational experiments.

Magnetic field topology
After 1.4 Gyr, our simulated galaxy displays a relatively ordered magnetic field that closely follows the gas density (see Fig. 3).The central star-forming region has the highest magnetic field strengths around 100 G, while in the disk, dense clumps with low temperatures  ≲ 10 2 K and strong magnetic fields  ≳ 10 G accumulate along the spiral arms.The magnetic field strength generally decreases with increasing radius.Our magnetic energy is roughly in equipartition with the thermal energy but not with the turbulent energy (see Fig. 4).Moreover, when averaged over cylindrical bins, the magnetic energy is dominated by the fluctuating radial and tangential components, while the mean radial and tangential fields are almost zero.This is consistent with our very weak initial magnetic field and indi-10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8  6.Left panel: mass-weighted (blue) and volume-weighted (red) gas fractions as a function of temperature, radially binned in 100 bins across a disk profile around the centre with maximum radius of 10 kpc and height of 3.2 kpc.Different ISM phase regions are annotated, along with the total mass-weighted (blue) and volume-weighted (red) mass fractions of each phase.They are the cold neutral medium (CNM;  < 250 K), the thermally unstable (UNM; 250 K ≤  < 3 × 10 3 K) and warm neutral medium (WNM; 3 × 10 3 K ≤  < 5 × 10 4 K), the warm ionised medium (WIM; 5 × 10 4 K ≤  < 4 × 10 5 K), and the hot ionized medium (HIM; 4 × 10 5 K ≤ ).Right panel: Mass-weighted gas fraction as a function of number density over the same binned cylindrical region as the left panel.We show this density PDF for all of the gas in the region (black), the cold neutral medium phase (blue), the thermally unstable phase (green), the warm neutral medium phase (light blue), the warm ionised phase (orange), and the hot ionised phase (red).cates that our magnetic field evolution converges towards a saturated small-scale turbulent dynamo.
In Figure 5, we have plotted the mass-weighted, radial, tangential and vertical components of the magnetic field, seen both face-on and edge-on.Although the toroidal component appears to be the strongest, we see also a strong radial component and a weak vertical component.We obtain main field reversals as revealed by the multiple changes of sign of both the radial and the tangential component, explaining why the average values vanish.Along the spiral arms,   ∼   ∼ 5 G, while in the regions between spiral arms, the magnetic field component magnitudes slightly decrease to ∼ 2 G.

Structure of the ISM
We distinguish the various phases of our simulated galaxy ISM through the mass-weighted and volume-weighted density distributions of the gas, visualised as a function of temperature in the left panel of Figure 6.Cold, dense gas (which we demarcate via a temperature criterion of  < 250 K) clearly dominates at about 73% of the mass fraction, consistent with our expectations for a star forming disk.However, there is not a clear distinction between this cold neutral medium (CNM) and a warm neutral medium (WNM) in our ISM, aside from a local peak in the gas fraction that occurs around  ∼ 2 × 10 4 K.In this final simulation snapshot the galaxy is still relatively gas-rich (with a total gas fraction  gas ∼ 20%) and actively star-forming, which suggests the presence of a thermally unstable neutral medium (UNM).Furthermore, our limited numerical reso-  lution prevents proper spatial resolution of this transitional medium between the CNM and WNM.Hence we consider in the left panel of this figure the UNM combined with the WNM, as opposed to their fractions individually, to account for gas within the temperature range of 250 K ≤  < 5 × 10 4 K. Gas at these temperatures make a non-negligible contribution to the mass fraction of ∼ 26%.Additionally, we label an apparent warm ionised medium (WIM) with temperatures 5 × 10 4 K ≤  < 4 × 10 5 K, and a hot ionised medium (HIM) with temperatures 4 × 10 5 K ≤ .These phases make up the least of the gas mass, while being the most volume filling: the HIM fills 51.6% of the volume, and the WIM fills 22.3%.The UNM/WNM phase also makes a non-neglible contribution to the volume with a 24% volume-filling fraction, and the smallest (1.7%) volume-filling fraction can be attributed to the CNM.These values align with our expectations and previous results on the ISM structure of star forming, gas-rich disks (e.g.Hopkins et al. 2012).
Differences in our denoted ISM phases are further illustrated in the right panel of Fig. 6, where we plot mass-weighted gas fraction as a function of number density.Here we have differentiated the UNM from the WNM, using gas temperature cutoffs of 250 K ≥  < 3 × 10 3 K and 3 × 10 3 K ≥  < 5 × 10 4 K for the UNM and WNM respectively.The HIM and WIM correspond to lower density gas, with their probability density functions (PDFs) centred at number densities of 10 −4 cm −3 and 10 −3 cm −3 respectively.The UNM and WNM lie largely in the expected intermediate density range of 0.01 cm −3 -1 cm −3 , while the CNM is concentrated at high densities (> 1 cm −3 ).
Figure 7 shows 2D mass-weighted histograms of temperature and sub-grid velocity dispersion as a function of gas density.Similar to Fig. 6, the temperature histogram reflects how the bulk of gas mass lies in a cold phase with number densities ∼ 100cm −3 and temperatures ≲ 50 K.Warmer, intermediately-dense gas ( ∼ 1000 K,  H ∼ 0.1 cm −3 makes up the second largest portion of the mass, while the most diffuse and hot gas ( > 10 5 K,  H < 0.01 cm −3 makes up the least.We note that a seemingly concentrated portion of mass at temperatures  ∼ 10 6 K and densities  H ∼ 10 −6 cm −3 are again numerical artefacts, attributable to the lack of resolution in our galaxy CGM.The sub-grid velocity dispersion reaches its minimum around  1D ∼ 1 km s −1 for a moderate density around 1cm −3 , close to the average gas density in the disk.Inside dense clumps, however, the sub-grid velocity dispersion increases following the expected adiabatic relation  1D ∝  1/3 (Semenov et al. 2016).

Why is the Alfvén speed so constant?
The 2D mass-weighted histograms of magnetic field strength and Alfvén speed as a function of gas density, shown in Figure 8, reflect typical magnetic field strengths and densities in the disk on the order of 10 G and 1cm −3 , respectively.Maximum strength of 100 G is particularly evident in the densest regions, namely cold clumps along the spiral arms and the nuclear region of the disk.A striking result visible in this figure is that the Alfvén speed  A appears as mostly constant, with  A ∼ 10 km s −1 for a large range of number densities from 0.1 to 100cm −3 .This reflects the same scaling  ∝  1/2 , predicted by theory and simulations (Padoan & Nordlund 1999) on small scales, deep inside molecular clouds.Note that there is however significant scatter that covers the classical adiabatic compression relation  ∝  2/3 (Crutcher et al. 2010).
Why is the Alfvén speed roughly constant across scales?We propose here two mechanisms working in tandem.First, on very large scales, the magnetic energy saturates at a fixed fraction of the turbulent kinetic energy.This implies  2 ∝  2 1D , and is the classical argument of equipartition proposed by Chandrasekhar & Fermi (1953); Fiedler & Mouschovias (1993); Cho & Vishniac (2000); Groves et al. (2003).Since the gas velocity dispersion is roughly constant throughout the disk (see Fig. 4), this leads to both the expected scaling  ∝  1/2 and the correct normalisation when using the average disk density  0 ≃ 1 cm −3 and the average magnetic field  0 ≃ 10 G.
What is less obvious is why this scaling relation persists inside dense clumps.The second, less classical explanation we propose here is the 2D magnetised collapse of these dense clumps from a razor thin (Toomre-unstable) disk.Assuming a tangential (or radial) magnetic field  embedded in a collapsing thin cloud of constant height H and shrinking radius R, conservation of mass writes  =  2 and conservation of magnetic flux writes  = 2.Combining these two relations while removing the radius, we obtain assuming that the typical Toomre-unstable clump has initial radius  ≃ , initial density  =  0 and initial magnetic field  =  0 .This two-step model explains why the Alfvén speed remains constant from the largest scales in the disk down to collapsing molecular clouds.Furthermore it remains constant at even smaller scales, as advocated by the work of Padoan & Nordlund (1999), although these scales are for us unresolved and captured only by our sub-grid model.

Effect of our new star formation model
Our sub-grid dynamo methodology allows the emergence of a strong magnetic field characterised by a roughly constant Alfvén speed of magnitude  A ≃ 10 km/s.This translates into a plasma  ≃ 10 −3 inside dense star-forming molecular clouds in which the sound speed is roughly  s ≃ 0.2 − 0.4 km/s.The nature of the turbulence inside the clouds is clearly dominated by the Alfvénic Mach number  A =  1D / A ≃ 1 − 2 as opposed to the sonic Mach number  s =  1D / s ≃ 10 − 20, therefore we expect the magnetic field to be the dominant component in regulating star formation at small scales.
In particular, because the large-scale turbulent velocity dispersion and the Alfvén speed are both roughly constant, as well as the sound speed inside molecular clouds, we expect an overall constant offset of the star formation efficiency.
The first difference between our magnetised star formation model and the non-magnetised model of Kretschmer & Teyssier (2020) can be seen in the star formation history of our isolated galaxy shown in Figure 9.In both cases, we observe an overall slow exponential decline with a prompt initial burst.Note that this initial burst is considerably smaller than in other simulations of isolated galaxies (see e.g.Stinson et al. 2006) because we have included decaying turbulence in our initial conditions, preventing the formation of an unrealistic razor thin disk after the initial vertical collapse.
Including the magnetic field reduces the initial SF rate from 5 − 6  ⊙ yr −1 to 2 − 3  ⊙ yr −1 .As a result, the gas depletion time scale has been reduced by a factor of 2, resulting in a flatter SF history.At the end of the simulation, both models predict a SF rate around 1  ⊙ yr −1 , but the magnetised case contains almost twice the amount of gas of the non-magnetised model by the end of the simulation.
In order to further understand the effect of the magnetic field on our star formation recipe, we have plotted in Figure 10 the local star formation efficiency in each cell for both the non-magnetised case (left panel) and the magnetised case (right panel) as a function of the local gas density.Most notably our magnetised model results in the distribution of star-forming cells clearly shifting to higher densities, from ∼ 10 cm −3 to ≳ 100 cm −3 .Overall however, the efficiency ranges themselves are quite similar in both cases, namely because the typical gas density is a factor of ∼ 2 larger in the magnetised case.This particular snapshot at 1.4 Gyr appears to have quite low efficiencies per free-fall time, with the non-magnetised and magnetised models reaching peak efficiencies of 0.7% and 0.35% respectively.However, we emphasise this is only a snapshot and hence not representative of the larger star forming history; preceding snapshots demonstrate more active star formation, with the peak efficiency lying at 1% for the densest cells in both model cases.
The total gas Schmidt law is represented in the left panel of Figure 11, where we plot the gas surface density Σ gas against the SFR surface density Σ SFR .Σ SFR is computed within a disk of radius 10 kpc and height 3 kpc, by azimuthally summing the mass of young stars (formed over the last 100 Myr) in radial bins of size 100 pc.There is a slight difference in the power law index, where  ∼ 2.5 and  ∼ 2 for the non-magnetised and magnetised star formation models respectively.The power law in the magnetised model however shows better agreement with observed values for individual spiral galaxies, which range from  ∼ 1.2 − 2.1 (Wong & Blitz 2002).
The most striking difference in the left panel of Figure 11 is the much higher gas surface density observed in the nuclear region.We reach 500  ⊙ pc −2 using the magnetised star formation model against only 100  ⊙ pc −2 .The particularly strong magnetic field in the nuclear region reduces the local star formation efficiency.This in turn requires more and denser gas to accumulate, in order to produce enough supernovae feedback to regulate the global star formation efficiency.This self-regulation mediated by feedback explains why our overall SF history is not strongly affected by our new model, while locally it is strongly modified.Interestingly, the magnetised star formation model allows also smaller values of the SFR at low surface densities.We see in the left panel of Figure 11 a tail of low SFR for surface densities below 10  ⊙ pc −2 .This tail is reduced for the non-magnetised star formation model, suggesting that the magnetic field in the outer regions is high enough with  ≃ 3 − 4 G to also quench star formation there.

DISCUSSION AND COMPARISON TO OBSERVATIONS
To assess the consistency of our simulation results with observations, we compare the star formation rate surface density Σ SFR , the total gas surface density Σ gas , and the volume-weighted average magnetic field strength  to observed relations found in literature.The left panel of Figure 11 displays the average total gas Schmidt law as fitted by Bigiel et al. (2008), for a sample of 18 nearby galaxies at sub-kpc resolution.Similar to our analysis, they provide values for Σ SFR and Σ gas across azimuthally-averaged radial profiles.While our non-magnetised and magnetised star formation models appear to under-predict their averaged power law, they both lie well within the observed spread around this power law.Bigiel et al. (2008) emphasise that Σ SFR and Σ gas differ largely across the galaxies in their sample, with the average Σ SFR possessing an RMS scatter of ∼ 0.3 dex.In fact, where some galaxies are well described by a single power law, others are not.Hence we wish to make the key point here that our simulations necessarily align with the generally observed range of values in SFR-gas space.
Furthermore, we note that the total gas Schmidt relations of both simulation models appear to possess a "knee" around the value of Σ gas ∼ 10  ⊙ pc −2 .Both Bigiel et al. (2008) and previous studies (Morris & Lo 1978;Martin & Kennicutt 2001;Wong & Blitz 2002) have demonstrated how the HI radial profiles in spirals reflect a remarkably consistent saturation at high column densities, with a threshold of Σ gas ≈ 9  ⊙ pc −2 .H 2 surface density, on the other hand, generally reflect no similar limitations.Therefore the apparent difference between the Σ gas − Σ SFR relation at lower and higher column densities, as opposed to a single power law, arguably reflects a transition between an HI-dominated ISM to an H 2 -dominated ISM that is consistent with observations.
We additionally compare the relation between magnetic field strength and Σ SFR produced by both models, to that inferred from radio continuum observations and reported in the recent work of Manna  Figure 11.Left: Gas surface density Σ gas and star formation rate surface density Σ SFR relation, within simulated galaxies using either the non-magnetised (blue) or magnetised (orange) sub-grid star formation models.The high gas surface density tail, which we attribute to the nuclear regions, is emphasised by thicker markers.Overplotted in the red dotted line is the average Σ gas − Σ SFR law empirically found by Bigiel et al. (2008).Right: Star formation rate surface density Σ SFR and magnetic field strength  relation, within simulated galaxies using either the non-magnetised (blue) or magnetised (orange) sub-grid star formation models.Similarly, the high SFR surface density tail is additionally emphasised by thicker markers.Overplotted in the red dotted line is the Σ gas −  law empirically found by Manna & Roy (2023).In both subplots, we consider mass-weighted quantities within a disk around the centre of radius 10 kpc and height 3.2 kpc, split into 100 cylindrical bins.Σ SFR is computed assuming a stellar age of 100 Myr.
& Roy (2023).In the right panel of Figure 11, we show their average relation found using a sample of 11 nearby galaxies is in striking agreement with our simulations.Interestingly, our non-magnetised star formation model appears to result in central magnetic fields that are too weak.Previous semi-analytical models of the  − Σ SFR relation, such as that proposed by Schleicher & Beck (2013), similarly predict central magnetic field values that are too weak by a factor of 2 to 3. As Manna & Roy (2023) discuss, a key parameter of their model is the fraction of turbulent kinetic energy that is converted to magnetic energy, otherwise known as the saturation level of the turbulent dynamo (  sat ).While Schleicher & Beck (2013) assume an  sat ∼ 5%, in our case (see Figure 4) we find a significantly larger value of  sat ≃ 25%.This could partly explain why our simulations are in better agreement with observations.

CONCLUSIONS
In this paper, we have performed ideal MHD simulations of an isolated, Milky Way-like galaxy using a new sub-grid model for our star formation recipe that accounts for the presence of magnetic fields.
For this, we used a Large Eddy Simulation approach to compute the velocity dispersion of the sub-grid turbulence, with the same LES model used to power a sub-grid turbulent dynamo.This mean-field  dynamo allows us to amplify our vanishingly small initial magnetic field to a value close to (but smaller than) equipartition.Our star formation model is based on the Sub-grid Scale (SGS) multi-freefall model proposed originally by Federrath & Klessen (2012), using both the sub-grid turbulence and the magnetic field computed selfconsistently by the RAMSES code, and extends numerical methods described in Kretschmer & Teyssier (2020).As discussed in Section 2, the effect of the magnetic field is to both reduce the width of the log-normal gas density PDF and make the collapse criterion for molecular cores more difficult to achieve.Overall, as shown in this paper, the efficiency of star formation at small, sub-grid scales, can be significantly reduced if the magnetic field strength is high enough.Interestingly, we find that the Alfvén speed is remarkably constant from galactic scales down to molecular cloud scales.It value is of the order of 10 km/s, corresponding exactly to the average sound speed of the warm gas, but smaller that the average velocity dispersion, which in our case is closer to 20 km/s.The origin of this remarkable uniformity of the Alfvén speed is due to 2 combined effects: 1-the saturation of the turbulent dynamo at a constant ratio of magnetic to turbulent energy, 2-the 2D collapse of razor thin molecular cloud with a magnetic field aligned with the disk.
Including the magnetic field in the multi-free fall star formation model reduces the overall initial star formation rate by nearly a factor of 2 when compared to the non-magnetised reference case.However, star formation somehow self-regulates so that the final SFR is around 1  ⊙ yr −1 after ≳ 1 Gyr in both cases.The most striking effect of reducing the local SF efficiency at small scales is to compress the gas to a much higher density (and field strength) by a factor of up to 5 in the nuclear region.As a result, the local SF efficiency does not change much compared to the non-magnetised case, and lies close to the value adopted in most galaxy formation simulations,  ff ∼ 0.01 (see Krumholz & Tan 2007;McKee & Ostriker 2007).The dispersion of the local efficiency is also in good agreement with observations of nearby galaxies, where  ff = 0.003 − 0.026 with a median value around 0.007 (Utomo et al. 2018).
Our results compare well with the mean SFR surface densitymagnetic field relation observed in nearby galaxies (Manna & Roy 2023).Although encouraging, our results are still lacking many realistic aspects of galaxy formation and cosmic magnetism.Because our initial magnetic field was vanishingly small, only the turbulent dynamo is at work here.As a result, our mean toroidal field is very small and the magnetic field geometry is dominated by the fluctuating component.We plan to explore in the future different initial field configurations, including ones with a strong pre-existing toroidal field.This could have an additional impact on the overall star formation efficiency.
Another important caveat of this model is that we consider the effect of the magnetic field only through the magnetic pressure, completely neglecting anisotropic effects.However as shown in the recent work of Beattie et al. (2021), the gas density distribution in sub-Alfvènic turbulence is highly anisotropic.Their theoretical model could be used to refine our model and account for this anisotropy.
Furthermore, as this study focuses on isolated galaxy simulations, we cannot account for the effect of the cosmological environment (inflows) on the galactic magnetic field topology and star formation.Past studies agree that on average, isolated galaxies have higher contents of HI and lower star formation rates than galaxies residing within clusters (see Boselli & Gavazzi 2006, and references therein).Alternatively, motions in a forming galaxy cluster combined with turbulence in the intracluster gas is believed to amplify the magnetic field in cluster galaxies still undergoing major mergers (Roettiger et al. 1999a,b;Subramanian et al. 2006).In future work, we intend to examine our star formation model with these factors in mind, starting with more comprehensive simulations of isolated galaxies that can be followed up by those of larger-scale galaxy clusters.

Figure 1 .
Figure1.Contours of the efficiency per free-fall time  ff , computed for various combinations of sonic virial parameter  vir,  and sonic Mach number M s .The left panel displays the shape of these contours in the non-magnetised sub-grid star formation recipe byKretschmer & Teyssier (2020), while the right panel displays their shape as per our magnetised model in the case of  = 0.5 (  = 3).At the low M s limit, star formation efficiency depends on the entire sonic virial parameter of the cell, occuring at  vir,  ∼ 1 for the non-magnetised model and  vir,  ∼   vir,  ∼ 1 for the magnetised model.For M s ≳ 1 we see efficiency increasing in both models as  vir,  decreases, though the efficiency is generally dampened in the magnetised model as seen in the shift leftwards of efficiency contours.

Figure 2 .
Figure2.Face-on and edge-on maps of gas surface density Σ gas (top row) and temperature (bottom row) of our simulated galaxy.Both quantities are densityweighted, and the gas surface density is averaged along the line of sight within a box 320 kpc wide in each direction.For the temperature map, face-on and edge-on views are taken from the  = 0 and  = 0 plane respectively.

Figure 3 .
Figure 3. Face-on and edge-on views of the density-weighted magnetic field strength, within the  = 0 plane (left) and  = 0 plane (right).

Figure 4 .
Figure 4. Total resolved plus sub-grid velocity dispersion  tot (red line), resolved velocity dispersion  res (green line), sound speed  s (orange line), and Alfvén velocity  A (purple line) as varying with radius  in the central region of the disk.The disk profile has a maximum radius of 10 kpc and height of 3.2 kpc.Quantities are computed as mass-weighted averages within 100 cylindrical bins.

Figure 7 .
Figure 7. Mass-weighted histograms of temperature (left) and sub-grid velocity dispersion   (right) versus number density  H .

Figure 8 .
Figure 8. Mass-weighted histograms of magnetic field strength and Alfvén velocity  A versus number density  H .

Figure 9 .
Figure9.A comparison of the star formation history in both simulated galaxies, using either the non-magnetised sub-grid star formation recipe from Kretschmer & Teyssier (2020) (blue) or our magnetised sub-grid recipe (orange).

Figure 10 .
Figure10.Mass-weighted 2D histograms of efficiency per free-fall time  ff versus number density  H , taken from a snapshot at 1.4 Gyr.The left panel is from our control simulation with the non-magnetised sub-grid star formation recipe, while the right panel is from our simulation that incorporates the magnetised model.