# Self-Generated Cosmic-Ray Turbulence Can Explain the Morphology of TeV Halos

Payel Mukhopadhyay<EMAIL_ADDRESS>ORCID:
orcid.org/0000-0002-3954-2005 SLAC National Accelerator Laboratory, Stanford
University, Stanford, CA 94039, USA  Physics Department, Stanford University,
Stanford, CA 94305, USA  Tim Linden<EMAIL_ADDRESS>ORCID:
orcid.org/0000-0001-9888-0971 Stockholm University and The Oskar Klein Centre
for Cosmoparticle Physics, Alba Nova, 10691 Stockholm, Sweden

###### Abstract

Observations have shown that spatially extended “TeV halos” are a common (and
potentially generic) feature surrounding young and middle-aged pulsars.
However, their morphology is not understood. They are larger than the
“compact” region where the stellar remnant dominates the properties of the
interstellar medium, but smaller than expected in models of cosmic-ray
diffusion through the standard interstellar medium. Several explanations have
been proposed, but all have shortcomings. Here, we revisit a class of models
where the cosmic-ray gradient produced by the central source induces a
streaming stability that “self-confines” the cosmic-ray population. We find
that previous studies significantly underpredicted the degree of cosmic-ray
confinement and show that corrected models can significantly inhibit cosmic-
ray diffusion throughout the TeV halo, especially when similar contributions
from the coincident supernova remnant are included.

##  I Introduction

TeV halos are a distinct class of galactic $\gamma$-ray emission sources
characterized by their hard $\gamma$-ray spectrum, spatially extended and
roughly spherically symmetric morphology, and coincidence with middle-aged
pulsars. Building on initial observations by the High Altitude Water Cherenkov
(HAWC) telescope of two TeV halos surrounding the Geminga and Monogem pulsars
Abeysekara  _et al._ (2017, 2017), subsequent observations by HAWC, the High-
Energy Steroscopic System (H.E.S.S) and Large High-Altitude Air Shower
Observatory (LHAASO) have identified at least 8 TeV halo systems Linden  _et
al._ (2017); Riviere  _et al._ (2017); Brisbois  _et al._ (2018); Fang  _et
al._ (2019); Sudoh  _et al._ (2019); Fang  _et al._ (2021); Aharonian  _et
al._ (2021). Dozens of other systems have been discussed as possible TeV
halos, or TeV Halo/Pulsar Wind Nebulae (PWN) composite systems Linden  _et
al._ (2017); Linden and Buckman (2018); Profumo  _et al._ (2018); Sudoh  _et
al._ (2019) , a distinction which primarily depends on the definition used to
distinguish standard PWNe from the more extended and diffuse TeV halos Sudoh
_et al._ (2019); Giacinti  _et al._ (2020).

The size and morphology of TeV halos is unexpected. Early theoretical models
by Aharonian and collaborators Aharonian (1995); Aharonian  _et al._ (1995)
predicted that energetic pulsars could produce an extended halo of TeV
electrons. However, these studies utilized standard values for galactic
diffusion, leading them to conclude that the resulting halos would be
extremely large and have a surface brightness too dim to observe.111See also,
later studies based on Milagro observations of Geminga, which noted source
emission Abdo  _et al._ (2007, 2009); Yuksel  _et al._ (2009). Current
observations indicate that the morphology of TeV halos is diffusive – matching
predictions from Aharonian (1995); Aharonian  _et al._ (1995) However, the
diffusion coefficient that was calculated from these observations lies
approximately two orders of magnitude below the standard value for galactic
diffusion Hooper  _et al._ (2017); Abeysekara  _et al._ (2017); Fang  _et al._
(2018); Profumo  _et al._ (2018); Di Mauro  _et al._ (2020, 2019). This
puzzling feature is not explained, and produces halos with a high surface
brightness and typical extension between 20–50 pc at TeV energies.

A number of models have been posited to explain the relatively compact
morphology of TeV halos – all of which have significant drawbacks: The first
assumes that TeV halos occupy regions of space with unusually small pre-
existing diffusion coefficients Abeysekara  _et al._ (2017); Tang and Piran
(2019); Giacinti and López-Coto (2019); López-Coto and Giacinti (2018) that
are not necessarily correlated with the pulsar or associated supernova
remnant. While such a model was reasonable when only Geminga and Monogem had
detected halos – it becomes less credible as more TeV halos are detected. The
second utilizes anisotropic diffusion models which indicate that diffusion is
primarily confined to one dimensional flux tubes on small distance scales that
orient particle diffusion in the direction along the line-of-sight between the
pulsar and Earth Liu  _et al._ (2019). This restricts the spatial extension of
$\gamma$-ray emission perpendicular to the source, but should produce highly
asymmetric sources that are not observed in most cases. A third model exploits
the transition between rectilinear propagation and diffusion to temporarily
oriented high-energy cosmic-rays towards Earth so that they can begin cooling
before diffusing Recchia  _et al._ (2021). However, the efficiency of
diffusion in these scenarios requires the pulsars be significantly more
energetic than indicated by radio observations Bao  _et al._ (2021).

Finally, some models argue that the pulsar (or associated supernova explosion)
generates turbulence that drives down the local diffusion coefficient through
a streaming instability Evoli  _et al._ (2018a); Fang  _et al._ (2019). This
model mirrors similar studies of inhibited diffusion near supernova remnants
D’Angelo  _et al._ (2018), which indicate that the net contribution of all
supernova remnants may even dominate cosmic-ray diffusion throughout the Milky
Way halo Evoli  _et al._ (2018b). In particular, the steep cosmic-ray gradient
produced by a bright source generates Alfvén waves that propagate outward
along the cosmic-ray gradient. Once excited, these Alfvén waves dominate the
turbulence spectrum at the scattering scale because they are naturally
resonant with the injected cosmic rays. The idea of self-generated turbulence
by cosmic ray streaming is not new one and is often taken into account to
model cosmic ray acceleration in supernova remnant shocks and propagation in
the galactic halo Skilling (1970); Lagage and Cesarsky (1983); Aloisio and
Blasi (2013).

Unfortunately, early models of cosmic ray propagation around pulsars with
self-generated turbulence indicated that the pulsar is not energetic enough to
inhibit diffusion on such large scales Evoli  _et al._ (2018a); Fang  _et al._
(2019). Additional contributions from the coincident supernova could
potentially boost this effect. However, mature pulsars such as Geminga have
often moved far from their parent SNR making it difficult to explain why the
TeV emission then remains centered on the quickly moving pulsar Fang  _et al._
(2019).

In this paper we revisit models of cosmic-ray self-confinement. We first
correct an error that affected previous results from Ref. Evoli  _et al._
(2018a). In our corrected models, both the amplitude and duration of cosmic-
ray self-confinement is significantly increased, making pulsars more than
capable of powering cosmic-ray self-confinement in 1D simulations. We show
that these models remain robust even when: (1) additional turbulence damping
terms are included, (2) the pulsar spectrum is softened, or (3) the radius of
the 1D flux tube is increased. In simulations where the initial diffusion
coefficient is unsuppressed in three dimensions, the effect of cosmic-ray
self-confinement decreases too rapidly with increasing radius. However, models
that include additional contributions from the coincident SNR confine cosmic-
rays on 10 pc scales. While fully consistent models of cosmic-ray self-
confinement lie beyond the scope of this work, our models indicate that
pulsars are energetically capable of confining their own cosmic-ray
populations on Myr timescales.

##  II Model of the pulsar source

We model the pulsar as a point-source that begins injecting electrons at time
$t=0$. Throughout this paper we use ‘electrons’ to refer to $e^{+}e^{-}$ pairs
unless otherwise specified. The spindown power ($L(t)$) of the pulsar is
assumed to be:

$L(t)=L(t=0)(1+t/\tau)^{-2}$ (1)

where $L(t=0)$ is the spindown power at $t=0$, $t$ is the time and $\tau$ is
the spindown timescale of the pulsar. The pulsar source function for 1D or 3D
spherically symmetric propagation of cosmic rays is taken to be, respectively:

$\displaystyle
Q^{1D}_{e}(z,p,t)=Q_{0}(t)\left(\frac{p}{m_{e}}\right)^{-\eta}e^{-\frac{p}{p_{m}}}\frac{\exp({-z^{2}/2\sigma^{2}})}{\pi
R^{2}\sqrt{2\pi\sigma^{2}}}$ (2) $\displaystyle
Q^{3D}_{e}(r,p,t)=Q_{0}(t)\left(\frac{p}{m_{e}}\right)^{-\eta}e^{-\frac{p}{p_{m}}}\frac{\exp({-r^{2}/2\sigma^{2}})}{(2\pi\sigma^{2})^{3/2}}$

