re-forming the clouds. The net effect is that the momentum of the hot wind is
transferred to the cold phase that survives the exposure to the hot wind
(Gronke & Oh, 2018, 2020; Sparre et al., 2020; Li et al., 2020a; Girichidis et
al., 2021; Farber & Gronke, 2021).
As mentioned above, cold clouds can be accelerated by CR pressure gradients.
While an accelerating cold cloud exposed to a CR gradient can also undergo
shattering, the tendency of a cloud to disrupt can be suppressed if the Alfvén
speed profile across the cloud has sharp edges at the cloud interface (Thomas
et al., 2021). Furthermore, even when the Alfvén speed profiles do not
resemble inverted top hat profiles, two-dimensional simulations with radiative
cooling demonstrate that the clouds can survive, and the survival is even
possible despite CR heating operating at cloud-ISM interfaces (Wiener et al.,
2019). In the simulations of this process without CR heating, the cloud
destruction is primarily driven by CR streaming on the sides of the cloud,
which results in a significant shear at the boundary layer (Brüggen &
Scannapieco, 2020). However, the cloud destruction rates can be slow enough to
allow for cloud survival and for cold cloud acceleration to hundreds of km s-1
as seen in the observed galactic winds.

##### Dynamical impact of cosmic rays in porous multiphase media.

In realistic situations, CRs interact with multiphase gas characterized by a
range of densities, temperatures, and ionization fractions. Additional
physical processes shape the dynamical interactions of CRs with such porous
media. Specifically, in the mostly neutral gas phase, the effective CR
streaming speed can vastly exceed the speed in the fully ionized gas because,
in the self-confinement model in the fully coupled regime, CR transport
proceeds at the ion Alfvén speed, i.e.,
$\varv_{\mathrm{a}}=B/\sqrt{4\pi\rho}$, where $B$ is the magnetic field
strength and $\rho$ is the ion mass density. Furthermore, in the weakly
ionized phase, ion-neutral friction can very efficiently damp self-generated
Alfvén waves on which CRs scatter (Kulsrud & Pearce, 1969). Consequently, in
the mostly neutral gas, CR can propagate much faster than in the fully ionized
medium. These processes have interesting implications not only for the CR
transport but also for the dynamical coupling of CRs to the gas and heating of
the gas by CRs.
The inclusion of ion-neutral damping leads to subtle yet important
modifications to the bottleneck process described above. As CRs stream down
their pressure gradient toward an overdense cloud, they propagate in a
progressively denser medium, and as long as the magnetic field remains
constant, the Alfvén speed decreases. As before, this creates a bottleneck
effect and the CR pressure flattens in the upstream region and drops
significantly past the minimum in the transport speed. However, the transport
speed eventually increases significantly as CRs propagate deeper into the
cloud because the ionization fraction drops and the ion-neutral friction damps
the confining waves. This effect can overcompensate for the initial drop in
the transport speed. Consequently, this very fast transport inside the cloud
leads to the reduction in the CR pressure. As the CRs traverse the cloud, they
eventually encounter the far side of the cloud, where the density decreases,
ionization fraction increases, and ion-neutral damping becomes ineffective
again. The latter of these three factors is decisive in controlling the
effective CR transport speed and CRs reestablish strong coupling with the gas
on the far side of the cloud. This results yet again in the reduction of the
CR transport speed and, thus, a second bottleneck (Bustard & Zweibel, 2021).
Consequently, the reduced CR pressure remains flat inside the cloud just as it
did in front of the leading edge of the cloud, which leads to a stair-step
distribution of the CRs across the cloud.
Substantial CR collisionless losses can occur at the locations of sharp CR
pressure gradients at cloud edges. These losses are further increased because
of the ionization fraction dependence of the collisionless heating,
$|\boldsymbol{\varv}_{\mathrm{a}}\boldsymbol{\cdot}\boldsymbol{\nabla}P_{\rm
cr}|$. In addition to collisionless losses, CR also experience substantial
hadronic and Coulomb losses. In objects that are CR calorimeters, by
definition, these energy losses are comparable to the CR energy production
rate in SNe (Lacki, 2011; Yoast-Hull et al., 2015, 2016; Krumholz et al.,
2020). Despite the complexities introduced to CR transport by the varying
ionization fraction, the overall CR losses may be only weakly dependent on
this factor. In the low ionization case, the decrease in the CR residence time
in the clouds, and the associated reduction in the collisional losses, can be
compensated for by the CR probing denser medium in which the losses are larger
and by the fact that collisionless losses can occur on both cloud-ISM
interfaces rather than just one (Bustard & Zweibel, 2021). These
considerations have implications for the $\gamma$-ray emission which probe
hadronic interactions of CRs with the gas, which we comment on in Section 4.3.
The dynamical impact of CRs can also be influenced by the ionization structure
of the multiphase medium. In terms of the effective momentum imparted on to
the gas by CRs, the outcomes can be bracketed by two extreme limits. When the
medium is very close to being neutral, ion-neutral damping can be very
effective throughout the clouds thus completely flattening CR pressure
profiles in the cold and dense phase and making CR driving of outflows from
the galaxy ineffective (Everett & Zweibel, 2011). On the other hand, if the
medium is fully ionized, CR bottlenecks form and, as mentioned above, the
clouds can be accelerated to astrophysically relevant speeds (Wiener et al.,
2017a, 2019; Brüggen & Scannapieco, 2020; Thomas et al., 2021). Interestingly,
in the intermediate case of partial ionization, cloud acceleration rates
remain largely unaffected by the level of ionization and the clouds can be
accelerated to velocities still comparable to those found for full ionization
cases (Bustard & Zweibel, 2021).
Controlled and highly idealized experiments of individual clouds are helpful
for developing physical intuition and interpreting simulations in which the
formation of the multiphase gas is computed from first principles. Simulations
of thermal instability in magnetized media with CRs naturally predict phase
space distributions of such cold clouds (e.g., Sharma et al., 2010b; Huang et
al., 2022, see also Fig. 19, which shows density, velocity, and CR pressure
distributions in thermally unstable gas exposed to CR flux).

Figure 19: Top panel: thermally unstable gas density distribution exposed to
CR flux. Bottom panel: gas density, CR pressure, and projected gas velocity
along a magnetic field line (magnetic field line is denoted as a solid line in
the top panel). Formation of flat CR plateaus and bottlenecks coincident with
gas density maxima is evident in the bottom panel. Image from Huang et al.
(2022) reproduced with permission from ApJ.

Such simulations reveal that CRs reduce the density contrast between the
clouds and the hot ambient medium and that the clouds can have very
filamentary morphology. In a purely hydrodynamical cases, when clouds are
exposed to an incoming wind, they experience the strongest acceleration on the
side exposed to the wind. However, when the clouds are instead exposed to CR
flux, bottleneck formation can help to accelerate the clouds. Assuming an
Alfvén-speed profile that follows an inverted Gaussian, acceleration tends to
be stronger on the far side of the clouds with respect to the CR source, thus
producing a differential force stretching the clouds, which can make the
clouds filamentary. The net acceleration of the multiphase gas formed via
thermal instability tends to be suppressed compared to the expectations based
on idealized experiments described above. This can be attributed to two
factors. First, because in such thermally unstable environments magnetic
fields are tangled, and CR transport is anisotropic, the forces associated
with the bottlenecks tend to push the gas both in the direction parallel and
perpendicular to the overall wind direction. Second, the strength of the
bottlenecks is now reduced because, in the process of thermal instability
development, the decrease of the Alfvén speed in the clouds due to elevated
density is partially compensated for by the increased magnetic field strength
in the cold phase. While the momentum imparted to the gas via CR pressure
forces is reduced, the velocity dispersion of the gas increases.
Interestingly, the impact of the CR flux on the acceleration of the cold gas
is nevertheless significant because CR heating in the hot gas increases
thermal pressure that in turn is able to accelerate the cold gas to the
observed velocities (Huang et al., 2022).

####  3.2.3 Parker instability and cosmic rays

We now turn our attention to larger length scales and discuss the impact of
CRs on the stability of gravitationally stratified patches of the ISM. These
considerations are highly relevant for a range of topics including regulation
of star formation, magnetic dynamos, or polluting of the CGM with magnetic
fields and CRs, where the latter may dominate pressure support. While we defer
the discussion of these topics to later sections, here we focus specifically
on the triggering and development of the Parker instability and the role that
CR physics plays in it.
In the solar neighborhood, typical values of the energy densities of CRs,
magnetic fields, radiation field, and kinetic energy of the gas are comparable
and at the level of 1 eV cm-3 (e.g., Boulares & Cox, 1990; Badhwar & Stephens,
1977). One consequence of this apparent equipartition is that CRs and magnetic
fields can play a dynamically important role in setting up gravitational
equilibrium in the gaseous component of galactic disks (Badhwar & Stephens,
1977; Ghosh & Ptuskin, 1983; Chevalier & Fransson, 1984). We note in passing
that, while equipartition may be a reasonable assumption in normal star
forming galaxies (but see the discussion in Section 4.3.2), it may not hold in
starburst galaxies (Yoast-Hull et al., 2016). In the latter case, high gas
densities and mass loss due to winds, can lead to significant CR energy
density losses. Consequently, magnetic fields need to be much stronger in
order for the synchrotron emission to match the observed radio luminosities.
This also implies that magnetic fields must be strongly correlated with the
sites of significant star formation. In such situations, magnetic energy
density may in general be out of equipartition and may significantly exceed
the CR energy density.
Since the magnetic fields and CRs are essentially weightless components of the
ISM, a gravitationally stratified ISM may become unstable to a buoyancy
instability that is analogous to the Rayleigh-Taylor instability that is
triggered when a heavier fluid is placed on top of a lighter one. In the
context of magnetized medium with CRs, this instability is known as the Parker
instability (Parker, 1966). In a Parker unstable system, the gas begins to
flow down along a perturbed magnetic field. In the process of infall,
gravitational potential energy is released and spent on compressing the gas
and counteracting magnetic tension. Parker’s instability has been demonstrated
to play a role in the magnetic field amplification by dynamo action and in the
transport of magnetic fields and CRs out of the disk (Hanasz & Lesch, 2000;
Hanasz et al., 2004a; Hanasz et al., 2004b), which could aid magnetic field
growth via the small-scale dynamo (Brandenburg & Subramanian, 2005; Pakmor et
al., 2017; Pfrommer et al., 2022) and a large-scale $\alpha$-$\Omega$ dynamo
(Ntormousi et al., 2020). While it was originally envisioned that Parker’s
instability could be responsible for the formation of molecular clouds, in the
modern picture, supersonic ISM turbulence driven by gravitational collapse and
stellar winds or supernova explosions regulate star formation (Mac Low &
Klessen, 2004; Krumholz & McKee, 2005), and GMCs are only transient objects
that form as a result of converging turbulent flows and only survive a few
crossing time scales, which may explain the small star formation efficiency.
In the classic treatment of the Parker instability, Parker (1966) assumed that
the CR pressure is constant along all magnetic field lines. This assumption is
justified when the speed of sound of the composite fluid is much larger than
any characteristic velocities in the problem and the CR fluid is not affected
by gravity. These assumptions can be relaxed in several ways by considering
general CR pressure perturbations, which leads to several interesting
consequences that we discuss next. In brief, while treating CRs as a
compressible fluid described by an adiabatic index $\gamma_{\rm cr}>0$ makes
the gas less susceptible to the Parker instability, while including non-
adiabatic effects renders the gas more unstable.
First, the classic treatment can be extended by approximating the CR fluid as
a relativistic fluid with $\gamma_{\rm cr}=4/3$, which makes the gas less
compressible. This in turn implies that more work needs to be done by gravity
to compress the gas, which makes it less susceptible to the development of the
instability (Zweibel & Kulsrud, 1975; Boettcher et al., 2016; Heintz &
Zweibel, 2018). This could explain why, in simulations of stratified disks
involving adiabatic CRs, the atmospheres can remain puffed up by CR pressure
rather than become unstable to the Parker instability (e.g., Uhlig et al.,
2012; Ruszkowski et al., 2017a, see also Section 3.2.4 below).
Second, including non-adiabatic effects, such as radiative cooling, CR
diffusion, streaming, and thermal conduction have destabilizing effect. The
impact of radiative cooling is easiest to grasp. Radiative losses make the gas
more compressible, which makes it easier for the gas to flow deeper into the
valleys formed by the magnetic field and accelerate the instability. CR
transport via anisotropic diffusion along the magnetic fields also increases
the growth rate of the Parker instability. As the coupling of CRs to the gas
is reduced (i.e., diffusion is faster), so is the steepness of the CR pressure
gradient that opposes gravity, which allows the gas to fall down more easily.
Consequently, the growth rate of the instability increases with the parallel
diffusion coefficient as demonstrated both by linear stability analysis and
numerical simulations (e.g., Ryu et al. 2003; Kuwabara et al. 2004; Rodrigues
et al. 2016; but see Hanasz & Lesch 2003, who reported the opposite trend).
This is evident in Fig. 20 that illustrates the stark differences between the
simulations with and without CR physics. In particular, the emergence of
prominent magnetic “Parker loops” is clearly visible in the bottom panel that
shows the results from CR simulations that include anisotropic CR diffusion.

Figure 20: Iso-density surfaces in a gravitationally stratified patch of a
galactic disk (yellow surfaces), density slices (vertical planes), and
magnetic fields (green lines) in simulations with and without CRs (top and
bottom panel, respectively; the simulation with CRs includes anisotropic
diffusion). Note that including CR physics makes the gas more unstable and
leads to prominent magnetic “Parker loops.” Image from Rodrigues et al.
(2016); reproduced with permission from ApJ.

