# Mirage Sources and Large TeV Halo-Pulsar Offsets: Exploring the Parameter
Space

Yiwei Bao1,2 Ruo-Yu Liu1,3<EMAIL_ADDRESS>Gwenael Giacinti2,4
<EMAIL_ADDRESS>Hai-Ming Zhang1 Yang Chen1,3<EMAIL_ADDRESS>1Department of Astronomy, Nanjing University, 163 Xianlin Avenue, Nanjing
210023, China
2Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 201210,
China
3Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University,
Ministry of Education, Nanjing, China
4School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China

###### Abstract

We investigate the asymmetric propagation of 100 TeV electrons (whose
radiation mainly concentrates on 20–30 TeV) in turbulent magnetic fields
around pulsars, using GPU-accelerated simulations to explore their
trajectories and interactions within pulsar wind nebulae and the interstellar
medium. Key results include the identification of “mirage” sources indicating
significant offsets in high-energy emissions from their originating pulsars,
challenging the results of traditional symmetric diffusion models. By varying
parameters like source distance, magnetic field strength, and electron
injection spectral index, the study delineates their effects on observable
phenomena such as the probability that a source has at least one mirage around
it, as well as the source separation. Our results offer insights into some
puzzling sources observed recently by the Large High Altitude Air Shower
Observatory (LHAASO), and shed light on the cosmic-ray transport mechanism in
the interstellar medium.

††preprint: arXiv:0000.00000

##  I Introduction

Pulsars are rapidly rotating neutron stars. They power relativistic pulsar
winds composed of $e^{\pm}$ pairs [1, 2, 3] at the expense of their rotational
energies. A termination shock is formed where the pulsar wind encounters the
ambient medium and reaches a pressure balance. Either at the termination shock
[4], or far upstream of the termination shock [5], magnetic reconnection may
occur and some of the $e^{\pm}$ pairs be accelerated up to $\sim$ PeV [3].
Alternatively, these particles may be accelerated at the termination shock,
especially in its equatorial region [6, 7]. The radiation of these
relativistic pairs in the shocked wind gives rise to the pulsar wind nebula
(PWN).

Recent observations of extended TeV emissions around middle-aged pulsars by
HAWC [8] and LHAASO [9] have sparked significant interest among the community.
It can be inferred from the surface brightness profiles of these extended
$\gamma$-ray emission that $e^{\pm}$ pairs are escaping their PWNe and
diffusing in the ambient interstellar medium (ISM) with a significantly
suppressed diffusion coefficient, by a factor of $\sim 100-1000$, compared to
the Galactic mean value. It may be interpreted as an enhancement of the
interstellar turbulence [10, 11, 12], or a reasonably large level of pre-
existing turbulence [13], or ascribed to the cross-field transport of
particles if the local magnetic field and the observer’s line of sight is
approximately aligned [14]. The size of the region with suppressed particle
diffusion is suggested to be at least 20 pc [15]. The extended $\gamma$-ray
emission region beyond the PWN is usually referred to as “TeV halos”.

In principle, each middle-aged pulsar could power a TeV halo [16, 17]. Indeed,
we see many sources identified by LHAASO [18] and HAWC [19] are associated
with middle-aged pulsars. However, many of them present significant spatial
offsets between positions of the pulsars and the centroid of the TeV sources,
which is not expected if pairs are injected from the position of the pulsar
(or from the compact PWN around the pulsar) and isotropically diffuse in the
surrounding ISM [20]. Actually, isotropic diffusion is a good approximation of
the real transport of particles only if the physical scale considered in the
problem is larger than multiple coherence lengths of the turbulent magnetic
field [21] (about 1/5 of the outer scale of the turbulence for Kolmogorov
spectra), such as e.g. in the modelling of the global transport of cosmic rays
(CRs) in the Galaxy.

However, if we focus on the particle transport around accelerators, especially
those located above/below the Galactic disk or in the outer Galaxy region, the
isotropic diffusion approximation could be problematic, given that the average
coherent length of the ISM is of the order of $\sim 100$ pc [22]. This is
particularly important for pulsar halo candidates detected by LHAASO-KM2A,
because these sources are due to the emission from electrons/positrons with
energies above 100 TeV. Given the rapid cooling rate of these extremely
energetic leptons, they may probably cool down before crossing one magnetic
coherence length. Therefore, the distribution of high-energy particles could
be highly asymmetric, following the direction of the local magnetic field. In
the companion paper [23] (Paper I), we have shown that the asymmetric
propagation of electrons/positrons around their accelerators may cause
multiple “mirage” sources and significant spatial offset between the
accelerator and the source. The main goal of the present paper is to give a
more quantitative description of our method, and to explore extensively the
influence of different model parameters on the obtained source morphology.

The rest of the paper is organized as follows: in Section II, we describe our
method in detail; we explore the properties of the simulated sources with
different parameters in Section III; In Section IV, we discuss several
possible applications of our findings, and we conclude in Section V.

##  II Method Description

