# Reacceleration of Galactic Cosmic Rays Beyond the Knee at the Termination
Shock of a Cosmic-Ray-Driven Galactic Wind

Payel Mukhopadhyay SLAC National Accelerator Laboratory, Stanford University,
Menlo Park, CA 94025, USA<EMAIL_ADDRESS>Physics Department, Stanford
University, Stanford, CA 94305, USA  Department of Physics, University of
California, Berkeley, CA, 94720, USA  Enrico Peretti Niels Bohr International
Academy, Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17,
DK-2100 Copenhagen, Denmark  Noémie Globus Department of Astronomy and
Astrophysics, University of California, Santa Cruz, CA 95064, USA
Astrophysical Big Bang Laboratory, RIKEN, Wako, Saitama, Japan  Paul Simeon
Independent Researcher  Roger Blandford Physics Department, Stanford
University, Stanford, CA 94305, USA  Kavli Institute for Particle Astrophysics
and Cosmology, 452 Lomita Mall, Stanford, CA 94305, USA

###### Abstract

The origin of cosmic rays above the knee in the spectrum is an unsolved
problem. We present a wind model in which interstellar gas flows along a non-
rotating, expanding flux tube with a changing speed and cross-sectional area.
Cosmic rays from Galactic sources, such as supernova remnants, which are
coupled to the plasma via Alfvén waves, provide the main pressure source for
driving this outflow. These cosmic rays are then subject to diffusive shock
reacceleration at the Galactic wind termination shock, which is located at a
distance $\sim 200\,{\rm kpc}$. Some of the highest-energy reaccelerated
particles propagate upstream against the wind and can contribute to the
PeV–EeV range of the spectrum. We analyze the conditions under which efficient
reacceleration can occur and find that rigidities $\sim$ 10–40 PV can be
obtained and that the termination shock may account for half of the proton
spectrum measured in IceCube/IceTop experiment. The highest-energy particles
that escape downstream from our termination shock, and similar shocks
surrounding most galaxies, can be further accelerated by intergalactic shock
fronts.

cosmic rays, galactic winds, termination shock

##  1 Introduction

The origin of cosmic rays (CRs) is a long-standing open question since their
discovery. The most established interpretation ascribes a Galactic origin to
CRs populating the spectrum up to the “knee”, with rigidity (momentum per unit
charge) $R=R_{15}\,{\rm PV}$ satisfying $R_{15}\sim 3$, where the spectrum
steepens (Tibet AS$\gamma$ Collaboration et al., 2021; Cao et al., 2021;
Albert et al., 2020; Aharonian et al., 2007). By contrast, cosmic rays with
energy above the “ankle” in the spectrum at $R_{15}\gtrsim 3000$, where the
spectrum flattens, are widely believed to have an extragalactic origin because
the average Galactic magnetic field is not strong enough to confine such
energetic particles. Furthermore, the arrival flux does not feature any
relevant anisotropy correlating with the Galactic disk (see e.g., Aab et al.,
2018).

Similar to extragalactic cosmic rays, the origin of cosmic rays between the
knee and the ankle, also referred to as the “shin” of the cosmic-ray spectrum,
is highly uncertain. We call these “intermediate-energy” cosmic rays. The most
widely discussed sources are pulsars, magnetars, circumgalactic and
extragalactic shock fronts, and unidentified “pevatrons” (Abeysekara et al.,
2021a; Cao et al., 2021; Amenomori et al., 2019; HESS Collaboration et al.,
2016; Gottlieb & Globus, 2021).

In this paper, we focus on one particular type of source of high energy CRs
that could possibly populate the shin region of the spectrum: the Galactic
wind termination shock (GWTS). We consider Galactic cosmic rays as being
injected by sources located in the Galactic disk from the lowest energies up
to the knee. These cosmic rays are assumed to escape the disk and be advected
in the Galactic wind on timescales shorter than the Galactic rotation period.
The spectral index of the escaping cosmic rays is essentially that produced by
the Galactic accelerators like supernova remnants (SNRs), not what is observed
at Earth, assuming the distribution of Galactic cosmic rays is approximately
stationary, with the time-averaged source flux equal to the escaping flux. The
Galactic wind advects cosmic rays in an accelerating flow that passes through
a critical point and becomes supersonic. Eventually, this flow will pass
through a strong shock front, behind which it will decelerate and eventually
become incorporated into the intergalactic medium. The shock reaccelerates
Galactic cosmic rays by diffusive shock (re)acceleration (DSA), contributing
to the shin region of the spectrum. Some of these higher-energy particles will
propagate upstream back to the Galactic disk, affecting the spectrum,
composition, and anisotropy observed at Earth.

Of course, all other galaxies, especially disk galaxies, where there is active
star formation and a relatively high supernova rate, should behave similarly
to our Galaxy. The high-energy cosmic rays that are advected downstream from
the termination shock escape their original galaxy and provide an input for
larger shock fronts. Galaxies with active nuclei, notably starburst galaxies,
likely emit an outsized flux of these intergalactic cosmic rays (see e.g.,
Zhang et al., 2020; Peretti et al., 2022). The largest intergalactic shock
fronts, especially those associated with clusters of galaxies and the
filaments that connect them in the “cosmic web,” may be the source of the
extragalactic cosmic rays. This possibility will be discussed in a future
publication (Simeon et al. in preparation).

Galactic wind termination shocks were first proposed as a cosmic ray
acceleration site in the papers by Jokipii & Morfill (1985) and Jokipii &
Morfill (1987). Since then, the idea has been explored by different authors
(Völk & Zirakashvili, 2004; Zirakashvili & Völk, 2006; Thoudam et al., 2016;
Bustard et al., 2017; Merten et al., 2018). The injected cosmic rays that are
accelerated at the GWTS can either come from particles injected at the shock
front itself (Bustard et al., 2017; Merten et al., 2018) or from Galactic
cosmic rays that are transported along the Galactic wind (Zirakashvili & Völk,
2006; Thoudam et al., 2016). Zirakashvili & Völk (2006) and Thoudam et al.
(2016) argued that cosmic rays reaccelerated at the GWTS that propagate back
to the disk can fully explain the observed all-particle spectrum at Earth
between the knee and the ankle region. These studies were limited in their
treatment of the cosmic-ray transport equation. Specifically, the maximum
energies achieved in the reacceleration process were not self-consistently
computed from the transport equation but were prescribed in the solution.
Thoudam et al. (2016) treated the maximum energies as free parameters chosen
such that the spectrum of reaccelerated particles propagating back to the disk
reasonably agrees with the measured all-particle spectrum. Zirakashvili & Völk
(2006) treated the GWTS as a reflecting boundary and defined the maximum
energy as the highest energy beyond which the particles can cross the
termination shock diffusively. Additionally, these studies neglected the
transport equation downstream of the shock and did not explore the possibility
of a fraction of the reaccelerated cosmic rays escaping into the intergalactic
medium, seeding it with intermediate-energy cosmic rays.

In this paper, we consistently solve the CR transport equation both upstream
and downstream of the GWTS. We find the maximum energies of the reaccelerated
cosmic rays self-consistently as a solution to the transport equation, instead
of prescribing it externally. We also calculate the flux of reaccelerated
cosmic rays that propagate back to the Galaxy, as well as the flux of
particles escaping into the intergalactic medium. Additionally, we provide an
analysis of the conditions under which efficient CR reacceleration can happen
at the GWTS. To this end, we explore a physically motivated model of the
Galactic wind where the wind is driven by cosmic rays and the pressure of
self-excited Alfvén waves streaming away from the disk. This model follows the
treatment presented by Breitschwerdt et al. (1991), albeit with an updated
treatment of the dark-matter potential of the Milky Way. We will call this
model, the “Cosmic Ray–Alfvén wave driven” wind model. We then solve the
transport equation by prescribing the velocity profiles obtained from the wind
model. The transport equation is computed for two different prescriptions of
the upstream diffusion coefficient to test the conditions under which
efficient reacceleration can occur and to understand what conditions are
suitable for efficient reacceleration of cosmic rays at the GWTS. In this
model, we neglect the rotation of the disk for computing the wind. In a
follow-up paper (Blandford et al. 2023 in preparation), we will develop our
model, by including rotation through a magnetocentrifugally-driven wind that
becomes a significant part of the mass, energy, and angular momentum budget of
the interstellar medium and which may accelerate cosmic rays more efficiently
to higher energy than in the present model while being subject to significant
observational constraints.

We review the current state of measurements in Sec. 2. We describe a
relatively simple, though quite general, formalism for describing the outflow
in Sec. 3, emphasizing the choice that we argue is most appropriate for the
solar neighborhood. We also describe the transport of Galactic CRs in the wind
and the reacceleration of these CRs at the GWTS. In Sec. 4, we present our
results. Finally, in Sec. 5, we summarize and present our conclusions

##  2 Cosmic-Ray Observations

It is widely believed that cosmic rays below $\sim 10^{15}$ eV are accelerated
by Galactic accelerators, possibly SNRs, and that the ankle of the cosmic-ray
spectrum marks the end of the transition from Galactic to extragalactic cosmic
rays (e.g., Globus et al., 2015). The explanation of the intermediate-energy
cosmic rays is still unsolved but increasingly constrained by the many
detailed observations of the past decade. We summarize some of these
observations below.

##### Spectrum

The AMS-02, CALET, DAMPE, HAWC, and ISS-CREAM experiments have measured the
individual spectra of various intermediate energy species with unprecedented
accuracy (Aguilar et al., 2015, 2017, 2020; Adriani et al., 2019, 2021;
Alemanno et al., 2021; Abeysekara et al., 2021b). These spectra show ankles
and breaks. A spectral softening, or a “knee,” can be caused by the maximum
energy of the acceleration process, a leakage from the Galaxy, or a
combination of both. An “ankle” is a hardening of the cosmic-ray spectrum, a
natural feature marking the transition between two components of the spectrum.

Above a few GeV/nucleon, the energy spectrum follows a power law proportional
to $E^{-2.7}$. A first softening in the spectrum occurs at around $1.2\times
10^{13}$ eV (Alemanno et al., 2021). It is then followed by a hardening at
$\sim 10^{14}$ eV (Abeysekara et al., 2021b). The cosmic-ray spectrum then has
another stronger knee at around $3\times 10^{15}$ eV — a well-known feature of
the cosmic-ray spectrum known since the 1950s. Beyond that energy, the cosmic-
ray spectrum then follows a power law proportional to $E^{-3.0}$. At $\sim
10^{17}$ eV, the spectrum shows another softening only in the heavy component,
the so-called “heavy knee” (Apel et al., 2011). At the same energy of about
$\sim 10^{17}$ eV, the light component (proton and helium) shows a hardening,
a “light ankle” (Apel et al., 2013), that has been interpreted as the
beginning of the light extragalactic component of cosmic rays (Globus et al.,
2015). At $4\times 10^{18}$ eV, the total cosmic-ray spectrum gets harder
again. This final ankle marks the end of the transition between Galactic and
extragalactic cosmic rays.

##### Composition

Several experiments, e.g., ARGO-YBJ, KASCADE-Grande, Tunka, IceCube, and
IceTop (e.g., Andeen & Plum, 2019; Kang et al., 2021), have measured the
composition of cosmic rays between the knee and the ankle of the spectrum. The
average composition is light (predominantly protons and helium nuclei) at the
knee, i.e., at energy $\sim 10^{15.5}$ eV. Above the knee, the average
composition becomes heavier as a function of energy up to about $10^{17}$ eV,
and then, although the statistical errors become significant at these
energies, the average mass seems to become lighter up to the ankle at $\sim 4$
EeV. The different interaction models predict a different average mass, and
the differences increase with primary energy. This trend is naively expected
when the maximum energy attained during the acceleration process is
proportional to the nuclear charge; if the proton knee occurs at few PeV, then
the iron knee would occur at around $10^{17}$ eV, which is consistent with the
evolution of the composition. The heavy knee reported by KASCADE-Grande at
$10^{17}$ eV could mark the end of the heavy component from the SNR
contribution.

##### Anisotropy

Ahlers (2019) reports some hints of intermediate-scale anisotropy at 33 PeV,
although the level of anisotropy in the PeV–EeV energy range is low. The upper
limit to the amplitude of the first harmonic (dipole) is $10^{-3}$ at $2\times
10^{14}$ eV and increases to $\sim 10^{-2}$ at $10^{17}$ eV. KASCADE-Grande,
ESA-TOP, IceCube, and IceTop observed a change of the phase of the first
harmonic in the direction of the Galactic center at energies of $2\times
10^{14}$ eV. After that, the phase remains roughly flat until $\sim 10^{17}$
eV and then changes again (to $\sim 100$ degrees). The first phase change
could indicate a transition between local sources located in the neighboring
arms to a more global contribution of the whole Galaxy, while the second phase
change could be the sign of an extragalactic origin of cosmic rays.

##### Impact of reaccelerated CRs on the observables

In this paper, we will compute the flux of Galactic cosmic rays that are
reaccelerated at the GWTS and propagate back to the Galactic disk. We will
compare how the spectrum of the reaccelerated CRs compare with the measured
spectrum at the disk beyond the knee. We will see that these backstreaming
particles can form spectral bump features as a mark of transition between the
Galactic and the GWTS components. We will find that these bump features are
consistent with existing measurements of different elements and that they can
be further constrained in future by more precise experiments.

