# Simulating Radio Synchrotron Morphology, Spectra, and Polarization of Cosmic
Ray Driven Galactic Winds

H.-H. Sandy Chiu Department of Astronomy, University of Michigan, Ann Arbor,
MI 48109, USA  Mateusz Ruszkowski Department of Astronomy, University of
Michigan, Ann Arbor, MI 48109, USA  Max-Planck-Institut für Astrophysik (MPA),
Karl-Schwarzschild-Str. 1, 85748 Garching, Germany  E-mail<EMAIL_ADDRESS>(MR) Timon Thomas Leibniz-Institute for Astrophysics Potsdam (AIP), An der
Sternwarte 16, 14482 Potsdam, Germany  Maria Werhahn Max-Planck-Institut für
Astrophysik (MPA), Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
Christoph Pfrommer Leibniz-Institute for Astrophysics Potsdam (AIP), An der
Sternwarte 16, 14482 Potsdam, Germany

###### Abstract

The formation of galaxies is significantly influenced by galactic winds,
possibly driven by cosmic rays due to their long cooling times and better
coupling to plasma compared to radiation. In this study, we compare the radio
observations of the edge-on galaxy NGC 4217 from the CHANG-ES collaboration
catalog with a mock observation of an isolated galaxy based on the arepo
simulation that adopts the state-of-the-art two-moment cosmic ray transport
treatment and multiphase interstellar medium model. We find significant
agreement between the simulated and observed images and spectroscopic data.
Specifically, we find that (i) while the radio flux is sensitive to the
magnitude of the magnetic field, the shapes of intensity profiles are only
weakly dependent on the normalization of the magnetic field and are robustly
predicted by our galaxy simulations; (ii) the multi-wavelength spectrum above
0.1 GHz is in agreement with the radio observations and its slope is also only
weakly sensitive to the magnetic field strength; (iii) the magnetic field
direction exhibits an X-shaped morphology, often seen in edge-on galaxies,
which is consistent with the observations and indicates the presence of a
galactic-scale outflow. Since the normalization of the intensity profiles and
spectrum are degenerate with respect to the magnetic field strength and the
normalization of the electron spectrum, this qualitative match is possible
provided that the simulated magnetic field is only mildly boosted, which could
be attributed to the missing galactic disk interaction with the CGM in our
idealized simulation. Our results highlight the importance of incorporating
advanced cosmic ray transport models in simulations and provide a deeper
understanding of galactic wind dynamics and its impact on galaxy evolution.

cosmic rays – galaxies: magnetic fields – radio continuum: galaxies – methods:
numerical – radiation mechanisms: non-thermal

#1 #1

##  1 Introduction

Most galaxies exhibit baryon fractions significantly below the average
cosmological value. The missing baryons are either not captured by the
developing protogalaxies’ potential wells or are expelled due to feedback
mechanisms during the galaxy formation process (Anderson & Bregman, 2010;
Somerville & Davé, 2015). In halos with masses comparable to that of the Milky
Way, which have the most effective star formation efficiency, only around 20%
of the available baryons are converted into stars (Bregman, 2007; Moster et
al., 2010; Tumlinson et al., 2017). In galaxies exceeding the mass of the
Milky Way, active galactic nuclei are believed to play a significant role in
supplying the necessary feedback to suppress excessive star formation (Croton
et al., 2006), while in less massive galaxies, starburst activities are widely
accepted to be the main mechanism of generating outflows and thus limit the
amount of gas mass available for star formation (Veilleux et al., 2005; Bland-
Hawthorn et al., 2007).
In the traditional framework of galactic winds powered by supernovae
(Chevalier & Clegg, 1985), galactic outflows are initiated by thermal energy,
propelling the gas out of the galaxy ballistically. Although these models have
demonstrated the ability to generate superwinds observed in starburst galaxies
(Bustard et al., 2016), they are still facing obstacles such as the
overcooling problem (Walch et al., 2015; Naab & Ostriker, 2017), where the
thermal energy released by supernovae is rapidly dissipated, leading to a
decrease in the effectiveness of supernova explosions in regulating star
formation. Feedback from supernovae can also be driven by radiation pressure
(Murray et al., 2005, 2011; Hopkins et al., 2012). However, this model also
faces challenges in effectively coupling stellar radiation to the gas (Rosdahl
et al., 2015), resulting in reduced momentum transferred to the gas and
allowing radiation to escape through a path of least resistance. Furthermore,
reprocessing of UV and optical radiation from massive stars is not expected to
have a significant dynamical effect on the surrounding gas due to the
reprocessing of this radiation to far-IR, where the coupling of radiation to
the gas is weak (Reissl et al., 2018; Menon et al., 2022). These processes
reduce the efficiency of stellar feedback and make both the thermal energy-
driven and radiation-pressure-driven models insufficient in explaining
observations (Steidel et al., 2010).
Cosmic rays (CRs) have recently gained significant attention for their role in
galactic feedback (Ruszkowski & Pfrommer, 2023). These relativistic particles
accelerated in supernova remnants (SNRs) (Blandford & Eichler, 1987; Caprioli
& Spitkovsky, 2014) or by winds of massive young stars (Bykov, 2014), are
relatively sparse in the interstellar medium (ISM). However, observations from
galaxies with distinct outflow characteristics such as M82 and NGC 253
(Ackermann et al., 2012; Bolatto et al., 2013; Yoast-Hull et al., 2013)
suggest that the CR energy density is comparable to, or even exceeds, that of
thermal gas and magnetic fields (Acciari et al., 2009; Paglione & Abrahams,
2012). These observations indicate that CRs can drive, shape, and alter large-
scale galactic outflows.
In contrast to traditional feedback models, where thermal energy is injected
at SNR sites and quickly radiated away, injection of a fraction of SN energy
in the form of CRs can result in a prolonged dynamical impact of CRs on the
gas due to relatively long CR cooling times. Thus, CR feedback offers a
possible solution to the “overcooling" problem (Jubelgas et al., 2008).
Following their injection at the SNR sites, CRs can be transported via
diffusion, advection, and streaming to a low-density region above the galactic
disk, where they can accelerate this gas more effectively than radiation due
to better coupling of CRs to the gas compared to radiation (Socrates et al.,
2008; Huang & Davis, 2022). Consequently, CRs can be a critical feedback
agent, the efficiency of which depends sensitively on the properties of CR
transport.
This motivates the need to constrain CR transport processes. Edge-on galaxies
serve as good laboratories for studying CR transport models because in these
objects non-thermal emission generated by CRs can be directly observed above
and below galactic disks. One-dimensional CR transport models (Heesen et al.,
2016, 2018) have been developed to constrain transport processes and quantify
the effective CR diffusion coefficient. However, these models rely on
simplifying assumptions regarding the magnetic field and suffer from
uncertainties associated with modeling of thermal emission. While one-
dimensional models are valuable tools for constraining physical parameters and
understanding individual observations, their reliance on these assumptions
restricts predictive capability of this approach.
Magneto-hydrodynamical simulations of supernova feedback including simplified
prescriptions for CR transport (i.e., spatially-constant diffusion coefficient
or constant streaming velocity) have demonstrated that CR pressure forces can
drive galactic winds (e.g., Breitschwerdt et al., 1991; Uhlig et al., 2012;
Salem & Bryan, 2014; Pakmor et al., 2016). Models that include this physics
can reproduce electron and proton spectra in the local environment observed by
Voyager-I and AMS-02 (Werhahn et al., 2021b), FIR–$\gamma$-ray correlation
(Pfrommer et al., 2017; Werhahn et al., 2021c), and the FIR-radio correlation
(Werhahn et al., 2021a; Pfrommer et al., 2022). However, in a realistic
situation, CR transport is likely to include a combination of diffusion and
streaming depending on the CR scattering frequency with low and high
scattering rates corresponding to diffusion and streaming, respectively.
Simulations that include these effects have been performed (e.g., Ruszkowski
et al., 2017; Wiener et al., 2017), but have not faced the same level of
scrutiny in the context of specific observational predictions and comparisons
to the observational data as those based on the above simplified CR transport
prescriptions. Furthermore, CR transport is generally expected to be highly
spatially variable (e.g., Farber et al., 2018; Hopkins et al., 2021). By
comparing the streaming instability growth rate to the dominant wave damping
rates – ion neutral and non-linear Landau damping rate – one can infer the
relative contributions to CR flux due to diffusion and streaming (Armillotta
et al., 2021, 2022). The resulting effective CR transport speeds are very
large in the predominantly neutral ISM, which imposes severe limitations on
the numerical scheme requiring prohibitively small timesteps for numerical
stability.
The above considerations motivate development of more efficient computational
methods that can tackle more realistic spatially variable and fast CR
transport without sacrificing numerical accuracy, and thus enable more
meaningful comparisons to observations. This challenge can be successfully
addressed by employing two-moment methods, which draw on the analogy between
CR and radiative transport and extend the set of governing equations to
include explicit time evolution of the CR flux rather than assuming that it is
given by its steady-state form (Jiang & Oh, 2018). This treatment can be
further extended to include explicitly evolving the energy density in the
Alfvén waves responsible for scattering CRs (Thomas & Pfrommer, 2019; Thomas
et al., 2023). The two-moment method has been tested against observations of
radio harps near the Galactic center (Heywood et al., 2019). These radio harps
are likely created by CRs emitted by a moving source such as a pulsar wind
nebula or stellar wind and are injected into the ambient magnetic field. These
CRs follow magnetic field lines and encode their transport history.
Comparisons of the observed brightness profiles to those simulated using the
two-moment method, demonstrate that CR transport is best described by a
combination of streaming and diffusion (Thomas et al., 2020).
Given the successful validation of the two-moment method on small scales, this
study aims to determine if this method, when applied to CR transport governed
by the streaming instability, diffusion in the presence of both the non-linear
Landau and ion-neutral damping, and state-of-the-art multiphase ISM model, can
produce galactic outflow properties that agree with the observational data. To
this end, we generate mock observations from a state-of-the-art galactic wind
simulation (Thomas et al., 2024) and compare them with the observational data
from an edge-on galaxy NGC 4217 (Stein et al., 2020). Our work complements
that of Ponnada et al. (2024), who presented comparisons of mock radio
observations based on cosmological simulations of Milky Way-mass galaxies
using a CR magneto-hydrodynamical approach with on-the-fly computation of the
CR spectrum, but differs in several respects. Our two-moment transport method
includes explicit evolution of the energy density in Alfvén waves responsible
for scattering CRs and includes a different approach to modeling the ISM
physics. Our ISM modeling includes a thermochemistry model incorporating
important ionization and low-temperature cooling and heating processes. This
allows us to faithfully evolve the thermodynamic state of the diffuse ISM and
also to model the effects of ion-neutral damping based on atomic collision
data and its subsequent impact on CR transport. We also focus on different
types of observables such as the spatially-resolved profiles of radio
intensity, morphology of the polarization fraction and field orientation
(including Faraday rotation), and radio spectrum.
This paper is structured as follows: In Section 2, we present the simulation
setup and describe how we calculate steady-state CR spectra, radio emission,
and Stokes parameters. In Section 3, we compare mock observations to
observational data focusing on the intensity profiles, radio emission
morphology, polarization, and multi-wavelength spectrum of the galaxy.
Finally, we outline possible caveats and summarize key findings in Section 4.
In Appendix A, we scrutinize the steady-state assumption of our approach of
modeling CR electron spectra and show the Faraday rotation measure of our
galaxy in Appendix B. Finally, we show the vertical synchrotron profiles in
five strips in Appendix C, which are lined up horizontally along the disk and
probe different outflow regions of the galaxy vertically, ranging from the
outskirts to the nuclear region.