To characterise the asymmetric propagation of particles, we conduct test-
particle simulations by calculating the trajectories of 100 TeV electrons, the
radiation spectrum of which typically peaks at approximately 20 TeV via the
inverse Compton (IC) scattering off the cosmic microwave background (CMB). We
employ a first-principle approach by injecting electrons into a turbulent
magnetic field and determining their motion based on the Lorentz force [24,
25, 26]. We employ GPU-acceleration methods, allowing the code to operate in a
highly parallelized manner. To preserve GPU on-board memory, we calculate the
magnetic field anew at every time step for every particles’ position. While
this approach is time-intensive, it offers a greater accuracy compared to the
method which specifies fields on a discrete spatial lattice in advance and
interpolates while computing. Here, we thoroughly explain in detail the
process of conducting the simulation in terms of particle motion and field
generation.

###  II.1 Turbulent magnetic fields

We follow the same method in Ref.[24] to generate the turbulent magnetic
field. The field is written as a superposition of $N_{m}$ plane waves with

$\bm{B_{\rm
t}(r)}=\sum_{n=1}^{N_{m}}A_{n}\epsilon_{n}\exp({i\bm{k}_{n}\cdot\bm{r}+i\beta_{n}}),$
(1)

where the explanation of $A_{n}$, $\epsilon_{n}$, $\bm{k}$ and $\beta_{n}$ can
be found in [24]. In addition, in the formula

$A^{2}(k_{n})=\sigma^{2}G(k_{n})\left[\sum\limits_{n=1}^{N_{m}}G(k_{n})\right]^{-1},$
(2)

we use $G(k_{n})\propto k^{-5/3}\Delta k$ instead of the one used in [24].
Technically, a value of 100 waves per decade [see e.g., 25] is used to balance
the speed and accuracy. In practice, we generate the random constants on CPU
which determine the turbulent magnetic field and store them in the “constant
memory” of GPU to accelerate the computation. There are 600 waves in total,
which covers a range of 6 magnitudes. We have calculated 100 distinct
configurations of the turbulent magnetic field and show a subset of these
configurations, which exhibit intriguing properties. $T_{\rm age}$, the age of
the pulsar, is adopted to be 110 kyr.

###  II.2 Energy losses

For the synchrotron loss (in terms of $\gamma$), we utilize

$\dot{\gamma}_{\rm sync}=-3.25\times 10^{-8}\frac{B^{2}}{8\pi}\gamma^{2},$ (3)

where $B$ is the magnetic field strength (here we use $B^{2}=B_{\rm
r}^{2}+B_{\rm t}^{2}$, where $B_{\rm r}$ is the regular component and $B_{\rm
t}$ the root mean square of the turbulent component). For IC losses we use the
approximation given by Fang et al. [27]:

$\dot{\gamma}_{\rm IC}=-\frac{20c\sigma_{\rm
T}\omega\gamma^{2}}{\pi^{4}m_{e}c^{2}}Y(\gamma,T),$ (4)

where $\omega$ is the grey-body energy density, $T$ the grey-body temperature
and $Y(\gamma,T)$ the parameterised approximation function [27]. For the IC
seed photons, we use the interstellar radiation field averaged over a 2 kpc
region surrounding the Sun as a reference, which is composed of five gray-body
components with temperatures of 23209 K, 6150.4 K, 3249.3 K, 313.3 K, and 33.1
K [27]. These components have corresponding energy densities of 0.12, 0.23,
0.37, 0.055, and 0.25 eV cm-3 [27], respectively.

###  II.3 Particle motions

We employ the original Boris pusher which avoids the phase error [see e.g.,
Boris A pusher in 28] to preserve phase space volume. Since the electric field
in the turbulent plasma may be neglected [24], the Boris pusher can be written
as

$\displaystyle\bm{x}^{\prime}:$
$\displaystyle=\bm{x}_{n-1}+\frac{1}{2}\bm{v}_{n-1}\cdot\Delta t$ (5)
$\displaystyle\bm{t}:$ $\displaystyle=\tan\frac{\theta}{2}\hat{\bm{b}}$ (6)
$\displaystyle\gamma\bm{v}^{\prime}:$
$\displaystyle=\gamma\bm{v}_{n-1}+\gamma\bm{v}_{n-1}\times\bm{t}$ (7)
$\displaystyle\gamma\bm{v}_{n}$
$\displaystyle=\gamma\bm{v}^{\prime}+\frac{2}{1+|\bm{t}|^{2}}\left(\gamma\bm{v}^{\prime}\times\bm{t}\right)$
(8) $\displaystyle\bm{x}_{n}$
$\displaystyle=\bm{x}^{\prime}+\frac{1}{2}\bm{v}_{n}\cdot\Delta t$ (9)

where $\hat{\bm{b}}:=\bm{B}/|\bm{B}|$ is the unit vector, $\gamma$ the
electron Lorenrz factor, $n$ the time step, $\bm{x}^{\prime}$, $\bm{t}$,
$\bm{v}^{\prime}$ are three intermediate variables. A time step of 1/100
gyration period is utilized in the simulation. We also take synchrotron/IC
losses into consideration. The spin-down luminosity of the pulsar is
$L(t)=L_{0}/[1+(t/\tau_{0})^{2}]$, where $L_{0}$ is the initial spin-down
luminosity and $\tau_{0}$ the initial spin-down timescale. Electrons are
injected in a power-law form as ${\rm d}N/{\rm
d}\gamma\propto\gamma^{-\alpha}$, where $\alpha$ the power-law index. To mimic
the continuous injection of electrons, we inject $\approx 10^{4}$ numerical
particles continuously over time with energy $\gamma(t_{\rm inj})=\int^{t_{\rm
inj}}_{t_{\rm age}}\left[\dot{\gamma}_{\rm sync}(\gamma)+\dot{\gamma}_{\rm
IC}(\gamma)\right]+{\rm 100\,TeV}/{m_{e}c^{2}}$ so that all of them are cooled
down to 100 TeV at time $t_{\rm age}$, where $t_{\rm inj}$ is the injection
time of the particle and $t_{\rm age}$ the current age of pulsar. The value of
$t_{\rm inj}$ is set such that every numerical particle represents roughly the
same number of real electrons.