This study was further extended to include self-gravity of the gas and it
corroborated earlier findings regarding the impact of diffusion on the growth
rates (Kuwabara & Ko, 2006, 2020). This work also demonstrated that CRs can
effectively oppose gravity and suppress the Jeans instability. For low
parallel diffusion, large CR gradients can be established that prevent the gas
from collapsing, which leads to dense filaments aligned with the local
direction of the magnetic field. In the opposite case of fast diffusion, the
gas slides along the magnetic field and forms filaments perpendicular to it.
Including CR transport via streaming further destabilizes the system. By
invoking similarities to overstability of acoustic waves driven by CR
streaming (Drury & Falle, 1986; Begelman & Zweibel, 1994), and considering
non-radiative simulations of the Parker instability with streaming, Heintz &
Zweibel (2018) identified the CR streaming instability heating as an important
destabilizing process. They noted that capturing the instability in this case
requires (volume-filling) numerical resolution better than 100 pc, which is
higher than in the case of classic Parker instability. However, since the
instability growth rates can be comparable for both the diffusion and
streaming cases (Heintz et al., 2020), the destabilizing impact of CR
streaming heating and the exact nature of transport may be secondary for the
development of the Parker instability, i.e., the instability can develop more
easily as long as some form of transport is present that allows CRs to escape
from the compressing gas regions.
The ISM may be even more unstable to the Parker instability in the presence of
efficient anisotropic thermal conduction. The origin of this instability can
be understood by considering a horizontal magnetic field in a gravitationally
stratified medium in which temperature increases in the direction of gravity.
As the field is perturbed in the vertical direction, the heat begins to flow
along the field in the upward direction. If the process is isobaric, then the
displaced parcel of gas becomes hotter and less dense, and thus more buoyant.
While this instability was considered by Balbus (2000, 2001) in the context of
high plasma $\beta$ medium, such as the ICM, it has been used to generalize
the Parker instability in the ISM by Dennis & Chandran (2009), who
demonstrated that the ISM is more prone to the Parker instability in the
presence of anisotropic conduction.
We note that the regime of applicability of the above models is limited by the
fact that they neglect driving of supersonic turbulence by SNe and
superbubbles in the ISM. In the environments where this driving is important,
slow buoyant motions are subdominant compared to fast advection, and the
Parker instability is unlikely to operate. However, in the coronal hot phase
regions of the ISM, the gas may be Parker unstable. The outcome will therefore
depend on a combination of uncertain factors such as the ratio of the time
between superbubble explosions and the instability growth timescale, spatial
separation between superbubble locations, and the properties (temperature,
velocity dispersion) of the phase ISM relatively unaffected by SN activity. If
these complications can be neglected, in the saturated state of the Parker
instability, large magnetic loops extending to significant heights above the
disk midplane can be formed. These loops could transport both the magnetic
fields and CRs away from the disk, possibly leading to CR-dominated galactic
halos. These processes are closely connected to the issue of galactic wind
launching that we discuss in the next section.

####  3.2.4 Wind launching and simulations of disk patches and stratified
boxes

In this section we begin the discussion of the impact of CRs on the scales
intermediate between the ISM scales responsible for star formation in
molecular clouds and the large scales of the CGM. Specifically, we discuss
what role CRs can play in shaping the smooth transition between the disk-like
distribution of the ISM and the more spherically distributed CGM surrounding
the galaxy. One of the key processes shaping this transition is the process of
galactic fountain and wind launching. We start by discussing the basic
mechanisms behind CR wind driving and the crucial role that CR transport plays
in this process. We then examine the consequences of the physics of wind
launching for star formation, mass and energy loading in the wind, and the
phase space distribution of the gas in the ISM and the wind immediately above
the disk.

##### CR Eddington limit and CR staircases.

As discussed above, CRs can exchange momentum with the gas via pitch angle
scattering off of Alfvén waves and thus accelerate the gas. In a pioneering
work, Socrates et al. (2008) suggested that if the CR energy generation rate
in a galaxy (i.e., CR luminosity) exceeds a critical value, then the gas can
be unbound. This argument is analogous to the Eddington luminosity argument
with the difference being that photons are replaced by different relativistic
particles, i.e., by CRs. The critical CR luminosity $L_{\rm Edd,cr}$ above
which CR forces exceed gravitational forces can be easily obtained in the
diffusion limit (cf. Huang & Davis, 2022, similar expressions to the ones
presented below can be obtained for the streaming-dominated case). In this
case, balancing CR flux with the critical flux needed to offset gravity, and
assuming spherical symmetry, we have

$\frac{L_{\rm Edd,cr}}{4\pi
r^{2}}\boldsymbol{\hat{r}}=-\frac{\kappa}{\gamma_{\rm
cr}-1}\boldsymbol{\nabla}P_{\rm cr}=-\frac{\kappa}{\gamma_{\rm
cr}-1}\rho\boldsymbol{g},$ (64)

where $\kappa$ is the CR diffusion coefficient, $P_{\rm cr}$ is CR pressure,
$\rho$ is the gas density, and $\boldsymbol{g}$ is the gravitational
acceleration. For a point source of gravity, the Eddington luminosity is

$L_{\rm Edd,cr}=\frac{4\pi GM\kappa\rho}{\gamma_{\rm cr}-1}.$ (65)

Analogous expression can be derived for radiation, in which case
$\kappa\rho/(\gamma_{\rm cr}-1)=m_{\rm p}c/\sigma_{\rm T}$, for Thomson
scattering. If $\kappa$ is constant (as is often assumed in simple CR
propagation models), then equation (65) suggests that lower density gas should
be more easily accelerated. This is in contrast to the radiation case, whereby
photons tend to transfer little momentum to the gas as they escape the
atmosphere through low density channels. This comparison can be put on a more
quantitative footing by expressing the Eddington limit in terms of mean free
path $\lambda_{\rm mfp}$ of CRs or photons. In both cases, respective $L_{\rm
Edd}\propto\rho\lambda_{\rm mfp}$, and consequently the ratio of CR and
radiative Eddington limits is $L_{\rm Edd,cr}/L_{\rm
Edd,\gamma}\sim\lambda_{\rm mfp,cr}/\lambda_{\mathrm{mfp,}\gamma}$. Based on
the spallation argument presented in the Introduction (and detailed in Section
4.1.1), the relatively long residence time $\tau_{\rm esc}\sim 3\times 10^{7}$
yr of CRs near the disk midplane (compared to the light crossing time $\sim
H/c\sim 2\times 10^{4}$ yr; assuming a height of the CR scattering halo of
$H\sim 5$ kpc) translates to a short mean free path of CRs
$\lambda_{\mathrm{mfp,cr}}\sim 3H^{2}/c\tau_{\rm esc}\sim 8$ pc. Since this
mean free path is much shorter than the Thomson mean free path for photons
($\sim 500$ kpc for typical ISM density of 1 cm-3), the CR Eddington limit is
much lower than the usual limit for radiation. Above this CR limit, the gas
cannot exist in hydrostatic equilibrium and a wind should develop. Since the
production rate of CRs is closely linked to the star formation rate, this
limit corresponds to a critical star formation rate above which winds should
be launched. We note that, just as in the case of the classic derivation of
the Eddington limit for radiation, the CR Eddington limit is not strict due to
a number of simplifying assumptions such as spherical symmetry, neglect of
inelastic energy losses, spatial variation in opacity that can flatten
pressure profiles, etc.
Whether the wind is actually launched depends on a variety of factors. Using
one-dimensional steady state calculations, Heintz & Zweibel (2022) generalized
the above argument to include gravitational potential from a spherical mass
distribution at the galactic center and the dark matter in the halo. By
resorting to the solar wind analogy (Parker, 1958), where the solar wind is
expected to be launched when the asymptotic pressure in the solar corona
exceeds that of the interstellar space, they compared asymptotic CR pressure
to the CGM pressure to determine the conditions for CR wind launching. By
examining the role of different CR transport mechanisms, Heintz & Zweibel
(2022) found that solutions were strongly dependent on the physics of CR
transport and that winds were only likely to develop in the presence of
streaming. Moreover, while they confirmed that the Eddington limit does exist,
it was not likely to be exceeded in realistic galaxies unless gas densities
were reduced below the observed average values.
Nevertheless, in realistic multiphase environments, it is plausible that CRs
could be transported to lower density regions and drive the winds for a wider
range of parameters. More broadly, the work of Heintz & Zweibel (2022)
therefore underscores the need for considering the role of multidimensional
and time-dependent effects and the rich physics of CR transport in launching
the winds. The work of Huang & Davis (2022) makes a step in that direction by
performing two-dimensional MHD simulations of stratified atmospheres including
radiative cooling and CR transport. In their simulations, CRs are injected at
the midplane of the atmosphere and propagate either via diffusion or
streaming. They find that the Eddington limit can be reached for a broad range
of star formation environments and that winds can be launched when CR
transport is governed by diffusion or streaming. In simulations with
diffusion, streaming heating is absent and all of the CR energy can be used to
increase the kinetic energy of the gas, thus making the gas acceleration more
efficient. In simulations with streaming, they observed growth and saturation
of the CR acoustic instability. This instability has the potential to
considerably alter the picture of CR wind acceleration.
In general, the CR acoustic instability arises as a result of a relative
displacement due to streaming or diffusion of the CR fluid with respect to the
gas during the propagation of a sound wave. In the diffusion-dominated regime,
this instability only occurs in rare cases when the CR pressure scaleheight is
shorter than the CR diffusion length (Drury & Falle, 1986) and, consequently,
it is not seen in the simulations of Huang & Davis (2022). In the streaming-
dominated regime, the instability can occur even in the presence of weak CR
background gradient (Begelman & Zweibel, 1994). In this case, the relative
shift between CRs and the gas, combined with the CR streaming heating, leads
to effective forcing of the gas that may either amplify the oscillations in
the gas density, if the forcing is in phase with the waves, or damp them in
the opposite case. As demonstrated by Tsung et al. (2022) via one-dimensional
time-dependent simulations, as the waves steepen, they form sequences of
shocks in the gas. The density jumps associated with these shocks may result
in the formation of CR bottlenecks. This implies that the CR distribution
develops plateaus interspersed by sharp jumps at the locations of the shocks.
Consequently, CRs are dynamically decoupled throughout most of the volume but
exert strong forces at the locations of the jumps in CR pressure. While the
gas is not heated by CRs in the plateau regions, it is subject to intense
heating at the shock locations. Furthermore, since radiative cooling is not
offset by CR heating at the plateaus, the gas can be severely out of thermal
equilibrium, which can facilitate the formation of additional shocks. Overall,
the transfer of energy and momentum from CRs to the gas is expected to be
highly spatially and temporally dependent. The formation of the bottlenecks
and CR decoupling at stair-like plateaus in the distribution of CRs was also
observed by Huang & Davis (2022). They also determined that the CR heating of
the gas is dominated by that occurring near the shocks rather than the more
volume-filling regions characterized by shallower CR gradients. Despite these
intricacies of CR interactions with the gas, the numerical experiments by
Huang & Davis (2022) show that CR and the gas are nevertheless well coupled
dynamically – the efficiency of momentum transfer was found to be nearly
perfect irrespectively of whether CR transport occurred via diffusion or
streaming. Overall, these simulations suggest that CR fluxes can achieve
super-Eddington levels. The impact of CR staircase formation and CR transport
on launching highly variable outflows is further discussed in Section 3.2.1.
The level of sophistication of the models presented above can be increased by
considering a number of factors. First, the ability of bottlenecks to impede
CR transport is likely to be reduced in three-dimensional simulations compared
to those performed in one and two dimensions because the volume-filling factor
of the plateaus would be significantly reduced. Second, the severity of the
impact of the bottlenecks can be further affected by considering self-
consistent simulations of thermal instability in which the magnetic field is
allowed to increase with the gas density. Allowing for this effect would
reduce the frequency with which the bottlenecks occur because Alfvén speed
minima would be less deep. Third, the level of decoupling of CRs from the gas
can be further affected by, e.g., ion-neutral damping. Further work is needed
to ascertain the impact of these factors on CR propagation and triggering of
outflows from multiphase gravitationally stratified atmospheres. All of these
issues fall under the broad category of the impact of CR transport on wind
launching, which is the topic we discuss in the next section. We return to the
issue of the CR Eddington limit in Section 3.2.4, where we more closely
examine the link between this limit and star formation.

##### Impact of CR transport on wind launching.

In this section, we discuss in more detail the role of CR transport in
triggering the winds in gravitationally stratified atmospheres. We begin by
considering the dependence of transport on environment and then discuss how CR
pressure builds up in the disk and launches winds. From the perspective of
wind launching, the key reason why CR transport is important is that it
regulates the magnitude of the outward-pointing gas acceleration. It does so
either via (i) shaping the CR pressure gradient and/or (ii) by allowing for
the CR momentum to be deposited in lower density gas that is easier to
accelerate away from the disk. These issues are central to the question of
whether an adequate CR pressure gradient can be built up to accelerate the gas
and launch the wind.
As mentioned in Section 2.3.6, CRs can either scatter on Alfvén waves
generated via the streaming instability or on extrinsically generated waves
that are due to turbulent cascade. For CRs with energies less than $\sim 10$
GeV, which contribute most of the CR pressure support, self-confinement due to
the streaming instability is expected to dominate over extrinsic confinement
(Evoli et al., 2018). The effective drift speed of CRs with respect to the
bulk flow of the gas, either dominated by streaming or diffusion, can then be
determined by comparing, under the assumption of steady state, the streaming
instability growth rate of the scattering waves to their damping rates. The
damping can be due to ion-neutral damping, turbulent damping, linear or non-
linear Landau damping. Depending on the environment, the effective speed of CR
transport will be dominated by advection, streaming, or diffusion. By
neglecting CR-gas dynamical coupling, Armillotta et al. (2021) focused on the
physics of transport and considered CR propagation in the presence of non-
linear Landau damping and ion-neutral damping by post-processing TIGRESS
simulations (Three phase ISM in Galaxies Resolving Evolution with Star
formation and Supernova feedback, Kim & Ostriker, 2017, 2018) of galactic disk
patches typical of the solar neighborhood. They found that transport is highly
spatially variable and the scattering coefficient varies over four orders of
magnitude. Near the disk plane, where most of the mass of the ISM is in
neutral form, ion-neutral damping is extremely important, and consequently,
diffusion is much faster than advection. While in this environment streaming
is also expected to be very fast due to very low ionization, and thus large
ion Alfvén speed, diffusion is nevertheless more important than streaming for
$\sim$ GeV CRs. The coupling of CRs to the gas is restored in the layers
adjacent to the disk midplane, where ionization is much higher so that non-
linear Landau damping is the dominant damping process, significantly slowing
down CR transport.
Due to the sudden drop in the effective CR transport speed in that region, CRs
can be effectively trapped in the cold gas. As a result, their actual
transport speed may depend more strongly on the conditions in the volume
filling phase above the disk rather than in the region where ion-neutral
damping leads to very fast transport (Hopkins et al., 2021c; Bustard &
Zweibel, 2021; Huang & Davis, 2022; Armillotta et al., 2021). While the CR
pressure can build up near the midplane, the actual impact of CRs on the gas
acceleration depends on the gradient of pressure rather than its magnitude. In
the cold gas, where transport is very fast, pressure gradients along the
magnetic field will be very small and the ability of CR to accelerate the gas
will be diminished (e.g., Armillotta et al., 2021; Bustard & Zweibel, 2021).
Consequently, the impact of CRs on suppressing star formation on very small
scales or on launching winds in such environments may be reduced. However,
above the disk midplane, or perpendicular to the local magnetic field
orientation, or at the interfaces between individual cold clouds and the
hotter ISM, the coupling between CR and the gas can be restored, and large CR
pressure gradients can be present.
Significant CR pressure gradients can affect the gas distribution and the
impact of these gradients can be larger when they act on less dense gas.
Figure 21 illustrates schematically the impact of CR transport on triggering
the development of a galactic outflow (Salem & Bryan, 2014). The blue and
yellow lines correspond to hypothetical gas density and CR energy density
profiles, respectively. Gravitational field is in