##  3 Cosmic ray–Alfvén wave driven wind model

###  3.1 General Considerations

Galactic winds driven by thermal pressure, cosmic rays, or magnetic fields
have been studied for more than 50 years (Burke, 1968; Johnson & Axford, 1971;
Ipavich, 1975; Chevalier & Clegg, 1985; Breitschwerdt et al., 1991; Everett et
al., 2008; Recchia et al., 2016; Mao & Ostriker, 2018). If such a wind for the
Galaxy were to be thermally driven, the sound speed at the solar circle would
be $v_{\rm s}\sim\sqrt{\frac{2}{3}}v_{\rm esc}\sim 470~{}{\rm km\,s}^{-1}$,
where $v_{\rm esc}$ is the local escape speed, which is $\sim$580 km s-1 at
the solar circle (Monari et al., 2018; Necib & Lin, 2022). This sound speed
corresponds to a temperature higher than $10^{7}$ K, which is very high
considering the warm interstellar medium (ISM) temperature of the Galaxy is
between $10^{5}$ K and $10^{6}$ K. Therefore, the scenario of a purely
thermally driven wind in the Galaxy is in tension with the observations of the
ISM.

One alternative is that cosmic rays assist in driving the wind, a possibility
that has been discussed by a number of authors (Ipavich, 1975; Breitschwerdt
et al., 1991; Everett et al., 2008; Dorfi & Breitschwerdt, 2012; Recchia et
al., 2016; Dorfi et al., 2019; Mao & Ostriker, 2018). One plausible scenario,
first described by Breitschwerdt et al. (1991), is that the cosmic-ray sources
in the Galactic disk produce energetic particles that cannot freely escape
from the Galaxy but rather amplify Alfvén waves (Wentzel, 1974). Such waves
lead to an efficient coupling of the thermal gas to energetic particles
through the cosmic-ray resonant streaming instability (Lerche, 1967; Kulsrud &
Pearce, 1969), and the pressure gradient of cosmic rays drives a wind outflow.
These winds start at subsonic speeds of a few km s-1 at the base of the wind
and become supersonic at $\mathcal{O}(10~{}\mathrm{kpc})$ distances. The flow
forms a wind termination shock at $\mathcal{O}(10^{2}~{}\mathrm{kpc})$
distances (Breitschwerdt et al., 1991). We adopt the simple model of the wind
used by Breitschwerdt et al. (1991) and modify it by including an updated
dark-matter potential of the Galaxy, consistent with GAIA observations.

###  3.2 Wind geometry and dynamics

Our reference wind model treats cosmic rays as a rarefied plasma in dynamic
interaction with the magnetic field in the interstellar medium. This
approximation holds because rapid pitch-angle scattering of cosmic rays on the
fluctuations (waves) of the field makes the cosmic-ray momentum distribution
isotopic to the lowest order in the wave frame. The residual cosmic-ray
diffusion thorough the gas, along the mean magnetic field, adds itself to any
convective motions of the thermal plasma. Propagation of waves along the mean
magnetic field will, in addition, give rise to a drift speed of cosmic rays
equal to the wave phase velocity. The cosmic ray motion is therefore composed
of convection, wave drift, and diffusion.

Figure 1: Sketch of the problem under consideration. The wind geometry is
shown as a magnetic flux tube. Cosmic rays produced by Galactic sources in the
disk travel with the wind up to the Galactic wind termination shock (GWTS),
where they get reaccelerated to higher energies.

When the streaming speed of the cosmic rays through the ambient medium exceeds
the Alfvén speed, $v_{\mathrm{A}}={B}/\sqrt{{4\pi\rho}}$ ($B$ is the magnetic
field and $\rho$ is the ion density of the medium), there is a quasi-linear
instability in which the amplitude of the resonant Alfvén waves grows
exponentially with time until their amplitude is strong enough to limit the
cosmic-ray streaming speed to slightly more than the Alfvén speed (Lerche,
1967). See Blandford & Eichler (1987) for a review. The cosmic rays,
scattering off these resonantly excited Alfvén waves, inevitably transfer
momentum to the background medium. This ensures that the pressure gradient of
the cosmic rays accelerates the background fluid, which accelerates the wind
in our model. This mechanism only works if the growth rate of these Alfvén
waves, $\Gamma_{\rm res}$, is faster than $\frac{v}{z}$, the inverse of the
dynamical time scale of the problem. The growth rate of these waves is

$\Gamma_{\rm res}=\Omega\left[\frac{v_{\rm
stream}}{v_{\mathrm{A}}}-1\right]\frac{n_{\rm cr}}{n_{\mathrm{i}}},$ (1)

where $n_{\rm cr}\sim 4\pi fp^{3}$ is an estimate of the total number of
resonant comsic rays, and $n_{\mathrm{i}}$ is the background ion density
(Blandford & Eichler, 1987). $v_{\rm stream}$ is the cosmic-ray streaming
speed, $v_{A}$ is the Alfvén speed, $\Omega$ is the non-relativistic
gyrofrequency of the cosmic rays, and $n_{\rm cr}$ is the cosmic-ray number
density, which is $\sim 10^{-12}$ cm-3 for GeV cosmic rays, the primary
drivers of the wind. Using $v_{\rm stream}=2~{}v_{\rm A}$, $n_{\rm i}\sim
10^{-4}$ cm-3, and $B\sim 2~{}\mu{\rm G}$, which are order-of-magnitude
numbers from our calculation at $z\sim 10$ kpc, we find that the condition
$\Gamma_{\rm res}\times z/v\sim 10^{6}\gg 1$ is easily satisfied for $O($1
GeV) CRs. Accordingly, the cosmic-ray pressure is sufficient to excite Alfvén
waves to launch the wind, and this instability can grow faster than the
dynamical time ($z/v$) of the problem for cosmic rays of energies up to
$O(1~{}\mathrm{TeV})$. We note these waves could be damped by processes such
as ion–neutral or nonlinear Landau damping. We neglect the former because a
hot outflow plasma with $T\gtrsim 10^{5}$ K can maintain a high level of
ionization. The latter is expected to become important when $\frac{\delta
B}{B}\sim 1$, when quasilinear theory breaks down. In this paper, we neglect
these damping processes and acknowledge the caveat that these processes could
be important.

We assume the magnetic field, which is also the geometry of the flow, has a
“mushroom-type” geometry (Breitschwerdt et al., 1991), which emphasizes the
concept of isolated open magnetic flux tubes, and can be conveniently
represented by a one-dimensional flux tube model, as shown in Fig. 1. Edge-on
observations of spiral galaxies like NGC 5775 reveal large scale magnetic
fields of similar topology in their galactic halo (Soida et al., 2011). 3D
numerical simulations of the magnetic field evolution in barred galaxies under
the influence of cosmic ray driven dynamo find that the halo magnetic field
can have strong vertical components, possibly transported by mass outflow from
the galactic disk (Kulpa-Dybeł et al., 2011).

This formulation accounts for a vertical flow geometry near the disk, which
transitions to a spherical geometry adequate to describe the flow at altitudes
greater than the disk radius, $R_{\mathrm{G}}\sim 15$ kpc. Following this
discussion, the geometry for the system is physically represented by outwardly
directed magnetic lines of force:

$\displaystyle A(z)=A_{0}\left[1+\left(\frac{z}{z_{0}}\right)^{2}\right],$ (2)

where $z$ is the altitude from the disk, and $z_{0}$ is the vertical scale
beyond which the geometry becomes radial. The parameter $z_{0}$ is typically
expected to be of the order of the Galactic radius (Breitschwerdt et al.,
1991; Dorfi & Breitschwerdt, 2012). $A(z)$ represents the area cross section
of the flow as a function of $z$.

The following dynamical equations describing the overall balance of mass,
momentum, and energy are solved (Breitschwerdt et al., 1991):

$\displaystyle\nabla\cdot\left(\rho\textbf{v}\right)=0$ (3)
$\displaystyle\nabla\cdot\left(\rho\textbf{v}\otimes\textbf{v}+\left[P_{\mathrm{g}}+P_{\mathrm{c}}+\frac{\langle\left(\delta\textbf{B}\right)^{2}\rangle}{8\pi}\right]\cdot\textbf{I}\right)=-\rho\nabla\Phi$
(4)
$\displaystyle\nabla\cdot\left(\rho\textbf{v}\left[\frac{1}{2}v^{2}+\frac{\gamma_{\mathrm{g}}}{\gamma_{\mathrm{g}}-1}\frac{P_{\mathrm{g}}}{\rho}+\Phi\right]+\frac{\langle\left(\delta\textbf{B}\right)^{2}\rangle}{4\pi}\left[\frac{3}{2}\textbf{v}+\textbf{v}_{\rm
A}\right]\right.$
$\displaystyle~{}~{}\left.+\frac{1}{\gamma_{\mathrm{c}}-1}\left[\gamma_{\mathrm{c}}P_{\mathrm{c}}(\textbf{v}+\textbf{v}_{\rm
A})-\overline{D}\nabla P_{\mathrm{c}}\right]\right)=0$ (5)
$\displaystyle\nabla\cdot\left(\frac{\gamma_{\mathrm{c}}}{\gamma_{\mathrm{c}}-1}(\textbf{v}+\textbf{v}_{\rm
A})P_{\mathrm{c}}-\frac{\overline{D}}{\gamma_{\mathrm{c}}-1}\nabla
P_{\mathrm{c}}\right)=(\textbf{v}+\textbf{v}_{\rm A})\nabla P_{\mathrm{c}}$
(6)
$\displaystyle\nabla\cdot\left(\frac{\langle\left(\delta\textbf{B}\right)^{2}\rangle}{4\pi}\left[\frac{3}{2}\textbf{v}+\textbf{v}_{\rm
A}\right]\right)=\textbf{v}\nabla\left(\frac{\langle\left(\delta\textbf{B}\right)^{2}\rangle}{8\pi}\right)-\textbf{v}_{\rm
A}\nabla P_{\mathrm{c}}$ (7) $\displaystyle\nabla\cdot\textbf{B}=0.$ (8)

Here, $\rho$ denotes the gas density, v denotes the outflow gas velocity
vector, $\textbf{v}_{\rm A}$ is the Alfvén speed, B is the mean magnetic
field, $P_{\mathrm{c}}$ and $P_{\mathrm{g}}$ are the pressures of the cosmic
rays and gas, respectively, and $\Phi$ is the gravitational potential.
$\langle\left(\delta\textbf{B}\right)^{2}\rangle/(8\pi)$ is the pressure
contribution of the magnetic field fluctuations. These fluctuations are
generated via the resonant streaming instability of cosmic rays streaming at
roughly the Alfvén speed with respect to the background medium. In the rest of
the text, the pressure contribution due to magnetic field fluctuations will be
denoted as
$P_{\mathrm{w}}=\langle\left(\delta\textbf{B}\right)^{2}\rangle/(8\pi)$. The
adiabatic index of the thermal gas is $\gamma_{\mathrm{g}}=5/3$, and the
adiabatic index of the cosmic rays is $\gamma_{\mathrm{c}}=4/3$. The first
three equations (3, 4, and 5) are the overall conservation of mass, energy,
and momentum. Eq. 6 is the equation for the energy balance of the cosmic ray
component. Eq. 7 is the energy balance of the mean-squared fluctuating waves
that are assumed to propagate down the cosmic-ray gradient with the Alfvén
speed, and finally Eq. 8 states that the divergence of the mean magnetic field
is zero. $\overline{D}$ in Eqs. 5 and 6 represents the effective diffusion
coefficient of the cosmic rays. We neglect the diffusion term for finding the
wind solutions and assume that the transport is advection dominated for the
GeV cosmic rays, providing most of the pressure support for the wind. Finally,
the symbol $\otimes$ is the tensor product, and I is the unit tensor.

With a flux tube geometry, the divergence operator becomes $\nabla\cdot
x=\frac{1}{A}\frac{\differential(A\cdot x)}{\differential z}$, where $A$ is
the cross-section area function and $z$ is altitude. We define the Alfvén Mach
number as $M_{\mathrm{A}}=\frac{v}{v_{\mathrm{A}}}$. We simplify the above set
of equations as follows:

$\displaystyle\rho vA=\mathrm{const}$ (9) $\displaystyle\frac{\differential
P_{\mathrm{g}}}{\differential
z}=\gamma_{\mathrm{g}}\frac{P_{\mathrm{g}}}{\rho}\frac{\differential\rho}{\differential
z}$ (10) $\displaystyle\frac{\differential P_{\mathrm{c}}}{\differential
z}=\gamma_{\mathrm{c}}\frac{P_{\mathrm{c}}}{\rho}\left(\frac{M_{\mathrm{A}}+\frac{1}{2}}{M_{\mathrm{A}}+1}\right)\frac{\differential\rho}{\differential
z}$ (11) $\displaystyle\frac{\differential P_{\mathrm{w}}}{\differential
z}=\frac{1}{2(M_{\mathrm{A}}+1)}\left[(3M_{\mathrm{A}}+1)\frac{P_{\mathrm{w}}}{\rho}\frac{\differential\rho}{\differential
z}-\frac{\differential P_{\mathrm{c}}}{\differential z}\right]$ (12)
$\displaystyle BA=\mathrm{const}$ (13) $\displaystyle\frac{\differential
v}{\differential z}=v\frac{~{}c^{2}_{*}\frac{1}{A}\frac{\differential
A}{\differential z}-\frac{\differential\Phi}{\differential
z}}{v^{2}-c^{2}_{*}},$ (14)