##  2 Methods

In this section, we describe the simulation used in this study. This is
followed by the description of our method to calculate via post-processing the
CR spectra, radio emission, and polarized intensity.

###  2.1 Simulation

The simulation used in this study was performed using the moving-mesh code
arepo (Springel, 2010; Pakmor et al., 2016) with its two-moment CR
magnetohydrodynamics module (CRMHD, Thomas & Pfrommer, 2019; Thomas et al.,
2021; Thomas & Pfrommer, 2022) and the crisp ISM physics module (Thomas et
al., in prep., see also Thomas et al. 2024). Below, we first describe our
simulation settings, then provide an overview of the two-moment CR method and
the crisp framework.
We study an isolated Milky Way-type galaxy with a dark matter halo mass of
$M_{\rm 200}=10^{12}~{}\mathrm{M}_{\odot}$. The simulation starts from an
idealized setup that includes a doubly exponential disk for both the gaseous
and stellar components. The gaseous disk has a radial scale length of 5 kpc, a
scale height of 0.5 kpc, and a total mass of $8\times
10^{9}~{}\mathrm{M}_{\odot}$. The stellar disk has the same scale length and
height, but a total stellar mass of $3.2\times 10^{10}~{}\mathrm{M}_{\odot}$.
Star particles are initialized with individual particle masses of
$M_{*}=1000~{}\mathrm{M}_{\odot}$ placed randomly following the exponential
disk profile. Gas particles are initialized in the same way and are
(de-)refined based on the more stringent of two criteria: (i) the mass of the
gas particle should differ by a factor of at most 2 (at least 0.5) from
$1000~{}\mathrm{M}_{\odot}$ and (ii) the volume of the gas particle should
differ by a factor of at most 2 (at least 0.5) from $V=(4\pi/3)r^{3}$ where
the effective resolution is better than or equal to $r=100$ pc in our region
of interest. The latter criterion ensures that our numerical resolution is
approximately an order of magnitude finer than the beam size of the
observation.
Star formation is modeled following a Schmidt-type approach that converts gas
in cells whose number density is above the threshold of 100 hydrogen particles
per cm3 into new stars within a free-fall time, which is 4.5 Myr at this
threshold. These newborn stars participate in the stellar feedback process by
ejecting mass, metals, and energy into the ISM, with respective rates based on
STARBURST99 calculations (Leitherer et al., 1999). We assume that each
supernova injects $1.06\times 10^{51}$ erg of mechanical energy and that 5% of
the mechanical energy supplied by SNe and stellar winds is injected as CRs.
The CR energy density is then evolved using the two-moment CRMHD model of
Thomas & Pfrommer (2019). The initial magnetic field consists of a toroidal
magnetic field with a strength of $B_{\rm tor}=10^{-1}\sqrt{\rho/\rho_{\rm
max}}~{}\mu{\rm G}$, where $\rho$ is the gas density and $\rho_{\rm max}$ is
the maximum gas density at the galaxy center, and a constant vertical magnetic
field with a strength of $10^{-3}~{}\mu$G. The final magnetic field is
unaffected by the initial conditions as the field strength is significantly
amplified by the small-scale dynamo process until the system reaches
equilibrium.
Our isolated galaxy simulation aims to test the two-moment CRMHD method and
the multi-phase ISM model crisp. The two-moment CRMHD method (Thomas &
Pfrommer, 2019) is used to evolve the CR energy density and CR flux density.
The CR transport velocity is determined by how the CRs interact with small-
scale Alfvén waves, which can be damped by non-linear Landau effect (Miller,
1991) or ion-neutral damping (Kulsrud & Pearce, 1969). In addition to the
energy loss due to the interactions of CRs with Alfvén waves, CRs can also
lose energy through hadronic and Coulomb interactions (Pfrommer et al., 2017).
The CR energy density is evolved assuming the gray approach (i.e., we do not
explicitly track CR energy density corresponding to different momenta of CR
particles) and the CR energy loss rates are averaged over the spectrum of CR
momentum.
We use the crisp module (Thomas et al., 2024) to incorporate the multi-phase
physics of the ISM. crisp tracks the ionization stages of 12 different
species, including $\rm H_{2}$, H I, H II, all ionization stages of He, and
the first two ionization stages of C, O, and Si. The package evolves the
ionization stages considering heating from the far ultraviolet (FUV) photons
emitted by newly born stars and various cooling processes including Ly$\alpha$
cooling by H, rotation-vibrational lines from $\rm H_{2}$ (Moseley et al.,
2021), bremsstrahlung cooling at high temperatures (Cen, 1992), and metal line
cooling calculated directly from collision rates (Abrahamsson et al., 2007;
Grassi et al., 2014) at low temperatures ($T<10^{4}$ K) or by interpolating
tables pre-compiled using the chianti code (Dere et al., 1997) at high
temperatures ($T>10^{4}$ K).

###  2.2 Modeling Cosmic-Ray Spectra

We followed the procedure described in Werhahn et al. (2021b) to model CR
spectra using the Crayon+ code. As this approach adopted a spatially constant
diffusion coefficient, here we outline the key steps in the calculation and
discuss the modifications needed to accommodate spatially varying CR transport
due to a combination of diffusion and streaming in the two-moment model.
In order to compute the energy distribution of CR particles, we consider a
cell-based steady-state approximation, which assumes that the timescale over
which the CR energy density changes is longer than the proton cooling
timescale, which is generally longer than the electron cooling timescale. This
assumption is validated in both Werhahn et al. (2021b) and Appendix A. We then
solve the diffusion-loss equation separately in each cell for protons, primary
electrons, and secondary electrons. Under the steady-state assumption, the
source and loss term in the diffusion-loss equation are in balance such that
the total energy of the CR population in each cell remains constant in time.
We model the source term $q_{i}=q_{i}(p_{i})$ as a power-law function of
momentum with an exponential cutoff:

$q_{i}(p_{i})\mathrm{d}p_{i}=C_{i}p_{i}^{-\alpha_{\rm
inj}}\exp{[-(p_{i}/p_{\rm cut,i})^{n}]}\mathrm{d}p_{i}$ (1)

where $p_{i}=P_{i}/(m_{i}c)=\sqrt{[E_{i}/(m_{i}c^{2})]^{2}-1}$ is the
normalized particle momentum and $E_{i}$ is the CR total particle energy,
subscript $i$ specifies CR species ($={\rm p,e}$ for protons and electrons,
respectively), $m_{i}$ is the proton/electron rest mass, and $c$ is the speed
of light. The normalization of the source function $C_{i}$ is determined
through subsequent re-normalization processes (see below). We adopt the
injection spectral index $\alpha_{\rm inj}=2.1$ for both protons and primary
electrons, cut-off momentum of protons $p_{\rm cut,p}=1\ {\rm PeV}/m_{\rm
p}c^{2}$ (Gaisser, 1990) and electrons $p_{\rm cut,e}=20\ {\rm TeV}/m_{\rm
e}c^{2}$ (Vink, 2011), and $n=1$ for protons and $n=2$ for primary electrons
(Zirakashvili & Aharonian, 2007; Blasi, 2010). The source function of
secondary electrons and positrons is calculated from the steady-state proton
spectrum and the parametrization of the differential cross section of
secondary particle production from Kelner et al. (2006) at large proton
energies ($T_{\rm p}>100$ GeV). For $T_{\rm p}<10$ GeV, we use the
parametrization from Yang et al. (2018) and interpolate between the two in the
intermediate range (see equations B1, B5 and B6 in Werhahn et al., 2021b).
The loss term in the diffusion-loss equation consists of cooling and escape
losses. For the cooling losses, we adopt the same approach as that described
in Werhahn et al. (2021b). Specifically, we include hadronic and Coulomb
interactions for protons, as well as synchrotron and inverse Compton (IC)
losses for primary electrons (as detailed in sections 3.1.2–3.1.3 in Werhahn
et al., 2021b). In contrast to Werhahn et al. (2021b) where only advection and
diffusion escape losses were included, we include escape losses due to
advection, diffusion, and streaming losses. In Werhahn et al. (2021b), the
diffusion timescale is calculated by $\tau_{\rm diff}=L_{\rm CR}^{2}/D(E)$,
where $E$ is total energy of the CR spectral bin, $L_{\rm CR}=e_{\rm
CR}/|\nabla e_{\rm CR}|$ is the diffusion length in each cell, $e_{\rm CR}$ is
the energy density of the CR fluid, and $D(E)=D_{0}\ (E/E_{0})^{0.3}$ is the
energy-dependent diffusion coefficient, where $D_{0}=10^{28}~{}{\rm
cm^{2}s^{-1}}$ and $E_{0}=3~{}{\rm GeV}$ are free parameters. This particular
choice of a typical CR energy $E_{0}$ is motivated by the fact that CRs that
contribute most to the pressure of the CR fluid have energies close to the
peak in the CR spectrum that occurs at that characteristic energy (e.g.,
Ruszkowski & Pfrommer, 2023). Unlike in the simulations in Werhahn et al.
(2021b), CR energy flux $f_{\rm CR}$ and the streaming velocity are explicitly
modeled in this work. This allows us to define an effective diffusion
coefficient as $D_{\rm 0,eff}=|\varv_{\rm CR}|L_{\rm CR}$, where the velocity
of the CR fluid is given as