where $m_{e}$ is the electron mass, $\sigma$ is the source size, $\eta$ is the
spectral index, $p_{m}$ is the momentum cutoff beyond which the source
injection exponentially falls, $R$ is the flux tube radius for the case of 1D
diffusion and $Q_{0}$ is a normalization constant. For the reference case
considered in this paper, we set $\sigma=$ 1 pc in both models and $R=1$ pc in
the 1D case.

The normalization $Q_{0}$ is determined such that the total power in electrons
integrated over momenta is equal to the simultaneous spindown power of the
pulsar. In reality, this term is modified by some efficiency factor that
determines the fraction of the spindown power that is channeled to electrons.
Observations of TeV halo luminosities indicate that both young and middle-aged
pulsars convert a sizable fraction, $\mathcal{O}$(9–30%), of their spindown
power into electrons Hooper  _et al._ (2017); Profumo  _et al._ (2018);
Johannesson  _et al._ (2019); Evoli  _et al._ (2021). For our reference case
of the Geminga, we take $L(340~{}\mathrm{kyr})=3.7\times 10^{34}$ erg/s
Manchester  _et al._ (2005) setting the age to 340 kyr and $\tau=10$ kyr. We
set the efficiency parameter to $\alpha$ = 0.1, corresponding to a 10%
electron injection efficiency.

##  III Cosmic Ray Transport Modeling

The coherence length of the galactic magnetic field is expected to be
somewhere between $L_{c}\sim$ 1–100 pc Iacobelli  _et al._ (2013); Haverkorn
_et al._ (2008); López-Coto and Giacinti (2018). It likely varies
substantially throughout the galaxy, with shorter coherence lengths in regions
with higher stellar and gas densities. Cosmic ray diffusion on distance scales
smaller than the coherence length can be approximated by 1D diffusion in a
flux tube. On larger scales, it is more appropriate to utilize 3D diffusion
models. In this paper, we separately consider cosmic ray transport for both 1D
and 3D geometries, ignoring the transition from 1D to 3D propagation for
diffusion distances near $L_{c}$. In this section, we first detail a 1D
treatment of diffusion, and then subsequently discuss several adjustments
needed to consider particle propagation in 3D.

###  III.1 Diffusion Equation in 1D

Particle diffusion in a cylindrically symmetric 1D flux tube can be described
by the following transport equation:

$\displaystyle\frac{\partial f}{\partial t}=Q_{e}(p,z,t)-u_{A}\frac{\partial
f}{\partial z}+\frac{\partial}{\partial z}\left(D(p,z,t)\frac{\partial
f}{\partial z}\right)$ $\displaystyle+\frac{\partial u}{\partial
z}\frac{p}{3}\frac{\partial f}{\partial
p}-\frac{1}{p^{2}}\frac{\partial}{\partial p}\left(p^{2}\frac{dp}{dt}f\right)$
(3)

where $f$ is the phase space distribution function, $u_{A}$ is the advection
speed, which we assume to be the Alfvén velocity given by
$v_{A}=B_{0}/\sqrt{4\pi m_{i}n_{i}}$ ($B_{0}$ is the background magnetic
field, $m_{i}$, $n_{i}$ represent the mass and number density of ions in the
plasma), and $D(p,z,t)$ diffusion coefficient, which varies in momentum, space
and time.

###  III.2 Energy Loss Terms

The last term, which includes $\frac{dp}{dt}$, takes into account the energy
losses of $e^{+}e^{-}$ pairs. The energy loss rate is due to a combination of
synchrotron radiation and inverse Compton scattering Blumenthal and Gould
(1970) and is given by:

$\displaystyle-\frac{dp}{dt}=\sum_{i}\frac{4}{3}\sigma_{T}\rho_{i}S_{i}(E)\left(\frac{E_{e}}{m_{e}}\right)^{2}+\frac{4}{3}\sigma_{T}\rho_{mag}\left(\frac{E_{e}}{m_{e}}\right)^{2}$
(4)

where $\sigma_{T}$ is the Thomson cross section, $\rho_{mag}$ is the magnetic
field energy density and the sum over in the expression with $\rho_{i}$ is
carried out over the various radiation backgrounds consisting of the cosmic
microwave background (CMB), infrared radiation (IR), ultraviolet emission (UV)
and starlight. We choose $\rho_{CMB}=0.260$ eV/cm3, $\rho_{IR}$ = 0.6 eV/cm3,
$\rho_{star}=0.6$ eV/cm3 and $\rho_{mag}$ = 0.025 eV/cm3 for a background
magnetic field of $B_{0}$ = 1 $\mu G$ Hooper  _et al._ (2017). At very high
energies, inverse Compton scattering is further suppressed by Klein-Nishina
effects (encoded in $S_{i}(E)$) which have also been taken into account
following the analytic formalism of Schlickeiser and Ruppel (2010), the
accuracy of which has recently proven insufficient for the 1% measurements of
electron spectra by AMS-02 Fang  _et al._ (2021); Di Mauro  _et al._ (2021),
but is more than sufficient for our purposes here.

###  III.3 Self-Generated Turbulence and Non-Linear Damping

The confinement of cosmic rays in our model depends on the generation of
magnetic turbulence by the sharp gradient produced in the cosmic-ray
population. Specifically, the diffusion coefficient ($D$) is related to the
spectral power $W$ as given in Ref.Skilling (1971)

$D(p,z,t)=\frac{4}{3\pi}\frac{r_{L}(p)v(p)}{k_{res}W_{res}}$ (5)

where $r_{L}(p)$ is the Larmor radius and $v(p)$ is the particle speed, which
we assume to be equal to the speed of light $c$. $k_{res}=1/r_{L}(p)$ is the
resonant wavenumber and $W_{res}$ is the spectral power computed at the
resonant wavenumber. The evolution of $W$ can then be calculated by:

$\frac{\partial W}{\partial t}+v_{A}\frac{\partial W}{\partial
z}=(\Gamma_{CR}-\Gamma_{D})W$ (6)

where the wave growth rate $\Gamma_{CR}$ for 1D propagation is given by:

$\Gamma_{CR}=\frac{16\pi^{2}}{3B_{0}^{2}}\frac{v_{A}}{kW(k)}\left(p^{4}v(p)\frac{\partial
f}{\partial z}\right)_{p_{res}}$ (7)

with $\Gamma_{CR}$ being the growth rate of Alfvén waves due to the resonant
streaming instability Skilling (1971).

The turbulence will also be damped by a combination of several different
processes. Our default models focus on non-linear wave damping (NLD) terms,
which have been found to be the dominant damping term in regimes with similar
magnetic field strengths and gas densities as observed near Geminga Evoli  _et
al._ (2018a). However, in the next section we will also discuss the impact of
additional damping terms. For the NLD term, we follow the prescription of
Ptuskin and Zirakashvili (2003) and set:

$\displaystyle\begin{split}\Gamma_{NLD}&=(2c_{K})^{-3/2}|v_{A}|k^{3/2}W^{1/2}\;\;(\mathrm{Kolmogorov})\\\
\\\
&=(2c_{K})^{-3/2}|v_{A}|k^{2}W\;\;\;\;\;\;\;\;\;\;(\mathrm{Kraichnan})\end{split}$
(8)

where $c_{K}=3.6$. Eq.(8) is equivalent to having a Kolmogorov or Kraichnan
description of the cascade given by $\frac{\partial}{\partial
k}\left(D_{kk}\frac{\partial W}{\partial k}\right)$ in the limit where the
cascade damping timescale is faster than any other damping timescales Brunetti
_et al._ (2004). We note that it was this term in Ref. Evoli  _et al._ (2018a)
which included a numerical error, where cK was accidentally set to be 0.21
(though the value was correctly listed in the paper). The error in this value
significantly increased (by more than an order of magnitude) the rate of non-
linear damping, significantly decreasing both the duration and the maximum
inhibition of cosmic-ray diffusion throughout that work.

###  III.4 Additional Damping Processes

