The global structure of magnetic fields and gas in simulated Milky Way-analogue galaxies

ABSTRACT

We simulate an isolated, magnetized Milky Way-like disc galaxy using a self-consistent model of unresolved star formation and feedback, evolving the system until it reaches statistical steady state. We show that the quasi-steady-state structure is distinctly layered in galactocentric height z, with a broken power-law structure in Alfven Mach number and plasma beta. Magnetic pressure exceeds turbulent and thermal pressures after the gas is depleted to levels below that of the present-day Galaxy, but is subdominant at higher gas fractions and star formation rates. We find field strengths, gas surface densities, and star formation rates that agree well with those observed in the Solar neighbourhood. The most significant dynamical effect of magnetic fields on the global properties of the disc is a reduction of the star formation rate by a factor of 1.5–2 with respect to an unmagnetized control simulation. At a fixed star formation rate of approximately |$2 \, {\rm M}_{\odot }$| yr⁻¹, there is no significant difference in the mass outflow rates or profiles between the magnetized and non-magnetized simulations. Our results for the global structure of the magnetic field have significant implications for models of cosmic ray-driven winds and cosmic ray propagation in the Galaxy, and can be tested against observations with the forthcoming Square Kilometre Array and other facilities. Finally, we report the discovery of a physical error in the implementation of neutral gas heating and cooling in the popular gizmo code, which may lead to qualitatively incorrect phase structures if not corrected.

1 INTRODUCTION

The role of magnetic fields in the Galaxy has been long debated. The discovery of polarized starlight by Hiltner (1949) and Hall & Mikesell (1949) led to the hypothesis that the observed polarization was caused by dust grains aligned with the magnetic field in the interstellar medium, with the immediate implication from the observed polarization direction that Galactic magnetic fields are aligned parallel to the Galactic plane, i.e. are toroidal (Davis & Greenstein 1949, 1951). At nearly the same time, Fermi (1949) theorized that fluctuating magnetic fields in the interstellar medium were the origin of cosmic rays (via ‘second-order Fermi acceleration’). Soon afterward, the Galactic radio emission first observed by Jansky and Reber was proposed to be due to the radiation produced by the gyration of cosmic rays around magnetic fields (Kiepenheuer 1950). These theoretical advances were followed by radio observations of polarized Galactic emission and Faraday rotation in the interstellar medium in the 1950s and 1960s (see Wielebinski 2012 for a detailed historical review of polarized radio observations).

In parallel with this observational progress, theorists began to consider the possible dynamical role of magnetic fields. This topic was first considered by Alfvén (1942), who argued that magnetic fields may be dynamically important on the surface of the sun through their wave interactions with partially ionized gas, and Fermi (1949), who recognized that the same considerations applied to the interstellar medium, further proposed that the energy of the Galactic magnetic field should be of order the energy of the turbulent motions of the gas, implying a magnetic field strength of order ∼|$\mu$|G. However, the topology of the field, and in particular its structure as one moves away from the Galactic plane, remain uncertain. One possible picture is provided by Parker (1966), who discovered an MHD instability of toroidal fields that leads to the magnetic field buckling into loops above the plane. However, Parker’s instability was suppressed in the magnetohydrostatic model of Boulares & Cox (1990), which proposed a significant vertical (i.e. poloidal) component of the magnetic field at ∼kpc heights above the Galactic plane and allowed for the diffusion of cosmic rays within a region of tangled (but mean-toroidal) magnetic fields near the plane, yielding a solution with approximate energy equipartition between turbulent motions, cosmic rays, and magnetic fields, and an extended distribution of H ii gas at heights |z| ≳ 1 kpc (the so-called Reynolds layer; Reynolds 1989). This basic picture of the magnetohydrostatic steady-state structure of the gas, cosmic rays, and magnetic fields of the Galaxy has survived to the present.

The uncertainty in the field topology is only one of several outstanding questions about the origin and dynamical role of the magnetic field. One major area of uncertainty is the origin of the field, and whether it is governed by a mean field dynamo process (e.g. Gent et al. 2021). A second is how magnetic fields interact with galactic winds driven by either supernova breakout (e.g. Tomisaka 1998) or by cosmic rays (e.g. Breitschwerdt, McKenzie & Voelk 1991; Everett et al. 2008; Mao & Ostriker 2018; Quataert, Thompson & Jiang 2022). There have been some attempts to address these questions in magnetohydrodynamic (MHD) cosmological simulations (e.g. Pakmor & Springel 2013; Pakmor et al. 2017; Su et al. 2017; Pakmor et al. 2018; Hopkins et al. 2020; we note that Pakmor & Springel 2013 carried out isolated simulations but with comparable resolution and physics to cosmological simulations) and isolated disc MHD simulations (e.g. Ntormousi 2018; Körtgen et al. 2019), but these efforts have been somewhat limited by resolution. Non-zoom-in cosmological simulations, such as the highest resolution Illustris-TNG simulations (Pillepich et al. 2019; Nelson et al. 2019), have a resolution of Δx ∼ 150 pc at the mean density of the Milky Way’s interstellar medium (n_H ∼ 1 cm⁻³), meaning the gas scale height of the Milky Way (h ∼ 100 pc) is resolved with less than one gas particle in the diffuse interstellar medium. Zoom-in cosmological simulations do only slightly better, reaching Δx ∼ 60 pc at n_H ∼ 1 cm⁻³, thereby resolving the scale height with ≲ 2 gas particles (e.g. Pakmor et al. 2017; Hopkins et al. 2020). This is clearly insufficient to capture the topology of the field and its changes in and out of the plane.

On the other end of the resolution spectrum lie MHD simulations of local patches of the Milky Way, such as those of Kim et al. (2016), Kim, Choi & Flauger (2019), Kim et al. (2020), Kim & Ostriker (2017), Kim & Ostriker (2018), and Rathjen et al. (2021). These simulations benefit from uniform parsec-scale resolution of all of the phases of the interstellar medium, but they cannot resolve the global structure of the disc in radius and height, which means that cannot address questions of the field topology. Additionally, as emphasized by Martizzi et al. (2016), any local simulation does not allow streamlines to diverge in outflows that would otherwise be spherical, thus preventing the outflow from crossing the sonic point of a classical hot superwind (Chevalier & Clegg 1985). This limits their ability to study the interaction of magnetic fields with outflows.

These limitations in previous work motivate us to consider an intermediate resolution regime. We present a new dynamical model of the Galaxy, evolved for ∼1 Gyr until it has reached a quasi-steady-state. We reach roughly an order of magnitude higher mass resolution than even zoom-in cosmological simulations of Milky Way-like galaxies, and, though our resolution is still substantially smaller than that in local patch simulations, we retain the full geometry of the problem, so that we can study field topology and outflows. Relatively few simulations of this type have been published, and those have largely been concerned with studying the magnetization of the neutral medium (e.g. Wang & Abel 2009) or attempting to explain and interpret the observed Faraday rotation sky (e.g. Kulpa-Dybeł et al. 2015; Butsky et al. 2017). There have been no previous efforts to use simulations of this type to map out the vertical structure and topology of the magnetic field, or to study how fields interact with galactic winds. These questions are the focus of our study.

