Time evolution of the galactic B − ρ relation: the impact of the magnetic field morphology

,


Introduction
Magnetic fields are fundamental for several crucial processes in the interstellar medium (ISM).This influence naturally translates into measurable effects in galaxy evolution (e.g., Pakmor & Springel 2013;van de Voort et al. 2021;Martin-Alvarez et al. 2020;Whittingham et al. 2021).On galactic scales, magnetic fields have the ability to guide the propagation of cosmic rays (Fermi 1949;Kulsrud & Pearce 1969;Cesarsky 1980;Desiati & Zweibel 2014;Shukurov et al. 2017) and determine the locations where molecular clouds form, through the interplay between gravity and magnetic buoyancy (Parker 1966;Mouschovias et al. 1974).Very importantly, magnetic fields can significantly influence star formation (see, e.g.Hennebelle & Inutsuka 2019;Pattle et al. 2023, for recent reviews).When the mass-to-flux ratio (M g /Φ B ) exceeds a critical threshold, the magnetic support against gravity can slow down or even stop the collapse of a cloud (e.g., Mestel & Spitzer 1956;Mouschovias & Spitzer 1976).
To quantify the importance of magnetic fields with respect to other forces several metrics are used, such as the plasma beta (β) to assess the significance of magnetic fields relative to the ther-⋆ akonstantinou@physics.uoc.gr⋆⋆ evangelia.ntormousi@sns.itmal pressure.Plasma β represents the ratio of thermal pressure, P th , to magnetic pressure, P mag , and if it is greater than 1 the gas is thermally dominated (e.g., Pattle et al. 2023).It is essential to note that accurately measuring the magnetic field strength and the thermal pressure in the ISM can be challenging.
Similarly, when considering the importance of magnetic fields in relation to gravity, measuring the M g /Φ B in astrophysical systems poses challenges because the complexity of astrophysical environments presents significant obstacles in accurately determining the magnetic flux (e.g., Pattle et al. 2023).
To overcome the challenges in directly measuring the M g /Φ B and to quantify the importance of magnetic fields, astronomers often resort to the B − ρ relation, which is typically described by a power law.This relation provides insights into how the magnetic field reacts to condensations and offers information about the collapse geometry of molecular clouds (Mestel 1965;Mouschovias & Spitzer 1976;Crutcher et al. 2010;Tritsis et al. 2015;Pattle et al. 2023).Mestel (1965) was the first to derive a theoretical B − ρ relation for a spherically symmetric, isotropically collapsing cloud, based on the mass and magnetic flux conservation: where M g is the gas mass, R the cloud radius, Φ B the magnetic flux and B the magnetic field, assumed uniform.With R ∝ ρ −1/3 from Eq. 1, Eq. 2 gives B ∝ ρ 2/3 .Some years later, Mouschovias (1976a,b) derived a B − ρ relation for isothermal self-gravitating clouds.He assumed that the ratio of magnetic to thermal pressure in the deep interior of clouds tends to remain constant and close to unity, β = 1: where c s the isothermal sound speed.Clearly, from Eq. 3, B ∝ ρ 1/2 .Summarizing, the two slopes, 2/3 (Mestel 1965) and 1/2 (Mouschovias 1976a,b), imply different kinds of compressions for the clouds: the former suggests an isotropic spherical geometry, while the latter implies a slab-like or filamentary geometry; in the 1/2 case, the magnetic field lines are either perpendicular to the slab or inclined relative to the primary axis of the filament.More generally, if the equation of state is P ∝ ρ γ , then B ∝ ρ γ/2 (Mouschovias 1991).In the following, we will be referring to the slope κ of the B − ρ relation by assuming B ∝ ρ κ .Since Zeeman measurements of the magnetic field in the atomic and molecular ISM became accessible, the B − ρ relation has also been thoroughly studied in observations.For example, Verschuur (1969) reports that κ = 2/3 based on Zeeman measurements of nine HI clouds, however, without giving an explicit fit.Later, Crutcher (1999) supported κ ≃ 1/2 by fitting a larger number of cloud data with number densities greater than 100 cm −3 , and finding κ = 0.47 ± 0.08.In a subsequent analysis with a larger number of atomic and molecular regions, Crutcher et al. (2010) displayed a log(B)-log(n) diagram that is flat at low densities, indicating no correlation between them, while for higher densities the data follow a trend close to B ∝ ρ 2/3 .However, Tritsis et al. (2015) revisited the Crutcher et al. (2010) data by examining different cloud morphologies; they found that the preferred slope is κ = 0.5 when the entire density range is considered.Li et al. (2015) fitted an even flatter slope, κ = 0.41 to observations of the massive star-forming region NGC 6334.By extending the Bayesian analysis from Crutcher et al. (2010) and by using recent observational and theoretical developments, Jiang et al. (2020) showed that κ cannot be reliably estimated when the observational uncertainty exceeds 2 due to the errors-in-variables bias.Further, Myers & Basu (2021) examined 17 dense cores using the Davis-Chandrasekhar-Fermi (DCF) method to estimate plane-of-sky field strengths B pos ; they obtained a best-fit value of κ = 0.66 for their dataset.Additionally, Liu et al. (2022) compiled DCF data and obtained a best-fit value of κ = 0.57, but they noted that the relation contains substantial uncertainties.
Despite the great progress of magnetic field observations in the ISM, we are far from fully understanding the connection between the magnetic field and galaxy evolution, mainly because it is inherently difficult to measure the magnetic field strength and to characterize its topology.In parallel with observational efforts, theorists have been trying to gain insights from numerical simulations.These simulations can provide the B − ρ relation and allow for a comparison with observational data.
Early on, Fiedler & Mouschovias (1992, 1993) performed numerical MHD simulations of isothermal, axially symmetric, gravitational contracting cores, predicting a slope of κ = 0.47.Kudoh et al. (2007) developed the first fully 3D simulation to study molecular cloud fragmentation; they found that for higher densities, κ tends to be 0.5.Collins et al. (2011) simulated supersonic, super-Alfvénic, and self-gravitating turbulence: throughout the collapsing gas, κ = 0.5 was found.Mocz et al. (2017) simulated supersonic, turbulent, isothermal, self-gravitating gas with a range of magnetic mean-field strengths; they showed that when the kinetic energy dominates over the magnetic pressure, κ = 2/3 implying an isotropic collapse while for a dominant large-scale magnetic field the collapse is anisotropic with κ = 0.5.
In more recent studies, Seta & Federrath (2022) investigated the turbulent dynamo in a two-phase medium, obtaining κ = 0.22 and κ = 0.27 for a solenoidal turbulent driving at low and high temperatures, respectively; for the compressive case, the slope was found to be κ = 0.71 for low temperatures and κ = 0.51 for higher temperatures, in all cases with a large scatter around the relation.In simulations of an isothermal turbulent medium collapse, Brandenburg & Ntormousi (2022) found that the scatter around the B − ρ relation was high enough to accommodate various exponents.In simulations of decaying turbulence, Auddy et al. ( 2022) simulated turbulent molecular clouds and they observed a distinct break density that separates a relatively flat, low-density regime from a power-law regime at higher densities.They reported that the transition density increases with increasing values of the Alfvén Mach number.
On galactic scales, numerical models so far reach contradicting results regarding the B − ρ relation.Wang & Abel (2009) modeled the formation and early evolution of disk galaxies with a uniform magnetic field and discovered a flat relation for low densities and a slope of 1/2 for higher densities.However, Pardi et al. (2017), focusing on kpc-scaled regions within a galactic disk, did not observe an increase in magnetic field with density.Meanwhile, Girichidis et al. (2018), adopting a similar setup, reported a slope of 2/3 for low-density gas and 1/4 for highdensity gas.Furthermore, conducting cosmological simulations of Milky Way-like galaxies with an initial uniform magnetic field, Ponnada et al. (2022) found results consistent with both slopes of 2/3 and 1/2.
Until now, studies investigating the B − ρ relation have primarily assumed a uniform magnetic field, with little exploration of different magnetic field topologies.Existing research focusing on kpc-and pc-sized regions consistently indicates that increasing magnetic field strength suppresses the dense gas fraction and overall star formation rate within the model (Myers et al. 2014;Federrath 2015;Pardi et al. 2017;Girichidis et al. 2018;Wurster et al. 2019).However, recent simulations have revealed that a stronger random magnetic field can actually lead to a faster collapse in dense regions (Brandenburg & Ntormousi 2022).
In this paper, we investigate the impact of a galaxy's initial magnetic field topology on the B − ρ relation over time, specifically exploring the possibility that different field morphologies induce different compression modes (κ values).To achieve this, we perform two galaxy simulations, each initialized with a different initial magnetic field: one with an ordered toroidal field and the other with a random magnetic field.
We describe the numerical code and the setup in Sect. 2. The results of our investigations are presented in Sect.3, a discussion in Sect.4, and the conclusions in Sect. 5.

