# Study of maximum electron energy of sub-PeV pulsar wind nebulae by
multiwavelength modelling

Jagdish C. Joshi1,2 , Shuta J. Tanaka3,4 , Luis Salvador Miranda5 , Soebur
Razzaque2,6,7
1 Aryabhatta Research Institute of Observational Sciences (ARIES), Manora
Peak, Nainital 263001, India
2 Centre for Astro-Particle Physics (CAPP) and Department of Physics,
University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa
3 Department of Physical Sciences, Aoyama Gakuin University, 5-10-1 Fuchinobe,
Sagamihara, Kanagawa 252-5258, Japan
4 Graduate School of Engineering, Osaka University, 2-1 Yamadaoka, Suita,
Osaka 565-0871, Japan
5 Department of Physics, The Chinese University of Hong Kong, Shatin, Hong
Kong China
6 Department of Physics, The George Washington University, Washington, DC
20052, USA
7 National Institute for Theoretical and Computational Sciences (NITheCS),
South Africa
jagdish<EMAIL_ADDRESS>
###### Abstract

Recently, the Large High Altitude Air Shower Observatory (LHAASO) reported the
discovery of 12 ultrahigh-energy (UHE; $\mathrm{\varepsilon}\geq 100$ TeV)
gamma-ray sources located in the Galactic plane. A few of these UHE gamma-ray
emitting regions are in spatial coincidence with pulsar wind nebulae (PWNe).
We consider a sample of five sources; two of them are LHAASO sources (LHAASO
J1908+0621 and LHAASO J2226+6057) and the remaining three are GeV-TeV gamma-
ray emitters. In addition, X-rays, radio observations or upper limits are also
available for these objects. We study multiwavelength radiation from these
sources by considering a PWN origin, where the emission is powered by spin-
down luminosity of the associated pulsars. In this leptonic emission model,
the electron population is calculated at different times under the radiative
(synchrotron and inverse-Compton) and adiabatic cooling. We also include the
onset of the reverberation phase for the PWN, by assuming radially symmetric
expansion. However, in this work, we find that multiwavelength emission can be
interpreted before the onset of this phase. The maximum energy of the
electrons based on the spectral fit is found to be above 0.1 PeV and close to
1 PeV. For LHAASO J2226+6057, using its observations in radio to UHE gamma-
rays, we find that UHE gamma-rays can be interpreted using electrons with
maximum energy of 1 PeV. We estimate the upper limits on the minimum Lorentz
factor of the electrons and it also infers the minimum value of the pair-
multiplicity of charged pairs.

###### keywords:

radiation mechanisms: non-thermal, relativistic processes, gamma-rays: stars,
ISM: cosmic rays, supernova remnants, stars: pulsars

##  1 Introduction

The PWN structure is energetically supported by the spin-down luminosity of
the central pulsar and its composition is dominated by $e^{\pm}$ pair-plasma
coupled with the magnetic field, as well as nuclei (Pacini & Salvati, 1973;
Rees & Gunn, 1974; Reynolds & Chevalier, 1984a; Arons & Tavani, 1994; Bednarek
& Protheroe, 1997; Amato et al., 2003; Chevalier, 2004; Gaensler & Slane,
2006; Volpi et al., 2008; Kirk et al., 2009; Bühler & Blandford, 2014;
Kashiyama, 2017; Gelfand, 2017; Torres, 2017; Amato & Olmi, 2021; López-Coto
et al., 2022). A strong pulsar wind makes their nebula brighter in gamma-rays,
and in general spin-down luminosity greater or equal to $4\times 10^{36}$
erg/s is sufficient (Gotthelf, 2004). In our Galaxy, PWNe is the dominant
class of very high energy (VHE; $100~{}{\rm GeV}\leq\mathrm{\varepsilon}<100$
TeV) gamma-ray sources, detected in the Galactic plane survey by the High
Energy Stereoscopic System (H.E.S.S.) telescope (H. E. S. S. Collaboration et
al., 2018b, a), also known as multiwavelength emitters (Reynolds et al.,
2017). Mattana et al. (2009) found that the production of TeV emission in PWNe
is not correlated with the spin-down luminosity and the characteristic age of
the pulsar. Further, for the VHE emission, the target photon field can be a
combination of synchrotron photons, cosmic microwave background (CMB) photons,
dust infrared (IR) and stellar photons (Tanaka & Takahara, 2010; Torres et
al., 2014; Zhu et al., 2018; Mares et al., 2021). For older PWNe, the magnetic
field is weaker, which makes them fainter in X-rays but due to the inverse-
Compton (IC) scattering in the CMB, IR radiation they remain brighter in
gamma-rays (de Jager et al., 2009). In older PWNe, the gamma-ray emission is
mostly due to the up-scattering of CMB photons or IR photons by relativistic
electrons (Torres et al., 2014); and for younger ($t_{\rm age}<300$ yr) PWNe,
this emission is dominated by the upscattering of the synchrotron photons
(Tanaka & Takahara, 2011). The detailed dynamical and radiative models of PWNe
are useful to understand the physical parameters of the progenitor supernova
(SN), energetics of the pulsar and its wind, properties of the surrounding
environment, etc. (Chevalier, 2005; Gelfand et al., 2009; Martín et al., 2016;
Bandiera et al., 2020). The time-dependent leptonic spectral evolution model
for PWNe can provide us with details about the electron population and
magnetic field as functions of their age. The cooling of these electrons in
the magnetic and radiation fields leads to a multiwavelength spectrum from
radio to gamma-rays. For example, in the Crab Nebula, synchrotron self-Compton
(SSC) mechanism and IC scattering off the CMB photons have been used to
explain the multiwavelength radiation (Tanaka & Takahara, 2010). Their
magnetic field evolution in time is also consistent with the rate of flux
decrease in radio wavelengths.

In general, the multiwavelength emission from PWNe can be modelled using
leptonic models (Zhang et al., 2008; Tanaka & Takahara, 2010, 2011; Martín et
al., 2012; Torres et al., 2013, 2014; Zhu et al., 2018) or lepto-hadronic
models (Atoyan & Aharonian, 1996; Bednarek & Bartosik, 2003; Zhang & Yang,
2009; Li et al., 2010). The IceCube collaboration used the gamma-ray flux
levels of 35 Galactic PWN sources and used stacking analysis to find the
neutrino signal (Aartsen et al., 2020). They found neutrino flux in the TeV-
PeV range is less than or equal to 4 $\%$ from these classes of objects.
Recently, 12 UHE, Galactic gamma-ray sources have been discovered by the
LHAASO and few of these gamma-ray sources are in spatial coincidence with some
of the PWNe (Cao et al., 2021; The LHAASO Collaboration, 2021). In particular,
the 1.1 PeV gamma-ray event was found to be associated with the Crab Nebula
(The LHAASO Collaboration, 2021) and the 1.4 PeV maximum energy photon was
correlated with the Cygnus OB2 region (Cao et al., 2021). The radio to UHE
gamma-ray emission from the Crab Nebula is consistent with the SSC+IC model
and constrains the size of the electron pevatron in between 0.025 to 0.1 pc
(The LHAASO Collaboration, 2021). They found that the luminosity in the PeV
electrons is approximately 0.5 $\%$ of the pulsar’s spin-down luminosity. This
opens up a new domain of UHE gamma-ray Astronomy and powerful pulsars in our
Galaxy play an important role in their origin (Albert et al., 2021). The UHE
gamma-ray detection is useful to test the theoretical models of electron
acceleration in the PeV range (Giacinti & Kirk, 2018; Breuhaus et al., 2021a).
Further, PeV gamma-ray detection in PWNe can constrain accelerator size, the
minimum acceleration rate, magnetic field, and the maximum Lorentz factor of
electrons, etc., (The LHAASO Collaboration, 2021). Also, the pulsar wind can
carry electrons of maximum energy and these are injected by the polar cap
potential regions in the pulsars (Bucciantini et al., 2011). The UHE gamma-ray
sources and maximum photon energy detected from them provide some hint that
particles are accelerating with maximum efficiency (de Oña Wilhelmi et al.,
2022).

Here, we study multiwavelength emission from PWNe powered by their associated
pulsars using the interactions between the non-thermal population of
relativistic electrons with magnetic and radiation fields. The paper is
organized as follows. In Section 2, we discuss the evolution of the PWN
radius, cooling timescales that affects the electron distributions, magnetic
field evolution, etc. In Section 3, we describe the PWN sources that have data
or upper limits in radio, X-rays, and gamma-rays and perform modelling using a
one-zone model. In Section 4, we discuss and conclude our results.

##  2 The Model

In this section, we have described our model, based on Tanaka & Takahara
(2010) and further included the impact of SN reverse shock on the PWN radius
inside a non-radiative supernova remnant (SNR) (Gelfand et al., 2009). The
compression of the PWN radius can enhance the magnetic field and as a result,
it affects the non-thermal radiation from relativistic electrons (Reynolds &
Chevalier, 1984b; van der Swaluw et al., 2001; Gelfand et al., 2009; Bandiera
et al., 2020). These effects are most important if the PWN age is greater than
10 kyr (López-Coto et al., 2022).

The spin-down luminosity of a pulsar at a given time can be estimated from the
observed quantities, i.e., period $P$ of the pulsar, it’s derivative $\dot{P}$
and moment of inertia $I$ of the neutron star (NS) (Gaensler & Slane, 2006).
Further, it enables a continuous supply of energy into particles and fields,
that changes according to the relation (Gaensler & Slane, 2006)

$L(t)=L_{0}\left(1+\frac{t}{\tau_{0}}\right)^{-(n+1)/(n-1)},$ (1)

where $L_{0}$ is the initial spin-down luminosity and braking index $n$ is set
equal to 3 for pulsars (PSRs) in our calculations. The total energy injected
by a pulsar to its nebula in its lifetime approximately lies in the range of
$\sim(1-5)\times 10^{-2}\,{E_{\rm SN}}$, where $E_{\rm SN}=10^{51}$ erg is the
SN kinetic energy (Gaensler & Slane, 2006; Bucciantini et al., 2011). The
pulsar age $t_{\rm age}$, characteristic timescales, $\tau_{c}=P/2\dot{P}$ and
$\tau_{0}$, are related by the relation (Tanaka & Takahara, 2011)

