∎ #1 #1

11institutetext: ∗shared first authorship with the authors contributing
equally to this work
1Department of Astronomy, University of Michigan, 1085 S. University Ave., 323
West Hall, Ann Arbor, MI 48109-1107, USA
2Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482
Potsdam, Germany
3Max-Planck-Institut für Astrophysik Karl-Schwarzschild-Str. 1 85748 Garching,
Germany
11email: Mateusz Ruszkowski<EMAIL_ADDRESS>Christoph Pfrommer:
<EMAIL_ADDRESS>
# Cosmic ray feedback in galaxies and galaxy clusters

A pedagogical introduction and a topical review of the acceleration,
transport, observables, and dynamical impact of cosmic rays

Mateusz Ruszkowski1,3,∗ Christoph Pfrommer2,∗

(Received: date / Accepted: date)

###### Abstract

Understanding the physical mechanisms that control galaxy formation is a
fundamental challenge in contemporary astrophysics. Recent advances in the
field of astrophysical feedback strongly suggest that cosmic rays (CRs) may be
crucially important for our understanding of cosmological galaxy formation and
evolution. The appealing features of CRs are their relatively long cooling
times and relatively strong dynamical coupling to the gas. In galaxies, CRs
can be close to equipartition with the thermal, magnetic, and turbulent energy
density in the interstellar medium, and can be dynamically very important in
driving large-scale galactic winds. Similarly, CRs may provide a significant
contribution to the pressure in the circumgalactic medium. In galaxy clusters,
CRs may play a key role in addressing the classic cooling flow problem by
facilitating efficient heating of the intracluster medium and preventing
excessive star formation. Overall, the underlying physics of CR interactions
with plasmas exhibit broad parallels across the entire range of scales
characteristic of the interstellar, circumgalactic, and intracluster media.
Here we present a review of the state-of-the-art of this field and provide a
pedagogical introduction to cosmic ray plasma physics, including the physics
of wave-particle interactions, acceleration processes, CR spatial and spectral
transport, and important cooling processes. The field is ripe for discovery
and will remain the subject of intense theoretical, computational, and
observational research over the next decade with profound implications for the
interpretation of the observations of stellar and supermassive black hole
feedback spanning the entire width of the electromagnetic spectrum and multi-
messenger data.

###### Keywords:

Cosmic rays plasmas galactic winds interstellar medium circumgalactic medium
AGN feedback

††journal: The Astronomy and Astrophysics Review

###### Contents

  1. 1 Introduction
    1. 1.1 Discovery of cosmic rays
    2. 1.2 Basic properties of cosmic rays
    3. 1.3 Role of cosmic rays in galaxy evolution
    4. 1.4 Scope and context of the review
    5. 1.5 Outline of the review
  2. 2 Physics
    1. 2.1 Cosmic ray interactions with electromagnetic waves
      1. 2.1.1 Estimates of cosmic ray number densities
      2. 2.1.2 Cosmic ray-wave scattering and diffusion
      3. 2.1.3 Cosmic ray driven plasma instabilities
    2. 2.2 Cosmic ray acceleration and escape from sources
      1. 2.2.1 Cosmic ray sources
      2. 2.2.2 Particle acceleration and magnetic amplification
      3. 2.2.3 Particle escape from supernova remnants to the interstellar medium
    3. 2.3 Cosmic ray spatial transport
      1. 2.3.1 Theoretical background
      2. 2.3.2 One-moment cosmic-ray hydrodynamics
      3. 2.3.3 Two-moment cosmic-ray hydrodynamics
      4. 2.3.4 Cosmic ray energy vs. entropy methods
      5. 2.3.5 Streaming instability, wave damping mechanisms and cosmic ray self-confinement
      6. 2.3.6 Cosmic ray scattering on MHD turbulence – external confinement by turbulence
    4. 2.4 Radiative and non-radiative cosmic ray processes and their cooling times
      1. 2.4.1 Overview
      2. 2.4.2 Cosmic ray ion interactions
      3. 2.4.3 Cosmic ray lepton interactions
      4. 2.4.4 Equilibrium electron distribution
    5. 2.5 Cosmic ray spectral transport
      1. 2.5.1 Momentum-dependence of spatial cosmic ray transport
      2. 2.5.2 Numerical methods for evolving the cosmic ray momentum spectrum in space
  3. 3 Astrophysical Systems
    1. 3.1 Cosmic ray ionization
      1. 3.1.1 Direct and indirect constraints on low-energy cosmic rays
      2. 3.1.2 Low-energy cosmic ray ionization rates
      3. 3.1.3 Penetration of molecular clouds by low-energy cosmic rays
    2. 3.2 Cosmic ray-driven galactic winds
      1. 3.2.1 One-dimensional models
      2. 3.2.2 Dynamical impact of cosmic rays on multiphase medium
      3. 3.2.3 Parker instability and cosmic rays
      4. 3.2.4 Wind launching and simulations of disk patches and stratified boxes
      5. 3.2.5 Role of cosmic rays in global models of isolated galaxies
    3. 3.3 Cosmological effects of cosmic ray-driven winds
      1. 3.3.1 Analytic considerations
      2. 3.3.2 Structural properties of galaxies
      3. 3.3.3 CGM thermodynamics, kinematics, and cosmic ray transport
    4. 3.4 Thermal instability and cosmic rays in the CGM and ICM
      1. 3.4.1 Onset of thermal instability in the presence of cosmic rays
      2. 3.4.2 Non-linear outcome of thermal instability with cosmic rays
      3. 3.4.3 Impact of magnetic fields on thermal instability in hot halos
      4. 3.4.4 Cosmic ray heating of the cold thermally unstable gas
    5. 3.5 Impact of cosmic rays from AGN in massive hot halos
      1. 3.5.1 Thermalization of the AGN energy and heating of the CGM and ICM
      2. 3.5.2 Limitations of the heating mechanisms
      3. 3.5.3 The role of AGN-injected cosmic rays in the CGM and ICM in massive halos
      4. 3.5.4 Impact of cosmic rays on CGM/ICM instabilities
  4. 4 Observational signatures
    1. 4.1 Cosmic ray propagation in the Milky Way
      1. 4.1.1 Cosmic ray residence time and the boron-to-carbon ratio
      2. 4.1.2 Leaky box model: derivation and limitations in light of new cosmic ray data
      3. 4.1.3 (Semi-)analytical cosmic ray propagation models
      4. 4.1.4 Numerical integration of cosmic ray transport in the static Galaxy
      5. 4.1.5 MHD simulations of active cosmic ray transport in evolving galaxies
    2. 4.2 Cosmic ray driven/aided outflows in the Milky Way
      1. 4.2.1 Observational evidence for Galactic center outflow
      2. 4.2.2 Hadronic star formation–driven wind model
      3. 4.2.3 Leptonic AGN jet model
    3. 4.3 Non-thermal emission from galaxies
      1. 4.3.1 Radio emission
      2. 4.3.2 Energy equipartition between magnetic fields and cosmic rays
      3. 4.3.3 Gamma-ray emission
    4. 4.4 Observational signatures of cosmic ray feedback in the CGM
      1. 4.4.1 Gamma rays and constraints on cosmic ray transport in the inner CGM
      2. 4.4.2 Interplay between cosmic ray feedback and magnetic fields in the CGM
      3. 4.4.3 Cosmic ray feedback and integrated X-ray emission from the CGM
      4. 4.4.4 CGM line diagnostics of cosmic ray feedback: thermodynamics
      5. 4.4.5 CGM line diagnostics of cosmic ray feedback: kinematics
    5. 4.5 Observational evidence for cosmic rays in galaxy groups and clusters
      1. 4.5.1 Radio relics
      2. 4.5.2 Radio halos
      3. 4.5.3 Hadronic and re-acceleration models of radio halos
      4. 4.5.4 Gamma-ray emission limits
    6. 4.6 Cosmic rays from AGN jets in galaxy groups and clusters
      1. 4.6.1 Non-thermal content of lobes associated with FRI and FRII sources
      2. 4.6.2 Gamma-ray constraints on the contents of AGN bubbles
      3. 4.6.3 Probing non-thermal AGN bubble content using the Sunyaev-Zeldovich effect
    7. 4.7 Current and future multi-messenger observatories
  5. 5 Open questions and future directions
    1. 5.1 Plasma physics and cosmic ray transport challenges
      1. 5.1.1 Plasma physics challenges
      2. 5.1.2 Building a self-consistent cosmic ray transport theory
      3. 5.1.3 Incorporating complete theory of cosmic ray transport in MHD simulations
      4. 5.1.4 Cosmic ray transport and non-thermal signatures in the Milky Way and other galaxies
    2. 5.2 Astrophysical challenges
      1. 5.2.1 Effective cosmic ray transport near sources
      2. 5.2.2 Interactions of low-energy cosmic rays with molecular clouds
      3. 5.2.3 Launching galactic winds
      4. 5.2.4 Impact of cosmic rays on the CGM
      5. 5.2.5 Cosmic rays and AGN jet feedback
      6. 5.2.6 Impact of cosmic rays on the reionization of the Universe
    3. 5.3 Concluding remarks

##  1 Introduction

###  1.1 Discovery of cosmic rays

Cosmic rays (CRs) are relativistic charged particles that originate from
outside the Solar System. The highest energy CRs are the fastest nuclei in the
Universe, moving close to the speed of light. Most CRs are protons and alpha
particles. Collisions between these primary CRs and atoms or ions in space
produce secondary CRs. In particular, collisions of primary CRs with the atoms
in Earth’s atmosphere produce CR muons that can be directly detected at sea
level (thousands of these particles will pass through your body in the few
minutes it takes to read the introduction to this review). CRs were discovered
over a century ago by an Austrian-American physicist Victor Hess in Bohemia,
which was then part of the Austro-Hungarian Empire. In a series of high-
altitude balloon flights, Victor Hess measured the rate of ionization as a
function of height above ground using an electroscope. In 1912, he determined
that the ionization rate increased with height above a certain altitude up to
four times compared to the value measured at sea level, and this conclusion
was unaffected even when the experiments were conducted during a partial solar
eclipse (Hess, 1912b).111Original publication in German: Hess (1912); for
English translation published in 2018 by Angelis and Schultz, see Hess
(1912a). These experiments firmly established that the source of ionization
was extraterrestrial and that the Sun was not the main contributor to the
ionizing flux.222It is worth noting that, in 1912 Domenico Pacini reported on
the discovery that the ionizing rate decreased by $\sim$$20\%$ when
measurements were made three meters below the surface of water, and concluded
that “a significant proportion of the pervasive radiation that is found in air
had an origin that was independent of direct action of the active substances
in the upper layers of the Earth’s surface” (Pacini, 1912, for English
translation by De Angelis from 2017, see Pacini 1912). The first evidence for
cosmic “rays” being charged particles was provided by a Dutch physicist Jacob
Clay, who observed a decrease in the ionizing flux closer to the equator
during his sea voyages between Italy and Indonesia (Clay, 1927, 1928). This
finding was later interpreted as the result of CRs deflection by the Earth’s
magnetic field, demonstrating that CRs are charged particles rather than
ionizing radiation. In 1936, Victor Hess received a Nobel Prize for the
discovery of CRs.333Victor F. Hess – Nobel lecture.

###  1.2 Basic properties of cosmic rays

A variety of modern instruments have been used to determine the spectrum of
CRs across many orders of magnitude in particle energy. Figure 1 combines the
results of various measurements including data for the extragalactic
$\gamma$-ray and neutrino backgrounds for comparison. The spectrum in this
figure is multiplied by the square of particle energy in order to provide a
proxy for the energy density per decade in particle energy.

Figure 1: Cosmic ray spectrum obtained using variety of instruments. Original
figure adopted from Lenok (2022, see references therein for data sources) and
extended to include electron and positron data detected by Voyager 2 (Stone et
al., 2019). Note that Voyager 2 data has been measured outside the
magnetopause and is thus free from the effects of solar modulation (shown with
a fading color gradient towards small energies below 10 GeV for protons and
below higher energies for heavier ions). Thin colored lines connecting the
Voyager to Earth-bound low-energy CR data indicate the shape of the
interstellar CR spectra in the Local Bubble. Notably the Voyager 2 data show
an electron flux that dominates over that of protons at the lowest energies.

As stated above, the energy density of CRs is dominated by contributions from
hydrogen and helium, followed by that of heavier nuclei, electrons, positrons,
and antiprotons. In particular, the electron-to-proton energy density ratio at
energies of 10 GeV (which are large enough so that those CR fluxes are not
affected by the magnetized solar wind) is about $K_{\mathrm{ep}}\approx 0.02$.
Elemental abundances of CR nuclei closely track those seen in the Solar
System, with the notable exception of lithium, beryllium, and boron (as well
as the elements surrounding titanium, chromium, and manganese), which are more
abundant in CRs. This deviation can be attributed to the spallation process in
which CR carbon nuclei (for the light elements) and CR iron nuclei (for the
heavy element groups) collide with interstellar medium (ISM) hydrogen atoms to
form these more abundant elements.
The spectrum of CRs exhibits a number of important features. At proton
energies well below several GeV (and for heavier nuclei below energies that
are larger by a factor equal to the ion charge), it is a power-law
characterized by a spectral index $\alpha\sim 0$
($\mathrm{d}n/\mathrm{d}E\propto E^{-\alpha}$). When observed close to Earth,
the spectrum in this energy range is severely attenuated by solar modulation –
a process in which CRs interacting with the turbulent and magnetized solar
wind are effectively scattered and partially prevented from reaching Earth.
The low-energy data shown in Fig. 1 comes from the Voyager probes after they
had passed the heliopause, which is the boundary between the solar wind and
the ISM, and is thus unaffected by the modulation. The maximum contribution to
the CR energy density is provided by particles with energies of a few GeV. At
larger particle energies, the spectrum declines as a power-law in energy,
$\mathrm{d}n/\mathrm{d}E\propto E^{-2.7}$, where $n$ is the number density of
particles of energy $E$. This power-law is commonly attributed to a
combination of CR acceleration in supernova remnant (SNR) shock waves yielding
a spectrum close to $\mathrm{d}n/\mathrm{d}E\propto E^{-2.2}$ and diffusive
escape losses that steepen the spectral index of CRs remaining within the disk
by $0.5$. Given the finite duration of the shock waves and the decreasing
shock speed in the energy conserving Sedov-Taylor phase, the maximum energy to
which CR protons can be accelerated via this process is limited to about a few
$\times~{}10^{14}$ eV possibly reaching as far as the “knee” ($\sim 3\times
10^{15}$ eV) depending on the physics of magnetic field amplification,
particle transport, and ambient ISM conditions. Past the knee, the composition
of CRs includes substantial contributions from heavier ions and the all-
particle spectrum steepens to a power law with spectral index $\alpha\approx
3.1$ that extends up to the “ankle” ($\sim 10^{18}$ eV), where it flattens
again.444There have been more spectral features detected between the knee and
ankle such as the second knee (Dembinski et al., 2017; Schröder et al., 2019),
which we do not further discuss here.
At the location of the ankle, the sources of CRs likely change from Galactic
to extragalactic accelerators because, at this energy, the CR gyroradius is
already comparable to the thickness of the Galactic disk (e.g., Kampert,
2007). Collisions of CR protons with cosmic microwave background (CMB) photons
can generate electron-positron pairs resulting in a pair production dip in the
CR proton spectrum between $1\times 10^{18}$ and $4\times 10^{19}$ eV, which
could be responsible for (some of) the observed spectral flattening
(Berezinsky et al., 2006). Following the flattening, the spectrum exhibits a
sharp cutoff at around $5\times 10^{19}$ eV, the reason for which is still the
subject of active research. Initially, this cutoff was thought to be due to
the Greisen–Zatsepin–Kuzmin (GZK, Greisen, 1966; Zatsepin & Kuz’min, 1966;
Aloisio, 2013) limit, which arises because of the collisions of CR protons
with the CMB photons that lead to photo-production of pions, a process that
drains energy from protons. Additional data from the Pierre Auger
Collaboration (Aab et al., 2017) established a change in trend from a proton-
dominated CR composition at the ankle to an increasingly heavier composition
with oxygen-group elements dominating the region at the cutoff (Aab et al.,
2020). This casts doubt on the GZK interpretation and suggests that the
maximum energies of heavier CR nuclei could be limited by photo-disintegration
due to collisions with the CMB and extragalactic background light.
Alternatively, the cutoff in the CR spectrum could signal a limiting voltage
of the acceleration process because it is very challenging for astrophysical
accelerators to reach the highest CR energies (Kotera & Olinto, 2011). The
largest detected CR energy is $\sim 3\times 10^{20}$ eV, which coincidentally
is the energy of a standard tennis ball travelling at a typical serve speed of
100 miles per hour (corresponding to $\sim 160~{}\mathrm{km}$ per hour), and
exceeds the highest particle energy achievable at the world’s largest particle
accelerator, the Large Hadron Collider at CERN, by about $3\times 10^{3}$
times.555Assuming proton-proton collisions, the center-of-mass energy
$\sqrt{s}$ and the beam energy $E_{\mathrm{lab}}$ in a fixed target experiment
(measured in the laboratory frame) are related via
$\sqrt{s}=(2E_{\mathrm{lab}}m_{\mathrm{p}}c^{2})^{1/2}$ so that the energy
reachable by the Large Hadron Collider $\sqrt{s}\approx 14$ TeV corresponds to
a CR particle energy $E_{\mathrm{lab}}\approx 10^{17}$ eV, which falls short
by a factor of $3\times 10^{3}$ in comparison to the highest CR energies ever
measured. Here, $m_{\mathrm{p}}$ and $c$ denote the proton rest mass and light
speed, respectively.