Figure 21: Impact of CR transport on wind launching. The blue and yellow lines
correspond to hypothetical gas density and CR energy density profiles,
respectively. The top panel shows these distributions early in the evolution,
and the bottom panel corresponds to a time when CR got transported away from
the disk midplane. Image from Salem & Bryan (2014) reproduced with permission
from MNRAS.

the horizontal direction and points toward the density peak. The top panel
shows the gas and CR distributions in the initial state. In the absence of CR
transport relative to the gas such as diffusion and streaming, CR pressure
forces can thicken the disk (e.g., Girichidis et al., 2016b; Simpson et al.,
2016). When active CR transport is present, CRs can diffuse or stream away
from the midplane and their pressure gradients can begin to act on less dense
gas at larger heights above and below the disk. This effectively breaks the
global hydrostatic balance212121In this argument, we neglect turbulent
pressure and velocities that, by definition, are not hydrostatic. and allows
the gas to begin to flow out of the potential well. In the context of zoom-in
simulations of galactic disk patches, this was seen in the simulations of
Girichidis et al. (2016b) and Simpson et al. (2016), who considered
anisotropic CR diffusion with a constant coefficient and a detailed chemical
network that followed the abundances of key ions and molecules. In terms of
the treatment of CR transport, this was extended by Farber et al. (2018), who
considered a spatially-dependent anisotropic transport prescription in which
diffusion was boosted in the low-temperature ISM. This approach was motivated
by the desire to emulate the effects of ion-neutral damping combined with the
increased Alfvén speed in the cold gas. They demonstrated that despite the
effective CR decoupling from the gas in the low ionization environment, CRs
can still drive strong winds. While this work was performed under the
assumption that diffusion dominates transport, in a similar setup, Holguin et
al. (2019) showed that the results carry over to the case of CR streaming.
They showed that streaming regulated by turbulent damping dependent on local
ISM magnetization leads to effectively super-Alfvénic CR drifts. Stronger
damping resulted in more extended CR and gas distributions and stronger winds.
Due to their ability to remove gas from galaxies, CRs have a potentially
significant bearing on star formation, mass and energy loading of the
outflows, and the gas phase-space distribution in galactic winds. Below we
continue the discussion of these topics and focus on role of CRs in star
formation feedback.

##### CR feedback and star formation.

Stellar feedback plays a critical role in regulating star formation. In order
to place the impact of CRs in the appropriate context, we first list well
established stellar feedback mechanisms. The most important contributions to
the feedback come from massive stars in their late evolutionary stages. These
stars interact with the surrounding ISM via: (i) radiation pressure associated
with ultra-violet and IR photons interacting with dust (e.g., Krumholz &
Matzner, 2009; Fall et al., 2010; Murray et al., 2010; Raskutti et al., 2016;
Thompson & Krumholz, 2016), (ii) photoionization heating (e.g., Geen et al.,
2016; Gavagnin et al., 2017; Dale et al., 2012; Walch et al., 2013; Dale et
al., 2014; Haid et al., 2019; Kim et al., 2021), (iii) stellar winds (e.g.,
Castor et al., 1975; Weaver et al., 1977; Lancaster et al., 2021), and (iv)
SNR shock waves (e.g., Rogers & Pittard, 2013; Calura et al., 2015).
Each of these processes comes with its limitations in terms of the ability to
regulate star formation. For example, the question of the importance of
radiation pressure feedback is actively debated. While it has been argued that
radiation pressure can be very effective in regulating star formation in high
density environments (Krumholz & Matzner, 2009; Murray et al., 2010; Kim et
al., 2018b), claims to the contrary have also been presented (Rosdahl &
Teyssier, 2015; Rosdahl et al., 2015; Tsang & Milosavljević, 2018), arguing
that radiation can escape low-density channels without depositing a large
fraction of the photon momentum. In particular, Menon et al. (2022) suggested
that IR radiation pressure should have only very mild impact on star formation
and is unlikely to drive galactic winds. Similarly, Reissl et al. (2018)
demonstrated that UV and optical photons from young stars do not impart
significant momentum to the ambient gas as these photons are efficiently
absorbed and re-emitted in the far-IR, where the gas is optically thin.
Heating resulting from photoionization of the surrounding gas due to ultra-
violet photon absorption is also an important feedback mechanism, but its
importance is dependent the environment. For example, Fukushima & Yajima
(2021) argue that photoionization feedback is inefficient in giant molecular
clouds with surface densities above 100 M⊙ pc-2. Similarly, Haid et al.
(2018), demonstrated that photoionizing radiation is important in the cold
neutral medium, but it is subdominant compared to stellar winds in the warm
ionized medium. Stellar winds can help to disperse gas clouds surrounding star
formation sites (e.g., Castor et al., 1975; Wünsch et al., 2008, 2011; Toalá &
Arthur, 2011; Dale et al., 2012; Rogers & Pittard, 2013). However, the
efficiency of this form of feedback can be limited by the tendency of the hot
wind gas to escape through low-density channels and by fast cooling at the
mixing layers between the hot and cold gas (Rogers & Pittard, 2013; Rosen et
al., 2014; Lancaster et al., 2021). The last of the feedback processes
mentioned above involves deposition of the momentum and energy of expanding
SNR shocks in the ISM. SN are the most violent of these feedback processes.
While energetically very important, SN explosions are delayed with respect to
the onset of star formation and other processes can disperse GMCs before SN
feedback can act (Fall et al., 2010; Geen et al., 2016; Grudić et al., 2022).
Nevertheless, SN have a dramatic impact on driving turbulence in the ISM
(McKee & Ostriker, 1977; Mac Low & Klessen, 2004; Padoan et al., 2016; Kim &
Ostriker, 2017), especially, if the SNe are clustered in space (Girichidis et
al., 2016a).
CRs injected by SNe add a new “dimension” to the feedback problem and possess
a number of appealing features that help to maximize the impact of SN
feedback. In most numerical simulations of SN feedback, the resolution is
insufficient to resolve individual SNe and the SN energy is quickly radiated
away. This catastrophic energy loss is known as the overcooling problem (e.g.,
Walch et al., 2015; Naab & Ostriker, 2017). From the numerical point of view,
several approaches have been proposed to avoid this problem including:
temporarily switching off radiative cooling (e.g., Stinson et al., 2006;
Teyssier et al., 2013), distributing SN energy injection sites stochastically
(e.g., Dalla Vecchia & Schaye, 2012; Simpson et al., 2016), including
injection of momentum in addition to thermal energy (e.g., Martizzi et al.,
2015; Keller et al., 2022), and including natal kicks (e.g., Steinwandel et
al., 2022). As discussed in Section 3.2.2, CRs are injected by SNe and the
pressure inside SNRs may be dominated by CRs later in their evolution because
the adiabatic index of CRs is lower than that of the thermal gas. Moreover, in
situ injection of CRs at the SN shock fronts tends to result in larger
momentum transfer to the ambient ISM (Diesing & Caprioli, 2018). This effect,
combined with the relatively slower CR energy losses compared to radiative
losses of the thermal gas (see e.g., Figure 5 in Jubelgas et al., 2008), makes
CRs an attractive element of the feedback process in which the issues of
overcooling could be mitigated.
The impact of CR feedback can be amplified in superbubble environments, where
multiple shocks can accelerate CRs. Because most massive stars form in giant
molecular clouds and live for relatively short times before they explode as
SNe, the feedback due to SNe is expected to be strongly clustered in space and
time. Multiple SN explosions then lead to the formation of hot and low density
superbubbles. Consequently, rather than expending most of the SN energy in
high density environments, where catastrophic radiative cooling losses would
quickly radiate the energy away (Gatto et al., 2015), the SN energy can be
injected in lower density environment of overlapping SNRs (e.g., Mac Low &
McCray, 1988), where the effect of SN explosions can be sustained over longer
periods. Indeed, the impact of SN feedback is enhanced if they detonate in
low-density environments both in terms of their momentum transfer to the
surrounding ISM (Walch et al., 2015; Martizzi et al., 2015; Gentry et al.,
2017) and the ability to drive winds (Creasey et al., 2013; Fielding et al.,
2017b; Fielding et al., 2018). The dynamical impact of CRs could be further
increased in clustered SNe as the re-acceleration of pre-existing CR
populations could result in an elevated CR acceleration efficiency (e.g.,
Caprioli et al., 2018; Vieu et al., 2022).

##### CRs and energy and mass loading of the winds.

As discussed above, CRs can play an important role in launching galactic winds
and regulating star formation. The efficiency of this feedback can be
quantified in terms of either energy $\eta_{\mathrm{e}}$ or mass
$\eta_{\mathrm{m}}$ loading factors. The energy and mass loading factors are
typically defined as the ratio of the kinetic plus thermal luminosity in the
outflow to the energy injection rate from SNe, and the ratio of the mass
outflow rate to the star formation rate, respectively. In simulations, star
formation and outflow rates are typically measured simultaneously. However,
the instantaneous star formation rate does not directly control the outflow at
some height, because the outflow was launched by momentum deposition
associated with earlier star formation events. The situation is even more
complicated in observations. The star formation rate is indirectly inferred
through various estimators such as the rest-frame UV or FIR emission, which
probe an instantaneous star formation event, while the outflow rate measured
at a finite height above the disk is typically inferred from spectral line
shifts excited some time after the corresponding star formation event that
triggered the wind near the disk midplane. These are difficult to localize in
space and hence, cannot be unambiguously attributed to outflowing gas at a
specific height above the disk, which complicates direct comparison to
simulations.

Observations suggest that the mass loading factor $\eta_{\mathrm{m}}$ varies
widely and ranges from $\sim$ 0.01 to $\sim$ 10 (e.g., Veilleux et al., 2005)
Values broadly falling within this range have been seen in numerical
simulations of small stratified boxes. For example, strong outflows with mass
loading factor values up to 10 were reported from small-box SILCC simulations
(SImulating the Life Cycle of molecular Clouds, Girichidis et al., 2016a) and
values ranging from 0.01 to 10 were seen in TIGRESS simulations (Kim &
Ostriker, 2018). Furthermore, Li et al. (2017a) found that $\eta_{\mathrm{m}}$
decreases as a function of galactic surface density with 50
$\mathrm{M}_{\odot}~{}$pc-2 marking the transition where
$\eta_{\mathrm{m}}=1$. However, particular mass loading values depend on the
distance from the galactic midplane where the measurements of gas flows are
made and on a particular gas phase. An insightful unification of these type of
small-box simulations in terms of the wind loading factors has been presented
by Li & Bryan (2020). They conducted a systematic literature search and showed
that hot outflows are much more powerful than cool outflows. This analysis
revealed a unifying trend that the specific energy density of the hot phase in
the wind is $\propto\eta_{\mathrm{e}}/\eta_{\mathrm{m}}$ over four orders of
magnitude in star formation surface density. However, it has been argued
(e.g., Martizzi et al., 2016) that simulations of Cartesian boxes with a uni-
directional gravitational field may not fully capture the outflow physics as
they neglect the correct global geometry and galactic gravitational potential.
This omission tends to suppress the mass loading factors and reduces wind
speeds. We return to this question in Section 3.2.5, where we discuss the role
of CRs in global models of isolated galaxies.
How do CRs change this picture? The fundamental difference between the SN-
powered winds and those also assisted by the injection of CRs is that the
latter can experience continuous acceleration due to CR pressure gradients
while the former are more episodic in nature. Consequently, a generic feature
of CR-driven outflows is that their speeds increase as a function of the
distance from the disk, which is seen not only in stratified box simulations
(e.g., Girichidis et al., 2018) but also in one-dimensional models (e.g.,
Dorfi & Breitschwerdt, 2012) discussed in Section 3.2.1 and global models
(e.g., Pakmor et al., 2016c) discussed in Section 3.2.5. Importantly, this
offers a natural way to explain observational trends (Steidel et al., 2010).
However, other models not invoking CRs can also explain these trends (see
Section 4.4.5).
The impact of CRs on wind driving is best illustrated via controlled
experiments designed to isolate differences in results that arise solely due
to CR physics. For example, Simpson et al. (2016)

Figure 22: Gas density distributions (top row) in the plane perpendicular to a
galactic disk and the corresponding phase space gas distributions (bottom
row). From left to right shown are the cases corresponding to no CRs, CRs with
advection only, CRs with advection and anisotropic diffusion, and random
placement of SNe without CR physics. Dashed horizontal lines indicate the
height that contains half the gas mass in the initial conditions. Image from
Simpson et al. (2016) reproduced with permission from ApJ.