$\tau_{c}=\frac{n-1}{2}(\tau_{0}+t_{\rm age}).$ (2)

As we have taken $n=3$, in this case, the pulsar loses energy via its spin-
down by the magnetic dipole radiation (Gaensler & Slane, 2006). This value
also implies that $\tau_{c}>t_{\rm age}$ to get a positive value of
$\tau_{0}$. Note that the uncertainty in the age of the pulsar $t_{\rm age}$
affects the model parameters (Tanaka & Takahara, 2013). Also, $\tau_{0}$ and
$t_{\rm age}$ play important roles in shaping the spectral energy distribution
(SED) of a PWN (Torres et al., 2014). As $\tau_{c}$ is known from pulsar
observations and for a known or assumed value of $t_{\rm age}$, we can
estimate $\tau_{0}$. Further, these calculation depends on $n$, however, in
our calculations we have taken $n=3$.

###  2.1 Evolution of PWN Radius

The radius of the PWN can be estimated analytically and is defined as (Blondin
et al., 2001)

$\begin{split}R_{\rm PWN}(t)\simeq 0.5~{}{\rm pc}~{}\left(\frac{E_{\rm
SN}}{10^{51}~{}\rm erg}\right)^{\frac{3}{10}}\left(\frac{M_{\rm
ej}}{8~{}M_{\odot}}\right)^{-\frac{1}{2}}\\\
\times\left(\frac{L_{0}}{10^{38}~{}{\rm
erg/s}}\right)^{\frac{1}{5}}\left(\frac{t}{{\rm
500~{}yr}}\right)^{\frac{6}{5}}.\end{split}$ (3)

Figure 1: The SN reverse shock radius (shown in grey color) and the evolution
of the PWN radius. Once the SN reverse shock radius collides with the PWN
radius the compression and later the re-expansion takes place. For better
visibility of the curves, we have divided the radius values for HESS
J1640-465, HESS J1813-178, LHAASO J2226+6057, and HESS J1303-631 by factors
50, 30,10, and 3 respectively. The circled point in each curve represents the
radius at the current age of the PWN for which multiwavelength emission is
calculated.

This analytical solution is only valid up to pulsar age $t\ll\tau_{0}$.
$M_{\rm ej}$ is the amount of mass ejected during the SN explosion. For the
constant density $n_{\rm ISM}$ of the interstellar medium (ISM), the ejecta
phase remains dominant up to a time $t_{0}=196~{}{\rm yr}~{}(M_{\rm
ej}/M_{\odot})^{5/6}E_{\rm SN,51}^{-1/2}n_{\rm ISM}^{-1/3}$ (Truelove & McKee,
1999). Afterward, during the Sedov-Taylor (ST) phase, the SNR shock radius is
given by $R_{\rm SNR}=1.17\times\left(E_{\rm SN}t^{2}/\rho\right)^{1/5}$,
where the density $\rho=\mu_{\rm H}n_{\rm ISM}$ and the mean molecular weight
of ISM is $\mu_{H}=1.4\times 1.67\times 10^{-24}~{}{\rm g/cm^{3}}$ (Shu,
1992). Using these initial conditions we study the evolution of the PWN radius
inside a non-radiative SNR (Gelfand et al., 2009). Based on this, the first
compression starts when the PWN radius is equal to the SN reverse shock
radius. After this stage, the pressure balance between the PWN volume and the
ejecta material decides the PWN radius evolution. The details of the radius
evolution are modelled based on the formalism discussed by Gelfand et al.
(2009) and the model is only valid up to the time when the ST phase ends. In
Figure 1, the dynamics of the PWN radius are shown w.r.t. the pulsar or PWN
age. Initially, the PWN radius expands freely and after its collision with the
SN reverse shock radius, the compression phase starts under the condition of
sub-dominant pressure inside the PWN volume compared to the SN ejecta. In an
opposite scenario, for later times, the PWN volume re-expands, and so on. The
compression and re-expansion of the PWN volume are termed the reverberation
phase. For our input parameters, in Table 1, we have listed the collision time
$t_{\rm coll}$, after which compression starts. The onset of this phase
dominates during the end of the ejecta phase or onward the ST phase. However,
the details of the evolution also depend on the pulsar energetics and
parameters. Based on Figure 1, we also infer that for lower ejecta mass the
compression starts early in time. We have kept the same value of SN energy and
ISM density for all the objects, however, the onset of the compression phase
is also sensitive to these parameters. The circled point in each curve
represents the radius of the PWN at its current age $t_{\rm age}$, and in this
work, we have investigated the multiwavelength radiation from all sources in
the pre-contraction phase.

Figure 2: The magnetic field inside the PWN and its evolution with time.
Initially, the B-field decreases with the expansion of the pulsar wind and
once the reverberation phase starts then the B-field undergoes expansion and
re-compression phases, respectively. The circled point in each curve
represents the magnetic field at the current age of the PWN for which
multiwavelength emission is calculated. The magnetic field is estimated using
equation 12.

###  2.2 Particle Distribution Under Radiative and Adiabatic Cooling

In our one-zone model, the time-dependent evolution of electron population
$N(\gamma,t)$ at a time $t$, in the spherically expanding PWN can be
determined by their energy losses in the magnetic field, photon fields, and in
addition via adiabatic losses. The particle distribution $N(\gamma,t)$ under
energy losses $\dot{\gamma}(\gamma,t)$ and source term $Q(\gamma,t)$ can be
found using the formalism discussed in Tanaka & Takahara (2010). The
continuity equation for $N(\gamma,t)$ is defined as

$\frac{\partial N(\gamma,t)}{\partial
t}+\frac{\partial}{\partial\gamma}[\dot{\gamma}(\gamma,t)N(\gamma,t)]=Q(\gamma,t),$
(4)

where $\gamma$ is the electron Lorentz factor. The second term on the left-
hand side of equation 4, describes the cooling of relativistic electrons. We
have considered the cooling rate due to adiabatic expansion $\dot{\gamma}_{\rm
ad}(\gamma,t)$, synchrotron $\dot{\gamma}_{\rm Sy}(\gamma,t)$ and IC
$\dot{\gamma}_{\rm ic}(\gamma)$ scattering, respectively (Blumenthal & Gould,
1970; Tanaka & Takahara, 2010). The adiabatic losses affect the low energy
part of the electron distribution while the synchrotron and IC losses are
important for larger values of the electron Lorentz factor (Tanaka & Takahara,
2010; Vorster et al., 2013).

We use the adiabatic cooling time as $t_{\rm ad}=|\gamma/\dot{\gamma}_{\rm
ad}|=R_{\rm PWN}/V_{\rm PWN}$, and for our modelling this choice is okay as we
are mostly in the expansion stage of the PWN radius (Martín et al., 2012). The
values of PWN radius $R_{\rm PWN}$ and velocity $V_{\rm PWN}$ at the current
age are listed in Table 1 and same used for the estimation of $t_{\rm ad}$, as
shown in Figures 10 to 12. The electrons are cooled by synchrotron radiation
in the magnetic field and the cooling rate for a single electron in the
magnetic field is $\dot{\gamma}_{\rm
Sy}(\gamma,t)=-4c\sigma_{T}\gamma^{2}U_{B}(t)/3m_{e}c^{2}$, where $c$ is the
speed of light, $\sigma_{T}$ is the Thomson cross section, $U_{B}$ is the
magnetic energy density and $m_{e}$ is the electron mass. The synchrotron
cooling time-scale is (Blumenthal & Gould, 1970)

$t_{\rm Sy}=\frac{6\pi m_{e}c^{2}}{c\sigma_{T}B^{2}\gamma}=2.5\times
10^{3}{\rm yr}\left(\frac{B(t)}{100~{}\mu{\rm
G}}\right)^{-2}\left(\frac{\gamma}{10^{6}}\right)^{-1}.$ (5)

Further, the cooling rate of electrons in the target photon field including
the Klein-Nishina (KN) scattering regime is $\dot{\gamma}_{\rm
IC}(\gamma,t)\approx-4c\sigma_{T}\gamma^{2}U_{\rm
ph}(t)/3m_{e}c^{2}(1+4\gamma\epsilon_{0}/m_{e}c^{2})^{3/2}$, and corresponding
time-scale is (Moderski et al., 2005)

$\begin{split}t_{\rm IC}=\frac{3}{4}\frac{m_{e}c^{2}}{c\sigma_{T}U_{\rm
ph}(t)\gamma}\left[1+\frac{4\gamma\epsilon_{0}}{m_{e}c^{2}}\right]^{3/2}\\\
=6.1\times 10^{5}{\rm yr}\left(\frac{U_{\rm ph}}{1~{}{\rm
eV\,cm^{-3}}}\right)^{-1}\left(\frac{\gamma}{10^{6}}\right)^{-1}\left[1+\frac{4\gamma\epsilon_{0}}{m_{e}c^{2}}\right]^{3/2}.\end{split}$
(6)

The target photons for the IC mechanism are synchrotron photons, CMB, IR
photons from dust grains and stellar optical photons. In the IC scattering,
the Thomson regime is valid if the target photon of energy $\epsilon_{0}$
interacts with an electron with Lorentz factor $\gamma$ satisfying the
condition $4\gamma\epsilon_{0}/m_{e}c^{2}\ll 1$ and later the KN effects are
important when $4\gamma\epsilon_{0}/m_{e}c^{2}\gg 1$ (Jones, 1968; Blumenthal
& Gould, 1970; Fang et al., 2021). While interacting with the graybody photons
with temperature $T$, KN effects become important when the energy of the
relativistic electron is larger than $0.27m_{c}^{2}c^{4}/k_{B}T$ (Schlickeiser
& Ruppel, 2010). Thus KN process reduces the energy loss for the electron and
a harder spectrum for non-thermal electrons is expected. These electrons
interact with the CMB radiation field with energy density $U_{\rm
CMB}=0.26~{}{\rm eV/cm^{3}}$ but the IR and stellar photon density depend on
the source location in the Galaxy (Porter et al., 2006). This should create a
harder electron spectrum above the electron Lorentz factor $\gamma_{\rm
KN}=0.27m_{e}c^{2}/k_{B}T$ (Schlickeiser & Ruppel, 2010), i.e., $6\times
10^{8},8\times 10^{7}$ and $3\times 10^{5}$, for the pair-plasma interactions
with the CMB, IR and stellar photons, respectively.

