Synchrotron Emission on FIRE: Equipartition Estimators of Magnetic Fields in Simulated Galaxies with Spectrally-Resolved Cosmic Rays

Synchrotron emission is one of few observable tracers of galactic magnetic fields ( B ) and cosmic rays (CRs). Much of our understanding of B in galaxies comes from utilizing synchrotron observations in conjunction with several simplifying assumptions of equipartition models, however it remains unclear how well these assumptions hold, and what B these estimates physically represent. Using FIRE simulations which self consistently evolve CR proton, electron, and positron spectra from MeV to TeV energies, we present the first synthetic synchrotron emission predictions from simulated L ∗ galaxies with "live" spectrally-resolved CR-MHD. We find that synchrotron emission can be dominated by relatively cool and dense gas, resulting in equipartition estimates of B with fiducial assumptions underestimating the "true" B in the gas that contributes the most emission by factors of 2-3 due to small volume filling factors. Motivated by our results, we present an analytic framework that expands upon equipartition models for estimating B in a multi-phase medium. Comparing our spectrally-resolved synchrotron predictions to simpler spectral assumptions used in galaxy simulations with CRs, we find that spectral evolution can be crucial for accurate synchrotron calculations towards galactic centers, where loss terms are large.


INTRODUCTION
Magnetic fields (B) and relativistic charged particles (cosmic rays, hereafter CRs) are of significant importance in astrophysics.Both are known to provide a significant source of non-thermal pressure support in the interstellar medium (ISM) (Boulares & Cox 1990), influence the physics of giant molecular clouds (for relevant reviews, see Crutcher 2012;Beck et al. 2020), and magnetic fields determine the transport of the dynamically coupled CRs through the ISM and into the CGM (Zweibel 2013), therefore influencing the structure of galactic outflows (see Zhang 2018;Owen et al. 2023;Ruszkowski & Pfrommer 2023, for relevant reviews).
A major difficulty in our understanding of galactic magnetic fields, and subsequently on our understanding of CR propagation and ensuing effects, is the measurement of magnetic field strengths and geometries.A unique extragalactic probe of B and CRs comes from the synchrotron emission radiated by CRs as they gyrate around magnetic field lines (Ginzburg & Syrovatskii 1965).Most estimates of magnetic field strengths in dwarf and spiral galaxies (Fitt & Alexander 1993;Beck et al. 2000;Chyzy et al. 2011;Fletcher et al. 2011; Basu ★ E-mail: sponnada@astro.caltech.edu& Roy 2013; Beck 2015) are derived indirectly from the intensity of this non-thermal synchrotron emission, often at radio frequencies, with the key assumption of energy equipartition between CR protons and B, used to resolve their formal strict degeneracy.
As discussed in Beck & Krause (2005); Stepanov et al. (2014); Seta & Beck (2019), the energy equipartition method is motivated by the intertwined dynamics of cosmic rays and magnetic fields, in approximate pressure equilibrium in the local, warm ISM.While this technique has been applied to several radio observations of galaxies across a wide range of spatial scales (Beck et al. 2020), it is not clear whether B and CRs are in approximate pressure equilibrium at all spatial scales, and across the large dynamic range of ISM gas densities and temperatures.Additionally, the method is subject to further assumptions regarding the spectrum, the spatial distribution of CR electrons (CRe) and protons (CRp), as well as the effective emitting volume, which are potentially order-of-magnitude uncertain given the multi-phase nature of the ISM.
Given these assumptions and caveats, it is unclear how to interpret the equipartition estimates of B. B is known to vary by orders of magnitude with gas density and with ISM phase (Crutcher 2012;Han et al. 2017), and while equipartition-inferred estimates of B may inform us about some volumetric and LOS-averaged estimator of B, it remains unclear how this maps to measuring B within a given ISM phase or gas density.In other words, what phase of the ISM dominates the synchrotron emission, and how well does the equipartition estimate of B trace it?
Exploring questions about synchrotron emission and equipartition in numerical simulations has been limited in scope by the computational complexity of self-consistently evolving B and CRs, in the context of the multi-phase ISM, which exhibits great spatial and temporal variation in gas properties across cosmological time.Simulations have only recently been able to evolve magnetic fields (Peng & Tom 2009;Pakmor et al. 2014;Marinacci et al. 2014;Su et al. 2017;Rieder & Teyssier 2017;Butsky et al. 2017;Martin-Alvarez et al. 2018;Ntormousi et al. 2020) including CRs as a coupled fluid term (Booth et al. 2013;Girichidis et al. 2016;Butsky & Quinn 2018;Chan et al. 2019;Buck et al. 2020;Werhahn et al. 2021a,b,c;Hopkins et al. 2020;Farcy et al. 2022), though it varies whether the simulations in these studies are cosmological, in addition to differences in resolution, physics prescriptions, and numerical methods.Furthermore, these simulations which incorporate a fluid treatment of CRs often take the so-called "single-bin" approximation, evolving solely the ∼1-10 GeV CR proton energy density (or a constant underlying spectrum, with transport coefficients corresponding to the 1-10 GeV CR protons) and by construction fail to capture the CRe spectra, which are required for predictive synchrotron calculations.
Recent algorithmic advances have pushed this boundary of explicitly modeling the CR spectrum (Yang & Ruszkowski 2017;Girichidis et al. 2020;Ogrodnik et al. 2021), allowing us to evolve the spectra of various CR species in live-kinetic magneto-hydrodynamic (MHD) simulations of galaxy formation (Hopkins et al. 2022a).Additionally, detailed study of B within this family of simulations (which share the same physical prescriptions, up to the novel spectral treatment of CRs) has shown that they produce realistic B strengths and geometries in resolved ISM phases (Ponnada et al. 2022).
In this work, we forward-model synchrotron emission from highresolution simulations of galaxy formation with cosmological initial conditions.We present our methodology for computing synchrotron emission from a set of L * galaxies from the Feedback in Realistic Environments project (FIRE)1 simulation suite (Hopkins et al. 2018(Hopkins et al. , 2022b) ) in Section 2. In Section 3, we compare our results for the spectrally-resolved FIRE runs to relevant observations, and explore how equipartition estimates of magnetic field trace the underlying magnetic fields within these simulations.Through resolving multiphase ISM structure in our simulations, we disentangle how different phases of the ISM contribute to the resulting synchrotron intensity, and physically interpret equipartition measurements of B, discussing how well each of the fundamental assumptions in empirical equipartition models does or does not apply.In Section 4, we present a toy model for understanding equipartition magnetic field estimates in the context of a multi-phase ISM motivated by our results in the prior section.Finally, we discuss our results and conclusions and summarize our findings in Section 5.