performed small stratified box experiments in which SNe and CRs were injected
at density peaks and CRs were allowed to diffuse out either isotropically or
along the magnetic fields. These cases were then compared to non-CR runs in
which SNe were placed randomly in space near the disk midplane (see Fig. 22).
These experiments revealed that in the former case, despite the fact that SN
feedback coincided with the densest regions where the radiative losses were
the strongest, the mass loading of the wind was comparable to that
corresponding to the latter case in which SNe were injected at random
locations and CRs were absent. However, unlike the peak-driving cases (with
and without CRs) that produced a self-regulated ISM, the random-placement
model resulted in a unphysical cooling-induced formation of overdense ISM
clouds. In an independent study of non-CR SN feedback in stratified boxes,
Girichidis et al. (2016a) confirmed that indeed only when SNe were placed
randomly with respect to the density peaks could the SN energy be injected in
sufficiently low density regions that the excessive radiative losses could be
avoided and the winds launched with significant mass loading factors.
Furthermore, by comparing non-CR runs to those that incorporate CR pressure
and diffusion, Girichidis et al. (2016b) showed that including CR physics
results in loading factors around unity, which were much larger compared to
those seen in non-CR simulations. In the CR simulations, the dominant source
of acceleration above a height of $\sim 1$ kpc was CR pressure gradients
(Girichidis et al., 2018). Moreover, emulating the effect of ion-neutral
decoupling of CRs from the gas in low temperature ISM regions reduces the mass
loading factors of the winds only moderately (Farber et al., 2018). Similarly,
Holguin et al. (2019) demonstrated that the increase in turbulent damping led
to faster CR streaming that helped to drive winds with progressively larger
mass outflow rates while allowing for increasingly larger star formation rates
(due to reduced CR pressure support in the disk). Overall, these small-box
studies demonstrate that CRs can play an important dynamical role in launching
galactic winds.

##### Phase space structure.

The observed galactic outflows are multiphase (e.g., molecular: Leroy et al.
2015; Chisholm & Matsushita 2016; neutral: Heckman et al. 2000; Rupke et al.
2005; Contursi et al. 2013; ionized: Steidel et al. 2010; Erb et al. 2012;
Heckman et al. 2015). As discussed in Section 3.2.2, the origin of the
multiphase nature of these outflows is uncertain, but it is either due to (i)
in situ formation of the cold phase via local thermal instability in hot
outflows or (ii) acceleration of multiphase gas from the vicinity of the
galactic disks.
Given that the outflows are launched from the disk that is occupied by
multiphase ISM, it is not surprising that the winds themselves could be
multiphase. The phase-space distribution of the gas in the midplane and in the
outflow regions above it depend on a number of factors including the details
of stellar feedback briefly discussed in Section 3.2.4 above. In this regard,
CRs in particular have the ability to significantly alter these phase space
distributions. For example, by systematically varying the assumptions about
the environments in which SNe explode, Simpson et al. (2023) showed in
stratified box experiments that the ISM either evolves to a state where warm
gas is volume filling, when SN are preferentially placed at the locations of
gas density peaks, or to a “thermal runaway” state, which occurs when SNe are
placed randomly in space so that their energy deposition cannot prevent the
mass-carrying cold and warm phases from collapsing to very high densities and
the hot gas becomes volume filling. Interestingly, dynamically important CR
pressure and CR diffusion have the ability to transform one solution to
another. In a related work, Commerçon et al. (2019) considered CR propagation
in a bi-stable ISM and discovered a critical diffusion level that separates
the regime where CRs can propagate freely from the one in which CRs are
trapped to produce strong CR pressure gradients that suppress thermal
instability and thus, affect star formation rates. Interestingly, some of the
most complete models of the ISM that in addition to including a myriad of
processes mentioned above (e.g., SN clustering, pre-SN feedback, radiation
pressure) also incorporate CR physics (Rathjen et al., 2021), result in the
best matches to the observed gas phase space structure, energy density
partition in the ISM, and star formation rates.
A generic feature of CR-driven winds is the tendency of the CR pressure to
become progressively more important compared to the kinetic, thermal, and
magnetic pressures as a function of the distance from the disk (e.g., Hanasz
et al. 2009a; Simpson et al. 2016; Girichidis et al. 2018; Armillotta et al.
2021; see also Section 3.2.5 below). One consequence of this dynamically
important CR pressure support is that thermal gas can be significantly cooler
because thermal pressure is replaced by CR pressure and different gas phases
no longer have to be in pressure equilibrium with each other. Indeed,
simulations demonstrate that CR driven outflows are denser and much cooler
compared to those that are thermally driven (e.g., Girichidis et al. 2016b,
2018; see also bottom row in Fig. 22). The addition of CRs also generally
tends to make the outflows slower and smoother in terms of density
distribution (cf. third and fourth panel in the top row of Fig. 22) while
reducing the volume filling factor of the hot gas phase provided half of the
SNe are placed randomly in the ISM (see figure 18 of Simpson et al., 2023),
though details of CR transport may affect these findings (Farber et al.,
2018). These overall trends may also help to explain the properties of the
warm ionized gas layers above the disks of the Milky Way and other disk
galaxies (Haffner et al., 2009). This warm gas is most likely not expelled
from the disk by very hot SN-driven thermal outflows, and cooler CR-driven
outflows may offer a solution to the problem of lifting the gas without
overheating it. While the warm ionized medium is most likely heated by
photoionization due to massive OB stars (Reynolds, 1990b, a), additional
heating mechanisms may be required to provide the necessary heating to explain
the [S II]/H$\alpha$ and [N II]/H$\alpha$ line intensity ratios (Reynolds et
al., 1999). Interestingly, this supplementary heating could be at least
partially provided by the streaming CRs (Wiener et al., 2013b) that are
already present in the wind.

##### Impact of CRs on star formation and winds in varying galactic
environments.

We now turn our attention to the role of CRs in shaping star formation across
a range of galactic environments. As discussed in Section 3.2.4, whether a
particular galaxy, or region within a galaxy, can launch a CR-driven wind
depends on the relative value of (i) the CR flux supplied by SNe compared to
(ii) the critical CR flux needed to unbind the gas (Socrates et al., 2008).
Under the assumption of negligible CR losses, the former is proportional to
the surface density of the star formation rate, $\dot{\Sigma}_{\star}$. The
latter increases with the strength of disk gravity and the speed of CR
transport, and both of these factors help to stabilize the disk against
disruption. In particular, for faster transport of CRs away from the disk,
larger CR fluxes are needed to unbind the gas in order to compensate for
weaker coupling of CRs to the gas. Depending on the details of CR transport
physics, this competition between outward-pointing CR pressure forces and
downward-pointing disk gravity leads to a relationship between the surface
density of the disk $\Sigma_{\rm gas}$ and the critical star formation rate
density $\dot{\Sigma_{\star}}$ above which the gas can be unbound by CR
forces. Predictions from such models can be compared to the classic Kennicutt-
Schmidt law (Kennicutt, 1989; de los Reyes & Kennicutt, 2019; Kennicutt & De
Los Reyes, 2021), i.e., the observed relationship between
$\dot{\Sigma_{\star}}$ and $\Sigma_{\rm gas}$. Based on this argument, CRs
could drive winds over a wide range of star formation environments (Huang &
Davis, 2022) as long as CR losses can be neglected. However, the impact of CR
losses is generally expected to be higher for larger $\Sigma_{\rm gas}$. There
is an emerging consensus that the observed star formation rate may be
sufficient to drive winds at least in lower $\Sigma_{\rm gas}$ environments
(Crocker et al., 2021a, b; Heintz & Zweibel, 2022), but at higher surface
densities the impact of CRs, and their pressure compared to the gas pressure,
may be reduced due to a combination of catastrophic hadronic losses (Lacki et
al., 2010; Lacki, 2011; Krumholz et al., 2020; Crocker et al., 2021b) and/or
fast transport away from the disk (Armillotta et al., 2021, 2022).
The relative importance of CR hadronic losses is closely tied to the details
of modeling of the physics of CR transport. Consequently, a particular shape
of the relationship between critical star formation rate needed to launch
winds and the surface density of the gas likely depends on this physics. At
high surface densities, galaxies may be good CR calorimeters (and powerful
$\gamma$-ray sources), meaning that most of the energy injected by SNe in the
form of CRs is radiated away. However, where exactly the transition to
calorimetry occurs as a function of $\Sigma_{\rm gas}$ depends on CR transport
physics and the ability to reliably model the small-scale porosity and
magnetic connectivity of the different phases of the ISM. Moreover, it has
been suggested that CR transport in the diffuse medium surrounding the cold
ISM needs to be very fast in order to limit hadronic losses and prevent
overproduction of $\gamma$-ray emission (Chan et al., 2019; Hopkins et al.,
2021c; Bustard & Zweibel, 2021) but claims to the contrary have also been
presented (Nuñez-Castiñeyra et al. 2022; Simpson et al. 2023; see also
Girichidis et al. 2018; Heintz & Zweibel 2022). Further detailed and fully
self-consistent modeling of these processes for a wide range of relevant ISM
parameters is needed to establish the conditions for the onset of CR wind
driving and the consistency of the models with the observed far-infrared
(FIR)–radio correlation (Condon 1992; see Section 4.3 for further discussion).

####  3.2.5 Role of cosmic rays in global models of isolated galaxies

Certain aspects of the feedback process cannot be fully captured in small-
scale ISM simulations (e.g., Girichidis et al., 2016a) or zoom-in simulations
of galactic patches (e.g., Simpson et al., 2016) discussed above. Most
notably, zoom-in simulations are, by construction, performed in planar
geometry and cannot include realistic galactic potentials. One consequence of
adopting planar geometry is that the wind solutions cannot transition from
subsonic to supersonic regime as is expected in the standard Chevalier & Clegg
(1985) galactic wind theory (Martizzi et al., 2016). As a result, outflows are
slower, typical mass loading factors appear to be substantially below those
seen in rapidly star-forming galaxies, and fountain flows are more likely to
develop

Figure 23: Slices of mass flux, thermal gas pressure, CR pressure, and the
ratio of the CR-to-thermal pressures in the plane perpendicular to the disk;
from a global simulation of a disk galaxy from Salem & Bryan (2014);
reproduced with permission from MNRAS.

(cf. zoom-in TIGRESS simulations of Kim & Ostriker 2018 and CGOLS global
simulations of Schneider et al. 2020). By contrast, the effect of the
“dilution” of the gravitational potential in global models of galactic winds
naturally predicts conical outflows and larger mass loading factors. Global
models can also be used to study various aspects of the interactions of the
outflows with the ambient medium, e.g., by including radiatively cooling and
collapsing gaseous halos, that source turbulence and drive a small-scale
magnetic dynamo (e.g., Pfrommer et al., 2022), and gas streams, that feed the
disk and interact with the galactic wind (Peschken et al., 2021), or
considering the interplay between winds and ram pressure associated with the
bulk motions of galaxies (e.g., Bustard et al., 2020; Farber et al., 2022). In
this section, we mostly focus on those aspects of CR feedback that are either
unique to the global models or that can be more faithfully modeled in the
context of global models of galaxies in non-cosmological settings. We discuss
several intricately inter-related processes such as the magnetic field
amplification, CR losses due to hadronic and Coulomb processes, and CR
transport in global galactic models.

##### Impact of CR transport on wind driving.