Figure 1: Four sets of mirages sources (with distance at 250 pc). The color
shows the relative intensity. The electrons are injected at (0,0), and the
black, green and orange dots show the detected source locations.

Figure 2: Four sets of large offsets (with distance at 1 kpc). The color shows
the relative intensity. The black dots show the electron injection site
locations, and the green dots correspond to the detected source locations.

##  III Simulation data analysis

###  III.1 Methods

Figure 3:  Dependence of $\xi$ on other parameters.

Figure 4: Scatter plot of the flux ratio $\chi$ and separation distance for
different $L_{\rm c}$ (top panel), magnetic field strength $B^{2}$ (middle
panel) and the strength ratio between regular field and turbulent field
$B_{r}/B_{t}$ (bottom panel). The distance is fixed at 1 kpc. The dashed line
in the histogram mark the average value of the separation.

To mimic observations, we envisage that the halo is measured by LHAASO-KM2A,
the point spread function (PSF) of which can be approximated as a 2D Gaussian
function with the standard deviation $\sigma=0.3^{\circ}$ around 20 TeV [29].
To present the simulation result more conveniently, we consider the
measurement of the pulsar halo solely at 20 TeV and set a reference luminosity
of the pulsar halo to be $L_{\rm ref}=6\times 10^{32}\,$erg/s at this energy.
This luminosity is about 6 times that of the observed luminosity of Geminga’s
halo in the range of $8-40$ TeV [8], which is not unreasonable. We choose this
high value in order to make the relevant features of simulated pulsar halos
more significant. According to LHAASO’s effective area, the reference
luminosity corresponds to a total event rate of $1.6\times 10^{6}(250\,{\rm
pc}/d)^{2}$ counts per year for our reference scenario, where $d$ is the
distance of the pulsar to us. We set our region of interest (ROI) to a region
with projected size of 100 pc $\times$ 100 pc. For the case of $L_{\rm c}=400$
pc or $d=6$ kpc, we use a ROI of 400 pc $\times$ 400 pc. In any case, the ROI
is divided into many elementary bins of $0.1^{\circ}\times 0.1^{\circ}$ or
correspondingly $0.436\,{\rm pc}(d/{\rm 250\,pc})\times 0.436\,{\rm pc}(d/{\rm
250\,pc})$. We then can calculate the expected event count number in each bin
based on the simulated intensity distribution of the pulsar halo.

To mimic the background, we generate a Poisson noise map (noise A) based on
the background event rate measured from the Crab Nebula region, which is 1
event per hour within a 1∘ cone [3], or 0.318 count/hour in each bin. We then
add the count number of the noise to the event number in each bin so that we
obtain a simulated counts map. We then generate the model, composed of another
Poisson noise (noise B) based on the measured background rate and a
symmetrical 2D Gaussian template (with the normalisation, extension, and
positron free) to fit the obtained counts map, and calculate the test
statistic (TS) value (TS1) [8, 29]. If the TS value exceeds 25, we deem that
we identify one source. Next, we set a model composed of noise B plus two
symmetrical 2D Gaussian template to fit the obtained counts map again, and
calculate the TS value (TS2). If TS2 exceeds TS1 by 25, we deem that we
identify two sources. The chain continues by adding another hypothetical
source until the increment of TS value is less than 25. Practically, we would
identify one or two sources in most cases, and in a few cases more mirage
sources may be also found in the identification process. In Figure 1 and
Figure 2, we show four configurations of mirage sources and four of large
offsets, respectively, for reference. In all panels in these two Figures, the
electrons are injected at (0,0), in the center. The colors show the level of
the relative gamma-ray intensity as calculated from our simulations. In Figure
1, the black, green and orange dots correspond to the locations of the sources
we detect. One can see that, in each panel, two more sources are detected far
from the actual astrophysical source. These are mirage sources. In Figure 2,
the green dots display the locations of the detected sources. In each panel, a
considerable offset is clearly visible.

The formation of mirage sources is related to the magnetic field geometry.
After escaping accelerators, electrons spiral around extensive magnetic field
lines in the surrounding medium. They are scattered by irregularities in the
magnetic field but their trajectories roughly follow the mean direction of the
magnetic field over one coherence length $L_{\rm c}$. By definition, the mean
field direction changes dramatically in another coherence length. If the field
line within certain coherence length roughly orients with our line of sight, a
visual projection effect leads to the apparent enhancement of electron column
density in that direction, which is not necessarily toward the injection
point, creating a mirage source or an offset. Note that the turbulent magnetic
field is treated as the superposition of many plane waves. Therefore, apart
from the main field direction change at the scale of $\sim L_{\rm c}$, there
are also smaller bents in the magnetic field within one coherence length. Even
if particles do not cross over one coherence length before cooling, they may
still form mirage sources or offset due to these small bents.