This paper is organized as follows: Section 2 details the numerical methods used in our simulations, including the ‘subgrid models’ used for star formation and feedback; Section 3 describes the initial conditions and evolution of our simulations as they relax into a quasi-steady-state; Section 4 summarizes the zonal structure of the simulations in steady state, the dependence of outflow rate on star formation rate (SFR), and the observations to which our simulations may be compared. We conclude in Section 5 with a summary of our results and possible future directions for research.

2 METHODS

We solve the equations of ideal MHD using the gizmo code (Hopkins 2015), which implements the method of Hopkins (2016) and Hopkins & Raives (2016). This method significantly improves upon the divergence-cleaning approach of Dedner et al. (2002) by the addition of a local approximate Hodge projection that further suppresses numerical magnetic monopoles in the discrete magnetic field. Note that we use the MHD solver even when running simulations with zero magnetic field in order to ensure that our results are not biased by the use of a different solver in different simulations. Since we initialize the magnetic field in these simulations to zero, it remains exactly zero at all subsequent times.

In addition to MHD, we solve a time-dependent chemistry network for the abundance of H i, H ii, He i, He ii, He iii, and free electrons that we use to compute the atomic cooling rates for hydrogen and helium. Additionally, we interpolate from a table of metal line cooling rates as a function of density and temperature (assuming ionization equilibrium and solar metallicity) that was computed using cloudy (Ferland et al. 1998) following the method of Smith, Sigurdsson & Abel (2008), as implemented in the Grackle chemistry and cooling library (Smith et al. 2017).¹ We assume an optically thin, spatially uniform photoionizing background radiation field based on the redshift z = 0 tabulation from Haardt & Madau (2012) when computing the ionization state and cooling rates for both primordial species and metals. For gas at temperatures T < 2 × 10⁴ K, we also assume a photoelectric (volumetric) heating rate

where n_H is the sum of H i and H ii number densities, which is the default setting for photoelectric heating in Grackle and agrees at order of magnitude with the solar neighbourhood value of photoelectric heating calculated by Wolfire et al. (2003, their equation 19).

The Grackle cooling implementation results in a three-phase neutral ISM, with a warm phase (WNM), a cool phase (CNM), and an unstable phase at intermediate temperatures. We use Grackle rather than the default cooling implementation in gizmo, as used in, e.g. the FIRE-2 simulations (Hopkins et al. 2018b), because the default gizmo cooling implementation contains an error that prevents it from producing an unstable phase of the neutral ISM, as shown in Appendix  A.

Using the implementation in gizmo (originally based on that of Springel 2005), we form stars by stochastically converting gas particles into star particles (e.g. Katz, Weinberg & Hernquist 1996) such that the expectation value of the instantaneous SFR density satisfies

where p is the probability of star formation, ρ is the gas density, ϵ_⋆ is the star formation efficiency parameter, and

is the gas free-fall time-scale. We choose the critical gas density |$\rho _{\text{crit}} = 100 \, \text{H cm}^{-3}$|⁠, where this density is approximately the Jeans density for 50 K gas at our gas particle mass resolution, and we set the star formation efficiency parameter ϵ_⋆ = 0.01 to match the observed star formation efficiency per free-fall time in dense molecular clouds (see e.g. Fig. 10 of Krumholz, McKee & Bland-Hawthorn 2019 and references therein). For densities ≥100 times that of the critical density ρ_crit, we increase the star formation efficiency parameter to unity in order to avoid runaway collapse in Jeans-mass-unresolved dense regions with infinitesimal time-steps.² We note that this star formation criterion implies that over a time-step, the star formation probability p for a given gas particle takes the form

since we must integrate equation (2), which gives a probability rate, over a time-step Δt in order to obtain a probability p (Katz et al. 1996).

We implement photoionization feedback following the method of Armillotta et al. (2019), which re-implements the Stromgren volume method described by Hopkins et al. (2018b). However, because our resolution is lower than that of Armillotta et al. (2019), for this paper we do not use stochastic sampling from the IMF for each star particle. Instead, for simplicity, we assume that the stellar initial mass function is fully sampled, and we adopt a constant ionizing luminosity of 10³⁹ ionizing photons per second per 100 M_⊙ of stellar mass from birth to 5 Myr after formation, and zero thereafter. We note that in our implementation, while a gas particle is identified as being within an H ii region, the cooling and heating source terms are turned off for that particle.³

For supernova feedback, we use the method of Hopkins et al. (2018a) as implemented in gizmo. In simplified form, this method couples the momentum from the unresolved Sedov–Taylor phase to the faces of the neighbouring fluid elements. Then, after subtracting the resulting kinetic energy from a fiducial explosion energy of 10⁵¹ erg, the remaining energy is coupled to the neighbouring gas particles as a thermal energy source term, with a minimum thermal energy injection of one-half of the initial explosion energy. This method is very similar to the algorithm of Kimm & Cen (2014), except that the gizmo algorithm also adds the momentum of the pre-shock ejecta and has minor differences in the treatment of the difference between the simulation frame and the explosion frame. The fiducial normalization of the maximum injected momentum from the unresolved Sedov–Taylor phase used by gizmo is ⁴

where

and

This normalization is somewhat larger than the commonly used normalization of 3 × 10⁵ M_⊙ km s⁻¹ (Thornton et al. 1998), but is consistent with the wide range of values proposed in the literature (e.g. Kim & Ostriker 2015; Gentry et al. 2017, 2019; Gentry, Madau & Krumholz 2020). Following the default used in gizmo, we assume a constant supernova rate of 3 × 10⁻⁴ SNe M_⊙⁻¹ Myr⁻¹ for star particles with ages 0 < t_age < 30 Myr. Each event is assumed to have an explosion energy of 10⁵¹ erg. We neglect type Ia supernova explosions in our model.

3 SIMULATIONS

3.1 Initial conditions

We carry out two simulations: one with hydrodynamics only and one with MHD. The initial conditions of the gas and collisionless components in both simulations are identical to those used by the AGORA isolated disc galaxy comparison project in their ‘high-resolution’ case (Kim et al. 2016). These include a dark matter halo component of mass M₂₀₀ = 1.07 × 10¹² M_⊙ (defined as the mass enclosed within a mean density of 200 times the critical density) with a concentration parameter c = 10, a stellar disc of mass 3.4 × 10¹⁰ M_⊙, a bulge of mass 4.3 × 10⁹ M_⊙, and a gas disc of mass 8.6 × 10⁹ M_⊙. This implies a gas fraction of ∼20 per cent by mass in the initial conditions. The stellar disc has a scale height of approximately 350 pc and thus represents the observed stellar ‘thin disc’ in the Galaxy. The gas disc has a scale height of R₀ = 3.43218 kpc and scale length z₀ = 0.343218 kpc, with the gas disc initialized with the density profile