###  1.3 Role of cosmic rays in galaxy evolution

Switching the focus to a seemingly unrelated topic – galaxy formation and
evolution in a cosmologically expanding Universe – allows us to relate the
underlying physics of galaxy formation to CR transport and to demonstrate that
CRs could hold the key to solving some of the most important problems in this
field. The standard cosmological cold dark matter model of the Universe with a
cosmological constant $\Lambda$ predicts that baryons should follow dark
matter potential wells and sink to their centers as a result of very efficient
cooling. However, the baryon content of dark matter halos formed during the
process of hierarchical structure formation falls significantly below the
cosmological mean baryon density, which is known as the “missing baryon
problem” (e.g., Bregman, 2007; Tumlinson et al., 2017). In particular, the
distribution of the ratios of stellar-to-dark matter halo mass as a function
of halo mass (Moster et al., 2010, see bottom graph in Fig. 2) reveals that
even at halo masses around the Milky Way ($\sim 10^{12}$ M⊙), where star
formation is most efficient, only $\sim 20\,\%$ of the available baryons are
converted into stars, and the star formation efficiency falls off rapidly
toward both ends of the mass spectrum. These missing baryons either did not
fall into dark matter potential wells (or remain undetected) or were expelled
due to feedback processes, which can be ejective (when gas is expelled from
galaxies before it can form stars) and/or preventative (when gas is heated and
prevented from accreting onto galaxies). Star formation suppression is thought
to be due to stellar and supermassive black hole feedback below and above this
critical halo mass threshold, respectively (e.g., Somerville & Davé, 2015).

Figure 2: Top: From left to right the images show the starburst galaxy M82
(infrared in red: NASA/JPL-Caltech/Univ. of AZ/Engelbracht; optical in yellow:
NASA/ESA/STScI/AURA/The Hubble Heritage Team; X-ray in blue:
NASA/CXC/JHU/Strickland), Fermi Bubbles in the Milky Way (Platz et al., 2022,
reproduced with permission from A&A), and a LOFAR image of the giant
elliptical galaxy M87 (de Gasperin et al., 2012, reproduced with permission
from A&A). Bottom: Stellar-to-halo mass ratio as a function of dark matter
halo mass (Moster et al., 2010); reproduced with permission from ApJ.

Key challenges facing commonly employed thermal and radiation pressure–driven
stellar feedback models are (i) the overcooling problem, where the thermal
energy injected by supernovae (SNe) is quickly radiated away, thus diminishing
the impact of SN explosions on regulating star formation (or, in case of high-
resolution simulation models, causing the expelled gas to be recycled too fast
and re-accreted in fountain flows), and (ii) poor coupling of stellar
radiation to the gas and thus reduced momentum imparted to the gas due to
photoionization and/or radiation pressure that open up low-density channels in
the optically thick gas and dust enshrouding star forming environment, through
which radiation can escape following the path of least resistance (Rosdahl et
al., 2015). Both of these issues severely reduce the efficiency of stellar
feedback and galactic wind launching. However, as we will argue in this
review, CRs accelerated at SNR shocks (Blandford & Eichler, 1987) provide an
efficient feedback mechanism that addresses both challenges because CR cooling
times are typically much longer than those of the thermal gas and CRs are
generally well coupled to the plasma via particle-wave interactions
facilitated by CR driven plasma instabilities as well as direct inelastic
collisions with the gas. These processes result in CRs imparting momentum to
and heating of the plasmas, and can change in fundamental ways the properties
and evolution of astrophysical objects. In fact, the CR energy density ($\sim
1$ eV cm-3) is close to equipartition with the thermal, magnetic, and
turbulent energy densities in the ISM of our Milky Way (Boulares & Cox, 1990),
suggesting that CRs are dynamically very important in regulating star
formation.
Despite travelling at speeds close to the speed of light, $c$, arguments based
on the production of secondary CRs via spallation reactions and the observed
boron-to-carbon ratio reveal that CRs are confined to the Galactic disk over
timescales ($\tau_{\rm diff}\sim 3\times 10^{7}$ yr) much longer than the
light crossing time of the disk thickness ($\sim 3\times 10^{3}$ yr; assuming
disk thickness $h_{\rm disk}\sim 1$ kpc). When the escape of CRs from the
Galaxy over these long timescales is balanced by their acceleration in SNRs,
CR pressure can build up to the observed dynamically important levels
(Ginzburg & Syrovatskii, 1964). Since a typical mean free path of these CRs
is, on average, on the order of $\lambda_{\rm mfp,cr}\sim 3h_{\rm
disk}^{2}/(c\tau_{\rm diff})\sim$ 0.3 pc, i.e., much smaller than for photons,
CRs may be well coupled to the ISM, and their diffusion from the disk can
establish dynamically important pressure gradients that may drive strong
galactic winds in late type galaxies, as first realized by Ipavich (1975).
Transport of CRs into the circumgalactic medium (CGM) can affect its
observational properties and provide non-thermal pressure support in the CGM,
that could in turn affect thermal instability, precipitation, and thus feeding
of star formation in the disk as we will explain in this review. On even
larger scales characteristic of active galactic nuclei (AGN) jet outflows in
early type galaxies in galaxy groups and clusters, CRs accelerated in
relativistic jets may play a key role in addressing the classic cooling flow
(or cool core) problem by facilitating efficient heating of the intracluster
medium (ICM) and preventing excessive star formation (e.g., Guo & Oh, 2008).
Even when the CR pressure support is low in the CGM and ICM, CR heating of the
gas can be competitive with radiative cooling, thus helping to prevent
excessive mass accretion rates.
It is evident from the above considerations that CRs should be crucially
important for our understanding of cosmological galaxy formation and
evolution. The last decade has seen tremendous progress in our understanding
of the role of CRs in these feedback processes. On the theoretical side, our
ability to model CRs has improved dramatically thanks to a significant and
concerted effort in the community to model from first-principles the transport
of CRs in the context of galactic wind launching and AGN feedback. One of the
key challenges is to identify and fully explore the essential physics needed
to realistically model CR feedback. This will likely be accomplished via a
combination of research involving particle-in-cell (PIC) simulations of CR
transport, one-dimensional models needed to interpret results from multi-
dimensional and multi-physics simulations, and fully cosmological high-
resolution simulations incorporating the physics-driven parameterizations of
CR transport and momentum and energy deposition from stars and AGN. On the
observational side, new data from telescopes such as LOFAR, MeerKAT, Jansky
VLA, Fermi, H.E.S.S., MAGIC, VERITAS, HAWC, IceCube, and space missions such
as Voyager, AMS-02 and others have transformed our understanding of the role
of CRs in feedback processes, and new and upcoming observatories such as
Square-Kilometre Array (SKA; radio), James Webb Space Telescope (JWST; optical
to the mid-infrared), and Cherenkov Telescope Array (CTA; $\gamma$ rays) will
continue to drive progress in this very active field of research. This
conviction is shared by the Decadal Survey on Astronomy and Astrophysics 2020
panel which has voiced its appreciation of the importance of CR feedback in
galaxy formation: “The impact of CRs is one of the largest uncertainties in
understanding feedback in galaxy formation. The primary uncertainty is how CRs
are scattered by small-scale fluctuations in the magnetic field, which sets
whether CRs can escape a region or whether their pressure builds up to the
point where it can drive an outflow. […] It is remarkable that tiny solar-
system scale fluctuations in the galactic magnetic field are a key ingredient
in understanding how galaxies drive winds on scales of tens of kiloparsecs, or
that the large scale magnetic field properties or distant supernovae can
affect the formation of pre-stellar cores” (Decadal Survey on Astronomy and
Astrophysics 2020 (Astro2020) et al., 2021).
Perhaps the following closing statement from the Nobel lecture by Cecil F.
Powell in 1950 accurately captures not only the spirit of the heady years
after the original CR discovery but also the current sentiment and excitement
in the field of galaxy formation and CR feedback: “It will indeed be of great
interest if the contemporary studies of the primary radiation lead us […] to
the study of some of the most fundamental problems in the evolution of the
cosmos.”666Closing sentence of the Nobel lecture “The cosmic radiation” by the
British physicist and discoverer of the neutral pion – Cecil F. Powell,
December 11, 1950; The term “primary radiation” refers to primary CRs. Powell
shared the Nobel prize with Carl D. Anderson – an American physicist and
discoverer of the positron and muon.

###  1.4 Scope and context of the review

The overarching goal of this review is to present a comprehensive and
pedagogical review of the field of the astrophysical CR feedback. Our approach
encompasses an overview of the essential physical processes that play key
roles in CR feedback followed by a systematic review of the role of CRs over a
very wide range of physical scales emphasizing in the process the connections
between CR physics operating on the various scales, and putting CR feedback
processes in context of other stellar and supermassive black hole feedback
mechanisms.
A number of other reviews discussing CRs have appeared recently. While these
reviews mention CRs, our review occupies a distinct niche as it focuses almost
exclusively on the in-depth discussion of the physics of CR feedback across a
wide range of physical scales – from the scales of individual SNe to the
scales of galaxy clusters. Here we list some the other reviews stating also
their main focus and order them from small plasma kinetic to large
cosmological scales, starting with (i) Zweibel (2017), a concise overview of
CR physics with applications to feedback; (ii) Caprioli (2016); Pohl et al.
(2020), particle-in-cell simulations of CR acceleration; (iii) Marcowith et
al. (2020), simulations of particle acceleration in astrophysical systems;
(iv) Padovani et al. (2020); Gabici (2022), impact of low-energy CRs in the
ISM; (v) Grenier et al. (2015), CR impact on the ISM, focusing on interstellar
chemistry and CR propagation in molecular clouds; (vi) Bykov et al. (2020),
high-energy particles and radiation in star-forming regions; (vii) Veilleux et
al. (2020), in-depth review of the observational aspects of galactic outflows;
(viii) Amato & Blasi (2018), CR transport in the Milky Way; (ix) Recchia
(2020), a concise overview of CR-driven winds including discussion of one-
dimensional models; (x) Yang et al. (2018), origin of the Fermi Bubbles; (xi)
Faucher-Giguere & Oh (2023), physical processes in the CGM; (xii) Hlavacek-
Larrondo et al. (2022), AGN feedback in groups and clusters, (xiii) Yang &
Bourne (2023); macro- and micro-physics of AGN jet feedback in galaxy
clusters, (xiv) Kunz et al. (2022), plasma physics of the ICM; (xv)
Vogelsberger et al. (2020), short overview of recent advances in galaxy
formation; (xvi) Kotera & Olinto (2011), astrophysics of ultra-high energy
CRs; and (xvii) Hanasz et al. (2021), numerical methods for CR transport.

###  1.5 Outline of the review

We start the review with an in-depth discussion of the physics relevant to CR
interactions with plasmas (Section 2). In that section, we describe a number
of processes that are essential for the modelling of the various astrophysical
phenomena and the interpretation of multi-wavelength observations. These
processes include particle-wave interactions and CR driven instabilities
(Section 2.1), CR acceleration (Section 2.2), spatial and spectral CR
transport (Section 2.3 and 2.5), and CR energy loss mechanisms (Section 2.4).
In the process, we also present a brief overview of the approaches to study
these processes (kinetic, hybrid, and fluid descriptions).
Following this pedagogical physics introduction, we transition to the
discussion in Section 3 of the applications of these processes to
astrophysical situations. In that section, we organize the discussion by
relevant astrophysical scales, focusing first on CR impact on small scales
before moving onto the discussion of the phenomena relevant to larger physical
scales. We also separate the discussion of low-energy CRs (with energies
$\lesssim$ GeV) that are critically important for ISM ionization, our
understanding of non-thermal emission, and calibrating transport coefficients,
but do not contribute significantly to gas pressure, from the discussion of
CRs with energies above $\sim$ GeV that play a crucial role in dynamical
feedback. Specifically, we discuss the role of CR physics in the following
astrophysical settings: low energy CR ionization in the ISM (Section 3.1),
stellar feedback and CR-driven winds (starting with one-dimensional models and
moving on to progressively more sophisticated descriptions of the dynamical
interactions of CR with the ISM including the role of CRs in the dynamics of
SNe, CR interactions with cold clouds and multiphase ISM, impact of CRs on
star formation, and the physics of galactic wind launching; Section 3.2),
impact of CRs in cosmological galaxy formation simulations (Section 3.3), the
role of CRs in thermal instability in the CGM and ICM (Section 3.4), and the
impact of CR feedback on the CGM and ICM in massive hot halos of the largest
galaxies and galaxy clusters (Section 3.5).
As the next step, we discuss the observational signatures of CR feedback
(Section 4), covering CR propagation in the Milky Way (Section 4.1), CR aided
outflows in the Milky Way (Section 4.2), non-thermal emission from galaxies
(Section 4.3), observational signatures of CR feedback in the CGM (Section
4.4) and ICM (Section 4.5 and Section 4.6), and current and future
observatories (Section 4.7). We conclude in Section 5, where we identify open
questions and future directions. While we discuss the GeV-TeV CR connection as
a means of calibrating CR feedback, we refrain from discussing issues related
to maximum energy of Galactic accelerators (e.g., hadronic PeVatrons) and do
not cover the topic of ultra-high energy CRs.

##  2 Physics

Before we discuss the details and subtleties of CR feedback in galaxies and
clusters, we introduce the basics of CR acceleration and transport. We put
particular emphasis on the different physics models used to study CRs, from
plasma-kinetic simulations to hydrodynamical models that characterize the CR
population with a small number of thermodynamic variables. We start our
discussion with CR interactions with electromagnetic waves in Section 2.1 and
explain the concept of CR pitch-angle scattering and CR driven instabilities.
In Section 2.2, we review our knowledge on CR acceleration and escape from
SNRs to the ISM. In particular, we provide an overview of the general picture
of acceleration processes and discuss in detail diffusive shock acceleration
of ions and electrons. We then transition to larger scales and discuss CR
spatial transport in Section 2.3, including a detailed explanation of one-
moment and two-moment CR hydrodynamics. This is followed by a discussion of
new results on the CR streaming instability, wave damping mechanisms and CR
self-confinement as well as CR scattering with magneto-hydrodynamic (MHD)
turbulence, which gives rise to external CR confinement by turbulence. In
Section 2.4, we discuss radiative and non-radiative CR (ion and electron)
interactions and their cooling times. Finally, we review the CR spectral
transport in Section 2.5. Specifically, we discuss how CR transport acquires a
momentum dependence in the self-confinement and external-confinement picture
of CR transport and provide an overview of the various numerical methods
developed for evolving the CR momentum spectrum in space and time.