###  III.2 Mirage Ratio

Because magnetic fields are generated randomly, the resulting morphology of
the pulsar halo and properties of the emission in one run can be quite
different from another. Therefore, we need to study the problem in a
statistical way. One relevant issue is how likely we would identify multiple
sources around an accelerator (i.e., pulsar/pulsar wind nebula). This can be
assessed with what we will call the “mirage ratio”, defined here as the
fraction of all runs with more than one sources identified in the mock data.
We run the simulation with 100 randomly generated magnetic field
configurations, and then measure the mirage ratio, denoted by $\xi$. We also
vary parameters such as pulsar distance $d$, injection power-law indices
$\alpha$, $B_{\rm r}/B_{\rm t}$ and $B^{2}$ one by one to see their
influences. The luminosity of the pulsar halo is always fixed at $L_{\rm
ref}$. The resulting values of $\xi$ are displayed in the four panels of
Figure 3.

We see that mirage sources are commonly expected for nearby sources, i.e.,
$d<1\,$kpc. For more distant sources, the value of $\xi$ declines
significantly with increasing distance $d$. The decline can be attributed to
the blurring effect of the PSF, and to the reduction in the signal-to-noise
ratio at greater distances. Above 4 kpc, we almost do not expect any mirage
source around the accelerator. In addition to such a general trend, we see
that the mirage ratio between $1.5-3.5$ kpc also depends on various model
parameters.

We explore the influence of the coherence length in the top-left panel. As
$L_{\rm c}$ increase, the value of $\xi$ increases significantly. The
underlying reason is that a longer coherence length leads to a greater
separation between two mirage sources and hence less influence by the PSF, as
will be explained subsequently. If the coherence length is small, electrons
can cross multiple coherence magnetic fields before cooling, so their
distribution becomes largely isotropic and hence do not produce an enhancement
of electron column density toward certain directions. The injection spectral
index also affects the mirage ratio but not as strongly as $L_{\rm c}$ does,
as shown in the top-right panel. $\xi$ rises for a steeper spectrum (larger
$\alpha$). This can be explained as follows: particles with larger energies
propagate more isotropically [21], and a smaller $\alpha$ implies that more
particles were injected earlier (and cooled down from a higher energy), and
therefore propagate more isotropically. $\xi$ also increases slightly as
$B_{\rm r}/B_{\rm t}$ rises. This is intuitive because when $B_{\rm r}/B_{\rm
t}$ is larger (and the total magnetic field energy remains the same), less
energy is deposited in the turbulent energy, and the transport is closer to
that in the anisotropic diffusion scenario [14], preferentially resulting in
an elongated source with a gradually varying 2D intensity map. Meanwhile,
$\xi$ decreases as $B^{2}$ increases, because the distance on which particles
are transported depends on their energy losses, which are more severe when the
magnetic field is stronger, and it then becomes harder for the separation of
the sources to overcome the PSF. Comparing the cases of $B^{2}=9\mu{\rm
G}^{2}$ and $81\mu{\rm G}^{2}$, the impact of $B^{2}$ is relatively weak
because when the magnetic field is weak, the energy losses are dominated by
the IC losses, the energy density of which is equivalently a magnetic field of
$\sim 3.6\ \mu$G.

###  III.3 Source separation and flux ratio

In this subsection, we discuss the separation distance (defined as the maximal
separation distance between any two sources if two or more mirage sources are
identified) and the flux ratio ($\chi$, defined as main source/mirage source,
where the main source is the source with the smallest distance to the
injection point and the mirage source is the farthest one) of the sources. For
each chosen set of parameters, we carry out 100 runs and calculate the above
two quantities of each run, and show their distribution in the plane of
separation value and flux ratio. In all runs, the accelerator’s distance is
fixed at $d=1\,$kpc. To reveal the difference more clearly, we set here a
higher criterion TS$=25$ for the identification of a source. As such, the
mirage ratio will decrease compared to the case shown in Figure 3 with the
same parameters (i.e., not all runs can generate a mirage source at
$d=1\,$kpc).

Our results are shown in Figure 4. In the top panel of the figure, we can see
that the mean separation between the main source and the mirage source is
positively correlated with $L_{\rm c}$. Given that the mirage source is formed
due to the projection effect when particles propagate to another coherence of
the magnetic field as discussed before, the separation value is correlated
with the coherence length. On the other hand, we do not see a clear dependence
of the flux ratio $\xi$ on $L_{\rm c}$.

The influence of the magnetic field strength is noticeable on both the
separation value and the flux ratio, as shown in the middle panel of Figure 4.
The effect can be attributed to faster cooling of electrons in stronger
magnetic fields. Injected electrons may not cross one coherence length before
cooling if the magnetic field strength is strong, and therefore only reveal
the small bents within one coherence length of the magnetic field, leading to
a smaller separation between the main source and the mirage source. Due to the
same reason, the electron number consisting of the mirage source is also
smaller for a stronger magnetic field, leading to a lower flux ratio. From the
middle panel of Figure 4, we see that the flux of the mirage source may be
larger than the main source in a considerable fraction of runs with $B^{2}=1$
and 9$\,\mu\rm G^{2}$. For $B^{2}=81\,\mu\rm G^{2}$, the flux ratio is greater
than unity in most of runs.