The gas temperature and circular velocity are initially computed via solution of the Jeans equations to be in hydrostatic and centrifugal equilibrium using the method described by Springel, Di Matteo & Hernquist (2005). However, we override the hydrostatic temperature values with a uniform gas temperature of 10⁴ K within the disc in order to allow the disc to rapidly cool and collapse to form stars. After stars form, the disc is then partially re-inflated by the injection of momentum and energy from supernovae. Our simulations use gas particles with mass 859.3 M_⊙, dark matter particles with mass 1.254 × 10⁵ M_⊙, and stellar disc and bulge particles with mass 3.4373 × 10³ M_⊙. (Star particles formed during the simulation have the mass of the gas particle from which they were stochastically converted.)

In the MHD simulation, the initial magnetic field is

where B₀ = 10 |$\mu$|G. This field is analytically divergence free. However, when discretized on to the gas particles in the initial conditions it is not, but the residual numerical divergence is rapidly transported outside the domain via divergence-cleaning and local approximate Hodge projection (Hopkins 2016). This field geometry (purely toroidal) and strength are chosen to be in rough approximate agreement with the observed ‘ordered’ (or ‘regular’) magnetic field of the Galaxy (Beck 2015), with the expectation that the turbulent component of the field (as well as any non-toroidal component) would be generated by the gravitational collapse and stellar feedback in the simulation.

3.2 MHD simulation

We run the MHD simulation for a time t = 976 Myr. The face-on projected density, the face-on projected SFR, and the face-on, in-plane magnetic field of the final output of the simulation are shown in the left column of Fig. 1. We see morphology typical of an S0 galaxy, with relatively weak spiral arms and no visible bar. The lack of a bar or a grand design spiral pattern may be due to a lack of close-in orbiting satellites, such as the Large Magellanic Cloud, which are absent in our simulation. The surface density normalization is roughly consistent with observations of the Milky Way, with values in the range of 50−100 M_⊙ pc⁻².

The magnetic field strength generally ranges from 10 to 100 |$\mu$|G, which is also consistent with observations when compared with the total magnetic field strength, not just the so-called ordered field. We have illustrated the magnetic field with arrows to convey the sense of the magnetic field direction. There do not appear to be any magnetic field reversals in the azimuthal direction (i.e. along a circular orbit), which is in tension with the common interpretation of the Galactic Faraday rotation measurements that suggest a reversal of the azimuthal magnetic field direction between spiral arms (Beck 2015). However, since our initial conditions do not have any magnetic field reversals (the initial field is everywhere aligned in the +ϕ direction), nor any gas accretion from halo gas or cosmological sources, the lack of field reversals may not be unexpected. The mechanism for generating field reversals is unknown, although they are seen in some cosmological simulations of magnetic fields (Pakmor et al. 2018), and are suggested to form as a result of a mean-field dynamo process (Beck 2015).

In the top left panel of Fig. 2, we show the SFR as a function of elapsed simulated time. We see that it initially spikes at a value of around 10 M_⊙ yr⁻¹, as a the disc cools and there is insignificant turbulence (or any other feedback processes) to slow the collapse of gas and subsequent formation of stars. The SFR then slowly declines as collapse and feedback processes come into a quasi-equilibrium toward the end of the simulation, plateauing at a rate of ∼1−2 M_⊙ yr⁻¹. As star formation proceeds, the gas is necessarily consumed. In the top right-hand panel of Fig. 2, we show the gas fraction of the disc (computed as the mass of gas particles relative to the total non-dark-matter mass). From the initial gas fraction of ∼0.2, we see a slow decline over ∼1 Gyr to a value of ∼0.1 gas fraction. The decline in gas fraction is somewhat more rapid than the decline in SFR, so the total depletion time (bottom left-hand panel of Fig. 2), defined as the ratio of the gas mass to the SFR, very gradually increases as the simulation runs, reaching ≈6 − 8 Gyr at the final time. This is roughly consistent with the depletion times observed in nearby spiral galaxies, where molecular gas depletion times average ≈2 Gyr (e.g. Leroy et al. 2013), but molecular gas constitutes only ≈1/3−1/4 of the total gas mass (with the balance as H i; Saintonge et al. 2011), so the total gas depletion time is a ≈6−8 Gyr.

Finally, the lower right-hand panel of Fig. 2 shows the gas mass-weighted mean magnetic field of the MHD simulation. The initialization of the magnetic field according to equation (11) implies a mass-weighted mean field of ∼1 |$\mu$|G. The field is subsequently strongly amplified by the initial burst of star formation described previously, and then relaxes into a quasi-steady-state value of ∼7 |$\mu$|G.

In Fig. 3, we show the mass-weighted distribution of temperature and gas density in the MHD simulation (left-hand panel) and the mass-weighted distribution of thermal pressure and gas density in the MHD simulation (right-hand panel). For comparison, we also overplot as solid lines the equilibrium temperature and thermal pressure as a function of density for our adopted radiative heating and cooling physics (using Grackle version 2.2; Smith et al. 2017). We see that the gas in the simulation quite closely tracks the equilibrium curves, except at low densities, where the gas is far out of thermal equilibrium due to shock heating from supernova feedback, and in dense clouds that are heated to ∼10⁴ K due to photoionization from stars ≲ 5 Myr old, in good agreement with similar simulations of the Milky Way ISM (e.g. Goldbaum, Krumholz & Forbes 2016; Kim & Ostriker 2017). There is also a small systematic offset from thermal equilibrium in the warm neutral medium due to the ionization equilibrium time-scale generally exceeding the cooling time in this phase, an effect anticipated from theory (Wolfire et al. 2003). Some recent simulations (e.g. Hopkins et al. 2018b; Gurvich et al. 2020) show large deviations of cold neutral gas from its thermal equilibrium state. Such behaviour is almost certainly incorrect given the extremely short cooling times of cold neutral gas (e.g. Wolfire et al. 2003) and is more likely due to the error affecting the heating and cooling rates of dense gas in these simulations that we identify in Appendix  A.

3.3 Hydrodynamics-only simulation

We run an additional simulation without magnetic fields as a control in order to examine the differences between the magnetic and non-magnetic simulations. For this simulation, we simply set the initial magnetic field to zero everywhere, and leave all other properties of the initial conditions as they are in the previous simulation. We use the same code settings, including using the MHD solver, in order to ensure that the differences between the two simulations are entirely due to the presence (or absence) of magnetic fields, not the numerical properties of the hydrodynamic versus MHD solver.

Examining the bulk properties shown in Fig. 2, we see that the non-magnetized simulation has a greater initial burst of star formation, peaking just above 12 M_⊙yr⁻¹ (upper left-hand panel). This burst is reflected in a steeper initial decline of the gas fraction (upper right-hand panel). In contrast to the magnetized simulation, following the initial burst, the SFR rapidly declines, undergoes a second burst at a time ∼200 Myr, and then slowly declines over the next several hundreds of Myr to an SFR of ∼2 M_⊙yr⁻¹. The larger initial burst but slower decline means that the gas fraction of the hydro simulation is very similar to the MHD simulation after both simulations are stopped after t ∼ 1 Gyr (upper right-hand panel). As a result of the slightly higher quasi-steady-state SFR compared to the magnetized simulation, the gas depletion time-scale has a quasi-steady-state value of ∼4 Gyr when the simulation is stopped, which is 1.5–2 times lower than the value for the magnetized disc.