$\varv_{\rm CR}=\frac{f_{\rm CR}}{e_{\rm CR}+P_{\rm CR}}=\frac{3f_{\rm
CR}}{4e_{\rm CR}}$ (2)

where $P_{\rm CR}=(\gamma_{\rm CR}-1)e_{\rm CR}$ is the CR pressure and
$\gamma_{\rm CR}=4/3$ (Thomas et al., 2023). With such-defined effective
diffusion coefficient, the diffusion timescale (that indirectly includes
streaming) becomes $\tau_{\rm diff}=L_{\rm CR}^{2}/D_{\rm eff}(E)=[L_{\rm
CR}/|\varv_{\rm CR}|](E/E_{0})^{-0.3}$. The advection timescale is calculated
by $\tau_{\rm adv}=L_{\rm CR}/\varv_{z}$, where $\varv_{z}$ is the velocity in
the $z$-direction. In order to account for CR losses, we only include in the
calculation the computational cells that include gas leaving the disk, i.e.,
whose for which $\varv_{z}>0$ for $z>0$ and $\varv_{z}<0$ for $z<0$ (see
Werhahn et al., 2021b, for a justification of the neglect of velocities in the
disk).

After solving the diffusion-loss equation, we renormalize the steady-state
spectra of CR protons such that the total kinetic energy density of the CR
population $\int^{\infty}_{0}T_{\rm p}(p_{\rm p})f_{\rm p,0}(p_{\rm
p})\mathrm{d}p_{\rm p}$ equals the CR energy density in our CRMHD model in
each cell, where $T_{\rm p}(p_{\rm p})/(m_{\rm p}c^{2})=\sqrt{p_{\rm
p}^{2}+1}-1$ is the kinetic energy of each momentum bin and $f_{\rm p,0}$ is
the CR proton spectrum. In order to obtain the normalization of the CR
electron spectrum, we assume a constant electron-to-proton injection number
ratio $K^{\rm inj}_{\rm ep}\approx 0.037$ at 10 GeV. Specifically, $K^{\rm
inj}_{\rm ep}$ is defined as

$K^{\rm inj}_{\rm ep}(T_{\rm norm})=\frac{q^{\rm prim}_{\rm e}(p_{\rm
norm,e})}{q_{\rm p}(p_{\rm norm,p})}\frac{\mathrm{d}p_{\rm
e}}{\mathrm{d}p_{\rm p}},$ (3)

where $T_{\rm norm}=10$ GeV is a fixed normalization kinetic energy for both
electrons and protons, $q_{\rm p}$ and $q^{\rm prim}_{\rm e}$ are the
injection spectrum for protons and primary electrons as defined in Equation
(1), and $p_{\rm norm,e/p}=\sqrt{[T_{\rm norm}/(m_{\rm e/p}c^{2})+1]^{2}-1}$.
This shows that $K^{\rm inj}_{\rm ep}$ is independent of the kinetic energy of
CRs as long as $T_{\rm norm}\gg m_{\rm e/p}c^{2}$ and the exponential cutoff
term can be neglected in Equation (3), i.e., $p_{\rm norm,i}\ll p_{\rm
cut,i}$. The former criterion is approximately fulfilled for our chosen
parameters and the associated small correction is discussed below (see the end
of this section), while the latter criterion is very accurately met for our
chosen parameters. This value of $K^{\rm inj}_{\rm ep}$ can be further related
to the electron-to-proton energy injection fraction

$\zeta_{\rm prim}=\frac{\epsilon^{\rm inj}_{\rm e}}{\epsilon^{\rm inj}_{\rm
p}}=\frac{\int^{\infty}_{0}T_{\rm e}q_{\rm e}\mathrm{d}p_{\rm
e}}{\int^{\infty}_{0}T_{\rm p}q_{\rm p}\mathrm{d}p_{\rm p}}$ (4)

by

$K^{\rm inj}_{\rm ep}=\zeta_{\rm prim}\left(\frac{m_{\rm p}}{m_{\rm
e}}\right)^{2-\alpha_{\rm p}}$ (5)

as long as the cut-off term can also be neglected in the calculation of
$\epsilon^{\rm inj}_{\rm p/e}$. We find that this criterion is fulfilled as
the integral values are very similar (10 percent difference for $\alpha_{\rm
p}=2.2$ and 25 percent difference for $\alpha_{\rm p}=2.1$) regardless of
whether the exponential terms are included. Here, $\alpha_{\rm p}=2.1$ is the
injection spectrum slope (assumed to be the same for electrons and protons)
and $\epsilon^{\rm inj}_{\rm e/p}$ is the energy injected into electrons or
protons (Werhahn et al., 2021b).
While the value of $(m_{\rm p}/m_{\rm e})^{2-\alpha_{\rm p}}$ is relatively
certain, the value of $\zeta_{\rm prim}$ may fluctuate due to different CR
acceleration efficiency for CR protons and electrons. One-zone model
simulations (Lacki et al., 2010) reveal that $2$ percent of the kinetic energy
of supernova explosion should go to primary CR electrons and $\sim$10-20
percent to CR protons to fit the FIR-radio correlation, which gives
$\zeta_{\rm prim}=0.1-0.2$ and corresponds to $K^{\rm inj}_{\rm
ep}=0.024-0.048$ for $\alpha_{\rm p}=2.2$ and $K^{\rm inj}_{\rm
ep}=0.059-0.12$ for $\alpha_{\rm p}=2.1$. On the other hand, cosmological
simulations (Ponnada et al., 2024) adopt even smaller values as $\zeta_{\rm
prim}=0.02$, which corresponds to $K^{\rm inj}_{\rm ep}=0.0048$ for
$\alpha_{\rm p}=2.2$ and $K^{\rm inj}_{\rm ep}=0.012$ for $\alpha_{\rm
p}=2.1$. In conclusion, the values of $\zeta_{\rm prim}$ can vary from 0.02 up
to 0.2, and $K^{\rm inj}_{\rm ep}$ that we adopt corresponds to $\zeta_{\rm
prim}=0.03$ and falls within the reasonable range of $\zeta_{\rm prim}$.
In practice, we normalize the injection spectrum for protons and electrons at
their momentum bins that are closest to a kinetic energy of 10 GeV, which is
9.63 and 9.08 GeV for protons and electrons, respectively, due to the limited
resolution when modeling the spectrum in momentum space. This small offset
leads to a correction factor of $(T_{\rm e,code}/T_{\rm p,code})^{-\alpha_{\rm
p}}(1+\alpha_{\rm p}m_{\rm p}c^{2}/T_{\rm p,code})$, where the second term in
brackets comes from a Taylor expansion of the normalized CR momentum with
$m_{\rm p}c^{2}/T_{\rm p,code}$ serving as a smallness parameter. Therefore,
$K^{\rm inj}_{\rm ep}$ used in the code is $K^{\rm inj}_{\rm ep,code}=K^{\rm
inj}_{\rm ep}(9.08/9.63)^{-\alpha_{\rm p}}[1+\alpha_{\rm p}(0.94/9.63)]=0.05$.

###  2.3 Radio Emission

Radio emission comprises of non-thermal synchrotron emission and thermal free-
free emission. Here we outline how we model these radiative processes to
understand their contributions to the overall radio emission.

####  2.3.1 Non-thermal Synchrotron Emission

We compute the synchrotron emission following Rybicki & Lightman (1986) as

$j_{\nu}=\frac{\sqrt{3}e^{3}B_{\perp}}{m_{\rm e}c^{2}}\int^{\infty}_{0}f_{\rm
e}(p_{\rm e})F(\nu/\nu_{c})\mathrm{d}p_{\rm e},$ (6)

where $e$ is the elementary charge, $B_{\perp}$ is the magnetic field
perpendicular to the line of sight, $f_{\rm e}=\mathrm{d}N_{\rm
e}/(\mathrm{d}E_{\rm e}\mathrm{d}V)$ is the CR electron spectral density
(i.e., number of particles per unit volume and unit energy), $\nu$ is the
observed frequency, and $F(x)\equiv
x\int_{x}^{\infty}K_{5/3}(\zeta)\mathrm{d}\zeta$, where $K_{5/3}(\zeta)$ is
the modified Bessel function of order 5/3. We adopted an analytical
approximation to $F(\nu/\nu_{\rm c})$ provided by Aharonian et al. (2010). The
role of the function $F(\nu/\nu_{\rm c})$ is a synchrotron emission kernel of
an individual electron that is centered on the typical synchrotron emission
frequency $\nu_{\rm syn}\approx 2\nu_{\rm c}$ for a given energy of a CR
particle, where the critical frequency $\nu_{\rm c}$ is defined as

$\nu_{c}=\frac{3eB_{\perp}}{4\pi m_{\rm e}c}\gamma_{\rm e}^{2},$ (7)

where $\gamma_{\rm e}=\sqrt{p_{\rm e}^{2}+1}$ is the Lorentz factor
corresponding to non-dimensional CR electron energy. This implies that
observations of synchrotron emission at a fixed frequency corresponding to
lower magnetic field strength probe electrons with higher energy and thus
lower CR electron number density.

####  2.3.2 Thermal Free-Free Emission

In addition to synchrotron emission, we also expect thermal bremsstrahlung to
contribute to the radio band. We follow Rybicki & Lightman (1986) and
calculate the free-free emissivity from

$j_{\nu,\rm ff}=6.8\times 10^{-38}Z^{2}\frac{n_{\rm e}n_{\rm i}}{\rm
cm^{-6}}\left(\frac{T}{{\rm K}}\right)^{-0.5}e^{-h\nu/(k_{\rm B}T)}{\rm
ergs^{-1}cm^{-3}Hz^{-1}}$ (8)

where $Z$ is the atomic number, $n_{\rm e},n_{\rm i}$ is the thermal electron
and ion density respectively, and $T$ is the thermal gas temperature. Unlike
in Werhahn et al. (2021a), who assume constant gas temperature of $8000$ K, we
adopt the thermal gas temperature that is self-consistently calculated in the
CRMHD simulation.

####  2.3.3 Generating Mock Observation