###  2.1 Cosmic ray interactions with electromagnetic waves

In order to motivate the importance of collective interactions of CRs with
plasma waves, we first provide some order of magnitude estimates of CR number
densities. We then explain the various ways in which CRs interact with
electromagnetic waves and end by introducing CR driven resonant plasma
instabilities.

####  2.1.1 Estimates of cosmic ray number densities

CRs are charged particles that form a collisionless species that is embedded
in a magnetized background plasma, which is composed of various phases ranging
from cold clouds (with temperatures $T\sim 10$ K) to ionized and neutral warm
($T\sim 10^{4}$ K) and hot phases ($T\sim 10^{6}$ K in galaxies and $T\sim
10^{7}$–$10^{8}$ K in galaxy clusters), which dominate by volume (Cox, 2005;
Draine, 2011). In the midplane of the Milky Way, the energy densities of CRs,
magnetic fields, turbulence, and the thermal population are in equipartition
(Boulares & Cox, 1990; Cox, 2005; Naab & Ostriker, 2017), suggesting that
these components are crucial for maintaining the energy equilibrium of the
ISM. As discussed in Fig. 1, the Galactic CR population is dominated by
particles with $\sim$ GeV energies while the warm ISM is dominated by thermal
particles with energies around $\sim$ eV, which amounts to a ratio of $10^{9}$
in particle energy. Being in pressure equilibrium, we conclude that to order
of magnitude, CRs must be extremely rare and only one in $\sim 10^{9}$ ISM
particles is a CR ion (i.e., mostly a CR proton), implying a CR number density
of about $10^{-9}~{}\mathrm{cm}^{-3}$ for a mean density of the warm ISM phase
of $1~{}\mathrm{cm}^{-3}$. The intracluster plasma is characterized by
particle energies of $k_{\mathrm{B}}T\sim 1$–$10$ keV and densities $\sim
10^{-3}~{}\mathrm{cm}^{-3}$. The non-detection of $\gamma$ rays from clusters
translates into an upper limit of the CR-to-thermal pressure ratio of
$\lesssim 10^{-2}$ in the bulk of the ICM (Ackermann et al., 2010; Aleksić et
al., 2010, 2012; Arlen et al., 2012; Ackermann et al., 2014a) so that we
obtain a CR-to-background density ratio of $\lesssim 10^{-7}$, implying CR
number densities in clusters of $\lesssim 10^{-10}~{}\mathrm{cm}^{-3}$.
Because of the very low CR number density in the ISM and ICM, CRs almost never
collide with each other. The interactions of CRs with background ions via
Coulomb or hadronic particle collisions cool the CR population on too long
time scales to establish a collisional equilibrium between CRs and thermal gas
particles (see Section 2.4.2). Instead, CRs collectively interact with the
background plasma predominantly via particle-wave scattering, which
substantially reduces the effective mean free path of CRs (Wentzel, 1974).

####  2.1.2 Cosmic ray-wave scattering and diffusion

In the following we first explore the interaction of a single CR particle with
an MHD wave and then discuss collective wave-particle interactions of CR
populations. There are two types of CR interactions with MHD waves.
(i) A gyro-resonant interaction occurs when the Doppler-shifted rotation rate
$\omega_{\mathrm{r}}$ of a circularly polarized plasma wave is a multiple of
the CR gyro frequency (Schlickeiser, 2002),

$\displaystyle k_{\parallel}\mu\varv-\omega_{\mathrm{r}}=\pm
n\Omega,\qquad\mbox{where}\qquad\Omega=\frac{qB}{\gamma mc}$ (1)

is the relativistic gyro frequency of a CR (ion or electron) population with
characteristic Lorentz factor $\gamma=[1-(\varv/c)^{2}]^{-1/2}$, velocity
$\varv$, particle mass $m$ and charge $q=Ze$ (where $Z$ is the charge number
and $e$ is the elementary charge). Here, $c$ is the light speed,
$k_{\parallel}$ is the component of the wave number parallel to the mean
magnetic field, and
$\mu=\cos\theta=\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{B}/(pB)$ is the
cosine of the pitch angle between the magnetic field and momentum vectors,
$\boldsymbol{B}$ and $\boldsymbol{p}$, while $B=|\boldsymbol{B}|$ and
$p=|\boldsymbol{p}|$ are the magnitude of the magnetic field and the
relativistic CR momentum, respectively, all measured in the frame that is
comoving with the background plasma. For a positive sign of the wave vector,
the sign $\pm$ denotes right and left circular polarization of the CR. Note
that $n$ is a natural number and $n\geq 1$ denotes the order of the resonant
interaction. The case $n=1$ denotes the interaction with plasma waves
propagating parallel to the magnetic field and $n>1$ accounts for the
interaction with obliquely propagating waves, which can be seen by exploring
the geometry of the wave-particle interaction.
To understand the physics of the resonant wave-particle interaction, we
consider the motion of a particle in a uniform and constant magnetic field
along which a small-amplitude shear Alfvén wave777A shear Alfvén wave (or
Alfvén wave for short) is a type of plasma wave that arises due to
oscillations of background ions in response to a restoring force provided by
the tension of magnetic field lines. Evaluating the response of the plasma to
such perturbations results in the dispersion relation for Alfvén waves,
$\omega_{\mathrm{r}}=k\varv_{\mathrm{a}}$, where
$\varv_{\mathrm{a}}=B/\sqrt{4\pi\rho}$ is the Alfvén velocity in the Gaussian
system of units. The ion mass density $\rho$ provides the inertia for the
resulting electromagnetic wave and causes it to oscillate much slower than a
light wave, provided the plasma is non-relativistic. As such, Alfvén waves are
low-frequency (compared to the ion gyro frequency) oscillations and the
perturbations of the magnetic field are transverse to the direction of
propagation. In astrophysical plasmas, Alfvén waves propagate primarily
alongside or opposite to the direction of the mean magnetic field with the
Alfvén velocity because obliquely propagating Alfvén waves are quickly damped
via the Landau damping process (Kulsrud, 2005). travels. Because this
represents a small perturbation to the constant magnetic field, the resulting
particle motion is approximately described by a circular gyration in the plane
perpendicular to the constant field, while the particle moves at uniform
velocity along the field. Evaluating the two solutions of the resonance
condition (1), we find that for super-Alfvénic particle motions
($\mu\varv-\varv_{\mathrm{a}}>0$), the magnetic field perturbation as seen by
the particle has the same sense of polarization as the particle itself (i.e.,
the sense of rotation direction for a fixed coordinate along the field but for
varying time). For sub-Alfvénic motions, the particle also probes the same
wave polarization for each solution, but with an opposite polarization in
comparison to the super-Alfvénic case. Hence, a resonant wave-particle
interaction requires a Doppler-shifted Alfvén wave with the same rotation
direction and frequency as that of the particle gyration frequency in its rest
frame.888The process of pitch angle scattering of an electron interacting with
a circularly polarized electromagnetic wave is beautifully visualized in
videos prepared by the NASA scientific visualization studio.
(ii) The case $n=0$ in equation (1) (also called the “Landau resonance”)
describes to a non-resonant wave-particle interaction named “transit time
damping” with
$\omega_{\mathrm{r}}=k_{\parallel}\mu\varv=k_{\parallel}\varv_{\parallel}$.
This implies that the time it takes for a particle to traverse the confining
region (i.e., the “ transit time” of the particle),
$\tau=\lambda_{\parallel}/\varv_{\parallel}=2\pi/(k_{\parallel}\varv_{\parallel})$
matches the wave period, $T=2\pi/\omega_{\mathrm{r}}$. Physically, this means
that an electron or proton is confined by a magnetic mirror force, which
requires the presence of compressible electromagnetic perturbations. In MHD,
those compressible perturbations are caused by fast and slow wave
modes,999Fast and slow magnetosonic waves are compressible, longitudinal
plasma waves for which the restoring force is provided by (thermal and
magnetic) pressures. Fast magnetosonic waves propagate parallel and
perpendicular to $\boldsymbol{B}$ and are equivalent to sound waves in
high-$\beta$ plasmas, where
$\beta=P_{\mathrm{th}}/P_{B}=2c_{\mathrm{s}}/\varv_{\mathrm{a}}$ is the
thermal-to-magnetic pressure ratio and characterizes the magnetization of a
warm plasma (and $c_{\mathrm{s}}$ is the isothermal sound speed). Slow
magnetosonic waves propagate predominantly parallel to $\boldsymbol{B}$ and
are equivalent to compressible Alfvén waves in high-$\beta$ plasmas. with slow
modes containing most of the compressible energy in subsonic turbulence. The
particle surfs the wave (i.e., the particle experiences an accelerating
electrostatic field of the wave) and gains energy because head-on interactions
between particle and wave are more frequent than head-tail interactions.

Figure 3: The schematic drawing shows a CR proton (red) orbiting around a
magnetic field line (left) and resonantly interacting with an Alfvén wave
(center), i.e., when its wavelength equals the particle’s gyroradius. In this
drawing, we choose the phase of the Alfvén wave and gyrating CR such that the
resulting Lorentz forces act opposite to the parallel velocity component of
the proton along its entire orbit and decelerate it. Note that for a CR pitch
angle $<90^{\circ}$ the magnetic vector of the circularly polarized Alfvén
wave must be tilted such that the resulting Lorentz force decelerates the
parallel velocity component. Since there are no electric fields in the
reference frame of the moving wave, the proton energy (and total velocity
$\varv^{2}=(\varv_{\parallel}-\varv_{\mathrm{a}})^{2}+\varv_{\perp}^{2}$) is
conserved and the perpendicular velocity component must increase, increasing
the pitch angle between the momentum and magnetic field vectors. Switching the
magnetic field perturbations ($\delta\boldsymbol{B}$) by $180^{\circ}$ would
result in an accelerating Lorentz force along the CR orbit and a decreasing
pitch angle. Hence, for a CR population with random phase shifts between the
Alfvén wave and CR gyro orbits, the CR protons would experience random gyro-
averaged Lorentz forces with random net pitch angle changes. In consequence,
proton scattering by an Alfvén wave packet (right) can be described by a
diffusion process leading to an isotropization of the CRs in the reference
frame of the wave for frequent CR-wave scatterings (see Fig. 4). In this case,
CRs propagate on average at the Alfvén velocity, which in the ISM is about
$10^{4}$ times smaller than the speed of light at which individual
relativistic CR particles approximately travel (Jacob & Pfrommer).  Figure 4:
A visualization of the CR scattering process with an Alfvén wave in CR
velocity space, where velocities are measured in the frame that is comoving
with the Alfvén wave. An anisotropic CR distribution moving either rightward
(red) or leftward (blue) has initial values of the cosine of the pitch angle
$|\mu|=|\varv_{\parallel}/\varv|\lesssim 1$. Scattering off of an Alfvén wave
leads to smaller (larger) values of $\varv_{\parallel}$ if initially
$\varv_{\parallel}>\varv_{\mathrm{a}}$ ($\varv_{\parallel}<\varv_{\mathrm{a}}$
) while conserving the particle energy in the rest frame of the Alfvén wave
(see also Fig. 3). This process can be described as a diffusion process in
$\mu$ along the equal-energy circle shown here so that the initially
anisotropic distribution diffuses over time to assume a homogeneous
distribution in $\mu$.

Back to the resonant interaction of a CR particle with a shear Alfvén wave. As
demonstrated in Fig. 3, upon resonantly interacting with an Alfvén wave, CRs
change their pitch angle. This can be visualized by adopting the dispersion
relation of Alfvén waves, $\omega_{\mathrm{r}}=k\varv_{\mathrm{a}}$ and
evaluating the resonance condition (1) for wave propagation parallel to the
magnetic field, yielding
$\left(\varv_{\parallel}-\varv_{\mathrm{a}}\right)k_{\parallel}=\pm\Omega$.
This is a manifestation of gyroresonance and can be visualized as a circle in
CR velocity space, centered on $\varv_{\mathrm{a}}$. Shifting our frame of
reference so that it is moving with the Alfvén speed has interesting
consequences. In this frame, the Alfvén wave only retains a magnetic field
while the electric field vanishes. As a result, CRs can neither change their
energy nor their total velocity, because this would require the presence of an
electric field. In other words,
$\varv^{2}=(\varv_{\parallel}-\varv_{\mathrm{a}})^{2}+\varv_{\perp}^{2}=\mathrm{const.}$,
which is the mathematical representation of the velocity space circle in Fig.
4. The remaining magnetic field of the Alfvén wave is able to shift the pitch-
angle of the CRs in an energy conserving manner, which leads to motion of CRs
along the velocity space circle. The interaction of an Alfvén wave with a CR
population that is randomly distributed in space leads to random phase shifts
between the individual gyro motions and the wave (see Fig. 3). This causes a
random change in pitch angle for each scattering event so that we are
witnessing a random walk in $\mu$, implying that this can be described by a
diffusion process along the equal energy surface in velocity space (see Fig.
4). CR particles with $\varv_{\parallel}>\varv_{\mathrm{a}}$ decrease their
parallel velocity component (or equivalently, decrease their pitch angle
cosine $\mu=\varv_{\parallel}/\varv$, which corresponds to an increasing angle
$\theta$) and particles with $\varv_{\parallel}<\varv_{\mathrm{a}}$ increase
their $\mu$. Formally, scattering across pitch angles of $\sim 90^{\circ}$
(corresponding to $\mu\sim 0$) would require wave modes with vanishing small
wavelengths, $\lambda=2\pi/k_{\parallel}\to 0$. However, those modes do not
exist because of fast damping of these waves. It turns out that resonance
broadening as a result of non-linear effects of the interaction enables
populating the other hemisphere of CR velocity space so that the system
eventually approaches an isotropic CR pitch angle distribution (Shalchi,
2005). Mathematically, these collective CR pitch angle scatterings can
described by a diffusion process for the gyro-averaged CR distribution
$f=f(\boldsymbol{x},p,\mu)$ in $\mu$:

$\displaystyle\left.\frac{\partial f}{\partial t}\right|_{\rm
scatt}=\frac{\partial}{\partial\mu}\left[\frac{1-\mu^{2}}{2}\,\nu(p,\mu)\,\frac{\partial}{\partial\mu}f\right],$
(2)

where the time and the pitch angle are evaluated in the Alfvén wave frame,
$\nu(p,\mu)$ denotes the scattering frequency and the factor
$(1-\mu^{2})/2=(\sin^{2}\theta)/2$ derives from gyro averaging the Lorentz
force terms. A transparent derivation of quasi-linear theory of CR transport
and pitch-angle scattering of CRs starting from elementary physics
considerations can be found in Thomas (2022).
The effective spatial scattering coefficient along the magnetic fields can be
derived using the following heuristic random walk arguments. CRs can scatter
on Alfvén waves if the resonance condition (1) is met. Adopting the dispersion
relation for Alfvén waves $\omega_{\mathrm{r}}=\varv_{\mathrm{a}}k$ and
realizing that those propagate at typical speeds of tens of km s-1, i.e.,
$\varv_{\mathrm{a}}/c\ll 1$, we drop $\omega_{\mathrm{r}}\ll
k_{\parallel}\varv$ (because CRs move relativistically, $\varv\sim c$) and
obtain the simplified resonance condition $k_{\parallel}^{-1}=\mu r_{\rm L}$,
where $r_{\mathrm{L}}=pc/(qB)$ is the CR gyroradius. This defines a minimum CR
momentum that is in resonance with the Alfvén wave of a given $k_{\parallel}$,

$\displaystyle p_{\mathrm{min}}=\frac{qB}{ck_{\parallel}}.$ (3)