4 DISCUSSION

4.1 Global magnetic structure

One of the primary goals of our study is to determine the global structure of the magnetic fields in and around Milky Way-like disc galaxies, and the relationship of the field structure to other physical quantities. To begin this investigation, in Fig. 4 we show the azimuthally averaged (i.e. in the ϕ-direction) structure of the disc in a number of quantities; all quantities shown are mass-weighted means except for the mass flux, which is an intrinsically volumetric quantity. These projections in the radius-vertical height plane illustrate the distinct ‘atmospheric’ components present in the disc. The azimuthal averages for Alfven Mach number and plasma beta clearly show three distinct zones – a very thin disc ∼300 pc in size near the mid-plane, a thicker, ∼2−3 kpc-wide zone around that, and then a distinct third zone at larger heights. This structure is reminiscent of a multilayered cake. Henceforth, we will refer to this stratification as the ‘layer-cake structure’ of the disc.

In the upper left-hand panel of Fig. 4, we see that the gas density is approximately exponentially distributed in the vertical direction, as expected. There is some substructure within the disc, as well as a few kiloparsec-scale plumes of gas above the disc as a result of supernova-driven hot outflows.⁵ The scale height of the gas density increases with galactocentric radius, indicating a ‘flaring’ disc structure. In Appendix  B, we quantify this flaring with profiles of gas scale height versus galactocentric radius at both our intermediate simulation epoch (t ∼ 550 Myr) and our final simulation epoch (t ∼ 1 Gyr). We also show radial column density profiles at the two times for comparison.

In the upper right-hand panel, we see the magnetic field strength, with the direction of the field indicated with arrows. The magnetic field strength in the galactic centre is of order 100 |$\mu$|G, consistent with observations (Beck 2015). The field strength falls off with height and radius less steeply than the gas density, a phenomenon observed in other simulations and suggested by observations. However, unlike the gas density, the magnetic field does not vary smoothly with height. Instead, it is significantly tangled at heights of within ∼1 kpc, consistent with the schematic field geometry of Boulares & Cox (1990) that was suggested by observations at the time. At large heights (>3 kpc), the field becomes coherent on kiloparsec scales, with the projected field lines stretching nearly vertically out of the galaxy near zero galactocentric radius, whereas the projected field lines further away in radius (>5 kpc) curve towards the disc and may form toroidal flux tubes at heights of several kiloparsecs above the disc, as seen at R ≈ 10 kpc and z ≈ 7 kpc in the image.

The middle left-hand panel shows the temperature. We see a cold gas disc (50−100 K) in the inner z < 300 pc region. The cold disc is surrounded by a region of significant spatial extent, ranging from ∼300−500 pc to ∼3 kpc, with gas at temperatures between several thousand Kelvins and ∼10⁴ K. The outer part of this region can be identified with the so-called Reynolds layer (Reynolds 1989) of diffuse ionized gas surrounding the Galaxy. Further out in height (>3 kpc), the gas temperature is typically around 10⁶ K, originating from the hot outflows driven by supernovae, intermixed with ‘plumes’ of rapidly cooling intermediate temperature (∼10⁵ K) gas.

The middle right panel of Fig. 4 shows the vertical mass flux, which we define as |$\dot{M}_z \equiv \rho v_z\mbox{sgn}(z)$|⁠, i.e. positive values indicate mass flow away from the mid-plane, while negative values indicate flow towards it. We use a diverging logarithmic colour scale, so that colour indicates flow direction – outflows are purple and inflows are orange. The strength of outflows/inflows is largest at the mid-plane (z = 0) of the galaxy, with a sharp drop-off in magnitude with height, and the direction of inflow or outflow shows many local reversals. This is suggestive of a fountain-type of outflow. At large heights (>3 kpc), by contrast, we see mostly coherent outflow, with only small patches of inflow, suggesting that this region is not primarily a fountain, but instead represents a true wind of mass leaving the galaxy.

The dimensionless plasma beta (P_thermal/P_magnetic) is plotted in the lower left-hand panel. We see the same basic ‘layer-cake’ structure in plasma beta as a function of galactocentric height, but with a more pronounced intermediate region. At time t ∼ 550 Myr, the gas disc is supported by turbulent pressure, with magnetic pressure about one order of magnitude less significant, and thermal pressure about one order of magnitude smaller still. However, this is not a steady-state situation: as the gas depletes, the ratio of turbulent to magnetic pressure decreases, so by t ∼ 1 Gyr magnetic pressure is a factor of several larger at the mid-plane, as we discuss further below. The bubble-like substructure is likely a result of individual H ii regions and supernova remnants. Above this, there is a ‘corona’ of magnetically dominated gas extending to ∼3 kpc, with stronger magnetic dominance toward the galactic centre. Above this corona is a gas-pressure-dominated region that is the product of hot outflows. All of these features are likely related to the decreased turbulent input power as a result of the lower SFR (making the turbulent pressure significantly lower) and the smaller gas fraction (making the magnetic pressure slightly larger).

Finally, the lower right-hand panel shows the ratio of the toroidal component (B_ϕ) to the total magnetic field strength. Again, we see a zonally stratified structure near the plane. First focus on the outer disc, R ≳ 5 kpc. In this radial region, near the plane at z ≲ 300 pc, the mean toroidal fraction is ≈0.5 (indicated by white on the plot), with a great deal of substructure corresponding to the positions of individual supernova remnants. Above this, at z ∼ 0.3−3 kpc, is a zone where the field is |$\approx 80{{\ \rm per\ cent}}$| toroidal, while at even larger height, z ≳ 3 kpc, the field becomes much less toroidal. The same layering is evident closer to the galactic centre, at R ≲ 5 kpc, except that each of these zones is thinner, so the transitions between them happen at smaller z. At the very centre of the galaxy, R ≲ 100 pc, the field is almost purely poloidal at all heights.

While Fig. 4 provides qualitative evidence for the existence of distinct layers, to demonstrate it quantitatively, and to provide more accurate estimates for the thicknesses of the layers, we compute vertical profiles through the disc; we define the vertical profile q_R(z) for some quantity q at galactocentric radius R by

i.e. our vertical profiles are azimuthal averages computed at fixed radius. We show the vertical profiles of density and magnetic field at t ∼ 560 Myr and ∼1 Gyr in Figs 5 and 6, respectively. It is clear that these profiles are well described by a two-component broken power law, consistent with our visual inspection in Fig. 4 that density and total magnetic field strength do not show three distinct layers. We also see that the magnetic field has a significantly greater vertical extent than the mass. We carry out simple least squares fits to the profiles, which we show as dashed lines in the figures, leaving the break height, the inner and outer slopes, and an overall normalization as free parameters; the best-fitting break heights are printed in the figure legend. We find that the average break height for the density profiles is ≈215 pc at 1 Gyr, whereas the average break radius for the magnetic field strength profiles at the same time is ≈446 pc.