In addition to emission, we also consider free-free absorption and synchrotron
self-absorption (see Appendix A in Werhahn et al. (2021a)) when calculating
the final emission maps. We calculate the mock observation with the following
procedure: (i) rotate the whole simulation domain to the desired inclination,
(ii) calculate the magnetic field perpendicular to line of sight and the
resulting emissivity for non-thermal and thermal emission, (iii) divide the
simulation cube into thin slices perpendicular to the line of sight, and
calculate absorption coefficients for each slice, and finally (iv) calculate
the final emission modified by absorption using the radiative transfer
equation:

$4\pi I_{\nu}=\int_{\rm
s_{0}}^{0}j_{\nu}(s)\exp{\left[-\tau_{\nu}(s,0)\right]}\mathrm{d}s,$ (9)

where $s$ is the distance to the observer along the line of sight, $s_{0}$ is
the distance to the far end of the computational domain,
$j_{\nu}=j_{\rm\nu,syn}+j_{\rm\nu,ff}$ is the total emissivity,
$\tau_{\nu}(s,0)=\int_{s}^{0}(\kappa_{\rm\nu,syn}+\kappa_{\rm\nu,ff})\mathrm{d}s$
is the optical depth, and $\kappa_{\rm syn/ff}$ is the absorption coefficient
for synchrotron and free-free emission. We include synchrotron emission and
absorption from both primary and secondary CR electrons.
Finally, in order to account for the finite angular resolution of observations
and to compare our results to the data for NGC 4217 (Stein et al., 2020), we
convolve the resulting emission with the imsmooth function in the CASA package
(McMullin et al., 2007) adopting the beam of 13.4 arcsec in diameter (half
power beam width), which corresponds to $1.3~{}{\rm kpc}$ appropriate for this
galaxy at a distance of $20.3~{}{\rm Mpc}$.

###  2.4 Calculation of Stokes Parameters

In order to compare the polarized intensity and the observed magnetic field
orientation with the observation, we need to calculate the Stokes $Q$ and $U$
parameters including the Faraday rotation. In this section, we outline the
procedure for this calculation, which closely follows that in Waelkens et al.
(2009).
First, we calculate the polarized emissivity in each cell using the formula in
Rybicki & Lightman (1986), which is similar to Equation (6) for the total
emissivity except that it contains a different synchrotron kernel:

$j_{\nu,\rm pol}=\frac{\sqrt{3}e^{3}B_{\perp}}{m_{\rm
e}c^{2}}\int^{\infty}_{0}f_{\rm e}(p_{\rm e})G(\nu/\nu_{c})\mathrm{d}p_{\rm
e},$ (10)

where $G(x)=xK_{2/3}(x)$ is the dimensionless synchrotron kernel function and
$K_{2/3}$ is the modified Bessel function of order $2/3$. Second, we calculate
the Stokes $Q$ and $U$ emissivities from

$\displaystyle j_{Q}$ $\displaystyle=j_{\nu,\rm pol}\cos{(2\chi)},$ (11)
$\displaystyle j_{U}$ $\displaystyle=j_{\nu,\rm pol}\sin{(-2\chi)},$ (12)

where $\chi$ is the polarization angle of each computational cell given by

$\chi={\rm RM}\lambda^{2}+\chi_{0},$ (13)

where $\chi_{0}$ is the magnetic field direction in each computational cell,
$\lambda$ is the wavelength corresponding to the measured frequency $\nu$, and
the rotation measure RM is defined as

${\rm RM}=\frac{e^{3}}{2\pi m_{\rm e}^{2}c^{4}}\int_{s_{0}}^{0}n_{\rm
e}B_{\parallel}\mathrm{d}s,$ (14)

where $B_{\parallel}$ is the signed magnetic field parallel to the line of
sight and $s_{0}$ is the distance to the observer. Finally, we integrate the
emissivities to obtain Stokes parameters from

$\displaystyle Q$ $\displaystyle=\int_{s_{0}}^{0}j_{Q}\mathrm{d}s$ (15)
$\displaystyle U$ $\displaystyle=\int_{s_{0}}^{0}j_{U}\mathrm{d}s.$ (16)

The polarized intensity is given by $\sqrt{Q^{2}+U^{2}}$ and the observed
polarization angle $\phi$, which is the orientation of the magnetic field
projected onto the plane of the sky, is defined as

$\phi=\frac{1}{2}\arctan{\frac{U}{Q}}.$ (17)

Figure 1:  Edge-on view of the slice plot (line of sight thickness of 0.5 kpc)
of magnetic field magnitude (upper left), CR electron energy density at
$10~{}{\rm GeV}$ (upper right), synchrotron emissivity (lower left) and line-
of-sight integrated synchrotron intensity (lower right, projection thickness
of 44 kpc). The lower right panel assumes that the galaxy is located at the
same distance as NGC 4217 (20.3 Mpc). This plot shows how the magnetic field
and the CR electron spectrum contribute to the synchrotron emissivity and the
projected intensity, both calculated at 1.6 GHz (L-band) and including
contributions from both primary and secondary CR electrons. Figure 2: Stacked
profile comparison at 1.6 GHz (L band, left panel) and 6 GHz (C band, right
panel). Shown are observational data points (grey points with error bars)
along with profiles corresponding to the (a.) original simulation, (b.)
simulation normalized to the observed flux at the corresponding frequency,
respectively, and (c.) simulation with magnetic field boosted by a factor of
four, which is consistent with the observational data points (see legend
embedded in the figure). Although the magnitude of emission depends on the
magnetic field, the shape of the profile is relatively insensitive to it. We
adopt the same size of beam area as in Stein et al. (2020), who use a beam
size of 203.5 arcsec2.

##  3 Results

###  3.1 Radio Intensity Profiles

We first show in Fig. 1 how the magnetic field and CR electron spectrum
contribute to the synchrotron emissivity. Overall, the magnetic field strength
(upper left panel) has similar morphology to that of the CR electron spectrum
at 10 GeV (upper right panel). In particular, a bubble-like feature near
$x=-2$ kpc is seen in the distribution of both quantities. However, the
magnetic field exhibits more disk-like morphology compared to the spatial
distribution of the CR electron spectrum at 10 GeV. These features are
reflected in the lower-left panel, where we show a slice plot of synchrotron
emissivity. The dim bubble-like feature near $x=-2$ kpc corresponds to weak
magnetic fields and low normalization of the CR electron spectrum, while the
bright disk is mainly associated with the strong magnetic field strength.
Finally, we integrate the synchrotron emissivity over 44 kpc along the line of
sight (i.e., over the whole galactic disk) and present the resulting intensity
in the lower right panel. The bubble-like feature in the slice plot is covered
by synchrotron emission from the outflow at different depths along the line of
sight and the intensity distribution is relatively smooth and decreasing as a
function of the distance from the galactic midplane.
In order to quantitatively compare outflow properties to observations, we
adopt an approach similar to that in Stein et al. (2020) and compute the total
intensity profile as a function of height in vertical strips centered on the
disk. Specifically, five strips are set to be perpendicular to the disk with
20 boxes in each vertical strip (see Fig. 9). The size of the box is set to be
the same as in Stein et al. (2020), which has a height of 5 arcsec and width
of 35 arcsec, corresponding to 0.5 kpc and 3.4 kpc, respectively. The
intensity profile as a function of height $z$ is measured by taking the
average intensity within each box. The resulting intensity profiles are shown
in Appendix C (see Fig. 10 and Fig. 11). Finally, we construct the overall
profile for all strips together and plot the profile.
We show the mean intensity (solid blue line) as a function of height along
with the observational data points (gray points with error bars) in Fig. 2.
Both the simulated and observed profiles have a peak at $z=0$ where the
galactic disk is located and decrease by $\sim$1.5 orders of magnitude towards
higher latitudes ($\sim$5 kpc). However, the simulated profile systematically
underpredicts the observations by about 1 order of magnitude. We find that the
simulated profile can match the observations within the error bars when the
magnetic field is boosted by a factor of four and the CR spectrum is re-
calculated accordingly (especially at 1.6 GHz). Such-obtained profile is shown
as a solid red line. We suggest that the underestimation of the magnetic field
in the simulation may be due to the lack of a realistic CGM. In isolated
galaxy simulations, the CGM is completely provided by the galactic wind,
whereas in cosmological simulations, the CGM originates both from the gas
accreted from large distances and the gas expelled by galactic winds in the
past across cosmic time. The interactions between the winds and the accreting
gas can result in amplification of the magnetic field via a small-scale
fluctuating dynamo. This process could be further affected by merger events at
higher redshift. Both the accretion processes and merger events can drive
turbulence and therefore amplify the magnetic fields (Pakmor et al., 2014,
2017; Martin-Alvarez et al., 2018). For example, Ponnada et al. (2024) also
study Milky Way–mass galaxies but in the context of cosmological zoom-in
simulations. They obtain a magnetic field magnitude of $\sim$10 $\mu$G at the
disk and $\sim$1 $\mu$G at $z=4$ kpc, which is a factor of 2 to 3 larger than
our magnetic field magnitude.
We note that the normalization of the intensity profiles (and the radio
spectrum) depends not only on the strength of the magnetic field but also on
the normalization of the primary electron spectrum, which is a free parameter
in our model. For the assumed magnetic field boost, the model fits to the data
are consistent with $\zeta_{\rm prim}\sim 0.03$, while $\zeta_{\rm prim}$ can
be as large as 0.2 (Lacki et al., 2010). Since (i) the intensity profile
shapes are only very weakly sensitive to the B-field boost and (ii) the
profile normalization depends on both the boost and $\zeta_{\rm prim}$, the
intensity profile normalization is degenerate with respect to the these two
parameters. Specifically, reasonable model fits to the intensity profiles (and
the radio spectrum) may be obtained by adopting a combination of these two
parameters, and may require smaller B-field boost for higher $\zeta_{\rm
prim}$, but detailed fitting of the models to the data is beyond the scope of
this work.
The impact of the underestimation of the magnetic field on radio intensity can
be significant. In steady state, the dependence of the radio emissivity on the
magnetic and radiation field is (see, e.g., Pfrommer et al., 2008; Ruszkowski
& Pfrommer, 2023)

$j_{\nu}\propto\frac{B^{\alpha_{\nu}+1}}{B_{\rm ph}^{2}+B^{2}},$ (18)