In our default models, we consider only the NLD processes for turbulent
damping. However, several additional dynamic processes can counteract the
generation of magnetic turbulence and relax the diffusion coefficient to
standard ISM values. Most notable is the ion neutral damping (IND) mechanism
which occurs in partially ionized plasmas Piddington (1956); Kulsrud and
Pearce (1969). The IND mechanism dissipates wave energy through a viscosity
term produced by charge-exchange interactions between ions and neutral atoms
which cause neutral particles to begin participating in hydromagnetic
phenomena.

The damping rate depends on how the wave frequency $\omega_{k}=v_{A}k$
compares to the ion-neutral collisional frequency $\nu_{in}=8.4\times
10^{-9}\times(n_{n}/cm^{-3})\times(T/10^{4}~{}K)^{0.4}$ s-1, where $n_{n}$ is
the neutral atom number density and $T$ is the plasma temperature. The IND
damping rate also depends on the ratio of ion to neutral densities
$f_{i}=n_{i}/n_{n}$. In the limit $f_{i}>>1$, which best describes
interstellar regions that lie outside of dense molecular clouds, the IND
damping rate is well approximated by Ref. Zweibel and Shull (1982):

$\Gamma_{IND}\simeq\frac{\nu_{in}}{2}\frac{\omega_{k}^{2}}{\omega_{k}^{2}+f_{i}^{2}\nu_{in}^{2}}$
(9)

For our default model, we assume $n_{i}=1$ cm-3 and $n_{n}\sim 10^{-2}$ cm-3,
representing a neutral fraction of 1%. These values are well-motivated by
observations indicating the absence of H$\alpha$ emission from Geminga’s bow
shock nebula, which indicate that the neutral fraction surrounding the pulsar
is less than 1% Caraveo  _et al._ (2003). In this case, we find that the IND
damping timescale at 10 TeV is $\frac{1}{\Gamma_{IND}}>10^{6}$ kyr, a value
that significantly exceeds the NLD damping timescale of
$\frac{1}{\Gamma_{NLD}}\sim O(100~{}\mathrm{kyr})$. Based on this line of
reasoning, the IND term was ignored in Ref. Evoli  _et al._ (2018a) and is
similarly ignored in our default analysis.

However, we note that the IND term can be significantly more important at low
energies, or in regions with higher neutral gas densities. In order to
determine the diffusion coefficient over the full energy range, we thus
include IND damping for some simulations using neutral fractions of 1% and
10%, finding that such models significantly increase the diffusion coefficient
at low energies.

In addition to the NLD and IND damping mechanisms, additional damping terms
such as Nonlinear Landau damping (NLLD) and linear turbulent damping (also
known as the Farmer-Goldreich: FG damping), have been considered in the
literature, but are not included in this analysis. NLLD damping effectively
converts wave energy into heat Kulsrud (1978); Wiener  _et al._ (2013). If
this damping term dominates, the turbulent cascade is suppressed and the NLD
cascade damping does not develop. However, the effectiveness of this damping
term in a turbulent plasma is not clear. The functional form for the damping
Kulsrud (1978); Wiener  _et al._ (2013) is obtained from a simplified system
where two Alfvén waves (with slightly different wavelengths) travel in
opposite directions, forming beat waves that interact with the thermal ions
and absorb the wave energy. A critical assessment of the role of this
mechanism, using for example, Particle-in-cell (PIC) simulations, is still
missing Blasi (2019). Additionally, we note that if NNLD damping were dominant
through the Galactic volume, then the turbulence injected by Supernova
remnants would be inefficient to explain the diffusion coefficient observed in
the Galaxy Ptuskin and Zirakashvili (2003).

The turbulent or FG damping term is given by
$\Gamma_{FG}\sim\frac{v_{A}}{\sqrt{L_{inj}r_{L}}}$, where $v_{A}$ is the
Alfvén speed and $L_{inj}$ is the injection scale of the turbulence Farmer and
Goldreich (2004). For reasonable magnetic field strengths around Geminga,
$B\sim 1$ $\mu G$ Liu  _et al._ (2019), the timescale associated to this
damping term is $\sim 500$ kyr for $p\sim 10$ TeV particles, which is
generally longer than the damping timescales associated to the NLD damping.
Therefore, for the models considered in this paper, FG damping is expected to
be subdominant. We note, however, that FG damping may become relevant for TeV
halos formed in regions with larger interstellar magnetic fields.

###  III.5 Adjustments for Diffusion in 3D

As noted earlier, for cosmic ray diffusion on distance scales larger than the
galactic coherence length, the diffusion is expected to be better represented
by particle propagation in 3D – although diffusion can be more complicated in
the transition region where TeV halos may reside. In this subsection, we
consider several adjustments to our model which are necessary to consistently
model diffusion in a regime where particle propagation in all three spatial
dimensions is initially unsuppressed. We additionally note that TeV halo
observations indicate that the halo is roughly spherical (or at least, not a
linear tube), which strengthens the case for considering 3D diffusion models.
We first note that in Eq. II, we have provided the pulsar morphology in both
the 1D and 3D spherically symmetric case, but find that these differences do
not significantly impact the results of our study.

More importantly, spherically symmetric particle diffusion in 3D space is
governed by the 3D diffusion equation, given by:

$\displaystyle\frac{\partial f}{\partial
t}=Q^{3D}_{e}(p,r,t)-v_{A}\frac{\partial f}{\partial
r}+\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}D(p,r,t)\frac{\partial
f}{\partial r}\right)$ (10)
$\displaystyle+\frac{2u_{A}}{r}\frac{p}{3}\frac{\partial f}{\partial
p}-\frac{1}{p^{2}}\frac{\partial}{\partial p}\left(p^{2}\frac{dp}{dt}f\right)$

Finally, in addition to changes in the diffusion term, the wave growth rate in
a spherically symmetric 3D geometry is given by:

$\Gamma^{3D}_{CR}=\frac{16\pi^{2}}{3B_{0}^{2}}\frac{v_{A}}{kW(k)}\left(p^{4}v(p)\frac{\partial
f}{\partial r}\right)_{p_{res}}$ (11)

where essentially the $\frac{\partial f}{\partial z}$ in eq. 7 is replaced by
the gradient in 3D spherical geometry Ptuskin  _et al._ (2008).

###  III.6 Contributions from a Supernova Remnant

Finally, we study a scenario where, in addition to electrons injected by the
pulsar, a supernova remnant (SNR) also injects protons into the interstellar
medium (see Fang  _et al._ (2019)). In this case, the self-generation of waves
is governed by the combined distribution function of both protons and
electrons. The SNR source term is taken to be:

$Q^{SNR}(r,p,t)=Q_{0}(t)\left(\frac{p}{m_{p}}\right)^{-\eta_{SNR}}e^{-\frac{p}{p_{c}}}\\\
\times\frac{\exp{(-r^{2}/2\sigma_{SNR}^{2}})}{(2\pi\sigma_{SNR}^{2})^{3/2}}$
(12)

where $Q_{0}$ is given by:

$Q_{0}(t)=\frac{\xi_{CR}E_{SN}}{T_{SN}I}$ (13)

where $T_{SN}$ is the time over which the SNR provides injection power. For
this paper, we assume $T_{SN}=1000$ yr with the SNR injection stopping after
this time. $I$ is a constant obtained by equating the total energy in protons
to the fraction of SNR kinetic energy transmitted to protons. $\xi_{CR}$ is
the fraction of the supernova kinetic energy transmitted to cosmic rays.

Protons and electrons are treated separately and for each component, the
propagation equation (Eq. 10) is solved. Protons are assumed to have zero
energy losses while electrons are taken to lose energy according to Eq. 4.

##  IV Numerical Modeling

The diffusion of particles and growth/damping of turbulence represent a set of
coupled differential equations which we solve using an explicit finite-
difference scheme. The timestep $\Delta t$ and spatial gridstep $\Delta x$ are
chosen such that the system is stable. In this paper we choose $\Delta t$ =
0.05 yr and $\Delta x$ = 1 pc. The momentum grid has 120 logarithmic
gridpoints between $p_{min}=1$ GeV and $p_{max}=10^{6}$ GeV. The results have
been confirmed for a range of spatial and momentum grids. The validity of the
results has also been tested in an unconditionally stable semi-explicit Crank
Nicolson scheme.

Figure 1: The diffusion coefficient at a radius of 10 pc and an energy of 10
TeV as a function of the pulsar age in the 1D case. Two example scenarios, of
Kolmogorov and Kraichnan turbulence, are shown. The shaded bands show results
for magnetic field strengths ranging from 0.5–2 $\mu G$. In both cases, a
strong suppression of the diffusion coefficient ($\sim 2-3$ orders of
magnitude) is observed. The suppression is seen to persist over the Myr
timescale in our simulation.