As argued above in the context of one-dimensional models and zoom-in
simulations of galactic patches, CR transport plays central role in shaping
the properties of galactic outflows. These outflow properties can differ
substantially depending on whether CRs are advected or actively transported
via diffusion and/or streaming, and how exactly these processes are modeled.
Despite these differences, some generic features emerge from these models
including cooler and slower winds (compared to simulations where winds are
driven by thermal feedback), potentially significant CR pressure support in
the halo, substantial suppression of star formation, and thickening of the
galactic disks.
Global three-dimensional simulations of galactic winds including CR diffusion
were performed by a number of authors (e.g., Hanasz et al., 2013; Booth et
al., 2013; Pakmor et al., 2016c; Ruszkowski et al., 2017a). Below we contrast
two examples of early simulations in order to isolate the impact of non-
adiabatic CR energy losses. The first example is the global three-dimensional
hydrodynamical simulations of galactic winds including isotropic CR diffusion
presented by Salem & Bryan (2014). These simulations demonstrated that strong
winds with mass loading factors exceeding unity can be driven by CR feedback
in Milky Way mass halos. Figure 23 shows slices of mass flux, thermal gas
pressure, CR pressure, and the ratio of the CR-to-thermal pressures (lower
right panel) in the plane perpendicular to the disk. It is evident from these
results that the outflow is conical and that CR pressure dominates over
thermal pressure in the halo. As a result of the additional pressure support
due to CRs, the disk is inflated and star formation is reduced. Since the CR
content near the disk midplane is regulated by the speed of CR transport, the
CR pressure support decreases with increasing CR diffusion coefficient
$\kappa$, which leads to the trend for star formation rate to increase with
$\kappa$. The second example is the extension of the above work by Booth et
al. (2013), who incorporated the effect of CR energy losses due to hadronic
Coulomb collisions in dense regions and studied the impact of CR feedback in a
dwarf and Milky Way-mass galaxies. Despite CR losses, the impact of CR
pressure is still significant in these simulations because CR diffusion allows
for the energy to be transferred from the densest regions to more tenuous gas
before CRs cool completely. As a result, the buildup of CRs provides a robust
and gentle wind acceleration mechanism. In agreement with observations, this
gentle and continuous acceleration of the gas results in slower winds compared
to those accelerated by SN feedback without CRs. In the case of dwarfs, the
boost in the mass loading factor in these simulations is an order of magnitude
larger than in the case of non-CR SN feedback. In agreement with predictions
from smaller scale simulations of zoom-in patches of galactic disks discussed
earlier, global simulations also predict significantly cooler winds compared
to the cases where winds are powered by SNe only. This effect is due to the
fact that a significant fraction of the energy density in the wind goes to the
CR component rather than the thermal component. This shift reduces the gas
temperature, which further increases the efficiency of radiative cooling as
the gas is subject to cooling closer to the peak of the cooling function.
The impact of CR streaming is similar to that of CR diffusion in the sense
that the presence of this CR transport process facilitates the wind launching.
While in the absence of any active CR transport mechanism the disk only puffs
up and star formation is severely reduced, when CR streaming (or diffusion) is
present the wind can develop. First global three-dimensional hydrodynamical
simulations that included CR streaming were performed by Uhlig et al. (2012).
In addition to transport via streaming, CRs were subject to hadronic and
Coulomb losses and cooling associated with the excitation of Alfvén waves
(though the process was approximated by assuming that streaming occurs at the
local sound speed in these purely hydrodynamical simulations). The inclusion
of streaming resulted in powerful and sustained winds with wind speeds and
mass loading factors in agreement with simplified analytical estimates that
are based on the Bernoulli equation. Moreover, CR streaming losses led to the
excitation of Alfvén waves, and the associated heating of the halo. This work
was further extended by Ruszkowski et al. (2017a), who included streaming CR
in MHD simulations and also observed development of winds with substantial
mass loading factors and suppression of the star formation rates.
Despite these similarities, there are important differences between winds
enabled by streaming and diffusion of CRs. Generally speaking, the mass
loading factor of the wind launched by CR feedback depends on the speed of CR
transport. When the transport speed is vanishingly small, CRs are concentrated
near the disk midplane, reduce star formation, and are unable to lift the
dense gas. As the transport speed is increased, progressively larger CR
pressure gradients can build up in the lower density gas above the disk
midplane, which leads to the wind launching and the increase of the mass
loading factor with CR transport speed. In the other extreme, very fast
transport of CRs allows them to leave the midplane on a short timescale so
that they spend systematically less time accelerating this gas, which implies
a decreasing mass loading factor with increasing transport speed. While the
latter regime corresponds to the diffusion-dominated wind driving (e.g., Salem
& Bryan 2014), the former regime regime is applicable to winds driven by CR
streaming (e.g., Ruszkowski et al. 2017a). This argument can be made more
quantitative by considering characteristic transport speeds in the diffusion
and streaming cases. For typical conditions operating near star forming
regions (e.g., cloud size of $L=100$ pc, magnetic field of 5 $\mu$G, diffusion
coefficient $\kappa=10^{28}$ cm2 s-1, and gas density of $n=10$ cm-3), the
characteristic diffusion speed is $\kappa/L=320\,(\kappa/10^{28}~{}{\rm
cm}^{2}~{}{\rm s}^{-1})(L/100~{}\rm pc)^{-1}$km s-1, while the streaming speed
is the Alfvén speed $\sim 3.5\;(B/5~{}\mu{\rm G})(n/10\;{\rm cm}^{-3})^{-1/2}$
km s-1. This implies that the streaming speed is generally much smaller than
the diffusion speed. In fact, if CR streaming alone is to drive the winds,
then the streaming speed needs to be boosted by collisionless wave damping
mechanisms such as turbulent damping, non-linear Landau damping, or ion-
neutral damping, and/or by the increase in the Alfvén speed in the low
ionization regions, for the gas to be more easily accelerated (e.g., Farber et
al., 2018; Farcy et al., 2022). In a realistic situation, the effective CR
drift speed arises from a combination of streaming and diffusion and is
sufficiently fast to launch winds.
In addition to the above CR transport speed arguments, additional complexity
is introduced by the fact that, unlike in the extrinsic confinement picture
where CR transport is dominated by diffusion, in the self-confinement picture
the streaming CRs also excite Alfvén waves on which they scatter, and
collisionless damping of these waves effectively cools CRs and heats the
thermal gas, which subsequently quickly radiates this energy away.
Consequently, this energy is not efficiently used for driving the winds and
the winds are expected to be weaker as lower CR pressure gradients are
available to accelerate them (Wiener et al., 2017b). On the other hand, in the
purely diffusive regime, CRs only experience adiabatic changes of their energy
because spatial diffusion itself is energy conserving (neglecting hadronic and
Coulomb losses that are present in both transport regimes) and the effect of
CRs on wind driving can be overestimated. Thus, pure streaming and pure
diffusion cases likely bracket the range of outcomes in terms of wind driving.
While the exact CR transport mechanisms are still uncertain, for CR energies
below $\sim$ 100 GeV that contribute most significantly to the CR energy
density, CR scattering is likely to be dominated by self-generated turbulence
(e.g., Aloisio et al., 2015; Evoli et al., 2018), rather than scattering on
extrinsic turbulence (Yan & Lazarian, 2004, 2008). Hydrodynamical simulations
aimed at bracketing the range of possibilities in the scenarios where CR
transport is dominated by these two regimes (Wiener et al., 2017b) demonstrate
that diffusive transport leads to mass loss rates that are an order of
magnitude larger in comparison to those inferred in the streaming scenario.
These findings are broadly consistent with those obtained using MHD
simulations (Dashyan & Dubois, 2020). Specifically, while the simulations
employing isotropic and anisotropic diffusion with typical values of the
diffusion coefficient $\sim 3\times 10^{28}-10^{29}$ cm2s-1 result in
efficient wind driving, the simulations with streaming resemble those with
pure CR advection. Including streaming speed boost helps to launch the wind
but the outflow rates are nevertheless smaller than in the diffusion case and
smaller than those typically observed. This is consistent with the findings of
Chan et al. (2019) and Hopkins et al. (2020) performed in the cosmological
context (for cosmological simulations see Section 3.3) and including boosted
streaming, though significant differences exist in the typical levels of
magnetic field strength, and consequently the baseline Alfvén speed, between
these simulations and those of Dashyan & Dubois (2020) making these
comparisons uncertain (see below).
The relationship between the mass loading and transport speed discussed above
can be significantly altered in the presence of CR energy losses.
Specifically, the results of Dashyan & Dubois (2020) differ from those of
Ruszkowski et al. (2017a), who reported significant outflows in the case of
boosted streaming. A possible resolution to these differences lies in the fact
that Dashyan & Dubois (2020) include additional CR cooling due to hadronic and
Coulomb processes. Since the timescale for CR cooling can be shorter than that
associated with the CR escape timescale due to streaming at speeds moderately
boosted beyond the Alfvén speed, CR may lose energy before they can exert
sufficient pressure forces on the gas to launch the wind. However, larger
streaming boost due to damping processes could perhaps help to alleviate this
problem. Similarly, the trend for the mass loading to decrease with the value
of the diffusion coefficient (Salem & Bryan, 2014) discussed above is reversed
in Dashyan & Dubois (2020, see also ). In the latter work this dependence of
the mass loading on transport speed is also attributed to the fact that, for
higher $\kappa$, CRs spend less time in the denser regions where CR cooling
times are shorter.
Irrespectively of whether CR cooling losses are included, recent studies
report that faster transport speeds lead to more star formation (e.g.,
Jubelgas et al., 2008; Salem & Bryan, 2014; Ruszkowski et al., 2017a; Farber
et al., 2018; Chan et al., 2019). However, the corresponding wind mass loading
does not have to be monotonic with the transport speed. Based on global
hydrodynamical simulations that include diffusion and CR losses, Jacob et al.
(2018) report that while the mass loading decreases with the diffusion
coefficient for low halo masses, at larger masses the relationship is non-
monotonic. In the same regime of low halo masses, Farcy et al. (2022) find the
opposite trend but the comparison is made difficult by the fact that the
latter work also includes radiation pressure. Irrespectively of the halo mass,
total outflow rates increased with $\kappa$ up to a peak and then either
saturated or began to decrease near very high values of $\kappa\sim 3\times
10^{29}$ cm2 s-1, where CRs did not spend enough time accelerating the gas
before escaping from the galaxy. These trends for the mass outflow rate to
peak at a certain $\kappa$ were more evident for cold gas component. For the
most massive $10^{12}$ M⊙ halo, the increase in $\kappa$ led to larger total
outflow rates and hotter winds. On a qualitative level, this is consistent
with the trends seen in zoom-in simulations where the transport speed is
locally increased in cold gas due to the combination of ion-neutral friction
and the increase in the ion Alfvén speed (Farber et al., 2018).

##### Impact of non-equilibrium processes on CR transport.

Most models of CR feedback make simplifying assumptions regarding the steady
state CR spectra and transport. These assumptions are only justified when the
relevant timescales associated with the processes governing CR evolution are
shorter than other relevant timescales. Below we discuss the impact of
relaxing these assumptions on the evolution of outflows in global models of CR
feedback.
Despite the fact that the energies of CR particles span twelve orders of
magnitude and the CR transport speed depends on CR energy, most models of CR
feedback consider a single momentum-integrated CR fluid (in the “gray”
approximation). In this approximation, the evolution of CRs depends on
momentum-integrated CR cooling rates and transport coefficients. This
assumption is justified when the CR spectrum can reach a steady state on
sufficiently short timescales. The evolution of CR energy density at a
particular location depends on the competition between (i) gains due to CR
injection and advective, diffusive, and streaming transport into the region of
interest, and (ii) losses due to advective, diffusive, streaming transport out
of that region combined with hadronic and Coulomb cooling. Therefore, the
assumption of steady state, i.e., of balance between gains and losses, is
justified when the instantaneous timescales associated with these processes
are shorter than the timescale on which the overall net change of CR energy
occurs. The validity of the steady state assumption can be assessed a
posteriori by post-processing results obtained under the gray approximation.
This type of analysis reveals that, close to the disk, the steady state
approximation is valid, but it breaks down in very dynamical regions
especially above the disk (Werhahn et al., 2021a). This has important
implications of the modeling of galactic emission, which we defer to Section
4.3.
In order to relax the steady assumption, CR spectra need to be modeled on-the-
fly using multi-energy CR fluid approach coupled with MHD (Miniati, 2001; Yang
& Ruszkowski, 2017; Girichidis et al., 2019, 2022; Ogrodnik et al., 2021, see
also Section 2.5.2 for a detailed discussion). Irrespective of whether we
model the full CR spectrum or whether we treat them in the gray approximation,
star formation rates and galaxy morphologies are not significantly affected.
However, mass loading factors of galactic winds are smaller by up to a factor
of four in dwarf galaxies in the spectral CR model (Girichidis et al., 2023).
In addition, the spectral CR treatment produces substantially more cold gas in
fountain flows in Milky Way-mass galaxies than grey models. The limited
influence of spectral CR modeling on hydrodynamics can be attributed to the
dominant contribution of GeV CRs to CR pressure. Additionally, the proportion
of fast-diffusing high-energy CRs is too small to exert a significant impact
on fluid dynamics of galactic winds, particularly in larger galaxies. These
simulations show that a powerful CR-driven outflow in smaller galaxies (with
halo masses $M\lesssim 3\times 10^{11}~{}\mathrm{M}_{\odot}$) primarily
advects CRs with the gas. By contrast, in Milky Way-mass galaxies ($M\approx
1\times 10^{12}~{}\mathrm{M}_{\odot}$) the CR pressure gradient alone is not
able to power a sustained outflow, implying that CR transport is entirely
diffusive. Interestingly, the energy-weighted CR diffusion coefficient varies
spatially by up to two orders of magnitude. In the disk and in the outflow
region, it assumes the canonical values of GeV CRs of $1\times 10^{28}$ to
$3\times 10^{28}~{}\mathrm{cm}^{2}\mathrm{s}^{-1}$ and increases up to
$3\times 10^{29}~{}\mathrm{cm}^{2}\mathrm{s}^{-1}$ in the CGM, where CRs with
energies exceeding TeV dominate the spectrum. As a result, the CGM is
supported by the pressure from high-energy CRs, which enables a more
structured gas density and lower temperatures in comparison to the gray CR
models (Girichidis et al., 2023).
In general, both the CR flux, that determines how fast CRs are transported,
and the diffusion coefficient, that determines how fast CRs are scattered, are
time-dependent quantities. The steady state of the diffusion coefficient is
determined by the balance between the wave growth and damping, and is reached
on characteristic timescales corresponding to these processes. On the other
hand, the steady state of the CR flux is determined by both the CR scattering
rates and the hydrodynamic adjustments to the CR pressure. The timescale for
the establishment of the steady state of the CR flux is in general longer than
that needed to reach the steady state of the diffusion coefficient (Thomas et
al., 2023). Whether or not the steady state of CR flux can be reached depends
on the bulk CR speed. It is the CR pressure gradient along the magnetic field
that drives a CR flux to propagate along the field lines. If CRs move faster
than the Alfvén speed, they are able to excite gyroresonant Alfvén waves on
which they scatter. Sufficiently frequent scattering is a necessary condition
for reaching the steady state of the CR flux. Thus, in the opposite case,
i.e., when the bulk CR speed is smaller than the Alfvén speed, CRs damp the
waves, thereby reducing scattering and thus increasing the dynamical CR
timescales associated with the adjustments to the CR pressure so that a steady
state of the CR flux is not reached. In this case, the CR pressure gradient
along the magnetic field is small, which implies a small growth rate of the
gyroresonant instability. In addition, the CR pressure gradient can be
rearranged on MHD timescales in response to gas motions and magnetic forces,
which can also cause the CR flux to deviate from its steady state value.
Galactic wind models that include CR feedback and CR interactions with
gyroresonant Alfvén waves make it possible to self-consistently compute the CR
diffusion coefficient and CR transport speeds in the two-moment approximation
of CR transport (see Section 2.3.3). These models demonstrate that CR
transport does not reach steady state solutions describable solely in terms of
diffusion, streaming, or a combination thereof (Thomas et al., 2023).
Moreover, while CR diffusion indeed quickly reaches steady state in most
environments, there also exist regions where equilibrium diffusion coefficient
is not established. In these regions, the CR gradient is almost perpendicular
to the local magnetic field (and the CR pressure is almost constant along the
magnetic field). Consequently, the CR speed is smaller than the local Alfvén
speed and CRs are unable to excite gyroresonant Alfvén waves leading to very
high values of the diffusion coefficient at these “Alfvén wave dark regions.”
This is illustrated in Fig. 24 that shows cross sections perpendicular to the
disk of the Alfvén energy density (left panel) and CR diffusion coefficient
values (right panel).