The source term $Q(\gamma,t)$ is due to the pulsar wind that depends on time.
In a phenomenological model, the broadband emission from PWN is generally
explained using a broken power-law distribution of electrons. The broken power
law has a harder index $p_{1}$ in the range from 1 to 1.8, below the break
Lorentz factor $\gamma_{b}$ and above it a softer index $p_{2}$ can take
values in the range from 2 to 3.1 (Torres et al., 2013; Tanaka & Asano, 2017;
Zhu et al., 2018). The time-dependent injected spectrum of electrons from the
pulsar into the PWN can be either a power-law or a broken power-law type.

$\displaystyle Q(\gamma,t)$
$\displaystyle=Q_{0,e}\begin{cases}\left(\frac{\gamma}{\gamma_{b}}\right)^{-p1},~{}~{}~{}~{}~{}~{}~{}\gamma_{\rm
min}<\gamma\leq\gamma_{b}\\\
\left(\frac{\gamma}{\gamma_{b}}\right)^{-p2},~{}~{}~{}~{}~{}~{}~{}\gamma_{b}<\gamma<\gamma_{\rm
max}.\end{cases}$ (7)

The normalization $Q_{0,e}$ can be estimated using the relation,

$\eta_{e}L(t)=\int_{\gamma_{\rm min}}^{\gamma_{\rm max}}Q(\gamma,t)\gamma
m_{e}c^{2}d\gamma,$ (8)

where $\eta_{e}=L_{e}(t)/L(t)$ is the fraction of total energy into electrons.
Further, we consider that $\eta_{B}=L_{B}(t)/L(t)$ is the fraction of total
pulsar energy, available for the amplification of the magnetic energy.
Approximately $10\%$ of the pulsar spin-down energy radiates via pulsed
emission (Vorster et al., 2013). These fractional parameters are used to
estimate the value of the magnetization parameter $\sigma$, which is the ratio
between magnetic energy flux and the particle energy flux,
$\sigma=L_{B}(t)/(L_{e}(t)+L_{A}(t))\approx\eta_{B}$ (Kennel & Coroniti,
1984). Here, $L_{e}(t)+L_{A}(t)$ represents the total pulsar energy
distribution in $e^{\pm}$ pairs and nuclei of mass number $A$. In this work,
we assume $L_{A}(t)=0$ and hence all the injected energy is distributed in
between the magnetic field and the $e^{\pm}$ pairs. The value of $Q_{0,e}$ can
be estimated by using the following expression

$Q_{0,e}=\frac{\eta_{e}L_{0}}{m_{e}c^{2}}\left(1+\frac{t}{\tau_{0}}\right)^{-\frac{n+1}{n-1}}\left[\frac{\gamma_{b}^{p}\gamma_{\rm
max}^{2-p}}{2-p}-\frac{\gamma_{b}^{p}\gamma_{\rm
min}^{2-p}}{2-p}\right]^{-1}.$ (9)

While most of the particles inside the PWN are accelerated at the termination
shock region but the polar cap potential is useful to scale the maximum energy
of particles (Bucciantini et al., 2011). The polar cap (PC) potential injected
maximum energy of the particles is (Goldreich & Julian, 1969)

$E_{\rm max,PC}\approx 6\times 10^{12}~{}{\rm
eV}~{}\left(\frac{B_{p}}{10^{12}~{}{\rm G}}\right)\left(\frac{R_{\rm
NS}}{10~{}{\rm km}}\right)^{3}\left(\frac{P}{1~{}{\rm s}}\right)^{-2},$ (10)

where $B_{p}\approx 3.2\times 10^{19}(P\dot{P})^{1/2}~{}$G is the magnetic
field on the surface of the NS, the radius of the NS is $R_{\rm NS}\approx 10$
km. The second limit can be derived by comparing the synchrotron cooling time-
scale of $e^{\pm}$ pair-plasma in the magnetic field $B$ of the termination
shock with the acceleration time-scale $t_{\rm acc}=\pi\gamma
m_{e}c^{2}/eBc\eta_{\rm acc}$, i.e., time taken for completing half of a full
gyration (Giacinti & Kirk, 2018). This is given by the expression

$E_{\rm max,cool}=\left(\frac{6m_{e}^{2}c^{4}e\eta_{\rm
acc}}{\sigma_{T}B}\right)^{1/2}\approx 3.4~{}{\rm
PeV}\left(\frac{B(t)}{100~{}\mu G}\right)^{-1/2}\left(\frac{\eta_{\rm
acc}}{1}\right)^{1/2},$ (11)

where $\eta_{\rm acc}\leq 1$ is the acceleration efficiency parameter (de
Jager et al., 1996). We note that this is an upper limit on the maximum energy
due to the cooling of particles, as effects due to Compton cooling are not
considered. However, this is sufficient for our purpose of finding the
possible reservoir for the injection of maximum energy of particles inside the
PWN. The value of $\eta_{\rm acc}\leq 1$ can be found by the maximum value of
the synchrotron photon energy and $\eta_{\rm acc}=(\epsilon_{\rm
max,Sy}/230~{}{\rm MeV}$) (Amato & Olmi, 2021), where $\epsilon_{\rm
max,Sy}=h\nu_{\rm max}$, is listed in Table 1. The break energy for a given
$\gamma$ is calculated by using the expression $h\nu=17.4~{}{\rm
keV}~{}(\gamma/10^{8})^{2}(B/100~{}{\mu\rm G})$ (Blumenthal & Gould, 1970).

The magnetic energy of the PWN is supported by the spin-down luminosity $L(t)$
by an amount $\eta_{B}$. Also, adiabatic losses affect the total amount of
magnetic energy available inside the PWN. We use the following equation to
estimate the magnetic field, as also discussed earlier by (Martín et al.,
2016),

$\frac{dW_{B}}{dt}=\eta_{B}L(t)-\frac{W_{B}(t)}{R_{\rm PWN}(t)}\frac{dR_{\rm
PWN}(t)}{dt},$ (12)

where $W_{B}=B^{2}R_{\rm PWN}^{3}/6$. The evolution of the magnetic field is
shown in Figure 2 and during the compression phase magnetic field gets
amplified.

The average energy of the particles $E_{w}\approx\Gamma_{w}m_{e}c^{2}$ inside
the PWN is approximately $\gamma_{b}m_{e}c^{2}(\gamma_{b}/\gamma_{\rm
min})^{-p_{1}+1}$. We have also estimated the pair multiplicity $\kappa$ that
is defined as $\kappa=\int_{\gamma_{\rm min}}^{\gamma_{\rm
max}}Q(\gamma,t)d\gamma/\dot{n}_{\rm GJ}$, where $\dot{n}_{\rm
GJ}=(cI\Omega\dot{\Omega}/e^{2})^{1/2}$ is the Goldreich–Julian number flux
(Goldreich & Julian, 1969). The estimated values are listed in Table 1 and
provides insights on the pulsar environment.

Table 1: Model parameters for the PWN based on the observations and spectral
fitting. The symbols that are based on the observations are pulsar period $P$,
period derivative $\dot{P}$, distance to the pulsar $d$, characteristic age
$\tau_{c}$, surface magnetic field of the NS $B_{p}$, current spin-down
luminosity $\dot{E}$ respectively and wavelength dependent size of the PWN, if
available from observations. We have assumed a value for PSR age $t_{\rm
age}$, a fixed value for the input parameters; braking index $n$, interstellar
medium gas density $n_{\rm ISM}$, SN kinetic energy $E_{\rm SN}$, IR energy
density $U_{\rm IR}$ and optical energy density $U_{\rm Opt}$. However, for
HESS J1640-465, a large value of energy density in optical photons is
motivated from recent observations. We fit the pulsar or PWN age $t_{\rm
age}$, ejecta mass $M_{\rm ej}$ of the progenitor star, spectral index
$p_{1}$, $p_{2}$, electron distribution with Lorentz factors $\gamma_{\rm
min},\gamma_{\rm b},\gamma_{\rm max}$ and fractions of the total energy in
electrons and magnetic field $\eta_{e}$ and $\eta_{B}$ based on the
multiwavelength observations. The derived parameters from the above
information are: pulsar spin-down timescale $\tau_{0}$, initial spin-down
luminosity of the pulsar $L_{0}$, PWN magnetic field $B$, radius $R_{\rm PWN}$
and velocity $V_{\rm PWN}$, etc., at $t_{\rm age}$. $t_{\rm coll}$ is the
collision time of the PWN radius with the SN reverse shock. Further, the value
of average Lorentz factor $\Gamma_{w}$, pair multiplicity $\kappa$ is
estimated. The last four parameters are electron acceleration efficiency
$\eta_{\rm acc}$ based on maximum synchrotron energy and corresponding maximum
electron Lorentz factor allowed by the synchrotron cooling $\gamma_{\rm
max,cool}$ and $\gamma_{\rm max,PC}$ is the allowed electron Lorentz factor
due to the pulsar injection.