For our 1D fiducial model, we use $\alpha$ = 0.1, $B_{0}=1$ $\mu G$, $\sigma$
= 1 pc and set the background turbulence to have a Kolmogorov spectrum with
$D=3.466\times 10^{28}~{}p_{\mathrm{GeV}}^{1/3}$ cm2 s-1. We take the
background ion number density $n_{i}=1~{}\mathrm{cm^{-3}}$. We set the
momentum dependence of the electron injection spectrum to be $\propto
p^{-3.5}\exp(-\frac{p}{100~{}\mathrm{TeV}})$ over a range between 1 GeV and
107 GeV. The flux tube radius, $R$ is assumed to be 1 pc. The outer boundary
of the flux tube is taken to be $z_{b}=500$ pc, an unrealistically large
number chosen such that boundary conditions do not affect our result near the
pulsar. We set the electron distribution, $f$, to be 0 at the outer boundaries
of the flux tube. The initial condition is set to be $f(z,p)=0$ at $t=0$.

For the 3D fiducial model, we keep most parameters to be the same as the
fiducial 1D model, but set the default pulsar efficiency to be $\alpha$ = 1.0.
The spatial gridstep is set to be $\Delta r=1$ pc. The 3D boundary condition
is set so that the electron distribution $f$ is zero at the outer boundary of
the spherical volume considered, at $r_{b}=500$ pc. The inner boundary of the
simulation is set at $r_{in}=0.1$ pc, instead of $r=0$ to avoid singularities.
For the inner boundary condition, we set $f(r_{0},p)=f(r_{1},p)$, where
$r_{0}$ and $r_{1}$ denote the first two points of the radial grid.

For 3D models that include an SNR, we set the SNR kinetic energy to be
$E_{SN}=10^{51}$erg, the SNR size, $\sigma_{SNR}=1$ pc, SNR spectral index,
$\eta_{SNR}=4.2$ and fraction of energy transmitted to protons,
$\xi_{CR}=0.2$. The electron distribution function ($f_{e}$) and proton
distribution function ($f_{p}$) are evolved separately with identical boundary
conditions in each case.

##  V Results

###  V.1 1D model

Figure 2: The inhibition of the diffusion coefficient at an energy of 10 TeV
and a distance of 10 pc from the location of the pulsar as a function of time.
Curves correspond to different efficiencies in the conversion of the pulsar
kinetic energy into e+e- acceleration. Models with high efficiencies produce
diffusion coefficients that fall below the Bohm limit, at which point our
diffusion model may no longer be valid. In our models, we do not allow the
diffusion coefficient to drop below the Bohm limit at any gridpoint.

In Fig. 1, we show the time evolution of of the diffusion coefficient for
models utilizing the Kolmogorov and Kraichnan wave-damping phenomenology in
our 1D simulation. We note that the initial turbulence spectrum is set to the
Kolmogorov values in each case. We show results for values of the magnetic
field strength ranging between 0.5–2 $\mu G$. The injection efficiency assumed
for this plot $\alpha=0.1$. The result is plotted at an electron energy of 10
TeV and at a distance 10 pc from the pulsar, roughly corresponding to the
energy and distance scales corresponding to HAWC observations of the Geminga
TeV halo Abeysekara  _et al._ (2017). For both the Kolmogorov and Kraichnian
cases, the diffusion coefficient decreases by more than 2.5 orders of
magnitude compared to the background ISM value, reaching minimum values
between 100–300 kyr, consistent with the ages of the Geminga and Monogem
halos. Notably, the suppressed diffusion coefficient persists throughout the 1
Myr timescale of our simulation.

The suppression of the diffusion coefficient is governed by the growth of
Alfvén waves governed by $\Gamma_{CR}$, which induces a period of highly
inhibited diffusion. The subsequent increase in the diffusion coefficient at
late times is governed by the relaxation term $\Gamma_{D}$. Our models
indicate that the diffusion coefficient does not return back to ISM values
within the 1 Myr timescales of our simulation for either the Kolmogorov of
Kraichnan models, with Kraichnan models having an even longer relaxation
timescale Evoli  _et al._ (2018a).

Our results are markedly different from those obtained in Ref. Evoli  _et al._
(2018a) (See Fig. 1 of Ref. Evoli  _et al._ (2018a) for a direct comparison.)
While both studies found that cosmic-ray self-confinement could inhibit
diffusion throughout the halo by 2.5–3 orders of magnitude at 10 TeV, our
study requires only a 10% conversion of the pulsar spindown power into
electron pairs to achieve this effect, while Ref. Evoli  _et al._ (2018a)
required a 100% conversion of the pulsar power. More importantly, our models
indicate that the diffusion coefficient requires more than 1 Myr to rebound to
ISM values, while the relaxation time in Ref. Evoli  _et al._ (2018a) was much
shorter and the diffusion coefficient returned to ISM values within 100 kyr.
We have confirmed with the authors that the discrepancy between these results
is due to an error in the relaxation term computed in Ref. Evoli  _et al._
(2018a). The results presented here have been checked against a corrected
version of the code used in Evoli  _et al._ (2018a), finding good agreement.

Figure 3: The diffusion spectrum as a function of energy for different time
snapshots for our default model with an injection efficiency of $\alpha=0.1$.
The grey dashed line indicates the Bohm limit. The diffusion coefficient is
highly suppressed below $\sim$ 100 TeV and the suppression persists till 1
Myr. The diffusion coefficient rises rapidly for very high energy electrons
above $\sim$ 10 TeV due to a combination of the exponential cutoff in the
pulsar spectrum and very fast cooling of high-energy electrons.  Figure 4: A
heatmap of the diffusion coefficient as a function of distance from the pulsar
and age for our default 1D scenario with $\alpha=0.1$. The diffusion
coefficient is found to be suppressed by $\sim$ 2–3 orders of magnitude within
10–20 pc of the pulsar. The resulting morphology of the halo is consistent
with HAWC observations. However, the diffusion remains inhibited out to more
than 100 pc from the pulsar, where the 1D approximation likely breaks down.

In Fig. 2, we show the evolution of the diffusion coefficient as a function of
time at a constant distance of 10 pc, but for different source injection
efficiencies. In our fiducial model, using a 10 % efficiency, the diffusion
coefficient for 10 TeV electrons is suppressed by more than two orders of
magnitude at 340 kyr. For values of the injection efficiency as low as 3%, we
still obtain a significant (factor of $\sim$5) suppression in the diffusion
coefficient, indicating that even relatively weak pulsars could significantly
affect local cosmic-ray propagation. We note that these effects are non-
linear. A model with a 10% pulsar efficiency produces a decrease in local
diffusion that is 50$\times$ larger than a model with a 1% pulsar efficiency.
This is due to the fact that the growth of turbulence is resonant. Larger
cosmic-ray gradients lead to suppressed diffusion, which in turn causes more
cosmic-rays to pile-up (leading to larger gradients). As a result, even small
increases the pulsar power can significantly affect the resulting model.

For models with injection efficiencies exceeding $\sim$30%, i.e $\alpha=0.3$,
the diffusion coefficient at 10 pc and 10 TeV falls below the Bohm limit. We
note that the assumptions of resonant diffusion off of turbulent waves breaks
down in such an environment. Thus, in our default models we set the minimum
diffusion coefficient in any spatial and energy bin to be equal to the Bohm
value at that energy. We note that even if the curve Fig 2 does not show a
value at the Bohm limit, the simulation may still be affected by this
parameter choice if the Bohm limit is being reached in a radial bin closer to
the pulsar.

In Fig. 3, we show the diffusion spectrum generated by our default model
($\alpha$ = 0.1) over the full energy range of our simulation. We immediately
note several key results: (1) the diffusion coefficient is highly suppressed
below energies of $\sim$100 TeV, (2) the minimum value of the diffusion
coefficient is reached on timescales of $\sim$100 kyr across a wide-energy
range, (3) the diffusion coefficient rises rapidly above $\sim$10 TeV, due to
the significantly decreasing pulsar power (and highly efficient electron
cooling) at higher energies. This result is potentially testable with upcoming
HAWC and LHAASO observations, (4) the diffusion coefficient below $\sim$10 TeV
is nearly energy-independent between 100-300 kyr, closely matching
observations of energy-independent diffusion in the Geminga TeV halo that were
first noted by Hooper  _et al._ (2017). Notably, these spectral features may
help differentiate cosmic-ray self confinement models from alternative
explanations for inhibited diffusion within TeV halos.