Figure 24: Slices perpendicular to the disk of the Alfvén wave energy density
(left panel) and CR diffusion coeffcient (right panel). Notice the spatial
coincidence between the regions devoid of Alfvén wave energy and those
characterized by very large values of the diffusion coefficient. Image from
Thomas et al. (2023) reproduced with permission from MNRAS.

##### Halo mass dependence of galactic wind mass loading.

Observationally-determined wind mass loading factors can serve as a means of
constraining CR feedback models. Ultra-violet absorption line measurements
suggest that the mass loading depends weakly on galaxy mass, with studies
hinting at either no significant scaling with mass, consistent with
$\eta_{\mathrm{m}}\propto M_{\rm vir}^{-\alpha_{M}}$ with $\alpha_{M}\sim 1/3$
to $2/3$ within the large scatter (Heckman et al., 2015), or a weak decrease
with stellar mass as $M_{\star}^{-0.4}$ (Chisholm et al., 2017). This slope of
the relationship is consistent with the scaling of the mass loading factor
$\eta_{\mathrm{m}}$ with galaxy circular velocity
$\eta_{\mathrm{m}}\propto\varv_{\rm circ}^{-\alpha_{\varv}}$, where
$\alpha_{\varv}\sim 1$ to $2$, that is necessary to explain the low-mass end
of the stellar mass function of galaxies (e.g., Oppenheimer et al., 2010).
From the theoretical standpoint, the consensus that emerges from global
galactic models of stellar feedback aided by CR processes is that they can
drive stronger winds in smaller systems (Uhlig et al., 2012; Booth et al.,
2013; Jacob et al., 2018; Dashyan & Dubois, 2020; Farcy et al., 2022) but are
insufficient to drive winds in halos more massive than $\sim 10^{12}$ M⊙
(e.g., Fujita & Mac Low, 2018). However, there are considerable differences
between theoretical predictions for the dependence of the mass loading factor
on halo mass both in terms of the normalization and slope of the relation.
We now discuss various global (non-cosmological) galactic simulations of CR
feedback from the perspective of deriving galactic wind mass loading factors.
There are two popular setups for this type of simulations that potentially
differ in the resulting wind speeds and mass loading factors, especially if
measured at larger heights: (i) “disk models” that start from an initial
configuration of a stellar bulge and an exponential disk of stars and gas
embedded in a dark matter halo, and (ii) a “collapsing halo” setting that
starts from a rotating dark matter halo filled with gas, which starts to cool
at the beginning of the simulation. In the latter scenario, the densest gas at
the center cools fastest, conserves its specific angular momentum, and starts
to form a disk inside out. This disk finds itself embedded in a turbulent and
hot CGM that has been heated by the outwards propagating accretion shock
(e.g., Springel & Hernquist, 2003; Pfrommer et al., 2022). Moreover, the
infalling CGM provides inertia against which the galactic wind has to do
volume work. We caution that while many trends that emerge from these models
are broadly self-consistent, exact comparisons of models with different
physics is complicated by the fact that no unique definition of the mass
loading factor exist with various authors using custom definitions of this
quantity. Similarly, the adopted definitions of the mass loading factor are in
general not identical to those considered in observational studies.
Nevertheless, clear general trends emerge from the simulations and illustrate
the role of CR feedback in galactic wind driving.
Early hydrodynamical simulations of CR feedback using the “collapsing halo”
setup demonstrated that strong winds can develop in disk galaxies less massive
than $\sim\\!10^{11}$M⊙ (Uhlig et al., 2012). General patterns that transpire
from these simulations are that in dwarf galaxies CR feedback leads to
significant suppression of star formation (up to a factor of $\sim 5$) and
large mass loading factors (also up to a factor of $\sim 5$). While the wind
speed increases with the halo mass, the mass loading factor exhibits the
opposite trend becoming negligibly small at $\sim 10^{11}$ M⊙. However, even
in massive systems, where streaming did not lead to significant outflows,
strong fountain flows were present. In these simulations, CR transport is
modeled as streaming at the local sound speed and, similarly, the heating of
the gas associated with streaming is proportional to the CR pressure gradient
and local sound speed. Consequently, hotter regions are preferentially heated,
which leads to a significant sink of CR energy (cf., Wiener et al., 2017b)
that may limit the impact of CR wind driving due to fast radiative cooling of
the gas heated by CRs.

Figure 25: Projections of the vertical component of the gas velocity for halos
ordered according to mass (logarithm of mass is shown in the upper left corner
of each panel). Note that quasi-spherical CR-driven outflows exist in dwarf
galaxies; they transition to biconical outflows at intermediate mass galaxies
($\gtrsim 10^{11}$ M⊙) and fail to exist in halos more massive than that of
the Milky Way. Image from Jacob et al. (2018) reproduced with permission from
MNRAS.

Despite the above limitations, other models that do not rely on these
simplifying assumptions find broadly consistent trends albeit with significant
quantitative differences. For example, global hydrodynamical simulations using
the “disk model” that incorporate CRs via isotropic diffusion also lead to
wind velocities increasing with circular velocity of the halo Booth et al.
(2013). These simulations by construction do not include streaming losses as
CR transport is treated within the framework of the extrinsic CR confinement
model. In the simulations by Booth et al. (2013) that include CR feedback, the
wind velocity is substantially reduced (compared to the hydrodynamic-only
cases) and, unlike in the case of thermal SN feedback, it matches the slope
and normalization of the observed relation (Schwartz & Martin, 2004; Rupke et
al., 2005). Similarly to Uhlig et al. (2012), the mass loading in these
simulations was found to decrease with the halo mass. In dwarf galaxies, the
magnitude of the mass loading was significantly larger when CR feedback was
included as opposed to the case where only thermal SN feedback was considered.
Surprisingly, the mass loading factor of $\sim 0.3$ for the Milky Way-size
halo matched that reported by Salem & Bryan (2014), who considered similar
models but without hadronic and Coulomb losses. The above models were further
generalized to include magnetic field and anisotropic CR diffusion (Pakmor et
al., 2016c; Ruszkowski et al., 2017a; Pfrommer et al., 2017b; Jacob et al.,
2018). Projections of the vertical component of the gas velocity for halos
ordered according to mass are shown in Fig. 25 and reveal that only
intermediate mass galaxies show biconical outflows (Jacob et al., 2018).
Similarly to earlier studies, the mass loading factor in these simulations was
found to decrease with the halo mass with CRs unable to drive winds for halos
more massive than $\sim 10^{12}$ M⊙. However, the scaling of the mass loading
factor with the halo mass $M_{\mathrm{vir}}^{-\alpha_{M}}$, with
$\alpha_{M}\sim 1$ to 2 independently of whether CR diffusion was isotropic or
magnetic-field-aligned, was steeper than found in observations and other
galactic wind simulations (e.g., Muratov et al., 2015). However, despite the
above differences, the star conversion efficiency predicted by the model of
Jacob et al. (2018) agrees with the empirically derived one (Vogelsberger et
al., 2013).
In addition to the differences in the slope of the mass loading vs. halo mass
relation, Jacob et al. (2018) also predict different normalization of this
relation compared to other models. For example, while Salem & Bryan (2014) and
Booth et al. (2013) predict strong winds for halo mass $\sim\\!10^{12}$ M⊙,
Jacob et al. (2018) find negligible mass loading at this halo mass. This
difference is unlikely to be due to the neglect of CR cooling in Salem & Bryan
(2014) because both Jacob et al. (2018) and Booth et al. (2013) include CR
losses and yet Booth et al. (2013) and Salem & Bryan (2014) essentially find
the same mass loading factor value at this halo mass. For smaller mass halos
of $10^{10}$ to $10^{11}$ M⊙, MHD simulations suggest that while the mass
outflow rates increase significantly in the presence of diffusing and
streaming CRs, the mass loading factor falls short of that found in
observations (Dashyan & Dubois, 2020). Despite similar physics, mass loading
values reported by Dashyan & Dubois (2020) are an order of magnitude smaller
than those reported by Jacob et al. (2018) for this same mass range. However,
the exact comparisons of the mass loading factors inferred from the above
simulations is difficult due to (i) different mass loading factor definitions
adopted by various authors and, perhaps more importantly, (ii) because of
different ways of initializing the simulations either via a pre-existing “disk
models” (e.g., Booth et al., 2013) or a “collapsing halo” models (e.g., Jacob
et al., 2018)).
The first global simulations to simultaneously incorporate both CRs and
radiation feedback were presented by Farcy et al. (2022). In agreement with
earlier studies, CRs were found to substantially boost wind mass loading
factors and reduce the star formation rates with the efficiency decreasing
with halo mass. However, the combined CR and radiation feedback was still less
efficient than artificially boosted feedback designed to reproduce the
observed galaxy luminosity function in cosmological simulations (Rosdahl et
al., 2018).
Nevertheless, the true potential of CR feedback has not yet been fully
elucidated as current modelling of CR physics in simulations is incomplete.
Specifically, in addition to radiation pressure, several missing CR physics
factors that may boost the efficiency of feedback are required to increase the
realism of feedback modeling and need to be simultaneously incorporated in
future simulations. These factors include (i) self-consistent and non-
equilibrium computation of the CR diffusion coefficient and CR transport
speeds (Thomas et al., 2023), (ii) spectrally-resolved CR transport that not
only results in a qualitatively different wind driving (Girichidis et al.,
2022) but also enables detailed model comparisons to radio and $\gamma$-ray
data (Werhahn et al., 2023), (iii) CR transport in the multi-phase ISM with a
self-generated diffusion coefficient that responds to the prevalent wave
damping processes of the different phases of the ISM, (iv) in situ CR
acceleration boost at the location of SNRs (Diesing & Caprioli 2018; cf. Jana
et al. 2020, who do include in situ acceleration and find negligible effect of
CR feedback on wind driving in general, but consider only early stages of
feedback), and (v) cosmological effects (e.g., Chan et al. 2019; Hopkins et
al. 2020; Buck et al. 2020; see also Section 3.3).

##### Magnetic dynamo and CR transport.

Galactic winds and fountain flows discussed above play a critical role in
magnetic field amplification via the galactic dynamo process (Parker, 1992;
Hanasz & Lesch, 2000; Brandenburg & Subramanian, 2005; Shukurov et al., 2006).
Since CRs can launch the winds, they are an important ingredient in models
aiming to explain the magnetic field amplification. The galactic CR dynamo
model includes a number of interdependent elements such as differential
rotation, Coriolis force, outflow driven by CRs, and fast magnetic
reconnection. In the original Parker galactic dynamo model (Parker, 1992), CRs
are an essential ingredient of the model as they provide buoyancy to the gas.
Unlike in the case of buoyancy due to magnetic fields, CRs do not introduce
tension forces that tend to suppress the convective rise of the gas. Thus, the
total buoyant forces can exceed magnetic tension forces when the CR energy
density is sufficiently large compared to the equipartition energy density of
the magnetic field. These buoyant forces play an important role in the dynamo
process as they facilitate the production of the radial magnetic field
component from the azimuthal one. As the magnetic loops rise and twist, they
reconnect forming amplified magnetic structures that posses a radial magnetic
field component. For this process to operate efficiently, reconnection needs
to be fast in order to form large scale poloidal magnetic fields. Similarly,
the azimuthal component can be generated from the radial one via differential
rotation.
Several aspects of this process can be captured via local zoom-in simulations
of galactic disks. In particular, early work demonstrated that magnetic fields
can grow exponentially fast, reaching equipartition with the turbulent energy
on timescales comparable to 100 Myr even in the limit of vanishing explicit
magnetic resistivity (Hanasz et al., 2004b; Hanasz et al., 2004a).
Interestingly, these shearing box simulations showed that the mean radial
component of the magnetic field can be generated from the mean azimuthal one
and vice versa, consistent with the incoherent $\alpha-\Omega$ dynamo theory
(Vishniac & Brandenburg, 1997). However, other aspects of the problem are
better investigated using global simulations. For example, the fastest growing
modes of the Parker instability have characteristic sizes close to $\sim 1$
kpc (Giz & Shu, 1993; Rodrigues et al., 2016), which is comparable to the
sizes of zoom-in boxes (e.g., Simpson et al., 2016). This limitation may thus
lead to reduced instability growth rates. Limiting the geometry of the
computational domain to zoom-in boxes may also affect the distribution of CRs
that play key role in the dynamo process. For example, in the stratified box
cases with periodic boundary conditions in the horizontal directions, CRs will
find it harder to escape the disk, which may result in CR energy density
vastly exceeding that in the magnetic fields, which is contrary to
observations (e.g., Hanasz et al., 2009a). This can be avoided in global
models in which CRs can escape in all directions rather than only along the
vertical axis, and where disk spiral structure can develop enabling more
realistic CR distributions. As discussed below, general consensus that emerges
from global simulations is that the field amplification is fast, magnetic
fields reach saturation with (a fraction of) the kinetic energy of the gas,
and that the level of anisotropy of CR transport has an important positive
impact on the magnetic field amplification. While these conclusions are in
line with the findings from those based on zoom-in simulations, there are
important differences between the results based on these two approaches.
The dynamo process can generate and amplify large scale fields, but it
requires weak initial seed fields to operate. These fields can be supplied by
SNRs that can simultaneously inject CRs via diffuse shock acceleration. This
process could result in small yet astrophysically relevant fields before the
dynamo process could amplify them further. For example, explosions of $\sim
10^{6}$ pleirionic SNe could generate seed fields at the level of
$10^{-3}~{}\mu$G (Rees, 1987). Global simulations that emulate these processes
by injecting weak small-scale dipole magnetic fields associated with
individual SNe, rather than starting from a predefined initial field,
demonstrate that the volume-filling magnetic field can indeed be amplified to
equipartition with the gas energy and that small scale fields can be
efficiently reordered on timescales comparable to the galactic orbital
timescale forming large scale fields in the disk and galactic halo (Hanasz et
al., 2009b). In these simulations, the field reconfiguration is aided by
explicit Ohmic resistivity.
The tendency for the field amplification to reach saturation values in global
simulations is a general one and is robust with respect to changes in the
model initialization and assumptions. A number of global simulations
demonstrated that a fluctuating small-scale dynamo (Rieder & Teyssier, 2016)
and differential rotation (Dubois & Teyssier, 2010), can amplify the magnetic
fields up to typical values of several to several tens of $\mu$G in dense ISM
regions, a process that can be aided by CRs (e.g., Pakmor et al. 2016c;
Pfrommer et al. 2017a, 2022; Dashyan & Dubois 2020; Nuñez-Castiñeyra et al.
2022). This is in agreement with the observed values (Beck, 2000; Crutcher,
2012).
Moreover, the conclusions regarding the formation of large scale fields
extending to the halo also appear to be robust to changes in the magnitude and
topology of the initial field. In simulations starting with large scale
toroidal fields of $\mu$G strength, such fields develop despite continuous
small scale tangling of the field by exploding SNe, and the winds expel
magnetic energy into the galactic halo forming large-scale open magnetic field
lines (Hanasz et al., 2013). Compared to the CR advection-only case, the
formation of the vertical component of the magnetic field can be enhanced as
much as two orders of magnitude when active CR transport is also taken into
consideration (Dashyan & Dubois, 2020).
Important quantitative differences in the field amplification emerge when CR
transport anisotropy is taken into account. When diffusion is anisotropic (or
suppressed), CRs are initially confined to the disk because the dominant
component of the magnetic field is aligned with the plane of the disk due to
galactic rotation. This confinement drives strong gradients in the vertical
and radial magnetic fields close to the disk. These strong gradients are
responsible for the rapid rate of magnetic field amplification in the
anisotropic diffusion and purely advective cases. In the case of isotropic CR
diffusion, the outflow rate is larger and hence it also carries more magnetic
flux to the halo that has to be replenished by the dynamo in the disk, which
results in a effectively suppressed dynamo (Pakmor et al., 2016c; Nuñez-
Castiñeyra et al., 2022). Depending on the outflow rate in the isotropic
diffusion case, field amplification may not reach equipartition with
turbulence or may only saturate at equipartition after a Hubble time, while
the fields can be very quickly amplified (with a rate that depends on the
magnetic Reynolds number) when diffusion is anisotropic or absent. This
picture is consistent with the analytic galactic dynamo theory (Shukurov et
al., 2006) and the results from zoom-in simulations that suggest the existence
of a critical level of anisotropy necessary for the CR dynamo to efficiently
amplify magnetic fields (Hanasz et al., 2009a). However, contrary to
observations, for Milky Way conditions, the magnetic energy density in these
simulations is always significantly subdominant compared to the CR energy
density which may be either due to the omission of CR losses or the
geometrical effects mentioned above. An interesting by-product of the
suppression of CR diffusion that leads to the increase in both CR and magnetic
pressures is that it slows down star formation, and thus has consequences for
the net energy losses of CR due to hadronic and Coulomb effects as described
below.
Further complexity to the problem of the field amplification can be introduced
by considering cooling and collapsing galactic halos (Pfrommer et al., 2022).
In this scenario, the field amplification occurs as a result of a fluctuating
small-scale dynamo driven by (i) turbulence injected at the outward-
propagating corrugated accretion shock front, and (ii) Kelvin-Helmholtz body
modes that are excited at the interface of the supersonically (with respect to
the warm ISM phase) rotating galactic disk and the hot CGM. In tandem with the
CR energy injection at SNe, these modes source ISM turbulence throughout the
entire star-forming disk. While this mechanism exponentially amplifies the
magnetic field to reach saturation with the thermal and CR energy densities in
Milky Way galaxies, in smaller galaxies the magnetic field energy density
reaches sub-equipartition values with respect to these energy densities
(though it saturates at the kinetic energy density of the gas). Interestingly,
this model matches the global FIR-radio correlation across four orders of
magnitude in luminosity (see Section 4.3.1 for more discussion).

