# Gamma-ray Pulsar Halos in the Galaxy

Kun Fang 1,∗

###### Abstract

Pulsar halos are extended gamma-ray structures generated by electrons and
positrons escaping from pulsar wind nebulae (PWNe), considered a new class of
gamma-ray sources. They are ideal indicators for cosmic-ray propagation in
localized regions of the Galaxy and particle escape process from PWNe. The
cosmic-ray diffusion coefficient inferred from pulsar halos is more than two
orders of magnitude smaller than the average value in the Galaxy, which has
been arousing extensive discussion. We review the recent advances in the study
of pulsar halos, including the characteristics of this class of sources, the
known pulsar halos, the possible mechanisms of the extremely slow diffusion,
the critical roles of pulsar halos in the studies of cosmic-ray propagation
and electron injection from PWNe, and the implications on the problems of the
cosmic positron excess and the diffuse TeV gamma-ray excess. Finally, we give
prospects for the study in this direction based on the expectation of a larger
sample of pulsar halos and deeper observations for bright sources.

##  1 Keywords:

cosmic rays, cosmic ray propagation, supernova remnants, pulsar, pulsar wind
nebulae

##  2 Introduction

In 2007, the Milagro Gamma-ray Observatory reported an extended TeV gamma-ray
source around the Geminga pulsar (Abdo et al., 2007). At that time, the source
was considered a TeV pulsar wind nebula (PWN) associated with Geminga (Yuksel
et al., 2009). However, the $\sim$$$ extension of the source is hard to
understand, which is dozens of times larger than the bow-shock PWN of Geminga
observed at X-ray wavelengths (Caraveo et al., 2003). A decade later, the
High-Altitude Water Cherenkov Observatory (HAWC) updated the measurement of
this gamma-ray source with better sensitivity and angular resolution
(Abeysekara et al., 2017). The gamma-ray morphology shows radial symmetry and
can be well interpreted by the diffusion model, indicating that the source is
more likely to be generated by electrons and positrons111Electrons will denote
both electrons and positrons hereafter if not specified. that have escaped
from the central PWN and are freely diffusing in the interstellar medium (ISM)
(Hooper et al., 2017; Fang et al., 2018; Profumo et al., 2018; Tang and Piran,
2019; Giacinti et al., 2020). This new source type is known as the gamma-ray
pulsar halo, TeV halo, or ICS halo. We will use the name pulsar halo
throughout this paper.

Pulsar halos are ideal probes for cosmic-ray (CR) propagation in localized
regions of the Galaxy. High-energy electrons released by the PWNe inverse
Compton (IC) scatter homogeneous background photons to produce the gamma-ray
halos, so the gamma-ray morphologies unambiguously trace the propagation of
the parent electrons. The most intriguing result is that the inferred electron
diffusion coefficient is several hundred times smaller than the average CR
diffusion coefficient in the Galaxy (Abeysekara et al., 2017; Aharonian et
al., 2021). It complicates the image of Galactic CR propagation and has
significant impacts on some key issues of CRs, such as the origin of the
positron excess (Adriani et al., 2009; Ackermann et al., 2012b; Aguilar et
al., 2013) and the diffuse TeV gamma-ray excess (Atkins et al., 2005; Abdo et
al., 2008; Amenomori et al., 2021).

Pulsar halos are essential for the study of electron injection from PWNe (Fang
and Bi, 2022; Fang et al., 2022). As the gamma-ray halos are generated by the
escaping electrons from PWNe, the energy spectrum of the electrons injected
from PWNe into the ISM (injection spectrum for short) should be inferred from
the halo spectrum rather than the spectrum of PWNe themselves. The injection
spectrum could be further used to estimate the acceleration limit of PWNe
(Fleischhack et al., 2019). Unlike the high-energy spectral cutoff of X-ray
PWNe that is determined by both the maximum electron energy and the magnetic
field of PWNe, the gamma-ray spectrum of pulsar halos can unambiguously
indicate the electron spectral cutoff. Pulsar halos are also proposed to
search for invisible pulsars (Linden et al., 2017).

This work reviews the recent advances in the study of pulsar halos. In Section
3, we introduce the characteristics of this class of sources and the known
pulsar halos. In Section 4, current theories about the origin of pulsar halos
are reviewed. We present the significance of pulsar halos to the studies of CR
propagation and electron injection from PWNe in Sections 5 and 6,
respectively, before presenting the implications for interpreting the positron
excess in Section 7 and the diffuse TeV gamma-ray excess in Section 8. We look
forward to the future studies of pulsar halos in Section 9 and give the
summary points in Section 10. One may also refer to other recent reviews of
pulsar halos (López-Coto et al., 2022; Liu, 2022), which have different
focuses from the present work.

##  3 Characteristics of Pulsar halos

###  3.1 Criteria of a Pulsar Halo

We first present several criteria for a pulsar halo and then introduce the
characteristics of pulsar halos by explaining these criteria. Some of them
have been proposed in Fang et al. (2022).

  1. 1.

The TeV pulsar halo of a visible pulsar should have spatial coincidence with
that pulsar.

  2. 2.

The spin-down luminosity of the central pulsar should be large enough to
generate the pulsar halo.

  3. 3.

The extension of a gamma-ray pulsar halo should be significantly larger than
that of the X-ray PWN if the X-ray PWN is observable.

  4. 4.

The morphology of a pulsar halo should be interpreted well by the diffusion-
loss propagation of electrons with a reasonable diffusion coefficient.

After escaping from the PWN, the electron propagation in ISM is generally
described by the diffusion process, so the gamma-ray halo should be centered
at the host pulsar, which is a more specific description of criterion 1 (CT
1). Even if the diffusion environment is asymmetric, the flux peak of the
gamma-ray halo should be coincident with the pulsar. CT 1 holds for
$E_{\gamma}\gtrsim 1$ TeV (Di Mauro et al., 2020), where $E_{\gamma}$ is the
gamma-ray energy. For example, HESS J1026$-$582 is a TeV gamma-ray source once
considered a possible pulsar halo (Di Mauro et al., 2020). However, the offset
between the source centroid and the assumed host pulsar PSR J1028$-$5819 is
significantly larger than the Gaussian extension of the source (Abramowski et
al., 2011), which cannot pass CT 1. On the other hand, a GeV pulsar halo may
deviate from the position of the associated pulsar owing to the pulsar’s
proper motion (Johannesson et al., 2019; Di Mauro et al., 2019; Zhang et al.,
2021).

CT 2 is essential to judge the association between a gamma-ray source and a
nearby pulsar. For example, HESS J1632$-$478 was once connected with a nearby
middle-aged pulsar PSR J1632$-$4757 (Di Mauro et al., 2020). However, the
gamma-ray luminosity of HESS J1632$-$478 is significantly larger than the
spin-down energy of the pulsar, and this source is more likely to be a PWN
associated with a much younger pulsar (Balbo et al., 2010). CT 2 is also
helpful for constraining theoretical models of pulsar halos. A model will be
disfavored if the required injection energy is larger than the pulsar can
provide. An example is shown in Section 4.4.

According to the standard picture of PWN evolution (Gaensler and Slane, 2006),
a pulsar and its original PWN are located at the center of the parent
supernova remnant (SNR) in the early age. The expanding PWN encounters the SNR
reverse shock and is disrupted after a time of $\sim 7$ kyr, leaving a relic
PWN. Meanwhile, the pulsar leaves the birthplace and then the relic PWN owing
to the kick velocity (typically $400-500$ km s-1) and finally escapes into the
ISM after a time of $\sim$50 kyr. Only then can the gamma-ray halo morphology
indicates the CR propagation in the ISM. Thus, we may give a preference to the
pulsars older than $\sim$50 kyr for the study of pulsar halos. For a gamma-ray
source associated with a much younger pulsar, we should be more cautious as it
is likely to be a pure PWN or a source in the mixed state of a pulsar halo and
a relic PWN. Older pulsars with $\tau\gtrsim 1$ Myr may not hold an observable
pulsar halo owing to the low spin-down energy, so pulsar halos are generally
considered to be generated by middle-aged pulsars with $\tau\sim 100$ kyr.
However, a recycled millisecond pulsar with a very large characteristic age
($\tau\gtrsim 50$ Myr) can be as bright as middle-aged pulsars and may also
hold a pulsar halo (Hooper and Linden, 2018, 2022).

When a pulsar travels outside its SNR (or in the outer part of its SNR), the
pulsar motion can be supersonic and drive a bow-shock PWN that is
significantly different from its original PWN. The size of the bow-shock PWN
is limited by the ram pressure caused by the pulsar motion, which is $\lesssim
1$ pc and not expected to expand with time (Gaensler and Slane, 2006). The
bow-shock PWN in X-ray and the corresponding TeV halo are generated by
electrons with similar energies. Thus, the X-ray PWN can be regarded as the
source region of the TeV halo and should then be significantly smaller than
the TeV halo, resulting in CT 3. On the other hand, a TeV source that has a
comparable size to the X-ray PWN is more likely a gamma-ray PWN (Abramowski et
al., 2012). CT 3 is important in distinguishing a pulsar halo from a PWN.

It is also necessary to distinguish between a pulsar halo and a gamma-ray
source with a hadronic origin. A nearby hadronic accelerator and a reasonable
gas distribution in the source region are the essential conditions of a
hadronic origin. Moreover, these two scenarios predict different gamma-ray
morphologies (Yang and Liu, 2022). The morphology of a pulsar halo should be
interpreted well by the diffusion-loss propagation of electrons, which is CT
4. The derived diffusion coefficient may also help to constrain the origin of
the source. Bohm limit is considered the lower limit of the diffusion
coefficient and is written as $1/3r_{g}c$, where $r_{g}$ is the gyro-radius of
CR particles. If the derived diffusion coefficient is smaller than the Bohm
limit under the typical magnetic field strength in the ISM, the gamma-ray
structure may not be a pulsar halo, or the pulsar distance is significantly
underestimated. Besides, CT 4 could be valuable in finding pulsar halos of
invisible pulsars.

The standard diffusion scenario fits well with the current morphology
measurements of pulsar halos, and CT 4 may play a significant role in
identifying pulsar halos at present. However, we cannot rule out the
possibility that more sophisticated propagation models are needed to interpret
more accurate observations in the future.

###  3.2 Observed Pulsar Halos

At present, about ten gamma-ray sources are proposed to be pulsar halos or
candidates for pulsar halos (Albert et al., 2020; Aharonian et al., 2021).
Four of these sources have been well studied as pulsar halos: the Geminga
halo, the Monogem halo, LHAASO J0621$+$3755, and HESS J1831$-$0952\. Figure 1
summarizes the significance maps of these sources. Information on the
associated pulsars and two important parameters derived from the observations
of the pulsar halos are summarized in Table 1. The Geminga halo is the
brightest pulsar halo and is regarded as the canonical source of this class.
Geminga is one of the brightest gamma-ray pulsars (Abdo et al., 2013) and is
close to Earth. It has received much attention from the CR community as it is
a likely source of the cosmic positron excess. We continue this topic in
Section 7.

Figure 1: Significance maps of the known gamma-ray pulsar halos. Left: Geminga
and Monogem halos observed by HAWC, taken from Abeysekara et al. (2017) with
permission from Science. ©2017 AAAS. Middle: LHAASO J0621$+$3755 observed by
LHAASO-KM2A, taken from Aharonian et al. (2021) with permission from Physical
Review Letters. ©2021 APS. Right: Significance map of HESS J1831$-$098 given
by the H.E.S.S. Galactic plane survey (Abdalla et al., 2018) with the
integration radius of $None$, taken from Fang et al. (2022). One may refer to
the original papers for the meaning of the marks in the figures. Table 1:
Pulsar parameters of known pulsar halos and parameters inferred from pulsar
halos.

Name | Parameters of associated pulsars | Derived parameters
---|---|---
Gl | Gb | $T$ | $D$ | $L$ | $D_{100}$ | $\eta$
[deg] | [deg] | [kyr] | [kpc] |  $10^{34}$ [erg s-1] | [cm2 s-1] | [%]
Geminga halo | $195.13$ | $4.27$ | $342$ | $0.25^{a}$ | $3.2$ |  $4.6\times 10^{27}$ c | $5^{c}$
Monogem halo | $201.11$ | $8.26$ | $111$ | $0.288$ | $3.8$ |  $1.5\times 10^{28}$ d | $4^{d}$
LHAASO J0621$+$3755 | $175.88$ | $10.96$ | $208$ | $1.6^{b}$ | $2.7$ |  $2.3\times 10^{27}$ c | $40^{c}$
HESS J1831$-$098 | $21.90$ | $-0.13$ | $128$ | $3.68$ | $110$ |  $9.0\times 10^{27}$ e | $7^{e}$

  * •