The Alfvén Mach number and plasma beta vertical profiles have a more complicated structure, shown in Figs 7 and 8. (See Appendix  C for details on how we compute the former of these two.) Here, we clearly see the three-zone structure that is easy to pick out by eye in Fig. 4; this structure is present at both ∼560 Myr and ∼1 Gyr, confirming that it is a long-lived, quasi-steady feature. We fit these profiles with a three-component power-law model, with the inner and outer power slopes fixed to zero, so the free parameters are the inner and outer amplitudes and the heights of the two breaks. Since the profiles are constrained to be continuous, these four parameters fully parametrize the three-component model. We again show the best-fitting models as dashed lines, and print the best-fitting heights of the two breaks in the figure legend.

For the Alfven Mach number, we find sub-Alfvenic regions in the inner R = 1.0 kpc up to 10 kpc in height. However, we note that this may not be representative of the Galaxy, since this snapshot (at t ∼ 1 Gyr) has a substantially suppressed SFR compared to the typical Galactic SFRs and is also dominated by magnetic pressure rather than turbulent pressure. Excluding the R = 12.5 kpc profile, since it appears to be an outlier, the average break heights are at ≈440 pc and ≈2100 pc. The plasma beta profile has a similar structure, although with different best-fitting values for the inner and outer break heights. Again excluding the R = 12.5 kpc profile, we find inner and outer break heights of ≈390 pc and ≈3900 pc at t ∼ 1 Gyr. These profile fits motivate our approximate three-zone ‘layer-cake’ decomposition.

Finally, we examine the sources of vertical pressure support that keep the gas disc in approximate vertical dynamical equilibrium. In Fig. 9, we show the average turbulent, magnetic, and thermal pressures as a function of height above the disc at a representative galactocentric radius R = 5.0 kpc. (See Appendix  C for a discussion of how we compute the turbulent pressure.) We see the turbulent pressure dominates at almost all heights for the intermediate snapshot (t ∼ 550 Myr; left-hand panel) whereas magnetic pressure dominates for the final snapshot (t ∼ 1 Gyr; right-hand panel). Thermal pressure is subdominant to both turbulent and magnetic pressure for both simulation epochs, except at the largest heights (≳3 kpc), where thermal pressure becomes more important than magnetic pressure.

4.2 Comparison to observations

It is interesting to compare our vertical profiles to previous observational estimates of the vertical structure of the gas in the galactic disc. A challenge for any quantitative comparison is that previous estimates have implicitly assumed an exponential (or, closely related, sech²) vertical profile, and quoted properties in terms of a scale height, whereas we find more profiles that are better described as power laws. While we can compute heights at which the profile breaks, this is no unambiguous way to define a scale height from our measured profiles. Thus, a full comparison would be best accomplished by forward modelling the observations from our simulations, an effort we defer to future work. None the less, we can compare qualitatively here.

We find that the total density in our simulations is roughly flat to ∼200 pc, and declines as a relatively steep power law thereafter. Qualitatively, our estimate is consistent with the total mass scale height of 160 pc obtained by Boulares & Cox 1990, and the much more recent estimate of 170 pc from McKee, Parravano & Hollenbach (2015). The density of ionized gas at larger heights is also constrained by pulsar dispersion measures, which provide evidence for the Reynolds layer (Reynolds 1989). Reynolds finds that, at the Solar Circle, the free electron density is ≈10⁻² cm⁻³ at z = 1 kpc, falling to ≈10⁻³ at z = 4 kpc. A comparison of these estimates to the simulation profiles shown in Fig. 5 indicates good agreement, suggesting that the region we have identified as the Reynolds layer in our simulations based on its qualitative appearance has properties that are quantitatively consistent with those inferred from observations.

We can also compare the properties of the disc in our simulations to observed properties of the Solar neighbourhood (assumed to be located at a radius of ∼8 kpc in the simulations). In the simulations, we find an average Solar neighbourhood gas surface density of ∼20 M_⊙ pc⁻² at the final epoch (see Appendix  B), in good agreement with the ≈15 M_⊙ pc⁻² observed value (McKee et al. 2015). Similarly, we compute a mid-plane Solar neighbourhood field strength of |$\sim 5\, \mu$|G in the simulations, which is consistent with the magnetic field of |$\sim 6\, \mu$|G inferred from nearby diffuse H i clouds (Crutcher et al. 2010). We can also compare to the SFR in the Solar neighbourhood, but this requires some care because star formation is highly stochastic in time, and thus observational tracers that average over different timescales do not necessarily agree with one another. Misiriotis et al. (2006) find that the SFR surface density in the Solar neighbourhood is ≈2.5 × 10⁻³ M_⊙ pc⁻² Myr⁻¹ (using their values as re-scaled by Kennicutt & Evans 2012) based on the Galactic distribution of 100 μ|$\mu$|m emission, a technique that averages over ∼10−20 Myr timescales, since these are the lifetimes of the stars responsible for dust heating. Taking a 10 Myr average of our simulation, we find a Solar neighbourhood SFR ∼10⁻³ M_⊙ pc⁻² Myr⁻¹, again in reasonable agreement with observations.

4.3 Magnetic effects on supernova-driven winds

We next compare the vertical outflows of gas between the MHD and non-MHD versions of our simulations. We compute the mass outflow rate as

where ρ is the gas density, |$\boldsymbol {v}$| is the gas velocity vector, δV_z denotes the surface of a cylinder centred on the galaxy with radius R_max and heights ±z (above and below the disc, respectively), |$\hat{n}$| is the surface normal unit vector, and dS is the area element. We choose R_max to be 17.5 kpc, encompassing virtually all of the mass of the disc. In Fig. 10 (left-hand panel), we show the vertical mass flux going away from the disc as a function of height, comparing the hydrosimulation at t ∼ 1 Gyr with the MHD simulation at the same simulated time and also with the MHD simulation at the same SFR (at an earlier simulated time). When comparing the two versions at the same simulated time, we find a substantially lower mass outflow rate at all heights in the MHD simulation (suppressed by 1–2 orders of magnitude). However, when the MHD and non-MHD simulations are matched at the same SFR (approximately 1.65 M_⊙yr⁻¹), there is no appreciable difference between the mass outflow rates as a function of height. (The large bump at ∼4 kpc in the hydro outflow rate is due to a transient star formation event and has propagated outward as a traveling wave, an effect apparent in the time-dependent version of this figure.)

In Fig. 10 (right-hand panel), we examine the vertical energy flux as a function of height above the galactic disc for both the magnetic and non-magnetic simulations. We compute the total kinetic and thermal energy in the outflow, neglecting the magnetic energy for a consistent comparison between the magnetized and non-magnetized simulations. The results are qualitatively identical to the mass outflow results, with the MHD simulation having a suppressed energy outflow rate with respect to the non-magnetized outflow rate when the two are compared at the same time, but with this effect disappearing after controlling for the SFR.