where we have defined a so-called composite sound speed $c_{*}$ defined to be
the following (Breitschwerdt et al., 1991):

$c^{2}_{*}=\gamma_{\mathrm{g}}\frac{P_{\mathrm{g}}}{\rho}+\gamma_{\mathrm{c}}\frac{P_{\mathrm{c}}}{\rho}\frac{(M_{\mathrm{A}}+\frac{1}{2})^{2}}{(M_{\mathrm{A}}+1)^{2}}+\frac{P_{\mathrm{w}}}{\rho}\frac{3M_{\mathrm{A}}+1}{2(M_{\mathrm{A}}+1)}.$
(15)

The gravitational potential $\Phi$ is modeled by taking inputs from GAIA data
(Monari et al., 2018; Necib & Lin, 2022). The dark-matter halo potential is
taken to be

$\Phi_{\mathrm{halo}}=\frac{v_{0}^{2}}{2}\ln{\left(1+\left[\frac{z^{2}+R_{0}^{2}}{R_{\rm
c}^{2}}\right]\right)},$ (16)

with the characteristic speed $v_{0}\sim 180$ km s-1, a radial concentration
parameter $R_{\rm c}\sim 5$ kpc (Monari et al., 2018), and the galactocentric
distance $R_{0}$, which is $\sim$ 8.5 kpc near the solar circle. Since
$\Phi_{\mathrm{halo}}$ does not go to 0 for large $z$, the escape speed for
this potential is computed by Monari et al. (2018):

$v_{\mathrm{esc}}(r)=\sqrt{2|\Phi(r)-\Phi(3r_{340})|},$ (17)

where $r_{340}$ is the spherical radius within which the average density of
the whole Galaxy is 340 times the critical density at redshift 0
$\left(\rho_{\rm c}=\frac{3H^{2}}{8\pi G}\right)$. This definition of the
escape speed physically means that when a particle is at large distances
$r\gtrsim 3r_{340}$, the particle has essentially ‘escaped’ because the Galaxy
is not an isolated system, and at those radii the gravitational potentials of
other galaxies, like Andromeda, start to dominate (Piffl et al., 2014; Monari
et al., 2018). The escape speed as a function of disk height $z$ at
$R_{0}=8.5$ kpc is shown in Fig. 2. Near the solar circle, the escape speed in
our model is $\sim 580$ km s-1, which is consistent with GAIA data (Monari et
al., 2018; Necib & Lin, 2022).

Figure 2:  Milky Way escape-speed profile used in this study. The halo
potential assumed is given in Eq. (16). The bulge and disk potential are
modeled according to Eq. (18). The escape speed near the solar circle is
$\sim$ 580 km/s, consistent with GAIA data.

The gravitational potential of the bulge–disk component is modelled following
Miyamoto and Nagai prescription (Miyamoto & Nagai, 1975):

$\Phi_{\rm B,D}(R_{0},z)=-\sum_{i=1}^{2}\frac{GM_{\rm
i}}{\sqrt{R_{0}^{2}+\left(a_{\rm i}+\sqrt{z^{2}+b_{\rm i}^{2}}\right)^{2}}}.$
(18)

The mass of the disk and bulge are taken to be $M_{\rm D}=9\times 10^{10}$ M⊙
and $M_{\rm B}=2\times 10^{10}$ M⊙, respectively. The disk parameters are
($a_{1}$, $b_{1}$) = (7.2 kpc, 0.52 kpc), and the bulge parameters ($a_{2}$,
$b_{2}$) = (0.0 kpc, 0.495 kpc).

###  3.3 Wind solution

####  3.3.1 Method of solution — Boundary conditions and critical points

As the starting conditions of the problem, we set a cosmic-ray pressure
($P_{\rm c}^{0}$), a gas pressure ($P_{\rm g}^{0}$), and a starting value of
the wave pressure ($P_{\rm w}^{0}$) at the base of the wind. We also set an
initial value of the gas mass density, $\rho_{0}$, at the wind base. The wave
pressure can be conveniently prescribed by setting an initial value of
${\left<\delta\textbf{B}\right>}/{\textbf{B}}$ at the base of the wind. The
magnetic field geometry is also prescribed by setting the geometry parameter
$z_{0}$ defined in Eq. (2).

With these conditions, we solve Eqs. (9)–(14) to obtain a supersonic wind
solution. The set of equations has a critical point when $v=c_{*}$ in Eq.
(14), which leads to both the numerator and the denominator of the momentum
conservation equation to go to zero at the sonic point. The continuous
solution across the critical point is obtained using L’Hospital’s rule by
integrating inwards from the critical point to match to the CR and gas
pressure at the disk. Integration is also performed from the critical point
outwards to get the solution up to $O(100~{}\mathrm{kpc})$ distances. Details
on how to solve for a wind solution can be found in Section 4.1 of the paper
by Breitschwerdt et al. (1991).

####  3.3.2 Galactic wind termination shock

The supersonic wind forms a termination shock when the ram pressure of the
flow ($\rho v^{2}$) becomes comparable to the pressure of the ambient medium,
which in our case is the pressure of the intergalactic medium
($P_{\mathrm{IGM}}$). This condition holds as long as the other sources of
pressure like the thermal gas pressure ($P_{\mathrm{g}}$), magnetic pressure
($\frac{B^{2}}{8\pi}$), and cosmic-ray pressure ($P_{\rm c}$) are sufficiently
smaller than the ram pressure. In our model we make sure that the cosmic-ray
pressure (including the pressure of reaccelerated CRs, which is computed
separately as a solution to CR transport equation) at the GWTS never exceeds
the ram pressure. Our reference model forms a GWTS at an altitude of
$z_{\mathrm{s}}\sim 200$ kpc for $P_{\mathrm{IGM}}\sim 10^{-15}$ erg cm-3.
This value of $P_{\mathrm{IGM}}$ is consistent with observations of the
warm–hot IGM (Nicastro et al., 2018) and with $P_{\mathrm{IGM}}$ values
assumed in previous works on cosmic-ray acceleration at termination shocks
(Zirakashvili & Völk, 2006; Bustard et al., 2017; Merten et al., 2018). The
termination shock causes a discontinuity in the speed, temperature, and
density of the fluid, and the strength of the shock jump
($\frac{v_{2}}{v_{1}}$) is obtained using Rankine–Hugoniot jump conditions
(Zel’dovich & Raizer, 1967). The evolution of the flow beyond the termination
shock (downstream) is obtained by assuming the downstream flow to be gas
dominated with adiabatic index $\gamma_{\mathrm{down}}=\frac{5}{3}$.

####  3.3.3 A reference solution

We assume that the wind is launched at a height of $z_{\mathrm{base}}=1$ kpc
from the midplane of the disk. Our results are ultimately not very sensitive
to the launching height. The cosmic ray pressure at the base is taken to be
$P_{\mathrm{c}}^{0}=2\times 10^{-13}$ erg cm-3, consistent with the measured
values of the CR pressure at Earth. The gas number density is taken to be
$n_{0}=2.5\times 10^{-3}$ cm-3, consistent with the observations of Milky
Way’s halo with oxygen lines (Miller & Bregman, 2015). The gas pressure at the
base is set to be $P_{\mathrm{g}}^{0}=3.4\times 10^{-14}$ erg cm-3, which
corresponds to a temperature of $\sim 10^{5}$ K, characteristic of a warm
ionized medium. The geometry parameter $z_{0}$ in Eq. (2) is set to be 20 kpc,
with the expectation that the flux tube opens up from vertical to spherical
geometry at a height of the order of the Galactic radius (Breitschwerdt et
al., 1991; Dorfi & Breitschwerdt, 2012). We also set the magnetic field
$B_{0}$ at the base of the wind which is the field along the $z$ direction.
The Galactic magnetic field has a rich structure consisting of disk and out-
of-plane components (Jansson & Farrar, 2012). Jansson & Farrar (2012) estimate
an out-of-plane field component of $\sim~{}1$ $\mu$G near the Solar circle,
which is what we are interested in. Consistent with this result, the magnetic
field at the base is set to be $B_{0}=2~{}\mu$G. Additionally, at the base, we
set the wave pressure such that $\frac{\delta B_{0}}{B_{0}}=0.1$
(Breitschwerdt et al., 1991), corresponding to a small wave pressure of $\sim
2\times 10^{-15}$ erg cm-3. With these initial parameters, the wind solution
results in a total CR energy flux of $\sim 10^{-5}$ erg cm-2 s-1, consistent
with observations (Murase & Fukugita, 2019).

Figure 3: Top left: speed profiles of the gas, composite sound speed as
defined in Eq. (15), and Alfvén speed. The sonic point occurs at $\sim$ 20
kpc. Top right: baryon number density in the wind profile. Bottom left: mean
magnetic field vs. height from the disk. Bottom right: pressure contribution
from the cosmic rays, waves, and the ram pressure. These figures are shown up
to the location of the GWTS at 200 kpc.

The resulting wind profile as a function of the height from the disk, $z$, is
shown in the top-left panel of Fig. 3. In this model, the gas starts with a
very small speed of $\sim$ 6 km s-1 and increases to $\sim$ 870 km s-1 at
$z\sim 200$ kpc. At the wind base, the magnetic field is, $B_{0}\sim
2~{}\mu{\rm G}$, and the corresponding Alfvén speed is
$v_{\mathrm{A}}=B_{0}/\sqrt{4\pi m_{\rm i}n_{\mathrm{i}}}\sim 100$ km s-1
($m_{\rm i}$, $n_{\mathrm{i}}$ are the mass and number density of ions in the
plasma). This plot shows the velocity profile up to the GWTS. Additionally,
plots for the number density (top-right panel), mean magnetic field (bottom-
left panel) and pressure contributions of the different fluid components
(bottom-right panel) are shown. The number density decreases by four orders of
magnitude between the wind base up to the GWTS. The mean magnetic field
follows the flux tube geometry, decreasing from 2 $\mu{\rm G}$ at the wind
base to $\sim 0.02~{}\mu{\rm G}$ at the GWTS. For the pressure contributions
shown in the bottom-right panel of Fig. 3, we see that the wave pressure
starts out much smaller than the cosmic ray pressure at the base of the wind.
However, it quickly increases and dominates over the CR pressure for most of
the wind halo, where accelerates the gas. The mass outflow rate per unit area
in this model is $\dot{m}\sim 10^{-21}~{}\mathrm{g~{}cm^{-2}s^{-1}}$,
corresponding to an integrated mass outflow rate from the Galaxy of roughly
$\sim$ $0.1$ M⊙ yr-1. This is a non-negligible mass outflow rate, considering
that the star formation of the Galaxy is $\sim O(0.1-1~{}\rm
M_{\odot}yr^{-1})$ (Scoville et al., 2017), indicating that the cosmic-ray-
driven Galactic wind can have a relevant impact on the Galactic ISM and its
dynamical evolution.

Due to the pressure of the intergalactic medium ($\sim 10^{-15}~{}{\rm
g~{}cm}^{-3}$), the GWTS forms at $z_{\mathrm{s}}\sim 200$ kpc. At this
location, the jump in the speed of the gas is approximately
$\frac{u_{2}}{u_{1}}=\frac{\gamma_{\mathrm{gas}}-1}{\gamma_{\mathrm{gas}}+1}+\frac{2}{\gamma_{\mathrm{gas}}+1}\frac{1}{M^{2}}$,
where $M=\frac{v_{\rm gas}}{c_{*}}$.

Beyond the GWTS at 200 kpc, the gas is simply an adiabatically expanding flow
approximately leading to a velocity profile that falls like $v\propto z^{-2}$.

###  3.4 Cosmic-ray transport

The processes involving high-energy particles (most notably the
reacceleration) are taking place over typical timescales that are much shorter
than the dynamical time of the CR-driven Galactic wind. Therefore we assume
that the transport of such particles takes place in stationary condition.

We solve the following stationary, linear transport equation:

$\displaystyle\begin{aligned} \frac{1}{A(z)}\frac{\partial}{\partial
z}\left(A(z)D(z,p)\frac{\partial f}{\partial z}\right)-v(z)\frac{\partial
f}{\partial z}\\\
+\frac{1}{A(z)}\frac{d}{dz}(A(z)v(z))\frac{1}{3}\frac{\partial
f}{\partial\ln{p}}=-Q(z,p),\end{aligned}$ (19)