Notes. Columns 2-6 are the Galactic longitude, Galactic latitude,
characteristic age, distance, and spin-down luminosity, respectively. Pulsar
parameters are taken from the ATNF pulsar catalog (Manchester et al., 2005)
unless specified. Column 7 is the diffusion coefficient normalized at 100 TeV.
Column 8 is the needed conversion efficiency from the pulsar spin-down energy
to the electron energy of pulsar halo. $\eta$ is obtained by assuming an
injection spectrum with a ECPL form, except the case of Monogem, where $\eta$
is obtained by assuming a simple power-law injection. $\eta$ will be
significantly smaller if the ECPL injection is adopted for Monogem.

  * •

References: aFaherty et al. (2007); bParkinson et al. (2010); cBao et al.
(2022); dAbeysekara et al. (2017); eFang et al. (2022).

The Monogem halo is associated with PSR B0656$+$14, known as the Monogem
pulsar. The Monogem halo was reported by HAWC together with the Geminga halo
(Abeysekara et al., 2017) and is also a nearby source. However, it is
significantly dimmer than the Geminga halo, so the diffusion coefficient in
the halo is not accurately constrained.

LHAASO J0621$+$3755 is the pulsar halo of PSR J0622$+$3749 and is the first
pulsar halo reported by the Large High-Altitude Air Shower Observatory
(LHAASO) (Aharonian et al., 2021). Intriguingly, this source could be far away
from the Galactic plane ($\approx 300$ pc above the plane), which is distinct
from the other known pulsar halos.

HESS J1831$-$098 is the pulsar halo of PSR J1831$-$0952\. This source was
first discovered by the High Energy Spectroscopic System (H.E.S.S) and was
considered an old PWN at that time (Sheidaei et al., 2011). However, Fang et
al. (2022) point out that it is more reasonable to regard it as a pulsar halo
as it passes all the criteria introduced above. Although this source is more
than $10$ times farther than Geminga, the significantly stronger injection
power ensures its visibility.

The four observed pulsar halos all meet CT 1-4. The separations between the
centroids of the gamma-ray halos and the pulsars are all significantly smaller
than the extensions of the corresponding halos (Albert et al., 2020; Aharonian
et al., 2021; Fang et al., 2022), indicating good spatial coincidence.
Criteria about conversion efficiency and pulsar age are all satisfied, as
shown in Table 1. Geminga, Monogem, and PSR J1831$-$0952 all hold X-ray PWNe,
which are at least dozens of times smaller than the corresponding TeV halos
(Posselt et al., 2017; Bîrzan et al., 2016; Abichandani et al., 2019). No
clear x-ray PWN is found for PSR J0622$+$3749, possibly owing to the faintness
of the PWN (Aharonian et al., 2021). Finally, all the TeV pulsar halos can be
interpreted with the standard diffusion framework under the current
measurement accuracy. The measured diffusion coefficients are summarized in
Table 1.

For other possible pulsar halos, such as HAWC J0543+233 (Riviere et al., 2017)
and HAWC J0635+070 (Brisbois et al., 2018), the information reported is still
not comprehensive enough, and we cannot test these candidates with all our
criteria at present.

##  4 Possible Origins of Pulsar halos

A pulsar halo is not expected under the average CR diffusion coefficient in
the Galaxy as the electrons escaping from the PWN will spread out rapidly.
Hooper et al. (2017) pointed out that the enhanced flux and small extension of
the gamma-ray source around Geminga are likely to be interpreted by a slow-
diffusion environment. This view was supported by the subsequent morphology
measurement of HAWC (Abeysekara et al., 2017). As the gamma-ray profile of the
Geminga halo can be fitted well by the diffusion model with a small diffusion
coefficient, the slow-diffusion picture has been widely accepted. However, it
is challenging to interpret the origin of the slow-diffusion environment.

###  4.1 Self-generated MHD Turbulence

The most straightforward interpretation may be the self-generated scenario,
which has been adopted to predict the slow-diffusion environment around SNRs
(Ptuskin et al., 2008; Malkov et al., 2013; D’Angelo et al., 2016; Recchia et
al., 2022) and spatially dependent diffusion of CRs at the Galactic scale
(Recchia et al., 2016; Evoli et al., 2018a). A large density gradient of CRs
near the sources can induce the growth of the small-scale magnetohydrodynamic
(MHD) turbulence of the background plasma, known as the resonant streaming
instability (RSI) (Skilling, 1971). CRs can then be trapped by the enhanced
MHD turbulence generated by themselves, corresponding to the suppression of
the diffusion coefficient. In the one-dimensional (1D) diffusion case,
numerical calculations show that the positron/electron pairs released by the
Geminga PWN may suppress the diffusion coefficient to the observed level
through RSI if the pulsar motion is ignored (Evoli et al., 2018b; Fang et al.,
2019a; Mukhopadhyay and Linden, 2022). However, the proper motion measurement
of Geminga indicates that the pulsar has been traveling about $None$ from its
birthplace (Faherty et al., 2007), significantly larger than the size of the
pulsar halo. It means that the positron/electron pairs released at the early
age of Geminga cannot contribute to the formation of the slow-diffusion zone.
Fang et al. (2019a) proved through a simple analytical derivation that if the
pulsar motion is taken into account, the diffusion coefficient cannot be
suppressed to the required level even under very optimistic conditions for the
RSI growth (see the left of Figure 2). Moreover, the 1D diffusion mode may
eventually evolve into the 3D mode as suggested by the symmetry of the Geminga
halo. Schroer et al. (2021) also pointed out that the 1D flux-tube assumption
will be broken by the transverse expansion caused by the transverse pressure
gradient. The RSI growth is further limited in the 3D mode (Mukhopadhyay and
Linden, 2022) as the injected electrons are diluted compared with the 1D case.

Another branch of the self-generated scenario is the turbulent growth through
the nonresonant streaming instability (NRSI) (Bell, 2004). The growth rate of
NRSI could be larger than the NRI (Schroer et al., 2021), while it is
proportional to the total CR current, which is zero for the case of the
positron/electron pairs. Although the possible asymmetry between positrons and
electrons during the acceleration or escape process may excite the NRSI
(Schroer et al., 2022), quantitative calculations are needed to test this
interpretation.

Figure 2: Left: Lower limit of the diffusion coefficient in the Geminga halo
as a function of the injection spectral index under the self-generated
scenario. The RSI is assumed to account for the turbulent growth. One-
dimensional diffusion is assumed. Energy losses of electrons and turbulent
damping are ignored, which are optimistic conditions for turbulent growth.
Electrons injected during the last third of the age of Geminga are all assumed
to contribute to the turbulent growth. Right: Wave-number spectrum of the MHD
turbulence and the corresponding diffusion coefficient for the Geminga halo,
assuming the turbulence is provided by the SNR shock. $B_{0}$ is the
background magnetic field strength, and $\eta$ is the conversion efficiency
from the SNR initial energy to the MHD turbulenct energy. The gray dotted line
is the average diffusion coefficient in the Galaxy. Both figures are taken
from Fang et al. (2019a).

###  4.2 External Sources of MHD Turbulence

Besides being excited by pulsars themselves, the slow-diffusion environment
could also be generated by external sources. Stellar feedback, such as SN
explosion or stellar wind, is believed to be the main source of turbulence in
the ISM (Ferriere, 2020). For pulsars, their parent SNRs could be the most
straightforward sources for the turbulent environment. In the SNR shock frame,
the upstream plasma loses part of its kinetic energy when streaming through
the shock front, and this part of the energy is transferred into turbulence
and thermal energy behind the shock (Bell, 1978). Thus, the downstream region
could be highly turbulent (Giacalone and Jokipii, 2007), although the
turbulence will gradually decrease with the SNR evolution. Assuming the MHD
turbulence is injected at 10 pc and cascades to smaller scales, Fang et al.
(2019a) pointed out that the diffusion coefficient downstream of the SNR shock
could be suppressed by more than two orders of magnitude if 1-10% of the
initial energy of SNR is converted to the MHD turbulence (see the right of
Figure 2). If the pulsars are still inside their associated SNRs, the observed
slow-diffusion environment could be understood.

For Geminga, the parent SNR has not been observed, possibly owing to old age.
The proper motion of Geminga indicates that it has been $\approx 70$ pc away
from the SNR center, so there are works considering that Geminga has already
left behind the SNR. However, owing to the rarefied circumstance of Geminga
(Caraveo et al., 2003), the scale of the SNR may reach $\sim 100$ pc (Leahy
and Williams, 2017; Fang et al., 2019a), which means that Geminga may still be
downstream of the SNR shock. This scenario is consistent with the high-
ionization degree environment of Geminga indicated by the measurement of the
H$\alpha$ luminosity (Caraveo et al., 2003). Furthermore, as the speed of
sound in the outer part of the SNR is significantly smaller than that of the
motion of Geminga, this scenario is consistent with the bow-shock nature of
the Geminga PWN, meeting CT 3 given in Section 3.1. For the other known pulsar
halos, Monogem is likely inside its parent SNR, namely the Monogem ring
(Plucinsky et al., 1996; Knies et al., 2018), which is still observable in
X-ray owing to its younger age. However, the SNRs of PSR J0622$+$3749 and PSR
J1831$+$0952 have not been detected. Searching for traces of the associated
SNRs may be an important subject for this interpretation.

On the other hand, a slow-diffusion environment could also exist upstream of
the SNR shock, as mentioned in Section 4.1. This scenario is supported by the
observations of the interaction between the escaping CR nuclei from SNRs and
the molecular clouds near the SNRs (Fujita et al., 2009; Li and Chen, 2012).
Thus, pulsar halos may be generated if the host pulsars are embedded in the
turbulent medium upstream of the SNR shocks. Unlike positron/electron pairs
discussed in Section 4.1, CR nuclei could effectively induce the NRSI, and the
slow-diffusion zone may be retained for a long time owing to the large growth
rate of NRSI (Schroer et al., 2021, 2022).

Fang et al. (2019a) also proposed another possible scenario: Geminga is inside
a stellar-wind bubble, the Gemini H$\alpha$ Ring (Kim et al., 2007; Knies et
al., 2018). The shocked wind may provide a fresh MHD turbulence environment
for Geminga. The Gemini H$\alpha$ Ring is likely generated by several OB-type
stars that may have similar distances with Geminga (Knies et al., 2018).
Further quantitative calculations are needed to test this scenario.

###  4.3 Anisotropic Diffusion

Some intriguing models could interpret the pulsar halos without assuming
anomalously slow diffusion. The typical coherent length of the magnetic field
in the ISM is $\sim 50-100$ pc, and CRs should mainly propagate along the mean
magnetic field lines in the case of Alfvénic turbulence (Goldreich and
Sridhar, 1995). Meanwhile, CRs diffuse perpendicularly to the mean magnetic
field with a much smaller diffusion coefficient of
$D_{\perp}=D_{\parallel}M_{A}^{4}$, where $D_{\parallel}$ is the parallel
diffusion coefficient assumed to be the average value in the Galaxy, and
$M_{A}$ is the Alfvénic Mach number (Yan and Lazarian, 2008). The CR
propagation, in this case, is anisotropic, which seems to be inconsistent with
the symmetry of the pulsar halos. However, Liu et al. (2019b) pointed out that
if the local interstellar magnetic field of the Geminga pulsar is aligned with
our line of sight towards the pulsar, the morphology of the pulsar halo could
be interpreted by the inefficient perpendicular diffusion assuming
$M_{A}=0.2$. If so, sources of strong MHD turbulence discussed in Section 4.1
and 4.2 are not required. Another advantage of this model is that it could
naturally interpret the weak X-ray synchrotron emission in the vicinity of
Geminga (Liu et al., 2019a).

Obviously, this model has stringent requirements for the inclination angle
$\phi_{\rm incl}$ between the line of sight and the direction of the magnetic
field tube the pulsar embedding. A larger $M_{A}$ can reduce asymmetry but
will boost the extension of the source at the same time. The calculation of De
La Torre Luque et al. (2022) indicates that to consistently reproduce the
symmetry and extension of the Geminga halo, the permitted region in the
$M_{A}-\phi_{\rm incl}$ plane is very small, which means that the possibility
of finding multiple Geminga-type halos is extremely low. This scenario may
predict more pulsar halos with significant asymmetry, which could be tested by
future experiments (Yan et al., 2022).

Figure 3: Gamma-ray profiles of the Geminga halo (left) and LHAASO
J0621$+$3755 (right) after the relativistic correction to the propagation
equation, for both the fast- and slow-diffusion scenarios. $D_{0}$ is the
diffusion coeffcient at 1 GeV, and $\eta$ is the conversion efficiency from
the pulsar spin-down energy to the electron injection energy. Both figures are
taken from Bao et al. (2022).