Method and setup
We perform MHD simulations of two Milky Way (MW)-sized disk galaxies, accounting for self-gravity, star formation, supernova feedback, and chemistry.The initial conditions comprise a dark matter (DM) halo, a gaseous halo, a gas disk, and a stellar disk.The simulations achieve a spatial resolution of 24 pc in high-density regions, providing detailed information at that scale.

Magneto-hydrodynamics
The galaxy simulations were conducted with RAMSES, an Eulerian, magneto-hydrodynamics (MHD), Adaptive Mesh Refinement (AMR) code (Teyssier 2002;Fromang et al. 2006).The models include dark matter and stars, represented by collisionless particles, and a multiphase gaseous component, treated as a magnetized fluid.RAMSES solves the ideal MHD equations in the following form: where v is the velocity, ϕ the gravitational potential, Λ, and Γ the cooling and heating rates of the gas respectively.RAMSES solves the MHD equations using a Godunov scheme with Constrained Transport to guarantee that the solenoidality condition for the magnetic field (∇ • B = 0) is fulfilled.

Non-equilibrium chemistry
The complex chemistry of interstellar gas plays a vital role in its thermal evolution.Properly modeling the non-equilibrium chemistry of molecular hydrogen formation and dissociation also leads to a more accurate treatment of star formation (Valdivia et al. 2018;Decataldo et al. 2020).For this reason, we use a customized version of RAMSES that follows the nonequilibrium chemistry of H 2 formation and dissociation through the KROME package (Grassi et al. 2014).The implemented chemical network contains H, H + , H − , H 2 , H + 2 , He, He + , He ++ , and electrons and allows us to track the evolution and thermodynamics of ionized, atomic and molecular gas (Bovino et al. 2016), and it has been adopted both in galaxy (Pallottini et al. 2017) and molecular clouds simulations (Decataldo et al. 2020).Various heating and cooling processes regulate the thermal state of the gas (eq.6).In addition to typical atomic processes, our simulations also incorporate the cooling resulting from the presence of H 2 in the gas phase, and metal cooling.Here we have assumed solar metallicity for the ISM gas.
The radiation field can modify the chemical evolution by causing the ionization and dissociation of atoms and molecules.On-the-fly radiative transfer is not explicitly included in our simulations.Instead, we allow for a uniform radiation background of 1 G 0 where G 0 = 1.6 × 10 −3 erg cm −2 s −1 is the far Ultraviolet (UV) flux in the Habing band (6 -13.6 eV) normalized to the average MW value (Habing 1968).The spectral shape is that of a non-ionizing Draine (1978) distribution, which is appropriate to describe a MW-like interstellar radiation field, and we select ten photon bins to cover the energy range from 5 to 13.6 eV, Note that self-shielding of molecular hydrogen from the photodissociating Lyman-Werner radiation is included adopting a cellby-cell prescription following Richings et al. (2014).