The regular-to-turbulent ratio $B_{\rm r}/B_{\rm t}$ also affects the source
separation, although the impact is less important than the previous two
parameters. We see from the bottom panel of Figure 4 that, when $B_{\rm
r}/B_{\rm t}$ decreases, the separation distance (statistically) increases
slightly. Adding a regular magnetic field is equivalent to adding a wave of
infinitely long wavelength. As this regular component becomes stronger, the
coherence becomes less distinct. Correspondingly, there will be less
significant change of the mean field direction. In this sense, it is similar
to the influence of the coherence length. Indeed, in the extreme case of
$B_{\rm r}/B_{\rm t}\to\infty$, electrons just spiral around the regular
magnetic field, so no mirage source and subsequently no source separation is
expected.

###  III.4 Offsets

Figure 5: Scatter plot of offsets and total counts of the identified source
which is closest to the accelerator for different $L_{c}$ (top panels),
$B^{2}$ (middle panels), and $B_{r}/B_{t}$ (bottom panels). The distance is
fixed at 6 kpc (the pixel size 0.1∘ corresponding to $\approx 10$ pc; the PSF
( 0.3∘) corresponding to $\approx 30$ pc). In the left columns, the luminosity
is $\mathcal{L}=\mathcal{L}_{\rm ref}$, while in the right columns, the
luminosity is $\mathcal{L}=10\mathcal{L}_{\rm ref}$. Dashed lines in the
histogram mark the average value of the offset.

Now we investigate the spatial offsets of the identified sources from the
accelerator. We locate the source at a relatively large distance, 6 kpc, to
reduce the mirage ratio, so that the definition of the offset is clear.
Nevertheless, if more than one sources are identified, the offset is defined
as the distance between the accelerator and the source closest to it.

For $L=L_{\rm ref}$, as shown in the left columns of Figure 5, there are only
a few points sparsely distributed in the plane of the offset and the total
counts of the source (as a representative of the source flux). This is owing
to the large distance we chose, which results in a low flux of the source and
hence not many of them are identified. For those with a source being
identified, the obtained offsets are quite large for $L=L_{\rm ref}$. This is
because when the source is faint and the distance is far, the position of the
Gaussian centre is subject to the noise fluctuation and large PSF. If we
assume a higher lumionsity, i.e., $L=10L_{\rm ref}$, as shown in the right
column of Figure 5, the offsets become generally smaller. Meanwhile, as can be
expected, the offset increases with larger $L_{c}$ and is inversely correlated
to $B^{2}$. The reasons are similar to that for the separation distance, as
discussed in the last section.

The $\gamma$-ray luminosity $\mathcal{L}$ of the source has an impact on
source identification. This can be seen in the count maps displayed in Figure
6 and in Figure 7. The former shows a configuration of the turbulence with
$L_{\rm c}=400$ pc, and the latter another configuration with $L_{\rm c}=4$
pc. In the configuration with $L_{\rm c}=400$ pc, the emission is found to be
identified as two sources when the luminosity is equal to $10\mathcal{L}_{\rm
ref}$ (see the left panel in Figure 6), while it can only be identified as one
source when it is equal to $\mathcal{L}_{\rm ref}$ (see the right panel). In
the configuration with $L_{\rm c}=4$ pc, the centre of the emission is not
affected by the noise when the luminosity is set to $10\mathcal{L}_{\rm ref}$
(see the left panel in Figure 7), while it is strongly affected by the noise
when it is set to $\mathcal{L}_{\rm ref}$ (see the right panel).

Figure 6: The counts map (relative intensity) for one configuration with
$L_{\rm c}=400$ pc and distance of 6 kpc. The emission can be identified as
two sources if the $\gamma$-ray luminosity is $\mathcal{L}=10\mathcal{L}_{\rm
ref}$ (left panel), and it can only be identified as one source when
$\mathcal{L}=\mathcal{L}_{\rm ref}$ (right panel). The noise has been removed.

Figure 7: The counts map (absolute counts number) for one configuration with $L_{\rm c}=4$ pc. The centre of the emission is not affected by the noise when $\mathcal{L}=10\mathcal{L}_{\rm ref}$ (left panel), and it is severely affected by the noise when $\mathcal{L}=\mathcal{L}_{\rm ref}$ (right panel). The noise is kept.  Table 1: Possible applicable sources Source name | PSR name | Gl | Gb | Distance (kpc) |  $B_{\rm r}$ ($\mu$G) |  $B_{\rm t}$ ($\mu$G, $\beta=0.27$) | Offset (pc) | Offset (∘) | $\zeta$
---|---|---|---|---|---|---|---|---|---
1LHAASO J0216+4237u | J0218+4232 | 139.508 | -17.527 | 3.150 | 0.89 | 0.52 | 18.14 | 0.33 | 41%
1LHAASO J1837-0654u | J1838-0655 | 25.246 | -0.196 | 6.600 | 1.35 | 2.28 | 13.82 | 0.12 | 26%
1LHAASO J1912+1014u | J1913+1011 | 44.485 | -0.167 | 4.613 | 0.62 | 8.65 | 10.47 | 0.13 | 12%
1LHAASO J2005+3415 | J2004+3429 | 71.425 | 1.571 | 10.783 | 0.25 | 0.85 | 47.05 | 0.25 | 84%
1LHAASO J2005+3050 | J2006+3102 | 68.667 | -0.530 | 6.035 | 0.27 | 1.85 | 21.07 | 0.20 | 10%