###  4.4 Relativistic Diffusion

The commonly used diffusion model for CR propagation has the superluminal
problem. This problem can be ignored if the typical diffusion scale
($\sim\sqrt{Dt}$, where $D$ is the diffusion coefficient) is significantly
smaller than the distance traveled at the speed of light ($ct$), namely $t\gg
D/c^{2}$. Otherwise, the diffusion model will predict a significantly broader
CR distribution than the real case, and a relativistic correction to the
solution is essential (Dunkel et al., 2007; Aloisio et al., 2009). Another way
to describe this issue is that if $ct$ is smaller than the mean free path of
CRs, $\lambda$, that defines the diffusion coefficient with $D=1/3\lambda c$,
the CR propagation should be ballistic rather than diffusive. The critical
time is also $t\sim D/c^{2}$. It can be seen that the effect is more
significant for a larger $D$ as the distribution of CRs injected in a longer
period needs to be corrected in this case. Taking this effect into account,
the calculation of Recchia et al. (2021) showed that a steep gamma-ray profile
around Geminga is also predicted assuming the average diffusion coefficient in
the Galaxy, and a slow-diffusion environment may no longer be required. This
scenario seems attractive as it could be more general than the geometric
effect introduced in Section 4.3.

However, two problems remain to be solved for the fast diffusion model: the
goodness of fit and the conversion efficiency. For the Geminga halo, the
minimal reduced $\chi^{2}$ of the profile fit is $\approx 2$, and the required
conversion efficiency is $\approx 150\%$. As mentioned in Section 3.1, a
reasonable model of pulsar halo should interpret the observations with a
conversion efficiency smaller than 100%. These problems are more severe for
LHAASO J0621$+$3755 as shown by the calculation of Bao et al. (2022): the
minimal reduced $\chi^{2}$ of the profile fitting is $\approx 4$,
corresponding to exclusion with a confidence level of $99.996\%$, and the
required conversion efficiency is $\approx 350\%$222For LHAASO J0621$+$3755,
Recchia et al. (2021) obtained different results from Bao et al. (2022) as
they took a 1D Gaussian function as the point-spread-function (PSF). However,
the PSF should be defined in the 2D plane.. By comparison, the slow-diffusion
model is hardly affected by the relativistic correction. The reduced
$\chi^{2}$ for both the sources are around $1$, and the conversion
efficiencies are significantly smaller than $100\%$. Thus, the slow-diffusion
model is obviously the more reasonable one (Bao et al., 2022). The profile
fits of the two pulsar halos are shown in Figure 3 for both the fast- and
slow-diffusion models.

##  5 Pulsar Halos as Indicators of Galactic CR Propagation

Pulsar halos could be the most potent indicator of CR propagation in the
localized medium. Molecular clouds illuminated by escaping CRs from SNRs can
also measure the diffusion coefficient in the near-source region (Fujita et
al., 2009; Li and Chen, 2012). However, this method is difficult to measure
the CR spatial distributions in the ISM and investigate detailed mechanisms of
CR propagation. In comparison, pulsar halos’ morphology can straightforwardly
indicate the propagation of CR electrons in the ISM. We will show below that
the gamma-ray spectra of pulsar halos are also valuable for constraining
propagation models. The measurements with pulsar halos are complementary to
the global probes of Galactic CR propagation.

###  5.1 Local Propagation

Although the diffusion model is quite simple, it successfully describes CR
propagation in the Galaxy (Ginzburg and Syrovatskii, 1964). The Brownian
motion could simulate particle transport in the turbulent interstellar
magnetic field on microscopical scales. However, the normal diffusion model is
applicable only when the inhomogeneity of the ISM is negligible for the scale
of interest. The inhomogeneity of the ISM can have two types of effects on
diffusion properties (Uchaikin and Zolotarev, 1999):

  1. 1.

The diffusion coefficient is spatially dependent, which can be considered as a
superposition of normal diffusion processes.

  2. 2.

The particle motion is different from the Brownian motion, and the shape of
the diffusion packet is no longer Gaussian, known as the anomalous diffusion.

The extremely small diffusion coefficient measured with the pulsar halos
cannot be representative in the Galaxy. To reconcile the slow diffusion around
pulsars and the typical diffusion coefficient of the Galaxy, the two-zone
diffusion model is proposed (Hooper et al., 2017; Fang et al., 2018), where
the diffusion coefficient is only suppressed within $r_{\star}$ around pulsars
and recovers to the typical value outside $r_{\star}$. For the mechanisms
discussed in Sections 4.1 and 4.2, the slow diffusion could indeed be a near-
source phenomenon. The self-excited mechanism may only suppress the diffusion
coefficient within $\sim 50$ pc around the pulsar (Evoli et al., 2018b), as
the particle gradient at farther distances is too low to induce turbulent
growth. If the slow-diffusion environment is produced by the parent SNR,
$r_{\star}$ could be comparable to the SNR size and may also be $\sim 50$ pc.
The two-zone diffusion belongs to the first scenario listed above.

Gamma-ray spectra of pulsar halos suggest the possibility of two-zone
diffusion. For LHAASO J0621$+$3755, the TeV spectrum measured by the 1.3 km2
array of LHAASO (LHAASO-KM2A) and the GeV upper limits (ULs) given by the
Fermi Large Area Telescope (Fermi-LAT) cannot be interpreted consistently with
the one-zone normal diffusion model, or an unreasonably hard injection
spectrum is needed (Aharonian et al., 2021). We continue the topic of the
injection spectrum in Section 6. Under the two-zone diffusion assumption,
electrons with lower energies have a longer lifetime and have more chance to
escape from the slow-diffusion zone, which can suppress the low-energy
spectrum below the Fermi-LAT ULs, as shown in the left of Figure 4. Thus, the
two-zone diffusion model can interpret the wide-band gamma-ray spectrum of
LHAASO J0621$+$3755 with a reasonable injection spectrum if $r_{\star}\lesssim
50$ pc (Fang et al., 2021), consistent with the theoretical expectations. For
the Geminga halo and HESS J1831$-$098, the two-zone diffusion model may also
help to interpret the low flux or non-detection in the GeV band (Fang and Bi,
2022; Fang et al., 2022).

Figure 4: Left: Interpretation of the wide-band gamma-ray spectrum of LHAASO
J0621$+$3755 with two-zone diffusion models (Fang et al., 2021). Different
sizes of the slow-diffusion zone are adopted. The Fermi-LAT ULs vary with the
models. Right: Fitting results of the Geminga halo’s gamma-ray profile
measured by HAWC with superdiffusion models (Wang et al., 2021), where
$\alpha$ is the superdiffusion index. Superdiffusion degenerates to the normal
diffusion when $\alpha=2$.

On the other hand, the inhomogeneity of the ISM may arise on all scales. In
this case, the CR propagation cannot be described only by the spatially
dependent diffusion. The ISM is more likely to be a fractal type, and the
normal diffusion can be generalized to superdiffusion (Lagutin et al., 2001;
Lagutin and Uchaikin, 2003; Uchaikin and Sibatov, 2012, 2017).
Microscopically, the particle motion is described by Lévy flight instead of
the Brownian motion, where occasional very long steps are permitted. The
widening of the diffusion packet with time is faster for Lévy flight than the
normal diffusion case; hence superdiffusion is named. The fractional Laplacian
operator, $\Delta^{\alpha/2}$, is used for the superdiffusion equation. The
index $\alpha$ represents the degree of fractality of the ISM, which is
defined in $(0,2]$. When $\alpha=2$, the propagation degenerates to the normal
diffusion. Superdiffusion has been observed in interplanetary space (Perri and
Zimbardo, 2009b, a), the ISM (Perri et al., 2016), and the cluster of galaxies
(Ragot and Kirk, 1997), indicating that it could be common in astrophysics.
This picture is also adopted in the global propagation of Galactic CRs to
interpret features of the CR energy spectra (Lagutin et al., 2001; Volkov and
Lagutin, 2013).

As superdiffusion yields different CR spatial distributions from the normal
diffusion, the morphology of pulsar halos can be used to test the
superdiffusion model. The Geminga halo is the best target among the known
pulsar halos. Owing to the large extension of the Geminga halo, the morphology
measurements are affected little by the resolution of the experiments. Wang et
al. (2021) fitted the profile of the Geminga halo measured by HAWC with
different $\alpha$ and found that a model with $\alpha\lesssim 1.3$ is
disfavored at 95% confidence level. The fitting results are shown in the right
of Figure 4. Compared with the normal diffusion model, profiles predicted by
superdiffusion models are steeper near the source while flatter away from the
source. Models with $\alpha$ significantly smaller than $2$ give poor fits to
the near-source data, while the measurement still allows superdiffusion with
$\alpha\lesssim 2$. Fang et al. (2021) pointed out that superdiffusion is also
helpful in interpreting the non-detection of the GeV emission of LHAASO
J0621$+$3755.

###  5.2 Global Propagation

The CR boron-to-carbon ratio (B/C) is the major probe of the global
propagation of Galactic CRs (Strong et al., 2007). The CR B is entirely
secondary, produced by the interactions between the CR C, N, and O and the
ISM. Thus, B/C is directly proportional to the thickness of the Galactic
diffusion halo and inversely proportional to the diffusion coefficient. Ratios
between unstable secondary to primary CR nuclei, such as 10Be/9Be, can further
help to disentangle the diffusion coefficient from the halo thickness. The
average diffusion coefficient in the Galaxy is then obtained.

However, the diffusion coefficient is very likely spatially dependent in the
Galaxy as the distribution of turbulent sources is not homogeneous. This
picture is supported by observations. The measured turbulent magnetic field
energy in the inner Galactic disk could be one order of magnitude larger than
that in the high-latitude region of the outer Galaxy (Han, 2017). Spatially
dependent propagation models can also interpret various features of CR
measurements, such as the spectral hardening of the nuclei spectra, the radial
distributions of the CR proton densities and spectral indices, the diffuse
gamma-ray spectrum in different sky regions, and the CR anisotropy
(Tomassetti, 2012; Feng et al., 2016; Guo and Yuan, 2018; Zhao et al., 2021).
Among these measurements, the CR anisotropy is sensitive to the diffusion
coefficient in the nearby ISM. Based on the assumption that the 10 TeV
spectral bump of the CR proton spectrum is due to a nearby CR source (An et
al., 2019; Yue et al., 2020; Yuan et al., 2021), Fang et al. (2020) argued
that a slow-diffusion region could exist between this nearby source and the
solar system to consistently interpret the proton spectral feature and the low
CR anisotropy. Besides, spatially dependent diffusion is predicted by
different theoretical models (Evoli and Yan, 2014; Evoli et al., 2018a).

The discovery of pulsar halos is direct evidence of the spatially dependent
diffusion in the Galaxy if the morphology of pulsar halos is indeed due to
slow diffusion. For the mechanisms introduced in Sections 4.1 and 4.2 (except
for the stellar-wind bubble case), we may use the explosion rate of supernovae
to estimate the birth rate of the slow-diffusion zone. For an average size of
$50$ pc and a persisting time of $10^{6}$ yr, the slow-diffusion regions may
occupy $\sim 3\%$ of the volume of the Galactic disk (Hooper et al., 2017).
The diffusion coefficient is inversely proportional to the wave-number
spectrum of the MHD turbulence (Skilling, 1971). Thus, if the wave-number
spectrum is 300 times stronger in the slow-diffusion zones, the weight-
averaged spectrum in the Galactic disk will be amplified by $\sim 10$ times
($300\times 3\%$). It corresponds to 10 times suppression of the diffusion
coefficient in the Galactic disk on average, which may be consistent with the
results implied by the observations mentioned above. However, the universality
of pulsar halos is still in doubt, and we will discuss this topic in Section
9.2.1.

##  6 Pulsar Halos as Indicators of Electron Injection from Pulsar Wind
Nebulae

There are abundant measurements for the electromagnetic radiation spectrum of
PWNe (Reynolds et al., 2017), and we can infer the electron energy spectrum in
PWNe from these measurements. However, as the electron injection spectrum from
PWNe may be different from the spectrum in PWNe, we cannot only rely on the
observations of PWNe themselves to study the electron escape from PWNe. As
pulsar halos are generated by escaping electrons, their spectra are ideal
indicators of electron injection.