Star formation and supernova feedback
Our model includes star formation and supernova feedback.The star formation recipe is based on the Schmidt-Kennicut relation (Schmidt 1959;Kennicutt 1998) following the H 2 density (Pallottini et al. 2017).During star formation, a fraction of the molecular mass of a cell transforms into a new stellar particle with a stellar mass M * ,c .This process occurs with an efficiency ϵ S F that we set to ϵ SF = 1%, which is reasonable for MW-like galaxies, and we kept it constant between the two models to facilitate their comparison.The random sampling from a Poisson distribution was originally implemented by Rasera & Teyssier (2006).To prevent the spurious formation of stellar particles, a minimum stellar mass of M * ,c = 10 3 M ⊙ is set, which represents the smallest stellar mass that can form in a cell.
The supernova feedback is simulated by injecting thermal energy into the 27 neighboring cells around the star particle.These 27 cells consist of the central cell containing the supernova, which we can call C i , the 6 cells that share a face with C i , 12 cells that share an edge with C i , and 8 cells that share a vertex with C i .This arrangement ensures that the explosion is represented in a numerically stable way across varying resolution levels.Each supernova produces 10 51 erg.In our simulations, 20% of the stellar population in a particle undergoes supernova explosions 3 Myr after the star particle's generation.The initial stellar disk, which is composed of older stars, does not contribute to the supernova budget.Note that radiation or wind feedback from the stellar particles is not included in our simulations.

Initial conditions
We generate the initial conditions for the dark matter, stars, and gas using the DICE code (Perret 2016).DICE can combine the different components of a galaxy and sample three-dimensional (3D) distributions of particles using a Markov chain Monte Carlo (MCMC) algorithm.The dynamical equilibrium is reached by solving the Jeans equations.
We prepare MW-like galaxies at redshift z = 0 with a virial velocity of 200 km/s, which corresponds to a total mass of M tot = 10 12 M ⊙ .The galaxy is composed of a DM halo (97.5% of M tot ), a thin stellar disk (1.425%), a gaseous disk (0.075%), and a gaseous halo (1%).The DM and gaseous haloes follow a pseudo-isothermal density profile with a scale length of 3 kpc and a radial density cut at 100 kpc.Each DM particle has a mass of 4.4 × 10 5 M ⊙ .The thin stellar disk follows a Miyamoto-Nagai profile (Miyamoto & Nagai 1975) with a scale length of 3 kpc and a radial density cut at 12 kpc while all the initial stellar particles have the same mass.The gaseous disk follows an exponential density model with a scale length of 4 kpc, a radial density cut at 15 kpc and a constant density vertical profile cut at 0.75 kpc.The initial temperature is 8000 K without an initial turbulent velocity field.
For these simulations we use an AMR with a coarse resolution of 128 3 and five levels of refinement, by adopting the following strategy.We begin with a geometry-based refinement, where the two first levels are triggered in two cylindrical regions of decreasing size that are centered on the galaxy and fully encompass the gaseous disk with densities n/cm −3 > 10 −3 at all times.Inside the inner region, three additional levels of refinement are triggered by a Jeans-based criterion, i.e. we require the local Jeans length to be resolved with at least 10 cells.This strategy results in the highest effective resolution of 4096 3 in the densest regions.The simulation box is a cubic volume with a side length of 100 kpc, so the maximum resolution corresponds to 24 pc.
It is important to note that star formation is enabled at around 100 Myr for both models, to allow the gas to collapse and form dense clumps in the disk.The simulations run for a total duration of 500 Myr.
The initial conditions of the two galaxies are identical, except for the magnetic field morphology.One galaxy, referred to as model T, is initialized with a toroidal magnetic field, while the other galaxy, referred to as model R, starts with a random magnetic field.In both models, the magnetic field has an exponentially declining profile.The initial magnetic field topologies are shown in a 3D representation of magnetic field vectors in Fig. 1 for both models.
The initial magnetic field is generated by specifying the vector potential A, and then calculating the magnetic field by taking the curl of A, B = ∇ × A. This method ensures that ∇ • B = 0 to machine precision.For model T the definition of A follows: where êz is the unit vector along the z axis, r represents the cylindrical radius, and z is the vertical height from the galactic plane.
The parameters R 0 and H 0 represent the scale length and scale height of the magnetic field, respectively.In this work, the values of R 0 and H 0 are set to 1 kpc, and the magnetic field strength is normalized such that it is 10 µG at the galactic center.
For model R, we create a vector potential in Fourier space using a power spectrum: where i = x, y, z represents the component of A. Each component of A is then obtained in physical space through an inverse Fast Fourier Transform.The resulting magnetic field is normalized to match the central strength of model T and is also convolved with the same exponential profile of radius and height as model T.

Results
We study the evolution of the morphology of the galaxies, their gas and stellar content, and their star formation histories (Sec. 3.1).Then, we examine the thermodynamical and magnetic state of the gas, in particular investigating the plasma beta relation in different gas phases (Sec.3.2).Finally, we focus on the time evolution of the B-ρ relations (Sec.3.3).