Model parameters | LHAASO J1908+0621 | LHAASO J2226+6057 | HESS J1640-465 | HESS J1813-178 | HESS J1303-631 |
---|---|---|---|---|---|---
| (PSR J1907+0602) | (PSR J2229+6114) | (PSR J1640-4631) | (PSR J1813-1749 ) | (PSR J1301-6305) |
| From Past Observations |  |  |  |  |
$P$ (ms) | 106.6[1] | 51.6[3] | 206 [5] | 44.7[8] | 184 [10] |
$\dot{P}$ (s/s) |  $87.3\times 10^{-15}$[1] |  $78.3\times 10^{-15}$[3] |  $9.758\times 10^{-13}$[5] |  $1.265\times 10^{-13}$[8] | $2.65\times 10^{-13}[10]$ |
$d$ [kpc] |  $3.2\pm 0.6$ [2] | 3 [4] | 10 [6] | 6.2 [9] | 6.6 [11] |
$\tau_{c}$ (kyr) | 19.4 | 10.5 | 3.1 | 5.6 | 11 |
$B_{p}$ ($10^{12}~{}$G) | 3 | 2 | 14 | 2.4 | 7 |
$\dot{E}$[erg/s] |  $2.84\times 10^{36}$[1] |  $2.25\times 10^{37}$[3] |  $4.4\times 10^{36}$[5] |  $6.8\times 10^{37}$[8] |  $1.7\times 10^{36}$[10] |
Radio: Size | No counterpart [12] |  $200^{{}^{\prime\prime}}$ (3 pc)[16] |  $8^{{}^{\prime}}$ (23 pc)[19] |  $3^{{}^{\prime}}$ (5.4 pc)[22] | No counterpart [25] |
Xray:Size |  $12^{{}^{\prime\prime}}$ 0.2 pc [13] |  $200^{{}^{\prime\prime}}$ (3 pc)[17] |  $1.2^{{}^{\prime}}$ (3.5 pc) [20] |  $80^{{}^{\prime\prime}}$ (2.4 pc) [23] |  $2^{{}^{\prime}}$ (3.8 pc) [26] |
VHE(HESS):Size |  $0.34^{\circ}$ (19 pc) [14] |  |  $12^{{}^{\prime}}$ (35 pc) [21] |  $(2.2\pm 0.4)^{{}^{\prime}}$ (4 pc) [24] |  $(0.16\pm 0.02)^{\circ}$ (18.4 pc) [27] |
VHE(VERITAS):Size |  |  $0.27^{\circ}$ (14 pc)[18] |  |  |  |
UHE (LHAASO):Size |  $0.45^{\circ}$ (25 pc)[15] |  $0.49^{\circ}$ (25.6 pc) [15] |  |  |  |
| Assumed Parameters |  |  |  |  |
$t_{\rm age}$ [kyr] | 10 | 7 | 2.1 | 2.5 | 10 |
$n$ | 3 | 3 | 3 | 3 | 3 |
$n_{\rm ISM}$ [$\rm{cm^{-3}}$] | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 |
$E_{\rm SN}$ [$10^{51}\rm{erg}$] | 1 | 1 | 1 | 1 | 1 |
$U_{\rm IR}[{\rm eV/cm^{3}}]$ | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 |
$U_{\rm Opt}[{\rm eV/cm^{3}}]$ | 0.2 | 0.2 | 1300 [28] | 0.2 | 0.2 |
| Fitted Parameters |  |  |  |  |
$M_{\rm ej}$ [$M_{\odot}$] | 12 | 10 | 8 | 8 | 12 |
$p_{1},p_{2}$ | 1.3,2.32 | 1.9, 2.48 | 1.5, 2.1 | 1.5,2.15 | 1.3,2.1 |
$\gamma_{\rm min},\gamma_{b},\gamma_{\rm max}$ | ($5\times 10^{2},8\times 10^{5},1.4\times 10^{9}$) | ($2\times 10^{2},10^{4},1.8\times 10^{9}$) | ($10^{2},10^{5},1.2\times 10^{9}$) | ($4\times 10^{2},10^{4},1.2\times 10^{9}$) | ($10^{4},10^{5},6\times 10^{8}$) |
$\eta_{e},\eta_{B}$ | ($\sim 1$, $6\times 10^{-4}$) | ($\sim 1$, $2.5\times 10^{-3}$) | ($\sim 1$, $10^{-4}$) | ($\sim 1$, $5\times 10^{-4}$) | ($\sim 1$, $5\times 10^{-4}$) |
| Derived Parameters |  |  |  |  |
$\tau_{0}$ [kyr] | 9.4 | 3.5 | 1 | 3.1 | 1 |
$L_{0}$ [erg/s] | $1.2\times 10^{37}$ | $2.0\times 10^{38}$ | $4.2\times 10^{37}$ | $2.2\times 10^{38}$ | $2.1\times 10^{38}$ |
$B(t_{\rm age})[\mu$G] | 0.55 | 2.4 | 1.6 | 4 | 0.5 |
$R_{\rm PWN}(t_{\rm age})$[pc] | 8.5 | 11 | 2.1 | 4 | 13.8 |
$V_{\rm PWN}(t_{\rm age})$[km/s] | 1034 | 1750 | 1151 | 1733 | 1297 |
$t_{\rm coll}$[kyr] | 12.4 | 7.3 | 8.7 | 6.3 | 10.3 |
$\Gamma_{w}$ | $8.7\times 10^{4}$ | $1.6\times 10^{2}$ | $3.2\times 10^{3}$ | $2\times 10^{3}$ | $5\times 10^{4}$ |
$\kappa$ | $6.5\times 10^{4}$ | $5.4\times 10^{7}$ | $2.2\times 10^{6}$ | $1.5\times 10^{7}$ | $9\times 10^{4}$ |
$\eta_{\rm acc}$ | $7.4\times 10^{-5}$ | $6\times 10^{-4}$ | $1.7\times 10^{-4}$ | $4.3\times 10^{-4}$ | $1.4\times 10^{-5}$ |
$\varepsilon_{\rm max,Sy}$ [eV] | $1.7\times 10^{4}$ | $1.3\times 10^{5}$ | $4\times 10^{4}$ | $10^{5}$ | $3.1\times 10^{3}$ |
$\gamma_{\rm max,PC}$ | $3.2\times 10^{9}$ | $9\times 10^{9}$ | $4\times 10^{9}$ | $1.4\times 10^{10}$ | $2.5\times 10^{9}$ |
$\gamma_{\rm max,cool}$ | $8.1\times 10^{8}$ | $10^{9}$ | $7\times 10^{8}$ | $7\times 10^{8}$ | $3.5\times 10^{8}$ |

  * •

References: [1]Abdo et al. (2010a), [2]Abdo et al. (2010b), [3]Abdo et al.
(2010a), [4]Halpern et al. (2001a), [5]Gotthelf et al. (2014), [6]Lemiere et
al. (2009), [7] Archibald et al. (2016a), [8]Halpern et al. (2012), [9]Yao et
al. (2017) [10] Manchester et al. (2005), [11]Cordes & Lazio (2002). [12]
Duvidovich et al. (2020), [13] Abdo et al. (2010b), [14] Aharonian et al.
(2009), [15] Cao et al. (2021) [16] Halpern et al. (2001b), [17] Halpern et
al. (2001a), [18] Acciari et al. (2009), [19] Whiteoak & Green (1996), [20]
Gotthelf et al. (2014), [21] Abramowski et al. (2014a), [22] Brogan et al.
(2005), [23] Ubertini et al. (2005); Funk et al. (2007); Helfand et al.
(2007), [24] Aharonian et al. (2006), [25] Sushch et al. (2017), [26] H. E. S.
S. Collaboration et al. (2012), [27] Aharonian et al. (2005b), [28] Mares et
al. (2021).

##  3 multiwavelength Emission Modelling under Adiabatic and Radiative Cooling

In this Section, we describe source properties, their multiwavelength
radiation, and their known interpretation and then we discuss the source
parameters based on our modelling. The criteria for the source selection is,
the distances should be known, and most plausibly the X-ray and radio
observations or upper limits are available. In the list of five sources we
have studied, two of them are detected by the LHAASO detector in the TeV-PeV
band. In the gamma-ray band, we have also included data from the Fermi-Large
Area Telescope (LAT), H.E.S.S., Very Energetic Radiation Imaging Telescope
Array System (VERITAS), High Altitude Water Cherenkov Observatory (HAWC) and
MILAGRO, if available.

Based on our calculation, we have estimated the population of non-thermal
electrons at different epochs. These electrons radiate from radio to UHE
gamma-rays and these emissions are calculated using the python package NAIMA
(Zabalza, 2015). The synchrotron intensity is calculated using the formalism
discussed in Aharonian et al. (2010); Baring et al. (1999). Further, the IC of
synchrotron and thermal photons is calculated using the formalism discussed in
Khangulyan et al. (2014). We evolve the electron population for each object up
to its current age and the electron population inside the PWN at four random
epochs is shown in the top panel of Figures 3 $-$ 7\. In the respective bottom
panels, we have shown the spectral fit to the observational data at the
current age of the PWN. The breaks in the SED are either due to the injection
spectrum of electrons $\gamma_{\rm min},\gamma_{b},\gamma_{\rm max}$ or the
cooling break, defined as $\gamma_{c}(t)\sim 2.45\times 10^{6}(t/1{\rm
kyr})^{-1}(B(t)/100~{}{\mu\rm G})^{-2}$ (Tanaka & Takahara, 2011).

###  3.1 LHAASO J1908+0621

LHAASO J1908+0621 is one of the Galactic UHE gamma-ray sources reported by the
LHAASO collaboration with its possible connection to SNR G40.5-0.5, PSR
J1907+0602 and PSR J1907+0631 (Cao et al., 2021). This source was also
associated with the HAWC detected UHE gamma-ray source eHWC J1907+063
(Abeysekara et al., 2020). The spin-down luminosity of PSR 1907+0602 was
sufficient to support the TeV emission of MGRO J1908+06 (Abdo et al., 2010b)
and associated HESS detection of J1908+063 (Aharonian et al., 2009). In the
same region, the extended emission was reported by VERITAS; VER J1907+062 and
its radio emission properties are investigated in detail and no radio
counterpart was observed (Duvidovich et al., 2020). Recently, Crestan et al.
(2021) derived the X-ray upper limits for this source using the XMM-Newton
observations. The emission from this spatial region has been studied in the
case of MGRO J1908+06 as a Galactic pevatron for which the gamma-ray spectrum
is harder above 100 TeV (Crestan et al., 2021; Li et al., 2021). Further, a
PWN origin of the UHE gamma-rays has been investigated by Breuhaus et al.
(2021b); Crestan et al. (2021); Li et al. (2021).