This corresponds to individual scattering events each lasting $\sim r_{\rm
L}/c$, i.e., a gyration time $\sim\Omega^{-1}$. In the frame of the Alfvén
wave, Lorentz forces acting during a scattering event alter the CR pitch angle
by $|\delta\alpha|\sim\delta B/B\ll 1$, where $\delta B$ is the magnetic field
fluctuation associated with the wave and $B$ is the mean magnetic field, but
do not change CR energy. Multiple and uncorrelated wave-particle interactions
lead to the random walk of CRs in pitch angle, which approaches CR
isotropization when $N_{\rm scatt}(\delta\alpha)^{2}\sim 1$, where $N_{\rm
scatt}$ is the number of scattering events. The isotropization timescale
$t_{\rm i}\sim N_{\rm scatt}\Omega^{-1}$ of the CR is thus $t_{\rm
i}\sim\Omega^{-1}(\delta B/B)^{-2}$, and the corresponding spatial diffusion
coefficient along the mean direction of the magnetic field is $\kappa\sim
c^{2}t_{\rm i}\sim cr_{\rm L}(\delta B/B)^{-2}$ (Blandford & Eichler, 1987).
In the following sections, we will improve upon this heuristic argument and
distinguish between spatial CR transport owing to diffusion and streaming, the
latter of which is the limiting case of CR transport of an almost isotropic CR
distribution in the Alfvén frame.

####  2.1.3 Cosmic ray driven plasma instabilities

In the previous section, we discussed the effect that Alfvén wave interactions
with CRs have on their motion through phase space. We implicitly assumed that
the waves have been driven by some unspecified process and have not accounted
for a possible backreaction of the CRs on wave propagation. Here, we show that
this backreaction naturally results from various resonant CR-driven
instabilities while we defer a discussion of the non-resonant hybrid
instability (Bell, 2004) to Section 2.2.2.
The nature of CR driven plasma instabilities can be understood by studying the
response of a plasma to resonant driving provided by streaming (and gyrating)
CRs. Perturbing the Vlasov equation for the CR distribution function and
deriving a Fokker-Planck equation for CR transport leads to the picture of
quasi-linear diffusion of CRs in phase space (Schlickeiser, 2002). Assuming an
equilibrium configuration with a uniform background magnetic field and no net
electric field, we can analyze the behavior of small perturbations by
linearizing the system. The resulting coupled system of equations describes
electromagnetic modes in the plasma that are excited by anisotropically
drifting CRs and consists of the dispersion tensor applied to electric
perturbations, $T^{ij}\delta E^{j}=0$ with the boundary condition that these
electric perturbations obey the linearized set of Maxwell’s equations in
Fourier space,
$\boldsymbol{k}\boldsymbol{\times}\delta\boldsymbol{E}=\omega\delta\boldsymbol{B}$
and $\boldsymbol{k}\boldsymbol{\cdot}\delta\boldsymbol{B}=0$. The dispersion
tensor, $T^{ij}$, is complicated in general (see Chapter 8 in Schlickeiser,
2002) but simplifies for parallel propagating wave modes:101010For CR driven
instabilities, modes that propagate at an oblique angle to
$\boldsymbol{B}_{0}$ have typically a subdominant growth rate with respect to
that of parallel modes (Zweibel, 2017; Kulsrud & Pearce, 1969).

$\displaystyle T^{ij}=\left(\begin{array}[]{ccc}T^{11}&0&0\\\
0&T^{22}&T^{23}\\\ 0&-T^{23}&T^{22}\\\ \end{array}\right).$ (7)

The explicit form of the matrix elements are given in Shalaby et al. (2021,
Appendix B) and we assume that the magnetic field and wave modes are aligned
with the $x$ axis, i.e., $\boldsymbol{B}_{0}=B_{0}\boldsymbol{\hat{x}}$ and
$\boldsymbol{k}=k\,\boldsymbol{\hat{x}}$.
The stability of electromagnetic plasma modes is determined by the linear
dispersion relation in Fourier space, which assesses whether a perturbation of
a certain wave mode decays or grows over time. The dispersion relation is the
determinant of the dispersion tensor,
$|T^{ij}|=T^{11}\left[(T^{22})^{2}+(T^{23})^{2}\right]=0$. The solutions to
$T^{11}=0$ represent electrostatic wave modes for which only the parallel
electric field is finite while the solutions to $T^{22}\pm\mathrm{i}T^{23}=0$
represent electromagnetic wave modes that are characterized by transverse
electric and magnetic field components and a vanishing parallel electric
field. Physically, the linear dispersion relation represents a relation
between the complex wave frequency,
$\omega=\omega_{\mathrm{r}}+\mathrm{i}\Gamma$, and the wave number $k$ of
waves propagating along the magnetic field:

$\displaystyle D^{\pm}\equiv
T^{22}\pm\mathrm{i}T^{23}=1-\frac{k^{2}c^{2}}{\omega^{2}}+\sum_{s}\chi_{s}^{\pm}=0,$
(8)

where $\chi_{s}^{\pm}$ is the linear plasma response for species $s$, which
includes the background electrons and ions as well as the drifting CR
electrons and ions (e.g., Shalaby et al., 2021). The circularly polarized
eigenmodes of the background plasma relevant for CR scattering lie on the
electron and ion cyclotron branches: on large scales with wave numbers
$kd_{\mathrm{i}}\ll 1$, the electron cyclotron branch hosts forward-
propagating Alfvén waves, which rotate at a rate
$\omega_{\mathrm{r}}=k\varv_{\mathrm{a}}$. These scales are larger than the
ion skin depth, $d_{\mathrm{i}}=c/\omega_{\mathrm{i}}$ (where
$\omega_{\mathrm{i}}$ is the ion plasma frequency) and represent the MHD
limit. The Alfvén waves turn into whistler waves111111The term whistler
describes low-frequency plasma waves generated by lightning that emits for a
few seconds in the radio band with a descending tone at the end, reminiscent
of whistling. A whistler mode represents the continuation of the electron
cyclotron branch on scales $kd_{\mathrm{i}}>1$. Like the Alfvén waves on
larger scales, whistler waves are circularly polarized electromagnetic
perturbations with a dispersion relation for wavelengths $kd_{\mathrm{e}}\ll
1$ given by $\omega_{\mathrm{r}}=\pm
k_{\parallel}kd_{\mathrm{e}}^{2}\Omega_{\mathrm{e,0}}=\pm
k_{\parallel}kd_{\mathrm{i}}\varv_{\mathrm{a,i}}$, where
$d_{\mathrm{e,i}}=c/\omega_{\mathrm{e,i}}=\varv_{\mathrm{a,e,i}}/\Omega_{\mathrm{e,i,0}}$
are the collisionless electron and ion skin depths,
$\omega_{\mathrm{e,i}}=(4\pi e^{2}n_{\mathrm{e,i}}/m_{\mathrm{e,i}})^{1/2}$
are the electron and ion plasma frequencies, $\Omega_{\mathrm{e,i,0}}$ are the
non-relativistic electron and ion gyro frequencies, and
$\varv_{\mathrm{a,e,i}}=B_{0}/(4\pi n_{\mathrm{e,i}}m_{\mathrm{e,i}})^{1/2}$
are the electron and ion Alfvén speeds. at smaller scales for
$kd_{\mathrm{i}}>1$ (with a dispersion relation
$\omega_{\mathrm{r}}=k^{2}d_{\mathrm{i}}\varv_{\mathrm{a}}$ for parallel
whistlers) in order to approach electron cyclotron waves at scales smaller
than the electron skin depth,
$kd_{\mathrm{e}}=kd_{\mathrm{i}}\sqrt{m_{\mathrm{e}}/m_{\mathrm{i}}}\gg 1$.
Those electron cyclotron waves show a constant wave rotation rate
$\omega_{\mathrm{r}}=|\Omega_{\mathrm{e}}|$ (see the upper dashed black line
in Fig. 5). The ion cyclotron branch hosts backward-propagating Alfvén waves
on scales $kd_{\mathrm{i}}\ll 1$, which directly turn into ion cyclotron waves
on smaller scales with $\omega_{\mathrm{r}}=-\Omega_{\mathrm{i}}$ (see the
lower dashed black line in Fig. 5). By contrast, a population of gyrating and
drifting CRs along the mean magnetic field excites a CR ion-cyclotron wave,
which is a rotating and propagating electromagnetic wave with the rotation
rate $\omega_{\mathrm{r}}=-\Omega_{\mathrm{i}}+k\varv_{\parallel}$ in the
background plasma frame, i.e., its rotation rate equals the Doppler-shifted CR
ion-cyclotron frequency. Hence, it is no coincidence that these gyrating CRs
also obey the gyro-resonance condition (1) for parallel wave modes, which in
this case is interpreted such that the Doppler-shifted plasma wave frequency
resonates with the CR gyro frequency.

A resonant CR driven instability emerges if the rotating electromagnetic field
of the CR ion-cyclotron wave as seen in the background frame equals one of the
circularly polarized eigenmodes of the background plasma, namely those
provided by the ion and electron cyclotron branches. As realized by Shalaby et
al. (2023), this deep insight has a simple graphical interpretation, which is
visualized in Fig. 5: the locations at which the rotation rate of the CR ion-
cyclotron wave $\omega_{\mathrm{r}}=-\Omega_{\mathrm{i}}+k\varv_{\parallel}$
intersects the wave frequency $\omega_{\mathrm{r}}$ of the electron and ion
cyclotron branches maximize the linear growth rate of a CR driven instability.
Interestingly, CR ion-cyclotron modes are driven unstable not only at the
points of intersection with the background circularly polarized eigenmodes but
also when the rotation frequency of a given CR ion-cyclotron mode is smaller
in magnitude than that of the closest background mode (in
$\omega_{\mathrm{r}}$) of either the ion or the electron cyclotron branch (see
Fig. 5; Shalaby et al., 2023).

Starting on large scales, the first two intersection points result from a
resonance of the CR ion-cyclotron wave with the backward- and then the
forward-propagating Alfvén wave, with $\omega_{\mathrm{r}}=\mp
k\varv_{\mathrm{a}}$, respectively. In consequence, the two peaks of the gyro-
resonant instability (Kulsrud & Pearce, 1969) emerge at wave numbers
$kd_{\mathrm{i}}=(\varv_{\parallel}/\varv_{\mathrm{a}}\pm 1)^{-1}$, where
$\varv_{\parallel}=\varv\mu$ is the parallel CR velocity component. This
instability can be seen to the left in the growth rate plot in the bottom
panel of Fig. 5. On smaller scales, streaming CRs excite the intermediate-
scale instability (Shalaby et al., 2021), with the shortest unstable wave mode
at $k\varv_{\parallel}=|\Omega_{\mathrm{e,0}}|=eB/(m_{\mathrm{e}}c)$, where
$\Omega_{\mathrm{e,0}}$ denotes the non-relativistic electron gyro frequency.
Interestingly, the steepening of the wave rotation rate in the whistler regime
causes the emergence of the intermediate-scale instability that extends to the
electron cyclotron scale, where the real part of the electron cyclotron branch
flattens again and approaches $\Omega_{\mathrm{e}}$.

Figure 5: Graphical visualization of resonant CR ion driven instabilities for
reduced parameters $\varv_{\mathrm{a}}=10^{-4}c$,
$n_{\mathrm{cr}}/n_{\mathrm{i}}=10^{-6}$,
$\varv_{\parallel}/\varv_{\mathrm{a}}=0.9\sqrt{m_{\mathrm{i}}/m_{\mathrm{e}}}=2.7$,
$\varv_{\perp}=\varv_{\mathrm{a}}$, and $m_{\mathrm{i}}=36\,m_{\mathrm{e}}$
for visual clarity. Shown are wave rotation rate $\omega_{\mathrm{r}}$ (top
panel) and instability growth rate in the linear regime (bottom panel) versus
wave number times the ion skin depth $d_{\mathrm{i}}$. Resonant CR driven
instabilities are excited right at the intersection points (vertical red
dashed lines) of the CR ion-cyclotron wave,
$\omega_{\mathrm{r}}=-\Omega_{\mathrm{i}}+k\varv_{\parallel}$ (red), and the
circularly polarized plasma waves (black dashed). At the gyroscale there are
two peaks, which correspond to the crossing of the CR ion-cyclotron wave with
backward and forward Alfvén waves (left to right), $\omega_{\mathrm{r}}=\mp
k\varv_{\mathrm{a}}$, respectively (zoom-in panel). The electron cyclotron
branch steepens at smaller scales to become a (parallel) whistler wave with
$\omega_{\mathrm{r}}\propto k^{2}$ until it levels off at
$kd_{\mathrm{e}}\gtrsim 1$ to turn into an electron cyclotron wave with
$\omega_{\mathrm{r}}=|\Omega_{\mathrm{e}}|$ (top panel). The bottom panel
shows the growth rate of the gyro-resonant instability on large scales and the
intermediate-scale instability toward smaller scales (larger $k$ values).
Image from Shalaby et al. (2023); reproduced with permission from JPP.

Physically, the gyro-resonant streaming instability describes the ability of
CRs to excite resonant Alfvén waves (Kulsrud & Pearce, 1969). While Fig. 5 was
computed for gyro-tropic CRs with a specific pitch angle, in practice the CR
population shows a distribution in $\mu$. In this case, we can approximate the
linear growth rate of the gyro-resonant instability for
$n_{\mathrm{cr}}/n_{\mathrm{i}}\ll 1$ (Kulsrud, 2005; Thomas, 2022):

$\Gamma_{\mathrm{gyro}}\approx\pi\int\mathrm{d}^{3}p\,\frac{1-\mu^{2}}{2}\frac{q^{2}\varv}{p}\left[p\frac{\partial
f}{\partial p}+\left(\frac{\varv}{\varv_{\mathrm{a}}}-\mu\right)\frac{\partial
f}{\partial\mu}\right]\delta_{\mathrm{D}}\left[k(\mu\varv-\varv_{\mathrm{a}})\mp\Omega_{\mathrm{i}}\right],$
(9)

where $\delta_{\mathrm{D}}$ denotes Dirac’s $\delta$ distribution and
represents the resonance condition (1) for CR-Alfvén wave interactions.
However, this condition is now interpreted from the viewpoint of the Alfvén
waves and solely selects those CRs that meet this condition of resonantly
interacting with parallely propagating Alfvén waves of wave number $k$ to
excite or damp them, depending on the overall sign of the momentum-derivatives
of the CR distribution function in the bracket of equation (9): a positive
sign of the bracket signals instability while a negative sign indicates wave
damping. Importantly, in the gyro-averaged CR rest frame, the CR gyro-motion
needs to match the wave rotation frequency so that the CRs resonate with the
Doppler-shifted wave oscillation. For a power-law momentum distribution of
CRs, $f\propto p^{-\alpha}$, this growth rate can be simplified to yield
(Kulsrud, 2005)

$\displaystyle\Gamma_{\mathrm{gyro}}\sim\frac{\pi}{4}\,\Omega_{\mathrm{i,0}}C_{\alpha}\,\frac{n_{\mathrm{cr}}(>p_{\mathrm{min}})}{n_{\mathrm{i}}}\,\left(\frac{\varv_{\mathrm{d}}}{\varv_{\mathrm{a}}}-1\right),$
(10)