Morphology and star formation rates of galaxies
Fig. 2 shows face-on total gas column density maps and projected magnetic field vectors for the two galaxies at different times (t = 200, 300, and 500 Myr).Model T is shown on the left-hand side, and model R is on the right-hand side.The magnetic field vectors are color-coded according to their strength.At early times, both models exhibit stronger magnetic fields in their central regions, as set in the initial conditions.However, as time progresses, we observe distinct differences in the magnetic field evolution between the two models.Specifically, model R displays a wider region of strong magnetic field, spanning approximately 0-5 kpc.In contrast, model T shows a narrower range, covering approximately 0-3 kpc, as visually inspected.
As we outlined in the introduction, star formation is sensitive to the strength and morphology of the magnetic field.Therefore, the sharp difference in magnetic energy distribution between the two models that we notice on visual inspection could in principle affect the star-forming properties of the galaxies.In Fig. 3 we plot the star formation rate (SFR, left-hand side) and the cumulative mass of the new stars (right-hand side) as a function of time.The SFR history is generated with a time binning of 10 Myr.We have checked that varying this time interval does not affect the resulting trend of SFR with time.The red dashed line indicates the moment where we switch on star formation in the simulation.
We observe some differences between the two galaxies.In model T, the peak of the star formation rate occurs at around 250 Myr and reaches approximately 5.5 M ⊙ /yr; instead, in model R, the peak of the star formation rate occurs later, at approximately 380 Myr, with a higher value of around 6.5 M ⊙ /yr.Despite the temporal offsets observed in the SFRs between the two models, both ultimately result in a similar total mass in stars, approximately around 10 9 M ⊙ while the total stellar mass from the initial conditions was 1.4 × 10 10 M ⊙ .Furthermore, both models exhibit a comparable average star formation rate, 2.19 ± 1.69 M ⊙ /yr for model T and 2.32 ± 1.95 M ⊙ /yr for model R.These SFR values align with observational estimates of the SFRs in Milky Waylike galaxies, (Licquia & Newman 2015;Fraser-McKelvie et al. 2019;Boardman et al. 2020;Elia et al. 2022); it is worth noting that Galactic observations also indicate temporal variations in the SFR (Lee et al. 2016).
In Fig. 4, we show the surface density of the SFR at the time when each model reaches the peak of its SFR, as obtained from Fig. 3 (250 and 380 Myr, for models T and R respectively).Overall, star formation is more intense in the central 0.5-1 kpc of the galaxies and in spiral structures at larger radii reaching peaks of 1.41 and 1.85 M ⊙ yr −1 kpc −2 for model T and model R at 250 and 380 Myr, respectively.The corresponding mean surface density of the star formation rate is 8 × 10 −4 and 7.3 × 10 −4 M ⊙ yr −1 kpc −2 .At 250 Myr, star formation in both models is mainly organized in clumps, however, model T includes a smooth central region of about 5 kpc of high star formation rate.

Thermodynamical and magnetic state of the gas
In order to study the dynamical state of the gas in different phases, we first examine the plasma β as a function of gas density and time.The relation between β and mean number density n is shown in the top panel of Fig. 5 for the final snapshot of each model at 500 Myr.In this figure, the dots represent the mean values of log β in each density bin.To produce these plots we take into account only the gas cells with |z| ≤ 1.5 kpc and radius R ≤ 13 kpc encompassing the entire disk. 1 We study the behavior of the gas by categorizing it into three distinct phases defined by number density: low-density (10 −2 < n/cm −3 < 1), intermediate-density (1 < n/cm −3 < 10 2 ), and high-density phase (10 2 < n/cm −3 < 10 4 ).
In model T, the gas in the low-density phase exhibits high values of β; at intermediate and high densities, the mean behavior suggests a transition toward equipartition.Furthermore, the regions with high mass bins are predominantly located in the magnetically dominated region (log(β) < 0).
In model R, the mean behavior in the intermediate and highdensity phases points to magnetically dominated gas.Interestingly, there is also a substantial fraction of very magnetically dominated gas (β < −2.5) in the density range from 10 −2 to 10 cm −3 .Note that the initial conditions consist of only thermally dominated gas with low-density a temperature of T = 8000 K.However, after 100 Myr, the gas is able to cool and condensate at intermediate densities, becoming magnetically dominated.
In order to understand the origin of this behavior, in the middle panel of Fig. 5 we plot the total pressure (log(P/k B )). Here, the blue and green dots represent the mean values of the thermal and magnetic pressure per density bin, respectively2 .
In the low-density medium of model T, thermal pressure dominates and the signature of the thermal instability (Field 1965) is evident from the characteristic bend of the average pressure line (Wolfire et al. 1995).The two contributions to pressure are equal in the higher-density regions.In model R, the magnetic pressure dominates over almost the entire density range, with the exception of very low-density gas (n < 10 −2 cm −3 ).As a result, the two atomic phases (warm neutral medium (WNM, −1 < log(n/cm −3 ) < 0) and cold neutral medium, (CNM, 1 < log(n/cm −3 ) < 2)) created by the thermal instability are both supported by magnetic pressure.
This magnetic support yields a visible signature in the gas temperature, which is shown 3 as a function of density in the bottom panel of Fig. 5. Model R contains a significant gas fraction at low densities (10 −1 < n/cm −3 < 10) with temperatures lower than 100 K.This phase is absent in the temperature-density plot of model T. The additional magnetic pressure in model R allows low-density gas to cool while still maintaining pressure equilibrium with its surroundings.
In order to study the time evolution of the thermal state of the gas, in Fig. 6 we plot ⟨log β⟩ (top panels) and the total thermal and magnetic energies (bottom panels) as a function of time for the low, intermediate, and high-density phases4 .
The blue line represents the low-density gas, the orange line represents the intermediate-density gas, and the green line represents the high-density gas for model T (left-hand side) and model R (right-hand side).In the top plots, the black dashed line corresponds to log(β) = 0, while the red dashed line represents the onset of star formation.The reason why we show the calculation of β and the different energies only after approximately 150 Myr is that, up to that point, the gas has not yet cooled enough for the intermediate-and high-density gas to be present in a stable state.In the initial condition, ⟨log β⟩ = 6.3, and the total thermal and magnetic energies are 10 56.6 and 10 54.4 ergs, respectively.
In both models, the low-density gas presents high values of ⟨log β⟩ indicating thermal dominance throughout the model evolution.However, there is a decrease of β with time because the gas is cooling.The plasma β in Model R is slightly lower than in model T, which is expected since we showed that the magnetic field of model R can be strong in regions with low-density gas (see Fig. 2).Calculating the temporal average and standard deviation, we find ⟨log β⟩ T = 3.1±1.8which is an order of magnitude higher than the temporaly averaged value of β in the low-density gas of model R, ⟨log β⟩ R = 2.3 ± 2.1.
The intermediate-density gas starts in a state of thermal dominance with a value of ⟨log β⟩ around 0.6 in both models, which is much lower than the initial values of the low-density gas.As the models evolve, ⟨log β⟩ reaches an equipartition value in model R, while in model T, it remains in the thermally dominant region, albeit close to equipartition.Calculating the temporal average and standard deviation, we find a higher value of β in model T ⟨log β⟩ T = 0.96±1.7.compared to model R (⟨log β⟩ R = 0.1±1.5).
The high-density gas is present starting from 180 Myrs in model T, while in model R, it forms earlier, at 120 Myrs.In both models, ⟨log β⟩ starts with lower values (-1.3 in both cases), indicating magnetic dominance.However, in model T, ⟨log β⟩ increases over time, finally reaching equipartition, whereas in model R, it remains in the magnetically dominated region throughout the evolution.However, they exhibit similar temporally-averaged values, ⟨log β⟩ T = −0.16± 1.2 for model T and ⟨log β⟩ R = −0.72 ± 1.1 for model R. We note the large scatter around all the above averages, which mainly comes from the large scatter in the individual β values.With this large scatter in mind, the average values indicate a stronger dynamical signif-  (blue), the intermediate-density (orange), and the high-density gas (green).In the initial conditions, the total thermal energy of the low-density gas in both models is two orders of magnitude higher than its total magnetic energy.However, by 500 Myr, these energies tend to become equal in model T. On the other hand, in model R, the total magnetic energy eventually surpasses the thermal energy.For both models, the total magnetic energy of the intermediate and high-density gas remains higher than the thermal energy throughout the entire evolution, except for the time period before 100 Myr, where the energies are approximately equal in the intermediate phase.
It is important to note that the ⟨log β⟩ values do not perfectly match the total energies due to several factors.Firstly, the ratio of total energies does not directly translate into the mean of log(β), as they represent different physical quantities and have distinct interpretations.Moreover, the high standard deviation observed in the data contributes to the tension.The variability in log(β) within the density ranges can lead to a wide spread of values, affecting the mean and making it less representative of the entire gas population.
We also present one-dimensional (1D) radial profiles of thermal and magnetic energies for the three different gas phases at three different times (200, 300, and 500 Myrs) in App.B. These profiles also illustrate that, model R consistently exhibits a higher degree of magnetic dominance than model T across different densities, and over time.They also show a slightly wider distribution of dense gas in model R compared to model T.