We have also used our comprehensive model by considering the role of PSR
1907+0602 in explaining the multiwavelength emission. We have chosen this PSR
compared to PSR 1907+0631 (another pulsar in the same region), due to its
higher ($\sim 5$ times) current spin-down luminosity. The distance $d$ to this
object, the current period $P$ and the period derivative $\dot{P}$ is $3.2\pm
0.6$ kpc, 106 ms and $87.3\times 10^{-15}$, respectively (Abdo et al., 2010b).
Using the values of the period and its derivative, the corresponding value of
$\tau_{c}$ is $\sim 19.4$ kyr. The source is located at
$(l,b)=[40.49^{\circ},-0.81^{\circ}]$ or at $(R,z)=[6.41,-0.05]$ kpc (by
taking the distance to the Galactic centre 8.3 kpc) (Abdo et al., 2013). If we
consider the age of the pulsar equal to the age of the SNR G40.5-0.5, that is
found to be in between 25 to 40 kyr (Downes et al., 1980), then from equation
2, a negative value of $\tau_{0}$ is obtained. However, this issue can be
resolved by taking a lower value of PWN age. To interpret the multiwavelength
observations we have taken a fiducial value $t_{\rm age}=10$ kyr. The emitting
regions, constrained by the observations are $12^{{}^{\prime\prime}}$ in
X-rays (Abdo et al., 2010b), the VHE size is $0.34^{\circ}$ (Aharonian et al.,
2009) and the UHE size of $0.45^{\circ}$ (Cao et al., 2021). As seen from the
radius evolution curve for LHAASO 1908+0621 in Figure 1, the radius is in the
pre-compression phase and its value is approximately $R_{\rm PWN}=8.5$ pc at
the current age as shown by the blue circle in Figure 1. That is greater than
the X-ray size of the nebula but lower by a factor of 2 to 3 compared to VHE
and UHE gamma-ray emitting regions. We will discuss in section 4, the possible
reasons for the difference in the size between the model and observation.

The time-dependent non-thermal electron distributions and estimated SED at
$t_{\rm age}=10$ kyr are shown in the top and bottom panels of Figure 3. The
model parameters are listed in Table 1. The electron spectrum inside the PWN
is similar to a standard PWN and the maximum energy of the electrons is $\sim
734$ TeV, that is lower compared to the 1-zone model but higher compared to
the 2-component model, as described in Crestan et al. (2021). Our estimated
value is different as this might be due to detailed modelling of the gamma-ray
data. In their model, they found that a single accelerator is unable to
interpret all the sets of observations. In that case, a 2-component model by
Crestan et al. (2021) would be useful, however, our interpretation of the data
is appropriate for our purpose of the estimation of the maximum energy of
electrons inside the PWN.

Figure 3: Top panel: The temporal evolution of electron distribution in the
PWN at different times for LHAASO J1908+0621. Bottom panel: The corresponding
SED at $t_{\rm age}=10$ kyr. The gamma-ray data points are taken from Cao et
al. (2021); Breuhaus et al. (2021b) and the X-ray upper limit is taken from
Crestan et al. (2021).

###  3.2 LHAASO J2226+6057

LHAASO J2226+6057 is another Galactic UHE gamma-ray source reported by the
LHAASO collaboration (Cao et al., 2021). It has a spatial association with SNR
G106.3+2.7 and its PSR J2229+6114, which supports wind nebula (also called
“Boomerang”) (Kothes et al., 2001). The MILAGRO collaboration also reported
detections in the Boomerang PWN associated with PSR J2229+6114 (Abdo et al.,
2009). PSR J2229+6114 was a bright gamma-ray pulsar in the first FERMI-LAT
catalog of gamma-ray pulsars (Abdo et al., 2010a). From the associated SNR,
multi-TeV gamma-ray emissions were detected by the VERITAS detector named as
VER J2227+608 (Acciari et al., 2009). In the spatial region with VER
J2227+608, Fermi-LAT detections observed GeV gamma-rays (Xin et al., 2019).
The SNR G106.3+2.7 was detected in radio (Pineault & Joncas, 2000) and also
X-ray data is available for this source (Fujita et al., 2021).

We consider that PSR J2229+6114 with current spin-down luminosity
$\dot{E}=2.3\times 10^{37}$ erg/s can power the PWN and UHE radiation from
LHAASO J2226+6057. The distance to PSR J2229+6114 is $\sim 3$ kpc based on
X-ray absorptions (Halpern et al., 2001a). Note that the distance is uncertain
and the atomic hydrogen and molecular velocity infers a distance $\sim 0.8$
kpc (Kothes et al., 2001). We have taken $d\sim 3$ kpc in our modelling of
this source. The current period $P=51.6$ ms and period derivative
$\dot{P}=78.3\times 10^{-15}$ s/s implies the value of $\tau_{c}=10.5$ kyr
(Abdo et al., 2010a). For PSR J2229+6114, we have assumed the age of PWN is
$t_{\rm age}=7$ kyr, similar to the value assumed in Yu et al. (2022). The
source location is $(l,b)=[106.28^{\circ},2.83^{\circ}]$ or at
$(R,z)=[9.7,0.15]$ kpc (Abdo et al., 2013) in our Galaxy.

The origin of UHE gamma-rays has been tested based on the leptonic emission in
PWN (Yu et al., 2022) and SNR G106.3+2.7 was investigated as a potential
pevatron candidate by Tibet AS$\gamma$ Collaboration et al. (2021). The GeV-
TeV gamma-ray data of VER J2227+608 with a hard gamma-ray spectral index
$1.90\pm 0.04$ and a cutoff in the proton spectrum above 400 TeV indicate a
plausible pevatron candidate (Xin et al., 2019). In their PWN model (Yu et
al., 2022), they have used a power law type electron spectrum with a spectral
index in the range of 2.3 to 2.5. In comparison, we have used a broken power
law type electron distribution with spectral index 1.9 and 2.48, before and
after the break energy.

The time-dependent electron distribution at four epochs and SED at $t_{\rm
age}=7$ kyr are shown in the top and bottom panels of Figure 4, respectively.
Our model constrains the average size of the PWN at the current age is
approximately 11 pc, which is again in between the X-ray emitting and VHE/UHE
gamma-ray emitting nebula. The emitting regions, constrained by the
observations are $200^{{}^{\prime\prime}}$ in radio and X-rays (Halpern et
al., 2001b, a), the VHE size is $0.27^{\circ}$ (Acciari et al., 2009) and the
UHE size of $0.49^{\circ}$ (Cao et al., 2021). The electron spectrum inside
the PWN is similar to a standard PWN and the maximum energy of the electrons
is $\sim 900$ TeV. For this source radio and X-ray observations are available
and this provides better constraints on the PWN magnetic field. We find
$B(t_{\rm age})=2.4\mu$G for this source. In Figure 9, we have fitted the SED
with maximum electron energy $E_{\rm max}=0.9,1$ and 3 PeV respectively. Based
on the UHE gamma-rays, we find that electrons of maximum energy 1 PeV are
available inside this PWN. For 3 PeV, we find excess flux of UHE gamma-rays.

Figure 4: Top panel: The temporal evolution of electron distribution in the
PWN at different times for LHAASO J2226+6057. Bottom panel: The corresponding
SED at $t_{\rm age}=7$ kyr. The data points are taken from, Radio (Pineault &
Joncas, 2000), X-ray (Fujita et al., 2021), gamma-ray (Xin et al., 2019;
Acciari et al., 2009; Abdo et al., 2009; Cao et al., 2021; Breuhaus et al.,
2021b).

###  3.3 HESS J1640-465

This source was discovered by the HESS telescope during their survey of the
inner Galaxy (Aharonian et al., 2006). Slane et al. (2010); Xin et al. (2018)
found that the gamma-ray observations for this source can be explained by
considering a PWN origin. Particle acceleration and radiation due to neutral
pion decay in the SNR G338.3-0.0 can also account for the TeV emission in HESS
J1640-465 (Abramowski et al., 2014b; Tang et al., 2015; Supan et al., 2016;
Mares et al., 2021). The TeV emission from HESS J1640-465 is spatially
correlated with a SNR G338.3-0.0 and PSR J1640-4631 (Supan et al., 2016). The
NuSTAR X-ray observations of the pulsar enable the estimation of the braking
index $n=3.15\pm 0.03$ for this source (Archibald et al., 2016b). However, for
simplicity, we have used $n=3$, which is very close to the observed value. The
characteristic age $\tau_{c}=3113$ yr corresponding to the period $P=206$ ms
and $\dot{P}=9.758\times 10^{-13}$ s/s (Gotthelf et al., 2014), that infers a
spin-down luminosity of $4.4\times 10^{36}$ erg/s. In our calculations if we
take the age of the pulsar or PWN the same as the SNR, which is found to be in
between 5-8 kyr (Slane et al., 2010) then the spin-down timescale found to be
negative based on equation 2, hence, we have chosen a lower value of $t_{\rm
age}=2.1$ kyr, that provides $\tau_{0}\sim 1$ kyr. The distance to this source
is in between 8.5-13 kpc based on the HI absorption (Lemiere et al., 2009) and
we have taken a reference value $d=10$ kpc in our modelling.

Recently, Mares et al. (2021) discussed the PWN origin of the multiwavelength
observations of this source and they found a rapidly rotating pulsar with
initial period $P_{0}\sim 10$ ms is needed. They also required an extremely
intense UV photon field, with energy density $1.3~{}{\rm keV/cm^{3}}$. HESS
J1640-465 is a composite SNR source, and an young pulsar PSR J1640-4631 is
also located in the shell-type SNR $\rm G338.3-0.0$ (Gotthelf et al., 2014).