where $A(z)$ is the area of the wind defined in Eq. (2), $v(z)$ is the wind
speed, $D(z,p)$ is the diffusion coefficient, $Q(z,p)$ is the source term, and
$f=f(z,p)$ is the phase space density of cosmic rays. We note that the
diffusion coefficient used in the transport calculation was not incorporated
into the wind equations (9)–(14). Previous calculations have shown that the
effect of the diffusion coefficient in the wind solutions starts to become
important only for mean diffusion coefficients greater than $\sim 10^{29}$ cm2
s-1 (Dorfi & Breitschwerdt, 2012). The mean diffusion coefficients in our
model varies between $10^{28}-10^{30}$ cm2s-1, depending on the location in
the halo. Therefore, the effect of diffusion in the wind profile could in
general, not be neglected. However, the primary effect of diffusion is to
produce lighter and faster winds with lower mass loss rates (Dorfi &
Breitschwerdt, 2012). Having faster winds will only increase the maximum
energies achievable in the reacceleration process, as will be described later.
Incorporating this effect will therefore only improve the results presented in
this paper. The $v(z)$ term that enters the transport equation contains
contributions from both the gas velocity and the Alfvén speed,
$v(z)=v_{\mathrm{gas}}(z)+v_{\mathrm{A}}(z)$ upstream of the GWTS. We note
that setting $v(z)=v_{\mathrm{gas}}(z)+v_{\mathrm{A}}(z)$ assumes that the CRs
are streaming at the Alfvén speed. This assumption, however, only holds for
CRs up to $\sim 1$ TeV for which streaming instability grows faster than the
dynamical timescale of the problem. For higher energy CRs, the wave growth
will not happen, and they will likely not stream at $v_{\mathrm{A}}(z)$. We
acknowledge this caveat, but we have checked that our results do not
sensitively depend on the value of $v_{\mathrm{A}}(z)$ that enters the
transport equation. Therefore, we keep using
$v=v_{\mathrm{gas}}+v_{\mathrm{A}}$ for the transport calculations of the
entire CR energy range under consideration. For downstream, we assume gas-
dominated transport and simply set $v(z)=v_{\mathrm{gas}}(z)$.

We assume a Gaussian source term for the Galactic CRs:

$Q(z,p)=Q_{0}(p)\frac{\exp\left({-\frac{z^{2}}{2\sigma_{z}^{2}}}\right)}{\sqrt{2\pi\sigma_{z}^{2}}}.$
(20)

Here, $\sigma_{z}\sim 100$ pc is taken as a proxy for the Galactic disk
thickness; however, we verify that the results do not strongly depend on the
value of $\sigma_{z}$ because the scale of the problem, given by the size of
the GWTS, is ${\cal O}(10^{2}~{}\mathrm{kpc})$, which is much larger than the
size of the the central source. The source normalization $Q_{0}(p)$ in Eq.
(20) is given by:

$Q_{0}(p)=\frac{\mathcal{F}_{\mathrm{SN}}\mathcal{R}_{\mathrm{SN}}}{\pi R_{\rm
d}^{2}},$ (21)

where $\mathcal{R}_{\mathrm{SN}}$ is the supernova rate and $R_{\rm d}$ is the
Galactic disk radius. $\mathcal{F}_{\mathrm{SN}}$ is the spectrum of cosmic
rays injected by a single supernova. We assume that supernova remnants (SNRs)
are the sources of Galactic cosmic rays, although the calculations presented
here do not depend crucially on such an assumption. The spectrum of an SNR can
be written as:

$\mathcal{F}_{\rm SN}(p)=\frac{\xi_{\mathrm{CR}}E_{\rm
SN}}{I(\alpha)}\left(\frac{p}{m}\right)^{-\alpha}\exp{\left(-\frac{p}{p_{\rm
c,SNR}}\right)},$ (22)

where $\xi_{\mathrm{CR}}$ is the cosmic-ray injection efficiency by the SNR,
$E_{\rm SN}$ is the kinetic energy released in a supernova explosion, and
$p_{\rm c,SNR}$ is the momentum cutoff of Galactic accelerators.

Finally, $I(\alpha)$ is a normalization factor chosen such that:

$\int_{p_{0}}^{\infty}\mathcal{F}_{\rm
SN}(p)T(p)d^{3}p=\xi_{\mathrm{CR}}E_{\rm SN},$ (23)

where $p_{0}$ is the minimum injection momentum for protons and $T(p)$ is the
kinetic energy of a particle of momentum $p$. For our reference model, we
choose $E_{\rm SN}=10^{51}$ erg, $\xi_{\mathrm{CR}}=0.1$, $p_{0}$ = 1 GeV,
$\mathcal{R_{\mathrm{SN}}}=2$ SN per century. We also take $\alpha=4.4$ and
$p_{\rm c,SNR}=3\times 10^{6}$ GeV for protons, consistent with standard
assumptions in the literature (Thoudam et al., 2016; Recchia et al., 2016,
2017). For heavier elements with atomic number $Z$, the maximum momentum
injected by Galactic sources is $Z$ times the cutoff momentum for protons.

As discussed in Sec. 3.3.3, the base of the wind is set at
$z_{\mathrm{base}}=$ 1 kpc. Below this altitude, the advection term in the
transport equation is taken to be zero, and the cosmic ray propagation is
assumed to be purely diffusive inside this region. Inside the wind base, in
agreement with measurements of the B/C ratio (Aguilar et al., 2016), a
Kolmogorov-type diffusion coefficient is assumed:

$D(p)=2\times
10^{28}~{}\left(\frac{p}{Z}\right)^{1/3}~{}\mathrm{cm^{2}~{}s^{-1}},~{}~{}~{}~{}~{}~{}z<z_{\mathrm{base}}.$
(24)

For $z_{\mathrm{base}}<z<z_{\mathrm{s}}$, i.e, in the upstream wind region up
to the GWTS, we study two different prescriptions for the upstream diffusion
coefficient. We briefly outline these prescriptions below.

####  3.4.1 Diffusion model 1: Quasi-linear theory

Quasi-linear theory (QLT) is a common approach for modeling diffusion
coefficients of magnetised particles in a turbulent medium (Jokipii, 1966;
Shalchi et al., 2004, 2009). Within QLT, the particle motion is assumed to be
a superposition of the gyromotion of the particle and stochastic motion along
the magnetic field lines. This approximation is valid only for weak
turbulence, $\delta B<B$. The scattering rate within the QLT formulation is
$\nu=2\pi^{2}\omega_{B}\frac{k_{\rm res}\epsilon(k_{\rm res})}{B^{2}}$
(Kulsrud & Pearce, 1969; Zweibel, 2013). Here, $\omega_{B}$ denotes the
synchotron frequency proportional to the resonant wavenumber $k_{\rm res}$. At
wavenumber $k_{\rm res}$, the wave energy contained is $k_{\rm
res}\epsilon_{\rm res}$. We assume isotropic turbulence with a Kolmogorov
spectrum and assume that energy is injected at an outer scale
$l_{\mathrm{max}}$ and dissipates to smaller scales after a cascade of energy
from large to small scales without energy loss. The coherence length of the
system, $l_{\mathrm{c}}$, is then given by,
$l_{\mathrm{c}}=l_{\mathrm{max}}/5$ for a Kolmogorov spectrum (Harari et al.,
2002).

As a plausible ansatz, we suppose that an isotropic magnetic turbulence
spectrum is maintained by essentially external processes (other than CRs
generating their own turbulence) that determine the diffusion coefficient,
which in turn controls the particle transport. In particular, we assume that
the wind itself can generate a turbulence at the level $\delta B$ on top of
the mean magnetic field $B$. We define a parameter, $\eta$, such that

$\eta\equiv\frac{\delta B^{2}}{4\pi\,\rho v^{2}}\,,$ (25)

implying $\eta$ is the fraction of energy used in generating turbulence out of
the ram pressure of the wind. The mean magnetic field is obtained from the
wind solution described in Sec. 3.3.3. This turbulence generates a diffusion
coefficient parallel to the mean magnetic field given by Reichherzer et al.
(2020):

$D(p,z)=D_{\mathrm{Bohm}}(p_{0})\left(\frac{p}{p_{0}}\right)^{\frac{1}{3}}\left(\frac{B}{\delta
B}\right)^{2}+D_{\mathrm{Bohm}}(p_{1})\left(\frac{p}{p_{1}}\right)^{2}\left(\frac{B}{\delta
B}\right)^{2},$ (26)

where $r_{L}(p)$ is the Larmor radius computed with respect to the mean
magnetic field. $p_{0}$ and $p_{1}$ are such that $2\pi
r_{L}(p_{0})=l_{\mathrm{c}}$ and $2\pi r_{L}(p_{1})=\frac{3}{2}l_{\mathrm{c}}$
(Globus et al., 2008). The Bohm diffusion coefficient is
$D_{\mathrm{Bohm}}=\frac{1}{3}r_{L}(p)c$. For momenta, $p>p^{*}$, where
$r_{L}(p^{*})=l_{\mathrm{c}}$, the diffusion coefficient changes its energy
dependence to $p^{2}$ due to lack of resonant perturbations (Globus et al.,
2008; Snodin et al., 2016; Subedi et al., 2017; Dundovic et al., 2020). We
prescribe $\frac{\delta B}{B}\sim 0.1$, a condition that can be easily
satisfied by converting $\eta<$ 5% of the wind ram pressure into generating
magnetic turbulence. Note that we keep $\frac{\delta B}{B}\sim 0.1$ fixed
throughout the halo in this case. We also prescribe a coherence length,
$l_{\mathrm{c}}\propto z$, i.e, the coherence length of the turbulent field
scales with height, as a reasonable guess. The prescription is such that
$l_{\mathrm{c}}$ increases from $\sim$ 5 pc to $\sim$ 1 kpc from the wind base
to the GWTS. We note that our results are not very sensitive to the details of
the prescription of the coherence length in the halo. With this prescription,
$D(p,z)$ continuously increases from the wind base up to the GWTS because the
coherence length increasing with $z$ implies an increasing mean free path of
the CRs with height, causing an increasing diffusion coefficient. We will
explain later that this model fails to generate any efficient reacceleration
of the galactic CRs at the GWTS because the diffusion coefficient is too high
close to the GWTS, and the maximum energies achieved in the acceleration
process goes as inverse power of the diffusion coefficient ahead of the shock.
In the following, we will see that in order to achieve more efficient
acceleration, we need to assume that a larger fraction of the ram pressure is
converted into magnetic turbulence close to the GWTS.

####  3.4.2 Diffusion model 2: QLT halo + fully turbulent field near GWTS

In this model, we prescribe a variant of the QLT model presented above that
allows CRs to reach energies up to 10–40 PV rigidities in the reacceleration
process. This is a phenomenological prescription where we allow the
possibility of an increased particle scattering close to the shock, leading to
a lower diffusion coefficient near the GWTS. As we will describe later, this
lower value of the diffusion coefficient close to the GWTS is crucial for an
efficient reacceleration. In this model, we assume the same QLT prescription
as outlined in Sec. 3.4.1 given by Eq. (26) up to $z=100$ kpc with
$\frac{\delta B}{B}=0.1$ and $l_{\mathrm{c}}\propto z$. Between 100 kpc and up
to the GWTS located at 200 kpc, we assume an enhanced turbulence induced by
converting $\eta\sim 20\%$ of the ram pressure of the wind into generating
turbulence. This criterion generates strong turbulence close to the GWTS such
that $\delta B>B$ is satisfied close to the GWTS between 100–200 kpc. We also
set a constant coherence length $l_{\mathrm{c}}=1$ kpc in this region of lower
diffusion coefficient in the halo.

In this strong scattering approximation, the diffusion coefficient near the
GWTS is given by Globus et al. (2008):

$D(p,z)=D_{\mathrm{Bohm}}(p_{0})\left(\frac{p}{p_{0}}\right)^{1/3}+D_{\mathrm{Bohm}}(p_{1})\left(\frac{p}{p_{1}}\right)^{2}.$
(27)

Note that the only difference between Eqs. (26) and (27) is the absence of the
$\left(\frac{B}{\delta B}\right)^{2}$ in the latter. We stress that the only
difference between this model and the previous one is the presence of a strong
scattering region near the GWTS, which as we will see later helps in efficient
reacceleration of the CRs. This strong scattering region is a phenomenological
prescription that we use to demonstrate the conditions under which particles
can be accelerated to up to $\sim$ 40 PV, which is about an order of magnitude
larger than the maximum achievable energies by Galactic CR sources. In
principle, the particle scattering close to the shock can be enhanced by a
number of mechanisms, such as the non-resonant streaming instability (Bell,
2004) and turbulent amplification of magnetic fields by CR-driven pressure
gradients (Beresnyak et al., 2009; Drury & Downes, 2012). These processes
might even decrease the diffusion coeffficient to the Bohm limit at the GWTS,
and such a strong turbulence limit has been assumed in various earlier studies
of CR acceleration at the GWTS (Zirakashvili & Völk, 2006; Dorfi &
Breitschwerdt, 2012; Dorfi et al., 2019). These processes, however, need
sophisticated numerical simulations to accurately assess whether the diffusion
coefficient can be decreased to the Bohm limit close to the GWTS. Our
prescription, instead provides a simple mechanism to generate a strongly
turbulent region close to the shock, which, as will find later, is important
in getting efficient CR reacceleration.

Finally, for the transport equation in the downstream region, we assume as a
benchmark scenario that the downstream diffusion coefficient is an order of
magnitude smaller than the diffusion coefficient right upstream of the GWTS.
This assumption is justified because one can expect that the GWTS will create
strong magnetohydrodynamic (MHD) turbulence downstream towards intergalactic
space (Völk & Zirakashvili, 2004; Donnert et al., 2018), and the diffusion
coefficient there should be small. The resulting diffusion coefficient is
small enough to make the diffusion term negligible compared to the advection
term in the downstream region. We note, however, that the maximum energies
reached in the DSA mechanism decreases with increasing downstream diffusion
coefficient because it increases escape of CRs from the downstream region. For
completeness, we will also consider the effect of having higher downstream
diffusion in Sec. 4.3.