B-ρ relations
Fig. 7 shows the magnetic field-number density (B-n) relation at t=470 Myr.Model T is on the left-hand side and the model R is on the right-hand side.As in Fig. 5 and 6, we take into account only the gas cells in the disk, with |z| ≤ 1.5 kpc and radius R ≤ 13 kpc.The B-n relations are given as 2D massweighted histograms in the log B-log n space.The colored dots indicate the mean magnetic field strength per density bin.The standard error of the mean is not noticeable here because it is very small.We fit the mean relation using a broken power law: where the two power-law slopes κ 1 and κ 2 are parameters of the fit, while the transition number density, n 0 , is fixed in this instance to 100 cm −3 .In this figure we show only one snapshot to showcase the break, but later in this section we will examine the impact of varying the parameter n 0 and show the time evolution of the slopes.The broken power-law fit was chosen to test the suggestion that the B-ρ relation should be flat below, and a power-law above, a critical density n 0 (Crutcher 1999;Pattle et al. 2023).However, even a simple visual inspection of Fig. 7 indicates that a single power law is a better fit to the B-ρ relations in both models.Therefore, in the following we will also present results for a single power-law fit with an index simply referred to as κ.
The slopes κ 1 and κ 2 are provided in the legend of each plot.Additionally, in the bottom right corner of each plot, we include the B-ρ relations proposed by various theories, which are based on different assumptions regarding conservation laws and geometry during the collapse (as discussed in Section 1).We recall here that B ∝ n 0 results from compression along magnetic field lines, B ∝ n 1 from compression perpendicular to the magnetic field lines, B ∝ n 1 2 implies a slab-like or filamentary geometry, where the magnetic field lines are either perpendicular to the slab or inclined relative to the primary axis of the filament, which collapses radially and B ∝ n 2 3 arises from spherical compression (see Fig. 1 from Tritsis et al. 2015, for a clear illustration).
In both models, there is a positive correlation between the magnetic field and density.In the intermediate-and high-density gas phase, we do not observe the flatness in the magnetic fields reported in observations for n < 300 cm 3 (Crutcher et al. 2010;Pattle et al. 2023).However, when we focus on the peak of the mass, represented by darker regions, at diffuse medium densities we observe a bimodal distribution.One branch appears to extend the high-density power law to lower densities, while the other branch shows a flattening of the magnetic field close to 10 µG.This behavior has been consistent across all the times we have examined.
Comparing the two panels, for both models we see that the B-n relation at intermediate densities has a power-law slope of κ 1 = 0.31 ± 0.01.For both models, in the high-density medium, the power law steepens with respect to the intermediate, with κ 2 = 0.61 ± 0.02 for model T and κ 2 = 0.53 ± 0.01 for model R.According to Mestel (1965); Mouschovias (1976a,b); Tritsis et al. (2015), the dense gas behavior of model T is closer to the theoretical prediction for isotropic spherical compression while the behavior of model R implies a slab or filament-like geometry.
However, it is crucial to highlight that this is for a single snapshot in time and that the unique aspect of model R does not follow the conventional assumption of a uniform magnetic field.Therefore, the existing theoretical predictions may not directly capture the specific dynamics and behaviors observed in model R. Using the same fitting procedure for the B-n relation, in Fig. 8 we plot the slopes (κ 1 , κ 2 ) for the broken power-law fit and the slope κ for the single power-law fit as a function of time for Model T (blue) and Model R (green).
In order to examine the dependence of the fitted exponents on the value of the transition density we consider three different values: n 0 = 50 cm −3 , shown in the top panel, n 0 = 100 cm −3 , shown in the middle panel, and n 0 = 300 cm −3 , shown in the bottom panel.We chose the particular transition densities because according to Crutcher et al. (2010) we expect the break of the power law at n 0 = 300 cm −3 .However, we investigate additional density ranges, as the break may occur at different values depending on time.The shaded areas correspond to the standard deviation (1-σ) of the slopes derived from the fit.The black and red dashed lines represent the slopes of 1/2 and 2/3, respectively.It is worth noting that both slopes vary significantly over time.
In Model T and for a transition density of n 0 = 50 cm −3 , the intermediate-density slope κ 1 fluctuates around 2/3 from 200 to 300 Myr.Later on, it decreases to values lower than 1/2, indicating compression predominantly along the magnetic field lines.Instead, in model R, κ 1 fluctuates around 0.4 for the entire duration of the simulation.Eventually, both models converge to similar values around 0.3-0.4.For higher transition densities (n 0 = 100 cm −3 and n 0 = 300 cm −3 ), the behavior of κ 1 is similar to that of n 0 = 50 cm −3 , but the variations are smoother.
For model R, and a transition density of n 0 = 50 cm −3 , the high-density slope κ 2 fluctuates around 1/2, implying a slab or filament-like geometry with a perpendicular or radial compression.At the same transition density, κ 2 for model T is slightly lower, staying slightly below 1/2 for almost the entire duration of the simulation.Calculating the temporal average of the slopes for n 0 = 50 cm −3 , we find ⟨κ 1 ⟩ T = 0.34 ± 0.19 for model T and ⟨κ 1 ⟩ R = 0.38 ± 0.09 for model R. Model T also exhibits ⟨κ 2 ⟩ T = 0.44 ± 0.07, and model R records an average slope of Model T Model R Fig. 6.Evolution of β⟩ (top) and total energies (bottom) for three different phases of the gas.The blue, orange, and green lines show the gas in a low (−2 < log(n/cm −3 ) ≤ 0), medium (0 < log(n/cm −3 ) ≤ 2) and high (2 < log(n/cm −3 ) ≤ 4) density phase, respectively, for model T (left-hand side) and model R (right-hand side).In the top plots, the black dashed line corresponds to ⟨log β⟩ = 0, which is the critical value.Above this value, the gas is thermally dominated, while below it, it is magnetically dominated.For the bottom plots, circles represent the thermal energy while triangles correspond to the magnetic energy.The red dashed line represents the onset of star formation.The shaded regions correspond to the standard deviation.
⟨κ 2 ⟩ R = 0.48 ± 0.08.The errors are attributed to the time evolution which is the primary source of variance.For higher n 0 values, κ 2 exhibits significant variations and differences between the two models.
For a transition density of n 0 = 100 cm −3 , κ 2 starts off near 0.2 for model T and 0.6 for model R.After about 250 Myr and for the remainder of the simulation, κ 2 in both models R and T fluctuates around 1/2 but with larger amplitude (between 0.4 and 2/3) than for the lower transition density.Calculating the temporal average of the slopes, ⟨κ 1 ⟩ T = 0.36 ± 0.16 for Model T and ⟨κ 1 ⟩ R = 0.40 ± 0.07 for Model R. For higher densities (> n 0 ), Model T exhibits a slope of ⟨κ 2 ⟩ T = 0.46 ± 0.15, and Model R ⟨κ 2 ⟩ R = 0.51 ± 0.10.
In the case of n 0 = 300 cm −3 , the initial values of κ 2 for Model T are close to zero, indicating compression along magnetic field lines, which is the theoretical expectation for an ordered field morphology initially.After 280 Myr, κ 2 in model T shifts to unity, the theoretical expectation for compression across the field lines.At later times, κ 2 fluctuates strongly between 0.5 and 1.Also for R κ 2 fluctuates strongly, although not following any of the theoretical expectations.The stronger fluctuations as the threshold density increases are due to the smaller fitting range, but also due to poorer statistics, because at this density range, we are occasionally depleting gas to star formation.
Overall, the B − ρ relation exhibits strong variations across different density ranges and over time.However, despite these variations, the average intermediate slope (⟨κ 1 ⟩) consistently appears lower than the average high-density slope ( ⟨κ 2 ⟩).Additionally, for higher densities, both slopes tend to converge toward 1/2 for lower transition densities (n 0 ) which gives better statistics.This consistency aligns with previous studies by Mouschovias (1976a,b).
The single power-law fit (bottom row of Fig. 8) is very similar in temporal behavior to that of κ 1 .The slopes in the two models both start near κ=0.5 and decrease over time, to a value close to κ=0.4.While there are time variations of the same order as what we noted for κ 1 in the previous cases, when averaged Model T Model R over time, the κ values in the two models are within one sigma from each other.Finally, it is worth noting that the χ 2 for the single power-law fit is lower than that of the double power law, for both models.We calculated the BIC (Bayesian Information Criterion, Schwarz 1978) for the two models as an indication for the preferred fit.The BIC: where k the number of parameters of the model, n the number of data points, and L the likelihood function, penalizes models with more fitting parameters, with lower values of BIC indicating a better fit.Here as L we have taken the value of χ 2 ).The BIC comparison favors the single power-law fit with respect to the broken power law for all snapshots, despite the small difference in the number of parameters between the two models.