where $\Omega_{\mathrm{i,0}}=qB/(m_{\mathrm{i}}c)$ is the non-relativistic ion
gyro frequency, $C_{\alpha}=(\alpha-3)/(\alpha-2)$ is a constant of order
unity, $p_{\mathrm{min}}$ is defined in equation (3),
$n_{\mathrm{cr}}(>p_{\mathrm{min}})\propto p_{\mathrm{min}}^{3-\alpha}$ is the
number density of resonant CRs, $n_{\mathrm{i}}$ is the number density of
thermal ions, and $\varv_{\mathrm{d}}$ is the CR drift velocity that defines a
frame with respect to which the CR distribution is isotropic. Equation (10)
shows that only super-Alfvénic CRs can drive Alfvén waves unstable with a
growth rate that depends on the resonant CR-to-thermal number density ratio,
which decreases towards larger CR momenta due to the diminishing spectra at
high energies for $\alpha>4$.
On smaller scales, CRs with a finite pitch angle resonantly excite parallel
electromagnetic waves on scales between the ion and electron gyro-resonances
through the intermediate-scale instability (Shalaby et al., 2021). As we can
infer from Fig. 5, CR ions drive background ion-cyclotron modes121212An
ion–cyclotron wave is a circularly polarized oscillation of ions in a
magnetized plasma that is excited by the gyro motion of background ions in
response to the Lorentz force with a dispersion relation
$\omega_{\mathrm{r}}=\pm\Omega_{\mathrm{i}}$. It represents the small
wavelength continuation of the backward propagating shear-Alfvén wave and
propagates along the magnetic field with the plasma. unstable, which are
comoving with the CR ions. While the electron modes have the wrong sense of
rotation to interact with the CR ions in the background frame, Lorentz
transforming to a frame that is comoving with the CR ions changes the sense of
rotation of these wave modes to adopt that of CR ions and hence to enable
resonance. Provided that $n_{\mathrm{cr}}\ll n_{\mathrm{i}}$ and
$\varv_{\mathrm{a}}\ll c$, Shalaby et al. (2021) approximate the peak growth
rate by

$\displaystyle\frac{\Gamma_{\mathrm{inter}}}{\Omega_{\mathrm{i,0}}}\approx\left(\frac{n_{\mathrm{cr}}}{n_{\mathrm{i}}}\right)^{3/4}+\left(\frac{n_{\mathrm{cr}}}{3n_{\mathrm{i}}}\right)^{1/3}\left(\frac{\varv_{\parallel}\varv_{\perp}}{\varv_{\mathrm{a}}^{2}}\right)^{2/3},$
(11)

where $\varv_{\perp}$ is the perpendicular CR velocity of a gyrotropic
distribution. This intermediate-scale instability typically grows faster by
more than an order of magnitude in comparison to the previously discussed
gyro-resonant streaming instability. The two peaks of the intermediate-scale
instability are approximately given by
$kc/\omega_{\mathrm{i}}=\varv_{\parallel}/\varv_{\mathrm{a}}$ and
$kc/\omega_{\mathrm{i}}=m_{\mathrm{r}}\varv_{\mathrm{a}}/\varv_{\parallel}-\varv_{\parallel}/\varv_{\mathrm{a}}$,
where $\omega_{\mathrm{i}}=(4\pi e^{2}n_{\mathrm{i}}/m_{\mathrm{i}})^{1/2}$ is
the ion plasma frequency, $m_{\mathrm{r}}=m_{\mathrm{i}}/m_{\mathrm{e}}$ is
the ion-to-electron mass ratio, and
$kc/\omega_{\mathrm{i}}=m_{\mathrm{r}}\varv_{\mathrm{a}}/\varv_{\parallel}$ is
the electron gyroscale. This amounts to a peak separation of

$\displaystyle\frac{\Delta kc}{\omega_{\mathrm{i}}}\approx
m_{\mathrm{r}}\frac{\varv_{\mathrm{a}}}{\varv_{\parallel}}-2\frac{\varv_{\parallel}}{\varv_{\mathrm{a}}}\qquad\Rightarrow\qquad\frac{\Delta
k\varv_{\parallel}}{\Omega_{\mathrm{i,0}}}\approx
m_{\mathrm{r}}-2\left(\frac{\varv_{\parallel}}{\varv_{\mathrm{a}}}\right)^{2},$
(12)

which is an excellent fit for low values of
$\varv_{\parallel}/\varv_{\mathrm{a}}$. For increasing
$\varv_{\parallel}/\varv_{\mathrm{a}}$, the peak separation decreases to the
point were both peaks merge at
$\varv_{\parallel}/\varv_{\mathrm{a}}=\sqrt{m_{\mathrm{r}}}/2$ (i.e., the
factor 2 is replaced by 4 in equation 12). For larger CR drift speeds, the
instability ceases to exist because there is no intersection of the CR ion-
cyclotron wave with the circularly polarized electron cyclotron branch at the
whistler scales. Using the graphical representation in Fig. 5, the slope of
the CR ion-cyclotron wave is too steep to intersect and resonate with the
overturning electron cyclotron branch towards small scales. While the gyro-
resonant instability is governing CR transport in interstellar,
circumgalactic, and intracluster media, the intermediate-scale instability may
also play an important role in CR scattering and transport (Shalaby et al.,
2021) as well as in pre-accelerating electrons at collisionless shocks to
enable them to participate in diffusive shock acceleration (Shalaby et al.,
2022, see also Section 2.2).

###  2.2 Cosmic ray acceleration and escape from sources

Provided the CR distribution exhibits only a small degree of anisotropy, we
can expand the CR distribution in a suitable system of basis functions
(Legendre or Taylor polynomials) and truncate the expansion at finite order
(e.g., after the dipolar or quadrupolar anisotropy) to derive a moment-based
description of CR transport. This is governed by the wave-particle scattering
rate, which itself depends on the amplitudes of resonant waves that result
from the interplay of wave growth and damping. The resulting system of
equations can be readily applied to the problem of galactic CR propagation, as
we will detail in Section 2.3. In situations where large CR fluxes (and
associated anisotropies) are present, which is realized in the vicinity of
shock fronts, one must solve the full kinetic equations and study the
excitation of plasma instabilities in the regime of large CR fluxes, which
scatter the particles and thereby reduce the CR flux, as will be discussed in
this section.

####  2.2.1 Cosmic ray sources

The sources of Galactic CRs are generally assumed to be SNRs in the Milky Way.
While there is an ongoing debate whether this source class is solely
responsible for all CRs up to the “knee” in the CR momentum spectrum (at
$E\approx 3\times 10^{15}$ GeV for CR protons while heavier CR ions have a
break at larger energies) or whether other sources contribute some flux,
energetic considerations clearly argue for a dominant contribution of SNRs for
generating the pressure carrying (GeV–TeV) CRs. This was suggested in a
visionary paper by Baade & Zwicky (1934) and made quantitative by Ginzburg &
Syrovatskii (1964); here we will sketch the argument by estimating the
luminosity required to supply all the Galactic CRs and balance their escape
losses:

$\mathcal{L}_{\mathrm{cr}}=\frac{V\varepsilon_{\mathrm{cr}}}{\tau_{\mathrm{esc}}}\sim
3\times 10^{40}~{}\mathrm{erg~{}s}^{-1},$ (13)

where $V$ is the volume of the thick galactic disk (“the CR scattering halo”)
with half-height $H\approx 2$ kpc and radius $R\approx 10$ kpc,
$\varepsilon_{\mathrm{cr}}\approx 0.8\,\mathrm{eV~{}cm}^{-3}$ is the CR energy
density, and $\tau_{\mathrm{esc}}\approx 3\times 10^{7}$ yr is the diffusive
escape time (as inferred from the boron-to-carbon ratio of CRs and assuming a
mean hydrogen density along the CR column in the scattering halo of
$\bar{n}_{\mathrm{H}}\approx 0.1~{}\mathrm{cm}^{-3}$, see Section 4.1.1). This
power can be delivered by Galactic SNe if 10% of their kinetic energy is
transferred to CRs:

$\mathcal{L}_{\mathrm{cr}}\sim
0.1\mathcal{L}_{\mathrm{sn}}\sim\frac{0.1E_{\mathrm{sn}}}{\tau_{\mathrm{sn}}}\sim\frac{0.1\times
10^{51}\mathrm{erg}}{100\,\mathrm{yr}}\sim 3\times
10^{40}~{}\mathrm{erg~{}s}^{-1},$ (14)

where $\tau_{\mathrm{sn}}$ is the average time scale between Galactic SNe.
Within the Galaxy, various other (less prominent) sources of CRs exist. These
include termination shocks found in stellar wind bubbles and superbubbles,
shocks resulting from the collision of stellar winds in binary systems, bow-
shocks generated by massive runaway stars, and pulsar wind nebulae. Time-
averaging the mass loss rate times the square of the terminal velocity of
stellar winds throughout different stages such as the main sequence, red
supergiant, and Wolf-Rayet phases, the collective stellar wind luminosity from
all massive stars in the Galaxy amounts to approximately
$\mathcal{L}_{\mathrm{w}}\approx 1\times 10^{41}~{}\mathrm{erg~{}s}^{-1}$ (Seo
et al., 2018), which is about 1/3 of the power of SN explosions,
$\mathcal{L}_{\mathrm{sn}}\approx 3\times 10^{41}~{}\mathrm{erg~{}s}^{-1}$.
Internal recollimation shocks in relativistic AGN jets as well as the
termination or back-flow shocks are other important sites of CR acceleration,
yet primarily for injecting CRs into the intergalactic and intracluster
plasmas.

####  2.2.2 Particle acceleration and magnetic amplification

##### General picture.

We now focus on the mechanism of diffusive shock acceleration that energizes
CR electrons and protons at astrophysical non-relativistic shocks (Krymskii,
1977; Axford et al., 1978; Bell, 1978a, b; Blandford & Ostriker, 1978) and
discuss the popular case of a SNR shock that causes the surrounding ISM to
expand radially. Since the mean free path in the hot phase of the ISM is
larger than the SNR, the shock wave cannot set the ISM in motion by momentum
transfer from individual particle-particle collisions, but only by particles
scattering at electromagnetic fluctuations in the surrounding plasma that are
driven unstable by the expanding flux of plasma (e.g., for negligible
magnetization through the Weibel instability as shown by Medvedev & Loeb 1999,
while other electromagnetic plasma instabilities dominate the magnetized
case). Such a shock wave is referred to as a collisionless shock wave.
Essential to the CR acceleration process is the amplification of standing
Alfvén waves by currents driven by high-energy ions (Bell, 2004) that outrun
the shock towards the upstream region to generate a foreshock precursor (see
Fig. 6 for a schematic drawing). This sources an (electron) return current in
the thermal plasma that exponentially grows a spectrum of helical Alfvén waves
with wavelengths smaller than the gyro radii of the high-energy protons. These
waves act as a mediator that absorbs energy and can effectively scatter
resonant protons and electrons at lower energies (Fig. 3) that are diffusing
into the precursor. The scattering rate is close to, or below, the Bohm limit
for a range of energies (Stage et al., 2006), i.e. they are scattered once per
gyro orbit, implying that scattering in self-excited turbulence significantly
reduces the diffusion coefficient (Reville et al., 2008).

Figure 6: Density structure of a quasi-parallel collisionless shock that
accelerates CRs. The plasma is entering from the right and is gradually slowed
down over the “foreshock precursor” where electron and ion driven plasma
instabilities drive plasma turbulence that scatters ionized particles. As the
plasma enters the shock transition of strong electromagnetic fluctuations that
stretches across several ion gyro radii $r_{\mathrm{gi}}$, frequent particle-
wave scattering ensues. This slows down the incoming flow and steeply
increases the plasma density in form of a “ramp” and an “overshoot”, after
which the plasma settles to a downstream state (Bohdan, private comm.).

This leads to isotropization of charged particles in the respective reference
frames (before and after the shock front) so that they pick up momentum
through head-on scatterings with the magnetic field without recoil to
(re)cross the shock front.131313The recoil of the magnetic field frozen into
the surrounding plasma can be neglected because of the large inertia of the
plasma. In each of these cycles (downstream, upstream and back downstream), a
part of the accelerated particles is advected with the plasma behind the shock
front, so that they do not participate in the further acceleration process.
Many acceleration cycles, i.e., repeated reflections on both sides of the
shock wave, generate relativistic CR particles. As long as the gyroradius is
smaller than the radius of the shock front, this scale invariant acceleration
process produces a (non-thermal) power law in the CR momentum distribution
with an efficiency that scales with $\varv_{\mathrm{s}}/c$, where
$\varv_{\mathrm{s}}$ is the shock velocity in the laboratory frame (which is
sometimes referred to as first-order Fermi acceleration, even though Fermi did
not discover this process). Effectively, all particle acceleration at shocks
is routed in the electrostatic field E of the shock, that slows down the
incoming flow and causes the global velocity divergence.

As originally pointed out by Fermi (1949, 1954), particles can also be
accelerated without a shock through interacting with externally driven
turbulence. Resonant scattering off of moving magnetic irregularities, with
$\textbf{{E}}=\mathbf{0}$ in the local rest frame causes isotropic and elastic
scattering in the scattering center rest frame. As a result there is a
momentum gain for head-on collision and a momentum loss for tail-on collisions
in the laboratory frame. On average, the mean energy is conserved but because
there are statistically more head-on than tail-on collisions due to the
particle motion, this causes the width of the particle distribution to
increase so that a small fraction of particles in the high-energy tail
experiences acceleration. In the quasi-linear picture of particle transport in
a bath of linear (Alfvén, magnetosonic) waves (Kennel, 1966; Skilling, 1971,
1975; Schlickeiser, 2002; Shalchi, 2009), this can be described by a diffusion
process in momentum space, which models the energy gain through resonant
interactions. This process is comparably inefficient because it scales as
$(\varv_{\mathrm{w}}/c)^{2}$, i.e. it is only second order in the
dimensionless wave velocity of the plasma wave scattering centers
$\varv_{\mathrm{w}}$ in the laboratory frame. In the following, we will not
consider this second-order Fermi process for particle acceleration (while it
may play an important role in particle transport).
The (microscopic) processes at shocks are governed by plasma kinetics and can
be studied by means of particle-in-cell (PIC) simulations that use macro
particles to represent electrons and ions of the thermal plasma and the
energetic CRs. According to Maxwell’s equations, moving charges source
currents that perturb electromagnetic fields. These generate Lorentz forces
that accelerate charged particles and hence modify the charge distribution and
currents. The system is evolved by numerically iterating this loop on a
fraction of the electron plasma timescale for macro-particles representing the
individual elementary particles of a plasma. This methodology is ideal for
exploring kinetic instabilities in the collisionless plasma around a shock.
Prominent PIC codes include Tristan-MP (Spitkovsky, 2005), VPIC (Bowers et
al., 2008), photo-plasma (Haugbølle et al., 2013), SHARP (Shalaby et al.,
2017; Shalaby et al., 2021), and Warp-X (Vay et al., 2018); see also Pohl et
al. (2020) for a review. However, this limits the physical length scale and
total simulation time of PIC simulations to microscopic dimensions. In order
to allow for longer run times of the simulations and/or multiple spatial
dimensions, the electron timescale is integrated out in hybrid-PIC
simulations, where the electron population is represented by an adiabatic
fluid and the ions are treated as macro particles in the kinetic PIC model.
Approaching even larger time and length scales (to study certain properties of
CR transport) requires to treat the entire background either as an MHD fluid
or as separate fluids for each particle species, which are coupled to the PIC
component representing CR particles.

##### Ion acceleration at shocks.

Hybrid-PIC simulations in large computational domains and long run times
demonstrate the formation of non-thermal tails at strong, non-relativistic
shocks. While the efficiency of ion acceleration far into the relativistic
regime approaches zero for quasi-perpendicular shocks (with magnetic fields
nearly perpendicular to the shock normal, $\theta_{Bn}\approx 90^{\circ}$), it
is maximized for quasi-parallel shock geometries (Caprioli & Spitkovsky,
2014a, b, see Fig. 7, bottom panels) because of efficient magnetic field
amplification through the Bell instability (Bell, 2004).