where $j_{\nu}$ is the synchrotron emissivity, $B$ is the magnetic field
strength, and $B_{\rm ph}^{2}=8\pi u_{\rm ph}$, where $u_{\rm ph}$ is the
energy density in the photon field with contributions from the interstellar
radiation field and the CMB. The observed radio spectral index $\alpha_{\nu}$
is related to the injected electron spectral index $\alpha_{\rm inj}$ as
$\alpha_{\nu}=\alpha_{\rm inj}/2=1.05$, yielding $j_{\nu}\propto B^{2.05}$ if
$B\ll B_{\rm ph}$. Before boosting, $B_{\rm ph}$ is larger than $B$. As we
increase the magnetic field, it becomes comparable to and smaller than $B$ at
$z\sim 4$ kpc and $z\sim 1$ kpc, respectively. Therefore, the dependence of
$j_{\nu}$ on $B$ becomes weaker than $j_{\nu}\propto B^{2.05}$. This explains
why the increase in the radio intensity near the disk midplane is slightly
weaker than in the wings of the profile as the field strength is increased. In
practice, we find that the dependence of the radio emission on the magnetic
field is still super-linear and hence the radio intensity can be boosted
significantly by increasing the magnetic field.

Figure 3: C-band (6 GHz) emission map based on the simulation (left panel) and
observation (right panel). The simulated emission map looks similar to the
observed image. Both images are smoothed with a beam size of 13.6 arcsec in
diameter, whose relative size is shown in the lower-left corner in both
panels. The simulated emission map is calculated with a boosted magnetic
field.

Although the magnitude of the intensity profile depends on the magnetic field,
the shape of the profile is relatively insensitive to it. For reference, in
Fig. 2 we show as a dashed blue line the simulated profile rescaled to match
the integrated intensity of the observed profile. This rescaled profile
closely follows the observed one (at 6 GHz) and is similar to the profile
computed by boosting the magnetic field, which demonstrates that changing the
magnetic field strength results in only a slight difference in the shape of
the intensity profile. Due to the relative sensitivity of the radio intensity
to the normalization of the magnetic field, it is challenging for the
simulated profile computed without any adjustments to the magnetic field
strength (or, alternatively, normalization of the electron spectrum) to
exactly match the observed level of emission. However, the emission profile
shapes match the observations very well and are a robust feature of our model.
Uncertainty in the electron-to-proton injection number ratio $K^{\rm inj}_{\rm
ep}$ only affects the normalization rather than the profile shape. Since
increasing $K_{\rm ep}^{\rm inj}$ leads to a higher normalization of the
spectrum of primary CR electrons, it only changes the synchrotron emission
from primary CR electrons. As the secondary synchrotron emission and free-free
emission are at least three times and three orders of magnitude smaller than
the primary synchrotron emission, respectively, the total simulated emission
is dominated by primary synchrotron emission. This implies that increasing the
value of $K_{\rm ep}^{\rm inj}$ will not change the shape of the total radio
emission profile but only increase its normalization.
The prediction of dominant synchrotron emission is consistent with the non-
thermal fraction estimated from observations (Stein et al., 2020), but the
exact synchrotron-to-free free emission ratio depends on detailed modeling of
the ISM properties. While the simulated intensity profile with boosted
magnetic field is generally consistent with the observational data points
within their error bars, the profiles slightly overestimate the intensity near
the profile peak at 1.6 GHz and the intensity in the wings at 6 GHz. There
could be several reasons for these deviations. If CR acceleration efficiency
is reduced, this would result in a decrease in the CR pressure gradient, thus
weakening the galactic outflow and making the intensity profile more
concentrated toward the midplane. This might solve the discrepancy at 6 GHz.
At the same time, this process would also steepen the profiles at 1.6 GHz
resulting in a worse match to observations. However, free-free absorption in
the disk can help to reduce the peak at $z=0$ kpc. Although the current
optical depth due to free-free absorption is $\lesssim 1$ and does not play an
important role in determining the final intensity, increasing the resolution
of the simulation at the disk may decrease the gas temperature (Hummels et
al., 2019; Nelson et al., 2020; Esmerian et al., 2021) and thus lead to an
increase of free-free absorption. Since the free-free absorption coefficient
$\propto\nu^{-2}$, this effect will be more important at lower frequencies
(and at low heights above the disk) and might thus help to cure the small
deviation from observations at 1.6 GHz near the peak of the profile while not
significantly affecting the emission at 6 GHz. Therefore, the exact shape of
the profiles may depend on various details of the model.

###  3.2 Morphology of Radio Emission

Figure 4:  Magnetic field direction calculated from Stokes Q and U maps
overlaid on Stokes I map (left and middle panel) and magnetic field direction
overlaid on the map of polarization degree (right panel). Faraday rotation is
included in the middle and right panel when calculating the Stokes Q and U
maps. The B field direction shows an X-shape pattern at $|z|>1$ kpc,
indicating the existence of an outflow. The polarization degree tends to be
higher in regions where magnetic fields are more ordered and vice versa.

Fig. 3 shows the total radio emission maps from both observations and
simulation at the frequency of 6 GHz with a boosted magnetic field. The right
panel shows the radio observation from Stein et al. (2020) but rotated to
align the disk with the horizontal direction for easier comparison to the
simulated mock image shown in the left panel. Both images are smoothed using
the same beam size (13.4 arcsec in diameter corresponding to $\sim$1.3 kpc at
the distance of interest). Observational noise is added to the simulation by
inverse-Fourier transforming the noise spectrum with different random phases
but the same amplitude as the observational background noise, which is
calculated using the observed intensity distribution at $|z|>6$. After
accounting for these observational effects, the simulated image appears
blurred to the same extent as the observed one and also exhibits some negative
flux regions at the corner of the image.
The simulated emission map displays a wind-like diffuse emission structure
that is distributed across the disk midplane, and closely resembles the radio
emission pattern observed in real galaxies. This is in contrast to previous
works (Werhahn et al., 2021a; Thomas et al., 2023), which exhibit a prominent
outflow feature in the central region of the disk at $x=0$. The main
difference between our simulation and the above earlier work lies in the
initial conditions and the adoption of the crisp ISM model in our work. First,
previous studies (Werhahn et al., 2021a; Thomas et al., 2023) start with a
gravitationally collapsing gas cloud embedded in dark matter halo and simulate
the disk formation process without pre-seeded star particles. Instead, we
post-process the simulations of Thomas et al. (2024) that start from an
idealized setup with star particles seeded randomly following the exponential
disk profile. This idealized disk setup allows the simulation to relax in a
shorter time and can be run at two orders of magnitude higher mass resolution
than the simulation shown by Thomas et al. (2023). Therefore, our simulation
represents a later evolutionary stage of the galaxy. Second, compared to
earlier studies that rely on the Springel & Hernquist (2003) ISM model (but
without adopting the phenomenological wind model), our implementation of the
much more realistic crisp ISM model enables gas cooling to temperatures as low
as several hundred Kelvin, leading to increased star formation activity and
preventing the formation of outflows that are too concentrated toward the
galactic center. As a result, the edge-on morphology of the radio intensity
closely resembles real observations.

###  3.3 Morphology of the Polarized Emission

The presence of the outflow in the simulated galaxy is also confirmed by the
X-shaped magnetic field orientation shown in Fig. 4. The magnetic field
direction without and with Faraday rotation considered are shown in the left
and middle panels, respectively. As shown in Equation (14), the rotation
measure is proportional to the thermal electron density and the magnetic
field, both of which have larger values in the disk midplane (see top panel of
Fig. 8). Therefore, the rotation measure is highest within the disk and
decreases toward the halo. As a result, comparing the left and middle panels,
the halo regions ($|z|>1$ kpc), which are relatively unaffected by Faraday
rotation, clearly show the X-shape pattern of the magnetic field indicating
the existence of an outflow. On the other hand, the magnetic field direction
near the disk is parallel to the disk midplane when the Faraday rotation is
not considered (left panel) and is strongly distorted by Faraday rotation as
the polarization angles can be Faraday rotated by at most 30 degrees at 6 GHz
in the disk (see bottom panel of Fig. 8). In the right panel of Fig. 4, we
present the field direction overlaid on top of the polarization degree with
Faraday rotation included. The polarization degree is smaller when the field
is more tangled, but it can reach the typical maximum synchrotron polarization
degree ($\sim 70-75$ percent for electron spectral index in between 2.1 and
3.1) when the field is more ordered, such as in the X-shape region and just
above the disk.

Figure 5:  Comparison of the X-shaped magnetic field morphology and the
polarized intensity between simulation and observation. The contours show
polarized intensity at 0.015, 0.050, 0.1, 0.5, and 1 mJy with beam
depolarization considered. The white-dashed lines show the magnetic field
direction, plotted on top of the total emission map in 6 GHz. The left panel
depicts the simulation results, while the right panel shows the real
observation (same data as in Figure 11 in Stein et al. 2020). The upper panel
corresponds to the field direction without the influence of Faraday rotation,
whereas the lower panel shows the field direction with Faraday rotation
included. The magnetic field direction exhibits a characteristic X-shaped
pattern indicative of an outflow. In lower panels, the field orientation
within the disk region is Faraday de-rotated, while it continues to display an
X-shaped pattern at $|z|>1$ kpc.

Fig. 5 shows mock observations (left column) and observations (right column).
The background shows total radio emission convolved to a beam size of 7.7
arcsec ($\sim 0.75$ kpc), while the contours show the polarized intensity at
0.015, 0.05, 0.1, 0.5, and 1 mJy beam-1 when applicable since the maximum
value of the four panels is different. Magnetic field orientations in regions
where the polarized intensity is larger than 0.015 mJy beam-1 are shown as
white dashed lines. Beam-depolarization is included by smoothing the Stokes Q
and U maps with a 7.7 arcsec beam in diameter before calculating the polarized
intensity and magnetic field orientation.

The upper (lower) panels show results without (with) Faraday rotation. In the
mock observation, observational noise generated by a Fourier transform
approach is added to the total intensity shown in the background (in the same
fashion as in Section 3.2). The upper-left panel shows a significantly larger
(maximum contour being 1 mJy/beam) polarized intensity than the upper-right
panel (maximum contour being 0.05 mJy beam-1). This is in part because the
observational image in the upper-right panel is obtained by Faraday de-
polarizing the lower-right panel, while the simulation data can be used to
directly calculate the polarized intensity without Faraday rotation and, thus,
without any information loss. When Faraday rotation is considered (lower left
panel), the simulated polarized intensity is still larger (maximum contour
being 0.1 mJy beam-1) than the observed polarized intensity because we do not
model observational noise for polarized intensity.
In Fig. 5, the field orientations in the simulation exhibit a pattern
consistent with the observational data. Both the simulation and observation
reveal that near the disk the field directions are parallel to the midplane
when the Faraday rotation is not considered (upper panel), and become
disordered when Faraday rotation is included (lower panel). At higher
latitudes, the B-field directions are not parallel, but instead slightly
tilted or even perpendicular to the midplane, suggesting the presence of an
outflow. This X-shaped feature persists even when Faraday rotation is
included, because it primarily affects the emission from the vicinity of the
midplane. Fig. 4 also shows that the field directions can still be
perpendicular to the midplane at even higher latitudes (than those shown in
Fig. 5), where the polarized intensity is too faint to be detected.