Figure 26: Face-on distribution of the hydrogen column density in a Milky Way-
like galaxy. Inclusion of CRs significantly smoothens the gas distribution
(cf. top/bottom panel for the case without/with CRs; Farcy et al. 2022);
reproduced with permission from MNRAS.

##### Disk structure and CR transport.

CR feedback can shape both the local and global structure of galactic disks
through driving of galactic winds and suppressing star formation. As pointed
out in several studies (e.g., Salem & Bryan 2014; Ruszkowski et al. 2017a;
Pfrommer et al. 2017a), the latter process occurs due to significant increase
in the midplane pressure provided by CRs. This additional pressure stabilizes
the disk, prevents the gas from collapsing to form very dense clouds, and thus
reduces the star formation rate. Figure 26 illustrates this effect and shows
face-on distributions of the hydrogen column density in a Milky Way-like
galaxy. Inclusion of CRs significantly smoothens the gas distribution (see
bottom panel) compared to the case without CRs (top panel, Farcy et al.,
2022).
As discussed in Section 3.2.2, observations of the enhanced $\gamma$-ray
emission in the vicinity of SNRs (Casanova et al., 2010; Hanabata et al.,
2014) suggests that CR transport speeds can be significantly suppressed near
these sources compared to the Galactic average. From the theoretical
standpoint, this suppression could be due to the excitation of non-resonant
modes that can locally amplify the magnetic field and drastically reduce the
diffusion coefficient by orders of magnitude (Bell, 2004; Reville & Bell,
2013; Caprioli & Spitkovsky, 2014b, c; Schroer et al., 2021, 2022). One
consequence of this suppression of transport is that the accumulation of CRs
near SNe and superbubbles further increases local CR pressure gradients that
prevent the formation of massive gas clumps and significantly suppress star
formation (Semenov et al., 2021). This is illustrated in Fig. 27, which shows
the impact of CRs on global galactic structure. Most notably, this figure also
demonstrates that when CR diffusion is locally suppressed in the vicinity of
SNe, the disk can maintain a “grand design” spiral structure. This is
consistent with recent discoveries of high-redshift galaxies with massive and
dynamically cold disks (Neeleman et al., 2020; Rizzo et al., 2020).

Figure 27: Impact of CRs on global galactic structure. From left to right,
columns show surface density of the gas, star formation rate, and CR pressure.
From top to bottom, different rows correspond to cases without CRs, isotropic
CR diffusion, and isotropic diffusion suppressed near the CR injection sites.
Inclusion of CRs prevents dense clump formation and allows the disk to
maintain a “grand design” spiral structure and the effect is particularly
strong when diffusion is locally suppressed. Image based on Semenov et al.
(2021); reproduced with permission from ApJ.

From the point of view of feedback modeling, some aspects of the impact of CR
diffusion suppression are qualitatively similar to those seen in the “delayed
cooling” approach that is sometimes employed in galaxy formation simulations
(e.g., Thacker & Couchman 2000; Agertz et al. 2011; see also Section 3.2.4).
In the latter approach, radiative cooling is switched off for some time
following SN explosions. Rather than getting quickly radiated away, the
feedback energy can locally accumulate and strong pressure gradients can
develop that are capable of dispersing dense gas clouds, which leads to the
suppression of star formation. This effect is also similar to that due to
increased star formation efficiency or including an ad hoc SN energy and
momentum boost. However, in the case of suppressed CR diffusion, this effect
is modeled more from first-principles, is better physically motivated and, in
contrast to the effective feedback models, CR feedback also affects the state
of the ISM and wind launching in other fundamental ways.

##### CR energy losses and $\gamma$-ray emission.

The issues of magnetic field amplification, multiphase structure of the ISM,
and CR transport all play important and interrelated roles in determining the
net CR energy losses and, thus, the $\gamma$-ray emission due to pion
production. Below we discuss the latest results from global simulations that
make predictions for the $\gamma$-ray luminosities of galaxies. Since there is
currently no consensus on the single most important factor responsible for the
significant differences between these predictions, we follow this discussion
by examining the various factors that may play important roles in determining
the level of $\gamma$-ray emission.
Global non-cosmological simulations of dwarfs and $L_{*}$ galaxies employing
the FIRE-2 model demonstrated that $\gamma$-ray luminosities are significantly
overproduced compared to observations when CR transport occurs only via
advection and streaming (Chan et al., 2019). However, in these simulations,
agreement with the observations can be achieved when very fast CR diffusion is
included with diffusion coefficients in the range $10^{29-31}$ cm2 s-1. In
these simulations the streaming transport is subdominant compared to
diffusion. One fundamental reason for the CR energy losses to increase for
slower CR diffusion, or when only streaming rather than diffusion is included,
is that in this case CRs spend more time in the densest regions, where the
hadronic and Coulomb CR cooling rates are the highest (e.g., Dashyan & Dubois
2020, who confirm that streaming losses are generally subdominant in dwarfs).
Even if CRs do not spend a significant amount of time and lose a lot of energy
in the densest gas due to efficient wave damping (Farber et al., 2018), large
values of the diffusion coefficient may be needed in the tenuous warm ionized
medium and CGM surrounding these clouds in order to reduce the losses and
excessive $\gamma$-ray emission as demonstrated in cosmological FIRE-2
simulations (Hopkins et al. 2021c, see also Section 3.3).
However, these results are in tension with other reports suggesting that
$\gamma$-ray luminosities are fully consistent with the observational Fermi
constraints even for typical diffusion coefficient values of $3\times 10^{28}$
to $10^{29}$ cm2s-1, equivalent to isotropic values comparable to or smaller
than $3\times 10^{28}$ cm2s-1 (Pfrommer et al., 2017b; Buck et al., 2020;
Nuñez-Castiñeyra et al., 2022). Specifically, these simulations suggest that
typical luminosities are an order of magnitude smaller than those based on the
above FIRE-2 simulations for the same level of diffusion, and about five times
smaller than those based on AREPO simulations when measured per star formation
rate (Werhahn et al., 2021b). Moreover, global simulations that include
suppression of CR diffusion near the sites of SN explosions also find
agreement with the observed star formation vs. $\gamma$-ray luminosity
relation (Semenov et al., 2021).
There are several factors that may be responsible for the discrepancies
between the above results. The key factors that may have the most important
impact on $\gamma$-ray emission are: (i) the level of magnetic field
amplification that affects the speed of CR transport and energy losses, (ii)
ISM structure and its relationship to the spatial distribution of CR feedback
injection sites, and (iii) modeling of CR spectra. The speed of CR streaming
and the associated streaming losses increase with the strength of the magnetic
field. As mentioned above, it has been suggested that streaming may be
subdominant compared to diffusion in terms of CR transport speed and that CR
streaming losses may be much smaller than Coulomb and hadronic CR cooling
rates. While the exact comparisons of the magnetic field strengths in the
simulations considering CR losses and the associated $\gamma$-ray emission are
difficult due to the non-uniformity of the data sets (i.e., due to the
different ways in which the fields are measured with groups reporting averages
vs. profiles, volume- vs. mass-weighted values, averages in volumes of
different sizes and shapes, or because the values are reported for halos of
different mass), the strengths of magnetic fields present in the simulations
likely vary substantially between simulations performed by different groups.
Consequently, the streaming speed may be underestimated when the magnetic
fields are relatively weak even if the streaming boost due to the efficient
wave damping processes, or the increase in the Alfvén speed due to gas
ionization, is very large (in this case, the streaming losses would also be
underestimated). Similarly, even if the magnetic fields are close to
equipartition but the streaming boost is only moderate, the effects of
streaming could be too weak. In addition to these factors, magnetic field
amplification, in conjunction with CR pressure support, affects the level of
star formation (e.g., Nuñez-Castiñeyra et al., 2022) and thus indirectly
reduces the $\gamma$-ray emission. We return to these issues in Section 3.3
below.
Setting aside the issues of the magnetic field amplification, $\gamma$-ray
luminosities can also be suppressed if the distribution of CRs correlates
weakly with the gas density and vice versa. This can be most easily
illustrated by comparing results from simulations performed for identical
initial conditions setups, resolutions, and numerical tools, but with
different prescriptions for the physics of the ISM. In this context, a
comparison of global simulations that follow the development of thermal
instability in the cold gas (Nuñez-Castiñeyra et al., 2022) and those that do
not lead to the multiphase structure of the ISM (Dashyan & Dubois, 2020)
reveals that the latter result in moderately reduced $\gamma$-ray emission per
star formation rate. For similar reasons, in the simulations where CR
diffusion is suppressed near feedback sites, $\gamma$-ray luminosites and CR
energy losses are significantly suppressed (Semenov et al., 2021) because CRs
primarily occupy underdense superbubbles, where collisions between CRs and
thermal ISM ions are less frequent. On a qualitative level, this is consistent
with the conclusions from controlled zoom-in numerical experiments that allow
for relative offsets between the locations of gas density peaks and the SNe
injecting CRs (Simpson et al., 2023). Such offsets could be due to SN
explosions in regions that have been evacuated by prior explosions, radiation
pressure, or stellar winds. These models are broadly consistent with
$\gamma$-ray observations.
Finally, incorporation of energy-dependent CR fluid can result in the
reduction of the level of $\gamma$-ray emission if steeper spectra are
considered (Nuñez-Castiñeyra et al., 2022) or when the spectra are allowed to
age due to CR losses (Werhahn et al., 2021c; Girichidis et al., 2022).

###  3.3 Cosmological effects of cosmic ray-driven winds

Complete understanding of CR feedback processes requires insights from
cosmological simulations capable of simultaneously resolving the structure of
the CGM and ISM. Such an extreme dynamical range is necessary because certain
aspects of CR-driven outflows cannot be fully captured in idealized global
models of isolated galaxies. Rather than relying on approximate prescriptions
for the distribution of mass and angular momentum in the CGM in global models
of isolated galaxies, cosmological simulations allow for self-consistent and
most realistic modeling of the accretion of the gas onto galaxies from the CGM
and the interaction of the winds with the accreting gas. The cosmological
approach also offers the most straightforward way to model the interaction of
a CR-dominated CGM with the accretion shocks, and to incorporate mergers and
continuous accretion of gas onto galaxies, which gives rise to extended (and
perhaps even bursty) star formation histories rather than a single episode of
accretion and one burst of CR injection as expected in global models of
isolated galaxies. Such a long-term cosmological evolution of galaxies
provides more robust predictions for the role of CR feedback in shaping
structural properties of galaxies and the surrounding CGM, which is discussed
next.
We caution that while the cosmological approach comes with clear advantages,
such as the possibility to realistically model boundary conditions (that is
not possible with zoom-in simulations of galactic patches or isolated global
models) to simulate the interaction of outflows with the CGM or the formation
of “galactic fountains,” it also comes with certain trade offs. These trade
offs are related to limited dynamical range of the simulations and involve the
need to resort to the effective modeling of the ISM physics as well as small
scale physics of mixing using sub-grid models (Fielding et al., 2020;
Weinberger & Hernquist, 2023).