Figure 7: Hybrid-PIC simulation of CR ion acceleration at a collisionless,
non-relativistic strong shock. The top panel shows the downstream ion energy
spectrum of a quasi-parallel shock, color coded by different times. The
thermal distribution can be accurately described by a Maxwellian distribution
with a temperature that is 80% of the expected temperature for a shock (with a
Mach number of 20) that does not accelerate particles (dashed line). The
remaining energy is shared among magnetic turbulence and CRs that follow a
power-law spectrum, which has an increasing maximum energy with time. To
demonstrate the consistency with the theoretically predicted energy scaling in
the non-relativistic regime of diffusive shock acceleration of test particles,
the spectrum is scaled by $E^{1.5}$. This corresponds to a momentum spectrum
$\propto p^{-4}$ (see inset). The bottom two panels show the magnitude of the
total magnetic field for $\mathcal{M}=50$ shocks and a magnetic obliquity of
$\theta_{Bn}=0^{\circ}$ and $80^{\circ}$, respectively (see gray arrows),
implying that magnetic field amplification and thus CR acceleration is only
efficiently at work in quasi-parallel shocks. Image from Caprioli & Spitkovsky
(2014a); reproduced with permission from ApJ.

The average distance that energetic ions travel before experiencing pitch-
angle scattering is similar to their gyro radii calculated in the turbulence
they generate. In the case of moderately strong shocks, the magnetic field
amplification can be characterized in the quasi-linear regime. This implies
that all particles diffuse with their self-generated diffusion coefficient and
undergo scattering once per gyro orbit (i.e., they diffuse in the Bohm limit).
By contrast, in the case of very strong shocks, the magnetic field undergoes
significant amplification, reaching non-linear levels. Most of the magnetic
energy is concentrated in modes with wavelengths similar to the gyroradii of
the highest-energy ions. As a result, only those particles undergo Bohm-like
diffusion, while others scatter less efficiently (Caprioli & Spitkovsky,
2014c). The evolution of Bell’s (2004) instability at high Mach number, quasi-
parallel shocks is qualitatively similar in two and three-dimensional
simulations (van Marle et al., 2019).
Interestingly, quasi-parallel shocks are not stationary but they reform quasi-
periodically on ion cyclotron timescales. This indicates that when ions
encounter the steepest gradient of the shock discontinuity, they undergo
specular reflection. This only happens for a quarter of the time. As a result,
those ions are energized via shock-drift acceleration141414Here, the term
“drift” refers to the $\boldsymbol{E\times B}$ drift of charged particles,
i.e., the constant motion of the particles’ gyro centers perpendicular to the
electric and magnetic fields. While shock-drift acceleration is most effective
at perpendicular and quasi-perpendicular shocks, it does not work at a
parallel shock (Ball & Melrose, 2001). For a quasi-perpendicular magnetic
guide field, B, the induced electric field
$\textbf{{E}}=-\boldsymbol{\varv\times B}/c$ (where $\boldsymbol{\varv}$ is
the gyro-averaged particle velocity) is preferentially oriented perpendicular
to the shock normal, so that the particles drift along the shock surface. As
the particles’ gyro orbits enter the shock transition, they experience an
accelerating force from the somewhat stronger electric field there. Depending
on the orientation of the gyro orbit, the field can either accelerate or
decelerate the particle. After passing through the shock, the stronger post-
shock magnetic field tightens the gyro orbits. If the gyro orbit is still
large enough to enter the shock transition region, the increased electron
field will continue to accelerate the particle. and only a fraction of those
gain enough energy to be injected into the process of diffusive shock
acceleration (Caprioli et al., 2015). Hence, the popular thermal leakage model
where all particles above a critical momentum in the tail of the Maxwell-
Boltzmann distribution are accelerated is likely not applicable and too
simplified to capture the correct physics. Energy transfer to a non-thermal
(relativistic) particle population causes shocks to be modified and to assume
an increased density jump at the shock transition due to an increased
compressibility brought about by the decrease in the adiabatic index in the
presence of the relativistic (CR) population. As a result, flatter CR spectra
should emerge at high energies in comparison to the test-particle prediction
$\propto p^{-4}$ (Drury & Voelk, 1981; Malkov & Drury, 2001). This is in
contradiction to the steeper slopes inferred from the non-thermal emission at
SNRs implying momentum spectra $\propto p^{-4.3}$. Possible solutions of this
spectral steepening include the insight that (i) the energy required to
turbulently amplify magnetic field during the particle acceleration at shocks
extracts energy from the CR population and steepens the CR energy spectrum
(Bell et al., 2019), (ii) that CRs potentially isotropize in the frame of
Alfvénic turbulence downstream of the shock in a region known as the
“postcursor” so that they would be advected downstream with this turbulence,
implying a steepening of the CR momentum spectrum despite the density jump
exceeding the canonical value for a non-relativistic gas of four (Caprioli et
al., 2020), and (iii) varying magnetic obliquity along the shock surface
(Hanusch et al., 2019).

##### Electron acceleration at shocks.

By contrast, the gyroradii of electrons are smaller by the mass ratio relative
to ions, $r_{\mathrm{e}}/r_{\mathrm{i}}=m_{\mathrm{e}}/m_{\mathrm{i}}$, so
that they random walk through the shock transition and do not see a coherent
electrostatic shock potential. This implies that thermal electrons cannot
directly participate in the diffusive shock acceleration process that
accelerates electrons to highly relativistic energies and requires a pre-
acceleration process, which increases their momentum by a factor
$\sim(m_{\mathrm{i}}/m_{\mathrm{e}})(c/\varv_{\mathrm{s}})\sim 3\times
10^{4}$, thereby constituting the famous “electron injection problem” at
shocks. Several processes have been suggested that differ depending on the
magnetic obliquity in the upstream. Quasi-perpendicular shocks (with an
obliquity $\theta_{Bn}\sim 90^{\circ}$) are characterized by a narrow shock
transition so that particles cannot escape the shock and are (inefficiently)
accelerated via shock-surfing (Shimada & Hoshino, 2000; Hoshino & Shimada,
2002; Bohdan et al., 2019a, b), magnetic reconnection (Matsumoto et al., 2015;
Bohdan et al., 2020), stochastic Fermi acceleration (Matsumoto et al., 2017;
Bohdan et al., 2019b), and compressional heating. In oblique shocks (with an
obliquity $\theta_{Bn}\approx 50^{\circ}-70^{\circ}$), electrons can escape
the shock transition region to form the electron foreshock where they generate
electrostatic electron-acoustic waves far upstream and whistler waves closer
to the shock that can accelerate electrons (Xu et al., 2020; Morris et al.,
2022; Bohdan et al., 2022). In addition, shock-surfing, magnetic reconnection,
and stochastic shock drift acceleration also energize electrons (Matsumoto et
al., 2017; Katou & Amano, 2019; Amano et al., 2020).

However, by far the most efficient electron accelerators are quasi-parallel
shocks ($\theta_{Bn}\lesssim 50^{\circ}$) where ions can escape the shock
transition region forming the ion foreshock (Park et al., 2015; Hanusch et
al., 2020; Arbutina & Zeković, 2021). In particular, the recently found
intermediate-scale instability (Shalaby et al., 2021) provides a natural way
to produce large-amplitude electromagnetic fluctuations in parallel shocks.
The instability drives ion–cyclotron waves unstable that are comoving with the
upstream plasma at the shock front (see Section 2.1.3). These unstable ion-
cyclotron modes scatter the electrons parallel to the magnetic field so that
some of them get accelerated to the required energies to be injected into
diffusive shock acceleration (Shalaby et al., 2022). However, these PIC
simulations are numerically very challenging, which limits the achievable
simulation time and/or dimensionality of the problem. Hence, we do not yet
have a complete picture of the kinetic plasma physics that is responsible for
particle acceleration.

##### Meso-scale models of shock acceleration.

There is an enormous range of scales between an entire SNR (with a typical
radius of $r\sim 3~{}\mathrm{pc}\approx 10^{19}$ cm) and typical scales of a
plasma. In a cold plasma, the ion and electron distributions oscillate at
their characteristic plasma frequencies, $\omega_{\mathrm{i,e}}=(4\pi
e^{2}n_{\mathrm{i,e}}/m_{\mathrm{i,e}})^{1/2}$, which implies corresponding
ion and electron skin depths of $d_{\mathrm{i}}=c/\omega_{\mathrm{i}}\sim
2\times 10^{7}~{}n_{\mathrm{i,}0}^{-1/2}\mathrm{cm}$ and
$d_{\mathrm{e}}=c/\omega_{\mathrm{e}}\sim 5\times
10^{5}~{}n_{\mathrm{e,}0}^{-1/2}\mathrm{cm}$, where
$n_{\mathrm{e,}0}=n_{\mathrm{i,}0}=1~{}\mathrm{cm}^{-3}$. Those length and
time scales need to be resolved by PIC models of non-relativistic shocks at SN
blast waves, which have typical shock transition widths of a few
$d_{\mathrm{i}}$. In consequence, these PIC simulations can only afford a
limited simulation run time typically corresponding to
$300~{}\Omega_{0}^{-1}=300\,m_{\mathrm{p}}c/(eB)\sim
9B_{\mu\mathrm{G}}^{-1}~{}\mathrm{hours}$ of physical time and rarely model
two or even three spatial dimensions, though hybrid PIC simulations are able
to significantly expand upon these constraints. To bridge this range in
scales, several different meso-scale approaches have been developed, which aim
at answering different open problems of the CR shock acceleration problem. The
answer to the question of whether SNR shocks are able to accelerate CRs to the
“knee” at $\sim 3$ PeV (for CR protons and greater energies for heavier ions)
depends on the time evolution of the CR diffusion coefficient and hence on the
magnetic amplification upstream the shock (Lagage & Cesarsky, 1983a, b; Bell,
2013). Three-dimensional MHD-kinetic simulations, which employ a spherical
harmonic expansion of the Vlasov–Fokker–Planck equation (Bell et al., 2013;
Reville & Bell, 2013) find that at the present time, historical SNRs such as
Cas A, Tycho and Kepler appear to be expanding too slowly to accelerate CRs to
the knee and argue for a larger maximum energy at an earlier stage of their
evolution. Coupling the hydrodynamics to the transport equations for CRs and
magnetic turbulence in one spherical dimension and zooming onto the shock
transition enables one to resolve the scales of the fast growing Bell
instability in the shock precursor and the associated Bohm diffusion while
following diffusive acceleration of CRs (Brose et al., 2020). This yields the
time evolution of the CR spectrum and facilitates accurate comparisons to
observational data. The simulated $\gamma$-ray emission spectra (Brose et al.,
2021) of the different stages of the remnant’s evolution such as free
expansion, energy-conserving Sedov-Taylor, and early post-adiabatic expansion
phases appear to match the observed time sequence of TeV $\gamma$-ray spectra
(Funk, 2015).

##### Three-dimensional MHD models of shock acceleration.

A complementary line of research uses three-dimensional MHD simulation models,
where a suitable algorithm finds the shock surface during the run-time of the
simulation and characterises its jump conditions, Mach number, and magnetic
obliquity (Pfrommer et al., 2017a; Pais et al., 2018; Dubois et al., 2019;
Böss et al., 2023). In the next step, CR ions are injected and transported
either in the grey approximation (Pfrommer et al., 2017a) or by accounting for
the Fokker-Planck evolution of the full CR ion spectrum (Yang & Ruszkowski,
2017; Girichidis et al., 2019, 2022), which enables us to account for the
dynamical back reaction of CRs on the MHD flow. Additionally, a full spectrum
of (primary, shock-accelerated) CR electrons can be injected according to a
“subgrid” model, which is inferred from PIC simulations or motivated by
theoretical considerations. The CR electron spectrum is transported with the
flow while accounting for all re-acceleration and cooling processes that
modify the spectrum (Winner et al., 2019). At every time step, we can compute
the hadronic proton-proton reaction that generates pions, which decay into
$\gamma$ rays and secondary electrons/positrons and neutrinos (see Section
2.4). Subsequently, the (primary and secondary) leptons radiate synchrotron
emission from the radio to the X-ray band (provided the shock generated
sufficiently energetic electrons and amplified magnetic fields). In addition,
they undergo IC interactions, which results in emission from the hard X-ray to
the $\gamma$-ray regime, which can be compared to observational data.
Comparison between observed and simulated emission maps, ranging from the
radio, to X-rays to the $\gamma$-ray regime and the detailed multi-frequency
spectrum of SN1006 argue in favor of a preferred quasi-parallel electron and
ion acceleration for TeV electrons (Winner et al., 2020, see Fig. 8). While
simulated radial profiles in the equatorial sectors match X-ray data very
well, there is some residual radio flux visible in the radio observations,
which is missed by those simulations. This may suggest that the acceleration
efficiency for GeV electrons is less strongly dependent on magnetic obliquity
while the acceleration efficiency for TeV electrons, which radiate synchrotron
X-rays, is suppressed by a factor of about 10 for quasi-perpendicular shocks
in comparison to quasi-parallel shock geometries. By incorporating the
obliquity-dependent CR ion acceleration model of Caprioli & Spitkovsky (2014a)
into these comprehensive three-dimensional MHD models, we can account for the
puzzling range of morphologies observed in TeV shell-type SNRs by adjusting
the magnetic morphology and attribute regions exhibiting high (low)
$\gamma$-ray intensity to quasi-parallel (quasi-perpendicular) shock
configurations. Moreover, this enables constraining the magnetic coherence
scale in the environment of several SNRs by relating the degree of
$\gamma$-ray patchiness of SNRs to the correlation scale in the magnetic field
(Pais et al., 2020; Pais & Pfrommer, 2020).

Figure 8: Multi-frequency maps of the SNR SN 1006: the bottom row shows
observations in the 1.4 GHz radio band (Dyer et al., 2009), the soft X-ray
band (Cassam-Chenaï et al., 2008), and at very-high energy $\gamma$ rays
(Acero et al., 2010). Those are compared to simulated maps derived from global
MHD simulations of Sedov explosions, which account for CR ion acceleration at
the shock and evolve the CR electron spectrum (Winner et al., 2019).
Initially, the magnetic field points from the bottom right to the top left.
These simulations account for an acceleration efficiency that depends on
magnetic obliquity, implying efficient acceleration at quasi-parallel shocks
(i.e., when the shock travels at a narrow angle relative to the orientation of
the magnetic field in the upstream region), and inefficient acceleration at
quasi-perpendicular shock configurations. Hence, they explain the polar cap
morphology seen in the maps. The simulated $\gamma$-ray map is convolved with
the point-spread function of the H.E.S.S. instrument and incorporates Gaussian
distributed noise, which matches the observed amplitude and correlation
structure. The maps have a side length of 21 pc or 42.5’ at a distance of 1.66
kpc. Image from Winner et al. (2020); reproduced with permission from MNRAS.

####  2.2.3 Particle escape from supernova remnants to the interstellar medium

Diffusive shock acceleration is a self-regulating process: provided the level
of upstream turbulence is low, many particles escape upstream, which sources a
large current and a very efficient amplification of resonant and non-resonant
Alfvén waves through streaming instabilities (Kulsrud & Pearce, 1969; Bell,
2004). In turn, this causes the level of Alfvénic upstream turbulence to
increase until the magnetic energy density in the amplified field is
comparable to the CR energy flux density at the shock divided by the light
speed, yielding a saturated value for the amplified magnetic energy density
(Bell, 2004; Niemiec et al., 2008; Riquelme & Spitkovsky, 2009; Ohira et al.,
2009; Gargaté et al., 2010; Blasi et al., 2015; Zacharegkas et al., 2022):

$\displaystyle\varepsilon_{B,\mathrm{sat}}\sim\frac{1}{2}\frac{\varv_{\mathrm{s}}}{c}\,\varepsilon_{\mathrm{cr}},$
(15)