As also mentioned by Mares et al. (2021) about its closeness with a nearby
young massive stellar cluster, which is $8^{{}^{\prime}}$ away from HESS
J1640-465 (Davies et al., 2012). This cluster can provide a higher stellar
photon density at the pulsar location $(l,b)=[338.28^{\circ},-0.04^{\circ}]$
or at $(R,z)=[3.78,-0.005]$ kpc (H. E. S. S. Collaboration et al., 2018b). We
also found that for this source stellar photon density 1.3 ${\rm keV/cm^{3}}$
is required. Further, Mares et al. (2021) has taken the age of the pulsar
equal to 3 kyr, which is a slightly larger value compared to our modelled
value $t_{\rm age}=2.1$ kyr. This creates a difference in the spin-down
timescale. The magnetization parameter proposed by Mares et al. (2021),
$\eta_{B}=10^{-1}$ is very large compared to our modeled value
$\eta_{B}=10^{-4}$ and this is also contrary to the convention that PWNe are
particle dominated. Our used value of gas density $0.1~{}\rm cm^{-3}$ is 10
times higher comparatively and the spin-down timescale is larger by a factor
of 250. The diversity in the model parameters for the interpretation of the
PWN emission is evident between these models and infers the degeneracy in
input model parameters. The electron distribution at four epochs and the SED
at $t_{\rm age}=2.1$ kyr are shown in the top and bottom panels of Figure 5,
respectively. At this age, the modelled value of the PWN radius is
approximately 2.1 pc which is a factor of 2 lower compared to the observed
Xray size of the nebula $3.5$ (Gotthelf et al., 2014) and approximately 10
times lower compared to HESS measured size (Abramowski et al., 2014a). The
maximum energy of electrons based on the spectral fit is $\sim$ 137 TeV.

Figure 5: Top panel: The temporal evolution of electron distribution in the
PWN at different times for HESS J1640-465. The concave shape of the electron
distribution at late times is due to the dominant IC cooling in the dense
stellar photon environment. Bottom panel: The corresponding SED at $t_{\rm
age}=2.1$ kyr. The radio flux upper limits and the remaining broadband data
are extracted from Abramowski et al. (2014b); Mares et al. (2021).

###  3.4 HESS J1813-178

HESS J1813-178 was discovered during the HESS survey of Galactic sources
(Aharonian et al., 2005a, 2006). Its multiwavelength observations are
available in radio, X-ray (Brogan et al., 2005), and gamma-rays (Araya, 2018).
This source has spatial overlap with the SNR G12.82-0.02 and at the centre of
the SNR, PSR J1813-1749 is located (Gotthelf et al., 2014). The current spin-
down luminosity of the pulsar is $\dot{E}=(6.8\pm 2.7)\times 10^{37}$ erg/s
(Gotthelf & Halpern, 2009; Camilo et al., 2021). Using the pulsar period
$P=44.7$ ms and its derivative $\dot{P}=1.265\times 10^{-13}$ s/s, we
estimated $\tau_{c}=5.6$ kyr (Halpern et al., 2012).

The radio observations of $\rm G12.82-0.02$ infers a compact size
$3^{{}^{\prime}}$ and also, independent of distance, the age of the SNR is
between 285-2,500 yrs (Brogan et al., 2005). The distance of the source is
quite uncertain and a value of 4.7 kpc was used to study the non-thermal
radiation (Fang & Zhang, 2010). In recent studies, a large dispersion measure
was found for the pulsar and, based on the electron distribution model of
Cordes & Lazio (2002), the source distance was found to be $12\pm 2$ kpc. We
caution that the dispersion in the distance will affect the model parameters.
We have taken its updated distance of 6.2 kpc based on the electron
distribution model of Yao et al. (2017). The source location is
$(l,b)=[12.81^{\circ},-0.03^{\circ}]$ or at $(R,z)=[2.81,-0.0032$] kpc in our
Galaxy (Aharonian et al., 2006). The electron distribution at four epochs and
the SED at $t_{\rm age}=2.5$ kyr are shown in the top and bottom panels of
Figure 6, respectively. The emitting regions for this source in the X-ray band
infer a size of value 2.4 pc (Ubertini et al., 2005; Funk et al., 2007;
Helfand et al., 2007), while the radio region is of size 5.4 pc (Brogan et
al., 2005). Also, in the case of this source, the HESS observations infer a
very compact size (4 pc) of the emission region (Aharonian et al., 2006). Our
modelled value of the PWN radius is approximately 4 pc. Further, for this
source radio and X-ray observations are available and this provides better
constraints on the PWN magnetic field. We find $B(t_{\rm age})=4\mu$G for this
source. The maximum energy of the electrons is $\sim 628$ TeV.

Figure 6: Top panel: The temporal evolution of electron distribution in the
PWN at different times for HESS J1813-178. Bottom panel: The corresponding SED
at $t_{\rm age}=2.5$ kyr. The MW data for HESS J1813-178 is taken from (Fang &
Zhang, 2010).

###  3.5 HESS J1303-631

The multiwavelength (radio, X-ray, and gamma-ray) observations of HESS
J1303-631 can be explained through the PWN emission (H. E. S. S. Collaboration
et al., 2012). They have used a stationary leptonic model with magnetic field
value 1.4 $\mu G$. In our revised version, the magnetic field is lower and in
our modelling we evolve the electron population under the adiabatic and
cooling losses plus reverberation and interpret the data. The distance to this
source is 6.6 kpc based on the Galactic electron distribution model by Cordes
& Lazio (2002) and the IR energy density is $1.3~{}{\rm eV/cm^{3}}$ at the
location $(l,b)=[304.21^{\circ},-0.33^{\circ}]$ or at $(R,z)=[7.25,-0.04]$ kpc
(Aharonian et al., 2006). We revisit the source modelling using time-dependent
PWN emission. The source is associated with PSR J1301-6305 having period
$P=184$ ms, period derivative $\dot{P}=2.65\times 10^{-13}$ s/s and spin-down
luminosity $\dot{E}=1.7\times 10^{36}$ erg/s (Manchester et al., 2005; H. E.
S. S. Collaboration et al., 2012). This provides the characteristic age of
this source to be approximately 11 kyr. The electron distribution at four
epochs and the SED are shown in the top and bottom panel of Figure 7 for
$t_{\rm age}=10$ kyr. There is no radio counterpart for this source (Sushch et
al., 2017) but the X-ray size is of value $\sim 4$ pc (H. E. S. S.
Collaboration et al., 2012) and a larger emitting region in the VHE gamma-rays
of value 18.4 pc (Aharonian et al., 2005b). Our model value for the PWN radius
is 13.8 pc and it’s closer to the VHE gamma-ray emitting nebula. For our
selected input parameters the full part of the VHE gamma-ray spectrum is not
explained however, it was explained by H. E. S. S. Collaboration et al.
(2012). This difference might be due to the contrast difference in our
modelling and IR field $>0.2~{}\rm eV/cm^{3}$ can be useful for reproducing
the full VHE spectrum. However, we focus on the estimation of the maximum
energy of the electrons based on the maximum energy of the VHE photons. The
maximum energy of the electrons from the spectral fit is $\sim 307$ TeV.

Figure 7: Top panel: The temporal evolution of electron distribution in the
PWN at different times for HESS J1303-631. Bottom panel: The corresponding SED
at $t_{\rm age}=10$ kyr. The upper limits for the radio, X-ray data, and GeV-
TeV data of HESS J1303-631 are taken from (H. E. S. S. Collaboration et al.,
2012).

##  4 Discussions and Conclusions

In recent years, the VHE gamma-ray spectrum of Galactic sources has been
extended to the UHE range by the HAWC and LHAASO detectors. The discovery of
1.1 PeV photon from the Crab PWN makes it one of the first Galactic sources of
PeV photons and also a Galactic leptonic pevatron source (The LHAASO
Collaboration, 2021). Liu & Wang (2021) have shown that the end part of the
UHE gamma-ray spectrum is harder for the CRAB Nebula and infers the
acceleration of cosmic ray (CR) protons up to PeV energies. Motivated by these
results we have explored the origin of the UHE gamma-ray spectrum in two of
the LHAASO detected sources. We have taken the cooling of the electrons in the
KN regime (Moderski et al., 2005) while calculating the electron distribution
and in the SED calculation these modifications are already present in the
NAIMA public code based on the formalism by Khangulyan et al. (2014). Further,
we compare these model parameters of two LHAASO detected PWN with a few other
H.E.S.S. detected objects studied by us. We find the UHE detected objects can
be interpreted using the spin-down luminosity of pulsars and their model
parameters are similar to modelled VHE detected sources. Hence, most of the
VHE gamma-ray sources in our Galaxy are powered by the pulsar spin-down
luminosity then they should also be detected by UHE gamma-ray detectors.
Recently, this is also shown by Albert et al. (2021), that the powerful
pulsars with $\dot{E}\geq 10^{36}$ erg/s, would have a UHE gamma-ray spectrum.

Figure 8: The multiwavelength fit based maximum energy $E_{\rm max}$ vs polar
cap potential  1 0 -.18 1$\phi$. We have also included these scaling values
from (Zhu et al., 2018) (BLJ et al.(2018)) and (Tanaka & Takahara, 2011) (TT
et al.(2011)) for the comparison.The young PWNe have age less than 2.5 kyr,
and for evolved PWNe the age is above 2.5 kyr but less than 5 kyr and for
mature PWNe it is above 5 kyr respectively as defined in (Zhu et al., 2018).

In Table 1, we have listed parameters: (i) From Past Observations (ii) Assumed
parameters (iii) Parameters based on spectral fit (iv) Derived from Parameters
in (i) and (ii) and (iii). Our current estimation of the PWN radius can’t
interpret the wavelength dependent size of the emitting nebula. The variation
of size in radio to gamma-ray wavelengths is very explicit from observations.
We think that age is also an important factor in deciding the radius of the
PWN. Wavelength dependent effects on PWN size requires the inclusion of the
particle particle transport scenarios (Tang & Chevalier, 2012). The evolution
of the PWN radius during the ejecta and ST phase is shown in Figure 1. The
injected electron spectrum used in our study follows a broken power law. The
spectral index $p_{1}$ and $p_{2}$ are consistent with the standard PWN
interpretation. Due to the very weak magnetic field at the current age of the
PWN, the cooling Lorentz factor is lower than the minimum Lorentz factor for
all the sources. For the breaks in the injected electron spectrum, we can find
their signatures in their SEDs. The magnetic field inside LHAASO J1908+0621 is
0.55 $\mu$G and for HESS J1303-631 is 0.5 $\mu$ G and the origin of these low
magnetic field values inside PWN is not very well known. In our case, these
two objects are older, compared to others and that leads to these small values
in the expansion phase. However, similar low values have been reported for the
PWN modelling, for example, in the modelling of UHE gamma-ray source, HAWC
J1826-128 by Burgess et al. (2022).