We further examine the mass outflow rate by decomposing it into thermal phases, including the three phases of the neutral ISM (cold neutral, unstable, and warm neutral; e.g. Wolfire et al. 2003), the warm ionized medium (WIM), the warm-hot ionized medium (WHIM), and a hot phase comprised of material at all higher temperatures, as shown in Fig. 11. The boundary between cold, unstable, and warm phases is determined by finding the zero-crossings of the derivative of the equilibrium thermal pressure with respect to gas density dP/dn, as explained in Appendix  A. The boundary between warm neutral and WIM is determined by finding the temperature at which the number of free electrons is one-half the number of hydrogen nucleons in thermal equilibrium. For our cooling function, we find this occurs at a temperature of 7105 K.⁶ We set a somewhat arbitrary temperature threshold between the WIM and WHIM at 2 × 10⁴ K, and set the boundary between the WHIM and the ‘hot’ phase at 5 × 10⁵ K. These choices are designed to allow the WHIM to encompass the material near the peak of the radiative cooling rate Λ(T) of interstellar gas in collisional ionization equilibrium at T ∼ 10⁵ K (e.g. fig. 8 of Sutherland & Dopita 1993).

We find that the phase structure of the outflow is not qualitatively different between the magnetized and non-magnetized simulations, and that all phases decline rapidly with height in both simulations, suggesting a fountain type of outflow. Neutral phases dominate in the fountain region z ≲ 3 kpc, while ionized gas dominates at larger heights. The localized peaks and troughs apparent in the profiles are of the same order of magnitude as the temporal variability in the profiles on ∼ Myr time-scales and we caution against overinterpreting differences between the MHD and non-MHD simulations.

In Fig. 12, we show the (dimensionless) mass loading factor η, defined as

where |$\dot{M}_{\text{wind}}$| is the mass outflow rate, and |$\dot{M}_{\text{SFR}}$| is the SFR, as a function of the SFR |$\dot{M}_{\text{SFR}}$|⁠. The mass outflow rate is the instantaneous value, while the SFR is averaged over the previous 10 Myr. Each point represents the galaxy-averaged value a simulation output, and we compute points for each simulation snapshot (output at Δt = 1 Myr intervals), excluding times t ≲ 0.2 Gyr. To show that the trend shown in Fig. 12 is not dominated by the transient burst of star formation at the start of the simulation, we show the mass outflow rate at |z| = 2.0 kpc and the 10 Myr-averaged SFR as a function of time in Fig. 13; this illustrates that the ratio of vertical mass flux to SFR undergoes a secular decline with time as the SFR decreases due to gas consumption, with no qualitative change in behaviour beyond the initial 0.2 Gyr period that we exclude. We have also confirmed that the qualitative trend shown in Fig. 12 is unchanged if we exclude the initial 0.4 Gyr rather than 0.2 Gyr, or if we change the averaging interval in the range of 10−100 Myr.

The colour of the points in Fig. 12 indicates at what galactocentric heights ±z we computed the mass outflow. We observe a linear trend, i.e. the mass loading factor is roughly proportional to the SFR. This trend is the same regardless of whether magnetic fields are present in the simulation (solid boxes correspond to the MHD simulation, dotted boxes correspond to the hydro simulation). This implies that there is an approximately quadratic dependence of mass outflow rate on SFR, regardless of the height at which the outflow is measured. Consistent with the picture of a fountain outflow, the mean trendline for each set of measurements for various galactocentric heights shows that net mass outflow rate declines with height. This sharply disagrees with the scaling found by Muratov et al. (2015), who found a linear relationship between mass outflow rate and SFR (implying a constant mass loading factor with SFR; their fig. B1), although they measured the mass outflow rate at 0.25R_vir ≈ 50 kpc (comoving) and found a mass outflow rate of ∼10 M_⊙yr⁻¹ at an SFR of 1 M_⊙yr⁻¹ at redshifts z ∼ 2−4 for progenitors of Milky Way-mass galaxies. The significance of this disagreement is difficult to interpret, since our simulations are both non-cosmological and much higher resolution than those of Muratov et al. (2015), but it is possible that the scaling properties of galactic outflows may strongly differ when measured near the galaxy versus near the virial radius, or may be a strong function of redshift (or gas fraction).

The qualitative trends of our fountain outflows broadly agree with those found in the simulations of Kim & Ostriker (2018) and Kim et al. (2020). However, we find that the structure of our outflows is significantly more extended in the vertical direction, and that at ∼1 kpc scales, all phases contribute a net mass outflow on average, whereas the simulations of Kim & Ostriker (2018) find that only the ‘ionized’ and ‘hot’ phases contribute to the net outflow at scales ≳ 1 kpc. At a fixed height of ∼2 kpc and a comparable SFR, we also find a mass loading factor that is approximately an order of magnitude greater than found by Kim & Ostriker (2018) and Kim et al. (2020). This is a significant discrepancy and more work is needed in order to determine its cause. One possible explanation is that our simulations include a model for ‘pre-supernova’ feedback in the form of photoionization, while Kim & Ostriker (2018) and Kim et al. (2020) include supernovae as their only form of spatially localized feedback. A number of authors have found that including pre-supernovae feedback can significantly alter the properties of galactic winds, by transforming the environment into which supernova energy and momentum are deposited (e.g. Agertz et al. 2013, Kannan et al. 2020; Jeffreson et al. 2021). Another possible explanation is that the local box geometry of Kim et al. (2020) does not allow streamlines to open up, which prevents the outflow from reaching the sonic point in the classical superwind solution of Chevalier & Clegg (1985), as emphasized by Martizzi et al. (2016).

Kim et al. (2020) additionally find a negative trend of mass loading factor with SFR, although this relationship was obtained by fitting models with significantly varying initial gas surface densities, and the same negative trend is obtained by fitting the mass loading factors and the initial gas surface densities of their models (their fig. C1). At fixed gas surface density, they likewise find a positive power-law relationship between mass loading factor and SFR for most of their models, with the strength of this relationship varying with initial gas surface density (their Fig. 8). We leave a more detailed exploration of the outflow properties and scalings with galaxy parameters to future work.

4.4 Implications for cosmic rays and radio observations

Our simulated Galactic magnetic field structure has implications for the transport of cosmic rays within and out of the Galaxy. In many analytic models of cosmic rays in the Galaxy, there is a distinction between an inner region of ‘tangled’ field lines and an outer region several kiloparsecs above the Galactic mid-plane with large-scale coherent magnetic fields (e.g. in the hydrostatic model of Boulares & Cox 1990, the wind model of Breitschwerdt et al. 1991; the ‘base radius’ of cosmic ray-driven winds in Quataert et al. 2022). Our results in Section 4.1 indicate that this transition occurs at ∼3 kpc. Due to the theoretical uncertainties in the cosmic ray diffusion coefficient as used in these models, the resulting predictions for mass outflow and gamma-ray luminosity (e.g. Lacki et al. 2011) cannot be used to directly test the magnetic field structure of our simulations, but our results provide a justification for a ∼3 kpc value of the launching radius in a cosmic ray-driven wind model of the Galaxy.