In Fig. 4, we show a two-dimensional representation of the diffusion
coefficient at our standard energy of 10 TeV as a function of both time (after
pulsar formation) and distance (from the pulsar) for our default simulation
(with $\alpha$=0.1). At time $t$=0, the diffusion coefficient is high
throughout the simulation volume, but it drops quickly in regions near the
pulsar, achieving a minimum value of $\sim$1027 cm2s-1 at an age of $\sim$100
kyr. The diffusion coefficient can continue to be significantly inhibited and
much larger radii – falling below a value of $\sim$1028 cm2s-1 (70$\times$
lower than the standard ISM) out to a radius of $\sim$35 pc. Such a model
matches (and even exceeds) the inhibition of diffusion observed near the
Geminga and Monogem TeV halos Abeysekara  _et al._ (2017), indicating the
potential impact of cosmic-ray self-confinement in these simulations.

###  V.2 Robustness of Results

Our 1D models indicate that the diffusion coefficient near pulsars: (1) can be
suppressed by 2–3 orders of magnitude at critical energies near 10 TeV, (2)
will remain suppressed on 1 Myr timescales consistent with the observation of
TeV halos around mature pulsars, (3) can be suppressed even if only $\sim$10%
of the pulsar spindown power is converted into electrons. The surprising
efficiency of cosmic-ray self-confinement indicates that these models are
likely to remain robust even if less optimistic input parameters are
considered.

Figure 5: The diffusion coefficient as a function of time at an energy of 10
TeV and a distance of 10 pc from the pulsar. The plot is shown for a softer
pulsar spectral index, $\eta=4.2$ than the one considered in Fig. 2, which
uses $\eta=3.5$. We find that even with an extremely soft spectral index, the
diffusion coefficient decreases by $\sim$ 2 orders of magnitude for a pulsar
efficiency of 30%, i.e $\alpha=0.3$. The pulsar efficiency needed to inhibit
the diffusion coefficient by 2 orders of magnitude is higher than our default
model with $\eta=3.5$ because the efficiency must increase to compensate for
the softer spectrum. Figure 6: The diffusion coefficient as a function of
energy in a model where we include IND damping. The result is shown at 300 kyr
and 10 pc from a pulsar with an efficiency of $\alpha$=0.1. The ion number
density, $n_{i}$ is set to 1 cm-3. We show results for neutral atom number
densities, $n_{n}$ = 0.1, 10-2 and 10-3 cm-3 (corresponding to a 10%, 1% and
0.1 % neutral fraction). For a 0.1% neutral fraction, IND damping has a
negligible impact on our model and the diffusion spectrum at 300 kyr very
similar to Fig. 3. When the neutral fraction increases to 1%, the low energy
($\lesssim 1$ TeV) flux begins to be significantly affected. However, the
diffusion coefficient at 10 TeV is not impacted for any of the neutral
fractions considered.

In Figure 5 we show the diffusion coefficient as a function of time in a model
where the pulsar electron injection spectrum has a momentum index p-4.2, which
is significantly softer than our default value of p-3.5. Our default value is
more consistent with observations of the Geminga and Monogem TeV halos Hooper
_et al._ (2017), while a softer value may be motivated by GeV observations of
Geminga Di Mauro  _et al._ (2019) as well as the spectrum of diffuse emission
from the TeV halo population Linden and Buckman (2018). However, we stress
that our choice of a p-4.2 momentum spectrum is even more pessimistic than
advocated by the results of Di Mauro  _et al._ (2019); Linden and Buckman
(2018), as we assume that the soft cosmic-ray injection spectrum continues all
the way to 1 GeV, while previous models have utilized a broken power-law to
decrease the electron power below energies of $\sim$500 GeV. Even in this
extremely pessimistic case, we find that the diffusion coefficient at 10 TeV
can be inhibited by 2–3 orders of magnitude so long as the electron efficiency
is increased to offset the softened spectrum.

In Figure 6, we show the diffusion spectrum for models that additionally
include ion-neutral damping (IND). The total damping term is taken as the sum
of the IND and NLD terms, as defined in Section III. We find that the IND term
generally has no effect at TeV energies, but can considerably affect diffusion
in the GeV range. This result is intriguing, as it indicates that TeV halos in
regions with low neutral gas densities (such as Geminga) may also produce GeV
halos Di Mauro  _et al._ (2019), while TeV halos in dense regions of the ISM
may lack GeV counterparts. This observation appears unique to cosmic-ray self-
confinement models and could be tested by future Fermi-LAT data.

Finally, in Figure 7 we investigate models with significantly wider flux tube
radii. Within a 1D model, such a change is degenerate with changing the
efficiency ($\alpha$) of the pulsar, because the electrons are evenly
distributed throughout the radius of the flux tube. However, this degeneracy
breaks down when the pulsar efficiency approaches unity. Further increases in
the flux tube radius can not be compensated by changing the pulsar energetics,
and the degree of cosmic-ray self-confinement begins to decrease. However,
even for flux tubes with a radius of 5 pc, we can obtain a diffusion
coefficient that is suppressed by more than an order of magnitude at a
distance 10 pc. Such a model provides a hint regarding the transition from 1D
to 3D models – indicating that in regimes where the flux tube approximation
begins to break down, the pulsar may continue to significantly inhibit local
diffusion.

Figure 7: The diffusion coefficient as a function of time for different flux
tube radii $R$ and source injection fractions $\alpha$. Increasing the flux
tube radius lowers the suppression of the diffusion coefficient by increasing
the propagation volume. This can initially be compensated for by increasing
the injection fraction $\alpha$. For example, the black and yellow solid lines
correspond to flux tube radii of 1 pc and 3 pc respectively, but have an
identical diffusion coefficient because the increasing radius is compensated
by increasing $\alpha$ from 0.1 to 0.9. Figure 8: Evolution of the diffusion
coefficient as a function of the pulsar age for 3D diffusion models. We assume
a 100% pulsar injection efficiency ($\alpha=1.0$). Results are shown at 5 pc
and 10 pc for two models: the fiducial 3D model and a ‘tuned’ model adjusted
to have extreme parameters. The fiducial model has $B=1\mu G$ and a spectral
index $\eta=3.5$. The tuned model has a spectral index of $\eta=3.0$, a
magnetic field of $B=0.5\mathrm{\mu G}$ and a minimum pulsar power of 10 GeV.
The 3D model produces a very low suppression (less than a factor of 3) of the
diffusion coefficient at 10 TeV for 10 pc distance. However, this model is
capable of producing very inhibited diffusion coefficients within the inner
$\sim$ 5 pc.

In view of the results above, we conclude that, within the context of 1D
simulations, the self-generation mechanism can be highly efficient in
suppressing the diffusion coefficient surrounding TeV halos by 2–3 orders of
magnitude at distances between 10–20 pc from the pulsar. Moreover, the
diffusion coefficient is maximally suppressed on timescales of $\sim$100 kyr,
consistent with the observations of extremely suppressed ISM diffusion in the
TeV halos that surround multiple middle-aged pulsars Linden  _et al._ (2017);
Riviere  _et al._ (2017); Brisbois  _et al._ (2018); Fang  _et al._ (2019);
Sudoh  _et al._ (2019); Fang  _et al._ (2021); Aharonian  _et al._ (2021). We
also note that diffusion remains relatively efficiency within the first
$\sim$10 kyr, which may explain the lack of an observed TeV halo around the
Crab nebula and the observation of only a dim halo around Vela (although we
note that the PWN, which is not modeled in this analysis, also plays a
significant role in electron propagation at such an early stage).

###  V.3 3D model

We note that in 1D models (where the volume increases linearly with distance),
the diffusion coefficient remains inhibited more than 100 pc from the pulsar.
Such a scenario appears unphysical, as the correlation length of galactic
magnetic field turbulence likely falls below 100 pc, even in particularly low-
density regions of the galaxy Iacobelli  _et al._ (2013); Haverkorn  _et al._
(2008). Additionally, the spherically symmetric nature of several TeV halo
observations motivates the investigation of 3D models.