###  3.4 Radio Spectrum

We compare the simulated and observed spectra of NGC 4217 in Fig. 6. The
references for the observational data are listed in Table 1. In order to
facilitate a meaningful comparison of the simulations and observations, we
include data points from various sources that provide measurements of the
total flux (i.e., solid angle-integrated) of NGC 4217 and exclude those that
only record a peak flux. The spectrum predicted from the simulation and
boosted magnetic field (solid red line) is consistent with the observations
for frequencies $\nu>0.1~{}{\rm GHz}$. The original simulated spectrum without
magnetic field boosting but renormalized to match the observations is shown as
a dashed blue line. The comparison of these simulation-based spectra shows
that the slope of the spectrum is very weakly sensitive to the normalization
of the magnetic field. Both the observed and simulated spectrum for
frequencies $\nu<0.1$ GHz are not reliable. On the observation side, the low
frequency (and hence low resolution) measurement from Marvil et al. (2015) at
74 MHz is determined by measuring the peak flux and multiplying it by the
ratio of the solid-angle-integrated flux to peak flux based on a measurement
at a higher frequency, which has better resolution and for which the target is
resolved. However, the integrated-to-peak flux ratio is not constant, but
increases with decreasing frequency. For example, the ratio of solid-angle-
integrated flux to the disk flux for NGC 4217 increases from 2.35 at 6 GHz to
5.4 at 0.15 GHz as measured in Stein et al. (2020). Therefore, using the ratio
of integrated to peak flux at a higher frequency to infer low-frequency
measurement will lead to an underestimate of the real total flux since the
actual ratio at low frequency is higher than it is at higher frequency. On the
simulation side, while the free-free absorption could significantly attenuate
the low-frequency spectrum, the modeling of free-free opacity depends
sensitively on the simulation resolution and the detailed ISM modeling. For
example, increasing resolution at the disk region may lead to increased
cooling and hence more absorption (Hummels et al., 2019; Nelson et al., 2020;
Esmerian et al., 2021). In conclusion, our simulation shows consistency with
observations for $\nu>0.1$ GHz, while the spectrum at low frequency ($\nu<0.1$
GHz) remains uncertain.

Figure 6: Observed (gray points with error bars) and simulated (solid red line) radio spectrum of NGC 4217. The simulated spectrum with a boosted magnetic field is consistent with the observations for $\nu>0.1~{}{\rm GHz}$. The original simulated spectrum, renormalized to observation at 1.58 GHz, is shown as the dashed blue line and demonstrates that the spectral slope is very weakly sensitive to the magnetic filed. Both the observational and simulation data for $\nu<0.1$ GHz (shaded region) are less reliable due to the low resolution of the observations at low frequencies and finite numerical resolution when modeling free-free absorption. Observational data and the associated references are given in Table 1. Table 1: Total flux of NGC 4217 Freq. (GHz) | Flux (Jy) | Reference
---|---|---
0.0575 | 0.3$\pm$0.2 | Israel & Mahoney (1990)
0.074 | 0.357$\pm$0.097 | Marvil et al. (2015)
0.15 | 0.608$\pm$0.121 | Shimwell et al. (2017)
0.325 | 0.291$\pm$0.009 | Marvil et al. (2015)
0.408 | 0.220$\pm$0.0205 | Ficarra et al. (1985)
1.4 | 0.123$\pm$0.0044 | Condon et al. (1998)
1.58 | 0.108$\pm$0.005 | Stein et al. (2020)
4.85 | 0.04$\pm$0.007 | Gregory & Condon (1991)
5 | 0.05$\pm$0.011 | Sramek (1975)
6 | 0.033$\pm$0.002 | Stein et al. (2020)

##  4 Conclusions

In this work, we model the radio emission from an edge-on galaxy using a two-
moment CR magnetohydrodynamical simulation (Thomas et al., 2024) and compare
the results to the edge-on galaxy NGC 4217 (Stein et al., 2020). We compute
the spectra of CR protons, primary electrons, and secondary electrons adopting
a cell-based steady-state approach (Werhahn et al., 2021b). We then calculate
thermal and non-thermal radio emission and absorption processes and the Stokes
parameters including Faraday rotation and beam depolarization.
We compare observations of the latitude- and longitude-dependent intensity
profiles obtained following a procedure similar to that adopted by Stein et
al. (2020). The intensity profiles at both 1.6 and 6 GHz are consistent with
the observations within the error bars when the magnetic field is moderately
boosted. We anticipate that the magnetic field could be increased by a factor
of a few in a cosmological simulation that self-consistently includes the
evolution of the CGM. Although the magnitude of radio emission depends on the
magnetic field, we find that the shape of the intensity profile is only very
weakly dependent on the strength of the magnetic field and is a robust
prediction of our model.
Our simulated galaxy exhibits a morphology similar to that of NGC 4217 in
terms of radio emission and magnetic field orientation. The emission map of
our simulated galaxy shows a spatially-extended wind-like structure with
emission evenly distributed across the mid-plane, which better resembles
galaxies that have no centrally biased star formation rate density. Such
centrally peaked star formation events can be realized in starburst systems or
systems with a strong stellar bar feature.

This extended wind-launching disk naturally arises in (i) our sophisticated
ISM crisp model that allows the gas to cool to lower temperatures and to
generate more stellar winds and supernovae to prevent star formation at the
center of unbarred galaxies, and (ii) because we effectively simulate later
evolutionary stages of the galaxy with higher mass resolution. Our simulations
also reveal X-shaped magnetic field patterns at heights $\sim$1 kpc above the
disk midplane irrespective of the presence of Faraday rotation. This feature
is also present in the observations and indicates the existence of an outflow
in both the simulation and observation.
Finally, the simulated spectrum of our galaxy matches that of NGC 4217 for
$\nu>0.1$ GHz, while the observed spectrum at lower frequencies remains
uncertain. The shape of the spectrum is also essentially insensitive to the
magnetic field strength. Our findings demonstrate the effectiveness of the
two-moment cosmic ray transport model in generating realistic outflows on
galactic scales and provide a deeper understanding of galactic wind dynamics
and its influence on the evolution of galaxies.

## Acknowledgements

HHC and MR thank Yelena Stein, Volker Heesen, and Jiangtao Li for their
valuable suggestions on observational data processing. MR acknowledges support
from the National Aeronautics and Space Administration grant ATP 80NSSC23K0014
and the National Science Foundation Collaborative Research Grant NSF
AST-2009227. TT and CP acknowledge support by the European Research Council
under ERC-AdG grant PICOGAL-101019746. This work was supported by the North-
German Supercomputing Alliance (HLRN) under project bbp00070. This research
was supported in part by grant NSF PHY-2309135 to the Kavli Institute for
Theoretical Physics (KITP). This work was performed in part at the Aspen
Center for Physics, which is supported by National Science Foundation grant
PHY-2210452, and by grants from the Simons Foundation (1161654, Troyer), and
Alfred P Sloan Foundation (G-2024-22395). MR acknowledges the support from the
Leinweber Center for Theoretical Physics which sponsored the organization of
the 7th ICM Theory and Computation Workshop during which HHSC, MR, and CP were
able to finalize key aspects of the project. MR thanks Volker Springel and the
Max Planck Institute for Astrophysics (MPA) in Garching for their hospitality
during his sabbatical stay at MPA. MR acknowledges Forschungsstipendium from
MPA.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

##  Appendix A Verifying the steady-state assumption

The steady state is defined by the requirement that

$\dot{e}_{\rm all}+\dot{e}_{\rm inj}\equiv\dot{e}_{\rm CR}=0,$ (A1)

where $\dot{e}_{\rm all}$ is the rate of CR energy change that includes
contributions from various cooling processes as well as the advection of CRs
in and out of the computational cells. That is, $e_{\rm CR}/|\dot{e}_{\rm
all}|\equiv\tau_{\rm all}$ includes all CR energy density changes other than
an imposed CR injection rate $\dot{e}_{\rm inj}$ that is necessary to keep the
CR distribution in a steady state (as assumed in the approach used to post-
process the simulations to generate steady-state radiation spectra). Thus, the
steady-state solution is exactly correct when $e_{\rm CR}/|\dot{e}_{\rm CR}|$
is infinite. In other words, the system is in a near steady state when
$|\dot{e}_{\rm all}|\sim|\dot{e}_{\rm inj}|$ or $|\dot{e}_{\rm
all}+\dot{e}_{\rm inj}|\ll(|\dot{e}_{\rm all}|+|\dot{e}_{\rm inj}|)/2$. In
that case, $e_{\rm CR}/|\dot{e}_{\rm all}|=\tau_{\rm all}\ll e_{\rm
CR}/|\dot{e}_{\rm CR}|\equiv\tau_{\rm CR}$, and so the system is close to a
steady state when $\tau_{\rm all}\ll\tau_{\rm CR}$ in which case the
application of the steady-state post-processing technique leads to near exact
results.
The CR energy change timescale $\tau_{\rm CR}$ is calculated by comparing the
CR energy in each computational cell of two adjacent snapshots separated in
time by $\mathrm{d}t\sim 0.1$ Myr. The energy loss timescale includes the
escape timescale $\tau_{\rm esc}$ and the cooling timescale $\tau_{\rm cool}$
and is given by

$\frac{1}{\tau_{\rm all}}=\frac{1}{\tau_{\rm esc}}+\frac{1}{\tau_{\rm cool}},$
(A2)

