© The Authors 2024

1 Introduction

Magnetic fields play fundamental roles in 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 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 thermal pressure. Plasma β represents the ratio of thermal pressure, P_th, to magnetic pressure, P_mag; 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 Mg/Φ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: $$M_{\text{g}}\propto \rho R^3=\text{const,}$$(1) $${\Phi }_{\text{B}}=B{R}^2=\text{ const, }$$(2)

where M_g is the gas mass, R the cloud radius, Φ_B the magnetic flux, and B the magnetic field, which is assumed to be 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: $$\frac{P_{\text{b}}}{P_{\text{th }}}=\frac{B^2}{8\pi \rho c_{\text{s}}^2}\simeq 1\text{, }$$(3)

where c_s the isothermal sound speed. Clearly, from Eq. (3), B ∝ ρ^1/2. To summarize, 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, when we refer to the slope, κ, of the B − ρ relation, we are assuming B ∝ ρ^K.

Since Zeeman measurements of the magnetic field in the atomic and molecular ISM became accessible, the B − ρ relation has been thoroughly studied in observations. For example, Verschuur (1969) reports that κ = 2/3 based on Zeeman measurements of nine HI clouds, though without giving an explicit fit. Later, Crutcher (1999) provided support for the κ ≃ 1/2 by fitting more cloud data with number densities greater than 100 cm⁻³ and finding κ = 0.47 ± 0.08. In a subsequent analysis with a larger number of atomic and molecular regions, Crutcher et al. (2010) present a log(B)-log(n) diagram that is flat at low densities, indicating no correlation, 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) fit 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) show that κ cannot be reliably estimated when the observational uncertainty exceeds 2 due to the errors-in-variables bias. Furthermore, 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 in 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 magnetohydrodynamic (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; for a dominant large-scale magnetic field, the collapse is anisotropic, with κ = 0.5.

In a more recent study, 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 observed a distinct break density that separates a relatively flat, low-density regime from a power-law regime at higher densities. They report that the transition density increases with increasing values of the Alfvén Mach number.

On galactic scales, numerical models have so far reached 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 kiloparsec-scale regions within a galactic disk, did not observe an increase in magnetic field with density. Meanwhile, Girichidis et al. (2018), adopting a similar setup, report a slope of 2/3 for low-density gas and 1/4 for high-density gas. Furthermore, conducting cosmological simulations of Milky Way (MW)-like galaxies with an initial uniform magnetic field, Ponnada et al. (2022) found results consistent with slopes of both 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 kiloparsec- and parsec-sized regions consistently indicates that increasing the magnetic field strength suppresses the dense gas fraction and overall star formation rate (SFR) within a 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 performed 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. We present the results of our investigations in Sect. 3, a discussion in Sect. 4, and the conclusions in Sect. 5.

2 Method and setup

We performed MHD simulations of two 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.

2.1 Magnetohydrodynamics

The galaxy simulations were conducted with RAMSES, an Eulerian MHD adaptive mesh refinement (AMR) code (Teyssier 2002; Fromang et al. 2006). The models include DM 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: $$\frac{\partial \rho }{\partial t}+\nabla (\rho \mathit{\upsilon })=0,$$(4) $$\frac{\partial (\rho \mathit{\upsilon })}{\partial t}+\nabla (\rho \mathit{\upsilon }\mathit{\upsilon }-\mathit{B}\mathit{B})+\nabla P_{\text{tot }}=-\rho \nabla \varphi ,$$(5) $$\frac{\partial E_{\text{tot}}}{\partial t}+\nabla \left(\left[\left(E_{\text{tot}}+P_{\text{tot}}\right)\mathit{\upsilon }-(\mathit{\upsilon }\cdot \mathit{B})\cdot \mathit{B}\right)\right]=-\mathit{\upsilon }\cdot \nabla \varphi -\rho \text{Λ}+\text{Γ,}$$(6) $$\frac{\partial \mathit{B}}{\partial t}-\nabla \times (\mathit{\upsilon }\times \mathit{B})=0,$$(7)

where υ is the velocity, ϕ the gravitational potential, and Λ 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.

2.2 Nonequilibrium chemistry

The complex chemistry of interstellar gas plays a vital role in its thermal evolution. Properly modeling the nonequilibrium 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 used a customized version of RAMSES that follows the nonequilibrium chemistry of H₂ formation and dissociation through the KROME package (Grassi et al. 2014). The implemented chemical network contains H, H⁺, H⁻, H₂, ${\text{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₂ 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 allowed for a uniform radiation background of 1 G₀, where G₀ = 1.6 × 10⁻³ erg cm⁻² s⁻¹ is the far-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 selected ten photon bins to cover the energy range from 5 to 13.6 eV. We note that self-shielding of molecular hydrogen from the photo-dissociating Lyman-Werner radiation is included adopting a cell-by-cell prescription following Richings et al. (2014).

2.3 Star formation and supernova feedback

Our model includes star formation and supernova feedback. The star formation recipe is based on the Schmidt-Kennicutt relation (Schmidt 1959; Kennicutt 1998) following the H₂ 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 ϵ_sf 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³ 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⁵¹ 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. We note that radiation or wind feedback from the stellar particles is not included in our simulations.

2.4 Initial conditions

We generated the initial conditions for the DM, stars, and gas using the DICE code (Perret 2016). DICE can combine the different components of a galaxy and sample 3D distributions of particles using a Markov chain Monte Carlo (MCMC) algorithm. The dynamical equilibrium is reached by solving the Jeans equations.

We prepared 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¹² 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⁵ 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 used an AMR with a coarse resolution of 128³ and five levels of refinement, adopting the following strategy. We began 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⁻³ > 10⁻³ at all times. Inside the inner region, three additional levels of refinement are triggered by a Jeans-based criterion, that is to say, we required the local Jeans length to be resolved with at least ten cells. This strategy results in the highest effective resolution of 4096³ 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: $$A\propto {\widehat{e}}_z\exp \left(-r/R_0\right)\exp \left(-z/H_0\right),$$(8)

where ê_ɀ is the unit vector along the ɀ axis, r represents the cylindrical radius, and ɀ is the vertical height from the galactic plane. The parameters R₀ and H₀ represent the scale length and scale height of the magnetic field, respectively. In this work, the values of R₀ and H₀ 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 created a vector potential in Fourier space using a power spectrum: $$P\left({\tilde{A}}_i\right)\propto k^{-4},$$(9)

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.

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.

3 Results

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

3.1 Morphology and star formation rates of galaxies

Figure 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). 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 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 upon visual inspection could in principle affect the star-forming properties of the galaxies. In Fig. 3 we plot the SFR (left) and the cumulative mass of the new stars (right) as a function of time. The SFR history is generated with a time binning of 10 Myr. We have verified that varying this time interval does not affect the resulting trend of SFR with time. The dashed red line indicates the moment where we switched on star formation in the simulation.

We observe some differences between the two galaxies. In model T, the peak of the SFR occurs at around 250 Myr and reaches approximately 5.5 M_⊙ yr⁻¹; instead, in model R, the peak of the SFR 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⁹ M_⊙ while the total stellar mass from the initial conditions was 1.4 × 10¹⁰ M_⊙. Furthermore, both models exhibit a comparable average SFR, 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 MW-like 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⁻¹kpc⁻² for model T and model R at 250 and 380 Myr, respectively. The corresponding mean surface density of the SFR is 8 × 10⁻⁴ and 7.3 × 10⁻⁴ M_⊙yr⁻¹ kpc⁻². 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 with a high SFR.

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

Fig. 3 SFR (left) and cumulative mass of the new stars ($M_{\star }^{\text{new }}$; right) as a function of time of model T (blue) and model R (green). The dashed red line represents the time when star formation is allowed in the simulation.

3.2 Thermodynamical and magnetic state of the gas

In order to study the dynamical state of the gas in different phases, we first examined 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 took only the gas cells with |z| ≤ 1.5 kpc and radius R ≤ 13 kpc encompassing the entire disk into account¹. We studied the behavior of the gas by categorizing it into three distinct phases defined by number density: low-density (10⁻² < n/cm⁻³ < 1), intermediate-density (1 < n/cm⁻³ < 10²), and high-density phase (10² < n/cm⁻³ < 10⁴).

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 high-density 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⁻² to 10 cm⁻³. We 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))².

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⁻² cm⁻³). As a result, the two atomic phases created by the thermal instability – the warm neutral medium (−1 < log(n/cm⁻³) < 0) and the cold neutral medium, (1 < log(n/cm⁻³) < 2) − are both supported by magnetic pressure.