The latest observations from the HAWC and H.E.S.S. experiments exhibit complex
features in the gamma-ray spectrum of the Geminga halo (Torres Escobedo et
al., 2021; Mitchell et al., 2021), which is likely attributed to the electron
injection process of the Geminga PWN (Fang and Bi, 2022). The HAWC spectrum
indicates a possible bump feature around 10 TeV. Meanwhile, the spectrum
unexpectedly climbs again below $\approx 3$ TeV. This low-energy feature can
also be clearly seen from the H.E.S.S. spectrum. These results cannot be
interpreted by a single power-law injection spectrum (Fang and Bi, 2022).
Noting that the X-ray PWN and the TeV gamma-ray halo are generated by
electrons with almost the same energy range, (Fang and Bi, 2022) propose a
two-population injection model based on the image and spectral measurements of
the Geminga X-ray PWN. One population is the freshly accelerated electrons
that escape from the PWN through rapid outflows, corresponding to the lateral
tails of the X-ray PWN. This population has an exponentially-cutoff power-law
(ECPL) injection spectrum and accounts for the high-energy bump structure of
the gamma-ray spectrum. The hard power-law term is suggested by the X-ray PWN
spectrum. The high-energy cutoff term corresponds to the acceleration limit of
the PWN, and the cutoff energy is $\approx 140$ TeV as determined by the data
fit. The other population is the electrons trapped longer in the PWN before
escaping, corresponding to the axial tail of the X-ray PWN. This population
has an additional spectral break in the injection spectrum owing to the
synchrotron energy loss of electrons inside the PWN, dominating the low-energy
gamma-ray spectrum. The two-population model can interpret the HAWC and
H.E.S.S. data well, as shown in the left of Figure 5.

Figure 5: Left: Best-fit gamma-ray spectra to the HAWC and H.E.S.S.
observations of the Geminga halo with the two-population injection model (Fang
and Bi, 2022). Pop A represents the freshly accelerated electrons that escape
from the PWN through rapid outflows, while Pop B represents the electrons
trapped longer in the PWN before escaping. The black solid (dotted) line is
the total spectrum within $None$ ($None$) around the pulsar, which fits the
HAWC (H.E.S.S.) spectrum. Right: Best-fit gamma-ray spectra to the H.E.S.S.
data of HESS J1831$-$098 with a single ECPL injection model (Fang et al.,
2022). The two sub-components of scattering the dust photons and CMB photons
are also separately shown. The Fermi-LAT ULs are shown for comparison but not
used in the fitting procedure.

The current gamma-ray spectrum measurements of the Monogem halo and LHAASO
J0621$+$3755 are mainly above $10$ TeV (Abeysekara et al., 2017; Aharonian et
al., 2021), where the spectra are probably dominated by the high-energy cutoff
term mentioned above. While for HESS J1831$-$098, the gamma-ray spectrum has
been measured by H.E.S.S. in the energy range of $\sim 0.3-30$ TeV (Sheidaei
et al., 2011), which may provide us a comprehensive understanding of its
electron injection spectrum. Fang et al. (2022) perform a combined fit to the
gamma-ray spectrum and morphology of HESS J1831$-$098, finding the electron
injection spectrum can be described by the ECPL form well. The best-fit
spectrum is shown in the right of Figure 5. The best-fit cutoff energy is $52$
TeV, while the best-fit power-law index is $0.88$, significantly smaller than
the typical value of PWNe ($1.5$). The Fermi-LAT flux ULs at lower energies
even requires a harder power-law term. Such a hard spectrum could not be the
general case as constrained by the CR positron spectrum (Bitter and Hooper,
2022), although simulations indicate that an electron spectral index of $1.0$
is possible for bow-shock PWNe (Bykov et al., 2017). The two-zone diffusion
model or the time-delayed slow-diffusion model may interpret the H.E.S.S. and
Fermi-LAT data with a milder power-law term of the injection spectrum (Fang et
al., 2022). These models can suppress the low-energy gamma-ray fluxes and
could be the more realistic scenarios considering the possible origins of the
slow-diffusion zone (Sections 4.1 and 4.2).

##  7 Positron excess

In early works, CR positrons were considered to be secondary products
generated by the collisions between CR nuclei and the interstellar material
(Moskalenko and Strong, 1998). However, with the measurement of the positron
energy spectrum breaking through $\sim 100$ GeV, secondary positrons alone can
no longer interpret the high-energy fluxes, known as the positron excess
(Adriani et al., 2009; Ackermann et al., 2012b; Aguilar et al., 2013). This
phenomenon has attracted a lot of attention as the high-energy positrons may
come from the annihilation or decay of dark matter (Bergstrom et al., 2008;
Cholis et al., 2009; Yin et al., 2009; Yuan et al., 2015). However, the
contribution of dark matter to the positron spectrum is strongly constrained
by the gamma-ray observation of Fermi-LAT on the dwarf galaxies and so on
(Huang et al., 2012; Ackermann et al., 2015; Lin et al., 2015; Liu et al.,
2017). The positron excess could also be originated from nearby pulsars
(Hooper et al., 2009; Yuksel et al., 2009; Yin et al., 2013; Della Torre et
al., 2015). Electron-positron pairs are produced in the magnetosphere of
pulsars and further accelerated in the PWNe.

###  7.1 Geminga: the Main Candidate

Geminga is a likely source of the positron excess owing to its proper age and
distance and the relatively large spin-down luminosity (Yin et al., 2013).
Moreover, the X-ray PWN and gamma-ray pulsar halo of Geminga indicate that the
Geminga PWN can indeed accelerate and release high-energy positrons.

If pulsar halos are interpreted by the slow-diffusion scenario discussed in
Sections 4.1 and 4.2, the positron flux produced by pulsars at Earth will be
significantly affected. Assuming the small diffusion coefficient around
Geminga can be extrapolated to the whole space between Geminga and the solar
system, the positrons produced by Geminga can hardly reach Earth, as shown by
the red dotted line in the left of Figure 6. If so, Geminga will be unlikely
to account for the positron excess (Abeysekara et al., 2017). However, as we
have mentioned in Section 5.1, the slow-zone diffusion zones around pulsars
may have a typical size of $\sim 50$ pc considering the possible mechanisms,
which is also supported by the observations of LHAASO J0621$+$3755 (Fang et
al., 2021). Under the two-zone diffusion model with $r_{\star}\sim 50$ pc,
Geminga can still interpret the positron excess well (see the left of Figure
6). The needed conversion efficiency is even more reasonable than the commonly
used one-zone fast diffusion model (Fang et al., 2018; Profumo et al., 2018;
Tang and Piran, 2019). Moreover, if superdiffusion happens, Geminga could
produce a much higher positron flux at Earth than the normal diffusion model,
owing to the nature of Lévy flight (Wang et al., 2021). In this case, Geminga
could significantly contribute to the positron excess even if the small
diffusion coefficient measured around Geminga is applied in the whole region
between Geminga and Earth.

The positron flux from Geminga is also determined by the positron injection
spectrum. As introduced in Section 6, the injection spectrum can be
constrained by the spectral measurements of the pulsar halo. Note that the
positron excess happens in the range of $\sim 10-500$ GeV. Positrons in this
energy range emit GeV photons through the IC scattering. Thus, only GeV
observations of the Geminga halo could effectively constrain Geminga’s
contribution to the positron excess. Xi et al. (2019) analyzed the 10-yr
Fermi-LAT data in the direction of Geminga and did not detect extended
emission under different input parameters. The flux ULs seriously limit the
positron injection power in the energy range of interest, which disfavors
Geminga as the dominant source of the positron excess. On the other hand, Di
Mauro et al. (2019) also analyzed the Fermi-LAT data and claimed the detection
of the GeV halo of Geminga. Compared with Xi et al. (2019), they adopted a
larger window for the analysis and took the proper motion of Geminga into
account. Based on this measurement, Di Mauro et al. (2019) claimed that
Geminga alone still cannot interpret the positron excess due to the mismatch
of the spectral shape. However, the calculation of Zhou et al. (2022)
indicates that Geminga is still likely to account for the positron excess
under the constraint of the same gamma-ray observation. The difference may
result from different parameter settings. As the Geminga halo is expected to
be very extended in GeV, the spectral measurement is challenging owing to the
large uncertainty of the gamma-ray background. More careful analyses are
needed to cross-check the current measurements.

Figure 6: Left: Interpretation of the positron excess with Geminga under the
two-zone diffusion model ($r_{\star}=50$ pc). The result of the one-zone slow-
diffusion model ($r_{\star}=\infty$) is also shown for comparison. This figure
is adapted from Fang et al. (2018). Right: Contours of the positron flux at
300 GeV as a function of pulsar age and distance. Pulsars inside the shaded
areas are those can contribute more than half of the measured flux of AMS-02
at 300 GeV. Three different models are presented with different colors. The
red band indicates the range of the cooling time corresponding to positron
spectral cutoff measured by AMS-02 (Aguilar et al., 2019). One may refer to
Fang et al. (2019b) for further details of this figure.

###  7.2 Other Candidate Pulsars

The latest positron spectrum measured by the AMS collaboration shows a sharp
spectral dropoff around 300 GeV with a significance of more than $3\sigma$
(Aguilar et al., 2019). This result provides important clues for the origin of
the positron excess. The spectral cutoff can hardly be interpreted by the
superposition of multiple pulsars (Manconi et al., 2020). Even old PWNe like
the Geminga PWN can accelerate electrons to $\sim 100$ TeV, so it is
unreasonable to assume a typical cutoff energy of $\lesssim 1$ TeV for the
injection spectrum of PWNe to interpret the data. Instead, the spectral cutoff
is more likely to correspond to the energy loss of the positrons released by a
single source. In other words, the positron excess could be dominated by a
single source (Fang et al., 2019b).

Although Geminga is still a likely source of the positron excess, it is
meaningful to investigate other possible dominant sources considering the
slow-diffusion phenomenon. Fang et al. (2019b) scanned all the identified
pulsars333Millisecond pulsars are not included in the analysis of Fang et al.
(2019b), while the nearest millisecond pulsar, PSR J0437$-$4715, is proposed
to be a likely source of the positron excess (Bykov et al., 2019). and
searched for the pulsars that can both interpret the positron intensity and
the spectral cutoff of the AMS-02 spectrum. The right panel of Figure 6 shows
the contours of the positron flux at 300 GeV as a function of pulsar age and
distance. Pulsars inside the shaded areas are those can contribute more than
half of the measured flux of AMS-02 at 300 GeV. Three different models are
presented with different colors: the one-zone slow-diffusion model, the one-
zone fast diffusion model, and a two-zone diffusion model with $r_{\star}=50$
pc. The red band indicates the range of the energy-loss time corresponding to
positron spectral cutoff. Thus, the candidate pulsars should lie in the
intersection region between a shaded area and the red band. As can be seen,
Geminga and PSR B1055$-$52 are the most likely candidates. Another interesting
finding is that the positron excess may not be interpreted by pulsars without
the assumption of slow-diffusion zones, since the blue shaded area has no
overlap with the red band.

PSR B1055$-$52 is a Geminga-like pulsar, while its distance was thought to be
larger than 1 kpc. The latest dispersion measure derives a close distance of
90 pc for PSR B1055$-$52 (Yao et al., 2017), and it becomes a very bright
positron source. The annual parallax measurement for PSR B1055$-$52 is
suggested to obtain a more precise distance. The Fermi-LAT observation for
extended emission is also essential to constrain its injection power.

##  8 Diffuse TeV Gamma-ray Excess

In this section, we briefly introduce the possible connection between pulsar
halos and the diffuse gamma-ray excess. The main origin of the diffuse gamma-
ray emission from the Galactic plane is generally believed to be the
$\pi^{0}$-decay process from the interaction between the propagating CR nuclei
and the ISM (Strong et al., 2007). The diffuse gamma-ray emission below 100
GeV measured by Fermi-LAT and the local CR nuclei spectrum can be consistently
interpreted under the homogeneous CR propagation model (Ackermann et al.,
2012a). At the TeV energy region, however, the $\pi^{0}$-decay emission
predicted by the local CR spectrum is significantly lower than the diffuse
gamma-ray measurements of Milagro and ARGO-YBJ experiments (Atkins et al.,
2005; Abdo et al., 2008; Bartoli et al., 2015), which is known as the diffuse
TeV gamma-ray excess (Prodanovic et al., 2007; Linden and Buckman, 2018).
Spatially dependent propagation models are constructed to solve this
contradiction (Gaggero et al., 2015).

Alternatively, the excess could be attributed to unresolved gamma-ray sources.
Linden and Buckman (2018) pointed out that pulsar halos and extended gamma-ray
PWNe could be promising candidates for this unknown component. If all the
young and middle-aged pulsars are associated with extended gamma-ray sources,
an average conversion efficiency of $\sim 10\%$ from the pulsar spin-down
energy to the high-energy electrons could be enough to account for the excess
(Linden and Buckman, 2018). This model can naturally avoid the constraint from
the local CR nuclei spectrum. Moreover, the gamma-ray spectra of pulsar halos
and young PWNe are much harder than the diffuse $\pi^{0}$-decay spectrum,
which means that this component may account for the TeV measurements while
keeping consistency with the Fermi-LAT observations. On the other hand, the
potential contribution of pulsar halos for the excess may vary with the
mechanism of pulsar halos introduced in Section 4. One may refer to Liu (2022)
for a comprehensive discussion on this topic.