####  3.4.3 Maximum energy achieved in reacceleration

Beyond a certain maximum momentum, $p_{\mathrm{max}}$, the distribution
function at the shock, $f_{s}(p)$, is exponentially suppressed (see Eq. (A23))
because the accelerator has a finite size, so it cannot accelerate particles
to infinitely large momenta. The maximum momentum achieved in reacceleration,
$p_{\mathrm{max}}$, is discussed below.

As seen from Eq. (A23), the exponential drop-off in the distribution function
at the shock, $f_{s}(p)$, is controlled by the two functions $\Gamma_{1}(p)$
and $\Gamma_{2}(p)$. The function $\Gamma_{1}$ represents the effects of the
geometry of the system and the adiabatic energy losses upstream. The
suppression term, e${}^{-\Gamma_{1}}$, is appreciably different from unity for
momenta close to $p_{\mathrm{max}}$, where the upstream diffusion length,
$D_{1}/v$, becomes comparable to the size of the system, $z_{\mathrm{s}}$
(Berezhko & Völk, 1997). For a spatially constant diffusion coefficient
upstream, $D_{1}(p)=\kappa_{1}~{}p^{2-\delta}$ cm2 s-1, and a uniform wind
velocity, $v$, the maximum energy achievable in the GWTS is (Morlino et al.,
2021; Peretti et al., 2022):

$p_{\mathrm{max}}\sim\left(\frac{v~{}z_{\mathrm{s}}}{\kappa_{1}}\right)^{1/(2-\delta)}.$
(28)

Note that Eq. (28) was obtained with the assumption of a spatially constant
upstream diffusion coefficient, but in our reference model, the diffusion
coefficient varies with $z$, and, therefore, $p_{\mathrm{max}}$ in our model
doesn’t necessarily agree with the equation above. However, some basic
features of the solution, like the dependencies of $p_{\mathrm{max}}$ on the
wind speed, shock radius, and diffusion coefficient, are the same as in Eq.
(28). Specifically, we find that $p_{\mathrm{max}}$ decreases with lower wind
speeds, smaller shock sizes, and higher upstream diffusion coefficients. Going
back to our prescription of the two upstream diffusion models in Secs. 3.4.1
and 3.4.2, we find that model 2 provides efficient CR acceleration because, as
we describe here, $p_{\mathrm{max}}$ increases with decreasing $D_{1}$.

The distribution function also depends on $\Gamma_{2}$, which is a result of
the transport in the downstream region. $\Gamma_{2}$ becomes important when
the diffusion length of the particles in the downstream region, $D_{2}/v_{2}$,
becomes comparable to the size of the shocked wind region,
$z_{\mathrm{b}}-z_{\mathrm{s}}$, where $z_{\mathrm{b}}$ is the location of the
outer boundary of the bubble. Typically, for the model scenarios considered
here, $\Gamma_{2}(p)\ll\Gamma_{1}(p)$, and $p_{\mathrm{max}}$ is effectively
set via $\Gamma_{1}(p_{\mathrm{max}})\sim 1$, which is equivalent to the
condition shown in Eq. (28). Therefore, the maximum energy achievable in the
GWTS in our model is set by the upstream plasma velocity, diffusion
coefficient, and the shock size. For more details on the interplay of
$\Gamma_{1}$ and $\Gamma_{2}$ for setting the maximum energy of accelerated
cosmic rays in the wind termination shock, see Morlino et al. (2021).

##  4 Results

###  4.1 Comparison of different diffusion prescriptions

Parameter | Value
---|---
$P_{\mathrm{CR}}^{0}$ | 2 $\times 10^{-13}$ erg cm-3
$P_{\mathrm{g}}^{0}$ |  $3.4\times 10^{-14}$ erg cm-3
$n_{0}$ | 2.5 $\times 10^{-3}$ cm-3
$B_{0}$ | 2 $\mu$G
$z_{0}$ | 20 kpc
$z_{\mathrm{s}}$ | 200 kpc
$v_{\mathrm{t}}$ | 870 km s-1
Table 1: Parameter values for the reference wind model. Values for the base
cosmic-ray pressure ($P_{\mathrm{CR}}^{0}$), gas pressure
($P_{\mathrm{g}}^{0}$), number density ($n_{0}$), magnetic field ($B_{0}$),
geometry parameter ($z_{0}$), shock location ($z_{\mathrm{s}}$), and terminal
gas speed ($v_{\mathrm{t}}=v_{\mathrm{gas}}$) are are all described in Sec.
3.3.3.

First, we summarize the key parameters used for the reference wind model in
Table 1. We also summarize the characteristics of the two upstream diffusion
models as follows:

  * •

Model 1 (QLT) — assumes isotropic turbulence at a length scale
$l_{\mathrm{c}}$. Prescription of a quasi-linear description with
$\frac{\delta B}{B}=0.1$ and $l_{\mathrm{c}}\propto z$ following Eq. (26). The
diffusion coefficient increases starting from the wind base up to the GWTS
(Sec. 3.4.1).

  * •

Model 2 (QLT + strongly turbulent GWTS) — we use QLT prescription following
Eq. (26) up to 100 kpc. Beyond 100 kpc till the shock location at 200 kpc, we
allow for strong scattering and a smaller diffusion coefficient, with
$l_{\mathrm{c}}=1$ kpc and $\delta B>B$ (Sec. 3.4.2).

Figure 4: Spectrum of CR protons at the disk for the different diffusion
prescriptions. Model 1 gives negligible contribution from the GWTS component,
while model 2 is capable of accelerating CRs up to 40 PeV, an order of
magnitude energy gain compared to the CRs produced within the disk.

Fig. 4 shows the proton spectrum at the disk for the two upstream diffusion
models. The spectrum at the disk is expected to be composed of two different
components. One is the Galactic CR contribution with a maximum energy cutoff
of 3 PeV, which marks the knee of the CR spectrum at the disk. Another
contribution is expected from the backstreaming CRs that have been
reaccelerated at the GWTS. The transition between the two components is
expected to lead to a ”bump” feature when the reacceleration is efficient.

For model 1, the spectrum at the disk consists only of the Galactic component
with negligible contribution from the GWTS reacceleration. This is because,
for model 1, the diffusion coefficient near the GWTS, which sets the maximum
energy (Sec. 3.4.3), is $D_{1,\mathrm{GWTS}}^{\mathrm{model~{}1}}\sim
10^{31}\times p_{\mathrm{GeV}}^{1/3}$ cm2 s-1. This diffusion coefficient,
combined with a speed of $\sim 870$ km/s and $z_{s}\sim$ 200 kpc gives a
maximum momentum, $\sim 10^{2}$ GV/c (Eq. 28), which is much smaller than the
$\sim$ 3 PeV energy reached in the Galactic CR sources of protons.
Additionally, these CRs have a harder time propagating back to the disk owing
to their lower rigidities. The result is a negligible contribution of the
reaccelerated particles to the disk spectrum for model 1.

On the other hand, model 2 can reaccelerate the Galactic protons to up to
$\sim$ 40 PeV energies, which is an order of magnitude higher than the knee,
as is shown in the figure. The lower diffusion coefficient of model 2,
$D_{1,\mathrm{GWTS}}^{\mathrm{model~{}2}}\sim(3\times 10^{28}\times
p_{\mathrm{GeV}}^{1/3}+10^{17}p_{\rm GeV}^{2})$ cm2 s-1, allows for a higher
maximum energy. In this case, $D_{1}/v$ becomes equal to the shock size
$z_{s}$ for an energy $\sim O(10~{}\rm PeV)$, which sets the maximum energy of
the system. In this paragraph and the rest of the paper, “maximum energy” in
reference to a plot (such as Fig. 4) refers to the energy at which the
$E^{3}\times$ intensity of the reaccelerated component is roughly smaller by a
factor of $1/e$ than the peak of the bump feature. The maximum rigidity is
that maximum energy divided by charge, $Ze$.

The main difference between model 1 and model 2 is the orders of magnitude
lower diffusion coefficient for the latter close to the GWTS. This leads to
much greater maximum energies achieved for model 2 in accordance with Eq.
(28). We note that the values for the diffusion coefficients in the halo at
$O(100~{}\mathrm{kpc})$ distances is completely unconstrained, thereby
rendering comparison between these extremely different assumptions reasonable.
In the disk, the two models have diffusion coefficient that converge to the
value constrained by the B/C ratio. In the following, we choose diffusion
model 2 to demonstrate other properties of our results, such as the elemental
composition of the CRs, observational implications, parameter space
investigation etc.

###  4.2 Elemental spectrum

Figure 5: Calculated cosmic-ray spectra at Earth for different species
(proton, helium, carbon, iron) compared with data. The reference model (model
2) results are plotted as a thick blue line. The two different components
(Galactic sources and backstreaming) are shown separately only for the proton
case because their behavior is qualitatively identical for the other chemical
elements. The reaccelerated component forms a bump feature due to the
transition from Galactic cosmic rays to the GWTS component. The reaccelerated
component of the proton spectrum at Earth is enough to explain the KASCADE
(2005) data. Additionally, these protons can contribute to the measured flux
at KASCADE-GRANDE and IceCube/IceTop at $\gtrsim$ 10% between 10-40 PeV, with
a maximum contribution of $\sim$ 50% at the peak of the second bump (marked in
green) at $\sim 20$ PeV. The behavior of the spectrum at Earth for heavier
elements is qualitatively similar to protons except the maximum energies
achieved at the Galactic sources and the GWTS are shifted by a factor of $Z$.

In Fig. 5, we plot the spectra, $E^{3}\times\rm intensity$ as a function of
energy for some individual elements and compare them to observed elemental
spectra at the disk from different experiments. The results for proton,
helium, carbon, and iron are shown. The Galactic source injection for the
elements is taken to be: $f_{Z}(p)\propto
p^{-\alpha}\exp{\left(-\frac{p}{Z\times p_{\rm c,SNR}}\right)}$, with
$\alpha=-4.4$ and $p_{\rm c,SNR}=$ 3 PeV. The chemical composition of the
elements is chosen in agreement with the composition inferred from detailed
studies of Galactic cosmic-ray transport problem (see Table 1 of the paper by
Thoudam et al., 2016). The resulting composition consists of $\sim$ 90%
protons, 9% helium, and 1% of remaining heavier nuclei. Note that we do not
take into account the additional complications of Galactic transport like
transitions in the diffusion coefficient (Génolini et al., 2017), particle
losses due to inelastic interactions with interstellar matter, or stochastic
reacceleration of cosmic rays in the interstellar medium which may impact the
observed spectrum at Earth (Thoudam et al., 2016).

The flux produced by the Galactic sources and those contributed by the
reaccelerated backstreaming particles are shown separately for protons in Fig.
5. For protons, we find that the reaccelerated component can contribute
between $\sim 10-50\%$ of the measured proton spectrum at IceCube (Aartsen et
al., 2019) and KASCADE-GRANDE experiments (Arteaga-Velázquez et al., 2018) for
$E\sim 10-40$ PeV, while it fully explains the flux measured by KASCADE
(Antoni et al., 2005). This shows the viability of this mechanism to
substantially contribute to the observed spectrum. The reaccelerated component
is also in agreement with the helium measurements at IceCube (Aartsen et al.,
2019) at $E\sim 30-60$ PeV, although these measurements have large
uncertainties as indicated by the large error bars. Carbon nuclei are
reaccelerated up to $\sim 100-200$ PeV, but their overall contribution to the
carbon spectrum at the disk is largely unconstrained due to the lack of data
in this energy range. Iron nuclei are also reaccelerated to up to $\sim 700$
PeV, and contributes at a level of $\sim$ 10% to the spectra measured by
KASCADE-GRANDE at the disk.

Figure 6: Fluxes of cosmic rays escaping out into the intergalactic medium for
different elements. The plot shows that reacceleration at the GWTS can be a
source of seeding the intergalactic medium with intermediate-energy cosmic
rays, reaching up to $O(1~{}\mathrm{EeV})$ for iron.

In Fig. 6, we plot the escaping flux into the IGM, $E^{3}\times
D\left.\frac{\partial f}{\partial
z}\right|_{z_{\mathrm{b}}}=\frac{E^{3}\times\phi_{esc}(E)}{A(z_{\mathrm{b}})}$
(Eq. (A21)) for different elements. The plot shows that protons and helium can
seed the IGM with particles having energies up to $\sim$
$O(100~{}\mathrm{PeV})$. Heavier nuclei, such as iron, can escape with
energies as high as $\sim$ $O(1~{}\mathrm{EeV})$. Therefore, particles
reaccelerated at the GWTS can play an important role in releasing
intermediate-energy particles into the intergalactic medium. Active star-
forming galaxies, with faster winds will reaccelerate CRs up to even higher
energies of more than $\sim 100$ PV (see e.g., the paper by Peretti et al.,
2022), and can dominantly contribute to seeding the IGM with intermediate-
energy CRs.

###  4.3 Parameter space investigation