where $\varv_{\mathrm{s}}$ is the shock speed, which is equal the CR drift
velocity relative to the upstream plasma. The strong Alfvénic turbulence
scatters particles more efficiently and suppresses particle escape in the
upstream (Schure & Bell, 2014; Cardillo et al., 2015). Besides regulating the
shock itself, this phenomenon also has a nonlinear impact on the escape and
confinement times of CRs in the immediate vicinity of their sources (e.g.,
Bell et al., 2013; Marcowith et al., 2021). The self-generated magnetic
turbulence increases the scattering rate of CRs, resulting in a reduction of
the diffusion coefficient. Consequently, CRs accumulate within the environment
of the sources and start to build up a significant pressure gradient. This
gradient excavates a cavity around the source and cause to the formation of a
CR-dominated bubble, so that the CRs can accumulate significant grammage
$X(t)=\rho ct$ in the bubbles before moving into the ISM (Reville et al.,
2008; Malkov et al., 2013; Schroer et al., 2021, 2022; Recchia et al., 2022);
the level of the increased CR scattering depends on the ISM properties in the
immediate environment of the sources (Reville et al., 2007, 2021; Nava et al.,
2016, 2019).
A crucial point to consider when discussing CR feedback in galaxies is the
presence of a spatially varying and reduced diffusion coefficient within the
source regions in comparison to its galactic average value. This modification
affects the phase-space structure of star-forming regions, leading to the
accumulation of CRs. As a result, even in situations where the disk is
unstable, the accumulation of CRs prevents local fragmentation and preserves a
well-defined spiral structure (Semenov et al., 2021). Provided the level of
electromagnetic fluctuations is small enough so that the quasi-linear theory
of CR transport applies, advancing to the next level calls for the two-moment
model of CR transport, which accounts for a temporarily and spatially varying
diffusion coefficient in the self-confinement picture (Jiang & Oh, 2018;
Thomas & Pfrommer, 2019, 2021; Thomas et al., 2021). These hold the promise to
elucidate the picture of how exactly CRs migrate from the sites of their
increased confinement around sources to the galactic environment where
transport is faster.

###  2.3 Cosmic ray spatial transport

####  2.3.1 Theoretical background

To model CR transport in macroscopic systems such as galaxies, galaxy
clusters, or AGN jets, and to study dynamical impact of CRs on the evolution
of these systems, we need to coarse grain the kinetic physics and develop a
fluid description for CRs that is coupled to MHD. Here, we sketch the
derivation and various approximations of the different approaches. The CR
distribution is defined in phase space, which is spanned by the momentum and
spatial coordinates $\boldsymbol{p}$ and $\boldsymbol{x}$, respectively, and
reads

$\displaystyle f\equiv
f(\boldsymbol{x},\boldsymbol{p},t)=\frac{\mathrm{d}^{6}N}{\mathrm{d}{x^{3}}\mathrm{d}{p^{3}}}.$
(16)

Its time evolution is governed by the Vlasov equation in the semi-relativistic
regime in the frame that is comoving with the fluid,

$\displaystyle\frac{\partial f}{\partial
t}+(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{cr}})\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}f+\boldsymbol{F}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}f=0,$
(17)

where the distribution function is advected by the total velocity that
consists of the mean gas velocity $\boldsymbol{\varv}$ and the CR velocity
$\boldsymbol{\varv}_{\mathrm{cr}}$. Note that $\boldsymbol{\varv}$ is measured
in the laboratory frame while CR velocity and momentum $\boldsymbol{p}$ are
measured in the comoving frame. The total force is represented by
$\boldsymbol{F}$. Introducing the comoving frame introduces pseudo forces,
which are denoted as $\boldsymbol{F}_{\mathrm{pseudo}}$ and which come about
because the momentum measured by a comoving observer changes as a result of a
change of the frame velocity $\boldsymbol{\varv}$. Consequently, the momentum
of CRs is altered by the Lorentz force, which can be divided into
contributions from large-scale and small-scale electromagnetic fields denoted
as $\boldsymbol{F}_{\rm macro}$ and $\boldsymbol{F}_{\rm micro}$, respectively
(as derived with a covariant formalism by Thomas & Pfrommer 2019, Appendix C):

$\displaystyle\boldsymbol{F}$ $\displaystyle=\boldsymbol{F}_{\rm pseudo}$
$\displaystyle~{}+~{}$ $\displaystyle\boldsymbol{F}_{\rm marco}$
$\displaystyle~{}+~{}$ $\displaystyle\boldsymbol{F}_{\rm micro}$ (18)
$\displaystyle=-m\frac{\mathrm{d}\boldsymbol{\varv}}{\mathrm{d}t}-\boldsymbol{(\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla})}\boldsymbol{\varv}$
$\displaystyle~{}+~{}$ $\displaystyle
q\frac{\boldsymbol{\varv}_{\mathrm{cr}}\boldsymbol{\times}\boldsymbol{B}}{c}$
$\displaystyle~{}+~{}$ $\displaystyle
q\left(\delta\boldsymbol{E}+\frac{\boldsymbol{\varv}_{\mathrm{cr}}\boldsymbol{\times}\delta\boldsymbol{B}}{c}\right),$
(19)

where $m$ and $q$ are the CR mass and charge, $\boldsymbol{B}$ is the local
mean magnetic field, and ${\rm d}/{\rm d}t=\partial/\partial
t+\boldsymbol{\varv}\boldsymbol{\cdot}\boldsymbol{\nabla}$ denotes the
Lagrangian time derivative. The first pseudo force arises from the
accelerating comoving frame itself. From the perspective of a comoving
observer, a CR at rest in the laboratory frame appears to be accelerating. The
second pseudo force stems from spatial variations in the velocity field of the
background plasma. If a CR moves in the laboratory frame, any change in its
position leads to a change in the frame velocity due to the inhomogeneities.
This relationship between the laboratory and comoving frames results in
acceleration in the comoving frame. The small-scale fluctuations, denoted as
$\delta\boldsymbol{E}$ and $\delta\boldsymbol{B}$, represent electric and
magnetic fluctuations, respectively. These fluctuations are generated by
plasma waves, especially Alfvén waves on MHD scales and ion-cyclotron waves
that are comoving with the CRs on smaller scales, and they serve as a source
of scattering for CR (see Fig. 3):

$\displaystyle\left.\frac{\partial f}{\partial t}\right\rvert_{\rm
scatt}=\boldsymbol{F}_{\rm
micro}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}f.$ (20)

The waves are either provided by external turbulence and cascading from large
scales or generated by the CR-driven gyroresonant instability.

The CR proton gyroradius is
$r_{\mathrm{g}}=p_{\perp}c/(eB)=0.22~{}\mathrm{AU}\,E_{\mathrm{GeV}}\,B_{\mu\mathrm{G}}^{-1}$
and is much smaller than any macroscopic scale in galaxies or galaxy clusters.
Hence, we can project out the full phase dynamics of CRs by taking their
gyroaverage to arrive at the focused transport equation (Skilling, 1971; Zank,
2014), which is expressed in the reference frame that moves with the average
velocity of the gas

$\displaystyle\frac{\partial f}{\partial t}$
$\displaystyle+(\boldsymbol{\varv}+\mu\varv_{\mathrm{cr}}\boldsymbol{b})\boldsymbol{\cdot}\boldsymbol{\nabla}f+\left[\frac{1-3\mu^{2}}{2}(\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varv}\boldsymbol{\cdot}\boldsymbol{b})-\frac{1-\mu^{2}}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varv}\right]p\frac{\partial
f}{\partial p}$
$\displaystyle+\left[\varv_{\mathrm{cr}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{b}+\mu\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varv}-3\mu(\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varv}\boldsymbol{\cdot}\boldsymbol{b})\right]\frac{1-\mu^{2}}{2}\frac{\partial
f}{\partial\mu}=\left.\frac{\partial f}{\partial t}\right\rvert_{\rm scatt}.$
(21)

Here, the mean gas velocity $\boldsymbol{\varv}$ and the direction of the mean
magnetic field on large scales, $\boldsymbol{b}=\boldsymbol{B}/B$ are measured
in the laboratory frame while the magnitude of the CR momentum,
$p=|\boldsymbol{p}|$, and the cosine of the pitch angle,
$\mu=\boldsymbol{\varv}_{\mathrm{cr}}\boldsymbol{\cdot}\boldsymbol{b}/\varv_{\mathrm{cr}}$
are given in the comoving gas frame. As a result, the CR distribution function
depends on 6 variables, $f=f(\boldsymbol{x},p,\mu,t)$. If the CR distribution
function is characterized by a high degree of anisotropy (which is the case
for CR escape ahead of shocks), we would have to directly solve the focused
transport equation (21) without making further approximations. Here, we define
the CR anisotropy as the deviation of the CR momentum-space distribution from
an isotropic distribution, $\delta f/f_{0}$, where
$f_{0}=f_{0}(\boldsymbol{x},p,t)$ is the isotropic part of the CR distribution
function.

####  2.3.2 One-moment cosmic-ray hydrodynamics

In our Galaxy, the CR anisotropy at particle energies in the range of GeV to
TeV is observed to be very small of order $\mathcal{O}(10^{-4})$ (Kulsrud,
2005). The simplest approximation for CR transport is the one-moment method,
which approximates $f$ by the isotropic part of the CR distribution function
and is derived in quasi-linear theory of CR transport that assumes
electromagnetic fluctuations with a small amplitude relative to the mean
field. Taking the zeroth $\mu$-moment of the focused transport equation (21),
yields a Fokker-Planck equation for CR transport (Skilling, 1971, 1975;
Schlickeiser, 2002),

$\frac{\partial f_{0}}{\partial
t}+(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}})\boldsymbol{\cdot}\boldsymbol{\nabla}f_{0}=\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\kappa\boldsymbol{b}\left(\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}f_{0}\right)\right]+\frac{1}{3}p\frac{\partial
f_{0}}{\partial
p}\,\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}})+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}\Gamma_{\mathrm{p}}\,\frac{\partial f_{0}}{\partial p}\right],$
(22)

where $p=|\boldsymbol{p}|$ is the magnitude of the momentum,
$\Gamma_{\mathrm{p}}$ is the momentum diffusion rate related to second-order
Fermi acceleration, and we have omitted CR sources and non-adiabatic CR
losses. $\boldsymbol{\varv}_{\mathrm{st}}$ and $\kappa$ are the CR streaming
speed and spatial diffusion coefficient along the mean magnetic field, which
will be introduced in detail below. This equation holds true only for time
periods significantly longer than the relaxation time for pitch angle
scattering $\tau\sim\mathcal{O}(D_{\mu\mu}^{-1})$ where
$D_{\mu\mu}=D_{\mu\mu}(\boldsymbol{x},p,\mu)$ is the Fokker-Planck coefficient
representing the frequency of pitch angle scattering of CRs by hydromagnetic
waves.

Multiplying the Fokker-Planck equation (22) by the CR kinetic energy and
integrating it over momentum space results in the evolution equation for the
CR energy density, $\varepsilon_{\mathrm{cr}}$:

$\displaystyle\frac{\partial\varepsilon_{\mathrm{cr}}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left[P_{\mathrm{cr}}\boldsymbol{\varv}_{\mathrm{st}}+\varepsilon_{\mathrm{cr}}(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}}+\boldsymbol{\varv}_{\mathrm{di}})\right]$
$\displaystyle=$ $\displaystyle-
P_{\mathrm{cr}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}}\boldsymbol{\cdot}\boldsymbol{\nabla}P_{\mathrm{cr}},$
(23) $\displaystyle\frac{\partial\varepsilon}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left[(\varepsilon+P_{\mathrm{th}}+P_{\mathrm{cr}})\boldsymbol{\varv}\right]$
$\displaystyle=$
$\displaystyle\phantom{-}P_{\mathrm{cr}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varv}-\boldsymbol{\varv}_{\mathrm{st}}\boldsymbol{\cdot}\boldsymbol{\nabla}P_{\mathrm{cr}}.$
(24)

Note that the CR equation (23) resembles the evolution equation for the
thermal and kinetic energy density,
$\varepsilon=\varepsilon_{\mathrm{th}}+\rho\varv^{2}/2$ shown in equation
(24). The equations of state relate the thermal and CR pressures,
$P_{\mathrm{th}}$ and $P_{\mathrm{cr}}$, to the corresponding energy densities
via

$\displaystyle P_{\mathrm{th}}$ $\displaystyle=(\gamma_{\rm th}$
$\displaystyle-1)\varepsilon_{\rm th},\hskip 40.0pt$ $\displaystyle\gamma_{\rm
th}$ $\displaystyle=\frac{5}{3},$ (25) $\displaystyle P_{\mathrm{cr}}$
$\displaystyle=(\gamma_{\rm cr}$ $\displaystyle-1)\varepsilon_{\mathrm{cr}},$
$\displaystyle\gamma_{\rm cr}$ $\displaystyle=\frac{4}{3}.$ (26)

CR streaming and diffusion are characterized by characteristic velocities,
$\boldsymbol{\varv}_{\mathrm{st}}$ and $\boldsymbol{\varv}_{\mathrm{di}}$
(defined below). Note that here we neglected second-order Fermi acceleration
and any non-adiabatic gain and loss terms for the thermal plasma and CRs for
sake of transparency. Combining equations (23) and (24) causes the terms on
the right-hand side to vanish identically. This means that in the absence of
non-adiabatic gain and loss terms, the total energy in form of CRs, thermal
plasma and kinetic energy is conserved. This can be readily inferred by
integrating the sum of equations (23) and (24) over volume and applying Gauss’
theorem. The physics of the CR energy equation (23) can be more transparently
discussed by expanding the $P_{\mathrm{cr}}\boldsymbol{\varv}_{\mathrm{st}}$
term in the brackets to arrive at

$\displaystyle\frac{\partial\varepsilon_{\mathrm{cr}}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\varepsilon_{\mathrm{cr}}(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}}+\boldsymbol{\varv}_{\mathrm{di}})\right]=-P_{\mathrm{cr}}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}}).$
(27)

This equation states that $\varepsilon_{\mathrm{cr}}$ is advected with the
total CR velocity,
$\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}}+\boldsymbol{\varv}_{\mathrm{di}}$.
CRs can either adiabatically gain or lose energy, depending on the sign of the
total velocity divergence: for adiabatic compression,
$\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varv}+\boldsymbol{\varv}_{\mathrm{st}})<0$,
so that CRs gain energy.

This surprising simplicity of CR transport can be understood by considering
the kinetic physics of CR-wave interactions. Because the magnetic field is
flux frozen into the ionized plasma in the MHD approximation, CRs are advected
with the field lines as they are themselves advected with the fluid velocity
$\boldsymbol{\varv}$. When CRs propagate along the magnetic field, their
microphysical transport depends crucially on the scattering frequency. A small
but non-negligible scattering rate implies diffusive transport along the
magnetic field. In addition to diffusion, CR particles can collectively drift
at an average velocity: when CRs move at a speed exceeding the local Alfvén
speed, they generate gyroresonant Alfvén waves through a phenomenon known as
the streaming instability (Kulsrud & Pearce, 1969). In turn, these waves
increase the CR scattering (Fig. 3), which decelerates the CR population along
the mean magnetic field, effectively transferring CR energy to further growing
resonant Alfvén waves. This self-reinforcing feedback loop is stopped when CRs
flow on average at the speed of Alfvén waves because there is no more net
energy transfer from CRs to Alfvén waves.
Wave damping processes transfer wave energy originally borrowed from the CRs
to the background thermal plasma. This means that CRs exert pressure on the
thermal plasma by scattering of off Alfvén waves. In the presence of weak wave
damping, the high (remaining) wave amplitudes enable a strong coupling,
causing CRs to stream alongside Alfvén waves. Conversely, when wave damping is
strong, wave amplitudes are reduced to the level that they cannot maintain a
frequent scattering that isotropizes CRs in the wave frame, so that the CRs
preferentially diffuse. This self-limiting picture of CR transport is
confirmed by recent simulations using the PIC technique (Holcomb & Spitkovsky,
2019; Shalaby et al., 2021), hybrid PIC (Weidl et al., 2019; Haggerty et al.,
2019) as well as MHD PIC models (Lebiga et al., 2018; Bai et al., 2019;
Plotnikov et al., 2021; Bambic et al., 2021). The highly idealized pure CR
transport modes advection, diffusion and streaming are visualized in Fig. 9.