Summary
In this work, we presented galaxy-scale MHD simulations that incorporate gravity, star formation, supernova feedback, and chemistry.For the first time, we explored different magnetic field morphologies, specifically an ordered toroidal magnetic field (model T) and a random magnetic field (model R).Our study focuses on star formation, the plasma beta parameter, and the B − ρ relation in the two models.We observe that the star formation rates (SFRs) of the two models peak at different times, and there are some differences in the instantaneous SFR at different instants.However, the timeaveraged SFR is similar between the two models, especially when taking into account the large scatter around the temporal mean.They also produce the same stellar mass after 500 Myrs of evolution.This implies that a different magnetic field morphology may result in a slightly different evolution of the SFR due to the complex nature of the problem, but this behavior does not affect the average star formation history of the galaxy.
We found that in both models, the low-density gas is thermally dominated, with model R showing slightly lower ⟨log β⟩ values.The intermediate-density gas tends toward equipartition in model R but remains thermally dominant in model T, while the high-density gas forms earlier in model R and exhibits magnetic dominance throughout its evolution.Concerning the total thermal and magnetic energies, initially, the total thermal energy of low-density gas significantly exceeds its total magnetic energy in both models.However, by 500 Myr, these energies approach equality in model T, whereas in model R, the total magnetic energy eventually surpasses the thermal energy.Throughout the evolution, both models maintain higher total magnetic energy for the intermediate and high-density gas.However, it is important to note that the scatter around the average β values is always very large, so the above differences are within the errors.Regarding the B − ρ relation, we observe a scaling that follows a power law with exponents κ ranging from 0.2 to 0.8 and varying over time.The temporal evolution of κ differs between the two models, but these disparities fall within the wide temporal scatter.