##  IV Possible applications

In the LHAASO catalog [18], some sources have been discovered in the vicinity
of middle-aged pulsars, but the centres of these sources show significant
offsets from the positions of the pulsars. See in particular those listed in
Table 1. These offsets may be linked with the mirage halo phenomenon discussed
here. We checked this possibility by simulating particle propagation around
these pulsars, based on the Jansson-Farrar model of the Galactic magnetic
field [30, 31]. This model provides the interstellar magnetic field strength,
and the ratio between the regular and the turbulent magnetic field strengths
at the position of each pulsar, which are important parameters in our
simulations as discussed in the previous Section. We note that the turbulence
level in this model needs be re-scaled by a factor $\beta$ to fit the boron-
to-carbon ratio (B/C) [32, 22]. We find that $\beta=0.27$ gives a reasonable
residence time of cosmic rays in our Galaxy by performing the same simulation
as presented in Ref.[32], so that the predicted cosmic-ray grammage 111The
cosmic-ray grammage is given by $\langle X\rangle=1.4m_{\rm p}n_{\rm p}\langle
t_{\rm res}\rangle c$, where $t_{\rm res}$ is the residence time, $n_{\rm p}$
is the hydrogen density in the ISM, $m_{p}$ is the mass of proton, and the
prefactor 1.4 takes into account the contribution of nuclei (such as helium)
in the ISM. is consistent with the cosmic-ray B/C measurement of AMS-02 [34].
In the simulation, we employ $L_{\rm c}=220\,$pc following the observation of
the power spectrum of Galactic foreground rotation-measure maps [22]. For a
smaller coherence length, $\beta$ could be smaller, but still at the level of
$\sim 0.1$, see Ref. [32].

Then, we further assume an injection spectral index $\alpha=2$ for the
electrons, and we simulate electron propagation around each of the pulsars
that are mentioned above. To reduce the influence of the randomly generated
turbulence on the results, we perform 100 separate runs of simulations for
each pulsar and calculate the offset generated in each run. A threshold TS
value of 49 (as is used in the LHAASO catalog), rather than 25, is used for
source identification. The possibility of explaining the offsets with mirage
halos is evaluated with $\zeta$, which is defined as the fraction of the runs
showing a source offset larger than the observed one. The obtained values of
$\zeta$ are shown in the last column of Table 1. We can see that relatively
weak regular magnetic fields are present in the ISM around these pulsars,
which is in favor of producing mirage halos. Indeed, the probability of
showing large offsets are found to be above 10% for all the listed pulsars,
demonstrating that mirage halos may be a reasonable explanation for these
measured large offsets.

##  V Summary and Discussion

Propagation of high-energy cosmic rays in the ISM close to their accelerators
plays an important role in shaping the morphology of extended gamma-ray
sources. In the companion paper [23], we found that the propagation of
electrons is highly asymmetric, and that this may lead to the appearance of
so-called mirage halos, as well as significant spatial offsets between the
positions of the measured gamma-ray sources and that of their parent
accelerators. In the present paper, we explored in detail the formation
mechanism of these mirage halos and spatial offsets, and we investigated
quantitatively the influence of the model parameters on the appearance of
these two phenomena.

We performed test-particle simulations of the propagation of electrons above
100 TeV in the interstellar magnetic fields around their accelerators, taking
radiative cooling into account. With the obtained distributions of electrons,
we then calculated their IC radiation and obtained the mock count maps as
would be seen by LHAASO, as well as the identifiable sources following the
standard likelihood analysis. We found that the mirage sources and the large
spatial offsets between identified sources and particle accelerators can be
ascribed to the asymmetric propagation of particles in the magnetic field and
to projection effects. More specifically, particles propagate preferentially
along the direction of their local magnetic field lines, which changes
dramatically between two coherence lengths. If the field line within a certain
magnetic coherence length is approximately aligned with our line of sight, an
apparent enhancement of the electron column density appears in our direction
due to a projection effect. A mirage source or a source offset is then
created.

Furthermore, based on the properties of identified sources, the influence of
model parameters on the propagation of electrons have been evaluated with
three indicators, namely, the mirage ratio, the maximum separation distance,
and the offset. The mirage ratio is the proportion of runs with identification
of two or more sources, and the maximum separation distance is the maximum
distance between any two sources identified in one run. The offset is the
distance between the pulsar’s position and the centroid of the closest
identified source. Our main findings include:

  * •

The mirage ratio is negatively correlated with $d$, $B_{\rm r}/B_{\rm t}$ and
$B^{2}$, and is positively correlated with $L_{\rm c}$ (strong dependence) and
$\alpha$ (weak dependence). The mirage ratio is almost 100% at 1 kpc with the
reference $\gamma$-ray luminosity. Even at a distance of 3 kpc, the mirage
ratio could still reach more than 10% with reasonable combinations of
parameters.

  * •