In Figure 8, we show the evolution of the diffusion coefficient as a function
of pulsar age for our 3D case with two different model parameters. One is the
3D fiducial parameters outlined in IV where the pulsar injection efficiency,
$\alpha=1.0$, $B=1~{}\mu G,\sigma=1$ pc, $\tau=10$ kyr and the background ISM
diffusion coefficient is $D=3.466\times 10^{28}p^{1/3}$ cm2 s-1. The momemtum
dependence of the pulsar injection for this fiducial model is $\propto
p^{-3.5}\exp{(-\frac{p}{100~{}\mathrm{TeV}})}$. Additionally, we consider a
“tuned” model where we adjust several parameters to optimistic values that
give the maximum suppression in the diffusion coefficient, namely
$B=0.5\mathrm{\mu G}$, $\tau=5$ kyr, a minimum injection momentum $p_{min}=10$
GeV and a momentum dependence of the pulsar injection $\propto
p^{-3.0}\exp{(-\frac{p}{100~{}\mathrm{TeV}})}$. For both the cases, we plot
the diffusion coefficient as a function of time at both 5 pc and 10 pc from
the pulsar. We note that this plot is similar to Figure 2 in the 1D case.

Figure 9: The diffusino coefficient as a function of pulsar age at 10 pc and
10 TeV for the case of 3D propagation with SNR proton injection. Results for
different SNR injection power into protons (denoted by SNRp) are shown. We set
$B=1\mu G$ and the SNR spectral index, $\eta_{SNR}=4.2$. In order to suppress
diffusion at 10 pc by more than one order of magnitude, the SNR must convert 4
$\times 10^{50}$ erg (or 40% of the SNR kinetic energy) into protons, which
exceeds standard SNR models.

Focusing on the black solid and dashed plots that correspond to the fiducial
model, we find that the diffusion coefficient at 10 pc is suppressed by less
than a factor of $2$ over the 1 Myr timescale of our simulation. However, at 5
pc, the diffusion coefficient is suppressed by more than 1.5 orders of
magnitude until 1 Myr. We find that even our extremely “tuned” model does not
significantly inhibit cosmic-ray diffusion at 10 pc (decreasing the diffusion
coefficient by only a factor of $\sim$3). However, it drastically decreases
the diffusion coefficient (by more than two orders of magnitude) at 5 pc.
Thus, our simulations indicate that in 3D models, the pulsar can significantly
inhibit cosmic-ray diffusion on $\sim$5 pc scales, but is insufficiently
powerful to inhibit diffusion on the 10–20 pc scales that are consistent with
TeV halo observations Abeysekara  _et al._ (2017).

###  V.4 3D model with Supernova Remnant

Because cosmic-rays from the pulsar appears to be insufficient to
significantly decrease the diffusion coefficient on 10–20 pc scales, we
consider the possibility that protons from the coincident supernova may
additionally contribute to lowering the diffusion coefficient through a
similar process. Utilizing an identical pulsar model as described for the
‘default” 3D model ($\alpha$ = 1.0), we add an SNR with a total kinetic energy
of 1051 erg. The momentum dependence of the proton injection spectrum by the
SNR is taken to be $\propto p^{-4.2}\exp(-\frac{p}{10^{7}~{}\mathrm{GeV}})$.
Protons from the SNR are treated identically as the electrons in our default
model, with the exception that they do not cool (because hadronic losses are
inefficient on these time- and density scales).

In Figure 9, we show the results for this model at our default radius of 10 pc
and energy of 10 TeV for different assumptions regarding the fraction of the
SNR kinetic energy that is converted into TeV protons. In all cases, we find
that the SNR significantly increases the suppression of the local diffusion
coefficient. However, a relatively high conversion (40% of the total SNR)
power must be converted into cosmic-ray protons to produce the factor of 30
suppression in the diffusion coefficient that begin to make these models
consistent with TeV halo observations. The full two-order of magnitude
suppression of the diffusion coefficient is only possible if the entire
kinetic energy of the SNR is converted into protons.

Figure 10: Colormap for $D$ as a function of distance from the pulsar and age.
This plot assumes an optimistic $\xi_{CR}=0.4$. The diffusion coefficient can
be highly suppressed in the inner $\sim$ 5 pc. The suppression is about an
order of magnitude at $\sim$ 10 pc. The diffusion coefficient quickly rises
back to its ISM values after 10 pc.

Finally, in Figure 10, we show a two-dimensional heatmap of the diffusion
coefficient as a function of time and distance from the pulsar in our 3D model
with a supernova contribution. We set the conversion efficiency of proton
kinetic energy into protons to $\xi_{p}=0.4$. These results demonstrate two
key (and intuitive) differences between our 1D and 3D models: (1) in 3D models
that include a significant SNR components, the diffusion coefficient near the
source becomes extremely small, saturating the Bohm limit for regions that are
within 5 pc of the pulsar on timescales exceeding 100 kyr. In addition to
demonstrating that SNR activity can significantly affect local diffusion -
this indicates that our simulation may not be robust, as it is difficult to
utilize our diffusive model under such extreme conditions, and (2) unlike the
case of 1D diffusion, the quickly expanding volume of 3D diffusion localizes
the effect of cosmic-ray induced self-confinement to a region within $\sim$10
pc of the central source. The diffusion coefficient rapidly returns to
standard ISM values by $\sim$20 pc from the pulsar.

##  VI Discussion and Conclusions

In this paper, we have re-examined models of cosmic-ray induced self-
confinement. We have corrected an error in the previous study of Ref. Evoli
_et al._ (2018a), which utilized an unphysically large non-linear damping
term. The net effect of this change demonstrates that pulsars are
significantly more capable of inhibiting cosmic-ray diffusion by more than a
factor of 100, consistent with TeV halo observations. Moreover, once diffusion
is suppressed by pulsar activity, it remains inhibited on nearly Myr
timescales, consistent with the observation of TeV halos around mature
pulsars. These two findings paint a much more optimistic picture regarding the
potential for pulsars to confine their own cosmic-ray population, compared
with previous work.

Excitingly, our models indicate several testable features that are unique to
the self-confinement scenario. These may be used to distinguish self-
confinement models from alternative models with pre-existing low diffusion
coefficients Abeysekara  _et al._ (2017); Tang and Piran (2019); Giacinti and
López-Coto (2019); López-Coto and Giacinti (2018), anisotropic diffusion
models Liu  _et al._ (2019) and models incorporating the transition from
rectilinear to diffusive propagation Recchia  _et al._ (2021). First, our
model predicts that in the absence of IND damping, the diffusion coefficient
will be roughly energy independent for $p\lesssim 1$ TeV, and the diffusion
coefficient will rapidly rise near the exponential energy cutoff of the pulsar
(Fig. 3). Second, the presence of stronger IND damping in dense regions of the
ISM will prevent the observation of GeV emission from TeV halos (Fig. 6) —
producing a bifurcation of TeV halo observations that depends on the local gas
density. Third, our models indicate that standard TeV halos should not form
around systems much younger than $\sim$10 kyr. All of the above predictions
are exclusive to this mechanism and are potentially testable with future TeV
observations.

###  VI.1 Comparison of 1D and 3D Models

While our results indicate that cosmic-ray self confinement models are robust
to changes in the initial pulsar spectrum, the presence of IND damping, and
the width of the 1D flux tube radius, our results strongly depend on whether
particle diffusion is initially confined to 1D flux tubes, or initially
diffuses efficiently in all three spatial directions.

In 1D models, many details of the self-generation mechanism are consistent
with TeV halo observations: (1) the diffusion coefficient is suppressed by 2–3
orders of magnitude, (2) the diffusion coefficient reaches a minimum value at
$\sim$100 kyr after pulsar formation. Moreover, we note that our 1D framework
is a reasonable approximation of reality, as propagation through 1D flux tubes
likely dominates diffusion on distances less than the 1–100 pc coherence
length of the galactic magnetic field López-Coto and Giacinti (2018);
Iacobelli  _et al._ (2013); Haverkorn  _et al._ (2008). On the other hand,
these models may appear insconsistent with data, because observations indicate
that TeV halos are roughly spherical sources.

However, we _stress_ that our 1D diffusion models are very different (and more
consistent with data) than competing flux-tube models (e.g., Liu  _et al._
(2019)) that do not use cosmic-ray self-confinement. In previous models, the
flux tube must be oriented directly towards Earth (with an offset less than
$\sim$5∘), to ensure that the halo does not extend too far in the transverse
(observable) direction. Such models are difficult to rectify with observations
which indicate that many pulsars produce observable TeV halos.