Comparison with previous studies and related works
Previous computational works have studied the B − ρ relation using different approaches and assumptions.Pardi et al. (2017) focused on kpc-sized portions of a galactic disk, examining the chemical, thermal, and dynamical evolution of the gas with self-gravity and supernova feedback.They found a large scatter around the B − ρ relation due to the dominance of kinetic energy density in their simulations.Unlike our results, they did not find a clear increase in magnetic field strength with density across the entire density range, suggesting a possible flattening of the relation above a density of 0.6 cm −3 .This flattening might be a result of a limited dynamical range.
Using a similar setup, Girichidis et al. (2018) found that weakly magnetized low-density gas (n < 0.6 cm −3 ) exhibited a power-law scaling with a slope of κ = 2/3, while higher densities (n > 0.6 cm −3 ) showed a flatter of κ = 1/4.Our results do not align with these studies since we observe scaling for densities higher than 10 cm −3 , and moreover, the power law becomes steeper at higher densities.However, there are some differences in the setups between our study and the previous ones by Pardi et al. (2017); Girichidis et al. (2018).Notably, the previous studies focus on smaller-scale regions, do not simulate entire galaxies, and utilize initially uniform magnetic fields.Additionally, variations in the heating and cooling mechanisms as well as the star formation recipe employed could contribute to the differences observed.Seta & Federrath (2022) modeled the turbulent dynamo in a two-phase medium considering solenoidal and compressive driving.They initialized the simulation with a weak random seed field and explored a density range of 10 −2 to 10 4 cm −3 , which varied depending on the turbulence driving and temperature.They find that the slope of the B − ρ relation is shallower in the solenoidal (κ = 0.22 for T < 10 3 K, κ = 0.27 for T > 10 3 K ) than in the compressive case (κ = 0.71 for T < 10 3 K, κ = 0.51 for T > 10 3 K ).In our study, we focus on fitting the B − ρ relation for temperatures lower than 10 4 K. Within the uncertainties due to the temporal variations, our models lie between these two cases.This fact is not surprising because turbulence in our simulation is injected both from supernova events, which inject a mixture of compressive and solenoidal modes (Pan et al. 2016), and the differential rotation of the disk, which preferentially injects solenoidal modes.However, a simple calculation indicates that, in our models, the dominant contribution to turbulence comes from the differential rotation.The number of supernova explosions can be estimated from the number of stars formed per unit time, indicatively with a ⟨S FR⟩ ≃ 1 M ⊙ /yr.With 20% of this mass exploding as SNe, E SN = 10 51 erg the energy per SN and 25% of E SN going to kinetic energy injected in the ISM, we estimate the rate of turbulent kinetic energy from SNe to be approximately E turb,S N ≃ 10 49 erg/yr.For the energy contribution from the differential rotation, we calculate the rotational velocity for each cell and use the relation K rot = 1/2 m v 2 rot to calculate the kinetic energy.We then divide this energy by the timestep, which is 10 3 yr, to obtain an estimate of E turb,rot ≃ 10 54 erg/yr.It is important to note that cooling effects are not included in this calculation, and cooling could counteract some of the energy injection from supernovae.
Finally, in collapse simulations of an isothermal turbulent medium, Brandenburg & Ntormousi (2022) employed both a random and a guide field with low and high magnetic field strengths.They consistently observed a large scatter around the relation, such that both κ = 2/3 and κ = 1/2 could fit the data.This finding is consistent with our study, where we also observe significant scattering in the relation.