Figure 9: The SED of LHAASO J2226+6057 is shown for maximum electron energy
$E_{\rm max}=0.9,1$ and 3 PeV, respectively. The other parameters are the same
as listed in Table 1 for this source. The UHE gamma-rays infer that electrons
of PeV energy can be inside the PWN.

The value of $\eta_{e}$ is dominant compared to $\eta_{B}$ for all sources and
implies that the PWN plasma is dominated by the $e^{\pm}$ pair-plasma. Similar
conclusions about the PWN composition were reported in earlier studies (Torres
et al., 2013). Further, using fitted parameters, we have estimated the average
energy per particle in the pulsar wind $\Gamma_{w}$ and the number of pairs
produced per photon, i.e., pair multiplicity $\kappa$. The values of $\kappa$
are approximately in between $\sim 10^{4}-10^{7}$.

Our assumed values for external radiation field energy density for the IR and
stellar photons are within the standard values as known from the Galactic
radiation field models (Popescu et al., 2017; Breuhaus et al., 2021b; Zhu et
al., 2018), however, these values can be location dependent. Surprisingly, For
the HESS J1640-465 source, a very large photon density is supported by a
nearby source as discussed in Section 3. Further, the temperature of the IR
radiation field can affect the IC radiation and lower values are useful to
produce UHE radiation (Breuhaus et al., 2021a), we have taken $T_{\rm IR}=20$
K in this work and this value is consistent with the dust temperature (Bernard
et al., 2010; Zhu & Huang, 2014). With more UHE source detection in the future
by LHAASO and CTA, the IR radiation environments can be tested.

The $\gamma-\gamma$ absorption in the ISM can affect the VHE to UHE part of
the gamma-ray spectrum (Moskalenko et al., 2006). However, for the two LHAASO
objects these effects on the TeV-PeV gamma-ray spectrum are negligible (Cao et
al., 2021). For HESS J1813-178 and HESS J1303-631, we found that attenuation
is not important. In the case of HESS J1640-465, these effects are dominant
due to the high density of the target photons. To minimize it, we have reduced
the target photon temperature to 5000 K and an exponential cut-off in the
gamma-ray spectrum is used above $\sim$ 50 TeV based on the pair production
condition $\epsilon_{o}\epsilon_{\gamma}\geq(2m_{e}c^{2})^{2}$, where
$\epsilon_{o}$ and $\epsilon_{\gamma}$ are the energy of optical and gamma-ray
photons, respectively.

Using the maximum synchrotron photon energy we calculate acceleration
efficiency $\eta_{\rm acc}$ of the electrons and we found its values in
between $10^{-4}-10^{-5}$, also listed in Table 1. Using, equation 11, we
obtained values of $\gamma_{\rm max,cool}$. These values are approximately 10
times lower than $\gamma_{\rm max,PC}$. The required values of $\gamma_{\rm
max}$ from the spectral fit are in between $\gamma_{\rm max,cool}$ and
$\gamma_{\rm max,PC}$. Hence, the maximum energy of the electrons is supported
by the polar cap potential (de Oña Wilhelmi et al., 2022). In Figure 8, we
have shown the scaling of the maximum energy of electrons based on the
multiwavelength fit, vs the polar cap potential $\mbox{\hskip 5.95834pt{1 0
-.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color$\phi$}}}=E_{\rm
max,PC}/e$. The cooling timescales for particles inside PWN are shown in
Figure 10-12, for the current age of the pulsar and input model parameters. It
is clear that within the acceleration timescale, adiabatic and synchrotron
cooling mechanisms dominate and limit the maximum energy of electrons at the
termination shock. To interpret the multiwavelength radiation we need another
source of electrons above these energies and hence, the injection of maximum
energy particles inside PWN must be due to the polar cap potential regions.
Our estimated values of the maximum electron energy from the spectral fit is
in the range 0.1-1 PeV. In Figure 9, we have shown the SED modelling of LHAASO
J2226+6057 at 1 and 3 PeV. We have selected it as for this object radio,
X-ray, and UHE observations are available). This shows that UHE gamma-ray
spectrum is a key to probing the PeV electron presence inside the PWN and
makes LHAASO J2226+6057 one of the potential leptonic pevatron candidates
compared to other sources investigated in this work.

Figure 10: LHAASO J1908+0621 and LHAASO J2226+6057 cooling timescales for the
model parameters listed in Table 1.

Figure 11: HESS J1640-465 and HESS J1813-178 cooling timescales for the model
parameters listed in Table 1. Figure 12: HESS J1303-631 cooling timescales for
the model parameters listed in Table 1.

##  5 Acknowledgements

JCJ is thankful to J. Carlos, N. Jingade, P. Kushwaha, R. Liu, V. Chand, E.
Amato, N. Gupta, and X.-Y. Wang for helpful discussions and S. Crestan for
providing the multiwavelength data files for the source MGRO J1908+06. SJT
would like special thanks to T. Tezuka for providing the base of the numerical
code used in this study. SJT is supported by Aoyama Gakuin University-
Supported Program “Early Eagle Program”. SR was supported by a grant from the
University of Johannesburg Research Council. Further, we are very thankful to
A. Joshi, N. Fraija and M. Cardillo for reading our work and for insightful
comments.

##  6 DATA AVAILABILITY

The data used in this work are available in the article and the code developed
for this work can be shared on reasonable request to the corresponding author.