METHODS
In this section, we will discuss the details of the simulations used in this study, and describe our methodology to compute synchrotron emissivities from the simulations in post-processing.
We investigate three zoom-in simulations of galaxies roughly similar in mass and size to the Milky Way, named m12i, m12f, and m12m.
These simulations were shown in previous works to all produce CR spectra consistent with those local ISM (LISM) constraints (Hopkins et al. 2022a).
All three simulations here have been presented and extensively studied in Hopkins et al. (2022a) (which did not, however, model their synchrotron properties), to validate that many other properties (e.g.stellar and gas masses, sizes, kinematics, etc) are plausibly consistent with observed galaxies of similar masses.The simulations all have a Lagrangian mass resolution of 56000 M ⊙ for gas cells, and the typical spatial resolution ranges from ∼ 1 − 10 pc in dense gas.

Simulations
The simulations studied here are all fully-dynamical, cosmological, magnetohydrodynamic-radiation-thermochemical-cosmic raygravitational star and galaxy formation simulations.This means they self-consistently follow galaxy formation from cosmological initial conditions at redshifts > 100 including both dark matter and baryons (in gas and stars), with magnetic fields grown self-consistently from arbitrarily small trace seed fields at  ≈ 100, with phase structure and thermo-chemistry in the galaxy emerging from cooling with temperatures  ∼ 1 − 10 10 K and self-gravity, fully-coupled to multi-band (EUV/FUV/NUV/OIR/FIR) radiation transport, with star formation in the most dense gas (≳ 1000 cm −3 ), and those stars influencing the medium in turn via their injection of radiation, stellar mass-loss, and both Type Ia and core-collapse SNe explosions (followed selfconsistently according to standard stellar evolution models).The cosmic ray physics is itself coupled directly to the dynamics, with CRs propagating along magnetic field lines according to the fully general CR transport equations, and interacting with the gas via scattering, Lorentz forces, and losses.
This means that when we model CRs, all quantities needed to compute the fully non-equilibrium CR dynamics and losses are captured in-code, except for the microphysical CR scattering rate .This arises physically from CR pitch-angle scattering off magnetic field fluctuations on gyro-resonant scales -far smaller than we can possibly resolve (∼ 0.1 au in the warm ISM, for ∼ GeV CRs).As such, the simulations must insert some assumed "sub-grid" scattering rate.We stress that the more familiar CR "diffusion coefficient" and/or "streaming speed" arise self-consistently from  in the appropriate limits of the CR dynamics equations (for example, if the CR distribution function is sufficiently close to isotropic and in local flux steady-state, the effective parallel diffusion coefficient is simply  ∥ ∼ v2 cr /3).As mentioned earlier, the simulations presented in Section 3 are the same as those presented in Hopkins et al. (2022a), and so we refer the reader for further details therein.Here, we summarize the most pertinent details.These simulations are run with GIZMO 2 , in the mesh-free finite-mass mode.All simulations include MHD as treated in (Hopkins & Raives 2016;Hopkins 2016), and fully-anisotropic Spitzer-Braginskii conduction and viscosity (Hopkins 2017;Su et al. 2017).These simulations were restarted at  ≈ 0.05 from a highresolution, cosmological, FIRE-2 'single-bin' CR-MHD simulation of the same galaxy, using the self-consistently evolved CR energy densities in each cell to populate the CR distribution function, and then run for ∼500 Myr.Starting from this low redshift and evolving for this runtime are sufficient for the galaxies CR spectra to reach quasi steady-state behavior in the disk and inner CGM at  = 0.This is true also of secondary leptons, which reach a steady-state on their dominant loss timescale (Coulomb/ionization, diffusive, and radiative losses for ∼MeV, ∼GeV, and ≳ 50 GeV leptons respectively), which never exceed ∼10 Myr in the disk and inner halo, as shown in numerical tests of (Hopkins et al. 2022a).The details of prescriptions for star formation, stellar feedback, CR-MHD spectral evolution, and coupling to gas follow those of Hopkins et al. (2022b).
Our implementation of the CR physics is described in detail in (Hopkins et al. 2022a), and self-consistently includes adiabatic, diffusive re-acceleration, gyro-resonant loss, Coulomb, ionization, hadronic and other collisional, radioactive decay, annihilation, Bremstrahhlung, inverse Compton (IC), and synchrotron loss terms.Hadronic losses are assumed to be dominated by the proton-proton interaction, with total pion loss rates following Mannheim & Schlickeiser (1994); Guo & Oh (2008) and those of Evoli et al. (2017) for antimatter.CRs are injected via a power-law spectrum in momentum at SNe (Types Ia & II) and stellar winds (OB/WR) to neighboring gas cells with fixed fractions  inj CR = 0.1 and  inj e = 0.002 of the initial ejecta kinetic energy going into CRs (protons) and leptons.These injection fractions are well motivated by theoretical work on nonlinear diffusive shock acceleration and inverse modeling of observed CR spectra (Caprioli 2012;Yuan et al. 2012).Most notably, the fiducial simulations here differ in their treatment of CRs from those presented and analyzed in Hopkins et al. (2020) and Ponnada et al. (2022) by explicitly evolving each bin of the injected CR spectra, following the method of Girichidis et al. (2020) and Ogrodnik et al. (2021) as presented in Hopkins et al. (2022a).These simulations follow standard practice and assume a spatially and temporally constant scaling for the scattering rate as a function of CR rigidity  ∼ 10 −9 (R/GV) −1/2 s −1 (calibrated explicitly therein to fit all of the observations of CRs in the Milky Way from observations such as Voyager, AMS-02, Fermi, and others).