The impact on galactic physics and observations
Variations in the B − ρ relation with time and with magnetic field morphology could have significant implications for our interpretation of the dynamical state of the gas and the star formation process.Our analysis shows that the B − ρ relation is not universal but rather context-dependent, reflecting the interplay between magnetic fields and the evolutionary stage and magnetic field morphology of the observed galaxy.
In principle, a galaxy's magnetic field will change morphology as the galaxy evolves.Mergers and accretion, which are particularly relevant for early galaxies (e.g.Kohandel et al. 2020), drive turbulence, which in turn can generate a turbulent dynamo.This would lead to a more random magnetic field morphology.
Instead, a mean-field dynamo could take over in more quiescent phases, leading to more coherent large-scale fields (e.g., Brandenburg & Ntormousi 2023, for a review on galactic dynamos).However, our results imply that any influence of the galactic magnetic field's morphology on the B-ρ relation must be smaller than the scatter introduced by the non-linear processes in the ISM, since the differences between the two models always fall within a very broad temporal scatter.However, we should stress that our experiments do not follow the complex, long-term evolution of the magnetic field from primordial seeds to the present day.Both models contain strong initial fields and none of them shows dynamo action.The imprints of the complex co-evolution of the galaxy and its magnetic field will be the subject of future work.Our results also suggest that interpreting the B − ρ relation in observations should be done cautiously, given the large scatter and the complex dependencies discovered in our models.While our study acknowledges that we are currently lacking proper production of observable-like data, it still emphasizes the importance of considering various factors and specific galactic conditions in comprehensively understanding the role of magnetic fields in star formation.

Future improvements: early feedback
One of the innovations in our work is the inclusion of nonequilibrium chemistry in the simulations through the KROME package, allowing us to follow the formation and dissociation of H 2 .This leads to a more realistic representation of star formation because it is based on the H 2 fraction.However, it is important to note that our simulations currently lack pre-SN feedback.Including early feedback, namely radiation and stellar winds alter their impact (Pallottini et al. 2017(Pallottini et al. , 2019;;Decataldo et al. 2020), in particular increases the time variability of the SFR (Pallottini & Ferrara 2023), which in turn might enhance the scatter in the B − ρ relation for some gas phases.The inclusion of early feedback is clearly the next step for this study.

Conclusions
Our analysis of galaxy-scale MHD simulations with different magnetic field morphologies contributes to the understanding of the B − ρ relation as a probe of the ISM's dynamical state with several important findings: 1.The two models (with ordered and random magnetic field) have a similar average star formation rate and they end up forming roughly the same total mass in stars, implying that the evolution of galaxies with comparable formation histories can have unique ISM characteristics depending on the magnetic field morphology.2. The plasma beta values indicate different thermal and magnetic dominance behaviors across different gas density phases.Both models have thermally dominated low-density gas (10 −2 < n/cm −3 < 1) with ⟨log β⟩>3, but model R shows slightly lower ⟨log β⟩ values.In the intermediate-density gas (1 < n/cm −3 < 10 2 ), model R tends toward equipartition, while model T remains thermally dominant with ⟨log β⟩ ≈ 1.
For high-density gas (10 2 < n/cm −3 < 10 4 ), model R exhibits magnetic dominance throughout its evolution with ⟨log β⟩ ≈ −1, whereas in model T high-density gas forms later and reaches equipartition.3. The total magnetic energy remains higher than the thermal energy for the intermediate and high-density gas phases in both models throughout their evolution.However, in the case of low-density gas, the total thermal energy initially exceeds the total magnetic energy both models.By 500 Myr, these energies approach equality in model T as the galaxy cools down and the magnetic field strength slightly increases.On the other hand, in model R, the total magnetic energy eventually surpasses the thermal energy due to a slightly stronger magnetic field compared to model T. 4. The B − ρ relation is not universal even within the same galaxy, displaying bimodal distributions and significant variations over time with κ ranging from about 0.2 to 0.8.Despite the variations, the slope κ at densities n > 50, cm −3 tends to converge to 1/2 for both models.5.Even considering the large scatter, a single power-law fit describes the B-ρ relation at intermediate and high densities (1 < n/cm −3 < 10 3 ) better than a broken power law, independently of the choice of the transition density, n 0 .6.The observed differences in the slope evolution between the two models are not large enough to indicate that the magnetic field morphology influences a galaxy's B−ρ relation because they fall within the very large scatter.
Overall, we see a small influence of the magnetic field morphology on a galaxy's B−ρ relation, which is, however, transient and below the level of the relation's fluctuations due to other stochastic phenomena.Eventually, both models tend to B ∝ ρ 0.5 .In general, the differences in between the two models, while measurable at each individual time, show significant scatter over time, highlighting the complex and context-dependent nature of the B − ρ relation.In closing, we should emphasize that, due to limited resources, our study only included two models of gaspoor, massive spiral galaxies.The effects of the magnetic morphology on the gas dynamics could be very different for gas-rich, more turbulent, or less massive galaxies.Further research and comprehensive observations will be crucial to fully understand the intricate role of magnetic fields in star formation processes across diverse galactic environments.

Fig. 1 .
Fig. 1. 3D representation of the initial conditions (t = 0) for the magnetic field vectors for the two simulations, showing the toroidal (T) and random (R) cases on the left and right side, respectively.For illustration purposes, we show only vectors in the central region of the simulation.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 2. Face on maps of the gas column density (N gas ) and the projected magnetic field vectors (B).Model T and model R are shown on the left-hand side right-hand side, respectively.From top to bottom we plot the maps at t = 200, 300, and 500 Myr.

Fig. 7 .
Fig. 7. |B| versus n for the two models at t = 470 Myr with the model T (left-hand side) and the model R (right-hand side).The B-ρ relations are given as 2D mass-weighted histograms in the log B-log n space.The dots indicate the mean magnetic field strength per density bin.The yellow and red lines are the power-law fits with slopes given in the legend of each plot.The B-ρ relations derived from Mestel (1965); Mouschovias (1976a,b); Tritsis et al. (2015) are reported with dashed black lines.

Fig. 8 .
Fig. 8. Evolution of the B − ρ slopes (κ 1 , κ 2 ) and the slope κ from the single power-law fit as a function of time for model T (blue) and model R (green) for intermediate (left-hand side) and high densities (right-hand side).The different transition densities for the broken power-law fit are indicated as labels on the left.The definition of the broken power law fit is given in Eq. 10.The shaded areas correspond to the standard deviation error on the slopes derived from the fit.The black and red dashed lines show the slopes of 1/2 and 2/3 respectively.