In our 1D models, the preferential direction for particle diffusion is
irrelevant. The cosmic-ray self-confinement mechanism efficiently inhibits
cosmic-ray propagation along this direction until it is _also_ inhibited.
Notably, the degree of inhibited diffusion (by a factor of $\sim$100) is
similar to the predicted efficiency of diffusion in directions perpendicular
standard 1D flux tubes Liu  _et al._ (2019). The sum of these effects produce
a cosmic-ray population whose diffusion is inhibited in every direction.
Therefore, it is _critical_ to note that the $\gamma$-ray emission morphology
from our 1D models is not yet known. It could potentially even be spherical,
if cosmic-ray self-confinement forces the diffusion coefficients in every
direction towards similar values. The interaction of cosmic-ray self-
confinement mechanisms with a multi-dimensional diffusion tensor, and the
comparison of such a model with observed TeV halos, is important work which we
leave for future publications.

Finally, while 1D models predict that diffusion continues to be inhibited on
distance scales that are much larger than the observed size of TeV halos, this
result is likely consistent with data for two reasons. First, the edges of TeV
halos are not well-defined and may be significantly affected by the flux
sensitivity of each $\gamma$-ray instrument. Indeed subsequent HAWC
observations of Geminga and Monogem found significantly larger TeV halos,
which was consistent with the improved exposure of each source Abeysekara  _et
al._ (2017). Second, we expect the flux tube approximation to break down on
larger distance scales, an effect which will quickly terminate the region
where diffusion is inhibited, as demonstrated from our 3D models.

In 3D models, on the other hand, the pulsar is not energetic enough to
significantly inhibit cosmic-ray diffusion – even when the model is strongly
tuned to maximize the inhibition of diffusion near 10 TeV. To produce
reasonable models, we added a contribution from the coincident SNR, which adds
an additional 4$\times$1050 ergs of cosmic-ray energy into the local ISM (for
our best-fit $\xi_{CR}=0.4$ value). In this case, we find that diffusion is
significantly inhibited on $\sim$5 pc scales, reaching the Bohm limit.
However, the large volume of 3D simulations still prevents our model from
significantly affecting cosmic-ray diffusion on 20 pc scales.

We note, however, that our 3D model is _maximally_ pessimistic for cosmic-ray
self-confinement models. The simulation begins with a diffusion coefficient
that is unsuppressed in every direction around the central source, allowing
cosmic-rays to quickly be diluted as they move outword in every direction.
Models that include inhibited diffusion in one direction (even on small
distance scales), will force the cosmic-ray flux to build up, producing larger
gradients (and thus significantly inhibited diffusion) throughout the
simulation volume. Notably, these effects can be strongly non-linear, as shown
by the significant changes in the diffusion coefficient which result from
relatively small changes in the pulsar efficiency. Finally, we note that these
models do not include SNR shocks or other molecular dynamics that may also
significantly increase magnetic turbulence and inhibit cosmic-ray diffusion.

To conclude, our 3D models where the self-confinement scenario is ineffective
at inhibiting cosmic-ray diffusion on 10–20 pc scales are models where the
initial diffusion coefficient is given by the standard ISM value in every
spatial direction surrounding the pulsar source. However, any adjustments to
this model, including diffusion coefficients that are initially suppressed in
one or more spatial directions, or contributions from SNR/PWN activity that
lower the initial diffusion coefficient near the pulsar source – will
unambiguously act to strengthen the self-confinement mechanism by allowing the
particle density to build up in compact regions. Thus, we conclude that in
realistic models, the cosmic-ray self-confinement mechanism is capable of
driving the inhibited diffusion coefficients observed out to 10–20 pc from
pulsar sources.

###  VI.2 Future Improvements

Finally, we note a number of potential improvements which lie beyond the scope
of this model, but which may be considered in upcoming work. First, we note
that some pulsars (including Geminga), have large proper velocities which have
significantly removed them from their natal SNR by the $\sim$300 kyr
observation time of their TeV halos. Such models make the inclusion of
turbulence from the initial SNR less realistic – though the SNR may stil act
as important seed to drive down the initial diffusion coefficient. Future
models must consider the impact of pulsar proper motion through the ISM.
Secondly, in our models, we have not yet considered the potential interplay
between 1D and 3D diffusion on the boundary between different flux tubes.
However, our results (which indicate that pulsars produce diffusion that is
too inhibited in 1D simulations but not inhibited enough in 3D simulations),
strongly motivate such an extension.

## Acknowledgements

We would like to thank Felix Aharonian, Elena Amato, John Beacom, Pasquale
Blasi, Ilias Cholis, Dan Hooper, Alexandre Marcowith, Pierrick Martin, Sarah
Recchia and Takahiro Sudoh for helpful comments which greatly improved the
quality of this manuscript. We would especially like to thank Carmelo Evoli
and Giovanni Morlino for a number of comments, code comparison tests, and
discussions that motivated and significantly improved the quality of this
manuscript. TL is partially supported by the Swedish Research Council under
contract 2019-05135, the Swedish National Space Agency under contract 117/19
and the European Research Council under grant 742104. This project used
computing resources from the Swedish National Infrastructure for Computing
(SNIC) under project Nos. 2020/5-463 and 2021-1-24 partially funded by the
Swedish Research Council through grant no. 2018-05973.