Synchrotron Forward Modeling
Our methodology to compute the synchrotron specific emissivities from our simulations follows the equations of synchrotron emission as derived in Ginzburg & Syrovatskii (1965), and summarized again in Padovani et al. (2021).
For each gas cell in our simulations, we calculate synchrotron emissivities as follows.First, we extract the internally evolved CRe and positron spectra,   () and components of the magnetic field perpendicular to the line of sight, B ⊥ .
Then, we compute the critical frequency of emission for each spectral bin of CRe, where B ⊥ is the magnetic field strength perpendicular to the line of sight, m e is the electron mass, and c is the speed of light.
We also calculate the power emitted per unit frequency by CRe in each spectral bin from MeV to TeV energies both parallel and perpendicular to the LOS, (2) where  ⊥ = |B ⊥ | and x = /  .
The functions F(x) and G(x) are defined as and where K 5/3 and K 2/3 are the modified Bessel functions of order 5/3 and 2/3.For each bin, we use the relevant F(x) and G(x) from precomputed look-up tables provided by Padovani et al. (priv. comm.).
We then compute the linearly polarized specific emissivities by integrating the contributions over j  , and where v e is the electron velocity.
From these specific emissivities, we also compute the Stokes Q  and U  specific emissivities, and where  is the local polarisation angle given by the orientation of B ⊥ rotated by ± 90 • .When calculating LOS-integrated quantities, we project the galaxy face-on or edge-on using the angular momentum vector of the stars to define the direction perpendicular to the galactic disk, which we define as the ẑ direction.All cell vector fields (r, B, v) are transformed accordingly.To compute the corresponding images of the specific intensity I  , Q  , and U  , the respective emissivity terms are integrated along the line of sight using a projection routine first described in Hopkins et al. (2005).This routine appropriately computes the contribution to the emission from every gas cell along the line of sight by taking into account their spatial distribution and hydrodynamic smoothing lengths.

RESULTS: FIRE SIMULATIONS WITH RESOLVED
COSMIC RAY SPECTRA

Synthetic Observations
In this section, we present our synthetic observations and explore which ISM phases contribute to the resultant emission.We first show images of I  at  = 6.2 cm for face-on projections of the galactic disks for the three fiducial simulations, following the procedure outlined in Section 2.2, in Figure 1.The images are in their fiducial highresolution format at a scale of ∼120 pc/pix and not smoothed to an observational beam.
A few things are evident from visual inspection: first, our synthetic synchrotron images have typical values of I  in qualitative agreement with those observed in nearby spiral galaxies (Basu & Roy 2013;Beck 2015) with spatially resolved radio continuum observations; we defer quantitative comparisons to the observations to Section 3.2.Second, the synchrotron emission is broadly coincident with the spatial distribution of neutral gas, with stronger emission coming from dense gas near the galactic center and correlated with spiral structure, and weaker emission in inter-arm regions and towards the galactic outskirts.
In Figure 2, we present 2-D histograms of the gas density and temperature, for the gas cells producing the emission visualized in Figure 1, all weighted by each gas cell's contribution to the specific synchrotron intensity I  .In these intensity-weighted ISM phase diagrams, the bimodal nature of gas which dominates the emission becomes clear; the synchrotron emission primarily arises from the cold, neutral medium (CNM -T ∼100-300 K, n ∼3-30 cm −3 ) and warm, neutral medium (WNM -T ∼3×10 3 -10 4 K, n ∼0.1-1 cm −3 ).Inspecting similar intensity-weighted histograms of the gas cells' neutral hydrogen fractions confirms that this gas is primarily neutral.
This evinces a physical scenario in which relatively cool, dense gas, which comprises a small fraction of the total ISM volume, contributes significantly to synchrotron emission.While warm, more appreciably volume-filling ISM gas is synchrotron-bright as well, the emission arises largely from the denser and neutral warm gas as opposed to the more diffuse (n < 0.1 cm −3 ), ionized warm/hot gas, which fills the majority of the ISM volume.Put more quantitatively, neutral gas with (n > 0.1 cm −3 , T < 10 4 K) contributes ∼80% of the emission, despite filling only ∼20% of the volume.
We note that a scenario in which CRs are not modeled selfconsistently would predict a different intensity-weighted phase distribution: if one assumed a priori that e CR was locally in equipartition with a self-consistently evolved B, then the intensity-weighted distribution would be more biased towards higher density gas as B ∼   , with  ∼ 0.5-0.66(Ponnada et al. 2022).Whereas if one assumed the ISM to be a homogeneous slab with a constant volumetric B in equipartition with CRs, the distribution would then effectively look like the volume-weighted distribution and thus weighted towards the more diffuse phases of the ISM.This result subsequently has physical implications for what traditional equipartition models à la Beck & Krause (2005); Lacki & Beck (2013) infer in the form of the "mean B" from synchrotron observations, which we show and discuss in Sections 3.4, 4, and 5.

Radial intensity profiles and comparisons to observations
In Figure 3, we compare the radial profiles of synchrotron specific intensity at observing frequencies of 0.33 GHz from our fiducial simulations to a few observed face-on spiral galaxies of roughly similar mass.We have also calculated the simulations' emission profiles at 4.8 GHz, as one of the observed galaxies, IC342, was measured at this frequency.But since the relative comparison is similar, we simply plot all systems at 0.33 GHz, shifting IC342 according to the authors' assumed spatially constant spectral index of -1.The radial profile of the non-thermal emission as plotted in Beck (2015) for IC342 is also shown, and for the other observed galaxies, we show the radial profiles of the non-thermal emission (corrected to remove the free-free emission) out to the radii published in Basu & Roy (2013).
We find the predicted I ,  = 4.83GHz for our simulations generally ranging from 10 5.5 to ∼10 4 Jy/sr at the galactic center to the "Solar circle" (galacto-centric cylindrical radius R = 7-9 kpc) and I ,  = 0.33GHz ranging from 10 7 to ∼10 5 Jy/sr in Fig. 3.When comparing our predicted radial emission profiles to observations of nearby spiral galaxies, we see an order-of-magnitude agreement, though two of the four observed profiles (NGC5055 and NGC 6946) appear to fall more slowly at large R ≥ 5 kpc.We caution that none of these simulations are meant to be exact analogs of the observed galaxies, but are only order-of-magnitude similar in morphology and galaxy mass.Various galaxy-to-galaxy differences in the profiles partly reflect this set of galaxies' diversity in properties like gas surface density and star formation rate, and in Figure A1 we show similar specific intensity profiles normalized to the gas surface density, which are qualitatively more similar in the shape of the profiles.