In this section we perform a parameter space exploration in order to
understand how the contribution of the reaccelerated component at the disk
changes with changing the main physical parameters of the problem. We focus in
particular on 1) the wind speed and, 2) the diffusion coefficient downstream
of the GWTS. These parameters have no observational constraints and, as we
will see shortly, can sensitively impact the particle (re)acceleration at the
GWTS. In this section, all the parameter space results presented are only for
cosmic ray protons. The conclusions are valid for other heavier elements as
well, albeit with maximum momentum of $Z\times p_{\rm max}^{\rm proton}$.

Figure 7: Proton intensity at the disk as a function of energy as the
diffusion coefficient downstream of the GWTS is increased. The reference model
assumes the downstream diffusion ($D_{2}$) is an order of magnitude smaller
than the diffusion coefficient right upstream ($D_{1}$) of the shock. As the
downstream diffusion coefficient is increased, the maximum energies achieved
at the GWTS decreases, and the overall flux of reaccelerated particles back to
the disk also decreases.

In Fig. 7, we show how the proton flux at the disk changes by parametrically
changing the downstream diffusion coefficient $D_{2}$ from the reference
assumption. In the reference case (blue curve), the downstream diffusion
coefficient is assumed to be an order of magnitude smaller than the upstream
diffusion coefficient ($D_{2}=D_{1}/10$). We show a few more curves as the
downstream diffusion coefficient is increased. The red curve describes a very
pessimistic scenario in which the downstream diffusion coefficient is equal to
the upstream one. We find that increasing the downstream diffusion coefficient
systematically suppresses the overall flux of backstreaming particles. For a
given $D_{1}$, as the downstream diffusion coefficient ($D_{2}$) increases,
the maximum energies achieved at the GWTS decreases (this comes from the
exponential suppression factor $\Gamma_{2}$ in Eq. A23). Another effect of
increasing the downstream diffusion is an increase in the escaping flux of
particles, as defined in Eq. (A21), thereby suppressing the flux of the
particles backstreaming to the disk. If, instead of treating the downstream
diffusion coefficient parametrically (as we have done so far), one assumes
that the decrease in the downstream diffusion coefficient is due to a simple
magnetic field compression at the GWTS, then the compression factor is $\sim
4$, i.e., $B_{2}\sim 4B_{1}$ (typical for strong shocks), where $B_{1,2}$ are
the magnetic fields right ahead and behind the shock, respectively. For this
compression factor111The ratio $D_{2}/D_{1}$ depends on whether the cosmic
rays are in the diffusive regime or the weak-scattering regime in Eq 27. In
model 2 the upstream magnetic field at the shock is about 70 nG, and 10 PV
cosmic rays are in the middle of the two diffusion regimes, where $D$ is
roughly $\propto B^{-1}$., $D_{2}/D_{1}\sim B_{1}/B_{2}$, and the maximum
proton energy achieved is similar to our reference model, i.e, $\sim 30$ PV.
We additionally explored scenarios of extreme suppression of the downstream
diffusion coefficient, namely $D_{2}\lesssim D_{1}/100$. In these
configurations we noticed that the spectrum of backstreaming particles does
not change substantially from our reference scenario since the transport in
the downstream region is practically advection dominated.

Figure 8: Proton intensity at the disk as a function of energy as the terminal
speed ($v_{\mathrm{t}}=v_{\mathrm{gas}}$) is varied. The reference model has
$v_{\mathrm{t}}\sim 870$ km s-1. Decreasing this speed to 200 km s-1 makes the
contribution of the reaccelerated particles to be negligible at the disk.
Increasing the speed to 1000 km s-1 increases the maximum energy achieved in
the reacceleration mechanism.

Fig. 8 illustrates the change in the reaccelerated component at the disk as
the terminal speed (defined as the gas velocity at the shock) is varied by
changing the geometry parameter $z_{0}$. The reference model has a terminal
wind speed of $\sim$ 870 km s-1, as obtained in Sec. 3.3.3 (also see Table 1).
Decreasing the speed to 200 km s-1 decreases the maximum proton energy
achieved at the GWTS to $\sim 5$ PeV, in accordance with Eq. (28). By
contrast, increasing the wind speed to $1000$ km s-1 increases the maximum
energy and shifts the GWTS bump to higher energies. A precise measurement of
the spectrum above the knee can therefore constrain the wind speed in our
model. We have also checked the dependence of the overall flux at the disk on
other parameters, such as the shock radius and slope of the diffusion
coefficients. The effect of changing the shock radius can be qualitatively
understood as follows. As the radius of the shock decreases, the size of the
system decreases, and the maximum energy achieved at the GWTS also decreases
in accordance with Eq. (28). The results remain similar if, instead of the
Kolmogorov term (slope $\frac{1}{3}$) in Eq. (26), we assume Kraichnan
diffusion with a slope of $\frac{1}{2}$. Kraichnan diffusion also gives
acceleration up to $\sim 40$ PeV for protons (equivalently, 40 PV rigidity for
all elements) and can similarly contribute to the observed spectrum at the
disk.

##  5 Summary and Conclusions

We propose that $\sim$ GeV Galactic cosmic rays drive a wind out into the halo
which then passes through a Galactic wind termination shock (GWTS) at a radius
$\sim 200\,{\rm kpc}$, where the wind ram pressure matches the pressure of the
circumgalactic medium. We have applied the theory of diffusive shock
acceleration at this shock and considered the possibility that some of these
cosmic rays are reaccelerated to energies up to $\sim 40$ PeV and diffuse
upstream back to the disk, thereby contributing to the observed spectrum in
the shin region.

In carrying out this investigation, we solved the wind hydrodynamic equations
to find a cosmic-ray-driven wind solution and then used a semi-analytic
iteration technique developed for transport calculations by Morlino et al.
(2021) and Peretti et al. (2022) for the case of spherical termination shocks
in star clusters and starburst galaxies, respectively, to solve for cosmic ray
transport in a Galactic wind. Our procedure takes into account cosmic-ray
injection at the disk, diffusion, advection, adiabatic energy losses, and
reacceleration at the GWTS. For the reacceleration, we considered two models
of diffusive CR transport upstream of the GWTS: “model 1” where $\frac{\delta
B}{B}\sim 0.1$ throughout the halo (weak turbulence introduced at an outer
scale which increases with $z$, as described in Sec. 3.4.1), and “model 2”
where $\frac{\delta B}{B}\sim 1$ close to the GWTS (fully developed magnetic
turbulence at an outer scale $l_{\rm max}\sim 5~{}\rm kpc$, corresponding to a
coherence length, $l_{\rm c}\sim 1~{}\rm kpc$, as described in Sec. 3.4.2 ).

Some important points about the model and our main results are:

  * •

The existence of a Galactic wind is strongly suggested by the requirement that
local cosmic rays escape in less than a hundred million years (Lipari, 2014).
The cosmic-ray energy flux from the local Galactic disk is adequate to drive
the wind and most of the cosmic rays out of the Galactic gravitational
potential.

  * •

The presence of a wind carrying off mass, angular momentum, and energy from
the Galactic disk should be an important feature of models of the interstellar
medium because these winds can have mass loss rates of $O(0.1~{}\rm
M_{\odot}yr^{-1})$, comparable to the star formation rate of the Galaxy
(Scoville et al., 2017).

  * •

We tested our results with two different prescriptions for the upstream
diffusion coefficient. We found that if the diffusion coefficient ahead of the
GWTS is too high, we get negligible reacceleration at the shock. The diffusion
coefficient ahead of the GWTS needs to be sufficiently small to give enough
reaccelerationup to $O(10~{}\rm PV)$ range.

  * •

In the latter case, cosmic rays can be efficiently reaccelerated at the
Galactic wind termination shock up to rigidities ranging from 10 PV to 40 PV.
This is the first time, to the best of our knowledge, that the maximum energy
is derived self-consistently based on a transport model that derived from the
physical parameters of the wind in the context of reacceleration at the GWTS.

  * •

We note that even in the limiting case when the turbulent pressure near the
shock is comparable with the ram pressure, we do not reach rigidities higher
than $\sim$ 70 PV in our model. In particular, they cannot account for the
observed cosmic rays all the way up to the ankle in the spectrum at energy
$\sim 4\,{\rm EeV}$. If there are additional sources of turbulence ahead of
the GWTS, such as turbulence generated by CRs themselves via resonant or non-
resonant streaming instabilities, and if these can cause the diffusion
coefficient to decrease close to the Bohm limit ahead of the shock, one could
possibly get reacceleration up to higher rigidities.

  * •

Even though the wind strongly suppresses the absolute number of reaccelerated
particles that can diffuse back to the Galactic disk, we find that, for our
benchmark scenarios, the flux of backstreaming protons can be as high as $\sim
50\%$ of the observed flux in the 10–40 PeV region measured in IceCube/IceTop
and KASCADE-Grande experiments. The backstreaming flux at the disk is marked
by a characteristic ‘bump’ feature formed by the reaccelerated proton
component at $\sim$ 10 PeV that could be detected in future observations.

  * •

The location of the outer boundary of the downstream flow behind the
termination shock at radius $\sim$ 500 kpc is chosen so as to roughly mimic
the distance corresponding to three times the virial radius of the Galaxy from
which cosmic rays can escape freely. We have verified that our results are
relatively robust with respect to changes in the location of the outer
boundary as long as it is located sufficiently far from the Galactic wind
termination shock ($z_{\mathrm{b}}\gtrsim 2\,z_{\mathrm{s}}$).

  * •

A fraction of the reaccelerated cosmic rays can seed the intergalactic medium
with intermediate energy particles. Protons with energies up to $\sim$ 40 PeV
and heavier nuclei like iron with energies $\sim$ $O(1~{}\mathrm{EeV})$ can be
transmitted into the intergalactic medium. These reaccelerated cosmic rays
that escape downstream from termination shocks around active star forming
galaxies can undergo further reacceleration to rigidities $\sim$ 10 EV by even
larger intergalactic shock fronts, which may ultimately account for the
spectrum of ultra-high-energy cosmic rays seen at Earth.

Our investigation opens up several directions for future studies. In order to
fully assess the velocity profile of the Galactic wind and the characteristics
of the termination shock, one needs to perform full MHD calculations of the
wind along with the transport equation taking into account cosmic-ray
pressure, thermal pressure, and magnetic stresses. Moreover, the reaccelerated
cosmic rays can cause feedback effects on the GWTS, which can only be
realistically modeled in the MHD framework. Additionally, the reaccelerated
particles can produce multi-messenger signals during propagation, like gamma
rays and neutrinos, which will be investigated elsewhere.

## Acknowledgements

We thank Chad Bustard, Lukas Merten, Vladimir Zirakashvilli, and Satyendra
Thoudam for answering various questions related to particle acceleration at
termination shocks. We thank Patrick Reichherzer for answering questions on
values diffusion coefficients. We thank Giovanni Morlino for fruitful
discussions. We also thank Sarah Recchia for many insightful comments on the
nature of galactic winds and cosmic-ray transport therein. The research
activity of E.P. was supported by Villum Fonden (project n. 18994) and by the
European Union’s Horizon 2020 research and innovation program under the Marie
Sklodowska-Curie grant agreement No. 847523 ‘INTERACTIONS’. N.G.’s research is
supported by the Simons Foundation, the Chancellor Fellowship at UCSC and the
Vera Rubin Presidential Chair.