## References

  * Aartsen et al. (2020) Aartsen M. G. et al., 2020, ApJ, 898, 117
  * Abdo et al. (2010a) Abdo A. A. et al., 2010a, ApJS, 187, 460
  * Abdo et al. (2010b) Abdo A. A. et al., 2010b, ApJ, 711, 64
  * Abdo et al. (2013) Abdo A. A. et al., 2013, ApJS, 208, 17
  * Abdo et al. (2009) Abdo A. A. et al., 2009, ApJ, 700, L127
  * Abeysekara et al. (2020) Abeysekara A. U. et al., 2020, Phys. Rev. Lett., 124, 021102
  * Abramowski et al. (2014a) Abramowski A. et al., 2014a, ApJ, 794, L1
  * Abramowski et al. (2014b) Abramowski A. et al., 2014b, MNRAS, 439, 2828
  * Acciari et al. (2009) Acciari V. A. et al., 2009, ApJ, 703, L6
  * Aharonian et al. (2009) Aharonian F. et al., 2009, A&A, 499, 723
  * Aharonian et al. (2005a) Aharonian F. et al., 2005a, Science, 307, 1938
  * Aharonian et al. (2005b) Aharonian F. et al., 2005b, A&A, 439, 1013
  * Aharonian et al. (2006) Aharonian F. et al., 2006, ApJ, 636, 777
  * Aharonian et al. (2010) Aharonian F. A., Kelner S. R., Prosekin A. Y., 2010, Phys. Rev. D, 82, 043002
  * Albert et al. (2021) Albert A. et al., 2021, ApJ, 911, L27
  * Amato et al. (2003) Amato E., Guetta D., Blasi P., 2003, A&A, 402, 827
  * Amato & Olmi (2021) Amato E., Olmi B., 2021, Universe, 7, 448
  * Araya (2018) Araya M., 2018, ApJ, 859, 69
  * Archibald et al. (2016a) Archibald R. F. et al., 2016a, ApJ, 819, L16
  * Archibald et al. (2016b) Archibald R. F. et al., 2016b, ApJ, 819, L16
  * Arons & Tavani (1994) Arons J., Tavani M., 1994, ApJS, 90, 797
  * Atoyan & Aharonian (1996) Atoyan A. M., Aharonian F. A., 1996, MNRAS, 278, 525
  * Bandiera et al. (2020) Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2020, MNRAS, 499, 2051
  * Baring et al. (1999) Baring M. G., Ellison D. C., Reynolds S. P., Grenier I. A., Goret P., 1999, ApJ, 513, 311
  * Bednarek & Bartosik (2003) Bednarek W., Bartosik M., 2003, A&A, 405, 689
  * Bednarek & Protheroe (1997) Bednarek W., Protheroe R. J., 1997, Phys. Rev. Lett., 79, 2616
  * Bernard et al. (2010) Bernard J. P. et al., 2010, A&A, 518, L88
  * Blondin et al. (2001) Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 563, 806
  * Blumenthal & Gould (1970) Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics, 42, 237
  * Breuhaus et al. (2021a) Breuhaus M., Hahn J., Romoli C., Reville B., Giacinti G., Tuffs R., Hinton J. A., 2021a, ApJ, 908, L49
  * Breuhaus et al. (2021b) Breuhaus M., Reville B., Hinton J. A., 2021b, arXiv e-prints, arXiv:2109.05296
  * Brogan et al. (2005) Brogan C. L., Gaensler B. M., Gelfand J. D., Lazendic J. S., Lazio T. J. W., Kassim N. E., McClure-Griffiths N. M., 2005, ApJ, 629, L105
  * Bucciantini et al. (2011) Bucciantini N., Arons J., Amato E., 2011, MNRAS, 410, 381
  * Bühler & Blandford (2014) Bühler R., Blandford R., 2014, Reports on Progress in Physics, 77, 066901
  * Burgess et al. (2022) Burgess D. A. et al., 2022, ApJ, 930, 148
  * Camilo et al. (2021) Camilo F., Ransom S. M., Halpern J. P., Roshi D. A., 2021, ApJ, 917, 67
  * Cao et al. (2021) Cao Z. et al., 2021, Nature, 594, 33
  * Chevalier (2004) Chevalier R. A., 2004, Advances in Space Research, 33, 456
  * Chevalier (2005) Chevalier R. A., 2005, ApJ, 619, 839
  * Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, arXiv e-prints, astro
  * Crestan et al. (2021) Crestan S., Giuliani A., Mereghetti S., Sidoli L., Pintore F., La Palombara N., 2021, MNRAS, 505, 2309
  * Davies et al. (2012) Davies B., de La Fuente D., Najarro F., Hinton J. A., Trombley C., Figer D. F., Puga E., 2012, MNRAS, 419, 1860
  * de Jager et al. (2009) de Jager O. C. et al., 2009, arXiv e-prints, arXiv:0906.2644
  * de Jager et al. (1996) de Jager O. C., Harding A. K., Michelson P. F., Nel H. I., Nolan P. L., Sreekumar P., Thompson D. J., 1996, ApJ, 457, 253
  * de Oña Wilhelmi et al. (2022) de Oña Wilhelmi E., López-Coto R., Amato E., Aharonian F., 2022, arXiv e-prints, arXiv:2204.09440
  * Downes et al. (1980) Downes A. J. B., Pauls T., Salter C. J., 1980, A&A, 92, 47
  * Duvidovich et al. (2020) Duvidovich L., Petriella A., Giacani E., 2020, MNRAS, 491, 5732
  * Fang & Zhang (2010) Fang J., Zhang L., 2010, ApJ, 718, 467
  * Fang et al. (2021) Fang K., Bi X.-J., Lin S.-J., Yuan Q., 2021, Chinese Physics Letters, 38, 039801
  * Fujita et al. (2021) Fujita Y., Bamba A., Nobukawa K. K., Matsumoto H., 2021, ApJ, 912, 133
  * Funk et al. (2007) Funk S. et al., 2007, A&A, 470, 249
  * Gaensler & Slane (2006) Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17
  * Gelfand (2017) Gelfand J. D., 2017, in Astrophysics and Space Science Library, Vol. 446, Modelling Pulsar Wind Nebulae, Torres D. F., ed., p. 161
  * Gelfand et al. (2009) Gelfand J. D., Slane P. O., Zhang W., 2009, ApJ, 703, 2051
  * Giacinti & Kirk (2018) Giacinti G., Kirk J. G., 2018, ApJ, 863, 18
  * Goldreich & Julian (1969) Goldreich P., Julian W. H., 1969, ApJ, 157, 869
  * Gotthelf (2004) Gotthelf E. V., 2004, in Young Neutron Stars and Their Environments, Camilo F., Gaensler B. M., eds., Vol. 218, p. 225
  * Gotthelf & Halpern (2009) Gotthelf E. V., Halpern J. P., 2009, ApJ, 700, L158
  * Gotthelf et al. (2014) Gotthelf E. V. et al., 2014, ApJ, 788, 155
  * H. E. S. S. Collaboration et al. (2018a) H. E. S. S. Collaboration et al., 2018a, A&A, 612, A2
  * H. E. S. S. Collaboration et al. (2018b) H. E. S. S. Collaboration et al., 2018b, A&A, 612, A1
  * H. E. S. S. Collaboration et al. (2012) H. E. S. S. Collaboration et al., 2012, A&A, 548, A46
  * Halpern et al. (2001a) Halpern J. P., Camilo F., Gotthelf E. V., Helfand D. J., Kramer M., Lyne A. G., Leighly K. M., Eracleous M., 2001a, ApJ, 552, L125
  * Halpern et al. (2012) Halpern J. P., Gotthelf E. V., Camilo F., 2012, ApJ, 753, L14
  * Halpern et al. (2001b) Halpern J. P., Gotthelf E. V., Leighly K. M., Helfand D. J., 2001b, ApJ, 547, 323
  * Helfand et al. (2007) Helfand D. J., Gotthelf E. V., Halpern J. P., Camilo F., Semler D. R., Becker R. H., White R. L., 2007, ApJ, 665, 1297
  * Jones (1968) Jones F. C., 1968, Physical Review, 167, 1159
  * Kashiyama (2017) Kashiyama K., 2017, in Astrophysics and Space Science Library, Vol. 446, Modelling Pulsar Wind Nebulae, Torres D. F., ed., p. 279
  * Kennel & Coroniti (1984) Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 694
  * Khangulyan et al. (2014) Khangulyan D., Aharonian F. A., Kelner S. R., 2014, ApJ, 783, 100
  * Kirk et al. (2009) Kirk J. G., Lyubarsky Y., Petri J., 2009, in Astrophysics and Space Science Library, Vol. 357, Astrophysics and Space Science Library, Becker W., ed., p. 421
  * Kothes et al. (2001) Kothes R., Uyaniker B., Pineault S., 2001, ApJ, 560, 236
  * Lemiere et al. (2009) Lemiere A., Slane P., Gaensler B. M., Murray S., 2009, ApJ, 706, 1269
  * Li et al. (2010) Li H., Chen Y., Zhang L., 2010, MNRAS, 408, L80
  * Li et al. (2021) Li J. et al., 2021, ApJ, 913, L33
  * Liu & Wang (2021) Liu R.-Y., Wang X.-Y., 2021, ApJ, 922, 221
  * López-Coto et al. (2022) López-Coto R., de Oña Wilhelmi E., Aharonian F., Amato E., Hinton J., 2022, Nature Astronomy, 6, 199
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Mares et al. (2021) Mares A. et al., 2021, ApJ, 912, 158
  * Martín et al. (2016) Martín J., Torres D. F., Pedaletti G., 2016, MNRAS, 459, 3868
  * Martín et al. (2012) Martín J., Torres D. F., Rea N., 2012, MNRAS, 427, 415
  * Mattana et al. (2009) Mattana F. et al., 2009, ApJ, 694, 12
  * Moderski et al. (2005) Moderski R., Sikora M., Coppi P. S., Aharonian F., 2005, MNRAS, 363, 954
  * Moskalenko et al. (2006) Moskalenko I. V., Porter T. A., Strong A. W., 2006, ApJ, 640, L155
  * Pacini & Salvati (1973) Pacini F., Salvati M., 1973, ApJ, 186, 249
  * Pineault & Joncas (2000) Pineault S., Joncas G., 2000, AJ, 120, 3218
  * Popescu et al. (2017) Popescu C. C., Yang R., Tuffs R. J., Natale G., Rushton M., Aharonian F., 2017, MNRAS, 470, 2539
  * Porter et al. (2006) Porter T. A., Moskalenko I. V., Strong A. W., 2006, ApJ, 648, L29
  * Rees & Gunn (1974) Rees M. J., Gunn J. E., 1974, MNRAS, 167, 1
  * Reynolds & Chevalier (1984a) Reynolds S. P., Chevalier R. A., 1984a, ApJ, 278, 630
  * Reynolds & Chevalier (1984b) Reynolds S. P., Chevalier R. A., 1984b, ApJ, 278, 630
  * Reynolds et al. (2017) Reynolds S. P., Pavlov G. G., Kargaltsev O., Klingler N., Renaud M., Mereghetti S., 2017, Space Sci. Rev., 207, 175
  * Schlickeiser & Ruppel (2010) Schlickeiser R., Ruppel J., 2010, New Journal of Physics, 12, 033044
  * Shu (1992) Shu F. H., 1992, The physics of astrophysics. Volume II: Gas dynamics.
  * Slane et al. (2010) Slane P., Castro D., Funk S., Uchiyama Y., Lemiere A., Gelfand J. D., Lemoine-Goumard M., 2010, ApJ, 720, 266
  * Supan et al. (2016) Supan L., Supanitsky A. D., Castelletti G., 2016, A&A, 589, A51
  * Sushch et al. (2017) Sushch I., Oya I., Schwanke U., Johnston S., Dalton M. L., 2017, A&A, 605, A115
  * Tanaka & Asano (2017) Tanaka S. J., Asano K., 2017, ApJ, 841, 78
  * Tanaka & Takahara (2010) Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248
  * Tanaka & Takahara (2011) Tanaka S. J., Takahara F., 2011, ApJ, 741, 40
  * Tanaka & Takahara (2013) Tanaka S. J., Takahara F., 2013, MNRAS, 429, 2945
  * Tang & Chevalier (2012) Tang X., Chevalier R. A., 2012, ApJ, 752, 83
  * Tang et al. (2015) Tang Y., Yang C., Zhang L., Wang J., 2015, ApJ, 812, 32
  * The LHAASO Collaboration (2021) The LHAASO Collaboration, 2021, Science, 373, 425
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration et al., 2021, Nature Astronomy, 5, 460
  * Torres (2017) Torres D. F., ed., 2017, Astrophysics and Space Science Library, Vol. 446, Modelling Pulsar Wind Nebulae. Springer International Publishing
  * Torres et al. (2014) Torres D. F., Cillis A., Martín J., de Oña Wilhelmi E., 2014, Journal of High Energy Astrophysics, 1, 31
  * Torres et al. (2013) Torres D. F., Cillis A. N., Martín Rodriguez J., 2013, ApJ, 763, L4
  * Truelove & McKee (1999) Truelove J. K., McKee C. F., 1999, ApJS, 120, 299
  * Ubertini et al. (2005) Ubertini P. et al., 2005, ApJ, 629, L109
  * van der Swaluw et al. (2001) van der Swaluw E., Achterberg A., Gallant Y. A., Tóth G., 2001, A&A, 380, 309
  * Volpi et al. (2008) Volpi D., Del Zanna L., Amato E., Bucciantini N., 2008, A&A, 485, 337
  * Vorster et al. (2013) Vorster M. J., Tibolla O., Ferreira S. E. S., Kaufmann S., 2013, ApJ, 773, 139
  * Whiteoak & Green (1996) Whiteoak J. B. Z., Green A. J., 1996, A&AS, 118, 329
  * Xin et al. (2019) Xin Y., Zeng H., Liu S., Fan Y., Wei D., 2019, ApJ, 885, 162
  * Xin et al. (2018) Xin Y.-L., Liao N.-H., Guo X.-L., Yuan Q., Liu S.-M., Fan Y.-Z., Wei D.-M., 2018, ApJ, 867, 55
  * Yao et al. (2017) Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29
  * Yu et al. (2022) Yu H., Wu K., Wen L., Fang J., 2022, New A, 90, 101669
  * Zabalza (2015) Zabalza V., 2015, Proc. of International Cosmic Ray Conference 2015, 922
  * Zhang et al. (2008) Zhang L., Chen S. B., Fang J., 2008, ApJ, 676, 1210
  * Zhang & Yang (2009) Zhang L., Yang X. C., 2009, ApJ, 699, L153
  * Zhu et al. (2018) Zhu B.-T., Zhang L., Fang J., 2018, A&A, 609, A110
  * Zhu & Huang (2014) Zhu J., Huang M., 2014, A&A, 564, A111