Spectral variation and emission properties
Evolving the full CR spectrum is computationally expensive, and most simulations of galaxy formation including CRs utilize the 'single-bin' approximation, solely evolving the 1-10 GeV e CRp .We are therefore motivated to test whether explicitly evolving CR spectra makes significant quantitative or qualitative differences in the predicted synchrotron emission.
To compare what the emission would look like for the same simulated galaxies if we did not evolve the full spectral information of CRs, we illustrate a few different spectral assumptions in Figure 4 starting from our fiducial simulations of m12i and m12f.First, we see that excluding emission from positrons (light-blue, dot-dashed lines) deducts minimally, contributing ≲ 10% of the emission at most radii, and reaching ∼20% near the galactic center of m12i.
Secondly, we show the results corresponding to scaling the CR electron spectrum of Bisschoff et al. ( 2019) by e CRe, sim , thereby holding the shape of j e constant (purple dotted lines), as well as scaling according to e CRp in each gas cell, which additionally fixes K (pink dashed lines), the total energy density ratio e CRp /e CRe to the literature value (akin to what might be done in the "single-bin" scenario where only e CRp is known).Holding the spectral shape constant with the "right" K shows that one would tend to over-predict the emission by up to factors of ∼6 in regions where losses at the energies of interest are strong (towards the galactic center).Additionally fixing K, we see that evolving the CR spectral information can make a modest difference at a factor of ∼2 in I  for typical spiral galaxy conditions as shown in much of m12f and the outer radii of m12i, but can make large (factor of ∼10-50) differences towards the central regions of the galaxy.
We break this down in more detail, showing that these deviations owe to a radially-varying proton-to-electron ratio between the galactic center and outskirts compounding with changing spectral shape (comparing constant spectral shape and constant K + spectral shape lines, which differ only in normalization) to further over-predict the emission, rather than the varying proton-to-electron ratio alone (green dashed-dotted lines, constant K, j e, sim ).This effect arises in m12i due to a quasi-starburst in the circum-nuclear region R < 3 kpc, where Σ gas quickly rises from ∼30-250  ⊙ pc −2 in comparison to lower surface densities ∼30-70 and 30-100  ⊙ pc −2 in m12f and m12m respectively.Correspondingly, loss terms owing to higher gas densities being correlated with higher B, increased radiation field intensities, etc., result in significant cooling of the CRe spectra at the energies of interest and a changing K.
As expected, there is little difference at "Solar circle" radii, where our simulations have been shown to faithfully agree with the observed LISM spectrum (Hopkins et al. 2022a).However, these results illustrate the need to self-consistently model the CR spectral shapes and K to accurately predict synchrotron emission in galactic environments where loss terms become important, as simple spectral assumptions from a 'single-bin' treatment of CR protons can lead to order-ofmagnitude differences in the predicted emission.These differences may also vary with frequencies different from the ∼GHz observa-  tions probed here which would probe different energy intervals of the CRe spectrum.Furthermore, evolution in K has implications on key assumptions invoked in equipartition estimates of B, as we detail in the next section. .Our synthetic synchrotron images for the fiducial model agree to within an order-of-magnitude with spatially resolved synchrotron emission from nearby spiral galaxies.Variation in the shapes of the profiles towards the galactic centers and outskirts partially arises due to differences in the gas surface density (see Figure A1 for corresponding plot), though the simulations appear to produce steeper radial profiles in comparison to the observed systems compared.) to internally evolved e CRp (pink dashed), holding the spectral shape constant by re-scaling the fiducial spectrum by e CRe (purple dotted), and holding K constant by re-scaling the fiducial j  to the nominal value at that simulation's "Solar Circle" (teal).Emission from secondaries contributes fractionally to overall emission at all radii, though reaching ∼20% contribution at high Σ gas .As demonstrated in m12i, strong CRe cooling at relevant energies corresponds to over-predicting emission towards the galactic center, where Σ gas is relatively higher, for constant K and/or spectral shape assumptions.This effect driven by variation in K compounding with spectral shape effects as shown by the pink line.Thus, in galactic environments where losses become important, modeling the CR spectra and K explicitly is necessary for predictive calculations, while in environments where cooling is relatively weak at relevant energies (galactic outskirts and much of m12f), spectral shape/K assumptions make at most factor of ∼2 differences.