## References

  * Aab et al. (2018) Aab, A., Abreu, P., Aglietta, M., et al. 2018, ApJ, 868, 4, doi: 10.3847/1538-4357/aae689
  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, Phys. Rev. D, 100, 082002, doi: 10.1103/PhysRevD.100.082002
  * Abeysekara et al. (2021a) Abeysekara, A. U., et al. 2021a, Nature Astron., 5, 465, doi: 10.1038/s41550-021-01318-y
  * Abeysekara et al. (2021b) —. 2021b, PoS, ICRC2021, 811, doi: 10.22323/1.395.0811
  * Adriani et al. (2019) Adriani, O., Akaike, Y., Asano, K., et al. 2019, Phys. Rev. Lett., 122, 181102, doi: 10.1103/PhysRevLett.122.181102
  * Adriani et al. (2021) —. 2021, Phys. Rev. Lett., 126, 241101, doi: 10.1103/PhysRevLett.126.241101
  * Aguilar et al. (2015) Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Phys. Rev. Lett., 114, 171103, doi: 10.1103/PhysRevLett.114.171103
  * Aguilar et al. (2016) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2016, Phys. Rev. Lett., 117, 231102, doi: 10.1103/PhysRevLett.117.231102
  * Aguilar et al. (2017) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2017, Phys. Rev. Lett., 119, 251101, doi: 10.1103/PhysRevLett.119.251101
  * Aguilar et al. (2020) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2020, Phys. Rev. Lett., 124, 211102, doi: 10.1103/PhysRevLett.124.211102
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, 464, 235, doi: 10.1051/0004-6361:20066381
  * Ahlers (2019) Ahlers, M. 2019, Astrophys. J. Lett., 886, L18, doi: 10.3847/2041-8213/ab552f
  * Albert et al. (2020) Albert, A., Alfaro, R., Alvarez, C., et al. 2020, ApJ, 896, L29, doi: 10.3847/2041-8213/ab96cc
  * Alemanno et al. (2021) Alemanno, F., An, Q., Azzarello, P., et al. 2021, Phys. Rev. Lett., 126, 201102, doi: 10.1103/PhysRevLett.126.201102
  * Amenomori et al. (2019) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Phys. Rev. Lett., 123, 051101, doi: 10.1103/PhysRevLett.123.051101
  * Andeen & Plum (2019) Andeen, K., & Plum, M. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 172, doi: 10.22323/1.358.0172
  * Antoni et al. (2005) Antoni, T., Apel, W. D., Badea, A. F., et al. 2005, Astroparticle Physics, 24, 1, doi: 10.1016/j.astropartphys.2005.04.001
  * Apel et al. (2011) Apel, W. D., Arteaga-Velázquez, J. C., Bekk, K., et al. 2011, Phys. Rev. Lett., 107, 171104, doi: 10.1103/PhysRevLett.107.171104
  * Apel et al. (2013) Apel, W. D., Arteaga-Velàzquez, J. C., Bekk, K., et al. 2013, Phys. Rev. D, 87, 081101, doi: 10.1103/PhysRevD.87.081101
  * Arteaga-Velázquez et al. (2018) Arteaga-Velázquez, C. J., et al. 2018, PoS, ICRC2017, 316, doi: 10.22323/1.301.0316
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550, doi: 10.1111/j.1365-2966.2004.08097.x
  * Beresnyak et al. (2009) Beresnyak, A., Jones, T. W., & Lazarian, A. 2009, ApJ, 707, 1541, doi: 10.1088/0004-637X/707/2/1541
  * Berezhko & Völk (1997) Berezhko, E. G., & Völk, H. J. 1997, Astroparticle Physics, 7, 183, doi: 10.1016/S0927-6505(97)00016-9
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Blasi (2002) Blasi, P. 2002, Astroparticle Physics, 16, 429, doi: 10.1016/S0927-6505(01)00127-X
  * Blasi (2004) —. 2004, Astroparticle Physics, 21, 45, doi: 10.1016/j.astropartphys.2003.10.008
  * Blasi (2013) —. 2013, A&A Rev., 21, 70, doi: 10.1007/s00159-013-0070-7
  * Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79
  * Burke (1968) Burke, J. A. 1968, MNRAS, 140, 241, doi: 10.1093/mnras/140.2.241
  * Bustard et al. (2017) Bustard, C., Zweibel, E. G., & Cotter, C. 2017, ApJ, 835, 72, doi: 10.3847/1538-4357/835/1/72
  * Cao et al. (2021) Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33, doi: 10.1038/s41586-021-03498-z
  * Chevalier & Clegg (1985) Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44, doi: 10.1038/317044a0
  * Donnert et al. (2018) Donnert, J., Vazza, F., Brüggen, M., & ZuHone, J. 2018, Space Sci. Rev., 214, 122, doi: 10.1007/s11214-018-0556-8
  * Dorfi & Breitschwerdt (2012) Dorfi, E. A., & Breitschwerdt, D. 2012, A&A, 540, A77, doi: 10.1051/0004-6361/201118082
  * Dorfi et al. (2019) Dorfi, E. A., Steiner, D., Ragossnig, F., & Breitschwerdt, D. 2019, A&A, 630, A107, doi: 10.1051/0004-6361/201935940
  * Drury & Downes (2012) Drury, L. O., & Downes, T. P. 2012, MNRAS, 427, 2308, doi: 10.1111/j.1365-2966.2012.22106.x
  * Dundovic et al. (2020) Dundovic, A., Pezzi, O., Blasi, P., Evoli, C., & Matthaeus, W. H. 2020, Phys. Rev. D, 102, 103016, doi: 10.1103/PhysRevD.102.103016
  * Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, ApJ, 674, 258, doi: 10.1086/524766
  * Génolini et al. (2017) Génolini, Y., et al. 2017, Phys. Rev. Lett., 119, 241101, doi: 10.1103/PhysRevLett.119.241101
  * Globus et al. (2008) Globus, N., Allard, D., & Parizot, E. 2008, A&A, 479, 97, doi: 10.1051/0004-6361:20078653
  * Globus et al. (2015) —. 2015, Phys. Rev. D, 92, 021302, doi: 10.1103/PhysRevD.92.021302
  * Gottlieb & Globus (2021) Gottlieb, O., & Globus, N. 2021, ApJ, 915, L4, doi: 10.3847/2041-8213/ac05c5
  * Harari et al. (2002) Harari, D., Mollerach, S., Roulet, E., & Sánchez, F. 2002, Journal of High Energy Physics, 2002, 045, doi: 10.1088/1126-6708/2002/03/045
  * HESS Collaboration et al. (2016) HESS Collaboration, Abramowski, A., Aharonian, F., et al. 2016, Nature, 531, 476, doi: 10.1038/nature17147
  * Ipavich (1975) Ipavich, F. M. 1975, ApJ, 196, 107, doi: 10.1086/153397
  * Jansson & Farrar (2012) Jansson, R., & Farrar, G. R. 2012, ApJ, 757, 14, doi: 10.1088/0004-637X/757/1/14
  * Johnson & Axford (1971) Johnson, H. E., & Axford, W. I. 1971, ApJ, 165, 381, doi: 10.1086/150903
  * Jokipii (1966) Jokipii, J. R. 1966, ApJ, 146, 480, doi: 10.1086/148912
  * Jokipii & Morfill (1987) Jokipii, J. R., & Morfill, G. 1987, ApJ, 312, 170, doi: 10.1086/164857
  * Jokipii & Morfill (1985) Jokipii, J. R., & Morfill, G. E. 1985, ApJ, 290, L1, doi: 10.1086/184430
  * Kang et al. (2021) Kang, D., et al. 2021, PoS, ICRC2019, 306, doi: 10.22323/1.358.0306
  * Kulpa-Dybeł et al. (2011) Kulpa-Dybeł, K., Otmianowska-Mazur, K., Kulesza-Żydzik, B., et al. 2011, ApJ, 733, L18, doi: 10.1088/2041-8205/733/2/L18
  * Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445, doi: 10.1086/149981
  * Lerche (1967) Lerche, I. 1967, Apj, 147, 689, doi: 10.1086/149045
  * Lipari (2014) Lipari, P. 2014, arXiv e-prints, arXiv:1407.5223.  https://arxiv.org/abs/1407.5223
  * Mao & Ostriker (2018) Mao, S. A., & Ostriker, E. C. 2018, ApJ, 854, 89, doi: 10.3847/1538-4357/aaa88e
  * Merten et al. (2018) Merten, L., Bustard, C., Zweibel, E. G., & Becker Tjus, J. 2018, Astrophys. J., 859, 63, doi: 10.3847/1538-4357/aabfdd
  * Miller & Bregman (2015) Miller, M. J., & Bregman, J. N. 2015, ApJ, 800, 14, doi: 10.1088/0004-637X/800/1/14
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Monari et al. (2018) Monari, G., Famaey, B., Carrillo, I., et al. 2018, A&A, 616, L9, doi: 10.1051/0004-6361/201833748
  * Morlino et al. (2021) Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096, doi: 10.1093/mnras/stab690
  * Murase & Fukugita (2019) Murase, K., & Fukugita, M. 2019, Phys. Rev. D, 99, 063012, doi: 10.1103/PhysRevD.99.063012
  * Necib & Lin (2022) Necib, L., & Lin, T. 2022, Astrophys. J., 926, 189, doi: 10.3847/1538-4357/ac4244
  * Nicastro et al. (2018) Nicastro, F., Kaastra, J., Krongold, Y., et al. 2018, Nature, 558, 406, doi: 10.1038/s41586-018-0204-1
  * Peretti et al. (2022) Peretti, E., Morlino, G., Blasi, P., & Cristofari, P. 2022, MNRAS, 511, 1336, doi: 10.1093/mnras/stac084
  * Piffl et al. (2014) Piffl, T., Scannapieco, C., Binney, J., et al. 2014, A&A, 562, A91, doi: 10.1051/0004-6361/201322531
  * Recchia et al. (2016) Recchia, S., Blasi, P., & Morlino, G. 2016, MNRAS, 462, 4227, doi: 10.1093/mnras/stw1966
  * Recchia et al. (2017) —. 2017, MNRAS, 470, 865, doi: 10.1093/mnras/stx1214
  * Reichherzer et al. (2020) Reichherzer, P., Becker Tjus, J., Zweibel, E. G., Merten, L., & Pueschel, M. J. 2020, MNRAS, 498, 5051, doi: 10.1093/mnras/staa2533
  * Scoville et al. (2017) Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150, doi: 10.3847/1538-4357/aa61a0
  * Shalchi et al. (2004) Shalchi, A., Bieber, J. W., Matthaeus, W. H., & Qin, G. 2004, ApJ, 616, 617, doi: 10.1086/424839
  * Shalchi et al. (2009) Shalchi, A., Skoda, T., Tautz, R. C., & Schlickeiser, R. 2009, A&A, 507, 589, doi: 10.1051/0004-6361/200912755
  * Snodin et al. (2016) Snodin, A. P., Shukurov, A., Sarson, G. R., Bushby, P. J., & Rodrigues, L. F. S. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3975, doi: 10.1093/mnras/stw217
  * Soida et al. (2011) Soida, M., Krause, M., Dettmar, R. J., & Urbanik, M. 2011, A&A, 531, A127, doi: 10.1051/0004-6361/200810763
  * Subedi et al. (2017) Subedi, P., Sonsrettee, W., Blasi, P., et al. 2017, ApJ, 837, 140, doi: 10.3847/1538-4357/aa603a
  * Thoudam et al. (2016) Thoudam, S., Rachen, J. P., van Vliet, A., et al. 2016, A&A, 595, A33, doi: 10.1051/0004-6361/201628894
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration, Amenomori, M., Bao, Y. W., et al. 2021, Nature Astronomy, 5, 460, doi: 10.1038/s41550-020-01294-9
  * Völk & Zirakashvili (2004) Völk, H. J., & Zirakashvili, V. N. 2004, A&A, 417, 807, doi: 10.1051/0004-6361:20040018
  * Wentzel (1974) Wentzel, D. G. 1974, ARA&A, 12, 71, doi: 10.1146/annurev.aa.12.090174.000443
  * Zel’dovich & Raizer (1967) Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena
  * Zhang et al. (2020) Zhang, Z., Murase, K., & Mészáros, P. 2020, MNRAS, 492, 2250, doi: 10.1093/mnras/staa022
  * Zirakashvili & Völk (2006) Zirakashvili, V. N., & Völk, H. J. 2006, Advances in Space Research, 37, 1923, doi: 10.1016/j.asr.2005.06.013
  * Zweibel (2013) Zweibel, E. G. 2013, Physics of Plasmas, 20, 055501, doi: 10.1063/1.4807033

##  Appendix A Iteration Technique

To solve Eq. (19), we use an iteration technique first introduced in Blasi
(2002, 2004), which was subsequently applied for spherical termination shocks
in Morlino et al. (2021). We recast Eq. (19) into the following form:

$\displaystyle\frac{\partial}{\partial z}\left(A(z)\,D(z,p)\frac{\partial
f}{\partial z}-v(z)\,A(z)f\right)+\frac{\differential}{\differential
z}(A(z)\,v(z))\left(f+\frac{p}{3}\frac{\partial f}{\partial
p}\right)+A(z)\,Q(z,p)=0.$ (A1)

We define a variable $q$ such that

$q=-\left(3+\frac{\partial\ln{f}}{\partial\ln{p}}\right).$ (A2)

The wind velocity profile $v(z)$ is obtained from the solution to the wind
equations described in Section 3.3.3. Note, crucially, that we do not solve
the wind equations and the transport equation self-consistently. The wind
calculation is done separately, and from that the velocity profile is plugged
into the transport equation. It is assumed that the wind $v(z)$ is launched
from a height $z_{0}=1$ kpc. Within this boundary ($z<$ 1 kpc), the wind is
assumed to quickly drop to zero velocity so that the propagation within the
Galactic disk and the inner magnetized halo is assumed to be diffusion
dominated. This assumption is in agreement with standard CR transport modeling
in a thin disk of 100-parsec height surrounded by a kiloparsec-sized
magnetized halo where CRs are diffusing (Blasi, 2013).

The solution of the transport equation is obtained by first solving the
equation in the unshocked region (upstream) and the shocked region
(downstream) and then by joining the two solutions across the shock surface
where the solution is $f_{s}(z,p)$. The wind velocity and the diffusion
coefficient $D(z,p)$ are spatially varying quantities. We now describe the
iterative procedure in detail.

In the upstream region ($z<z_{\mathrm{s}}$), the symmetry of the problem
imposes the inner boundary condition at $z=0$, resulting in a simple zero-flux
requirement:

$\left(\frac{\partial f}{\partial
z}\right)_{z=0}=0~{}~{}~{}~{}~{}(\mathrm{inner~{}B.C}).$ (A3)

Integrating from $z=0$ to some upstream $z$, we obtain,

$\left(A(z)\,D(z,p)\frac{\partial f}{\partial
z}-v(z)\,A(z)f\right)-\int_{0}^{z}\frac{\differential}{\differential
z}(A(z)\,v(z))\frac{f}{3}\,q\,\differential
z+\int_{0}^{z}A(z^{\prime})\,Q(z^{\prime},p)\,\differential z^{\prime}=0.$
(A4)

We define a new quantity $G(z,p)$ such that,