where the escape timescale is a combination of diffusion and advection
timescale $\tau_{\rm esc}^{-1}=\tau_{\rm diff}^{-1}+\tau_{\rm adv}^{-1}$ (see
Section 2.2 for $\tau_{\rm diff}$ and $\tau_{\rm adv}$) and the cooling
timescale is given by the total energy divided by the proton energy loss rate
$\tau_{\rm cool}=e_{\rm CR}/|b_{\rm coul}+b_{\rm pion}|$, where $b_{\rm coul}$
and $b_{\rm pion}$ are the energy loss rates due to Coulomb and hadronic
interactions, respectively. We calculate the escape and cooling timescales at
CR energy of 10 GeV. In order to verify the cell-based steady-state
assumption, we show in the left panel of Fig. 7 the ratio of timescales
computed as $\left<\log\left(\tau_{\rm CR}/\tau_{\rm all}\right)\right>$,
where the average is intensity-weighted.
We motivate this choice as follows. While $\log(\left<\tau_{\rm all}/\tau_{\rm
CR}\right>^{-1})$ is biased toward small ratios and $\log(\left<\tau_{\rm
CR}/\tau_{\rm all}\right>)$ is biased toward large ratios,
$\left<\log\left(\tau_{\rm CR}/\tau_{\rm all}\right)\right>$ provides a
typical in-between value of the timescale ratio and is independent of how we
define the timescale ratios (i.e., $\left<\log\left(\tau_{\rm CR}/\tau_{\rm
all}\right)\right>=-\left<\log\left(\tau_{\rm all}/\tau_{\rm
CR}\right)\right>$). In the absence of intensity weighting, the above
logarithms of averages would translate into the logarithm of the harmonic
$\left<\tau_{\rm CR}/\tau_{\rm all}\right>_{\rm H}$, geometric
$\left<\tau_{\rm CR}/\tau_{\rm all}\right>_{\rm G}$, and arithmetic average
$\left<\tau_{\rm CR}/\tau_{\rm all}\right>_{\rm A}$, respectively, which are
strictly ordered as $\left<\tau_{\rm CR}/\tau_{\rm all}\right>_{\rm
H}\leq\left<\tau_{\rm CR}/\tau_{\rm all}\right>_{\rm G}\leq\left<\tau_{\rm
CR}/\tau_{\rm all}\right>_{\rm A}$. In practice, we retain intensity weighting
to ensure that the computational cells that contribute most to intensity can
properly shape the timescale ratio.

Figure 7:  Left column: Adequacy of the steady-state assumption. Upper left
panel: projection of the timescale ratio within the whole box ($y=-10$ kpc to
$y=10$ kpc; upper panel) and a thin slice perpendicular to the line of sight
anchored at $y=0$ (middle panel) both weighted by total radio emissivity;
Lower left panel: histogram of weighted timescale ratios. The projection plots
show that the steady-state assumption is adequate for heights from $z=0$ to
$z\pm 5$ kpc (Fig. 2), while the histogram shows that the timescale ratio is
greater than 1 in cells that contribute most to the radio spectrum (Fig. 6).
Right column: Adequacy of the steady-state assumption for electrons. Upper-
right panel: projection of the cooling timescale ratio between protons and
electrons within the whole box ($y=-10$ kpc to $y=10$ kpc; upper panel) and a
thin slice perpendicular to the line of sight anchored at $y=0$ (middle panel)
both weighted by total radio emissivity; Lower-right panel: histogram of
weighted cooling timescale ratios. The projection plot shows that the steady-
state assumption is adequate for CR electrons for our intensity profile and
the radio spectrum, assuming that the energy density of electrons is always a
fraction of that of protons.

We evaluate timescale ratios in each computational cell and plot (i) a
projected map weighted by total radio emissivity in a box (20 kpc on a side;
upper panel) and a slice (0.5 kpc thick and perpendicular to the line of
sight; middle panel) and (ii) a histogram weighted by CR energy density and
total radio emissivity. Although there are some clumpy regions in the slice
plot where the $\tau_{\rm CR}/\tau_{\rm all}$ is smaller than 1 (shown in
blue), the radio emission-weighted emission average projection shows that the
ratio is around $10^{0.5}\sim 3$, meaning that the steady-state assumption is
approximately valid in the region where we calculate the intensity profile as
a function of height (Fig. 2). The histogram, showing the ratio distribution
of the most luminous computational cells, verifies the adequacy of the steady-
state assumption for the purpose of calculating the radio spectrum (Fig. 6).
In order to check the adequacy of the steady-state assumption for electrons,
we show in the right panel of Fig. 7 the distribution of the ratio between the
energy loss timescale for CR protons, $\tau_{\rm all,CRp}$, and electrons,
$\tau_{\rm all,CRe}$, with the same procedure as in the left panel of Fig. 7.
Both the maps and histogram show that the $\tau_{\rm all,CRe}$ is shorter than
$\tau_{\rm all,CRp}$, which is generally true for relativistic protons
(Ruszkowski & Pfrommer, 2023). In our simulation, the cooling timescale for
CRe and CRp are dominated by advective, diffusive, and streaming timescales,
which is identical for CRe and CRp. Therefore, $\tau_{\rm
all,CRp}\approx\tau_{\rm all,CRe}$ corresponds to regions where advective,
diffusive or streaming dominate while the tail of the distribution ($\tau_{\rm
all,CRp}>\tau_{\rm all,CRe}$) implies that the radiative leptonic loss rates
dominate over losses associated with spatial transport.
Assuming that the energy density of CR electrons is always some fraction of CR
protons, $\tau_{\rm CR}$ is then identical for protons and electrons.
Therefore, we obtain $\tau_{\rm all,CRe}<\tau_{\rm all,CRp}<\tau_{\rm CR}$ and
verify that the steady-state assumption is also fulfilled for CR electrons.
All ratios calculated here adopt the original magnetic field strength in the
simulation and thus correspond to the lower limit of $\tau_{\rm CR}/\tau_{\rm
all,CRe}$ because increasing the magnetic field strength will also increase
synchrotron losses and decrease the cooling timescale.

##  Appendix B Rotation measure

We show the rotation measure (RM) of the simulated galaxy in Fig. 8. The upper
panel shows the RM value in units of radian per $\rm m^{2}$ with RM values
near the disk midplane reaching values of $\sim\pm 200\ {\rm rad\ m^{-2}}$.
The lower panel shows the total rotation angle of the magnetic field along the
integrated line of sight at a frequency of 6 GHz to enable easy comparison to
Fig. 5, which shows that the field direction is mainly rotated near the disk
midplane, where the polarization angles can be rotated at most by $\sim 30$
degrees. Both panels show a transition of the RM values from positive to
negative at the galactic center, which is also seen in the RM of the Milky Way
(Reissl et al., 2023). This transition corresponds to the reversal of the
magnetic field direction and hence the toroidal magnetic field structure
feature in our simulation.

Figure 8: Rotation measure of the simulated galaxy. The top panel shows the
Faraday rotation angle per wavelength square, and the bottom panel shows the
total rotation angle in degrees at 6 GHz that can be used to compare with Fig.
5.

##  Appendix C Intensity profiles for individual data strips

We show the raw data of the average intensity along five different strips in
our simulation to demonstrate the large variation in the direction parallel to
the disk. For example, in Fig. 10, the intensity at $x=0$ increases from
around 2 mJy/beam to 10 mJy/beam between strip I and strip III; Fig. 11 also
shows about half an order of magnitude variation between strips.

Figure 9: Simulated intensity with five strips defined to fit the intensity
profile. Figure 10: Average intensity profile of five different strips at 1.6
GHz. Gray dots present the observation; the blue dashed line shows the profile
from the original simulation normalized to the total flux; the red line shows
the result from the simulation with a boosted magnetic field by a factor of
four. The profiles vary significantly between different strips.  Figure 11:
Same as Fig. 10, but shown at a frequency of 6 GHz.