####  3.3.1 Analytic considerations

Considerable insight into the role of CRs in shaping the properties of
galaxies, the CGM, and the dependence of the solutions on key model parameters
can be gleaned from a simple analytic model that incorporates CR injection and
diffusion out of gaseous halos in dark matter potentials (Hopkins et al.,
2020; Ji et al., 2020). In this approximate model, it is assumed that CRs are
generated by SNe at the energy injection rate of $\dot{E}_{\rm cr}$ at the
center of an isothermal sphere characterized by a circular velocity $V_{\rm
c}$, scale radius $r_{\rm s}$, and gas fraction $f_{\mathrm{gas}}$. In steady
state, and for a constant diffusion coefficient $\kappa$, this leads to a
simple distribution of CR energy density
$\varepsilon_{\mathrm{cr}}=\dot{E}_{\rm cr}/(4\pi\kappa r)$ and the
corresponding CR pressure gradient $\boldsymbol{\nabla}P_{\rm
cr}=-\varepsilon_{\mathrm{cr}}\boldsymbol{\hat{r}}/(3r)$, which can be used to
compare the magnitude of the CR pressure force opposing the gravitational
force $|\boldsymbol{\nabla}P_{\rm
cr}|/|\rho\boldsymbol{\nabla}\Phi|\sim\dot{E}_{\rm cr}Gr_{\rm s}/(3\kappa
f_{\rm gas}V_{\mathrm{c}}^{4})$. Assuming that (i) the CR energy injection
rate is proportional to the star formation rate $\dot{M}_{\star}$, i.e.,
$\dot{E}_{\rm cr}=\zeta_{\rm cr}\epsilon_{\rm sn}\dot{M}_{\star}$ (where
$\zeta_{\rm cr}\sim 0.05$ is the CR acceleration efficiency and $\epsilon_{\rm
sn}\sim 10^{51}{\rm erg}/(100~{}\mathrm{M}_{\odot})$ is the SN yield per solar
mass of star formation), and (ii) for sub-$L_{*}$ galaxies
$\dot{M}_{\star}\sim\alpha M_{\star}/\tau_{\mathrm{H}}$ (where $M_{\star}$ is
the mass in stars (Mitra et al., 2017), $\tau_{\mathrm{H}}\approx
14\,(1+z)^{-3/2}$Gyr is the Hubble time, $\alpha\sim 1$ which depends very
weakly on halo mass), the ratio of the CR to gravitational force becomes

$\frac{|\boldsymbol{\nabla}P_{\rm
cr}|}{|\rho\boldsymbol{\nabla}\Phi|}\sim\frac{\alpha\zeta_{\rm
cr,0.05}}{f_{\rm
gas,0.1}\kappa_{28}(1+z)^{3/2}}\left(\frac{M_{\star}}{M(<r_{\mathrm{s}})}\right),$
(66)

where $\kappa_{28}$ is the diffusion coefficient in units of $10^{28}$ cm2
s-1, $\zeta_{\rm cr,\,0.05}$ is the CR acceleration efficiency in units of
0.05, and $f_{\rm gas,\,0.1}$ is the gas mass fraction in units of 0.1. Given
that $M_{\star}\propto M(<r_{\mathrm{s}})^{2}$ for sub-$L_{*}$ galaxies in the
mass range $\sim 10^{11-12}~{}\mathrm{M}_{\odot}$, equation (66) reveals the
following interesting trends: (i) the relative importance of CR pressure
forces can be high in Milky Way-mass galaxies. For a constant $\kappa$, the
force ratio declines for dwarfs where $M_{\star}/M(<r_{\mathrm{s}})$ is
smaller. (However, CR transport could self-regulate in dwarfs so that the
effective $\kappa$ decreases and compensates for the impact of smaller
$M_{\star}/M(<r_{\mathrm{s}})$ thus making CR feedback also important in
smaller systems.) (ii) Slower CR transport speed allows for more CR pressure
buildup (see equation 66), which helps to drive the wind. (However, slower
transport also increases the average residence time of CRs in the ISM and CGM.
This could either increase CR energy losses compared to the CR energy
injection rate, which would in turn weaken the impact of CRs, or cause the
region to adiabatically expand due to strong CR pressure buildup, which would
then lower CR collisional losses.) (iii) The impact of CRs on the dynamics is
maximized at lower redshifts provided the above simple relationship between
$\dot{M}_{\star}$ and $\dot{E}_{\rm cr}$ applies. If instead the star
formation history is burstier in dwarfs (for which there is observational
evidence; Weisz et al., 2012; Kauffmann, 2014), this would decrease the
timescale of a star formation event $\tau_{\mathrm{H}}$ and thus temporarily
increase $\varepsilon_{\mathrm{cr}}$.
We stress that while this simple model offers first insights into the impact
of CRs, it does have limitations that stem from making simplifying assumptions
related to the neglect of more realistic CR transport, CR energy losses due to
streaming heating, and neglect of the impact CR inertia on CR forces imparted
to the gas, all of which can affect the distribution of CRs or effective CR
force acting on the ISM and CGM and CR energy losses. This in turn can shift
the regime of parameters where the effect of CRs is maximized. Nevertheless,
the model provides a useful reference framework that reveals some of the
trends seen in the simulations. In particular, it suggests that CR feedback
should be very important for halo masses ranging from LMC to Milky Way mass
($\sim 10^{11-12}~{}\mathrm{M}_{\odot}$). For smaller halo masses, more work
is needed to study the self-regulation of CR transport and whether the small
star formation rates are sufficient to build up sufficiently large CR
pressure. Similarly, for halo masses much larger than that of the Milky Way,
star formation rates are too low for CRs from SNe to produce dynamically
important pressure gradients especially given generally larger CGM pressures.
We will come back to this halo mass regime in Section 3.5 and discuss CR
feedback associated with AGN jets.

####  3.3.2 Structural properties of galaxies

While the first cosmological simulations to include the effects of CRs focused
on CR production in galaxy clusters, their contribution to the extragalactic
isotropic $\gamma$-ray background, and the physics of radio halos and relics
(e.g., Miniati et al., 2001a, b; Miniati, 2003; Pfrommer et al., 2007;
Pfrommer et al., 2008), an important role of CR feedback in shaping the
structural properties of galaxies was recognized early on with the advent of
first cosmological galaxy formation simulations that included CR physics
(Jubelgas et al., 2008; Wadepuhl & Springel, 2010; Salem et al., 2014; Chen et
al., 2016). These early hydrodynamic simulations demonstrated that CR feedback
can have a dramatic impact on the baryon distribution and kinematics in
galaxies. Specifically, they showed that CRs may significantly reduce star
formation rates (Jubelgas et al., 2008) and thus help to address the “missing
satellites” problem (Klypin et al., 1999; Moore et al., 1999) – by
substantially suppressing the star formation efficiency in the smallest
galaxies, the relatively small number of observed satellites could be
reconciled with the theoretically predicted numbers of dark matter satellites
(Wadepuhl & Springel, 2010). Furthermore, CRs could help to solve the “angular
momentum problem,” in which the efficient loss of the baryonic angular
momentum leads to galactic disk sizes that are too small compared to
observations (Jubelgas et al., 2008). This positive impact of CR feedback on
producing thin and extended galactic disks persisting to low redshifts was
seen both in dwarf systems (Chen et al., 2016) and Milky Way analogs (Salem et
al., 2014). CR feedback effects also led to flatter and more realistic
rotation curves when CR transport was included in the simulations.
Interestingly, in dwarf galaxies, the inclusion of CR physics, and in
particular CR diffusion, resulted not only in flatter and more extended disks
but also in a good match of the models to the observed baryonic Tully-Fisher
relation (Chen et al., 2016).
Adding magnetic fields could in principle already modify the structural
properties of galaxies without CRs through angular momentum transport by
magnetic tension forces. While this scenario does not seem to be realized in
practice in cosmological simulations of the growth of isolated galaxies
(Pakmor et al., 2017; Hopkins et al., 2020, using the MHD approximation),
spiral galaxies reforming after a major merger event show more extended disks
and better developed spiral structure in comparison to their pure
hydrodynamical analogues (Whittingham et al., 2021). This different galactic
morphology is caused by efficient angular momentum redistribution by magnetic
fields, which increases the central baryonic concentration after the galaxy
merger. This prevents the formation of a bar, which forms in the simulation
without magnetic fields, and thus implies a different star formation surface
density and feedback strength when taking into account magnetic fields
(Whittingham et al., 2023).
Early hydrodynamic models with CR feedback (Chen et al., 2016) suffered from
too low gas fractions compared to the observed systems and from cuspy dark
matter profiles that are also in tension with observations (e.g., Moore, 1994;
Walter et al., 2008). Moreover, more recent CR MHD models with subgrid-scale
feedback tailored to solve the above mentioned problems found either a weaker
or even opposite trend regarding disk sizes (Buck et al., 2020; Hopkins et
al., 2020). In agreement with analytic trends presented above, these models
adopt an (effectively) constant CR diffusion coefficient and suggest that the
effect of CRs in dwarf galaxies should be minor while the impact of SN CR
feedback could be very important in the halo mass range $\sim 10^{11}$ to
$\sim 10^{12}~{}M_{\odot}$. Depending on the details of the implementation of
stellar feedback and CR physics, these models predict that global or
integrated properties of galaxies either vary little when CRs are included
(Buck et al., 2020) or have a significant effect on global properties, such as
star formation rates that can be reduced by as much as a factor of $\sim 5$,
or total stellar mass that can be suppressed by a factor of a few due to CR
pressure gradients that support the CGM and slow down gas accretion (Hopkins
et al., 2020). These effects can consequently substantially reduce central
peaks in the galactic rotation curves.

Figure 28: Face-on and edge-on projections of stellar density (top two rows;
K-, B- and U-band stellar luminosities, are shown in red, green, and blue,
respectively) and the corresponding projections of the ratio of CR-to-gas
pressure (bottom). From left to right, columns correspond to the following
cases: (i) no CRs, (ii) anisotropic CR diffusion including Alfvén wave heating
of the gas, (iii) anisotropic CR diffusion, and (iv) CR advection only. Notice
the dramatic impact of the CR physics on the properties of the galaxy. Image
from Buck et al. (2020); reproduced with permission from MNRAS.

Irrespective of the magnitude of the impact of CRs on the global properties of
galaxies, CRs substantially modify the structural properties of the CGM such
as its density, temperature, and kinematics and, consequently, alter the
angular momentum acquisition by the galactic disks during their formation in
cosmological simulations. This can influence galactic properties such as disk
size or the ratio of the disk mass to total stellar mass with outcomes
strongly dependent on the physics of CR transport. For example, in simulations
of Buck et al. (2020), when CR advection or anisotropic diffusion are
included, the resulting CR-driven outflows tend to be more spherically
symmetric compared to those that either include additional CR streaming losses
or neglect CR physics altogether in their baseline Auriga model (Grand et al.,
2017, 2019). This outflow geometry interferes with coherent accretion near the
disk midplane and thus suppresses accretion of high angular momentum gas
leading to more compact disks. This is evident from the two top rows in Fig.
28 that show the projections of stellar density for the Auriga model: (i)
baseline model without CRs; (ii) model that includes CR advection, anisotropic
CR diffusion, and CR streaming heating; (iii) model with CR advection and
anisotropic CR diffusion; (iv) model with CR advection only. The trend for the
disks to decrease in size when CR physics is included is also present in
FIRE-2 simulations (Hopkins et al., 2020) though the outflow geometry is
different. This illustrates that exact conclusions are sensitive to the
physics of transport and how it shapes the properties of the CGM, which we
discuss next.

####  3.3.3 CGM thermodynamics, kinematics, and cosmic ray transport

The CGM lies at the intersection of the inflowing pristine intergalactic
medium (e.g., Kereš et al., 2009; Ocvirk et al., 2008; van de Voort et al.,
2011) and the outflowing metal-enriched gas from galactic disks (e.g.,
Oppenheimer & Davé, 2008; Behroozi et al., 2010; Davé et al., 2012) thus
making it a useful probe of the processes governing galaxy evolution. First
cosmological simulations of galaxy formation focusing on the impact of CRs
from SNe demonstrated that the CGM pressure support can be dominated by CRs
(Jubelgas et al., 2008; Salem et al., 2016; Chen et al., 2016) just as in CR
simulations of global models of isolated galaxies (e.g., Booth et al., 2013;
Salem & Bryan, 2014; Butsky & Quinn, 2018) and zoom-in simulations of
stratified boxes and disk patches (e.g., Girichidis et al., 2016b; Simpson et
al., 2016). In agreement with these conclusions, in the cosmological
simulations, the CGM was also more diffuse and characterized by substantially
lower temperatures (smaller than $\sim 10^{5}$ K) when CRs were included.
While these cosmological simulations neglected Coulomb, hadronic, and
streaming CR losses, they were nevertheless in broad agreement with the
observations suggesting that a significant fraction of the CGM may be
relatively cool (Tumlinson et al., 2011; Werk et al., 2013; Werk et al.,
2014).

Figure 29:  Top row: temperature distribution in the CGM of a simulated Milky
Way-size galaxy without (left panel) and with CR physics (right panel).
Warm/hot ($T\gg 10^{5}$K; red), cool ($T\sim 10^{4}-10^{5}$K; green), and cold
neutral gas ($T\ll 10^{4}$K; magenta/white). Bottom row: phase-space diagrams