This magnetic support yields a visible signature in the gas temperature, which is shown³ as a function of density in the bottom panel of Fig. 5. Model R contains a significant gas fraction at low densities (10⁻¹ < n/cm⁻³ < 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 a 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 phases⁴. 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 erg, 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 in β 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.8, which is an order of magnitude higher than the temporally 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 Myr in model T, while in model R, it forms earlier, at 120 Myr. 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 significance of the magnetic with respect to thermal energy in model R compared to model T.

It is indeed crucial to consider the high standard deviation in the data, as it indicates significant variability within the gas included in each density phase. The temperature 2D histogram shown in the bottom plot of Fig. 5 illustrates that the medium in the same density range exhibits different temperatures, leading to diverse distributions of log(β) that can exhibit multiple peaks. This variability should be taken into account when analyzing and interpreting the results.

In the bottom panel of Fig. 6, we show the total thermal (circles) and magnetic (triangles) energy of the low-density (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 1D radial profiles of thermal and magnetic energies for the three different gas phases at three different times (200, 300, and 500 Myr) in Appendix 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.

Fig. 4 Maps of the star formation surface density (ΣSFR) for model T (left) and model R (right) at the peak of their SFR, which happens at t = 250 and 380 Myr, respectively.

Fig. 5 2D mass-weighted histograms of the logarithm of the plasma beta (β; top), total pressure (P; middle), and temperature (T; bottom) of the gas plotted against its density (n) for model T (left) and model R (right) at t = 500 Myr. The dots represent the mean values per density bin. The dashed white line in the top panel indicates the equipartition of thermal and magnetic pressure.

Fig. 6 Evolution of 〈log β〉 (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) and model R (right). In the top plots, the dashed black 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, and triangles correspond to the magnetic energy. The dashed red line represents the onset of star formation. The shaded regions correspond to the standard deviation.

3.3 B−ρ relations

Figure 7 shows the magnetic field-number density (B-n) relation at t = 470 Myr. As in Figs. 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 mass-weighted 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: $$\log (\text{B}/\mu \text{G})=\{\begin{array}{ll}{\kappa }_1\log \left(\text{n}/{\text{cm}}^{-3}\right)\hfill & n<n_0\hfill \\ {\kappa }_2\log \left(\text{n}/{\text{cm}}^{-3}\right)\hfill & n>n_0\hfill \end{array},$$(10)

where the two power-law slopes κ₁ and κ₂ are parameters of the fit, while the transition number density, n₀, is fixed in this instance to 100 cm⁻³. In this figure we show only one snapshot to showcase the break, but later in this section we examine the impact of varying the parameter n₀ 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₀ (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 also present results for a single power-law fit with an index simply referred to as κ.

The slopes κ₁ and κ₂ 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 Sect. 1). We recall here that B ∝ n⁰ results from compression along magnetic field lines, B ∝ n¹ from compression perpendicular to the magnetic field lines, $B\propto n^{\frac{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\propto n^{\frac{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³ (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 κ₁ = 0.31 ± 0.01. For both models, in the high-density medium, the power law steepens with respect to the intermediate, with κ₂ = 0.61 ± 0.02 for model T and κ₂ = 0.53 ± 0.01 for model R. According to Mestel (1965), Mouschovias (1976a,b), and 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 (κ₁,κ₂) 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 considered three different values: n₀ = 50 cm⁻³ (top panel), n₀ = 100 cm⁻³, (middle panel), and n₀ = 300 cm⁻³ (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₀ = 300 cm⁻³. However, we investigated 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₀ = 50 cm⁻³, the intermediate-density slope κ₁ 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, κ₁ 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₀ = 100 cm⁻³ and n₀ = 300 cm⁻³), the behavior of κ₁ is similar to that of n₀ = 50 cm⁻³, but the variations are smoother.

For model R, and a transition density of n₀ = 50 cm⁻³, the high-density slope κ₂ fluctuates around 1/2, implying a slab or filament-like geometry with a perpendicular or radial compression. At the same transition density, κ₂ 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₀ = 50 cm⁻³, we find 〈κ₁〉_T = 0.34 ± 0.19 for model T and 〈κ₁〉r = 0.38 ± 0.09 for model R. Model T also exhibits 〈κ₂〉_T = 0.44 ± 0.07, and model R records an average slope of 〈κ₂〉_R = 0.48 ± 0.08. The errors are attributed to the time evolution, which is the primary source of variance. For higher n0 values, κ₂ exhibits significant variations and differences between the two models.

For a transition density of n₀ = 100 cm⁻³, κ₂ 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, κ₂ 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, 〈κ₁〉_T = 0.36 ± 0.16 for model T and 〈κ₁〉_r = 0.40 ± 0.07 for model R. For higher densities (> n₀), model T exhibits a slope of 〈κ₂〉_t = 0.46 ± 0.15, and model R 〈κ₂〉_r = 0.51 ± 0.10.

In the case of n₀ = 300 cm⁻³, the initial values of κ₂ 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, κ₂ in model T shifts to unity, the theoretical expectation for compression across the field lines. At later times, κ₂ fluctuates strongly between 0.5 and 1. Also, for R, κ₂ 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.

Calculating again the temporal average, the slopes for intermediate densities are 〈κ₁〉_T = 0.38 ± 0.12 for model T and 〈κ₁〉_r = 0.42 ± 0.06 for model R. In the case of high densities, model T shows a slope of 〈κ₁〉_T = 0.53 ± 0.39, while model R 〈κ₁〉_R = 0.60 ± 0.19.

Overall, the B − ρ relation exhibits strong variations across different density ranges and over time. However, despite these variations, the average intermediate slope (〈κ₁〉) consistently appears lower than the average high-density slope (〈κ₂〉). Additionally, for higher densities, both slopes tend to converge toward 1/2 for lower transition densities (n₀), 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 κ₁. 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 κ₁ in the previous cases, when averaged over time, the κ values in the two models are within one sigma from each other. Finally, it is worth noting that the χ² for the single power-law fit is lower than that of the double power law, for both models. We calculated the Bayesian information criterion (BIC; Schwarz 1978) for the two models as an indication for the preferred fit. The BIC – $$\text{BIC}=k\ln (n)-2\ln (L),$$(11)

where κ 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 χ²). 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.

Fig. 7 |B| versus n for the two models at t = 470 Myr with model T (left) and the model R (right). 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), and Tritsis et al. (2015) are reported with dashed black lines.

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) and high densities (right). 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 dashed black and red lines show the slopes of 1/2 and 2/3, respectively.

4 Summary and discussion

4.1 Summary

In this work, we have 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 SFRs of the two models peak at different times, and there are some differences in the instantaneous SFR at different instants. However, the time-averaged 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 Myr 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.

4.2 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⁻³. 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⁻³) exhibited a power-law scaling with a slope of κ = 2/3, while higher densities (n > 0.6 cm⁻³) showed a flatter slope of κ = 1/4. Our results do not align with these studies since we observe scaling for densities higher than 10 cm⁻³, 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) and 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⁻²−10⁴ cm⁻³, which varied depending on the turbulence driving and temperature. They find that the slope of the B − ρ relation is shallower in the solenoidal (k = 0.22 for T < 10³ K, κ = 0.27 for T > 10³ K) than in the compressive case (k = 0.71 for T < 10³ K, κ = 0.51 for T > 10³ K ). In our study, we focus on fitting the B − ρ relation for temperatures lower than 10⁴ 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 〈SFR〉 ≃ 1 M_⊙ yr⁻¹. With 20% of this mass exploding as supernovae, E_SN = 10⁵¹ erg the energy per supernova and 25% of E_SN going to kinetic energy injected in the ISM, we estimate the rate of turbulent kinetic energy from supernovae to be approximately E_turb,SN ≃ 10⁴⁹ erg yr⁻¹ . For the energy contribution from the differential rotation, we calculate the rotational velocity for each cell and use the relation $K_{\text{rot}}=1/2{\text{ m v}}_{\text{rot}}^2$ to calculate the kinetic energy. We then divide this energy by the timestep, which is 10³ yr, to obtain an estimate of E_turb,rot ≃ 10⁵⁴ 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.

4.3 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 coevolution 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 dependences 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.

4.4 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₂. This leads to a more realistic representation of star formation because it is based on the H₂ fraction. However, it is important to note that our simulations currently lack pre-supernova feedback. Including early feedback, namely radiation and stellar winds alter their impact (Pallottini et al. 2017, 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.

5 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 an ordered or a random magnetic field) have a similar average SFR, 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⁻² < n/cm⁻³ < 1) with 〈log β〉 > 3, but model R shows slightly lower 〈log β〉 values. In the intermediate-density gas (1 < n/cm⁻³ < 10²), model R tends toward equipartition, while model T remains thermally dominant with 〈logβ 〉≈ 1. For high-density gas (10² < n/cm⁻³ < 10⁴), 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 for 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⁻³ 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⁻³ < 10³) better than a broken power law, regardless of the choice of the transition density, n₀.

• 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, however, is 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 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 gas-poor, 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.

Acknowledgements

We would like to thank A. Ferrara, P. Hennebelle, R. Teyssier, T. Ch. Mouschovias, R. Skalidis, N. Loudas, G. Korkidis, and K. Kovlakas for insightful discussions related to this project. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282 and from the Hellenic Foundation for Research and Innovation (HFRI), second call for post-doctoral researchers, project number 224. This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility - ARIS - under project ID pr013007_thin. Part of the work has been performed under the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme. The authors gratefully acknowledge the support of Scuola Normale Superiore di Pisa (SNS) and the computer resources and technical support provided by the CINECA supercomputing center. We also gratefully acknowledge the computational resources of the Center for High Performance Computing (CHPC) at SNS and at the University of Crete. We acknowledge usage of the Python programming language (Van Rossum & de Boer 1991; Van Rossum & Drake 2009), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and SciPy (Virtanen et al. 2020).

Appendix A Evolution of thermal and magnetic phase diagrams

Figure A.1 illustrates 2D mass-weighted histograms presenting the logarithm of the gas temperature plotted against the logarithm of gas density. The left side of the figure represents model T, while the right side corresponds to model R. The histograms exhibit snapshots at different time points: t = 200, 300, and 500 Myr. Each histogram features colored dots indicating the mean gas temperature per density bin. Over time, model R develops a significant gas fraction at low densities (10⁻¹ < n/ cm⁻³ < 10 cm⁻³) characterized by temperatures below 100K. In contrast, this specific temperature-density regime remains unoccupied in the corresponding plot for model T.

Figure A.2 displays 2D mass-weighted histograms showing the logarithm of the total gas pressure plotted against the logarithm of gas density. The left side of the figure corresponds to model T, while the right side corresponds to model R. The histograms are shown at different time snapshots: t = 200, 320, 450, and 500 Myr. In each histogram, the blue dots represent the mean values of the thermal pressure per density bin, while the green dots represent the mean values of the magnetic pressure per density bin. The initial conditions are indicated by gray stars.

At 200 Myr, both models are thermally dominated for densities lower than 10 cm⁻³ and magnetically dominated for higher densities. As time progresses, model T transitions to equipartition, with magnetic and thermal contributions becoming comparable in higher densities, while maintaining thermal dominance in lower densities. In contrast, model R continues to become magnetically dominated even in lower densities over time.

Fig. A.1 2D mass-weighted histograms of the logarithm of the temperature of the gas plotted against log(n) for model T (left) and model R (right) at t = 200, 300, and 500 Myr. The colored dots represent the mean temperature per density bin.

Fig. A.2 2D mass-weighted histograms of the logarithm of the total pressure of the gas plotted against log(n) for model T (left) and model R (right) at t = 200, 320, 450, and 500 Myr. The blue dots represent the mean thermal pressure per density bin and the green ones the magnetic pressure. The gray stars represent the initial condition.

Appendix B Radial profiles of thermal and magnetic energies for the three different phases of the gas

We plot 1D radial profiles of thermal and magnetic energies for the three different phases of the gas. The initial condition for both models comprises exclusively low-density gas, where thermal energy dominates over magnetic energy across the entire range of radii. However, in the central region of the galaxy, the magnetic energy approaches the thermal energy, and it rapidly decreases as we move to larger radii compared to the thermal energy.

At 200 Myrs, the low-density gas is thermally dominated across the entire range of radii for model T. In contrast, for model R, within 1 kpc of the center, the mean magnetic energy is slightly higher than the thermal energy. The medium-density gas is magnetically dominated for the first 4 kpc from the center in both models. Then, it becomes thermally dominated. On the other hand, high-density gas extends up to 4.5 kpc from the center in model T, where it remains magnetically dominated. In model R, it primarily extends to 8 kpc, with occasional occurrences as far as 12.5 kpc. In model R, the high-density gas is magnetically dominated until 5 kpc and then transits to thermal dominance.

At 300 Myrs, the low-density gas is magnetically dominated for the first 1.5 kpc in both models, while beyond this region, it becomes thermally dominated. It is worth noting that the radius of model R extends to approximately 11 kpc, whereas model T reaches around 13.5 kpc. The medium-density gas remains magnetically dominated for the initial 4 kpc from the center in both models, after which it transits to thermal dominance, consistent with its state at 200 Myrs. The high-density gas extends to larger radii compared to the 200 Myrs snapshot in model T. In both models, they maintain magnetic dominance up to 4.5 kpc, after which they transit to thermal dominance.

At 500 Myrs, there are notable differences in the behavior of low-density gas between the two models. In model T, there is a transition from equipartition in the central region to slight magnetic dominance extending up to 3 kpc. Beyond this point, thermal energy prevails. Conversely, in model R, the gas exhibits strong magnetic dominance from the center up to 3 kpc, with the magnetic energy being approximately two orders of magnitude higher than thermal energy. At 4 kpc, it reaches equipartition before transitioning to thermal dominance. The medium-density gas maintains magnetic dominance for the initial 4 kpc from the center in model T and up to 5 kpc in model R, with model R having higher magnetic energy compared to model T in this region. Subsequently, both models transition to thermal dominance. The high-density gas is magnetically dominated up to 4 kpc in model T and up to 6 kpc in model R.

Overall, the low-density gas is thermally dominated for the majority of the disk, but the central regions appear to be magnetically dominated, particularly in the case of model R, where the extent of magnetic dominance increases over time. In the medium-density gas, both models consistently maintain magnetic dominance for the initial few kiloparsecs from the center before transitioning to thermal dominance. Notably, model R often exhibits higher magnetic energy compared to model T in this region. For the high-density gas in both models, magnetic dominance is the initial state, followed by a transition to thermal dominance at certain radii. In summary, these findings, in line with the patterns observed in the previous plots, suggest that model R consistently demonstrates a higher prevalence of magnetic dominance across various density regimes and time snapshots compared to model T.

Fig. B.1 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 200 Myrs.

Fig. B.2 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 300 Myrs.

Fig. B.3 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 500 Myrs.

References

All Figures

	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.
In the text

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

	Fig. 3 SFR (left) and cumulative mass of the new stars ($M_{\star }^{\text{new }}$; right) as a function of time of model T (blue) and model R (green). The dashed red line represents the time when star formation is allowed in the simulation.
In the text

	Fig. 4 Maps of the star formation surface density (ΣSFR) for model T (left) and model R (right) at the peak of their SFR, which happens at t = 250 and 380 Myr, respectively.
In the text

	Fig. 5 2D mass-weighted histograms of the logarithm of the plasma beta (β; top), total pressure (P; middle), and temperature (T; bottom) of the gas plotted against its density (n) for model T (left) and model R (right) at t = 500 Myr. The dots represent the mean values per density bin. The dashed white line in the top panel indicates the equipartition of thermal and magnetic pressure.
In the text

Fig. 6 Evolution of 〈log β〉 (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) and model R (right). In the top plots, the dashed black 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, and triangles correspond to the magnetic energy. The dashed red line represents the onset of star formation. The shaded regions correspond to the standard deviation.

In the text

Fig. 7 |B| versus n for the two models at t = 470 Myr with model T (left) and the model R (right). 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), and Tritsis et al. (2015) are reported with dashed black lines.

In the text

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) and high densities (right). 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 dashed black and red lines show the slopes of 1/2 and 2/3, respectively.

In the text

	Fig. A.1 2D mass-weighted histograms of the logarithm of the temperature of the gas plotted against log(n) for model T (left) and model R (right) at t = 200, 300, and 500 Myr. The colored dots represent the mean temperature per density bin.
In the text

	Fig. A.2 2D mass-weighted histograms of the logarithm of the total pressure of the gas plotted against log(n) for model T (left) and model R (right) at t = 200, 320, 450, and 500 Myr. The blue dots represent the mean thermal pressure per density bin and the green ones the magnetic pressure. The gray stars represent the initial condition.
In the text

	Fig. B.1 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 200 Myrs.
In the text

	Fig. B.2 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 300 Myrs.
In the text

	Fig. B.3 Radial profiles of the mean of the logarithm of thermal and magnetic energies for low, medium, and high densities at 500 Myrs.
In the text