Our results about the vertical structure of the magnetic field may be more directly tested via spatially resolved synchrotron emission, which is primarily sensitive to the magnetic field strength. Radio synchrotron observations typically find disk galaxies (of all inclinations) to have a scale length of ∼3−5 kpc (see Beck 2015 and references therein). Additional comparisons may be made to Galactic polarized synchrotron measurements (e.g. Planck Collaboration XLII 2016; Krachmalnicoff et al. 2018).

Another probe of the magnetic field are the polarized dust emission maps observed with the Planck satellite (Planck Collaboration XX 2015). Assuming the standard radiative torque alignment model (Davis & Greenstein 1951; Lazarian & Hoang 2007) and uniform dust properties across the Galaxy, these maps probe the plane-of-sky magnetic field orientation integrated along the line of sight. MHD simulations of local volumes of the ISM have been compared to the Planck maps (e.g. Planck Collaboration XX 2015; Kim et al. 2019) but comparisons with MHD simulations with star formation, supernova feedback, and a global disc geometry are currently lacking. The latter may be relevant due to the large-scale correlations inferred in the Galactic magnetic field (e.g. spiral arm structure). As a test of the realism of our simulations, however, such comparisons are somewhat limited by uncertain dust physics.

Undoubtedly the best observational comparison to our results will be the dense forest of quasar and pulsar Faraday rotation measures observable with the forthcoming Square Kilometre Array (SKA) to be completed in Western Australia (Haverkorn et al. 2015). These measurements of the line-of-sight Galactic magnetic field will enable direct tests of MHD simulations of the Galaxy via comparison with Faraday rotation maps and, in the plane of the Galaxy, tomographic mapping of the line-of-sight magnetic field.

5 CONCLUSIONS

Our first main result is the ‘layer-cake’ vertical structure of the magnetic field, evident most prominently in the Alfven Mach number and plasma beta as a function of Galactocentric height (Fig. 4). This structure decomposes into three approximate zones, with an innermost region having a roughly constant Alfven Mach number and plasma beta, an outermost region having a significantly larger but also roughly constant Alfven Mach number and plasma beta, and an intermediate region between 300 pc ≲ |z| ≲ 3 kpc smoothly connecting the two with a power-law structure. A similar zonal structure has been noted in the Galactic synchrotron emissivity, usually with a two-component structure of characteristic heights ∼200 pc and ∼1.5 kpc (e.g. Boulares & Cox 1990). Our simulation provides a theoretical picture that may help explain these observations.

Our second main result is the order unity effect of magnetic fields on the SFR (Fig. 2). We observe that the SFR is suppressed by a factor of 1.5−2 when magnetic fields of strength comparable to those observed in the Galaxy are present. However, when controlling for SFR, the mass outflow rate (both total and decomposed by phase) is indistinguishable between the magnetized and non-magnetized simulations (Figs 10 and 11). The mass outflow decomposed by phase is highly stochastic and is difficult to compare with precision between the simulations. Future studies are needed in order to quantify the residual differences at better than order-of-magnitude between outflows of magnetized and non-magnetized Galactic discs.

Lastly, we obtain a positive linear correlation between mass loading factor and SFR in our simulations when measuring the mass outflow on kiloparsec scales above the disc (Fig. 12). Previous work using a similar supernova feedback model found no correlation between these two quantities for gas-rich discs when measuring the mass outflow at significantly larger scales (Muratov et al. 2015). The apparent discrepancy should be examined by future work. In searching for a theory of this enhanced mass loading, it may be productive to examine the full distribution of gas densities in which supernova explode as a function of SFR, since this enters as an explicit factor in our adopted feedback model (equation 5) as a result of the density dependence of the radiative cooling of supernova remnants (Thornton et al. 1998).

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank L. Armillotta for providing the code to compute the H ii regions in our simulations and thank N. McClure-Griffiths for useful discussions.

This research was supported by the Australian Research Council through its Discovery Projects and Future Fellowship Funding Schemes, awards DP190101258 and FT180100375. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

The simulations in this work were run at a datacentre using 100 per cent renewable electricity sources.⁷

Software: matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), pandas (pandas development team 2020; McKinney 2010), numpy (Harris et al. 2020), h5py (Collette 2014), Snakemake (M’older et al. 2021), lic (Brinkmann 2020), meshoid (Grudic 2019), gizmo (Hopkins 2015), grackle (Smith et al. 2017).

DATA AVAILABILITY

Due to the large volume of data products produced (approximately 4.6 terabytes), a limited subset of the raw simulation outputs (and their associated processed data products) are permanently archived as an open-access Zenodo data set (Wibking & Krumholz 2021). Two simulation snapshots at different simulated times (t ∼ 556 Myr and t ∼ 976 Myr) are included from the MHD simulation and one simulation snapshot from the hydrodynamic simulation (t ∼ 976 Myr) is included in gadget-2/gizmo HDF5 format, along with processed data products from each of the simulation snapshots in NPZ (numpy array archive) format and visualizations of the processed data products as images in PNG format. Contingent on the ongoing availability of storage resources, additional simulation snapshots are available upon reasonable request to the authors.

Footnotes

References

APPENDIX A: COOLING CURVE COMPARISON

In Fig. A1, we show the equilibrium temperature and pressure as a function of density produced by the Grackle cooling code (version 2.2; Smith et al. 2017), as used in this work, and the gizmo cooling module as used in, e.g. the FIRE-2 simulations (Hopkins et al. 2018b), respectively.⁸ The unstable neutral medium is the phase for which dP/dn < 0, i.e. where the slope is negative in the right-hand panel of Fig. A1; the stable warm and cold atomic phases correspond the regions with dP/dn > 0 on the low- and high-density sides of this region, respectively. We see from the figure that the Grackle cooling curve features an unstable neutral phase between 980 and 4126 K, corresponding to a gas density of 0.26 and 0.38 H cm⁻³ and pressures of 302 and 866 K cm⁻³. The FIRE-2 cooling curve features an unstable neutral phase between 1091 and 8964 K, between densities of 0.02 and 0.06 H cm⁻³ and pressures of 51 and 181 K cm⁻³ for solar neighbourhood FUV irradiation.

Based on Wolfire et al. (2003), we expect an unstable phase at pressures between 1960 and 4810 K cm⁻³ between densities of 0.86 and 6.91 H cm⁻³ and temperatures of 258 and 5040 K for solar neighbourhood ISM conditions. Clearly, neither cooling curve agrees quantitatively with expectations. However, while for Grackle the unstable pressure range is a factor of ∼5 below that computed by Wolfire et al. (2003), for FIRE-2 the discrepancy is a factor of 30−50. Indeed, with the FIRE-2 cooling code, the diffuse interstellar medium (largely shielded from ionization due to young stellar populations) of galactic discs will be entirely composed of cold neutral medium, with no stable warm phase at all – with the FIRE-2 cooling, such a phase exists only for pressures P/k_B ≲ 50 K cm⁻³, which is far below those pressures found in galactic discs and instead is more typical of the circumgalactic medium.