The maximum separation is negatively correlated with $B^{2}$, and positively
correlated with $L_{\rm c}$. The large maximum separation suggests that in the
inter-arm region, and in the Galactic halo (where the coherence length is
larger), two sources separated with 100 pc could still be linked with one
another and belong to the same halo.

  * •

The offset is also negatively correlated with $B^{2}$, and positively
correlated with $L_{\rm c}$. Our calculations suggest that in the inter-arm
region, and in the Galactic halo (where the coherence length is larger), the
offsets between pulsar and the halo centroid could be as large as 50 pc.

We have also discussed the possible application of this mechanism to explain
some LHAASO sources with pulsar associations. For example, 1LHAASO
J0216+4237u, J0212+4254u, and J0206+4302u could have a common origin, and be
explained as mirage halos powered by the millisecond pulsar PSR J0218+4232.
The large offsets between a few middle-aged pulsars and LHAASO sources may
also be explained by our model naturally.

In our simulations, we do not consider the feedback of relativistic particles
on the magnetic field, such as the streaming instability. This simplification
is reasonable for energetic electrons ($>100$ TeV), because the energy density
of such high-energy electrons is usually very small compared with that of the
interstellar medium, and hence the growth of such instabilities should be very
slow.

We also note that the synthetic turbulent magnetic fields used in our
simulations are generated following Ref. [24], assuming an isotropic
turbulence with Kolmogorov-type power spectrum. The nature of the turbulence
in the interstellar medium may display some differences. In particular,
studies have shown that Alfvénic turbulence is highly anisotropic at the gyro-
scale of cosmic rays and thereby inefficient at scattering cosmic rays, while
fast modes are isotropic and should play a dominant role in scattering cosmic
rays [35], as supported by the recent analysis of radio polarisation
observations [36]. Therefore, some differences must exist between the cosmic-
ray scattering rates for particles propagating in our synthetic turbulence and
that for particles propagating in magnetohydrodynamic turbulence. Although our
synthetic turbulence may not be accurate on the cosmic-ray gyration scales, it
still qualitatively captures the isotropic and filamentary nature of the
interstellar turbulence on large scales, which is due to modes with large
wavelengths. These modes are actually those which drive the filamentary
propagation of cosmic-rays along magnetic field lines and thereby the
formation of asymmetric gamma-ray emissions [21], mirage sources, and large
offsets. Therefore, our calculations provide a reasonably good approximation
of realistic interstellar turbulence, for the study of the appearance of
asymmetric extended gamma-ray emissions around cosmic-ray sources.
Furthermore, this method can take into account a wide range of length scales
in the turbulence from the gyroradius of relevant electrons up to the
turbulence injection scale with reasonable computing times. More realistic
simulations, such as particle-in-cell simulations, or particle propagation in
3D cubes from magnetohydrodynamic simulations can, however, only cover a
limited range of spatial scales, limiting the ranges of electron energies that
can be studied with them.

###### Acknowledgements.

We thank Chunkai Yu, Hao Zhou and Giovanni Morlino for discussion. This work
is supported by NSFC under grants No. 12393852, 12350610239, 12393853,
12173018, 12121003, and U2031105.