## References

  * Abrahamsson et al. (2007) Abrahamsson, E., Krems, R. V., & Dalgarno, A. 2007, ApJ, 654, 1171, doi: 10.1086/509631
  * Acciari et al. (2009) Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Nature, 462, 770, doi: 10.1038/nature08557
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164, doi: 10.1088/0004-637X/755/2/164
  * Aharonian et al. (2010) Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D, 82, 043002, doi: 10.1103/PhysRevD.82.043002
  * Anderson & Bregman (2010) Anderson, M. E., & Bregman, J. N. 2010, ApJ, 714, 320, doi: 10.1088/0004-637X/714/1/320
  * Armillotta et al. (2021) Armillotta, L., Ostriker, E. C., & Jiang, Y.-F. 2021, ApJ, 922, 11, doi: 10.3847/1538-4357/ac1db2
  * Armillotta et al. (2022) —. 2022, ApJ, 929, 170, doi: 10.3847/1538-4357/ac5fa9
  * Bland-Hawthorn et al. (2007) Bland-Hawthorn, J., Veilleux, S., & Cecil, G. 2007, Ap&SS, 311, 87, doi: 10.1007/s10509-007-9567-8
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Physics Reports, 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Blasi (2010) Blasi, P. 2010, MNRAS, 402, 2807, doi: 10.1111/j.1365-2966.2009.16110.x
  * Bolatto et al. (2013) Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450, doi: 10.1038/nature12351
  * Bregman (2007) Bregman, J. N. 2007, ARA&A, 45, 221, doi: 10.1146/annurev.astro.45.051806.110619
  * Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79
  * Bustard et al. (2016) Bustard, C., Zweibel, E. G., & D’Onghia, E. 2016, ApJ, 819, 29, doi: 10.3847/0004-637X/819/1/29
  * Bykov (2014) Bykov, A. M. 2014, A&ARv, 22, 77, doi: 10.1007/s00159-014-0077-8
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Cen (1992) Cen, R. 1992, ApJS, 78, 341, doi: 10.1086/191630
  * Chevalier & Clegg (1985) Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44, doi: 10.1038/317044a0
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693, doi: 10.1086/300337
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Fossi, B. C. M., & Young, P. R. 1997, A&AS, 125, 149, doi: 10.1051/aas:1997368
  * Esmerian et al. (2021) Esmerian, C. J., Kravtsov, A. V., Hafen, Z., et al. 2021, MNRAS, 505, 1841, doi: 10.1093/mnras/stab1281
  * Farber et al. (2018) Farber, R., Ruszkowski, M., Yang, H. Y. K., & Zweibel, E. G. 2018, ApJ, 856, 112, doi: 10.3847/1538-4357/aab26d
  * Ficarra et al. (1985) Ficarra, A., Grueff, G., & Tomassetti, G. 1985, A&AS, 59, 255.  https://ui.adsabs.harvard.edu/abs/1985A&AS...59..255F
  * Gaisser (1990) Gaisser, T. K. 1990, Cosmic rays and particle physics. Cambridge Univ. Press, Cambridge
  * Grassi et al. (2014) Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386, doi: 10.1093/mnras/stu114
  * Gregory & Condon (1991) Gregory, P. C., & Condon, J. J. 1991, ApJS, 75, 1011, doi: 10.1086/191559
  * Heesen et al. (2016) Heesen, V., Dettmar, R.-J., Krause, M., Beck, R., & Stein, Y. 2016, MNRAS, 458, 332, doi: 10.1093/mnras/stw360
  * Heesen et al. (2018) Heesen, V., Krause, M., Beck, R., et al. 2018, MNRAS, 476, 158, doi: 10.1093/mnras/sty105
  * Heywood et al. (2019) Heywood, I., Camilo, F., Cotton, W. D., et al. 2019, Nature, 573, 235, doi: 10.1038/s41586-019-1532-5
  * Hopkins et al. (2012) Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522, doi: 10.1111/j.1365-2966.2012.20593.x
  * Hopkins et al. (2021) Hopkins, P. F., Squire, J., Chan, T. K., et al. 2021, MNRAS, 501, 4184, doi: 10.1093/mnras/staa3691
  * Huang & Davis (2022) Huang, (黄小珊), X., & Davis, S. W. 2022, MNRAS, 511, 5125, doi: 10.1093/mnras/stac059
  * Hummels et al. (2019) Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156, doi: 10.3847/1538-4357/ab378f
  * Israel & Mahoney (1990) Israel, F. P., & Mahoney, M. J. 1990, ApJ, 352, 30, doi: 10.1086/168513
  * Jiang & Oh (2018) Jiang, Y.-F., & Oh, S. P. 2018, ApJ, 854, 5, doi: 10.3847/1538-4357/aaa6ce
  * Jubelgas et al. (2008) Jubelgas, M., Springel, V., Enßlin, T., & Pfrommer, C. 2008, A&A, 481, 33, doi: 10.1051/0004-6361:20065295
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018, doi: 10.1103/PhysRevD.74.034018
  * Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445, doi: 10.1086/149981
  * Lacki et al. (2010) Lacki, B. C., Thompson, T. A., & Quataert, E. 2010, ApJ, 717, 1, doi: 10.1088/0004-637X/717/1/1
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Martin-Alvarez et al. (2018) Martin-Alvarez, S., Devriendt, J., Slyz, A., & Teyssier, R. 2018, MNRAS, 479, 3343, doi: 10.1093/mnras/sty1623
  * Marvil et al. (2015) Marvil, J., Owen, F., & Eilek, J. 2015, AJ, 149, 32, doi: 10.1088/0004-6256/149/1/32
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, Astronomical Data Analysis Software and Systems XVI, 376, 127.  https://ui.adsabs.harvard.edu/abs/2007ASPC..376..127M
  * Menon et al. (2022) Menon, S. H., Federrath, C., & Krumholz, M. R. 2022, MNRAS, 517, 1313, doi: 10.1093/mnras/stac2702
  * Miller (1991) Miller, J. A. 1991, ApJ, 376, 342, doi: 10.1086/170284
  * Moseley et al. (2021) Moseley, E. R., Draine, B. T., Tomida, K., & Stone, J. M. 2021, MNRAS, 500, 3290, doi: 10.1093/mnras/staa3384
  * Moster et al. (2010) Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903, doi: 10.1088/0004-637X/710/2/903
  * Murray et al. (2011) Murray, N., Ménard, B., & Thompson, T. A. 2011, ApJ, 735, 66, doi: 10.1088/0004-637X/735/1/66
  * Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569, doi: 10.1086/426067
  * Naab & Ostriker (2017) Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59, doi: 10.1146/annurev-astro-081913-040019
  * Nelson et al. (2020) Nelson, D., Sharma, P., Pillepich, A., et al. 2020, MNRAS, 498, 2391, doi: 10.1093/mnras/staa2419
  * Paglione & Abrahams (2012) Paglione, T. A. D., & Abrahams, R. D. 2012, ApJ, 755, 106, doi: 10.1088/0004-637X/755/2/106
  * Pakmor et al. (2014) Pakmor, R., Marinacci, F., & Springel, V. 2014, ApJL, 783, L20, doi: 10.1088/2041-8205/783/1/L20
  * Pakmor et al. (2016) Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ApJ, 824, L30, doi: 10.3847/2041-8205/824/2/L30
  * Pakmor et al. (2016) Pakmor, R., Springel, V., Bauer, A., et al. 2016, MNRAS, 455, 1134, doi: 10.1093/mnras/stv2380
  * Pakmor et al. (2017) Pakmor, R., Gómez, F. A., Grand, R. J. J., et al. 2017, MNRAS, 469, 3185, doi: 10.1093/mnras/stx1074
  * Pfrommer et al. (2008) Pfrommer, C., Enßlin, T. A., & Springel, V. 2008, MNRAS, 385, 1211, doi: 10.1111/j.1365-2966.2008.12956.x
  * Pfrommer et al. (2017) Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017, Mon. Not. R. Astron. Soc., 465, 4500, doi: 10.1093/mnras/stw2941
  * Pfrommer et al. (2017) Pfrommer, C., Pakmor, R., Simpson, C. M., & Springel, V. 2017, ApJ, 847, L13, doi: 10.3847/2041-8213/aa8bb1
  * Pfrommer et al. (2022) Pfrommer, C., Werhahn, M., Pakmor, R., Girichidis, P., & Simpson, C. M. 2022, MNRAS, 515, 4229, doi: 10.1093/mnras/stac1808
  * Ponnada et al. (2024) Ponnada, S. B., Panopoulou, G. V., Butsky, I. S., et al. 2024, MNRAS, 527, 11707, doi: 10.1093/mnras/stad3978
  * Reissl et al. (2018) Reissl, S., Klessen, R. S., Mac Low, M.-M., & Pellegrini, E. W. 2018, A&A, 611, A70, doi: 10.1051/0004-6361/201731698
  * Reissl et al. (2023) Reissl, S., Klessen, R. S., Pellegrini, E. W., et al. 2023, Nature Astronomy, 7, 1295, doi: 10.1038/s41550-023-02053-2
  * Rosdahl et al. (2015) Rosdahl, J., Schaye, J., Teyssier, R., & Agertz, O. 2015, MNRAS, 451, 34, doi: 10.1093/mnras/stv937
  * Ruszkowski & Pfrommer (2023) Ruszkowski, M., & Pfrommer, C. 2023, A&A Rev., 31, 4, doi: 10.1007/s00159-023-00149-2
  * Ruszkowski et al. (2017) Ruszkowski, M., Yang, H. Y. K., & Zweibel, E. 2017, ApJ, 834, 208, doi: 10.3847/1538-4357/834/2/208
  * Rybicki & Lightman (1986) Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics
  * Salem & Bryan (2014) Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312, doi: 10.1093/mnras/stt2121
  * Shimwell et al. (2017) Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017, A&A, 598, A104, doi: 10.1051/0004-6361/201629313
  * Socrates et al. (2008) Socrates, A., Davis, S. W., & Ramirez-Ruiz, E. 2008, ApJ, 687, 202, doi: 10.1086/590046
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel & Hernquist (2003) Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289, doi: 10.1046/j.1365-8711.2003.06206.x
  * Sramek (1975) Sramek, R. 1975, AJ, 80, 771, doi: 10.1086/111810
  * Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289, doi: 10.1088/0004-637X/717/1/289
  * Stein et al. (2020) Stein, Y., Dettmar, R.-J., Beck, R., et al. 2020, A&A, 639, A111, doi: 10.1051/0004-6361/202037675
  * Thomas & Pfrommer (2019) Thomas, T., & Pfrommer, C. 2019, MNRAS, 485, 2977, doi: 10.1093/mnras/stz263
  * Thomas & Pfrommer (2022) —. 2022, MNRAS, 509, 4803, doi: 10.1093/mnras/stab3079
  * Thomas et al. (2020) Thomas, T., Pfrommer, C., & Enßlin, T. 2020, ApJL, 890, L18, doi: 10.3847/2041-8213/ab7237
  * Thomas et al. (2021) Thomas, T., Pfrommer, C., & Pakmor, R. 2021, MNRAS, 503, 2242, doi: 10.1093/mnras/stab397
  * Thomas et al. (2023) —. 2023, MNRAS, 521, 3023, doi: 10.1093/mnras/stad472
  * Thomas et al. (2024) Thomas, T., Pfrommer, C., & Pakmor, R. 2024, arXiv e-prints, arXiv:2405.13121, doi: 10.48550/arXiv.2405.13121
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389, doi: 10.1146/annurev-astro-091916-055240
  * Uhlig et al. (2012) Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374, doi: 10.1111/j.1365-2966.2012.21045.x
  * Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769, doi: 10.1146/annurev.astro.43.072103.150610
  * Vink (2011) Vink, J. 2011, A&AR, 20, 49, doi: 10.1007/s00159-011-0049-1
  * Waelkens et al. (2009) Waelkens, A., Jaffe, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A, 495, 697, doi: 10.1051/0004-6361:200810564
  * Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238, doi: 10.1093/mnras/stv1975
  * Werhahn et al. (2021a) Werhahn, M., Pfrommer, C., & Girichidis, P. 2021a, MNRAS, 508, 4072, doi: 10.1093/mnras/stab2535
  * Werhahn et al. (2021b) Werhahn, M., Pfrommer, C., Girichidis, P., Puchwein, E., & Pakmor, R. 2021b, MNRAS, 505, 3273, doi: 10.1093/mnras/stab1324
  * Werhahn et al. (2021c) Werhahn, M., Pfrommer, C., Girichidis, P., & Winner, G. 2021c, MNRAS, 505, 3295, doi: 10.1093/mnras/stab1325
  * Wiener et al. (2017) Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906, doi: 10.1093/mnras/stx127
  * Yang et al. (2018) Yang, R.-z., Kafexhiu, E., & Aharonian, F. 2018, A&A, 615, A108, doi: 10.1051/0004-6361/201730908
  * Yoast-Hull et al. (2013) Yoast-Hull, T. M., Everett, J. E., Gallagher, J. S., & Zweibel, E. G. 2013, ApJ, 768, 53, doi: 10.1088/0004-637X/768/1/53
  * Zirakashvili & Aharonian (2007) Zirakashvili, V. N., & Aharonian, F. 2007, A&A, 465, 695, doi: 10.1051/0004-6361:20066494