Figure 9: Left: The magnetic fields, which are frozen-in, are carried along
with the plasma at an average velocity, $\varv=\varv_{\mathrm{adv}}$. Hence,
CRs orbiting individual field lines are also advected alongside the plasma.
Middle: When the amplitude of Alfvén waves is strongly damped, CRs are weakly
scattered and (after an initial time) diffuse away from the source. They are
transported with a characteristic root-mean-square velocity of
$\sqrt{2\kappa/t}$, where $t$ is the time and $\kappa$ denotes the diffusion
coefficient along the magnetic field. Right: If Alfvén waves are weakly
damped, their amplitude remains large so that CRs are effectively scattered
and stream at the Alfvén speed, $\varv_{\mathrm{a}}$, along the orientation of
the magnetic field. These different transport modes of CR fluids are not
realized in their pure forms in Nature as they are instead superposed onto
each other. Image from Thomas et al. (2020); reproduced with permission from
ApJ.

To concretize this qualitative discussion, we define the pitch-angle averaged
CR scattering frequencies of right- and left-ward propagating Alfvén waves
using their energy densities $\varepsilon_{\mathrm{a},\pm}$ (Thomas &
Pfrommer, 2019):

$\displaystyle\bar{\nu}_{\pm}=\frac{3\pi\,\Omega}{16}\frac{\varepsilon_{\mathrm{a},\pm}}{\varepsilon_{B}}=\frac{c^{2}}{3\kappa_{\pm}},$
(28)

where $\varepsilon_{B}=B^{2}/(8\pi)$ is the magnetic energy density and
$\kappa_{\pm}$ denote the CR diffusion coefficients associated with the
corresponding scattering frequencies. Hence, a larger level of Alfvénic
fluctuations increases the scattering rate. CRs stream down their own pressure
gradient relative to the background plasma along the local direction of B with
a velocity

$\displaystyle\boldsymbol{\varv}_{\mathrm{st}}=\boldsymbol{\varv}_{\mathrm{a}}\,\frac{\bar{\nu}_{+}-\bar{\nu}_{-}}{\bar{\nu}_{+}+\bar{\nu}_{-}}\to-\boldsymbol{\varv}_{\mathrm{a}}\,\mathrm{sgn}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}P_{\mathrm{cr}}).$
(29)

The limit on the right-hand side is realized in the self-confinement picture
of CR transport where CRs streaming down the CR gradient excite resonant
Alfvén waves. Those waves are moving in the direction of the streaming CRs so
that their wave amplitudes dominate over the counter propagating type of
waves, which are efficiently damped through various collisionless damping
processes, see Section 2.3.5. If $\boldsymbol{B}$ points in the opposite
direction of $\boldsymbol{\nabla}P_{\mathrm{cr}}$, equation (29) implies
streaming along the magnetic field,
$\boldsymbol{\varv}_{\mathrm{st}}=\boldsymbol{\varv}_{\mathrm{a}}$ and a
dominating scattering rate $\bar{\nu}_{+}$. Efficient scattering of CRs
isotropizes them in the Alfvén wave rest frame (see Section 2.1). As discussed
above, CRs can also diffuse in the wave frame along the local direction of B
due to less frequent pitch angle scattering by MHD waves (Thomas & Pfrommer,
2019):

$\displaystyle\boldsymbol{\varv}_{\mathrm{di}}=-\kappa\textbf{{b}}\,\frac{\textbf{{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\varepsilon_{\mathrm{cr}}}{\varepsilon_{\mathrm{cr}}},\qquad\mbox{where}\qquad\kappa=\frac{c^{2}}{3(\bar{\nu}_{+}+\bar{\nu}_{-})}$
(30)

is the total CR diffusion coefficient, which deceases for an increased CR-wave
scattering rate. Because $\boldsymbol{\varv}_{\mathrm{di}}$ depends on the
gradient of the CR energy density, this type of transport process spreads the
initial CR distribution in time. As a result, the diffusion operator is
characterized by a Laplacian (combining equations 27 and 30 and assuming
$\kappa=\mathrm{const}$). In one dimension, the solution of the diffusion
equation is given by convolving the initial condition with a Gaussian
representing the Green’s function of the diffusion operator. This causes the
diffusive spread of CRs away from a local maximum of CRs along the magnetic
field.

The first numerical simulations to include dynamical CR feedback on the
hydrodynamics studied the problem of diffusive shock acceleration at a SN
blast wave (Dorfi, 1984, 1985) and employed the so-called “two-fluid” model
laid out in equations (23) and (24), which are closed by the equations of
state with the corresponding ratio of specific heats, (25) and (26), for the
thermal and CR populations (Drury & Voelk, 1981; Kang & Jones, 1990). These
early studies assumed simplified expressions for the spatial and momentum
dependence of the CR diffusion coefficient and neglected the CR streaming
term. Using a simplified analytic two-fluid model of the SN expansion, Drury
et al. (1989) account for dynamical CR feedback as well as the influence of
Alfvén wave heating on the diffusion coefficient. The first fully time-
dependent and non-linear numerical simulations of spherical SN blast waves
with the “two-fluid” model of CR-hydrodynamics find that the shock can
transfer $\sim 10\%$ of its kinetic energy to the CR population at the end of
the adiabatic, Sedov-Taylor phase (Jones & Kang, 1990; Dorfi, 1990). The first
CR hydrodynamical simulations coupled to an equation for the resonant Alfvén
wave energy demonstrate the emergence of CR-driven galactic winds with one-
dimensional flux-tube models (Breitschwerdt et al., 1991, 1993).

In the new century, new computational capabilities enabled the study of the
dynamical effect of CR pressure and heating in galaxy formation and galaxy
clusters, using three-dimensional (magneto-)hydrodynamical simulations.
Pioneering work was performed with the Eulerian MHD code PIERNIK that followed
anisotropic CR diffusion (Hanasz & Lesch, 2003) and succeeded by simulations
of galaxies, clusters and cosmological volumes with the smoothed particle
hydrodynamics code Gadget, which evolves a CR population described by a single
power-law momentum spectrum while accounting for isotropic CR diffusion
(Pfrommer et al., 2006; Enßlin et al., 2007; Jubelgas et al., 2008) and CR
streaming (Uhlig et al., 2012). The promising results of galactic winds driven
by CR streaming were confirmed by several other groups that modelled isotropic
CR diffusion using adaptive mesh-refinement codes Ramses (Booth et al., 2013),
Enzo (Salem & Bryan, 2014), and ART (Semenov et al., 2021). Studies of
anisotropic CR diffusion coupled to MHD were made possible by CR
implementations in the Eulerian Pencil code (Snodin et al., 2006), the
adaptive mesh-refinement codes Flash (Yang et al., 2012; Girichidis et al.,
2014) and Ramses (Dashyan & Dubois, 2020), and the unstructured moving-mesh
code AREPO (Pakmor et al., 2016b; Pfrommer et al., 2017a). Anisotropic CR
streaming in MHD simulations was modelled by Ruszkowski et al. (2017a, with
Flash), Butsky & Quinn (2018, with Enzo), and Dubois et al. (2019, with
Ramses) while employing the regularization of the anisotropic CR streaming
term proposed by Sharma et al. (2010a).

####  2.3.3 Two-moment cosmic-ray hydrodynamics

The previously discussed picture of one-moment CR hydrodynamics has the
disadvantage that the (direction of the) CR streaming is ill-defined at the
extremes of the CR distribution (see equation 29). This causes numerical
instabilities of the CR streaming term, which can be cured by regularizing the
sharp transition of the sign function and smoothing it out (Sharma et al.,
2010a), but at the expense of introducing a numerical parameter that dominates
the solution in case of a CR source on top of a CR background (see figure 6 of
Thomas & Pfrommer, 2019). Capitalizing on analogies of CR and radiation
hydrodynamics, Jiang & Oh (2018) propose a two-moment CR hydrodynamics scheme
that assumes steady-state CR scattering approximation and additionally solves
for the CR flux density as an independent and evolved quantity, thereby
solving the problem of non-uniqueness of the CR streaming flux at the extremes
of $f$ (see also Chan et al., 2019, for a similar development).

In an alternative approach, we can expand $f$ into a complete set of basis
functions in pitch-angle (such as the Legendre polynomials, see Section 2.1),
which has a long history in CR transport (see e.g. Klimas & Sandri, 1971;
Earl, 1973; Webb, 1987; Zank et al., 2000; Snodin et al., 2006; Litvinenko &
Noble, 2013; Rodrigues et al., 2019). Using the Eddington approximation for
the CR distribution, $f=f_{0}+3\mu f_{1}$, and requiring that the isotropic
term dominates over the anisotropy, $f_{0}\gg f_{1}$, Thomas & Pfrommer (2019)
compute the zeroth and first $\mu$-moments of the focused transport equation
(21) to arrive at evolution equations for the isotropic and anisotropic parts
of the CR distribution, $f_{0}$ and $f_{1}$. Multiplying those with the
particle energy and the energy flux, respectively, and integrating the
equations for $f_{0}$ and $f_{1}$ over momentum space yields the evolution
equations for the CR energy and momentum density, $\varepsilon_{\mathrm{cr}}$
and $\boldsymbol{f}_{\mathrm{cr}}/c^{2}$. In the laboratory frame, the
pressure, energy and momentum densities transform as follows (neglecting the
relativistic corrections to the CR energy and pressure):

$\displaystyle\textbf{{P}}_{\mathrm{cr},\mathrm{lab}}=P_{\mathrm{cr}}\textbf{{1}},\quad\varepsilon_{\mathrm{cr},\mathrm{lab}}=\varepsilon_{\mathrm{cr}},\quad\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=\boldsymbol{f}_{\mathrm{cr}}+\boldsymbol{\varv}\left(\varepsilon_{\mathrm{cr}}+P_{\mathrm{cr}}\right),$
(31)

so that we obtain the governing equations for the CR energy and momentum
density as measured in the laboratory frame (Thomas & Pfrommer, 2019, see
Appendix E):

$\displaystyle\frac{\partial\varepsilon_{\mathrm{cr},\mathrm{lab}}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=$
$\displaystyle-\boldsymbol{\varw}_{+}\boldsymbol{\cdot}$
$\displaystyle\frac{\boldsymbol{b}\boldsymbol{b}}{3\kappa_{+}}\boldsymbol{\cdot}\left[\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}-\boldsymbol{\varw}_{+}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})\right]$
(32) $\displaystyle-\boldsymbol{\varw}_{-}\boldsymbol{\cdot}$
$\displaystyle\frac{\boldsymbol{b}\boldsymbol{b}}{3\kappa_{-}}\boldsymbol{\cdot}\left[\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}-\boldsymbol{\varw}_{-}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})\right]$
$\displaystyle-\boldsymbol{\varv}\boldsymbol{\cdot}\boldsymbol{g}_{\mathrm{Lorentz}}$
$\displaystyle+S_{\varepsilon},$ $\displaystyle\frac{1}{c^{2}}$
$\displaystyle\frac{\partial\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\textbf{{P}}_{\mathrm{cr},\mathrm{lab}}=$
$\displaystyle\,-$
$\displaystyle\frac{\boldsymbol{b}\boldsymbol{b}}{3\kappa_{+}}\boldsymbol{\cdot}\left[\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}-\boldsymbol{\varw}_{+}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})\right]$
(33) $\displaystyle\,-$
$\displaystyle\frac{\boldsymbol{b}\boldsymbol{b}}{3\kappa_{-}}\boldsymbol{\cdot}\left[\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}-\boldsymbol{\varw}_{-}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})\right]$
$\displaystyle-\boldsymbol{g}_{\mathrm{Lorentz}}$
$\displaystyle+\boldsymbol{S}_{f},$

where the Alfvén wave velocity in the laboratory frame is given by
$\boldsymbol{\varw}_{\pm}=\boldsymbol{\varv}\pm\boldsymbol{\varv}_{\mathrm{a}}$,
$S_{\varepsilon}$ and $\boldsymbol{S}_{f}$ are non-adiabatic source terms of
CR energy and flux, respectively, and the Lorentz force density in the ultra-
relativistic limit is $c^{2}\boldsymbol{g}_{\rm
Lorentz}=\boldsymbol{f}_{\mathrm{cr}}\boldsymbol{\times}\boldsymbol{\Omega}$,
where $\boldsymbol{\Omega}=\Omega\boldsymbol{b}$.

These equations also contain the previously discussed three CR transport
modes: advection, streaming, and diffusion. In the absence of sources and
along the magnetic field (where $\boldsymbol{g}_{\rm
Lorentz}=\boldsymbol{0}$), we recover CR streaming and diffusion. In the
presence of efficient CR-wave scattering, i.e., for a large CR scattering
frequency $\bar{\nu}_{\pm}=c^{2}/(3\kappa_{\pm})$, we approach a homogeneous
CR flux with a steady-state CR energy density, implying
$\partial_{t}\varepsilon_{\mathrm{cr},\mathrm{lab}}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=0$.
In consequence, the CR enthalpy streams with the Alfvén wave frame,
$\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=\boldsymbol{\varw}_{\pm}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})$,
where the sign of $\boldsymbol{\varw}_{\pm}$ determines the CR propagation
direction. If CR scattering is inefficient but finite, we also approach a
steady state CR flux
($\partial_{t}\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=\boldsymbol{0}$) and
recover the CR diffusion velocity,
$\boldsymbol{\varv}_{\mathrm{di}}\varepsilon_{\mathrm{cr},\mathrm{lab}}=-\kappa_{\pm}\boldsymbol{b}\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\varepsilon_{\mathrm{cr},\mathrm{lab}}$,
where we have identified the non-vanishing bracket on the right-hand side with
$\boldsymbol{\varv}_{\mathrm{di}}\varepsilon_{\mathrm{cr},\mathrm{lab}}$ and
used the ultra-relativistic equation of state for CRs,
$P_{\mathrm{cr}}=\varepsilon_{\mathrm{cr}}/3$. Finally, we recover CR
advection in steady state,
$\boldsymbol{f}_{\mathrm{cr},\mathrm{lab}}=\boldsymbol{\varv}(\varepsilon_{\mathrm{cr},\mathrm{lab}}+P_{\mathrm{cr},\mathrm{lab}})$
due to the balance of forces parallel and perpendicular to the magnetic field:
while the perpendicular pressure force balances the Lorentz force acting on
CRs,
$\boldsymbol{\nabla}_{\perp}P_{\mathrm{cr},\mathrm{lab}}=-\boldsymbol{g}_{\mathrm{Lorentz}}$
(which may however act on kinetic time scales), we have forces due to CR
streaming and diffusion balancing the pressure forces parallel to the magnetic
field.

These equations closely resemble the lab-frame equations for radiation energy
and momentum density, $\varepsilon$ and $\boldsymbol{f}/c^{2}$ (Mihalas &
Mihalas, 1984; Lowrie et al., 1999)

$\displaystyle\frac{\partial\varepsilon}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}$
$\displaystyle=-\sigma_{\mathrm{s}}\,\boldsymbol{\varv}\,\boldsymbol{\cdot}\,$
$\displaystyle\left[\boldsymbol{f}-\boldsymbol{\varv}\boldsymbol{\cdot}(\varepsilon\textbf{{1}}+\textbf{{P}})\right]$
$\displaystyle+S_{\mathrm{m}},$ (34) $\displaystyle\frac{1}{c^{2}}$
$\displaystyle\frac{\partial\boldsymbol{f}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\textbf{{P}}$
$\displaystyle=-\sigma_{\mathrm{s}}$
$\displaystyle\left[\boldsymbol{f}-\boldsymbol{\varv}\boldsymbol{\cdot}(\varepsilon\textbf{{1}}+\textbf{{P}})\right]$
$\displaystyle+S_{\mathrm{m}}\boldsymbol{\varv},$ (35)

where P is the radiation pressure tensor, $\sigma_{\mathrm{s}}$ is the
scattering coefficient, and $S_{\mathrm{m}}$ describes photon absorption and
emission processes, thus coupling radiation to matter. The similarities of the
CR and photon transport equations are apparent and result from the same
underlying physics: while efficient photon scattering causes the photon
enthalpy to be advected with the gas, efficient CR-wave scatterings causes the