Recently, diffuse gamma-ray above 100 TeV was reported by the Tibet AS+MD
experiment (Amenomori et al., 2021). If the $\pi^{0}$-decay component
dominates the diffuse emission, there will be tension between the Tibet AS+MD
measurement and the local CR spectrum even under the spatially dependent CR
propagation model (Qiao et al., 2022). Liu and Wang (2021) pointed out that
the upper limit of the Galactic neutrino emission may constrain the
$\pi^{0}$-decay component, also indicating the existence of extra sources. The
LHAASO experiment will test whether pulsar halos could contribute above 100
TeV in the near future. The sub-PeV excess may also be interpreted by the
potential PeV accelerators, such as massive star clusters, hypernova remnants,
or young PWNe (Liu and Wang, 2021).

##  9 Future Works

Our understanding of pulsar halos is still in the preliminary stage. Future
studies on pulsar halos should be carried out in both depth and breadth.
Deeper observations of bright pulsar halos can provide insight into the CR
propagation mechanism at small scales and then the properties of the MHD
turbulence in the ISM. At the same time, a larger sample may also help us
understand the origin of pulsar halos and the relation between local and
global propagation of Galactic CRs.

###  9.1 Deeper Studies on Bright Pulsar Halos

Morphology measurements of bright and nearby pulsar halos like the Geminga
halo can provide the richest details about small-scale CR propagation. As
shown in Figures 3 and 4, the gamma-ray profiles predicted by different
propagation models have different characteristics, while the current
experimental results may not distinguish them. For the superdiffusion
scenario, only the models with $\alpha$ close to $1$ are strongly disfavored.
For the relativistic-modified fast-diffusion scenario, the goodness-of-fit
test disfavors the model at a confidence level of $98.6\%$ for the Geminga
halo case (Bao et al., 2022), which has not reached the $3\sigma$ level. The
ongoing LHAASO experiment (Ma et al., 2022) is expected to give more definite
judgments to these models with better sensitivity and angular resolution.
LHAASO may also be able to detect weak asymmetry in extensive halos, which is
expected by the anisotropic diffusion model. Thus, more accurate morphology
measurements are essential to test the mechanisms proposed in Sections 4.3 and
4.4 and hence determine whether pulsar halos originated from slow diffusion.

Energy-dependent morphology measurement of pulsar halos is not available at
present but will be an essential subject for future study. This measurement
determines the energy dependency of the CR diffusion coefficient, which
reflects the properties of the MHD turbulence in the ISM and may also indicate
the origin of pulsar halos. If the turbulence originated from the Alfvénic or
slow magnetosonic cascade, the energy dependency could be $D(E)\propto
E^{1/3}$ (Goldreich and Sridhar, 1995; Cho and Lazarian, 2002), while the fast
magnetosonic scenario predicts $D(E)\propto E^{1/2}$ (Cho and Lazarian, 2002).
For the self-generated scenario, the energy dependency of $D$ is more
complicated and may not be described by a power law (Evoli et al., 2018b;
Mukhopadhyay and Linden, 2022). Current observations suggest that the
injection spectra of pulsar halos have a cutoff at $\sim 100$ TeV. As the
turbulence growth rate is positively correlated with the injection power under
the self-generated mechanism, a sharp increase of $D$ is predicted at high
energies (Mukhopadhyay and Linden, 2022). Besides, $D(E)\propto E$ may
correspond to the Bohm diffusion, which means that the turbulent growth has
reached saturation.

The energy spectra of pulsar halos are measured in a certain angular range and
thus determined by both the injection spectra and $D(E)$. Energy-dependent
morphology measurements are indispensable to decouple the injection spectra
and $D(E)$. LHAASO-KM2A combined with the Water Cherenkov Detector Array of
LHAASO (LHAASO-WCDA) can measure the morphologies of pulsar halos in the
energy range from $\sim 1$ TeV to $\sim 1$ PeV, which could be wide enough to
understand the forms of the injection spectra and $D(E)$. Bright pulsar halos
are feasible targets for energy-dependent measurements as it is possible to
divide more energy bins for these sources.

We should note that the comparison between the localized diffusion coefficient
inferred from pulsar halos and the global diffusion coefficient inferred from
B/C is based on the energy extrapolation of the result measured by B/C. For
example, the global diffusion coefficient suggested by the diffusion-
reacceleration model of Yuan et al. (2017) is $\approx 700$ times larger than
the value in the Geminga halo if extrapolated to 100 TeV. However, the B/C
measurement is only up to $\sim 1$ TeV at present, and the energy dependency
of the global diffusion coefficient is not available at higher energies. On
the one hand, the ongoing DArk Matter Particle Explorer (Chang et al., 2017,
DAMPE) will provide multi-TeV measurement of B/C in the coming future. We may
also expect future space experiments like the High Energy cosmic-Radiation
Detection (Zhang et al., 2014, HERD) facility to boost the B/C measurement to
$\sim 100$ TeV. On the other hand, we may measure the sub-TeV diffusion
coefficient in pulsar halos by Fermi-LAT or H.E.S.S to give a straightforward
comparison with the Galactic-averaged value.

###  9.2 A Larger Sample of Pulsar Halos

When a large sample of pulsar halos is available, we can statistically
investigate the correlations between the parameters of pulsar halos. This is
helpful in judging the origin of pulsar halos. The self-generated model will
be disfavored if there is no correlation between the diffusion coefficient and
the electron injection power. For example, as noticed by Fang et al. (2022),
the diffusion coefficient in HESS J1831$-$098 is comparable to the other known
pulsar halos, while the spin-down luminosity of PSR J1831$-$0952 is dozens of
times larger than the other pulsars. This implies that the slow-diffusion zone
around PSR J1831$-$0952 may not be self-generated. Another possibility is that
the slow-diffusion zone is self-generated while the turbulent growth has
reached saturation for the known pulsar halos. The latter case will be
possible if $D(E)\propto E$.

At the same time, if the diffusion coefficient in pulsar halos is independent
of the position in the Galaxy, the slow-diffusion environment is more likely
to be self-generated or produced by the associated SNR rather than the
superposition effect of multiple external turbulent sources. An example is
LHAASO J0621$-$3755\. It is far from the Galactic plane, while its diffusion
coefficient is comparable to the other pulsar halos near the Galactic plane.
It indicates that the slow diffusion may not only occur in the Galactic disk
but as a localized phenomenon around some middle-aged pulsars. Moreover,
searching for extended gamma-ray emission around millisecond pulsars is
important to test the SNR-generated scenario as these pulsars have been far
away from their parent SNRs.

The LHAASO experiment is expected to discover dimmer pulsar halos owing to the
better sensitivity. This is crucial for inferring the statistical properties
of pulsar halos. Assuming a fixed electron injection power, pulsar halos with
a larger diffusion coefficient will be dimmer. Thus, if the sample only
includes the brightest pulsar halos, the inference will be seriously affected
by the selection effect.

In addition, the local environment of pulsars, such as the interstellar
radiation field and the magnetic field strength, can also affect the
brightness of pulsar halos. These factors should be considered in the
statistical analyses. If the magnetic field strength has a significant spatial
variation in the local environment, the electron spatial distribution will be
affected, hence the diffusion coefficient measurement. Observations of diffuse
synchrotron emission in the X-ray band may help to test this scenario (Li et
al., 2022).

####  9.2.1 Are Pulsar Halos Universal?

The universality of pulsar halos is often mentioned in the prospect of future
observations (Linden et al., 2017; Sudoh et al., 2019). It may answer whether
the localized slow-diffusion around pulsars is related to the inferred
spatially dependent CR diffusion at the Galactic scale as discussed in Section
5.2 and how pulsar halos contribute to the excess of the diffuse TeV emission
introduced in Section 8.

Current observations indicate that pulsar halos may not be universal. PSR
J1809$-$2332, B0906$-$49, and J1105$-$6107 are the brightest middle-aged
pulsars within the survey region of the H.E.S.S. Galactic plane survey (HGPS).
Their $L/D^{2}$ is more than $5$ times larger than that of PSR J1831$-$0952,
where $L$ and $D$ are the pulsar spin-down luminosity and distance,
respectively. Significant TeV emission is expected around these pulsars if
pulsar halos are universal; however, no signal is detected by HGPS (Abdalla et
al., 2018). Moreover, if all the nearby middle-aged pulsars have the same
slow-diffusion environment and injection efficiency as the Geminga halo, the
predicted positron flux at Earth will be much higher than observed (Fang et
al., 2019b; Martin et al., 2022). Based on this, Martin et al. (2022) point
out that the occurrence rate of pulsar halos could be $\sim 5-10\%$.

It should be noted that the non-universality of pulsar halos is not equal to
the non-universality of the slow-diffusion environment around middle-aged
pulsars. The brightness of a pulsar halo is also determined by the electron
acceleration efficiency and escape efficiency of the center PWN. The reason
for the invisibility of pulsar halos is not necessarily the same for different
objects. For example, among the three bright pulsars mentioned in the last
paragraph, PSR J1809$-$2332 has a significant bow-shock X-ray PWN (Van Etten
et al., 2012). It indicates that this source can indeed generate high-energy
electrons, and the invisibility of the pulsar halo is possibly due to the low
escape efficiency from the PWN or the lack of slow-diffusion environment.
However, no X-ray PWN is observed for the other two pulsars, which means that
their ability to generate high-energy electrons is questionable.

We should also note that although the factors of the slow diffusion and
acceleration efficiency may be degenerate for the visibility of a pulsar halo,
the implication could be very different for the interpretation of positron
excess or the diffuse gamma-ray emission. In the absence of slow diffusion,
the electrons/positrons generated by the pulsar may still contribute to the
positron excess or the diffuse emission, while the opposite is true in the
absence of high-energy electron acceleration.

##  10 Summary

We briefly summarize the key points of this review as follows:

  1. 1.

Pulsar halos are a new class of gamma-ray sources represented by the Geminga
halo. They are generated by electrons diffusing in the ISM around pulsars. It
is essential to distinguish them from PWNe: PWNe should be regarded as the
electron sources of pulsar halos.

  2. 2.

The diffusion inferred from the morphologies of pulsar halos is extremely
slow, while the origin of the slow-diffusion environment is still uncertain.
It may be self-generated or left by external sources like SNRs. Some other
models may interpret the pulsar halos without slow diffusion.

  3. 3.

Pulsar halos are ’microscopes’ for studying the Galactic CR propagation, which
can test sophisticated propagation models at small scales. They are
complementary to the global probes of Galactic CR propagation like B/C.

  4. 4.

Pulsar halos are essential for inferring the electron injection spectra from
PWNe, which help us understand the particle escape process from PWNe.

  5. 5.

The discovery of pulsar halos has a significant influence on the
interpretation of the positron excess. Geminga is still the most likely source
of the positron excess. Deeper studies of the Geminga halo at $\sim 10$ GeV
are crucial to this issue.

  6. 6.

Pulsar halos may contribute to the unknown excess of the diffuse TeV gamma-ray
emission.

  7. 7.

Pulsar halos are probably not universal, which may be attributed to the non-
universality of the slow-diffusion environment or the difference in the
electron injection efficiency of PWNe.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Funding

This work is supported by the National Natural Science Foundation of China
under Grants No. 12105292, and No. U1738209.

## Acknowledgments

The author thanks Prof. Xiao-Jun Bi for providing motivation for writing this
paper. The author also thanks the reviewers for their detailed and valuable
suggestions.