The source of this discrepancy in the FIRE-2 cooling curve is the unphysically low grain photoelectric heating rate produced by the FIRE cooling module in the neutral atomic ISM. Specifically, the photoelectric heating efficiency [as defined by equation (20) of Wolfire et al. (2003)] is 3–4 orders of magnitude too low when computed by gizmo. This is a result of the model used for attenuating the ionizing flux from the extragalactic UV background, which sharply cuts off the flux above gas densities of ∼0.0123 cm⁻³, resulting in an unphysically low free electron density of ∼10⁻⁶ cm⁻³ in the warm neutral ISM. In the non-public FIRE-2 radiative transfer code, there are additional local photoionizing sources from young stellar populations, which partially mitigates this effect on the free electron number (but overestimates the photoionizing flux due to the optically thin approximation used for radiative transfer). However, the dominant source of ionization in the neutral atomic ISM of the Galaxy is not young stellar populations, but the combination of stellar EUV emission from old stellar populations (e.g. low-mass X-ray binaries) and the soft diffuse X-ray background produced primarily by X-ray line emission in supernova remnants (Slavin, McKee & Hollenbach 2000), with an additional contribution from C⁺ ionization in regions of high FUV irradiation or high density (Wolfire et al. 2003). These sources are not included in the FIRE-2 ionization model. Including these photoionizing sources in models (e.g. Wolfire et al. 2003) yields a free electron number density in the warm neutral interstellar medium consistent with the observationally inferred free electron density in the solar neighbourhood of ≈0.047 cm⁻³ (assuming a hydrogen nucleon density n_H = 0.5 cm⁻³; section 8.1 of Jenkins 2013).

Grackle version 2.2 does not include a self-consistent attenuation model for UV background radiation, so we did not enable it in our simulations, and our Grackle cooling calculations are therefore in the optically thin limit. As shown by the detailed radiative transfer calculations of Rahmati et al. (2013) (c.f. their fig. 3), optically thin extragalactic ionization remains a good approximation for calculating the ionization state of interstellar gas up to densities of at least ∼10 cm⁻³, because the attenuation of the extragalactic background is nearly compensated by an increase in ionizing flux due to diffuse galactic ionizing sources, although a precise accounting may require a more detailed calculation that both attenuates the extragalactic background and also includes the Galactic soft diffuse X-ray background (which Grackle does not include and which we do not attempt here). The equilibrium number of free electrons per hydrogen nucleon as a function of gas density for both Grackle and FIRE is shown in Fig. A2. We note that the Grackle equilibrium ionization predictions are broadly consistent with observations, whereas the FIRE equilibrium ionization model is not.

However, the Grackle ionization model is not perfect, and the location of the unstable neutral phase remains discrepant with the standard model of Wolfire et al. (2003). This discrepancy is likely due to the imperfect cancellation of the attenuation of the extragalactic background and the addition of the Galactic soft X-ray background (as discussed above), plus two additional details of the Wolfire et al. (2003) and Grackle models. First, the Wolfire et al. model strongly attenuates the extragalactic photoionizing UV background, such that in the Solar Neighbourhood, the extragalactic UV background only provides ∼1 per cent of the total ionization rate, which is probably too small for low-density gas, as suggested by the radiative transfer calculations of Rahmati et al. (2013). Secondly, the Grackle model neglects C⁺ ionization (also neglected by the FIRE model), which is the dominant source of free electrons at densities ≳30 H cm⁻³ for Solar metallicity gas in the CNM (equation C23; Wolfire et al. 2003), causing the number of free electrons per H to plateau at a value of approximately 2 × 10⁻⁴ under Solar Neighbourhood conditions. As shown in Fig. A2, while Grackle (unlike FIRE) produces a qualitatively reasonable estimate for the free electron abundance at densities of ∼30 − 100 H cm⁻³, does not correctly recover the plateau behaviour of the real ISM, because it does not include C⁺.

Since the photoelectric heating efficiency scales as the free electron number density to the 0.73 power (equation 20 of Wolfire et al. 2003), the ionization state can therefore make a difference of many orders of magnitude in the photoelectric heating rate. Due to this effect, in the diffuse interstellar medium (which should be far from ionizing radiation produced by young stellar populations), the FIRE-2 cooling model transitions between the warm and cold neutral phases at pressures and densities that are several orders of magnitude too low. In light of this qualitatively incorrect thermal structure, all conclusions regarding the thermal state of the interstellar medium in the FIRE-2 simulations (e.g. Gurvich et al. 2020; Pandya et al. 2021) should be critically re-examined.

APPENDIX B: SIMULATION PROFILES

In Fig. B1, we show the gas and magnetic scale heights of the MHD simulation at the intermediate (t ∼ 550 Myr) and the final snapshot (t ∼ 1 Gyr). We compute the scale height z₀ as

which exactly coincides with the exponential scale height z₀ for the case of an exponential profile |$\rho (z) = \rho _0 {\rm e}^{-z/z_0}$|⁠. The magnetic scale height is computed in a precisely analogous way, using the magnetic field strength profile |$|\boldsymbol {B}(z)|$|⁠.

In Fig. B2, we show the total gas surface density of the MHD simulation at the intermediate (t ∼ 550 Myr) and final snapshots (t ∼ 1 Gyr).

APPENDIX C: DEFINITION OF THE GAS VELOCITY DISPERSION AND RELATED QUANTITIES

Two of the quantities we compute in Section 4.1 – the Alfvén Mach number and the turbulent pressure – depend on the rms turbulent velocity of the gas. This not an entirely straightforward quantity to compute, because it depends on how we separate the velocity field into mean field and turbulent components. In Cartesian geometry this separation is often accomplished via a decomposition in Fourier space, but such an approach would not work in the cylindrical geometry of our simulation. More generally, unlike in the periodic boxes analysed most often in the turbulence literature, our system possesses large-scale gradients – for example variation of the galactic rotation speed with radius and height – that preclude a completely clean separation of the variation into large-scale and turbulent components.

Given this complexity, we adopt a definition of turbulent velocity that is appropriate for regions near the disc of the galaxy, where the velocity field can reasonably be described as ordered orbital motion plus a turbulent perturbation on top of it. We therefore define the turbulent velocity |$\delta \boldsymbol {v}$| by examining the variation of the velocity in azimuth ϕ at a given galactocentric radius R and height z:

where i denotes a velocity component in the set {R, z, ϕ} and 〈 · 〉_ϕ denotes an average over azimuthal angle ϕ. The rms velocity fluctuation δv is then computed as

For the purposes of computing the Alfvén Mach number, we use this method to estimate the turbulent velocity vector |$\delta \boldsymbol {v}$| at every point. We then compute the ratio |$|\boldsymbol {v_A}| / |\delta \boldsymbol {v}|$| before averaging in space; the result is our estimate for the mean Alfvén Mach number. Similarly, we compute the turbulent pressure at every point as |$\rho |\delta \boldsymbol {v}|^2$|⁠, and average this quantity to obtain the mean turbulent pressure.