## References

  * Abeysekara  _et al._ (2017) A. U. Abeysekara _et al._ , Astrophys. J. 843, 40 (2017), arXiv:1702.02992 [astro-ph.HE] .
  * Abeysekara  _et al._ (2017) A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, J. D. Álvarez, R. Arceo, J. C. Arteaga-Velázquez, D. Avila Rojas, H. A. Ayala Solares, A. S. Barber, N. Bautista-Elivar, A. Becerril, E. Belmont-Moreno, S. Y. BenZvi, and D. t. Berley, Science 358, 911 (2017), arXiv:1711.06223 [astro-ph.HE] .
  * Linden  _et al._ (2017) T. Linden, K. Auchettl, J. Bramante, I. Cholis, K. Fang, D. Hooper, T. Karwal, and S. W. Li, Phys. Rev. D 96, 103016 (2017), arXiv:1703.09704 [astro-ph.HE] .
  * Riviere  _et al._ (2017) C. Riviere, H. Fleischhack, and A. Sandoval, The Astronomer’s Telegram 10941, 1 (2017).
  * Brisbois  _et al._ (2018) C. Brisbois, C. Riviere, H. Fleischhack, and A. Smith, The Astronomer’s Telegram 12013, 1 (2018).
  * Fang  _et al._ (2019) K. Fang, X.-J. Bi, and P.-F. Yin, MNRAS 488, 4074 (2019), arXiv:1903.06421 [astro-ph.HE] .
  * Sudoh  _et al._ (2019) T. Sudoh, T. Linden, and J. F. Beacom, Phys. Rev. D 100, 043016 (2019), arXiv:1902.08203 [astro-ph.HE] .
  * Fang  _et al._ (2021) K. Fang, S.-Q. Xi, and X.-J. Bi, (2021), arXiv:2107.02140 [astro-ph.HE] .
  * Aharonian  _et al._ (2021) F. Aharonian _et al._ (LHAASO), Phys. Rev. Lett. 126, 241103 (2021), arXiv:2106.09396 [astro-ph.HE] .
  * Linden and Buckman (2018) T. Linden and B. J. Buckman, Phys. Rev. Lett. 120, 121101 (2018), arXiv:1707.01905 [astro-ph.HE] .
  * Profumo  _et al._ (2018) S. Profumo, J. Reynoso-Cordova, N. Kaaz, and M. Silverman, Phys. Rev. D 97, 123008 (2018), arXiv:1803.09731 [astro-ph.HE] .
  * Giacinti  _et al._ (2020) G. Giacinti, A. M. W. Mitchell, R. López-Coto, V. Joshi, R. D. Parsons, and J. A. Hinton, Astron. Astrophys. 636, A113 (2020), arXiv:1907.12121 [astro-ph.HE] .
  * Aharonian (1995) F. A. Aharonian, Nucl. Phys. B Proc. Suppl. 39, 193 (1995).
  * Aharonian  _et al._ (1995) F. A. Aharonian, A. M. Atoyan, and H. J. Volk, Astron. Astrophys. 294, L41 (1995).
  * Abdo  _et al._ (2007) A. A. Abdo _et al._ , Astrophys. J. Lett. 664, L91 (2007), arXiv:0705.0707 [astro-ph] .
  * Abdo  _et al._ (2009) A. A. Abdo, B. T. Allen, T. Aune, D. Berley, C. Chen, G. E. Christopher, T. DeYoung, B. L. Dingus, R. W. Ellsworth, M. M. Gonzalez, J. A. Goodman, E. Hays, C. M. Hoffman, P. H. Hüntemeyer, B. E. Kolterman, J. T. Linnemann, J. E. McEnery, T. Morgan, A. I. Mincer, P. Nemethy, J. Pretz, J. M. Ryan, P. M. Saz Parkinson, A. Shoup, G. Sinnis, A. J. Smith, V. Vasileiou, G. P. Walker, D. A. Williams, and G. B. Yodh, ApJ 700, L127 (2009), arXiv:0904.1018 [astro-ph.HE] .
  * Yuksel  _et al._ (2009) H. Yuksel, M. D. Kistler, and T. Stanev, Phys. Rev. Lett. 103, 051101 (2009), arXiv:0810.2784 [astro-ph] .
  * Hooper  _et al._ (2017) D. Hooper, I. Cholis, T. Linden, and K. Fang, Phys. Rev. D 96, 103013 (2017), arXiv:1702.08436 [astro-ph.HE] .
  * Fang  _et al._ (2018) K. Fang, X.-J. Bi, P.-F. Yin, and Q. Yuan, Astrophys. J. 863, 30 (2018), arXiv:1803.02640 [astro-ph.HE] .
  * Di Mauro  _et al._ (2020) M. Di Mauro, S. Manconi, and F. Donato, Phys. Rev. D 101, 103035 (2020), arXiv:1908.03216 [astro-ph.HE] .
  * Di Mauro  _et al._ (2019) M. Di Mauro, S. Manconi, and F. Donato, Phys. Rev. D 100, 123015 (2019), arXiv:1903.05647 [astro-ph.HE] .
  * Tang and Piran (2019) X. Tang and T. Piran, Mon. Not. Roy. Astron. Soc. 484, 3491 (2019), arXiv:1808.02445 [astro-ph.HE] .
  * Giacinti and López-Coto (2019) G. Giacinti and R. López-Coto, J. Phys. Conf. Ser. 1181, 012046 (2019), arXiv:1811.00494 [astro-ph.HE] .
  * López-Coto and Giacinti (2018) R. López-Coto and G. Giacinti, Mon. Not. Roy. Astron. Soc. 479, 4526 (2018), arXiv:1712.04373 [astro-ph.HE] .
  * Liu  _et al._ (2019) R.-Y. Liu, H. Yan, and H. Zhang, Phys. Rev. Lett. 123, 221103 (2019), arXiv:1904.11536 [astro-ph.HE] .
  * Recchia  _et al._ (2021) S. Recchia, M. Di Mauro, F. A. Aharonian, F. Donato, S. Gabici, and S. Manconi, arXiv e-prints , arXiv:2106.02275 (2021), arXiv:2106.02275 [astro-ph.HE] .
  * Bao  _et al._ (2021) L.-Z. Bao, K. Fang, and X.-J. Bi, (2021), arXiv:2107.07395 [astro-ph.HE] .
  * Evoli  _et al._ (2018a) C. Evoli, T. Linden, and G. Morlino, Phys. Rev. D 98, 063017 (2018a), arXiv:1807.09263 [astro-ph.HE] .
  * Fang  _et al._ (2019) K. Fang, X.-J. Bi, and P.-F. Yin, Mon. Not. Roy. Astron. Soc. 488, 4074 (2019), arXiv:1903.06421 [astro-ph.HE] .
  * D’Angelo  _et al._ (2018) M. D’Angelo, G. Morlino, E. Amato, and P. Blasi, Mon. Not. Roy. Astron. Soc. 474, 1944 (2018), arXiv:1710.10937 [astro-ph.HE] .
  * Evoli  _et al._ (2018b) C. Evoli, P. Blasi, G. Morlino, and R. Aloisio, Phys. Rev. Lett. 121, 021102 (2018b), arXiv:1806.04153 [astro-ph.HE] .
  * Skilling (1970) J. Skilling, MNRAS 147, 1 (1970).
  * Lagage and Cesarsky (1983) P. O. Lagage and C. J. Cesarsky, A&A 118, 223 (1983).
  * Aloisio and Blasi (2013) R. Aloisio and P. Blasi, J. Cosmology Astropart. Phys 2013, 001 (2013), arXiv:1306.2018 [astro-ph.HE] .
  * Johannesson  _et al._ (2019) G. Johannesson, T. A. Porter, and I. V. Moskalenko, Astrophys. J. 879, 91 (2019), arXiv:1903.05509 [astro-ph.HE] .
  * Evoli  _et al._ (2021) C. Evoli, E. Amato, P. Blasi, and R. Aloisio, Phys. Rev. D 103, 083010 (2021), arXiv:2010.11955 [astro-ph.HE] .
  * Manchester  _et al._ (2005) R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, AJ 129, 1993 (2005), arXiv:astro-ph/0412641 [astro-ph] .
  * Iacobelli  _et al._ (2013) M. Iacobelli _et al._ , Astron. Astrophys. 558, A72 (2013), arXiv:1308.2804 [astro-ph.GA] .
  * Haverkorn  _et al._ (2008) M. Haverkorn, J. C. Brown, B. M. Gaensler, and N. M. McClure-Griffiths, Astrophys. J. 680, 362 (2008), arXiv:0802.2740 [astro-ph] .
  * Blumenthal and Gould (1970) G. R. Blumenthal and R. J. Gould, Rev. Mod. Phys. 42, 237 (1970).
  * Schlickeiser and Ruppel (2010) R. Schlickeiser and J. Ruppel, New J. Phys. 12, 033044 (2010), arXiv:0908.2183 [astro-ph.HE] .
  * Fang  _et al._ (2021) K. Fang, X.-J. Bi, S.-J. Lin, and Q. Yuan, Chinese Physics Letters 38, 039801 (2021), arXiv:2007.15601 [astro-ph.HE] .
  * Di Mauro  _et al._ (2021) M. Di Mauro, F. Donato, and S. Manconi, Phys. Rev. D 104, 083012 (2021), arXiv:2010.13825 [astro-ph.HE] .
  * Skilling (1971) J. Skilling, ApJ 170, 265 (1971).
  * Ptuskin and Zirakashvili (2003) V. S. Ptuskin and V. N. Zirakashvili, A&A 403, 1 (2003), arXiv:astro-ph/0302053 [astro-ph] .
  * Brunetti  _et al._ (2004) G. Brunetti, P. Blasi, R. Cassano, and S. Gabici, MNRAS 350, 1174 (2004), arXiv:astro-ph/0312482 [astro-ph] .
  * Piddington (1956) J. H. Piddington, MNRAS 116, 314 (1956).
  * Kulsrud and Pearce (1969) R. Kulsrud and W. P. Pearce, ApJ 156, 445 (1969).
  * Zweibel and Shull (1982) E. G. Zweibel and J. M. Shull, ApJ 259, 859 (1982).
  * Caraveo  _et al._ (2003) P. A. Caraveo, G. F. Bignami, A. De Luca, S. Mereghetti, A. Pellizzoni, R. Mignani, A. Tur, and W. Becker, Science 301, 1345 (2003).
  * Kulsrud (1978) R. M. Kulsrud, in  _Astronomical Papers Dedicated to Bengt Stromgren_ , edited by A. Reiz and T. Andersen (1978) pp. 317–326.
  * Wiener  _et al._ (2013) J. Wiener, E. G. Zweibel, and S. P. Oh, ApJ 767, 87 (2013), arXiv:1301.4445 [astro-ph.GA] .
  * Blasi (2019) P. Blasi, Galaxies 7, 64 (2019), arXiv:1905.11149 [astro-ph.HE] .
  * Farmer and Goldreich (2004) A. J. Farmer and P. Goldreich, Astrophys. J. 604, 671 (2004), arXiv:astro-ph/0311400 .
  * Liu  _et al._ (2019) R.-Y. Liu, C. Ge, X.-N. Sun, and X.-Y. Wang, ApJ 875, 149 (2019), arXiv:1904.11438 [astro-ph.HE] .
  * Ptuskin  _et al._ (2008) V. S. Ptuskin, V. N. Zirakashvili, and A. A. Plesser, Advances in Space Research 42, 486 (2008).