Equipartition magnetic field strengths: what do they really measure?
The most commonly used formalism to determine equipartition estimates of the magnetic field (B eq ) from radio continuum observations of non-thermal emission from galaxies is that of Beck & Krause (2005).The equation is as follows: where  is the synchrotron spectral index, E p is the proton rest energy, and c 1 , c 2 , and c 4 are combinations of physical constants which encapsulate dependencies on  and the inclination of the magnetic field.The equipartition formula also requires assumptions about the depth of the emitting material, L (here the implicit dependence of the volume-filling factor of emitting gas, f  , is written explicitly), and K 0 , the ratio of number densities of CRp and CRe at energies from E p up to some energy E syn beyond which synchrotron and IC losses dominate for CR electrons3 (Beck & Krause 2005).
As described in the Introduction, equipartition estimation of the magnetic field is subject to a few extra assumptions beyond that of e CR = u B .The first auxiliary assumption is that of a constant ratio of the number densities of CRe and CRp, K 0 ∼100, motivated by the injection spectrum of primary CRs from SNe via diffusive shock acceleration and measurements at the Milky Way Solar Circle (Bell 1978(Bell , 2004;;Beck & Krause 2005).Though there have been modifications to this assumption for galaxies where the synchrotron emission is expected to have significant contribution from secondary CRe and positrons generated from the resultant pion-decays of primary CR collisional losses (Lacki & Beck 2013), these are not generally applied to galaxies like those in Fig. 3.This assumption further implicitly requires that the CRp and CRe spectrum to have constant and equal power-law indices, which holds close to injection sites, but may not hold in a spatially and temporally independent manner across galaxies.It also assumes negligible contribution from positrons, and that the synchrotron is optically thin.
Secondly, CRe and positrons, which dominate the synchrotron emission at frequencies of a few GHz, and ∼1-10 GeV CR protons, which dominate overall the energy density, are assumed to have the same, spatially and temporally uniform distribution within and across each pixel, and as a function of height above the mid-plane up to some height L. These assumptions are subject to questioning on galactic scales, where the dominant energy loss terms and correspondingly the loss timescales for CRe and CR protons can differ significantly due to large local variations (along a line of sight and across an observational beam) in gas density, magnetic and radiation energy densities, and phase structure (Wolfire et al. 1995;Evans II 1999).
Thirdly, the size of the emitting region both along the line of sight and within a given pixel or beam must be assumed.While traditional applications of the equipartition formula to galaxies have assumed path lengths of L ∼1-2 kpc for face-on observations, it remains unknown what the typical volume filling factor of synchrotron emitting gas actually is in galaxies.Implicit within this assumption is a homogeneous and volume-filling field, which given the multi-phase nature of the ISM, may not reflect the true B that primarily contributes to the synchrotron emission.
From the synthetic synchrotron images and radial profiles shown in the previous sections, we can now explore equipartition magnetic field strengths (B eq ) from the forward-modeled specific synchrotron intensity and test the several assumptions invoked.
In Figure 5, we show estimators of B weighted by I ,  = 0.33GHz and volume in m12f, and how B eq compares.We show B eq resulting from assuming  = 1, L = 1 kpc, K 0 = 100, and f  = 1, which are typical model assumptions in observational literature.In this case, it becomes immediately clear that Equation 10 under-predicts the "true" I  -weighted B by ∼0.3-0.6 dex across all radii and e CR by ∼0.3-1.0 dex for R > 4 kpc (though interestingly, traces the volumeweighted B well, which we discuss in detail in Section 4).Examining the assumptions for K 0 and f  in this model, we find K 0 varies from ∼100 ± 10 near the galactic center to ∼60 ± 10 near the "Solar circle" with the exception of m12i, where K 0 rises rapidly to ∼10 3 near the galactic center due to strong leptonic losses at relatively higher gas surface densities ≳ 150 M ⊙ pc −2 .Varying the assumed value of K 0 according to this radial variation would only affect the equipartition-inferred values at the tens of percent level (B eq ∼K 0 1/4 ), and so we find that K 0 = 100 is a decent assumption in our simulations for most of the galaxy conditions sampled here, with the exception of the very inner (R < 2 kpc) region of m12i.
The assumption of the depth of the emitting material being ∼1 kpc is more suspect, however.Computing the volume filling fraction of the gas cells which contribute to the upper 50% of I  within radial bins of vertical thickness 1 kpc reveals that f  varies from ∼0.05 -0.2 over the radial range.This volume filling fraction is very similar to the ratio of the face-on scale height of emission to the path length, H   /L, implying that the emission is primarily arising from the mid-plane.This, in conjunction with Figure 2, further develops the picture of the emission primarily being dominated by the cold and warm neutral medium relatively confined to the thin disk.
When f  is subsequently corrected to representative values ∼0.05 -0.2, the equipartition formula gives values of B eq and e CR that are closer to the "true," I  -weighted values of B (though over-predicting I  -weighted e CR owing to deviations from physical equipartition with u B ).The I  -weighted B is generally higher (0.2-0.6 dex) than the volume-weighted B, while the I  -weighted e CR and volumeweighted e CR exhibit less difference.This indicates that the primary effect of correcting f  is to correct for the higher B in the denser mid-plane gas, rather than to correct for e CR , which more weakly scales with gas density in our simulations owing to CRs' ability to diffuse (Ponnada et al. 2022;Hopkins et al. 2022a).It is important to note that the volume filling gas at heights L ∼ 1 kpc above the mid-plane at R ∼ 5 kpc has typical B ∼2-4 G ≪ ⟨B synch ⟩, (see Ponnada et al. 2022, for details), so the traditional estimator also severely over-estimates typical B in more diffuse gas above the disk, via assumptions of a homogeneous, height-independent B.
We emphasize that this is consistent with state-of-the-art, highresolution radio continuum observations (Krause et al. 2018;Heesen et al. 2018); detailed studies of nearby edge-on spiral galaxies have revealed very bright thin-disk components with scale heights ranging from 10s-100s of pc, which are ≲ the observational beam size and thus subject to large errors.The best-fit two-component models to these observations all indicate that this thin, mid-plane component almost ubiquitously dominates the overall emission, with additional extended thick-disk components with nominal scales of L ∼kpc as invoked in equipartition models contributing fractionally to the overall emission.Thus, when viewed face-on, the emission largely traces the thin disk component, and so taking typical values for L associated with these extended synchrotron halos rather than the scale height of the medium which contributes most to the emission (which is still subject to a high degree of observational uncertainty), leads to vastly under-predicting the I  -weighted B.
We have checked that treating edge-on synchrotron images of our simulations in an observational manner, smoothing to a representative observational beam of 10", and fitting two-component exponential disk models akin to observational studies yield similar results for thin and thick disk scale heights (∼100 pc, ∼1 kpc respectively; see Figure A2).

A TOY MODEL FOR ESTIMATING B FROM SYNCHROTRON EMISSION IN A MULTI-PHASE MEDIUM
In this section, we present a toy model to characterize B in a multiphase medium to gain intuition as to what factors cause the traditional model to deviate from the "true" values, given the understanding from the previous section that the synchrotron emission is primarily arising from neutral mid-plane gas (up to variations in the CR transport physics).
Consider the mean, volume-weighted magnetic energy density in a vertical ISM slab (at a given cylindrical radius R) of some finite height to be described by an exponential function anchored to the mean mid-plane (||≤ 0.1 kpc) value of u B : where the radial dependencies of ⟨ B 0 ⟩ and H B (the magnetic scale height) are implicit.Correspondingly, we can describe the volumeweighted average of the CR energy density in each slab: where  is ⟨ CR 0 ⟩/⟨ B 0 ⟩ and accounts for differences in the energy density at the mid-plane 4 and  accounts for differences in the vertical scale height of u B and e CR and the local inter-dependence.
Within each vertical slab, we can make the assumption that the coherence length of the magnetic field is much less than the vertical scale height or the radius, and thus construct a volume-weighted PDF of given magnetic fluctuations (   =  B / B 0 ) describing the probability of a given unit volume having a magnetic fluctuation    .We further make the ansatz that this PDF is a log-normal distribution of fluctuations in    , motivated by descriptions of density fluctuations in MHD turbulence (Hopkins 2013;Beattie et al. 2022): where V tot is the total volume within a vertical slab, S  describes the volume-weighted variance of ln(   ) arising from turbulent, multiphase structure.With these equations in hand, we can build a model which relates the observable synchrotron intensity (averaged over some finite pixel resolution A) to the underlying volume-weighted B.
Starting with the same assumptions as in Beck & Krause (2005) and including spectral dependencies implicitly ( will remain a free parameter), one can work to an expression for ⟨I  ⟩ given Equations 11 and 12 as follows:  10to the radial profile presented in Figure 3 is shown, with L = 1, K 0 = 100, f V = 1 (purple dashed) and f V computed within each radial bin (teal dashed).Shaded regions show the approximate ± 1 scatter (32-68 percentile) of the emission-weighted quantity and fiducial equipartition model at a given radial bin.Unilaterally, we see that without correcting for the volume-filling factor of the neutral mid-plane gas which dominates I  , the equipartition formula under-predicts the "true" I  -weighted B by ∼0.3-0.6 dex and e CR by by ∼0.1-0.6 dex.Furthermore, due to the emission being dominated by mid-plane gas, the equipartition model with the right volume-filling factor is generally not representative of the volume-weighted u B , though the fiducial model surprisingly traces the volume-weighted quantities well owing to a confluence of factors (see Section 4).While not shown in this figure, m12i and m12m show generally the same behavior. .
where dV is a volume element, C() is a combination of physical constants with minor dependencies on the power-law index and  5 , and  = [ (+1) 2 + ].From this expression, we can utilize the PDF given by Equation 13 to find a closed form expression for ⟨  ⟩ as follows: which can be re-arranged to find an expression for ⟨ B 0 ⟩: This expression is analogous to the Beck & Krause (2005) formalism of Equation 10, but relaxes the assumptions of equipartition between u B and e CR and a homogeneous medium with uniform B and

𝑝
, where c 1 -c 4 are defined in Appendix A of Beck & Krause (2005) e CR ; instead, we parameterize our ignorance of equipartition into a simple power-law relation between the two quantities and variations due to inhomogeneity in the ISM in the form of the volume-weighted variance in    and the magnetic scale height H B .Our revised expression reduces to Equation 10 when S  = 0,  = 1, and  = 1.
Note that our Eq.16 for ⟨ B 0 ⟩ = ⟨ B ⟩  is identical to the Beck & Krause (2005) (BK05) formula in Eq. 10 with the replacement  V →   cl  B / (BK05 implicitly take  V = 1), where We have three terms here parameterizing three physical assumptions/uncertainties: (1)  represents the (mean) deviation from equipartition; (2)  cl represents the "clumping factor" which can boost emission (for a given volume-weighted set of properties) owing to substructure with larger magnetic and/or CR energy density; and (3)   / simply corrects the ad-hoc constant  = 1 kpc to the "correct" size of the emitting disk.Observationally, of course, these are largely unknown, hence taking   = 1 in BK05.But here, we can estimate the true values in the simulations of each parameter.
The value of  can be directly read off from Fig. 5, increasing from ∼ 0.7 − 3 at  ∼ 1 − 10 kpc.This radial trend is expected given the diffusive nature of CRs (see discussion of this in Hopkins et al. 2022a), since it means the CR energy density will fall less-rapidly than the gas pressure (and  B ) at large galacto-centric radii.Taking  B to be the atomic disk scale-height gives  B ∼ 200 pc, similar to canonical scale heights for the star-forming disks of the Milky Way and other observed galaxies (Kalberla & Kerp 2009;Yim et al. 2014;Patra 2020;Gensior et al. 2023).
Direct examination of the  CR − B correlations (shown in Hopkins et al. 2022a) gives  ∼ 0.2 − 0.3: this arises again because the CRs are diffusive, so (especially on small scales) always form a locally-smooth distribution compared to the magnetic fields (Butsky & Quinn 2018;Chan et al. 2019;Buck et al. 2020).And we can estimate   directly from the midplane scatter in ln  B giving   ∼ 2 − 4 -this is naturally expected from the same turbulence models which motivated our lognormal assumption in the first place, which predict   ≈ ln[1 + M 2  ] (where M  is the compressive Mach number of ISM turbulence, and M  ∼ a few in both observations and simulations; see Federrath et al. 2010;Hopkins 2013).
We can immediately follow the same procedure to calculate the intensity-weighted magnetic energy density, ⟨ B ⟩   = (/( + 1)) exp (   ) ⟨  0 ⟩, as well as the volume-weighted and intensity-weighted midplane CR energy densities This provides a natural explanation for several phenomena we saw in Fig. 3.The ratio ⟨ B ⟩   /⟨ B ⟩ V ∼ 30 − 40 in the outer disk, arises primarily from the "clumping factor" correction.The ratio of ⟨ B ⟩ V = ⟨ B 0 ⟩ in the disk mid-plane to ⟨ B ⟩ BK05 (the standard equipartition formula) is given by (  cl   /) −2/(+3) , which in the outer disk is ∼ 0.7 (and is ∼ 1.4 in the inner disk) so the BK05 formula will slightly over-estimate (under-estimate) ⟨ B ⟩ V in the outer (inner) disk, as we see.The difference is small -i.e.BK05 appears to "correctly" obtain roughly the correct ⟨ B ⟩ V , because the corrections  cl and  B / are both relatively large but tend to go in opposite directions (cancelling each other out), combined of course with the effect of the small power law 2/( + 3) which tends to suppress any differences.In other words, the true emitting region is smaller along the line of sight than assumed by BK05, but the emission is also boosted within that region by clumping.Of course, BK05 under-estimates ⟨ B ⟩   by a very large factor as predicted.
Using the same exercise/model we can compare the intensity and volume-weighted CR energy densities ⟨ CR ⟩   and ⟨ CR ⟩ V .Given the small  ≪ 1, we predict that although ⟨ B ⟩   ≫ ⟨  ⟩ V , ⟨ CR ⟩   only exceeds ⟨ CR ⟩ V by a factor ∼ 1.5 − 2. Thus for a given ,  and   or  cl , and  B , we can quantitatively explain all of the relative values of the different estimators in Fig. 3.

SUMMARY AND CONCLUSIONS
In this work, we have presented the first end-to-end predictions of synchrotron emission from MHD galaxy formation simulations which self-consistently evolve B and the CR(e) spectra from Hopkins et al. (2022a).These simulations utilize a constant in space and time scaling for the CR scattering rate , calibrated to reproduce both Solar System CR data (e.g.Voyager, AMS-02) and resolved -ray observations of the MW and nearby galaxies (e.g.Fermi).
We have found that synchrotron emission in L * galaxies arises not only from the volume-filling, warm/hot phases of the ISM, but can be dominated by cooler and denser WNM/CNM phases.This is not unexpected; the long discovered and well studied FIR-Radio correlation of galaxies (Voelk 1985;Ivison et al. 1985;de Jong et al. 1985;Helou et al. 1985;Condon et al. 1991) which exists also on sub-kpc resolved scales (Murphy et al. 2006) requires a connection between the FIR emission, which arises from dust re-radiating UV photons from star formation, and the synchrotron emission which also arises from star forming regions with high neutral gas densities and magnetic field strengths.Furthermore, recent synchrotron observations of edge-on galaxies (Krause et al. 2018;Heesen et al. 2018) have found bright thin disk components which would be spatially coincident with thin, mostly neutral mid-plane gas, and contribute most of the emission when viewed face-on.Moreover, recent observations of structures within our Galaxy have found synchrotron emission arising from cold, neutral gas as well (Bracco et al. 2023).
While this is now known, the conventional wisdom when applying equipartition models to observations of extra-galactic non-thermal radio continuum emission has been to implicitly assume that the emission arises from a volume-filling phase of the ISM which is far more extended in the vertical direction compared to the atomic gas disk.This assumption can lead to underestimating the "true" B of the synchrotron-emitting dense gas, and over-estimating B in the volume-filling, tenuous thick disk/halo gas at ∼kpc above the disk.
We have found that explicitly evolving the CR(e) spectra is important for accurate synchrotron predictions towards galactic centers, where loss terms are drastically different from typical spiral galaxy conditions.Comparing to the "single-bin" scenario shows that the resulting predicted emission changes relative to assuming a constant LISM spectrum by modest factors of ∼1.5-2 in typical spiral galaxy conditions at outer radii (R > 3 kpc), but can be particularly important by a factor ∼10 -50 towards the galactic center, where loss terms can be drastically different.This is consistent with the softer CRe spectrum seen in these simulations towards the galactic center in Hopkins et al. (2022a), and highlights the varying importance of different loss rates and CR source distributions in generating predictive spectra and synchrotron images.We show that the difference in synchrotron owes primarily to variation in the p/e − ratio compounding with changes in spectral slope, rather than e + contributions.
Finally, we formulate a toy model that accounts for clumping factors, a varied magnetic scale height, and deviations from equipartition in order to more robustly trace B weighted by different quantities.From our toy model calculations, we find that uncertainty in connecting the equipartition values to the "true" values of u B and e CR boils down to the assumption of energy equipartition and the size/volumefilling factor of the emitting regions, and less so on spectral effects.When estimating a volume-averaged mean ⟨ B ⟩, the fact that CRs are diffusive and so naturally form a smoother distribution, means that there are large local violations of equipartition.But this smoothness partially cancels the "clumping factor" from in-homogeneous B. This leads to surprisingly reasonable values of B. Future high-resolution radio observations at GHz frequencies (Murphy et al. 2018) may further constrain the volume filling and clumping factors of synchrotron emission, including the current observationally uncertain thin-disk scale heights.
We note that while we have only studied simulations with constant power-law scattering (or diffusivity) of CRs in this work, we have also studied a set of FIRE-2 simulations with varied CR transport motivated by extrinsic turbulence and self-confinement in Ponnada et al. (2023), though with the caveat that those are "single-bin" simulations.There, we find that the CR transport physics can drive differences in gas phase structure, morphology, and non-thermal properties, leading to markedly different synchrotron emission particularly for the self-confinement models relative to extrinsic turbulence or constant diffusivity models.Correspondingly, there are implications for what phase of the ISM dominates the emission and the distribution of gas in the disk and inner CGM, which would affect what one would infer with equipartition assumptions.In particular, in runs with selfconfinement motivated CR transport, the strong trapping of CRs in regions of high e CR can lead to non-linear CR-driven, magnetized winds and thus result in a more volume-filling, diffuse phase of the ISM dominating the emission, more in-line with fiducial equipartition assumptions, though those simulated galaxies differ in detail from observed star-forming L * galaxies, as we discuss in Ponnada et al. (2023).In future work, we will continue to investigate how varied CR transport physics can act to vary synchrotron properties using fully cosmological, spectrally-resolved CR-MHD galaxy simulations.
Synchrotron emission has also been used to estimate physical properties of CRe like their transport length and diffusion coefficient (with an assumed streaming or advection speed), or the CRe scale height (Heesen 2021;Heesen et al. 2023).In future work, we will explore these estimators using our simulations.Azimuthally averaged radial profiles of synchrotron specific intensity normalized by HI surface density for m12i, m12f, and m12m, in navy dashed, solid, and dot-dashed lines.Shaded regions show the approximate ± 1 scatter (32-68 percentile) at a given radial bin.Corresponding radial profiles for nearby spiral galaxies from Basu & Roy (2013) and Beck (2015) are shown in dotted lines, normalized by the Σ  from (Casasola et al. 2017).Normalized in this way, the radial profiles exhibit broader similarities when compared to the specific intensity profiles shown in Figure 3, though the simulations' radial profiles still remain steeper than the observations compared to here.

APPENDIX A: AUXILIARY FIGURES
This paper has been typeset from a T E X/L A T E X file prepared by the author.or ⟨ CR ⟩ V /⟨ B ⟩ V , in cylindrical annuli of varying height from the mid-plane for m12f.For a cylindrical annulus of a given height, in the inner disk where gas densities are high, magnetic energy densities can dominate CR energy densities at the factor of ∼ 2 level for cylindrical heights ≤ 1 kpc in this volume-averaged sense, with this trend increasing as more of the dense midplane gas is sampled (lower heights).In the outer-disk, where gas densities are lower, CRs start to dominate the relative energy density, with the effect accentuated by sampling more "halo" gas at larger heights above the disk where e CR tends to dominate u B by factors of ∼ 2-4 for heights ≤ 1 kpc at the very outer radii (Ponnada et al. 2022).

Figure 1 .
Figure1.Images of specific intensity (I  ) for m12i, m12f, and m12m at  = 6.2 cm, with the colorbar showing log 10 (I  ).All three simulations exhibit a range of I  values broadly consistent with radio continuum observations of nearby spiral galaxies.The synchrotron intensity closely traces the neutral gas spatial distribution, with enhanced emission coincident with spiral structure and relatively lower levels of emission in inter-arm regions.

Figure 2 .
Figure 2. ISM phase diagram, or 2-D histogram of temperature vs. density of gas cells, weighted by contribution to specific synchrotron intensity for m12f at  = 6.2 cm, within cylindrical R < 10 kpc.Color-bar shows the probability density.Consistent with visual inference from Figure 1, the synchrotron emission is largely dominated by neutral gas typical in CNM and WNM ISM conditions, rather than diffuse, volume-filling, and ionized phases.

Figure 3 .
Figure3.Azimuthally averaged (mean) radial profiles of synchrotron specific intensity at 0.33 GHz for m12i, m12f, and m12m, in navy dashed, solid, and dot-dashed lines.Shaded regions show the 25-75 percentile range at a given radial bin.Corresponding radial profiles of non-thermal synchrotron emission (free-free emission corrected) for nearby face-on spiral galaxies fromBasu & Roy (2013) andBeck (2015) are shown in dotted lines, with 4.83 GHz observations of IC342 scaled up by / 0  .Our synthetic synchrotron images for the fiducial model agree to within an order-of-magnitude with spatially resolved synchrotron emission from nearby spiral galaxies.Variation in the shapes of the profiles towards the galactic centers and outskirts partially arises due to differences in the gas surface density (see FigureA1for corresponding plot), though the simulations appear to produce steeper radial profiles in comparison to the observed systems compared.

Figure 5 .
Figure 5. Azimuthally averaged radial profiles of u (left) and e CR (right) weighted by volume and I  in comparison to u Beq for m12f.Lines represent u B (pink), e CR (navy), volume-weighted (solid) and I ,  = 0.33GHz -weighted (dot-dashed) averages of each quantity, within |z| ≤ 0.5 kpc.The result of applying Equation10to the radial profile presented in Figure3is shown, with L = 1, K 0 = 100, f V = 1 (purple dashed) and f V computed within each radial bin (teal dashed).Shaded regions show the approximate ± 1 scatter (32-68 percentile) of the emission-weighted quantity and fiducial equipartition model at a given radial bin.Unilaterally, we see that without correcting for the volume-filling factor of the neutral mid-plane gas which dominates I  , the equipartition formula under-predicts the "true" I  -weighted B by ∼0.3-0.6 dex and e CR by by ∼0.1-0.6 dex.Furthermore, due to the emission being dominated by mid-plane gas, the equipartition model with the right volume-filling factor is generally not representative of the volume-weighted u B , though the fiducial model surprisingly traces the volume-weighted quantities well owing to a confluence of factors (see Section 4).While not shown in this figure, m12i and m12m show generally the same behavior.
Figure A1.Azimuthally averaged radial profiles of synchrotron specific intensity normalized by HI surface density for m12i, m12f, and m12m, in navy dashed, solid, and dot-dashed lines.Shaded regions show the approximate ± 1 scatter (32-68 percentile) at a given radial bin.Corresponding radial profiles for nearby spiral galaxies fromBasu & Roy (2013) andBeck (2015) are shown in dotted lines, normalized by the Σ  from(Casasola et al. 2017).Normalized in this way, the radial profiles exhibit broader similarities when compared to the specific intensity profiles shown in Figure3, though the simulations' radial profiles still remain steeper than the observations compared to here.

Figure A2 .Figure A3 .
Figure A2.Edge-on image analysis: Left: 10" Gaussian beam-convolved image of m12f at 6 GHz with the observer at a distance of 20 Mpc (1" ∼100 pc) with bins overlaid akin to the method of(Krause et al. 2018).Center: Lines showing intensity profiles for the central strip of bins (pink solid) with best-fit two-component exponential fit overlaid (purple dashed).Right: Lines showing vertical profiles of volume-weighted magnetic (pink dot-dashed) and CRe (navy solid) energy densities, with shaded regions showing 25-75 percentile ranges.Our results are consistent with observational constraints estimates of thin and thick disk components of edge-on spiral galaxies, with the thin disc component with small scale height (∼120 pc when averaged across each vertical strip) dominating the emission and the thick disk component contributing little to the emission but extended out to scale heights of ∼kpc.This primarily owes to a small magnetic scale height, with u B falling off more quickly than e CRe , which owing to CRe diffusion havs a flatter profile.
Bisschoff et al. (2019)ean) radial profiles the of I  to I ,   under different assumptions for j e for m12i (top) and m12f (bottom) at 4.83 GHz.Lines show fiducial intensity using internally evolved j e including e + contribution (navy solid), fiducial intensity without e + (light blue dot-dashed), holding spectral shape and K constant by re-scaling j e ofBisschoff et al. (2019)(B19 ) Figure 4.