## References

  * Abdalla et al. (2018) Abdalla, H. et al. (2018).  The H.E.S.S. Galactic plane survey.  _Astron. Astrophys._ 612, A1.  10.1051/0004-6361/201732098
  * Abdo et al. (2007) Abdo, A. A. et al. (2007).  TeV Gamma-Ray Sources from a Survey of the Galactic Plane with Milagro.  _Astrophys. J. Lett._ 664, L91–L94.  10.1086/520717
  * Abdo et al. (2008) Abdo, A. A. et al. (2008).  A Measurement of the Spatial Distribution of Diffuse TeV Gamma Ray Emission from the Galactic Plane with Milagro.  _Astrophys. J._ 688, 1078–1083.  10.1086/592213
  * Abdo et al. (2013) Abdo, A. A. et al. (2013).  The Second Fermi Large Area Telescope Catalog of Gamma-ray Pulsars.  _Astrophys. J. Suppl._ 208, 17.  10.1088/0067-0049/208/2/17
  * Abeysekara et al. (2017) Abeysekara, A. et al. (2017).  Extended gamma-ray sources around pulsars constrain the origin of the positron flux at Earth.  _Science_ 358, 911–914.  10.1126/science.aan4880
  * Abichandani et al. (2019) Abichandani, R., Mathur, M. B., Drake, J. J., Fruscione, A., Lee, N. P., and Glotfelty, K. (2019).  On the Possible Pulsar Wind Nebula Counterpart of the TeV Source HESS J1831-098.  _The Astronomer’s Telegram_ 12463, 1
  * Abramowski et al. (2011) Abramowski, A. et al. (2011).  Revisiting the Westerlund 2 Field with the H.E.S.S. Telescope Array.  _Astron. Astrophys._ 525, A46.  10.1051/0004-6361/201015290
  * Abramowski et al. (2012) Abramowski, A. et al. (2012).  Probing the extent of the non-thermal emission from the Vela X region at TeV energies with H.E.S.S.  _Astron. Astrophys._ 548, A38.  10.1051/0004-6361/201219919
  * Ackermann et al. (2012a) Ackermann, M. et al. (2012a).  Fermi-LAT Observations of the Diffuse Gamma-Ray Emission: Implications for Cosmic Rays and the Interstellar Medium.  _Astrophys. J._ 750, 3.  10.1088/0004-637X/750/1/3
  * Ackermann et al. (2012b) Ackermann, M. et al. (2012b).  Measurement of separate cosmic-ray electron and positron spectra with the Fermi Large Area Telescope.  _Phys. Rev. Lett._ 108, 011103.  10.1103/PhysRevLett.108.011103
  * Ackermann et al. (2015) Ackermann, M. et al. (2015).  Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data.  _Phys. Rev. Lett._ 115, 231301.  10.1103/PhysRevLett.115.231301
  * Adriani et al. (2009) Adriani, O. et al. (2009).  An anomalous positron abundance in cosmic rays with energies 1.5-100 GeV.  _Nature_ 458, 607–609.  10.1038/nature07942
  * Aguilar et al. (2013) Aguilar, M. et al. (2013).  First Result from the Alpha Magnetic Spectrometer on the International Space Station: Precision Measurement of the Positron Fraction in Primary Cosmic Rays of 0.5–350 GeV.  _Phys. Rev. Lett._ 110, 141102.  10.1103/PhysRevLett.110.141102
  * Aguilar et al. (2019) Aguilar, M. et al. (2019).  Towards Understanding the Origin of Cosmic-Ray Positrons.  _Phys. Rev. Lett._ 122, 041102.  10.1103/PhysRevLett.122.041102
  * Aharonian et al. (2021) Aharonian, F. et al. (2021).  Extended Very-High-Energy Gamma-Ray Emission Surrounding PSR J0622 + 3749 Observed by LHAASO-KM2A.  _Phys. Rev. Lett._ 126, 241103.  10.1103/PhysRevLett.126.241103
  * Albert et al. (2020) Albert, A. et al. (2020).  3HWC: The Third HAWC Catalog of Very-High-Energy Gamma-ray Sources.  _Astrophys. J._ 905, 76.  10.3847/1538-4357/abc2d8
  * Aloisio et al. (2009) Aloisio, R., Berezinsky, V., and Gazizov, A. (2009).  Superluminal problem in diffusion of relativistic particles and its phenomenological solution.  _Astrophys. J._ 693, 1275–1282.  10.1088/0004-637X/693/2/1275
  * Amenomori et al. (2021) Amenomori, M. et al. (2021).  First Detection of sub-PeV Diffuse Gamma Rays from the Galactic Disk: Evidence for Ubiquitous Galactic Cosmic Rays beyond PeV Energies.  _Phys. Rev. Lett._ 126, 141101.  10.1103/PhysRevLett.126.141101
  * An et al. (2019) An, Q. et al. (2019).  Measurement of the cosmic-ray proton spectrum from 40 GeV to 100 TeV with the DAMPE satellite.  _Sci. Adv._ 5, eaax3793.  10.1126/sciadv.aax3793
  * Atkins et al. (2005) Atkins, R. et al. (2005).  Evidence for TeV gamma-ray emission from the galactic plane.  _Phys. Rev. Lett._ 95, 251103.  10.1103/PhysRevLett.95.251103
  * Balbo et al. (2010) Balbo, M., Saouter, P., Walter, R., Pavan, L., Tramacere, A., Pohl, M., et al. (2010).  HESS J1632-478: an energetic relic.  _Astron. Astrophys._ 520, A111.  10.1393/ncc/i2011-10893-x
  * Bao et al. (2022) Bao, L.-Z., Fang, K., Bi, X.-J., and Wang, S.-H. (2022).  Slow Diffusion is Necessary to Explain the $\gamma$-Ray Pulsar Halos.  _Astrophys. J._ 936, 183.  10.3847/1538-4357/ac8b8a
  * Bartoli et al. (2015) Bartoli, B. et al. (2015).  Study of the Diffuse Gamma-ray Emission From the Galactic Plane With ARGO-YBJ.  _Astrophys. J._ 806, 20.  10.1088/0004-637X/806/1/20
  * Bell (1978) Bell, A. R. (1978).  The Acceleration of cosmic rays in shock fronts. I.  _Mon. Not. Roy. Astron. Soc._ 182, 147–156
  * Bell (2004) Bell, A. R. (2004).  Turbulent amplification of magnetic field and diffusive shock acceleration of cosmic rays.  _Mon. Not. R. Astron. Soc._ 353, 550–558.  10.1111/j.1365-2966.2004.08097.x
  * Bergstrom et al. (2008) Bergstrom, L., Bringmann, T., and Edsjo, J. (2008).  New Positron Spectral Features from Supersymmetric Dark Matter - a Way to Explain the PAMELA Data?  _Phys. Rev. D_ 78, 103520.  10.1103/PhysRevD.78.103520
  * Bîrzan et al. (2016) Bîrzan, L., Pavlov, G. G., and Kargaltsev, O. (2016).  Chandra observations of the elusive pulsar wind nebula around PSR B0656+14.  _Astrophys. J._ 817, 129.  10.3847/0004-637X/817/2/129
  * Bitter and Hooper (2022) Bitter, O. M. and Hooper, D. (2022).  Constraining the Milky Way’s Pulsar Population with the Cosmic-Ray Positron Fraction.  _arXiv e-prints_ , arXiv:2205.05200
  * Brisbois et al. (2018) Brisbois, C., Riviere, C., Fleischhack, H., and Smith, A. (2018).  HAWC detection of TeV source HAWC J0635+070.  _The Astronomer’s Telegram_ 12013, 1
  * Bykov et al. (2017) Bykov, A. M., Amato, E., Petrov, A. E., Krassilchtchikov, A. M., and Levenfish, K. P. (2017).  Pulsar wind nebulae with bow shocks: non-thermal radiation and cosmic ray leptons.  _Space Sci. Rev._ 207, 235–290.  10.1007/s11214-017-0371-7
  * Bykov et al. (2019) Bykov, A. M., Petrov, A. E., Krassilchtchikov, A. M., Levenfish, K. P., Osipov, S. M., and Pavlov, G. G. (2019).  GeV-TeV Cosmic Ray Leptons in the Solar System from the Bow Shock Wind Nebula of the Nearest Millisecond Pulsar J0437-4715.  _Astrophys. J. Lett._ 876, L8.  10.3847/2041-8213/ab1922
  * Caraveo et al. (2003) Caraveo, P. A., Bignami, G. F., De Luca, A., Mereghetti, S., Pellizzoni, A., Mignani, R., et al. (2003).  Geminga’s Tails: A Pulsar Bow Shock Probing the Interstellar Medium.  _Science_ 301, 1345–1348.  10.1126/science.1086973
  * Chang et al. (2017) Chang, J. et al. (2017).  The DArk Matter Particle Explorer mission.  _Astropart. Phys._ 95, 6–24.  10.1016/j.astropartphys.2017.08.005
  * Cho and Lazarian (2002) Cho, J. and Lazarian, A. (2002).  Compressible sub-alfvenic MHD turbulence in low-Beta plasmas.  _Phys. Rev. Lett._ 88, 245001.  10.1103/PhysRevLett.88.245001
  * Cholis et al. (2009) Cholis, I., Goodenough, L., Hooper, D., Simet, M., and Weiner, N. (2009).  High Energy Positrons From Annihilating Dark Matter.  _Phys. Rev. D_ 80, 123511.  10.1103/PhysRevD.80.123511
  * D’Angelo et al. (2016) D’Angelo, M., Blasi, P., and Amato, E. (2016).  Grammage of cosmic rays around Galactic supernova remnants.  _Phys. Rev. D_ 94, 083003.  10.1103/PhysRevD.94.083003
  * De La Torre Luque et al. (2022) De La Torre Luque, P., Fornieri, O., and Linden, T. (2022).  Anisotropic diffusion cannot explain TeV halo observations.  _arXiv e-prints_ , arXiv:2205.08544
  * Della Torre et al. (2015) Della Torre, S., Gervasi, M., Rancoita, P. G., Rozza, D., and Treves, A. (2015).  Pulsar Wind Nebulae as a source of the observed electron and positron excess at high energy: the case of Vela-X.  _JHEAp_ 8, 27–34.  10.1016/j.jheap.2015.08.001
  * Di Mauro et al. (2019) Di Mauro, M., Manconi, S., and Donato, F. (2019).  Detection of a $\gamma$-ray halo around Geminga with the Fermi -LAT data and implications for the positron flux.  _Phys. Rev. D_ 100, 123015.  10.1103/PhysRevD.104.089903.  [Erratum: Phys.Rev.D 104, 089903 (2021)]
  * Di Mauro et al. (2020) Di Mauro, M., Manconi, S., and Donato, F. (2020).  Evidences of low-diffusion bubbles around Galactic pulsars.  _Phys. Rev. D_ 101, 103035.  10.1103/PhysRevD.101.103035
  * Dunkel et al. (2007) Dunkel, J., Talkner, P., and Hanggi, P. (2007).  Relativistic diffusion processes and random walk models.  _Phys. Rev. D_ 75, 043001.  10.1103/PhysRevD.75.043001
  * Evoli et al. (2018a) Evoli, C., Blasi, P., Morlino, G., and Aloisio, R. (2018a).  Origin of the Cosmic Ray Galactic Halo Driven by Advected Turbulence and Self-Generated Waves.  _Phys. Rev. Lett._ 121, 021102.  10.1103/PhysRevLett.121.021102
  * Evoli et al. (2018b) Evoli, C., Linden, T., and Morlino, G. (2018b).  Self-generated cosmic-ray confinement in TeV halos: Implications for TeV $\gamma$-ray emission and the positron excess.  _Phys. Rev. D_ 98, 063017.  10.1103/PhysRevD.98.063017
  * Evoli and Yan (2014) Evoli, C. and Yan, H. (2014).  Cosmic ray propagation in galactic turbulence.  _Astrophys. J._ 782, 36.  10.1088/0004-637X/782/1/36
  * Faherty et al. (2007) Faherty, J., Walter, F. M., and Anderson, J. (2007).  The trigonometric parallax of the neutron star Geminga.  _Astrophys. Space Sci._ 308, 225–230.  10.1007/s10509-007-9368-0
  * Fang and Bi (2022) Fang, K. and Bi, X.-J. (2022).  Interpretation of the puzzling gamma-ray spectrum of the Geminga halo.  _Phys. Rev. D_ 105, 103007.  10.1103/PhysRevD.105.103007
  * Fang et al. (2019a) Fang, K., Bi, X.-J., and Yin, P.-F. (2019a).  Possible origin of the slow-diffusion region around Geminga.  _Mon. Not. Roy. Astron. Soc._ 488, 4074–4080.  10.1093/mnras/stz1974
  * Fang et al. (2019b) Fang, K., Bi, X.-J., and Yin, P.-F. (2019b).  Reanalyze the pulsar scenario to explain the cosmic positron excess considering the recent developments.  _Astrophys. J._ 884, 124–128.  10.3847/1538-4357/ab3fac
  * Fang et al. (2020) Fang, K., Bi, X.-J., and Yin, P.-F. (2020).  DAMPE proton spectrum indicates a slow-diffusion zone in the nearby ISM.  _Astrophys. J._ 903, 69.  10.3847/1538-4357/abb8d7
  * Fang et al. (2018) Fang, K., Bi, X.-J., Yin, P.-F., and Yuan, Q. (2018).  Two-zone diffusion of electrons and positrons from Geminga explains the positron anomaly.  _Astrophys. J._ 863, 30.  10.3847/1538-4357/aad092
  * Fang et al. (2022) Fang, K., Xi, S.-Q., Bao, L.-Z., Bi, X.-J., and Chen, E.-S. (2022).  Features of the gamma-ray pulsar halo HESS J1831$-$098\.  _arXiv e-prints_ , arXiv:2207.13533
  * Fang et al. (2021) Fang, K., Xi, S.-Q., and Bi, X.-J. (2021).  Self-consistent interpretations of the multiwavelength gamma-ray spectrum of LHAASO J0621+3755.  _Phys. Rev. D_ 104, 103024.  10.1103/PhysRevD.104.103024
  * Feng et al. (2016) Feng, J., Tomassetti, N., and Oliva, A. (2016).  Bayesian analysis of spatial-dependent cosmic-ray propagation: Astrophysical background of antiprotons and positrons.  _Phys. Rev. D_ 94, 123007.  10.1103/PhysRevD.94.123007
  * Ferriere (2020) Ferriere, K. (2020).  Plasma turbulence in the interstellar medium.  _Plasma Phys. Control. Fusion_ 62, 014014.  10.1088/1361-6587/ab49eb
  * Fleischhack et al. (2019) Fleischhack, H., Albert, A., Alvarez, C., Arceo, R., Solares, H. A. A., Beacom, J. F., et al. (2019).  Pulsars in a Bubble? Following Electron Diffusion in the Galaxy with TeV Gamma Rays.  _Bulletin of the American Astronomical Society_ 51, 311
  * Fujita et al. (2009) Fujita, Y., Ohira, Y., Tanaka, S. J., and Takahara, F. (2009).  Molecular Clouds as a Probe of Cosmic-Ray Acceleration in a Supernova Remnant.  _Astrophys. J. Lett._ 707, L179–L183.  10.1088/0004-637X/707/2/L179
  * Gaensler and Slane (2006) Gaensler, B. M. and Slane, P. O. (2006).  The evolution and structure of pulsar wind nebulae.  _Ann. Rev. Astron. Astrophys._ 44, 17–47.  10.1146/annurev.astro.44.051905.092528
  * Gaggero et al. (2015) Gaggero, D., Urbano, A., Valli, M., and Ullio, P. (2015).  Gamma-ray sky points to radial gradients in cosmic-ray transport.  _Phys. Rev. D_ 91, 083012.  10.1103/PhysRevD.91.083012
  * Giacalone and Jokipii (2007) Giacalone, J. and Jokipii, J. R. (2007).  Magnetic Field Amplification by Shocks in Turbulent Fluids.  _Astrophys. J. Lett._ 663, L41–L44.  10.1086/519994
  * Giacinti et al. (2020) Giacinti, G., Mitchell, A., López-Coto, R., Joshi, V., Parsons, R., and Hinton, J. (2020).  Halo fraction in TeV-bright pulsar wind nebulae.  _Astron. Astrophys._ 636, A113.  10.1051/0004-6361/201936505
  * Ginzburg and Syrovatskii (1964) Ginzburg, V. L. and Syrovatskii, S. I. (1964).  _The Origin of Cosmic Rays_ (New York: Macmillan)
  * Goldreich and Sridhar (1995) Goldreich, P. and Sridhar, S. (1995).  Toward a theory of interstellar turbulence. 2. Strong Alfvenic turbulence.  _Astrophys. J._ 438, 763–775.  10.1086/175121
  * Guo and Yuan (2018) Guo, Y.-Q. and Yuan, Q. (2018).  Understanding the spectral hardenings and radial distribution of Galactic cosmic rays and Fermi diffuse $\gamma$ rays with spatially-dependent propagation.  _Phys. Rev. D_ 97, 063008.  10.1103/PhysRevD.97.063008
  * Han (2017) Han, J. L. (2017).  Observing Interstellar and Intergalactic Magnetic Fields.  _Annu. Rev. Astron. Astrophys._ 55, 111–157.  10.1146/annurev-astro-091916-055221
  * Hooper et al. (2009) Hooper, D., Blasi, P., and Serpico, P. D. (2009).  Pulsars as the Sources of High Energy Cosmic Ray Positrons.  _JCAP_ 01, 025.  10.1088/1475-7516/2009/01/025
  * Hooper et al. (2017) Hooper, D., Cholis, I., Linden, T., and Fang, K. (2017).  HAWC Observations Strongly Favor Pulsar Interpretations of the Cosmic-Ray Positron Excess.  _Phys. Rev. D_ 96, 103013.  10.1103/PhysRevD.96.103013
  * Hooper and Linden (2018) Hooper, D. and Linden, T. (2018).  Millisecond Pulsars, TeV Halos, and Implications For The Galactic Center Gamma-Ray Excess.  _Phys. Rev. D_ 98, 043005.  10.1103/PhysRevD.98.043005
  * Hooper and Linden (2022) Hooper, D. and Linden, T. (2022).  Evidence of TeV halos around millisecond pulsars.  _Phys. Rev. D_ 105, 103013.  10.1103/PhysRevD.105.103013
  * Huang et al. (2012) Huang, X., Yuan, Q., Yin, P.-F., Bi, X.-J., and Chen, X. (2012).  Constraints on the dark matter annihilation scenario of Fermi 130 GeV $\gamma$-ray line emission by continuous gamma-rays, Milky Way halo, galaxy clusters and dwarf galaxies observations.  _JCAP_ 11, 048.  10.1088/1475-7516/2012/11/048.  [Erratum: JCAP 05, E02 (2013)]
  * Johannesson et al. (2019) Johannesson, G., Porter, T. A., and Moskalenko, I. V. (2019).  Cosmic-Ray Propagation in Light of the Recent Observation of Geminga.  _Astrophys. J._ 879, 91.  10.3847/1538-4357/ab258e
  * Kim et al. (2007) Kim, I. J., Min, K. W., Seon, K. I., Park, J. W., Han, W., Park, J. H., et al. (2007).  Far-Ultraviolet Observations of the Monogem Ring.  _Astrophys. J. Lett._ 665, L139–L142.  10.1086/521441
  * Knies et al. (2018) Knies, J. R., Sasaki, M., and Plucinsky, P. P. (2018).  Suzaku observations of the Monogem Ring and the origin of the Gemini H $\alpha$ ring.  _Mon. Not. R. Astron. Soc._ 477, 4414–4422.  10.1093/mnras/sty915
  * Lagutin et al. (2001) Lagutin, A. A., Nikulin, Y. A., and Uchaikin, V. V. (2001).  The “knee” in the primary cosmic ray spectrum as consequence of the anomalous diffusion of the particles in the fractal interstellar medium.  _Nuclear Physics B Proceedings Supplements_ 97, 267–270.  10.1016/S0920-5632(01)01280-4
  * Lagutin and Uchaikin (2003) Lagutin, A. A. and Uchaikin, V. V. (2003).  Anomalous diffusion equation: Application to cosmic ray transport.  _Nuclear Instruments and Methods in Physics Research B_ 201, 212–216.  10.1016/S0168-583X(02)01362-9
  * Leahy and Williams (2017) Leahy, D. A. and Williams, J. (2017).  A Python Calculator for Supernova Remnant Evolution.  _Astron. J._ 153, 239.  10.3847/1538-3881/aa6af6
  * Li et al. (2022) Li, B., Zhang, Y., Liu, T., Liu, R.-Y., and Wang, X.-Y. (2022).  Prospect of detecting X-ray haloes around middle-aged pulsars with eROSITA.  _Mon. Not. Roy. Astron. Soc._ 513, 2884–2892.  10.1093/mnras/stac711
  * Li and Chen (2012) Li, H. and Chen, Y. (2012).  Gamma-rays from molecular clouds illuminated by accumulated diffusive protons. II: interacting supernova remnants.  _Mon. Not. Roy. Astron. Soc._ 421, 935.  10.1111/j.1365-2966.2012.20270.x
  * Lin et al. (2015) Lin, S.-J., Yuan, Q., and Bi, X.-J. (2015).  Quantitative study of the AMS-02 electron/positron spectra: Implications for pulsars and dark matter properties.  _Phys. Rev. D_ 91, 063508.  10.1103/PhysRevD.91.063508
  * Linden et al. (2017) Linden, T., Auchettl, K., Bramante, J., Cholis, I., Fang, K., Hooper, D., et al. (2017).  Using HAWC to discover invisible pulsars.  _Phys. Rev. D_ 96, 103016.  10.1103/PhysRevD.96.103016
  * Linden and Buckman (2018) Linden, T. and Buckman, B. J. (2018).  Pulsar TeV Halos Explain the Diffuse TeV Excess Observed by Milagro.  _Phys. Rev. Lett._ 120, 121101.  10.1103/PhysRevLett.120.121101
  * Liu (2022) Liu, R.-Y. (2022).  The Physics of Pulsar Halos: Research Progress and Prospect.  _arXiv e-prints_ , arXiv:2207.04011
  * Liu et al. (2019a) Liu, R.-Y., Ge, C., Sun, X.-N., and Wang, X.-Y. (2019a).  Constraining the Magnetic Field in the TeV Halo of Geminga with X-Ray Observations.  _Astrophys. J._ 875, 149.  10.3847/1538-4357/ab125c
  * Liu and Wang (2021) Liu, R.-Y. and Wang, X.-Y. (2021).  Origin of Galactic Sub-PeV Diffuse Gamma-Ray Emission: Constraints from High-energy Neutrino Observations.  _Astrophys. J. Lett._ 914, L7.  10.3847/2041-8213/ac02c5
  * Liu et al. (2019b) Liu, R.-Y., Yan, H., and Zhang, H. (2019b).  Understanding the Multiwavelength Observation of Geminga’s Tev Halo: The Role of Anisotropic Diffusion of Particles.  _Phys. Rev. Lett._ 123, 221103.  10.1103/PhysRevLett.123.221103
  * Liu et al. (2017) Liu, W., Bi, X.-J., Lin, S.-J., and Yin, P.-F. (2017).  Constraints on dark matter annihilation and decay from the isotropic gamma-ray background.  _Chin. Phys. C_ 41, 045104.  10.1088/1674-1137/41/4/045104
  * López-Coto et al. (2022) López-Coto, R., de Oña Wilhelmi, E., Aharonian, F., Amato, E., and Hinton, J. (2022).  Gamma-ray haloes around pulsars as the key to understanding cosmic-ray transport in the Galaxy.  _Nature Astron._ 6, 199–206.  10.1038/s41550-021-01580-0
  * Ma et al. (2022) Ma, X.-H. et al. (2022).  Chapter 1 LHAASO Instruments and Detector technology *.  _Chin. Phys. C_ 46, 030001.  10.1088/1674-1137/ac3fa6
  * Malkov et al. (2013) Malkov, M. A., Diamond, P. H., Sagdeev, R. Z., Aharonian, F. A., and Moskalenko, I. V. (2013).  Analytic Solution for Self-regulated Collective Escape of Cosmic Rays from their Acceleration Sites.  _Astrophys. J._ 768, 73.  10.1088/0004-637X/768/1/73
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., and Hobbs, M. (2005).  The Australia Telescope National Facility pulsar catalogue.  _Astron. J._ 129, 1993.  10.1086/428488
  * Manconi et al. (2020) Manconi, S., Di Mauro, M., and Donato, F. (2020).  Contribution of pulsars to cosmic-ray positrons in light of recent observation of inverse-Compton halos.  _Phys. Rev. D_ 102, 023015.  10.1103/PhysRevD.102.023015
  * Martin et al. (2022) Martin, P., Marcowith, A., and Tibaldo, L. (2022).  Are pulsar halos rare ?  _arXiv e-prints_ , arXiv:2206.11803
  * Mitchell et al. (2021) Mitchell, A. M. W., Caroff, S., Hinton, J., and Mohrmann, L. (2021).  Detection of extended TeV emission around the Geminga pulsar with H.E.S.S.  _PoS_ ICRC2021, 780
  * Moskalenko and Strong (1998) Moskalenko, I. V. and Strong, A. W. (1998).  Production and propagation of cosmic ray positrons and electrons.  _Astrophys. J._ 493, 694–707.  10.1086/305152
  * Mukhopadhyay and Linden (2022) Mukhopadhyay, P. and Linden, T. (2022).  Self-generated cosmic-ray turbulence can explain the morphology of TeV halos.  _Phys. Rev. D_ 105, 123008.  10.1103/PhysRevD.105.123008
  * Parkinson et al. (2010) Parkinson, P. M. S. et al. (2010).  Eight gamma-ray pulsars discovered in blind frequency searches of Fermi LAT data.  _Astrophys. J._ 725, 571–584.  10.1088/0004-637X/725/1/571
  * Perri et al. (2016) Perri, S., Amato, E., and Zimbardo, G. (2016).  Transport of relativistic electrons at shocks in shell-type supernova remnants: diffusive and superdiffusive regimes.  _Astron. Astrophys._ 596, A34.  10.1051/0004-6361/201628767
  * Perri and Zimbardo (2009a) Perri, S. and Zimbardo, G. (2009a).  Ion and electron superdiffusive transport in the interplanetary space.  _Advances in Space Research_ 44, 465–470.  10.1016/j.asr.2009.04.017
  * Perri and Zimbardo (2009b) Perri, S. and Zimbardo, G. (2009b).  Ion Superdiffusion at the Solar Wind Termination Shock.  _Astrophys. J. Lett._ 693, L118–L121.  10.1088/0004-637X/693/2/L118
  * Plucinsky et al. (1996) Plucinsky, P. P., Snowden, S. L., Aschenbach, B., Egger, R., Edgar, R. J., and McCammon, D. (1996).  ROSAT Survey Observations of the Monogem Ring.  _Astrophys. J._ 463, 224.  10.1086/177236
  * Posselt et al. (2017) Posselt, B., Pavlov, G., Slane, P., Romani, R., Bucciantini, N., Bykov, A., et al. (2017).  Geminga’s puzzling pulsar wind nebula.  _Astrophys. J._ 835, 66.  10.3847/1538-4357/835/1/66
  * Prodanovic et al. (2007) Prodanovic, T., Fields, B. D., and Beacom, J. F. (2007).  Diffuse gamma rays from the galactic plane: probing the gev excess and identifying the TeV excess.  _Astropart. Phys._ 27, 10–20.  10.1016/j.astropartphys.2006.08.007
  * Profumo et al. (2018) Profumo, S., Reynoso-Cordova, J., Kaaz, N., and Silverman, M. (2018).  Lessons from HAWC pulsar wind nebulae observations: The diffusion constant is not a constant; pulsars remain the likeliest sources of the anomalous positron fraction; cosmic rays are trapped for long periods of time in pockets of inefficient diffusion.  _Phys. Rev. D_ 97, 123008.  10.1103/PhysRevD.97.123008
  * Ptuskin et al. (2008) Ptuskin, V. S., Zirakashvili, V. N., and Plesser, A. A. (2008).  Non-linear diffusion of cosmic rays.  _Adv. Space Res._ 42, 486–490.  10.1016/j.asr.2007.12.007
  * Qiao et al. (2022) Qiao, B.-Q., Liu, W., Zhao, M.-J., Bi, X.-J., and Guo, Y.-Q. (2022).  Galactic cosmic ray propagation: sub-PeV diffuse gamma-ray and neutrino emission.  _Front. Phys. (Beijing)_ 17, 64501.  10.1007/s11467-022-1188-8
  * Ragot and Kirk (1997) Ragot, B. R. and Kirk, J. G. (1997).  Anomalous transport of cosmic ray electrons.  _Astron. Astrophys._ 327, 432–440
  * Recchia et al. (2016) Recchia, S., Blasi, P., and Morlino, G. (2016).  On the radial distribution of Galactic cosmic rays.  _Mon. Not. Roy. Astron. Soc._ 462, L88–L92.  10.1093/mnrasl/slw136
  * Recchia et al. (2021) Recchia, S., Di Mauro, M., Aharonian, F. A., Orusa, L., Donato, F., Gabici, S., et al. (2021).  Do the Geminga, Monogem and PSR J0622+3749 $\gamma$-ray halos imply slow diffusion around pulsars?  _Phys. Rev. D_ 104, 123017.  10.1103/PhysRevD.104.123017
  * Recchia et al. (2022) Recchia, S., Galli, D., Nava, L., Padovani, M., Gabici, S., Marcowith, A., et al. (2022).  Grammage of cosmic rays in the proximity of supernova remnants embedded in a partially ionized medium.  _Astron. Astrophys._ 660, A57.  10.1051/0004-6361/202142558
  * Reynolds et al. (2017) Reynolds, S. P., Pavlov, G. G., Kargaltsev, O., Klingler, N., Renaud, M., and Mereghetti, S. (2017).  Pulsar-Wind Nebulae and Magnetar Outflows: Observations at Radio, X-Ray, and Gamma-Ray Wavelengths.  _Space Sci. Rev._ 207, 175–234.  10.1007/s11214-017-0356-6
  * Riviere et al. (2017) Riviere, C., Fleischhack, H., and Sandoval, A. (2017).  HAWC detection of TeV emission near PSR B0540+23.  _The Astronomer’s Telegram_ 10941, 1
  * Schroer et al. (2021) Schroer, B., Pezzi, O., Caprioli, D., Haggerty, C., and Blasi, P. (2021).  Dynamical Effects of Cosmic Rays on the Medium Surrounding Their Sources.  _Astrophys. J. Lett._ 914, L13.  10.3847/2041-8213/ac02cd
  * Schroer et al. (2022) Schroer, B., Pezzi, O., Caprioli, D., Haggerty, C., and Blasi, P. (2022).  Cosmic-ray generated bubbles around their sources.  _Mon. Not. Roy. Astron. Soc._ 512, 233–244.  10.1093/mnras/stac466
  * Sheidaei et al. (2011) Sheidaei, F., Djannati-Ataï, A., and Gast, H. (2011).  Discovery of very-high-energy $\gamma$-ray emission from the vicinity of PSR J1831-952 with H.E.S.S.  In _3rd International Fermi Symposium_. vol. 7, 244–247.  10.7529/ICRC2011/V07/1225
  * Skilling (1971) Skilling, J. (1971).  Cosmic Rays in the Galaxy: Convection or Diffusion?  _Astrophys. J._ 170, 265.  10.1086/151210
  * Strong et al. (2007) Strong, A. W., Moskalenko, I. V., and Ptuskin, V. S. (2007).  Cosmic-ray propagation and interactions in the Galaxy.  _Ann. Rev. Nucl. Part. Sci._ 57, 285–327.  10.1146/annurev.nucl.57.090506.123011
  * Sudoh et al. (2019) Sudoh, T., Linden, T., and Beacom, J. F. (2019).  TeV Halos are Everywhere: Prospects for New Discoveries.  _Phys. Rev. D_ 100, 043016.  10.1103/PhysRevD.100.043016
  * Tang and Piran (2019) Tang, X. and Piran, T. (2019).  Positron flux and $\gamma$-ray emission from Geminga pulsar and pulsar wind nebula.  _Mon. Not. Roy. Astron. Soc._ 484, 3491–3501.  10.1093/mnras/stz268
  * Tomassetti (2012) Tomassetti, N. (2012).  Origin of the Cosmic-Ray Spectral Hardening.  _Astrophys. J. Lett._ 752, L13.  10.1088/2041-8205/752/1/L13
  * Torres Escobedo et al. (2021) Torres Escobedo, R. et al. (2021).  Follow-up Analysis to Geminga’s contribution to the Local Positron Excess with HAWC Gamma-Ray Observatory.  _PoS_ ICRC2021, 842.  10.22323/1.395.0842
  * Uchaikin and Sibatov (2017) Uchaikin, V. and Sibatov, R. (2017).  Calculations of CR energy spectra within the NoRD model.  _arXiv e-prints_ , arXiv:1703.06486
  * Uchaikin and Sibatov (2012) Uchaikin, V. V. and Sibatov, R. T. (2012).  On fractional differential models for cosmic ray diffusion.  _Gravitation and Cosmology_ 18, 122–126.  10.1134/S0202289312020132
  * Uchaikin and Zolotarev (1999) Uchaikin, V. V. and Zolotarev, V. M. (1999).  _Chance and Stability_ (Utrecht, Netherland: VSP)
  * Van Etten et al. (2012) Van Etten, A., Romani, R. W., and Ng, C. Y. (2012).  A Chandra Proper Motion for PSR J1809-2332.  _Astrophys. J._ 755, 151.  10.1088/0004-637X/755/2/151
  * Volkov and Lagutin (2013) Volkov, N. and Lagutin, A. (2013).  Cosmic ray positron to electron ratio in the Galaxy: results of the fractional diffusion approach.  In _Journal of Physics Conference Series_. vol. 409 of _Journal of Physics Conference Series_ , 012042.  10.1088/1742-6596/409/1/012042
  * Wang et al. (2021) Wang, S.-H., Fang, K., Bi, X.-J., and Yin, P.-F. (2021).  Test of the superdiffusion model in the interstellar medium around the Geminga pulsar.  _Phys. Rev. D_ 103, 063035.  10.1103/PhysRevD.103.063035
  * Xi et al. (2019) Xi, S.-Q., Liu, R.-Y., Huang, Z.-Q., Fang, K., and Wang, X.-Y. (2019).  GeV observations of the extended pulsar wind nebulae constrain the pulsar interpretations of the cosmic-ray positron excess.  _Astrophys. J._ 878, 104.  10.3847/1538-4357/ab20c9
  * Yan and Lazarian (2008) Yan, H. and Lazarian, A. (2008).  Cosmic Ray Propagation: Nonlinear Diffusion Parallel and Perpendicular to Mean Magnetic Field.  _Astrophys. J._ 673, 942.  10.1086/524771
  * Yan et al. (2022) Yan, K., Liu, R.-Y., Chen, S. Z., and Wang, X.-Y. (2022).  Prospect of detecting TeV halos with LHAASO: in the framework of the anisotropic diffusion model.  _arXiv e-prints_ , arXiv:2205.14563
  * Yang and Liu (2022) Yang, R.-z. and Liu, B. (2022).  On the surface brightness radial profile of the extended $\gamma$-ray sources.  _Sci. China Phys. Mech. Astron._ 65, 219511.  10.1007/s11433-021-1777-4
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., and Wang, N. (2017).  A New Electron-density Model for Estimation of Pulsar and FRB Distances.  _Astrophys. J._ 835, 29.  10.3847/1538-4357/835/1/29
  * Yin et al. (2013) Yin, P.-F., Yu, Z.-H., Yuan, Q., and Bi, X.-J. (2013).  Pulsar interpretation for the AMS-02 result.  _Phys. Rev. D_ 88, 023001.  10.1103/PhysRevD.88.023001
  * Yin et al. (2009) Yin, P.-f., Yuan, Q., Liu, J., Zhang, J., Bi, X.-j., and Zhu, S.-h. (2009).  PAMELA data and leptonically decaying dark matter.  _Phys. Rev. D_ 79, 023512.  10.1103/PhysRevD.79.023512
  * Yuan et al. (2015) Yuan, Q., Bi, X.-J., Chen, G.-M., Guo, Y.-Q., Lin, S.-J., and Zhang, X. (2015).  Implications of the AMS-02 positron fraction in cosmic rays.  _Astropart. Phys._ 60, 1–12.  10.1016/j.astropartphys.2014.05.005
  * Yuan et al. (2017) Yuan, Q., Lin, S.-J., Fang, K., and Bi, X.-J. (2017).  Propagation of cosmic rays in the AMS-02 era.  _Phys. Rev. D_ 95, 083007.  10.1103/PhysRevD.95.083007
  * Yuan et al. (2021) Yuan, Q., Qiao, B.-Q., Guo, Y.-Q., Fan, Y.-Z., and Bi, X.-J. (2021).  Nearby source interpretation of differences among light and medium composition spectra in cosmic rays.  _Front. Phys. (Beijing)_ 16, 24501.  10.1007/s11467-020-0990-4
  * Yue et al. (2020) Yue, C. et al. (2020).  Implications on the origin of cosmic rays in light of 10 TV spectral softenings.  _Front. Phys. (Beijing)_ 15, 24601.  10.1007/s11467-019-0946-8
  * Yuksel et al. (2009) Yuksel, H., Kistler, M. D., and Stanev, T. (2009).  TeV Gamma Rays from Geminga and the Origin of the GeV Positron Excess.  _Phys. Rev. Lett._ 103, 051101.  10.1103/PhysRevLett.103.051101
  * Zhang et al. (2014) Zhang, S. N. et al. (2014).  The high energy cosmic-radiation detection (HERD) facility onboard China’s Space Station.  _Proc. SPIE Int. Soc. Opt. Eng._ 9144, 91440X.  10.1117/12.2055280
  * Zhang et al. (2021) Zhang, Y., Liu, R.-Y., Chen, S. Z., and Wang, X.-Y. (2021).  Morphology of Gamma-Ray Halos around Middle-aged Pulsars: Influence of the Pulsar Proper Motion.  _Astrophys. J._ 922, 130.  10.3847/1538-4357/ac235e
  * Zhao et al. (2021) Zhao, M.-J., Fang, K., and Bi, X.-J. (2021).  Constraints on the spatially dependent cosmic-ray propagation model from Bayesian analysis.  _Phys. Rev. D_ 104, 123001.  10.1103/PhysRevD.104.123001
  * Zhou et al. (2022) Zhou, G.-Y., Yu, Z.-H., Yuan, Q., and Zhang, H.-H. (2022).  Geminga contribution to the cosmic-ray positron excess according to the gamma-ray observations.  _arXiv e-prints_ , arXiv:2205.07038