$G(z,p)=\int_{0}^{z}\frac{\differential}{\differential
z}(A(z)\,v(z))\frac{f}{3}q\,\differential z.$ (A5)

Simplifying Equation (A4), we obtain,

$\frac{\partial f}{\partial
z}=\frac{v(z)f(z,p)}{D(z,p)}+\frac{G(z,p)}{A(z)\,D(z,p)}-\frac{\int_{0}^{z}A(z^{\prime})\,Q(z^{\prime},p)\differential
z^{\prime}}{A(z)\,D(z,p)}.$ (A6)

We define two new quantities:

$V_{\rm
eff}(z,p)=v(z)+\frac{G(z,p)}{A(z)f},~{}~{}~{}I(z,p)=\int_{0}^{z}A(z^{\prime})\,Q(z^{\prime},p)\,\differential
z^{\prime}.$ (A7)

Using Eqs. (A6) and (A7), we get

$\frac{\partial f}{\partial z}=\frac{V_{\rm
eff}(z,p)}{D(z,p)}f-\frac{I}{D(z,p)\,A(z,p)}.$ (A8)

Eq. (A8) has the following implicit solution:

$f_{1}(z,p)\exp{\left(-\int_{0}^{z}\frac{V_{\rm
eff}}{D(z^{\prime},p)}\differential
z^{\prime}\right)}=C_{1}-\int_{0}^{z}\frac{I}{D(z^{\prime},p)\,A(z^{\prime})}\exp{\left(-\int_{0}^{z^{\prime}}\frac{V_{\rm
eff}}{D(z^{\prime\prime},p)}\differential
z^{\prime\prime}\right)}\,\differential z^{\prime},$ (A9)

where $C_{1}$ is a constant to be determined. Equation (A9) gives the implicit
solution for the upstream distribution $f(z,p)$. The condition that at the
shock ($z=z_{\mathrm{s}}$) that $f=f_{s}(p)$ fixes $C_{1}$. The final solution
for the upstream distribution after imposing this condition gives the upstream
solution in an implicit form:

$f_{1}(z,p)=f_{s}\exp{\left(-\int_{z}^{z_{\mathrm{s}}}\frac{V_{\rm
eff}(z,p)}{D(z,p)}\differential
z^{\prime}\right)}+\int_{z}^{z_{\mathrm{s}}}\frac{I}{D(z^{\prime},p)\,A(z^{\prime})}\exp{\left(-\int_{z}^{z^{\prime}}\frac{V_{\rm
eff}}{D(z^{\prime\prime},p)}\differential
z^{\prime\prime}\right)}\,\differential z^{\prime}.$ (A10)

In the downstream region, $z_{\mathrm{s}}<z<z_{\mathrm{b}}$, where
$z_{\mathrm{b}}$ is the escape height, the solution is found by integrating
the transport equation from $z_{\mathrm{s}}$ to $z<z_{\mathrm{b}}$:

$\left[A(z)\,D(z,p)\frac{\partial f}{\partial
z}\right]_{z_{s^{+}}}^{z}-v(z)\,A(z)f+v_{2}A_{s}f_{s}+G_{2}(z,p)-G_{2}(z_{\mathrm{s}},p)=0,$
(A11)

where we defined

$v_{2}=v(z_{\mathrm{s}}^{+}),~{}~{}~{}~{}G_{2}(z,p)=\int_{z}^{z_{\mathrm{b}}}\frac{\differential}{\differential
z^{\prime}}(A(z^{\prime})\,v(z^{\prime}))\frac{f}{3}q\,\differential
z^{\prime}.$ (A12)

The downstream solution immediately after the shock becomes:

$\left(D\frac{\partial f}{\partial
z}\right)_{z_{\mathrm{s}}^{+}}=v_{2}f_{s}-\frac{\phi_{\rm
esc}}{A(z)}-\frac{G_{2}(z_{\mathrm{s}}^{+},p)}{A(z)},$ (A13)

where we have defined the quantity $\phi_{\rm esc}$ at the boundary as

$\phi_{\rm esc}=-A(z_{\mathrm{b}})\left[D(z,p)\frac{\partial f}{\partial
z}-v(z)\,f\right]_{z=z_{\mathrm{b}}}.$ (A14)

Using Equations (A11) and (A13), we can write the general downstream solution
in an implicit form as

$f_{2}(z,p)=\int_{z}^{z_{\mathrm{b}}}\left(\frac{G_{2}(z^{\prime},p)+\phi_{\rm
esc}}{A(z^{\prime})\,D(z^{\prime},p)}\right)\times\exp{\left(-\int_{z}^{z^{\prime}}\frac{v(z^{\prime\prime})}{D(z^{\prime\prime},p)}\differential
z^{\prime\prime}\right)}\differential z^{\prime}.$ (A15)

In terms of $f_{\rm s}$, the above equation becomes

$\begin{split}f_{2}(z,p)=\left[f_{\rm
s}(p)-\int_{z_{\mathrm{s}}}^{z}\left(\frac{G_{2}(z^{\prime},p)+\phi_{\rm
esc}}{A(z^{\prime})\,D(z^{\prime},p)}\right)\times\exp{\left(-\int_{z_{\mathrm{s}}}^{z^{\prime}}\frac{v(z^{\prime\prime})}{D(z^{\prime\prime},p)}\differential
z^{\prime\prime}\right)}\differential
z^{\prime}\right]\exp{\left(\int_{z_{\mathrm{s}}}^{z}\frac{v(z^{\prime})}{D(z^{\prime},p)}\differential
z^{\prime}\right)}.\end{split}$ (A16)

Equations (A10) and (A16) give the general solutions for the upstream
($f_{1}$) and downstream ($f_{2}$) distribution function in an implicit form.
The quantity $\phi_{\rm esc}$ has to be obtained from the outer boundary
condition to be described later. The two solution for $f_{1}$ and $f_{2}$
needs to be joined at the shock front.

At the shock front, the plasma speed discontinuously changes $v(z)\simeq
v_{2}+(v_{1}-v_{2})\Theta(z_{\mathrm{s}}-z)$, where $\Theta$ is the Heaviside
theta function and $u_{1},u_{2}$ are the speeds just before and after the
shock front. The derivative of a Heaviside function is the Dirac delta
function, leading to $\left(\frac{\differential}{\differential
z}(A(z)\,v(z))\right)_{z_{\mathrm{s}}}=-A_{\rm
s}(v_{1}-v_{2})\,\delta(z-z_{\mathrm{s}})$, where $A_{\rm s}$ is the area at
the shock front. Integrating the transport equation (Eq. 19) between
$z_{\mathrm{s}}^{-}$ to $z_{\mathrm{s}}^{+}$, we get

$(v_{1}-v_{2})\frac{p}{3}\frac{\partial f_{\rm s}}{\partial
p}=\left(D\frac{\partial f}{\partial
z}\right)_{z_{\mathrm{s}}^{+}}-\left(D\frac{\partial f}{\partial
z}\right)_{z_{\mathrm{s}}^{-}}.$ (A17)

The terms on the RHS of the above equation can be obtained from Equations (A8)
and (A13). Equation (A17) can then be simplified to

$(v_{1}-v_{2})\frac{p}{3}\frac{\partial f_{\rm s}}{\partial
p}=-(v_{1}-v_{2})f_{\rm s}-\left(\frac{\phi_{\rm
esc}}{A(z)}+\frac{G_{2}(z_{\mathrm{s}}^{+},p)}{A(z)}\right)-\frac{G_{1}(z,p)}{A(z)}+\frac{I(z_{\mathrm{s}},p)}{A(z)}.$
(A18)

The idea of the iteration technique is to implicitly solve for Equations
(A10), (A16), and (A18) in an iterative procedure.

In our case, for the downstream solution, the speed goes like $v\propto
z^{-2}$ in the downstream region as obtained from solving the wind equations.
For such large radii $z>z_{\mathrm{s}}$, where $z_{\mathrm{s}}\sim 200$ kpc,
the area also approximately becomes quasi-spherical $A(z)\propto z^{2}$. This
implies that $\frac{\differential}{\differential z}(A(z)\,v(z))\sim 0$ for the
large downstream radii. Therefore, $G_{2}(z,p)=0$, which implies zero
adiabatic losses in the downstream region. In this limit, the downstream
solution becomes analytical. Since the plasma is heated up to high
temperatures by the shock, it is reasonable to assume that the downstream
magnetic field as well as the diffusion coefficient $D_{2}(p)$ are spatially
constant. The latter one has therefore only a $p$ dependence. We define

$\chi_{2}(z,p)\equiv\int_{z_{\mathrm{s}}}^{z}\frac{v(z^{\prime})}{D_{2}(p)}dz^{\prime}=\frac{v_{2}\,z_{\mathrm{s}}}{D_{2}(p)}\left(1-\frac{z_{\mathrm{s}}}{z}\right).$
(A19)

For the far boundary condition, we assume

$f_{b}<<f_{\rm s}~{}~{}~{}~{}~{}(\mathrm{outer~{}B.C}),$ (A20)

where $f_{\rm b}$ is the distribution function at the outer boundary of the
bubble which we take to be $z_{\mathrm{b}}=500$ kpc, which is approximately
the order of three times the virial radius of Milky Way according to Global
Astrometric Interferometer for Astrophysics (GAIA) and RAdial Velocity
Experiment (RAVE) survey estimates (Piffl et al., 2014; Monari et al., 2018).
This free-escape type outer boundary condition is justified because at the
outer boundary of the bubble, the diffusion coefficient in the IGM is expected
to increase by many orders of magnitude than the diffusion coefficient in the
downstream region. We check that the exact location of the outer boundary does
not impact our results as long as the outer boundary is sufficiently far away
from the termination shock. Using the far boundary condition and Eq. (A16),
setting $f_{2}=f_{\rm b}$ and using $f_{\rm b}<<f_{\rm s}$ we get

$\phi_{\rm esc}=\frac{v_{2}\,f_{\rm
s}(p)\,A(z_{\mathrm{s}})}{1-e^{-\chi_{2}(p,z_{\mathrm{b}})}}.$ (A21)

The distribution function at the shock, $f_{\rm s}$ is obtained by inserting
$\phi_{\rm esc}$ from Eq. (A21) to Eq. (A18) to obtain

$\begin{split}p\frac{\partial f_{\rm s}}{\partial
p}=-\frac{3v_{1}}{v_{1}-v_{2}}\left(\left(1-\frac{u_{2}/u_{1}}{1-e^{\chi_{2}(p^{\prime},z_{\mathrm{b}})}}+\frac{G(p)}{v_{1}A(z_{\mathrm{s}})f_{\rm
s}}\right)f_{\rm
s}-\frac{I(z_{\mathrm{s}},p)}{v_{1}A(z_{\mathrm{s}})}\right),\end{split}$
(A22)

where $G(p)=G_{1}(z_{\mathrm{s}},p)$. The solution of Eq. (A22) can be
expressed in an implicit form as

$\begin{split}f_{\rm
s}(p)=\left[\frac{sp^{-s}}{v_{1}\,A(z_{\mathrm{s}})}\int_{p_{0}}^{p}\frac{I(z_{\mathrm{s}},p^{\prime})}{p^{\prime}}p^{{}^{\prime}s}e^{\Gamma_{1}(p^{\prime})}e^{\Gamma_{2}(p^{\prime})}dp^{\prime}\right]\times~{}e^{-\Gamma_{1}(p)}e^{-\Gamma_{2}(p)},\end{split}$
(A23)

where $s=\frac{3v_{1}}{v_{1}-v_{2}}$ and $p_{0}$ is the injection momentum of
the central source. We also define the velocity jump factor across the shock
to be $\sigma=\frac{v_{1}}{v_{2}}$. The exponential terms $\Gamma_{1}(p)$ and
$\Gamma_{2}(p)$ are

$\begin{split}\Gamma_{1}(p)=s\int_{p_{0}}^{p}\frac{G(p^{\prime})}{v_{1}A(z_{\mathrm{s}})f_{s}(p^{\prime})}\frac{dp^{\prime}}{p^{\prime}},\\\
\Gamma_{2}(p)=\frac{s}{\sigma}\int_{p_{0}}^{p}\frac{1}{e^{\chi_{2}(p^{\prime},z_{\mathrm{b}})}-1}\frac{dp^{\prime}}{p^{\prime}}.\end{split}$
(A24)

Eqs. (A10), (A16), (A23), and (A24) form a closed set of equations that can be
solved using the iteration technique. The iteration algorithm is described as
follows. One starts with a guess function for $f_{1}$ and $f_{s}$. A
convenient guess function is the solution for a planar shock. Then in each
subsequent $k^{th}$ step, $q^{k}(z,p)$, $G_{1}^{k}(z,p)$, and
$\Gamma_{1}^{k}(p)$ are computed using Eqs. (A2), (A5), and (A24) in sequence.
$f^{k}_{1}(z,p)$ and $f^{k}_{s}(p)$ are then updated at each step using the
values of $\Gamma_{1}^{k}(p)$, $q^{k}(z,p)$, and $G_{1}^{k}(z,p)$ using Eqs.
(A10) and (A23). The iteration sequence is then repeated until the solution
converges. The downstream escaping flux is then obtained from Eq. (A21) using
$D\frac{\partial f}{\partial z}=\frac{\phi_{\rm esc}}{A(z_{\mathrm{b}})}$.