## References

  * Sironi and Spitkovsky [2014] L. Sironi and A. Spitkovsky, ApJL 783, L21 (2014), eprint 1401.5471.
  * Cerutti et al. [2020] B. Cerutti, A. A. Philippov, and G. Dubus, Astronomy and Astrophysics 642, A204 (2020), eprint 2008.11462.
  * Lhaaso Collaboration et al. [2021] Lhaaso Collaboration, Z. Cao, F. Aharonian, Q. An, Axikegu, L. X. Bai, Y. X. Bai, Y. W. Bao, D. Bastieri, X. J. Bi, et al., Science 373, 425 (2021), eprint 2111.06545.
  * Sironi and Spitkovsky [2011] L. Sironi and A. Spitkovsky, Astrophys. J. 741, 39 (2011), eprint 1107.0977.
  * Cerutti and Philippov [2017] B. Cerutti and A. A. Philippov, Astronomy and Astrophysics 607, A134 (2017), eprint 1710.07320.
  * Giacinti and Kirk [2018] G. Giacinti and J. G. Kirk, Astrophys. J. 863, 18 (2018), eprint 1804.05056.
  * Cerutti and Giacinti [2020] B. Cerutti and G. Giacinti, Astronomy and Astrophysics 642, A123 (2020), eprint 2008.07253.
  * Abeysekara et al. [2017] A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, J. D. Álvarez, R. Arceo, J. C. Arteaga-Velázquez, D. Avila Rojas, H. A. Ayala Solares, A. S. Barber, et al., Science 358, 911 (2017), eprint 1711.06223.
  * Aharonian et al. [2021a] F. Aharonian, Q. An, L. X. Axikegu, Bai, Y. X. Bai, Y. W. Bao, D. Bastieri, X. J. Bi, Y. J. Bi, H. Cai, J. T. Cai, et al., Phys. Rev. Lett. 126, 241103 (2021a), eprint 2106.09396.
  * Evoli et al. [2018] C. Evoli, T. Linden, and G. Morlino, Phys. Rev. D 98, 063017 (2018), eprint 1807.09263.
  * Fang et al. [2019] K. Fang, X.-J. Bi, and P.-F. Yin, Mon. Not. R. Astron. Soc. 488, 4074 (2019), eprint 1903.06421.
  * Mukhopadhyay and Linden [2022] P. Mukhopadhyay and T. Linden, Phys. Rev. D 105, 123008 (2022), eprint 2111.01143.
  * López-Coto and Giacinti [2018] R. López-Coto and G. Giacinti, Mon. Not. R. Astron. Soc. 479, 4526 (2018), eprint 1712.04373.
  * Liu et al. [2019] R.-Y. Liu, H. Yan, and H. Zhang, Phys. Rev. Lett. 123, 221103 (2019), eprint 1904.11536.
  * Schroer et al. [2023] B. Schroer, C. Evoli, and P. Blasi, Phys. Rev. D 107, 123020 (2023), eprint 2305.08019.
  * Linden et al. [2017] T. Linden, K. Auchettl, J. Bramante, I. Cholis, K. Fang, D. Hooper, T. Karwal, and S. W. Li, Phys. Rev. D 96, 103016 (2017), eprint 1703.09704.
  * Sudoh et al. [2019] T. Sudoh, T. Linden, and J. F. Beacom, Phys. Rev. D 100, 043016 (2019), eprint 1902.08203.
  * Cao et al. [2024] Z. Cao, F. Aharonian, Q. An, Axikegu, Y. X. Bai, Y. W. Bao, D. Bastieri, X. J. Bi, Y. J. Bi, J. T. Cai, et al., ApJS 271, 25 (2024), eprint 2305.17030.
  * Albert et al. [2020] A. Albert, R. Alfaro, C. Alvarez, J. R. A. Camacho, J. C. Arteaga-Velázquez, K. P. Arunbabu, D. Avila Rojas, H. A. Ayala Solares, V. Baghmanyan, E. Belmont-Moreno, et al., Astrophys. J. 905, 76 (2020), eprint 2007.08582.
  * Zhang et al. [2021] Y. Zhang, R.-Y. Liu, S. Z. Chen, and X.-Y. Wang, Astrophys. J. 922, 130 (2021), eprint 2010.15731.
  * Giacinti et al. [2012a] G. Giacinti, M. Kachelrieß, and D. V. Semikoz, Phys. Rev. Lett. 108, 261101 (2012a), eprint 1204.1271.
  * Beck et al. [2016] M. C. Beck, A. M. Beck, R. Beck, K. Dolag, A. W. Strong, and P. Nielaba, JCAP 2016, 056 (2016), eprint 1409.5120.
  * Bao et al. [2024] Y. Bao, G. Giacinti, R. Liu, H. Zhang, and Y. Chen, Phys. Rev. Lett. , submitted (2024).
  * Giacalone and Jokipii [1999] J. Giacalone and J. R. Jokipii, Astrophys. J. 520, 204 (1999).
  * DeMarco et al. [2007] D. DeMarco, P. Blasi, and T. Stanev, JCAP 2007, 027 (2007), eprint 0705.1972.
  * Giacinti et al. [2012b] G. Giacinti, M. Kachelrieß, D. V. Semikoz, and G. Sigl, JCAP 2012, 031 (2012b), eprint 1112.5599.
  * Fang et al. [2021] K. Fang, X.-J. Bi, S.-J. Lin, and Q. Yuan, Chinese Physics Letters 38, 039801 (2021), eprint 2007.15601.
  * Zenitani and Umeda [2018] S. Zenitani and T. Umeda, Physics of Plasmas 25, 112110 (2018), eprint 1809.04378.
  * Aharonian et al. [2021b] F. Aharonian, Q. An, Axikegu, L. X. Bai, Y. X. Bai, Y. W. Bao, D. Bastieri, X. J. Bi, Y. J. Bi, H. Cai, et al., Chinese Physics C 45, 025002 (2021b), eprint 2010.06205.
  * Jansson and Farrar [2012a] R. Jansson and G. R. Farrar, Astrophys. J. 757, 14 (2012a), eprint 1204.3662.
  * Jansson and Farrar [2012b] R. Jansson and G. R. Farrar, ApJL 761, L11 (2012b), eprint 1210.7820.
  * Giacinti et al. [2015] G. Giacinti, M. Kachelrieß, and D. V. Semikoz, Phys. Rev. D 91, 083009 (2015), eprint 1502.01608.
  * Note [1] Note1, the cosmic-ray grammage is given by $\langle X\rangle=1.4m_{\rm p}n_{\rm p}\langle t_{\rm res}\rangle c$, where $t_{\rm res}$ is the residence time, $n_{\rm p}$ is the hydrogen density in the ISM, $m_{p}$ is the mass of proton, and the prefactor 1.4 takes into account the contribution of nuclei (such as helium) in the ISM.
  * Aguilar et al. [2023] M. Aguilar, L. Ali Cavasonza, B. Alpat, G. Ambrosi, L. Arruda, N. Attig, C. Bagwell, F. Barao, L. Barrin, A. Bartoloni, et al., Phys. Rev. Lett. 130, 211002 (2023).
  * Yan and Lazarian [2002] H. Yan and A. Lazarian, Phys. Rev. Lett. 89, 281102 (2002), eprint astro-ph/0205285.
  * Malik et al. [2024] S. Malik, K. H. Yuen, and H. Yan, Astrophys. J. 965, 65 (2024), eprint 2307.13342.