# A 3D Diffusive and Advective Model of Electron Transport Applied to the
Pulsar Wind Nebula HESS J1825-137

T. Collins,1, G. Rowell1, S. Einecke1, F. Voisin1, Y. Fukui2 and H. Sano3
1School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia
2Department of Physics, University of Nagoya, Furo-cho, Chikusa-ku, Nagoya,
464-8601, Japan
3Faculty of Engineering, Gifu University, Yanagido 1-1, Gifu, 501-1193, Japan
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

HESS J1825-137 is one of the most powerful and luminous TeV gamma-ray pulsar
wind nebulae (PWNe), making it an excellent laboratory to study particle
transportation around pulsars. We present a model of the (diffusive and
advective) transport and radiative losses of electrons from the pulsar PSR
J1826-1334 powering HESS J1825-137 using interstellar medium gas (ISM) data,
soft photon fields and a spatially varying magnetic field. We find that for
the characteristic age of $21\,$\mathrm{kyr}$$, PSR J1826-1334 is unable to
meet the energy requirements to match the observed X-ray and gamma-ray
emission. An older age of $40\,$\mathrm{kyr}$$, together with an electron
conversion efficiency of $0.14$ and advective flow of $v=0.002c$, can
reproduce the observed multi-wavelengh emission towards HESS J1825-137. A
turbulent ISM with magnetic field of
$B=$20\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{G}60\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$
to the north of HESS J1825-137 (as suggested by ISM observations) is required
to prevent significant gamma-ray contamination towards the northern
$\mathrm{TeV}$ source HESS J1826-130.

###### keywords:

ISM: cosmic rays - ISM: evolution - gamma-rays: general - X-rays: general -
ISM individual (HESS J1825-137) - pulsars: individual (PSR J11826-1334)

††pubyear: 2023††pagerange: A 3D Diffusive and Advective Model of Electron
Transport Applied to the Pulsar Wind Nebula HESS J1825-137–E

##  1 Introduction

HESS J1825-137 is a luminous pulsar wind nebula (PWN) powered by pulsar PSR
J1826-1334 with spin-down power $\dot{E}=2.8\times
10^{36}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and characteristic age
$\tau_{c}=P/2\dot{P}=21.4\,$\mathrm{kyr}$$ (Manchester et al., 2005). The
distance to PSR J11826-1334 has been estimated to lie at $3.6\,$\mathrm{kpc}$$
based on dispersion measurements (Taylor & Cordes, 1993; Cordes & Lazio,
2002), however we will use a distance of $4\,$\mathrm{kpc}$$ in line with Van
Etten & Romani (2011) and H.E.S.S. Collaboration et al. (2019). The
$\mathrm{TeV}$ gamma-ray emission from HESS J1825-137 has a characteristic
($1/e$) radius of $$$\pm$$_{\text{stat}}\pm$$_{\text{sys}}$, implying a radius
of $\approx 46\,$\mathrm{pc}$$ based on a distance of $4\,$\mathrm{kpc}$$
(H.E.S.S. Collaboration et al., 2019). Owing to its brightness in TeV gamma
rays, HESS J1825-137 is an ideal laboratory to study relativistic particle
transport in and around middle-aged PWNe. Several studies (e.g. Porth et al.
(2016), Giacinti et al. (2020)) suggest that both diffusive and advective
transport mechanisms are required to explain the extended gamma-ray morphology
towards PWNe.

Situated $None$ north of HESS J1825-137 (see Figure 1), HESS J1826-130 is a
TeV gamma-ray source and possible accelerator of cosmic rays up to
$\mathrm{PeV}$energies (Abeysekara et al., 2020; Cao et al., 2021). Due to its
close proximity to HESS J1825-137, HESS J1826-130 was originally considered an
extension of HESS J1825-137 until it was revealed to be a separate source of
gamma rays (H.E.S.S. Collaboration et al., 2018a). The two nearby supernova
remnants (SNRs) SNR G018.1-0.1 and SNR G018.6-0.2 (Odegard, 1986; Brogan et
al., 2006) were deemed to be unlikely to be associated with HESS J1826-137 due
to their offset positions and small angular diameters (H.E.S.S. Collaboration
et al., 2020). Instead, the Eel PWN (PWN G18.5-0.4) and PSR J11826-1256 are
associated with HESS J1826-130 based on spatial coincidence (H.E.S.S.
Collaboration et al., 2018a). PSR J11826-1256 has a spin-down power of
$3.6\times 10^{36}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and
characteristic age of $14\,$\mathrm{kyr}$$, well within the range of pulsar
properties associated with $\mathrm{TeV}$ PWNe (Manchester et al., 2005;
H.E.S.S. Collaboration et al., 2018b).

Araya et al. (2019) revealed $\mathrm{GeV}$ gamma-ray emission $None$ to the
Galactic south of HESS J1825-137. The same study postulated that the
$\mathrm{GeV}$ emission from this region originates from cosmic rays
accelerated by the SNR or PWN associated with HESS J1825-137 or a star-forming
region such as the Cygnus Cocoon. Comprehensive modelling of the spectral
energy distribution (SED) towards the $\mathrm{GeV}$ region suggests that the
emission may be reflective of an earlier epoch of the PWN or a combination of
HESS J1825-137 and nearby compact object LS 5039 (Collins et al., 2021).

The PWN associated with HESS J1825-137 must be expanding within the progenitor
SNR. A large H$\alpha$ rim-like structure discovered by Stupar et al. (2008)
is present towards the south of HESS J1825-137. Voisin et al. (2016)
postulated a connection between this rim and another southern H$\alpha$ rim
and the progenitor SNR of HESS J1825-137. Both structures lie
$\approx\\!\\!$$$ away from PSR J1826-1334 ($\approx 120\,$\mathrm{pc}$$ for a
distance of $4\,$\mathrm{kpc}$$), which is consistent with the predicted SNR
radius of $130\,$\mathrm{pc}$$ as suggested by de Jager & Djannati-Ataï
(2009).

Electrons released by a pulsar are subject to varying transport processes such
as diffusion and/or advection. It has been proposed that advection dominates
the particle transport close to the pulsar while diffusion dominates the outer
reaches of the nebula (Tang & Chevalier, 2012; Porth et al., 2016). TeV halos
around PWN have been suggested to form when electrons escape the PWN into the
surrounding ISM where diffusion dominates particle transport (Giacinti et al.,
2020; Recchia et al., 2021). H.E.S.S. Collaboration et al. (2019) found that
the energy-dependent radial extent of the $\mathrm{TeV}$ PWN associated with
HESS J1825-137 is unlikely to be explained with a diffusion-only scenario and
requires an overall bulk flow towards lower Galactic longitudes.

Kennel & Coroniti (1984a, b) developed the first 1D magnetohydrodynamic model
of the Crab Nebula as an extension to the model first developed by Rees & Gunn
(1974) by considering the PWNe evolving in a slowly expanding SNR shell.
Khangulyan et al. (2018) applied this approach to HESS J1825-137 and was able
to reproduce the size of the PWN and the position of the termination shock
($r_{\text{ts}}\approx 0.03\,$\mathrm{pc}$$) by assuming a short initial
period of PSR J11826-1334 ($P\approx 1\,$\mathrm{ms}$$), small braking index
($n\leq 2$), birth spin-down power $\geq
10^{41}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and evolution in dense
environment ($n_{\text{ISM}}\geq 1\,${\mathrm{cm}}^{-3}$$). This is supported
by the presence of dense molecular clouds towards HESS J1825-137 as described
by Voisin et al. (2016). Alternatively,Van Etten & Romani (2011) treated the
transport of electrons from PSR J11826-1334 as a series of uniform, spherical
‘bubbles’ to study the inverse Compton and Synchrotron emission from PWN (e.g.
(Aharonian et al., 1997). They found that a combined diffusive and advective
model was able to predict the multi-wavelengh SED and radial profile towards
HESS J1825-137. Recently, Lu et al. (2023) investigated the gamma-ray emission
towards HESS J1825-137 by combining a 1D diffusion-advection particle
transport model with Markov chain Monte Carlo techniques to obtain the best-
fitting parameters. However, the observed asymmetric gamma-ray morphology
observed towards HESS J1825-137 suggests a similarly asymmetric electron
density and/or magnetic field. Moreover, dense molecular clouds would prohibit
the escape of electrons out of the PWN, leading to an irregular electron
number density distribution around the pulsar, which cannot be predicted by a
1D model and therefore requires a more complex model.

Figure 1: HESS excess counts towards HESS J1825-137 (H.E.S.S. Collaboration et
al., 2019) overlayed by the regions used to extract the gamma-ray spectra
towards HESS J1825-137 (black) and HESS J1826-130 (red). 4FGL J1824.5-1351e,
eHWC J1825-134 and LHAASO J1825-1326 are shown by the purple, brown and cyan
dashed circles respectively with the positions of PSR J11826-1334 (blue) and
PSR J11826-1256 (black). The position of the nearby binary system LS 5039 is
indicated by the black arrow.

The evolution of the cosmic-ray number density distribution can be described
by the Fokker-Planck equation for particle transport (e.g. Skilling (1975);
Cesarsky & Volk (1978)). Henceforth, this equation will be called the
transport equation. Analytical solutions of the transport equation can be
found for specific cases, e.g. isotropic diffusion in a homogeneous
environment (Blumenthal & Gould, 1970; Atoyan et al., 1995; Aharonian &
Atoyan, 1996). However, it can only be solved numerically for more complex
systems, e.g. anisotropic diffusion where the diffusion coefficient varies
with position.

In this paper, we present a model that incorporates 3D distributions of the
ISM hydrogen number density and magnetic field and solves the transport
equation numerically. This model assumes PSR J11826-1334 to be a source of
high-energy electrons and aims at reproducing the X-ray and gamma-ray
morphology, spectrum and surface brightness radial profiles towards HESS
J1825-137.

##  2 Particle Transport and Multi-wavelength Emission

Upon the release from an accelerator, such as a SNR or PWN, cosmic rays are
transported through the ISM and experience radiative losses. The evolution of
the number density distribution of cosmic rays, $n\equiv
n\quantity(\gamma,t,\vec{r})$, with Lorentz factor
$\gamma\equiv\gamma\quantity(\vec{r})$, at position $\vec{r}\equiv(x,y,z)$ and
time $t$ after the birth of the accelerator, can be described by (e.g.
Skilling (1975); Cesarsky & Volk (1978)):

$\displaystyle\partialderivative{n}{t}=$
$\displaystyle\partialderivative{\gamma}\quantity(\dot{\gamma}n)+\nabla\quantity(\bar{\bar{D}}\cdot\nabla
n)-\nabla\cdot\quantity(n\vec{v}_{A})-\frac{1}{3}\partialderivative{\gamma}\quantity(\gamma\quantity(\nabla\cdot\vec{v}_{A}))n$
(1)
$\displaystyle+\partialderivative{\gamma}\quantity(\gamma^{2}D_{\gamma\gamma}\partialderivative{\gamma}\quantity(\frac{n}{\gamma^{2}}))+S\quantity(\gamma,t,\vec{r})\text{
.}$

The first term in Equation 1 gives the evolution of cosmic-ray density due to
radiative losses. The second term considers the spatial diffusion of cosmic
rays as a second-rank tensor
($\bar{\bar{D}}\equiv\bar{\bar{D}}\quantity(\gamma,t,\vec{r})$), allowing
preferential direction of transport. The third term describes the evolution of
cosmic-ray density due to advection as a co-moving fluid with velocity
$\vec{v}_{A}\equiv\vec{v}_{A}\quantity(\gamma,t,\vec{r})$. The fourth term
considers losses due to adiabatic expansion. The fifth term represents the re-
acceleration of cosmic rays due to stochastic processes with
$D_{\gamma\gamma}$ being the acceleration rate. Finally,
$S\quantity(\gamma,t,\vec{r})$ is the cosmic-ray source/injection function.

To numerically solve Equation 1, explicit finite difference techniques forward
in time can be used after discretising a region of interest into a grid of
voxels with dimension $\Delta x\Delta y\Delta z$ and time step $\Delta t$:

$\displaystyle\frac{n\quantity(\gamma,t+\Delta
t,\vec{r})-n\quantity(\gamma+\Delta\gamma,t,\vec{r})}{\Delta t}$
$\displaystyle=\quantity(\partialderivative{n}{t})^{\prime}_{\text{diff}}+\quantity(\partialderivative{n}{t})^{\prime}_{\text{adv}}+\quantity(\partialderivative{n}{t})^{\prime}_{\text{adb}}$
(2) $\displaystyle+\quantity(\partialderivative{n}{t})^{\prime}_{\text{re-
acc}}+S\quantity(\gamma,t,\vec{r})\text{ ,}$

where prime represents the evolution of the number density distribution after
radiative losses. The following discussion describes how the implemented model
treats individual terms in Equation 2.

###  2.1 Radiation Losses

High-energy electrons interact with the ISM via inverse Compton interactions
on ambient photons, via Bremsstrahlung with interstellar gas and via
synchrotron interactions against magnetic fields (see Appendix A). The
evolution of the electron number density with Lorentz factor $\gamma$ due to
radiative losses is given by:

$\displaystyle\partialderivative{n}{t}=\partialderivative{\gamma}\quantity(\dot{\gamma}n)\text{
,}$ (3)

where $\dot{\gamma}$ is the cooling rate as given by Manolakou et al. (2007):

$\begin{aligned}
\dot{\gamma}=b_{s}\gamma^{2}+b_{c}\quantity(3\ln\gamma+18.8)+5.3b_{b}+\sum\limits_{j}b_{\text{IC}}^{j}\gamma^{2}F_{\text{KN}}^{j}\quantity(\gamma)\end{aligned}\text{
,}$ (4)

for the case of ionisation or Bremsstrahlung losses in neutral hydrogen. Here,
$j$ sums over all radiation fields (CMB, infrared and optical photons),
$b_{s}$, $b_{c}$, $b_{b}$ and $b_{IC}$ are the coefficients for synchrotron
losses, Coulomb losses, Bremsstrahlung losses and inverse Compton losses
respectively and $F_{\text{KN}}$ is the Klein-Nishina cross section (see
Equation 20). The photon fields were assumed to be constant across the 3D
grid. The general solution to Equation 3 is:

$\displaystyle n\quantity(\gamma,t+\Delta t)$
$\displaystyle=\frac{\dot{\gamma}_{0}}{\dot{\gamma}}n\quantity(\gamma_{0},t)\text{
,}$ (5)

where $\gamma_{0}\equiv\gamma_{0}\quantity(x,y,z)=\gamma+\Delta\gamma$ is the
Lorentz factor at time $t$ before electrons cool to Lorentz factor $\gamma$ at
time $t+\Delta t$.

###  2.2 Diffusion

Over distances smaller than the gyro-radius, $r_{g}$, electrons propagate
through the ISM via ballistic motion. In a medium with randomised magnetic
turbulence ($\delta B$), electrons scatter and the motion switches to a
diffusive regime for distances larger than the gyro-radius (e.g. Prosekin et
al. (2015)). For a simple case of isotropic diffusion in magnetic field
$B\quantity(\vec{r})$, the diffusion tensor in Equation 1 becomes a scalar;
$\overline{\overline{D}}\rightarrow D\quantity(E,\vec{r})$, where $E$ is the
energy of the cosmic ray.

Suppression of cosmic-ray diffusion (compared to the Galactic average) is to
be expected towards PWNe and SNRs where magnetic field turbulence is enhanced
and the diffusion coefficient, $D\quantity(E,\vec{r})$, can be parameterised
by (e.g. Gabici et al. (2007)):

$\displaystyle D\quantity(E,\vec{r})$ $\displaystyle=\chi
D_{0}\quantity(\frac{E/$\mathrm{GeV}$}{B\quantity(\vec{r})/3\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$})^{\delta}\text{
,}$ (6)

where $D_{0}=3\times 10^{27}\,${\mathrm{cm}}^{2}\text{\,}{\mathrm{s}}^{-1}$$
is the average Galactic diffusion coefficient at $1\,$\mathrm{GeV}$$,
$\delta=0.5$ following cosmic-ray observations (e.g. see Strong et al. (2007))
and the diffusion suppression factor, $\chi$, takes values $\leq 1$ depending
on the environment (Berezinskii et al., 1990). For example, Gabici et al.
(2007) found that highly suppressed diffusion ($\chi\sim 0.01$) in molecular
clouds can significantly affect the shape of the observed gamma-ray spectrum.
However, the diffusion suppression factor is not well constrained and a
variety of $\chi$ have been found, e.g. Li & Chen (2010), Giuliani et al.
(2010) and Gabici et al. (2010) found values of $\chi=0.1$, $0.01$ and $0.06$
towards SNR W28, respectively. Similarly, Protheroe et al. (2008) showed that
the suppression factor towards the star-forming region Sgr B2 takes values
$<0.02$ based on the radio synchotron flux. Lu et al. (2023) found a diffusion
coefficient of $1.4\times
10^{26}\,${\mathrm{cm}}^{2}\text{\,}{\mathrm{s}}^{-1}$$ at $1\,$\mathrm{GeV}$$
towards HESS J1825-137.

Assuming isotropic inhomogeneous diffusion, the diffusive component of
Equation 2 is given by:

$\displaystyle\quantity(\partialderivative{n}{t})_{\text{diff}}$
$\displaystyle=\frac{1}{\Delta i^{2}}\sum_{i=x,y,z}$ (7)
$\displaystyle\quantity[\frac{D\quantity(\gamma,i+\Delta
i)+D\quantity(\gamma,i)}{2}]\cdot\quantity[n\quantity(\gamma,t,i+\Delta
i)-n\quantity(\gamma,t,i)]$ $\displaystyle+$
$\displaystyle\quantity[\frac{D\quantity(\gamma,i-\Delta
i)+D\quantity(\gamma,i)}{2}]\cdot\quantity[n\quantity(\gamma,t,i-\Delta
i)-n\quantity(\gamma,t,i)]\text{ ,}$

where $D\quantity(\gamma,i)$ is the diffusion coefficient from Equation 6. The
central finite difference technique used in Equation 7 only considers the
transport of electrons to/from the surrounding voxels. If the time step is too
large, electrons travel across more than one voxel and are lost from the
system. The finite difference technique is then said to be numerically
‘unstable’. Using Von Neuman stability analysis (e.g. see Isaacson (1966)),
Equation 7 is stable when:

$\displaystyle\Delta t$ $\displaystyle\leq\frac{\Delta
i^{2}}{2D\quantity(i)}\bigg{|}_{\text{min}}\text{ .}$ (8)

###  2.3 Advection

For simplicity, the velocity due to the bulk flow of electrons
($\vec{v}_{A}\equiv[v_{\text{A,x}},v_{\text{A,y}},v_{\text{A,z}}]$) was
assumed to be spatially-independent and energy-independent across the region
of interest. This assumption is reasonable for a model describing HESS
J1825-137 (see section 3) as (H.E.S.S. Collaboration et al., 2019) implied an
overall bulk motion towards lower Galactic longitudes. Using explicit finite
difference techniques, the advective component of Equation 2 is given by
Equation 9:

$\displaystyle\quantity(\partialderivative{n}{t})_{\text{adv}}$
$\displaystyle=-\sum_{i=x,y,z}v_{\text{A,i}}\frac{1}{\Delta
i}\begin{cases}n\quantity(\gamma,t,i+\Delta
i)-n\quantity(\gamma,t,i)\text{,}&v_{\text{A,i}}<0\\\
n\quantity(\gamma,t,i)-n\quantity(\gamma,t,i-\Delta
i)\text{,}&v_{\text{A,i}}>0\end{cases}\text{,}$ (9)

where $v_{\text{A,i}}$ is the component of advective velocity in the ith
direction. Equation 9 uses the forward difference method to approximate the
derivative in Equation 2 when $v_{\text{A,i}}<0$ and the backward difference
method when $v_{\text{A,i}}>0$.

For Equation 9 to be numerically stable, the time step must be chosen so that
an electron does not travel across more than one voxel in time $\Delta t$:

$\displaystyle\Delta t\leq\frac{\Delta
i}{\absolutevalue{v_{\text{A,i}}}}\bigg{|}_{\text{min}}\text{ .}$ (10)

The time step must satisfy both Equation 8 and Equation 10 when modelling a
scenario including both diffusion and advection.

###  2.4 Adiabatic Expansion and Re-acceleration of Electrons

The spatially-independent advective velocity assumed in our model results in
zero adiabatic losses in Equation 1 ($\nabla\cdot\vec{v}_{A}=0$). Moreover,
studies such as Tanaka & Takahara (2010) and Porth et al. (2016) who
considered spherically symmetric advection concluded that adiabatic losses are
dominant over radiative losses for electrons $<1\,$\mathrm{TeV}$$ (equivalent
to gamma-ray emission $<20\,$\mathrm{GeV}$$). As we are interested in the VHE
gamma-ray range which is not dominated by adiabatic losses, adiabatic
expansion is not considered here but is left for future work.

The termination shock (TS) of pulsar wind has been proposed as a site for the
re-acceleration of electrons through diffusive shock acceleration (DSA). By
ensuring the voxel width ($\Delta x$,$\Delta y$,$\Delta z$) is larger than the
diameter of the TS ($0.2\,$\mathrm{pc}$$, (Gaensler & Slane, 2006)), electrons
are both injected and re-accelerated within the same voxel. Therefore, the
source term in Equation 2 treats the injected electron spectra as the spectra
obtained after re-acceleration due to the TS. Furthermore, magnetohydrodynamic
models (e.g. Lemoine & Pelletier (2010); Sironi et al. (2015)) suggest that
DSA at the TS is too suppressed for electron acceleration up to energies
responsible for the $\mathrm{TeV}$ emission seen towards PWNe. Hence, the re-
acceleration of electrons is left for future work.

###  2.5 Multi-wavelength Photon Production

The final electron number density distribution was obtained by solving
Equation 2 in discrete time steps $\Delta t$ until the desired age was
reached. Based on the obtained electron number densities, the multi-wavelength
photon emission was derived for each voxel and summed along the line of sight,
$z$, to obtain the 2D photon distribution. Equations 17, 19a and 22 gives the
flux from synchrotron, inverse Compton and Bremsstrahlung interactions
respectively.

##  3 Application to HESS J1825-137

The modelling described in section 2 was applied to the PWN HESS J1825-137
with the pulsar PSR J11826-1334 being the accelerator of high-energy
electrons. PSR J11826-1334 is located at $\ell=$$$ $b=-$$$ and has a proper
motion of $\approx 440\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ (assuming a
distance of $4\,$\mathrm{kpc}$$) approximately perpendicular to the extended
$\mathrm{TeV}$ emission (see Figure 3) (Manchester et al., 2005). Hence, the
proper motion of the pulsar is unlikely to be related to the asymmetric gamma-
ray emission and our model assumed that electrons are injected at the current
position of the pulsar for simplicity. Two different ages of PSR J11826-1334
were considered, the characteristic age of $21.4\,$\mathrm{kyr}$$ and the
older age of $40\,$\mathrm{kyr}$$ suggested by Van Etten & Romani (2011). The
presence of the TeV halo toward HESS J1825-137 indicates that the system is a
middle aged PWN where diffusive particle transport dominates the outer reaches
of the Nebula (Tang & Chevalier, 2012; Porth et al., 2016; Giacinti et al.,
2020; Recchia et al., 2021). However, H.E.S.S. Collaboration et al. (2019)
suggested that both diffusive and advective transport mechanisms are present
in HESS J1825-137.

Each voxel in the 3D grid had a volume of $\Delta x\Delta y\Delta z$, where
$\Delta z$ is the voxel length in the line of sight and $\Delta x$ and $\Delta
y$ are the voxel length along Galactic longitude and latitude respectively.
For the purposes of this study, we utilised a $200\,$\mathrm{pc}$\times
200\,$\mathrm{pc}$\times 200\,$\mathrm{pc}$$ grid consisting of voxels of size
$2\,$\mathrm{pc}$\times 2\,$\mathrm{pc}$\times 2\,$\mathrm{pc}$$
($\approx$$\times$$\times 2\,$\mathrm{pc}$$). The pulsar was located in the
centre of the grid with the central $z$ slice lying at distance
$4\,$\mathrm{kpc}$$. The time step used for the finite difference technique
was $\approx 8\,\text{yr}$.

###  3.1 Electron Injection

High-energy electrons were injected into the 3D grid by PSR J1826-1334 and
follow an exponential cutoff power-law:

$\displaystyle S\quantity(E,t)$
$\displaystyle=A\cdot\quantity(\frac{E}{1\,$\mathrm{TeV}$})^{-\Gamma}\exp(-\frac{E}{E_{c}})\text{
,}$ (11)

following the observed $\mathrm{TeV}$ gamma-ray emission (e.g. see Blumenthal
& Gould (1970)) as observed by H.E.S.S. (H.E.S.S. Collaboration et al., 2019),
where $E_{c}$ is the cutoff energy and $A$ is the normalisation factor such
that:

$\displaystyle L_{\text{inj}}\quantity(t)$
$\displaystyle=\int_{E_{\text{min}}}^{E_{\text{max}}}S\quantity(E,t)\differential{E}\text{
,}$ (12)

with $L_{\text{inj}}\equiv\eta\dot{E}$ being the electron injection
luminosity, $\eta<1$ is the conversion efficiency of the pulsar spin-down
power, $E_{\text{min}}=1\,$\mathrm{MeV}$$ and
$E_{\text{max}}=500\,$\mathrm{TeV}$$. The spin-down power,
$\dot{E}\quantity(t)$ at time $t$ is given by (Haensel et al., 2007):

$\displaystyle\dot{E}\quantity(t)$
$\displaystyle=\dot{E}\quantity(t=t_{\text{age}})\quantity[1+\quantity(n-1)\frac{\dot{P}\quantity(t-t_{\text{age}})}{P}]^{-\Gamma_{\mathscr{n}}}\text{
,}$ (13)

where $n$ is the braking index of the pulsar,
$\Gamma_{n}\equiv\quantity(n+1)/\quantity(n-1)$ and
$\dot{E}\quantity(t=t_{\text{age}})$, $P$ and $\dot{P}$ are the spin-down
power, period and spin-down period of the pulsar at the current age
$t_{\text{age}}$.

###  3.2 The Environment Towards HESS J1825-137

####  3.2.1 Magnetic Field

Figure 2: The radially symmetric magnetic field due to the pulsar at slice
$z=4\,$\mathrm{kpc}$\pm 1\,$\mathrm{pc}$$. Overlaid are the green H.E.S.S.
contours at 4, 5, 6, 7, 8 and 9$\sigma$ (H.E.S.S. Collaboration et al., 2019).

Following Van Etten & Romani (2011), the magnetic field due to the PWN was
assumed to follow a time-independent power-law with a decreasing magnetic
field strength varying with distance $r$ from the pulsar:

$\displaystyle B_{\text{PWN}}\quantity(r)$
$\displaystyle=B_{0}\quantity(\frac{r}{r_{\text{ts}}})^{-\beta}\text{ ,}$ (14)

where $r_{\text{ts}}=0.03\,$\mathrm{pc}$$ is the radius of the termination
shock, and $B_{0}$ and $\beta$ are free parameters optimised to match the
multi-wavelength SED of HESS J1825-137. Van Etten & Romani (2011) suggested
$\beta=-0.69$ and $B_{0}=400\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ for an
age of $40\,$\mathrm{kyr}$$. Note that Van Etten & Romani (2011) considered an
additional dependence on the spin-down energy of the pulsar which was not
considered in this study.

####  3.2.2 Interstellar Medium

Figure 3: (top) Nanten 12CO(1-0) integrated intensity in the velocity range
$40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ corresponding to
$3.5-4.5\,$\mathrm{kpc}$$ overlaid by green H.E.S.S. contours (at 4, 5, 6, 7,
8 and 9$\sigma$). The yellow dot represents the Suzaku region A as defined in
Uchiyama et al. (2009) and is used to extract the X-ray SED. The region used
to obtain the gamma-ray spectra towards HESS J1825-137 and HESS J1826-130 are
shown in black and red respectively. (bottom) Calculated ISM number density
across the 3D grid at $4\,$\mathrm{kpc}$\pm 1\,$\mathrm{pc}$$ where the voxels
within the PWN extent ($R<$$$) are set to a density of
$0.5\,${\mathrm{cm}}^{-3}$$ to represent the bubble that has been swept out by
the stellar wind from the progenitor star. The proper motion of the pulsar is
shown by the red arrow with the projected birthplaces indicated by the red-
dashed line. The projected birthplaces for ages $21\,$\mathrm{kyr}$$ and
$40\,$\mathrm{kyr}$$ are indicated by the magenta and brown dot respectively.
The direction of advective particle transport suggested by H.E.S.S.
Collaboration et al. (2019) is shown by the green arrow. The black and green
segments represent the regions used to extract the X-ray and gamma-ray surface
brightness radial profile respectively. Molecular clouds R1-R5 from Voisin et
al. (2016) are shown in yellow with the position of PSR J11826-1256 indicated
by the red star. The width and height of the voxel in the 2D slice is
$2\,$\mathrm{pc}$$ compared to the minimum Nanten resolution of
$5\,$\mathrm{pc}$$ (assuming a distance of $4\,$\mathrm{kpc}$$).

The Nanten 12CO(1-0) survey (Mizuno & Fukui, 2004) was used to trace the
column density of molecular hydrogen towards HESS J1825-137:

$\displaystyle N_{H_{2}}=X_{\mathrm{12CO}}W_{\mathrm{12CO}}$ (15)

where $W_{\mathrm{12CO}}$ is the integrated intensity of the gas. The scaling
factor $X_{\mathrm{12CO}}=1.5\times
10^{20}\,${\mathrm{cm}}^{-2}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{km}}^{-1}\text{\,}\mathrm{s}$$
is assumed to be constant over the Galactic plane but may vary with
galactocentric radius ($1.3-1.5$ per $\mathrm{kpc}$) (Strong et al., 2004).
The length of the 3D grid ($<1\,$\mathrm{kpc}$$) allows the assumption of a
constant $X_{\mathrm{12CO}}$ towards the region of interest. PSR J1826-1334
has a dispersion measure distance of $4\,$\mathrm{kpc}$$, corresponding to a
velocity of $50\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ using the Galactic
Rotation model (Brand & Blitz, 1993). As there may be local motion of the gas
unrelated to Galactic rotation, we considered a velocity range of
$40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
($3.5\text{\,}\mathrm{kpc}4.5\text{\,}\mathrm{kpc}$) consistent with Voisin et
al. (2016). Atomic hydrogen in the same velocity range contributes less than
$1\%$ to the total column density towards HESS J1825-137 and thus was not
considered (Voisin et al., 2016; Collins et al., 2021). Assuming that all the
gas in the $40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ velocity range
lies within the 3D grid and the density along the line of sight is constant,
the number density of a voxel with column density $N_{H_{2}}$ is given by:

$\displaystyle n_{H}$ $\displaystyle=\frac{N_{H}}{200\,$\mathrm{pc}$}\text{
,}$ (16)

where $N_{H}\equiv 2.8N_{H_{2}}$ considers a $20\%$ He component.

Stellar winds from the progenitor star of PSR J11826-1334 pushes out gas in
the nearby vicinity (Castor et al., 1975). The subsequent supernova explosion
creates a ‘bubble’ of hot dense gas around a low-density interior. A region of
low-density gas in the $40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
velocity range can be seen towards the centre of the $\mathrm{TeV}$ emission
in Figure 3. To include this, any voxels lying within the extent of the PWN
volume (a sphere centered on the pulsar with radius $$$\approx
35\,$\mathrm{pc}$$) was set to a density of $0.5\,${\mathrm{cm}}^{-3}$$ based
on the average densities expected within massive stellar wind bubbles (Weaver
et al., 1977). The Nanten 12CO(1-0) integrated intensity between
$40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and calculated ISM number
density for the central slice lying at distance $4\,$\mathrm{kpc}$$ can be
seen in Figure 3. Any difference between the bottom and top panel of Figure 3
was due to the different resolutions of the 3D grid and Nanten.

Turbulent motion in the ISM results in an amplification of the magnetic field,
suppressing the diffusion of electrons as they travel through the ISM as given
by Equation 6. Figure 6 of Voisin et al. (2016) shows a three-coloured image
of the CS(1-0) and $\text{NH}_{3}$ integrated intensity between
$40-60\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and the $\text{H62}\alpha$
integrated intensity between $45-65\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
towards the cloud defined as R1 (see Figure 3 for the position of clouds R1-R5
from Voisin et al. (2016)). This suggests that cloud R1 is highly turbulent
with a minimum magnetic field strength of
$21\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ based on the density of
$600\,${\mathrm{cm}}^{-3}$$ calculated by Voisin et al. (2016) (see Equation
24). The amplification of the magnetic field towards cloud R1 was considered
in subsubsection 3.4.3. Given the likely physical proximity to HESS J1825-137,
cloud R1 may act as a barrier for electrons escaping into HESS J1826-130 from
PSR J1826-1334 (Voisin et al., 2016).

####  3.2.3 Soft Photon Fields

The photon fields around HESS J1825-137 was estimated utilising the radiation
field model described by Popescu et al. (2017); the far-infrared field (FIR)
with temperature $T=40$\mathrm{K}$$ and energy density
$U=1\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$, near infrared field (NIR)
with temperature $T=500\,$\mathrm{K}$$ and energy density
$0.4\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$ and optical light with
temperature $T=3500$\mathrm{K}$$ and energy density of
$U=1.9\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$.

###  3.3 Multi-wavelength Observations

The modelled gamma-ray SED of HESS J1825-137 was optimised to the
$\mathrm{TeV}$ gamma-ray energy flux presented by H.E.S.S. Collaboration et
al. (2019) and the $\mathrm{GeV}$ spectrum from the 4FGL catalogue (Abdollahi
et al., 2020). To compare the modelled surface brightness radial profile to
Figure 6 from H.E.S.S. Collaboration et al. (2019), a collection area of
$0.25\,${\mathrm{km}}^{2}$$ (Benbow, 2005) and observation time of
$387\,$\mathrm{h}\mathrm{r}$$ was used. The X-ray SED and surface brightness
radial profile was optimised to the results presented by Uchiyama et al.
(2009) using a collection area of $0.029\,${\mathrm{m}}^{2}$$.

To investigate the gamma-ray contamination of HESS J18260-130, by HESS
J1825-137, we utilised the gamma-ray SED presented by H.E.S.S. Collaboration
et al. (2020) and the spectrum from the 4FGL catalogue (Abdollahi et al.,
2020). H.E.S.S. Collaboration et al. (2020) estimated that the gamma-ray
contamination to be $40\%$ for photon energies below $1.5\,$\mathrm{TeV}$$ and
$20\%$ above $1.5\,$\mathrm{TeV}$$. The modelled X-ray emission towards HESS
J1826-130 was constrained by the ROSAT X-ray upper limit calculated using the
ROSAT X-ray background tool (Sabol & Snowden, 2019). The regions used to
extract the X-ray and gamma-ray SED towards HESS J1825-137 and HESS J1826-130
are shown in Figure 3.

###  3.4 Results

Table 1: Model parameters used for the application towards HESS J1825-137.
Fixed parameters refer to those constrained by measurements and non-fixed
refers to those that are optimised to observations discussed in subsection
3.3.

Fixed Parameters | Value | Reference
---|---|---
$t$ |  $21\,$\mathrm{kyr}$$ & $40\,$\mathrm{kyr}$$ | a, b
$d$ | $4\,$\mathrm{kpc}$$ | a
$P$ | $101\,$\mathrm{ms}$$ | a
$\dot{P}$ | $7.5\times 10^{-14}\,$\mathrm{s}\text{\,}{\mathrm{s}}^{-1}$$ | a
$\dot{E}$ | $2.8\times 10^{36}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ | a
$\Delta x$, $\Delta y$, $\Delta z$ | $2\,$\mathrm{pc}$$ |
$\Delta t$ | $8\,$\mathrm{y}\mathrm{r}$$ |
$E_{\text{min}}$ | $1\,$\mathrm{MeV}$$ |
$E_{\text{max}}$ | $500\,$\mathrm{TeV}$$ |
$D_{0}$ | $3\times 10^{27}\,${\mathrm{cm}}^{2}\text{\,}{\mathrm{s}}^{-1}$$ | c
$r_{\text{ts}}$ | $0.03\,$\mathrm{pc}$$ |
$U_{\text{CMB}}$, $T_{\text{CMB}}$ |  $0.26\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$, $2.72\,$\mathrm{K}$$ | e
$U_{\text{NIR}}$, $T_{\text{NIR}}$ |  $1\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$, $500\,$\mathrm{K}$$ | e
$U_{\text{FIR}}$, $T_{\text{FIR}}$ |  $0.4\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$, $40\,$\mathrm{K}$$ | e
$U_{\text{Opt}}$, $T_{\text{Opt}}$ |  $1.9\,$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$$, $3500\,$\mathrm{K}$$ | e
Non-fixed Parameters | Value | Reference
$\eta$ | $<1$ |
$\chi$ | $<1$ | c
$\Gamma$ | - |
$E_{c}$ | - |
$B_{0}$ | - |
$\beta$ | - |
$n$ | 2-3 |
$\vec{v}_{A}(\ell,b,z)$ | ($<0.01c,0,0$) | d
$B_{1826}$ | $\star$ |

  * •

$\star$ See subsubsection 3.4.3

  * •

a. Manchester et al. (2005)

  * •

b. Van Etten & Romani (2011)

  * •

c. Berezinskii et al. (1990)

  * •

d. H.E.S.S. Collaboration et al. (2019)

  * •

e. Popescu et al. (2017)

Table 2: Model parameters that match the multi-wavelength SED and gamma-ray
morphology towards HESS J1825-137. See Appendix D for the $10\%$ and $20\%$
systematic variation of parameters.

Parameter | Model 1 ($21\,$\mathrm{kyr}$$) | Model 1 ($40\,$\mathrm{kyr}$$) | Model 2 ($0.002c$) | Model 3∗ ($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$) | Model 3 ($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$)
---|---|---|---|---|---
$\eta$ | $10.7$ | $0.14$ | $0.14$ | $0.14$ | $0.14$
$\chi$ | $0.25$ | $0.1$ | $0.1$ | $0.1$ | $0.1$
$\Gamma$ | $2.0$ | $1.9$ | $1.9$ | $1.9$ | $1.9$
$E_{c}$ ($\mathrm{TeV}$) | $40$ | $500$ | $500$ | $500$ | $500$
$B_{0}$ ($\mathrm{\SIUnitSymbolMicro}\mathrm{G}$) | $70$ | $450$ | $450$ | $450$ | $450$
$\beta$ | $-0.9$ | $-0.7$ | $-0.7$ | $-0.7$ | $-0.7$
$n$ | $2$ | $2$ | $2$ | $2$ | $2$
$v_{A}$ | - | - | $0.002c$ | - | $0.002c$
$B_{\text{J1826}}$ ($\mathrm{\SIUnitSymbolMicro}\mathrm{G}$) | - | - | - | $60$ |  $60$***1
$\dot{E}_{\text{birth}}$ ($\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$) | $2.1\times 10^{37}$ | $1.1\times 10^{40}$ | $1.1\times 10^{40}$ | $1.1\times 10^{40}$ | $1.1\times 10^{40}$

  * • 11footnotemark: 1

See subsubsection 3.5.3

The full list of model parameters is summarised in Table 1, including any
constraints based on measurements. A computationally quick single-zone model,
where the electron number density is derived using a uniform sphere, was
utilised to investigate a large range of parameters to gain insight into HESS
J1825-137. The results of the single-zone modelling are summarised in
Appendix. C. For an age of $21\,$\mathrm{kyr}$$, the single-zone model
required electrons to follow an exponential cutoff power-law with spectral
index $\Gamma=2.1$ and cutoff $E_{c}=40\,$\mathrm{TeV}$$ to match the observed
gamma-ray SED while an older age of $40\,$\mathrm{kyr}$$ required an index of
$\Gamma=2.1$ and cutoff of $E_{c}=50\,$\mathrm{TeV}$$.

In the following, we present three applications of our model towards HESS
J1825-137. All models incorporated a simple assumption of isotropic diffusion
and radiative losses as described in subsection 2.2. Model 1 considered both
ages of PSR J11862-1334, $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$. Model
2 introduced an additional advective component to Model 1 with velocity
$\vec{v}_{A}=[v_{A},0,0]$ as suggested by H.E.S.S. Collaboration et al. (2019)
to explain the asymmetric gamma-ray morphology. H.E.S.S. Collaboration et al.
(2019) constrained the total flow velocity to be $<0.01c$. Model 3 expanded on
Model 2 by including turbulent ISM towards HESS J1826-130 (see subsubsection
3.2.2) to reduce the contamination by HESS J1825-137. The model parameters
were chosen based on the observations discussed in subsection 3.3 with the
parameter list shown in Table 1. The parameters we found to match the multi-
wavelength SED and morphology are shown in Table 2.

####  3.4.1 Model 1 ($21$ & $40\,$\mathrm{kyr}$$) - Isotropic Diffusion

Figure 4 and 5 show the modelled gamma-ray morphology in different energy
bands, the multi-wavelength SED and the $1-9\,$\mathrm{keV}$$ X-ray and
$0.1-91\,$\mathrm{TeV}$$ gamma-ray surface brightness radial profiles for the
$21$ and $40\,$\mathrm{kyr}$$ models respectively.

Figure 4: Model 1 ($21\,$\mathrm{kyr}$$), see Table 2 for model parameters.
(top & top-middle) Modelled gamma-ray morphology towards HESS J1825-137 in
different energy bands overlaid by green HESS significance contours ($5$, $10$
and $15\sigma$ for $E<10\,$\mathrm{TeV}$$ and $3$, $5$ and $10\sigma$ for
$E>10\,$\mathrm{TeV}$$) and grey $40$, $50$ and $60\sigma$ Nanten
${}^{12}\text{CO}$ integrated intensity contours. The positions of PSR
J11826-1334 and PSR J11826-1256 are indicated by the empty blue and black
crosses respectively. (bottom-middle )$1-9\,$\mathrm{keV}$$ X-ray (left) and
$1-91\,$\mathrm{TeV}$$ gamma-ray (right) surface brightness radial profiles in
comparison to Suzaku (Uchiyama et al., 2009) and HESS observations (H.E.S.S.
Collaboration et al., 2019) respectively. (bottom) SED towards HESS J1825-137
with the orange Suzaku X-ray spectral fit, blue 4FGL J1824.5-135e flux
observations (Abdollahi et al., 2020) and red HESS J1825-137 flux
observations.

Figure 5: Model 1 ($40\,$\mathrm{kyr}$$), see Table 2 for model parameters.
Same panel layout as in Figure 4. Figure 6: SED towards HESS J1826-130 from
electrons accelerated by PSR J1826-1256 by Model 1 ($40\,$\mathrm{kyr}$$)
against the observed flux of HESS J1826-130 (H.E.S.S. Collaboration et al.,
2020) and 4FGL J1826.1-1256 (Abdollahi et al., 2020). The $\mathrm{GeV}$ and
$\mathrm{TeV}$ gamma-ray flux observations towards HESS J1826-130 are
represented by blue and red respectively.

Both models predicted that the gamma-ray morphology towards HESS J1825-137 is
symmetric around the powering pulsar with some gamma-ray contribution
$<1\,$\mathrm{TeV}$$ via Bremsstrahlung radiation toward the region between
HESS J1825-137 and HESS J1826-130 (see Figure 4 and 5). The
$40\,$\mathrm{kyr}$$ gamma-ray emission between $1-10\,$\mathrm{TeV}$$ extends
further from the pulsar than the $21\,$\mathrm{kyr}$$ emission. Both models
predicted a steep surface brightness radial profile for X-rays between
$1-9\,$\mathrm{keV}$$ (see the bottom-middle panels of Figure 4 and 5). The
$21\,$\mathrm{kyr}$$ model was able to replicate the HESS surface brightness
radial profile for gamma rays between $0.133-91\,$\mathrm{TeV}$$ (see the
bottom-middle right panel of Figure 4) while the $40\,$\mathrm{kyr}$$ model
over-predicted the gamma-ray emission for distances $>$$$ from the pulsar (see
the bottom-middle right panel of Figure 5).

The $21\,$\mathrm{kyr}$$ modelled gamma-ray SED predicted by the multi-zone
model was able to match observations with a slight over-prediction ($\approx
94\%$) of the HESS data between $1-10\,$\mathrm{TeV}$$. While able to predict
the normalisation of X-rays produced by synchrotron emission, the model was
unable to replicate the slope of the observed Suzaku SED. The multi-zone
$40\,$\mathrm{kyr}$$ SED was able to predict both the X-ray and gamma-ray SED
with a similar over-prediction in $1-10\,$\mathrm{TeV}$$ photons as seen in
the $21\,$\mathrm{kyr}$$ model. A slight ’bump’ is present in the SED for
photons around $50-100\,$\mathrm{TeV}$$ for both ages.

The $21\,$\mathrm{kyr}$$ model required electrons with spectral index
$\Gamma=2.0$ and cutoff $40\,$\mathrm{TeV}$$ to be injected into the ISM with
a spin-down conversion factor of $10.4$ to match the multi-wavelength SED. The
$40\,$\mathrm{kyr}$$ required a conversion factor of $0.14$ with a spectral
index and cutoff of $1.9$ and $500\,$\mathrm{TeV}$$ respectively. As $\eta<1$,
Models 2 and 3 only considered an age of $40\,$\mathrm{kyr}$$. The initial
birth spin-down power, $\dot{E}_{\text{birth}}$, of the PSR J11826-1334
predicted by the $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$ models are
$2.1\times 10^{37}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and $1.1\times
10^{40}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ respectively.

Figure 6 shows the modelled SED towards HESS J1826-130 due to electrons
escaping from HESS J1825-137 for the $40\,$\mathrm{kyr}$$ model. The SED
towards HESS J1826-130 as a result of HESS J1825-137 exceeds observations for
photons below $2\,$\mathrm{TeV}$$. In the model, too many low-energy electrons
have escaped into the region towards HESS J1826-130 before losing their energy
to radiative losses. It is clear that further refinement of the model is
required to accurately describe the region surrounding HESS J1825-137.

####  3.4.2 Model 2 - Isotropic Diffusion + Advection

The gamma-ray morphology in Figure 5 shows that Model 1 ($40\,$\mathrm{kyr}$$)
did not reproduce the extended $\mathrm{TeV}$ gamma-ray morphology towards
HESS J1825-137 at lower Galactic longitudes (see the top-middle left panel of
Figure 5). Thus, Model 2 introduced an additional advective component as
suggested by H.E.S.S. Collaboration et al. (2019) towards lower Galactic
longitudes. The modelled flux, surface brightness radial profiles and gamma-
ray morphology for Model 2 ($40\,$\mathrm{kyr}$$) with an advective flow of
$v=0.002$ are shown in Figure 7. A comparison between different advection
speeds ($v=0.001c$, $v=0.002c$ and $v=0.003c$) is shown in Figure 17. All
models otherwise have the same parameters as Model 1 ($40\,$\mathrm{kyr}$$)
(see Table 2).

Figure 7: Model 2 ($0.002c$) vs Model 1 ($40\,$\mathrm{kyr}$$), see Table 2
for model parameters. (top and top-middle) panels show the gamma-ray
morphology for Model 2 ($0.002c$). (bottom-middle) $1-9\,$\mathrm{keV}$$ X-ray
(left) and $1-91\,$\mathrm{TeV}$$ gamma-ray (right) surface brightness radial
profiles. (bottom-left) SED towards HESS J1825-137. (bottom-right) SED towards
HESS J1826-130.

Figure 8: The energy flux along Galactic longitude for Model 1
($40\,$\mathrm{kyr}$$, black) and Model 2 ($0.002c$, green) for energy bands
$E<1\,$\mathrm{TeV}$$ (top-left), $1<E<10\,$\mathrm{TeV}$$ (top-right) and
$E>10$\mathrm{TeV}$$ (bottom-left) vs HESS excess counts (H.E.S.S.
Collaboration et al., 2019). (bottom-right) Model 2 ($40\,$\mathrm{kyr}$$)
gamma-ray morphology for energies $<1\,$\mathrm{TeV}$$. The rectangular
regions used to take the slice profile are indicated in black while the region
used to extract the gamma-ray SED towards HESS J1825-137 is shown by the black
circle. The positions of PSR J11826-1334 and PSR J11826-1256 are depicted by
the blue and black empty crosses respectively. See Table 2 for model
parameters.

For further comparison of the gamma-ray morphology towards HESS J1825-137, the
energy flux was extracted from rectangular regions taken along Galactic
longitude centred on PSR J11826-1334 and are shown in Figure 8. An advective
velocity of $0.002c$ was chosen so that the peak in the modelled gamma-ray
morphology in energy range $E<1\,$\mathrm{TeV}$$ and
$1\,$\mathrm{TeV}$<E<10\,$\mathrm{TeV}$$ corresponds to the HESS data (see the
left-upper panel of Figure 17).

While an additional advective flow of $0.002c$ lowered the gamma-ray SED
towards HESS J1826-130 for energies less than $2\,$\mathrm{TeV}$$, the
emission still exceeds H.E.S.S. observations.

####  3.4.3 Model 3 - Isotropic Diffusion + Advection + Magnetic Field towards
HESS J1826-130

As discussed in subsubsection 3.2.2, the turbulent molecular gas between HESS
J1825-137 and HESS J1826-130 can act as a barrier for electrons escaping from
the PWN. As clouds R1-R5 from Voisin et al. (2016) are positioned in an
approximate semi-circle around PSR J11826-1256 (see Figure 3), Model 3
expanded on Model 2 ($0.002c$) by including a shell of increased magnetic
field strength, $B_{1826}$, centred on HESS J1826-130 with inner and outer
radii $0.17\,$\mathrm{\SIUnitSymbolDegree}$$ and
$0.33\,$\mathrm{\SIUnitSymbolDegree}$$ respectively. Model 3∗ refers to Model
1 ($40\,$\mathrm{kyr}$$) with the shell of increased magnetic field strength
with no advective component ($v_{A}=0$).

Figure 9 shows the SED, surface brightness radial profiles and the gamma-ray
flux along Galactic longitude of HESS J1825-137 and HESS J1826-130 for Model 3
with magnetic field strengths of $B=20$, $60$ and
$100\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$. A comparison between Model 1
($40\,$\mathrm{kyr}$$), Model 3∗
($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$) and Model 3
$(60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$)$ is shown in Figure 10 as well
as the gamma-ray morphology for Model 3
($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$). All models otherwise have the
same parameters as Model 1 ($40\,$\mathrm{kyr}$$) (see Table 2).

Figure 9: The energy flux along Galactic longitude profiles (top & top-middle
left), surface brightness radial profiles (top-middle right & bottom-middle
left) and SED towards HESS J1825-137 (bottom-middle right) and HESS J1826-130
(bottom) for Model 3 ($20\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$, black),
($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$, green) and
($100\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$, purple) around HESS
J1826-130. See Table 2 for model parameters.

Figure 10: Comparison of the $40\,$\mathrm{kyr}$$ isotropic diffusion - Model
1 ($40\,$\mathrm{kyr}$$, black), Model 3∗ (Model 1 +
$60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$, green) and Model 3
($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$, purple). The morphology plots
towards HESS J1825-137 for Model 3
($60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$) are shown in the top and top-
middle panels. See Table 2 for model parameters.

###  3.5 Discussion

####  3.5.1 Model 1 - Isotropic Diffusion

The $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$ models were unable to
reproduce both the X-ray and gamma-ray surface brightness radial profiles. For
example, the diffusion suppression coefficient, $\chi$, could be increased to
compensate for the steep X-ray surface brightness radial profile for the
$21\,$\mathrm{kyr}$$ model. Electrons would then escape the PWN at a higher
rate and the gamma-ray surface brightness radial profile will flatten. This
can be seen in the $40\,$\mathrm{kyr}$$ model, which assumed a lower value of
$\chi$ than the $21\,$\mathrm{kyr}$$ model. The shallow $40\,$\mathrm{kyr}$$
gamma-ray surface brightness radial profile indicates that lower energy
electrons have started to accumulate near the pulsar, while high-energy
electrons rapidly lose their energy through radiative cooling and do not
escape far from the pulsar. This is demonstrated in the upper right panel in
Figure 5 where the gamma-ray flux below $10\,$\mathrm{TeV}$$ is relatively
constant over the grid while the flux above $10\,$\mathrm{TeV}$$ is
constrained to the pulsar. This accumulation is not as apparent in the
$40\,$\mathrm{kyr}$$ X-ray surface brightness radial profile and the SED as
the regions used in extracting X-ray spectra and surface brightness radial
profile are smaller than the regions used for the gamma-ray analysis (see
Figure 3 and top-middle right panel of Figure 4). The accumulation of lower-
energy electrons is also a reflected as a bump in the $\mathrm{TeV}$ gamma-ray
SED. The bump occurs when synchrotron losses start to dominate at electron
energies $>9\,$\mathrm{TeV}$$, resulting in inverse Compton radiation
$>6\,$\mathrm{TeV}$$, and radiative energy losses are balanced by the electron
injection luminosity (Manolakou et al., 2007; Hinton & Hofmann, 2009). This
bump is not present for a slightly younger age of $36\,$\mathrm{kyr}$$ (with
the same parameters as Model 1 ($40\,$\mathrm{kyr}$$), see Figure 14), where
the gamma-ray SED $\gtrsim 10\,$\mathrm{TeV}$$ at age $36\,$\mathrm{kyr}$$
matches Model 1 ($40\,$\mathrm{kyr}$$).

The $21\,$\mathrm{kyr}$$ model required a spin-down conversion factor of
$10.7$. To compensate, a braking index of $3$ would inject a greater quantity
of electrons at earlier times (see Equation 13). However, this results in an
accumulation of electrons at lower energies, consequently increasing the
gamma-ray flux for photons with energies $<1$\mathrm{TeV}$$ (see Figure 15)
and the modelled SED no longer reproduces observations. This suggests that the
age of HESS J1825-137 lies between $21\,$\mathrm{kyr}$$ and
$40\,$\mathrm{kyr}$$. The $40\,$\mathrm{kyr}$$ magnetic field profile takes
values of $B_{0}=450\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ and
$\beta=-0.7$ (see Equation 14) in comparison to
$B_{0}=400\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ and $\beta=-0.69$ used in
(Van Etten & Romani, 2011). Van Etten & Romani (2011) considered an evolving
magnetic field $B\propto\dot{E}\quantity(t)$ where the magnetic field takes
larger values at earlier times. This could explain the larger $B_{0}$
normalisation used in our modelling.

####  3.5.2 Model 2 - Isotropic Diffusion + Advection

An advective component of $0.002c$ towards lower Galactic longitudes was
included into Model 1 ($40\,$\mathrm{kyr}$$). The SED and X-ray surface
brightness radial profile with an advective transportation component remains
unchanged to Model 1 ($40\,$\mathrm{kyr}$$). Electrons rapidly escape the
small ($r=$$$) X-ray region, hence the subsequent X-ray SED and surface
brightness radial profile depends more on the injected electron spectrum
rather than the method of transport. On the other hand, for both Model 1
($40\,$\mathrm{kyr}$$) and Model 2 ($0.002c$), the majority of electrons
remain within the large ($None$) HESS region leaving the gamma-ray SED
unchanged. However, the electrons in Model 2 have migrated further from the
pulsar while remaining within the HESS region. Subsequently, the gamma-ray
profile for Model 2 ($0.002c$) is flatter than Model 1 ($40\,$\mathrm{kyr}$$).
Figure 16 shows the distance that electrons are transported before losing
their energy to radiation. It can be seen that advection is the dominant
particle transport method for electron energies less than
$25\,$\mathrm{TeV}$$, resulting in IC emission below $7\,$\mathrm{TeV}$$.
Diffusion is dominant for electrons above $25\,$\mathrm{TeV}$$. However, these
high-energy electrons do not travel far from their birthplace before losing
their energy to radiative cooling.

At all energies, the gamma-ray energy flux along Galactic longitude for Model
1 (40 $\mathrm{kyr}$) appears symmetric around the pulsar position and shows
no preferential direction of transportation. However, the HESS uncorrelated
excess data indicates that electrons are preferentially transported to lower
Galactic longitudes. With the addition of an advective flow of $0.002c$, the
peak in the $40\,$\mathrm{kyr}$$ gamma-ray slices for photons less than
$1\,$\mathrm{TeV}$$ is now offset from the pulsar and follows the shape of the
uncorrelated excess. For the $1\,$\mathrm{TeV}$<E<10\,$\mathrm{TeV}$$ energy
band, both Model 1 and Model 2 show a flatter slice profile compared to HESS
observations. Electrons resulting from this emission appear to be contained
near their birthplace before escaping into the nebula. Our model assumed that
the diffusion lies within the Kraichnan regime with the index being fixed at
$\delta=0.5$. The top-right panel of Figure 8 shows that the modelled gamma-
ray slice morphology is broader than that observed by HESS, suggesting that
electrons are constrained within the PWN. This suggests that the diffusion
index inside the PWN may be somewhat less than the $\delta=0.5$ value we
assumed.

By assuming that diffusion was isotropic in Equation 2, any preferential
direction for particle transport was a result of magnetic field irregularities
and/or advective flow. The highly asymmetric morphology towards HESS J1825-137
could be explained if diffusion is anisotropic with preferential diffusion
towards lower Galactic longitudes. However, an anisotropic diffusion model can
be approximated by an isotropic diffusion + advection model (i.e. Model 2).

####  3.5.3 Model 3 - Isotropic Diffusion + Advection + Magnetic Field towards
HESS J1826-130

Here, a spherical shell of increased magnetic field strength around HESS
J1826-130 was considered to replicate the turbulent gas towards cloud R1 from
Voisin et al. (2016).

The bottom-right panel of Figure 8 shows that cloud R1 lies within the area
used to determine the gamma-ray SED of HESS J1825-137. The ratio of
synchrotron to inverse Compton flux is proportional to the magnetic field
(Aharonian et al., 1997). Hence, as the magnetic field around HESS J1826-130
increases, electrons lose more energy through synchrotron losses and the
inverse Compton flux decreases at equivalent gamma-ray energies. This can be
seen in the SED for HESS J1825-137 in the bottom-middle right panel of Figure
9. This has the effect of improving the match to HESS observations between
$1-10\,$\mathrm{TeV}$$ compared to Model 1 ($40\,$\mathrm{kyr}$$) as shown in
the bottom-middle-left panel of Figure 10.

The gamma-ray energy flux along Galactic longitudes are shown in the top and
top-middle panels of Figure 9. As the magnetic field around HESS J1826-130
increases, gamma-ray emission less than $1<\,$\mathrm{TeV}$$ and
$>10\,$\mathrm{TeV}$$ remains unchanged at lower longitudes, with a decrease
at higher longitudes. Between $1\,$\mathrm{TeV}$<E<10\,$\mathrm{TeV}$$, the
gamma-ray slice profile drops of at a shallower rate compared to the HESS data
at lower Galactic longitudes. However the gamma-ray emission at higher
longitudes, representing the area towards HESS J1826-130, increases
proportionally with the magnetic field.

Additionally, the bottom-left panel of Figure 9 indicates that increasing the
magnetic field strength around HESS J1826-130 lowers the contamination of HESS
J1826-130 by the PWN associated with HESS J1825-137 for energies
$<2\,$\mathrm{TeV}$$. Regions of high magnetic field strength experience a
slower rate of diffusion (see Equation 6) and high synchrotron losses. Hence
regions of high magnetic field tend to ‘block’ cosmic rays from passing
through. The model implies that a minimum magnetic field strength of
$60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ is required to successfully
lower the contamination of HESS J1826-130 according to H.E.S.S. observations.
The bottom-right panel of Figure 9 shows the multi-wavelength SED towards HESS
J1826-130. An upper-limit to the X-ray emission towards HESS J1826-130 can be
obtained using HEARSAC’s X-Ray background tool utilising ROSAT data (Sabol &
Snowden, 2019). The estimated synchrotron flux towards HESS J1826-130 combined
with the ROSAT X-ray upper limit (obtained from the same region used to
extract the SED of HESS J1826-130) implies a maximum magnetic field strength
of $\approx 20\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ around HESS
J1826-130. This constraint violation suggests that the model is not fully
encapsulating the transport of particles between PSR J1826-1334 and HESS
1826-130.

###  3.6 Future Work

Presently, our model only considers isotropic diffusion and does not account
for diffusion parallel and perpendicular to the magnetic field (Drury, 1983;
Lazarian et al., 2023). The magnetic fields of PWNe are believed to be
toroidal in nature (Schmidt et al., 1979; Kothes et al., 2006), hence
diffusion is expected to be suppressed perpendicular to the magnetic axis of
the pulsar. Additionally, the recent detection of TeV halos (Abeysekara et
al., 2017) implies that the region surrounding the PWN experiences a higher
diffusion suppression compared to the average Galactic diffusion coefficient
(Evoli et al., 2018; Di Mauro et al., 2020; Schroer et al., 2023). While
current models of particle transport suggest that advection dominates particle
transport within the PWN and diffusion dominates at the edges, this could be
described by two different regions of diffusion suppression. Model 2 and 3 in
this study considered spatially-independent advection towards lower Galactic
longitudes to explain the asymmetric $\mathrm{TeV}$ gamma-ray morphology
towards HESS J1825-137. As a result, losses due to adiabatic expansion were
neglected. Future modelling of PWNe, in particular HESS J1825-137, could
investigate the effects of inhomogeneous anisotropic diffusion and an
azimuthal and surface brightness dependent advective velocity on the gamma-ray
morphology and SED. This can then be applied to model the formation of the
$\mathrm{TeV}$ halo around the PWN (Principe et al., 2020).

Our model assumed a time-independent magnetic field with decreasing strength
from the distance to the pulsar (see subsubsection 3.2.1). However, the
average magnetic field of PWNe are expected to decrease over time from the
conservation of magnetic energy density (e.g. Tanaka & Takahara (2010)) and
the normalisation, $B_{0}$, obtained from the modelling (see Table 2) can be
considered as the time-averaged normalisation. Any future predictions of the
formation of the $\mathrm{TeV}$ halo around HESS J1825-137 must consider time-
dependency on the magnetic field.

The implementation of a time-dependent source position will not affect the
gamma-ray morphology $\gtrsim 13\,$\mathrm{TeV}$$ around PSR J11826-1334 due
to high synchrotron losses. However, lower energy photons will appear to
originate at a position offset to the current position of the pulsar as seen
in Principe et al. (2020). This does not explain the extended $\mathrm{TeV}$
gamma-ray morphology towards lower Galactic longitudes as discussed in section
3, but could affect the modelled formation of the TeV halo. Future work could
investigate the effects of an evolving source position on the gamma-ray
surface brightness radial profile and morphology towards PWNe and other
sources. For example, the application of our model to SNRs would require
cosmic rays to be injected by an expanding shell to model diffusive shock
acceleration by the expanding SNR.

The model presented in this study is not limited to HESS J1825-137 and can be
used to model the transport of cosmic-rays (electrons and protons) from other
PWNe and cosmic ray sources.

##  4 Summary

By modelling the transport of electrons across a 3D Cartesian grid of varying
ISM density and magnetic field, we are able to reproduce the main
characteristics of the multi-wavelength spectrum and morphology towards HESS
J1825-137. Three different models were considered. Model 1 assumed a simple
case of isotropic diffusion and radiative losses for the characteristic age of
$21\,$\mathrm{kyr}$$ and an older age of $40\,$\mathrm{kyr}$$ as suggested by
Van Etten & Romani (2011). Model 2 included an additional advective component
to Model 1 ($40\,$\mathrm{kyr}$$) and Model 3 introduced turbulent ISM towards
HESS J1826-130 to Model 2.

The best fit $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$ Model 1 consisted
of a pulsar injecting electrons into the surrounding medium with a spin-down
conversion factor of $10.7$ and $0.14\%$ respectively, indicating that the
true age of the system is older than the characteristic age of PSR
J11826-1334. While able to reproduce the multi-wavelength SED, neither model
was able to sufficiently reproduce the gamma-ray flux along Galactic longitude
described in H.E.S.S. Collaboration et al. (2019) for photons with energies
$1\,$\mathrm{TeV}$<E<10\,$\mathrm{TeV}$$. However, the morphological profile
could be matched for gamma-rays with energies $<1\,$\mathrm{TeV}$$ (with an
offset of $None$ towards higher Galactic longitudes compared to the HESS data)
and energies $>10\,$\mathrm{TeV}$$.

Applying an advective bulk flow (with velocity $v=0.002\,c$) of electrons
towards lower Galactic longitudes did not alter the photon SED predicted by
Model 1. By extracting the energy flux along Galactic longitude, we were able
to compare the energy-dependent morphology towards HESS J1825-137. Model 2
($0.002c$) was able to reproduce the energy flux for photons
$<1\,$\mathrm{TeV}$$ and $>10\,$\mathrm{TeV}$$, however photons with energies
$1<E<10\,$\mathrm{TeV}$$ experience a shallower drop-off compared to the
uncorrelated HESS excess slices as revealed by H.E.S.S. Collaboration et al.
(2019). This suggests that the parent electrons are constrained within the PWN
before escaping into the interstellar medium to form a $\mathrm{TeV}$ halo.
The asymmetric gamma-ray morphology towards HESS J1825-137 cannot be predicted
using a diffusion-only model and requires an advective component with bulk
flow of $0.02c$.

As the gamma-ray emission associated with PSR J1826-1334 cannot exceed the
observed emission towards HESS J1826-130, HESS J1826-130 can be used to
constrain the model. Model 1 and Model 2 were found to over-predict the SED of
HESS J1826-130 for photons $<1.5\,$\mathrm{TeV}$$. By placing a shell of
increased magnetic field strength of at least
$60\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ around HESS J1826-130,
representing the turbulent gas between the two HESS sources (Voisin et al.,
2016), the contamination was successfully lowered below the levels closer to
those estimated by HESS. By combining the modelled synchrotron flux with the
ROSAT X-ray upper limit towards HESS J1826-130, we were able to constrain the
magnetic field shell to have maximum strength of
$20\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$. This constraint violation
suggests that further modelling of the turbulent gas is needed to fully
disentangle the particle transport towards HESS J1825-137.

## Acknowledgements

T. Collins acknowledges support through the provision of Australian Government
Research Training Program Scholarship. The Nanten project is based on a mutual
agreement between Nagoya University and the Carnegie Institution of Washington
(CIW). This research has made use of the NASA’s Astrophysics Data System and
the SIMBAD data base, operated at CDS, Strasbourg, France.

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abeysekara et al. (2017) Abeysekara A. U., et al., 2017, Science, 358, 911
  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Phys. Rev. Lett., 124, 021102
  * Aharonian & Atoyan (1996) Aharonian F. A., Atoyan A. M., 1996, A&A, 309, 917
  * Aharonian et al. (1997) Aharonian F. A., Atoyan A. M., Kifune T., 1997, MNRAS, 291, 162
  * Aharonian et al. (2005) Aharonian F., et al., 2005, Science, 307, 1938
  * Araya et al. (2019) Araya M., Mitchell A. M. W., Parsons R. D., 2019, MNRAS, 485, 1001
  * Atoyan et al. (1995) Atoyan A. M., Aharonian F. A., Völk H. J., 1995, Phys. Rev. D, 52, 3265
  * Benbow (2005) Benbow W., 2005, in Aharonian F. A., Völk H. J., Horns D., eds, American Institute of Physics Conference Series Vol. 745, High Energy Gamma-Ray Astronomy. pp 611–616, doi:10.1063/1.1878471
  * Berezinskii et al. (1990) Berezinskii V. S., Bulanov S. V., Dogiel V. A., Ptuskin V. S., 1990, Astrophysics of cosmic rays. John Wiley & Sons, Inc.
  * Blumenthal & Gould (1970) Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics, 42, 237
  * Brand & Blitz (1993) Brand J., Blitz L., 1993, A&A, 275, 67
  * Brogan et al. (2006) Brogan C. L., Gelfand J. D., Gaensler B. M., Kassim N. E., Lazio T. J. W., 2006, ApJ, 639, L25
  * Cao et al. (2021) Cao Z., et al., 2021, Nature, 594, 33
  * Castor et al. (1975) Castor J., McCray R., Weaver R., 1975, ApJ, 200, L107
  * Cesarsky & Volk (1978) Cesarsky C. J., Volk H. J., 1978, A&A, 70, 367
  * Collins et al. (2021) Collins T., Rowell G., Mitchell A. M. W., Voisin F., Fukui Y., Sano H., Alsulami R., Einecke S., 2021, MNRAS,
  * Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, arXiv e-prints, pp astro–ph/0207156
  * Crutcher et al. (2010) Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ, 725, 466
  * Di Mauro et al. (2020) Di Mauro M., Manconi S., Donato F., 2020, Phys. Rev. D, 101, 103035
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Evoli et al. (2018) Evoli C., Linden T., Morlino G., 2018, Phys. Rev. D, 98, 063017
  * Gabici et al. (2007) Gabici S., Aharonian F. A., Blasi P., 2007, Ap&SS, 309, 365
  * Gabici et al. (2010) Gabici S., Casanova S., Aharonian F. A., Rowell G., 2010, in Boissier S., Heydari-Malayeri M., Samadi R., Valls-Gabaud D., eds, SF2A-2010: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. p. 313 (arXiv:1009.5291)
  * Gaensler & Slane (2006) Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17
  * Giacinti et al. (2020) Giacinti G., Mitchell A. M. W., López-Coto R., Joshi V., Parsons R. D., Hinton J. A., 2020, A&A, 636, A113
  * Giuliani et al. (2010) Giuliani A., et al., 2010, A&A, 516, L11
  * H.E.S.S. Collaboration et al. (2018a) H.E.S.S. Collaboration et al., 2018a, A&A, 612, A1
  * H.E.S.S. Collaboration et al. (2018b) H.E.S.S. Collaboration et al., 2018b, A&A, 612, A2
  * H.E.S.S. Collaboration et al. (2019) H.E.S.S. Collaboration et al., 2019, A&A, 621, A116
  * H.E.S.S. Collaboration et al. (2020) H.E.S.S. Collaboration et al., 2020, A&A, 644, A112
  * Haensel et al. (2007) Haensel P., Potekhin A. Y., Yakovlev D. G., 2007, Neutron Stars 1 : Equation of State and Structure. Springer, New York
  * Hinton & Hofmann (2009) Hinton J. A., Hofmann W., 2009, ARA&A, 47, 523
  * Isaacson (1966) Isaacson E., 1966, Analysis of numerical methods. Wiley, New York
  * Kennel & Coroniti (1984a) Kennel C. F., Coroniti F. V., 1984a, ApJ, 283, 694
  * Kennel & Coroniti (1984b) Kennel C. F., Coroniti F. V., 1984b, ApJ, 283, 710
  * Khangulyan et al. (2018) Khangulyan D., Koldoba A. V., Ustyugova G. V., Bogovalov S. V., Aharonian F., 2018, ApJ, 860, 59
  * Kothes et al. (2006) Kothes R., Reich W., Uyanıker B., 2006, ApJ, 638, 225
  * Lazarian et al. (2023) Lazarian A., Xu S., Hu Y., 2023, Frontiers in Astronomy and Space Sciences, 10, 1154760
  * Lemoine & Pelletier (2010) Lemoine M., Pelletier G., 2010, MNRAS, 402, 321
  * Li & Chen (2010) Li H., Chen Y., 2010, MNRAS, 409, L35
  * Lu et al. (2023) Lu F.-W., Zhu B.-T., Hu W., Zhang L., 2023, MNRAS, 518, 3949
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Manolakou et al. (2007) Manolakou K., Horns D., Kirk J. G., 2007, A&A, 474, 689
  * Mizuno & Fukui (2004) Mizuno A., Fukui Y., 2004, in Clemens D., Shah R., Brainerd T., eds, Astronomical Society of the Pacific Conference Series Vol. 317, Milky Way Surveys: The Structure and Evolution of our Galaxy. p. 59
  * Odegard (1986) Odegard N., 1986, AJ, 92, 1372
  * Popescu et al. (2017) Popescu C. C., Yang R., Tuffs R. J., Natale G., Rushton M., Aharonian F., 2017, MNRAS, 470, 2539
  * Porth et al. (2016) Porth O., Vorster M. J., Lyutikov M., Engelbrecht N. E., 2016, MNRAS, 460, 4135
  * Principe et al. (2020) Principe G., Mitchell A. M. W., Caroff S., Hinton J. A., Parsons R. D., Funk S., 2020, A&A, 640, A76
  * Prosekin et al. (2015) Prosekin A. Y., Kelner S. R., Aharonian F. A., 2015, Phys. Rev. D, 92, 083003
  * Protheroe et al. (2008) Protheroe R. J., Ott J., Ekers R. D., Jones D. I., Crocker R. M., 2008, MNRAS, 390, 683
  * Recchia et al. (2021) Recchia S., Di Mauro M., Aharonian F. A., Orusa L., Donato F., Gabici S., Manconi S., 2021, Phys. Rev. D, 104, 123017
  * Rees & Gunn (1974) Rees M. J., Gunn J. E., 1974, MNRAS, 167, 1
  * Sabol & Snowden (2019) Sabol E. J., Snowden S. L., 2019, sxrbg: ROSAT X-Ray Background Tool, Astrophysics Source Code Library, record ascl:1904.001 (ascl:1904.001)
  * Sano et al. (2017) Sano H., et al., 2017, ApJ, 843, 61
  * Schmidt et al. (1979) Schmidt G. D., Angel J. R. P., Beaver E. A., 1979, ApJ, 227, 106
  * Schroer et al. (2023) Schroer B., Evoli C., Blasi P., 2023, Phys. Rev. D, 107, 123020
  * Sironi et al. (2015) Sironi L., Keshet U., Lemoine M., 2015, Space Sci. Rev., 191, 519
  * Skilling (1975) Skilling J., 1975, MNRAS, 172, 557
  * Strong et al. (2004) Strong A. W., Moskalenko I. V., Reimer O., Digel S., Diehl R., 2004, A&A, 422, L47
  * Strong et al. (2007) Strong A. W., Moskalenko I. V., Ptuskin V. S., 2007, Annual Review of Nuclear and Particle Science, 57, 285
  * Stupar et al. (2008) Stupar M., Parker Q. A., Filipović M. D., 2008, MNRAS, 390, 1037
  * Tanaka & Takahara (2010) Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248
  * Tang & Chevalier (2012) Tang X., Chevalier R. A., 2012, ApJ, 752, 83
  * Taylor & Cordes (1993) Taylor J. H., Cordes J. M., 1993, ApJ, 411, 674
  * Uchiyama et al. (2009) Uchiyama H., Matsumoto H., Tsuru T. G., Koyama K., Bamba A., 2009, PASJ, 61, S189
  * Van Etten & Romani (2011) Van Etten A., Romani R. W., 2011, ApJ, 742, 62
  * Voisin et al. (2016) Voisin F., Rowell G., Burton M. G., Walsh A., Fukui Y., Aharonian F., 2016, MNRAS, 458, 2813
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * de Jager & Djannati-Ataï (2009) de Jager O. C., Djannati-Ataï A., 2009, in Becker W., ed., Astrophysics and Space Science Library Vol. 357, Astrophysics and Space Science Library. p. 451 (arXiv:0803.0116), doi:10.1007/978-3-540-76965-1_17

##  Appendix A Non-thermal emission

This section will provide an overview of leptonic interactions and the
subsequent photon emission via synchrotron, Bremsstrahlung and inverse Compton
Processes.

Synchrotron radiation occurs when an electron interacts with background
magnetic fields. The resulting photon emission from a single electron with
Lorentz factor $\gamma$ with pitch angle $\alpha$ to the magnetic field $B$ is
given by:

$\displaystyle\derivative{N}{E}$
$\displaystyle=\frac{\sqrt{3}e^{3}B}{mc^{2}}\frac{\nu}{\nu_{c}}\int_{\frac{\nu}{\nu_{c}}}^{\infty}K_{\frac{5}{3}}(x)\differential{x}\text{
,}$ (17)

where $e$ and $m$ are the charge and mass of an electron respectively,
$K_{\frac{5}{3}}$ is the modified Bessel Function, $\nu$ is the frequency of
the gamma ray and $\nu_{c}$ is the critical frequency of the emission:

$\displaystyle\nu_{c}$ $\displaystyle=\gamma^{2}\frac{3eB\sin\alpha}{4\pi
mc}\text{ .}$ (18)

The inverse Compton gamma-ray emission from an electron with energy $E_{e}$
scattering off a target photon with energy in range
($\epsilon+\differential{\epsilon}$) and number density $n\quantity(\epsilon)$
can be found using:

$\displaystyle\frac{\differential{N}}{\differential{E_{\gamma}}}$
$\displaystyle=\frac{3\sigma_{T}mc^{3}}{4\gamma}\int_{E_{\gamma}/4\gamma^{2}}^{E_{\gamma}}\frac{n\quantity(\epsilon)\differential{\epsilon}}{\epsilon}f\quantity(q,\Gamma)$
(19a) $\displaystyle f\quantity(q,\Gamma)$ $\displaystyle=2q\ln
q+\quantity(1+2q)\quantity(1-q)+\frac{1}{2}\frac{\quantity(\Gamma
q)^{2}}{1+\Gamma q}\quantity(1-q)$ (19b) $\displaystyle q$
$\displaystyle=\frac{E_{\gamma}}{\Gamma\quantity(E_{e}-E_{\gamma})}\text{,}\quad\Gamma=\frac{4\epsilon\gamma}{m_{e}c^{2}}$
(19c)

where $\sigma_{T}=(3/8\pi)r_{0}^{2}$ is the Thompson cross section, $r_{0}$ is
the classical electron radius and $F_{\text{KN}}$ takes account the full
Klein-Nishina cross section for inverse Compton scattering (Manolakou et al.,
2007):

$\begin{split}F_{\text{KN}}&=\frac{1}{u_{0}}\int_{0}^{\infty}\mathscr{f}\quantity(\gamma,\epsilon)u_{\epsilon}\differential{\epsilon}\text{,}\quad\mathscr{f}\quantity(\gamma,\epsilon)=\quantity(1+4\gamma\epsilon)^{-\frac{3}{2}}\text{
.}\end{split}$ (20)

For a Planck distribution of photon energies, $F_{\text{KN}}$ can be
approximated by:

$\displaystyle F_{\text{KN}}$
$\displaystyle=\quantity(1+4\gamma\epsilon_{\text{eff}})^{-3/2}\text{,}\quad\epsilon_{\text{eff}}=\frac{2.8kT}{m_{e}c^{2}}\text{
.}$ (21)

Finally, the photon emission from Bremsstrahlung interactions is given by:

$\displaystyle\derivative{N}{E_{\gamma}}=nc\int\differential{\sigma}\quantity(E_{e},E_{\gamma},Z)\differential{E_{e}}\text{
,}$ (22)

where $Z$ is the atomic number of the target material and
$\differential{\sigma}$ is the Bremsstrahlung differential cross section as
defined in Blumenthal & Gould (1970).

The coefficients for leptonic losses in Equation 4 are:

  * –

$b_{s}\equiv 1.292\times
10^{-15}\quantity(B/10^{3}$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$)^{2}\,${\mathrm{s}}^{-1}$$
is the synchrotron loss coefficient

  * –

$b_{c}\equiv 1.491\times
10^{-14}\quantity(n_{H}/1${\mathrm{cm}}^{-3}$)\,${\mathrm{s}}^{-1}$$ is the
Coulomb loss coefficient

  * –

$b_{b}\equiv 1.37\times
10^{-16}\quantity(n_{H}/1${\mathrm{cm}}^{-3}$)\,${\mathrm{s}}^{-1}$$ is the
Bremsstrahlung loss coefficient

  * –

$b_{\text{IC}}\equiv 5.204\times
10^{-20}\quantity(u_{0}/$\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$)\,${\mathrm{s}}^{-1}$$
is the IC loss coefficient with the energy density of photons given by $u_{0}$

  * –

$n_{H}$ is the density of the ambient hydrogen gas

The diffusion length for electrons (Atoyan et al., 1995):

$\begin{aligned}
R_{\text{diff}}&=\sqrt{\frac{4D\quantity(\gamma)}{b_{s}\gamma\quantity(1-\delta)}\quantity[1-\quantity(1-\gamma
b_{s}t)^{1-\delta}]}\end{aligned}\text{ .}$ (23)

##  Appendix B Magnetic field due to turbulent ISM gas

The magnetic field due to the ISM gas with number density $n$ is given through
Crutcher’s relation (Crutcher et al., 2010):

$\begin{aligned}
B_{\text{gas}}\quantity(n)&=\begin{cases}B_{0,\text{gas}}\text{,}&n<300\,${\mathrm{cm}}^{-3}$\\\
B_{0,\text{gas}}\quantity({n}/{300\,${\mathrm{cm}}^{-3}$})^{\alpha}\text{,}&n>300\,${\mathrm{cm}}^{-3}$\end{cases}\end{aligned}\text{
,}$ (24)

where $B_{0,\text{gas}}=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ and
$\alpha=0.65$.

##  Appendix C Single-Zone Modelling

Here we considered a ‘single-zone’ model, where electrons are injected into a
spherical region of constant number density and magnetic field (Sano et al.,
2017; Collins et al., 2021). The final electron number density is calculated
by solving Eq. C7 from Collins et al. (2021) over the age of the system where
electrons escape the region at a rate dependent on diffusion. The multi-
wavelength SED from this region is then calculated. While unable to
encapsulate the complexity towards HESS J1825-137, a general insight of the
system was gained before more detailed modelling of the morphology and time
evolution.

###  C.1 Method

| ${t=21\,$\mathrm{kyr}$}$ | ${t=40\,$\mathrm{kyr}$}$
---|---|---
Parameter | HESS | Suzaku | HESS | Suzaku
$\dot{E}$ ($\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$) | $2\times 10^{38}$ | $8\times 10^{35}$ | $1\times 10^{38}$ | $4\times 10^{35}$
$r$ ($\mathrm{\SIUnitSymbolDegree}$) | $0.70$ | $0.025$ | $0.70$ | $0.025$
$n$ (${\mathrm{cm}}^{-3}$) | $0.5$ | $0.5$ | $0.5$ | $0.5$
$B$ ($\mathrm{\SIUnitSymbolMicro}\mathrm{G}$) | $5$ | $40$ | $5$ | $40$
$\chi$ | $0.25$ | $0.25$ | $0.1$ | $0.1$
$\Gamma$ | $2.1$ | $1.9$ | $2.1$ | $1.9$
$E_{c}$ ($\mathrm{TeV}$) | $40$ | $1000$ | $50$ | $1000$

Figure 11: SED for leptonic interactions towards HESS J1825-137 using single-
zone modelling for ages $21\,$\mathrm{kyr}$$ (top) and $40\,$\mathrm{kyr}$$
(bottom). The X-ray and gamma-ray spectra are fitted separately due to the
different coverage areas of HESS and Suzaku. The orange line shows the Suzaku
observations of X-rays between $1-9\,$\mathrm{keV}$$ towards the pulsar
associated with HESS J1825-137 (Uchiyama et al., 2009). Blue data points
represent the spectrum from the _Fermi_ -LAT 4FGL source catalogue towards
HESS J1825-137 while the red data shows the H.E.S.S. energy flux towards HESS
J1825-137 (H.E.S.S. Collaboration et al., 2019). The corresponding model
parameters are shown in the table.

The X-ray and gamma-ray emission was modelled separately using two spheres
with radii $r_{\text{X-ray}}$ and $r_{\text{gamma}}$ (see Figure 3) following
the extraction regions used by Uchiyama et al. (2009) and H.E.S.S.
Collaboration et al. (2019). Electrons were injected into the spherical region
at a constant rate $\dot{E}$ and followed a power-law spectrum with an
exponential cutoff: $\derivative{N}{E}\propto E^{-\Gamma}\cdot\exp(E/E_{c})$,
where $\Gamma$ is the spectral index and $E_{c}$ is the cutoff energy. Two
different ages were modelled, $21\,$\mathrm{kyr}$$ based on the characteristic
age of the pulsar and $40\,$\mathrm{kyr}$$ based on modelling conductede by
(Van Etten & Romani, 2011).

The HESS region adopted a uniform magnetic field of
$5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{G}$$ as suggested by Principe et al.
(2020) from comparing the estimated synchrotron emission to the Suzaku X-ray
emission. Subsequently, it was assumed that the smaller X-ray region has a
higher magnetic field strength than the HESS region due to the proximity of
the pulsar and was left as a free parameter. Both the gamma-ray and X-ray
region assumed a constant ISM density of $0.5\,${\mathrm{cm}}^{-3}$$

The fits to the SED towards HESS J1825-137 as well as the modelled parameters
can be seen in Figure 11.

###  C.2 Discussion

Figure 11 shows the modelled SED with corresponding parameters to the gamma
and X-ray spectra towards HESS J1825-137 . The majority of gamma rays in this
model are predicted to come from inverse Compton interactions from the
infrared and CMB photon fields. An electron injection luminosity of $2\times
10^{38}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and $1\times
10^{38}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ is needed to match the
gamma-ray spectra at ages $21$ and $40\,$\mathrm{kyr}$$ respectively. This is
a factor ten times greater than the spin-down power of PSR J1826-1334
($\dot{E}=2.8\times 10^{36}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$). The
single-zone model assumes a time-independent injection luminosity, whereas the
spin-down power of the pulsar decreases over time. The spin-down power of PSR
J1826-1334 could have been as high as
$10^{39}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ at a pulsar age of
$1\,$\mathrm{kyr}$$ (see subsection 3.1. Therefore the modeled injection
luminosities represents the average electron injection luminosity over the age
of the pulsar.

The X-ray emission towards PSR J1826-1334 can be predicted with an injection
luminosity of $8\times 10^{35}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ and
$4\times 10^{35}\,$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$ for the $21$ and
$40\,$\mathrm{kyr}$$ model respectively. The single-zone model can reasonably
predict both the X-ray and gamma-ray SED, yet the X-ray and gamma-ray photon
models require different injection spectra for both ages of the system. The
single-zone model assumes constant density and magnetic field strength across
the region of interest. However, the magnetic field structures towards PWNe
have been suggested to be toroidal in nature but the viewing angle results in
magnetic fields appearing radially dependent or tangled (Kothes et al., 2006).
If the dense clouds towards HESS J1826-130, as seen in Figure 3, lie at the
same distance as the pulsar, diffusion will be suppressed towards this region
with electrons losing their energy to bremsstrahlung losses. As previously
mentioned, the spin-down power of the pulsar decreases over time which has an
effect on the injection luminosity of electrons in the ISM. While the single-
zone model is able to predict the X-ray and gamma-ray SED towards HESS
J1825-137, it is unable to encapsulate the complexity of the PWN.

##  Appendix D Systematic Variation of Multizone Parameters

Figure 12: Model 1 ($21\,$\mathrm{kyr}$$) as in Figure 4 but with $10\%$ (grey
shaded band) and $20\%$ (pink shaded band) variation in parameters. Figure 13:
Model 2 ($40\,$\mathrm{kyr}$$) as in Figure 5 but with $10\%$ (grey shaded
band) and $20\%$ (pink shaded band) variation in parameters.

Figure 12 and Figure 13 shows the $10\%$ and $20\%$ systematic variation of
the free parameters $\beta$, $B_{0}$, $\chi$, $\Gamma$, $E_{c}$ and $\eta$ for
the $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$ models. These figures show
that the spectral index of injection electrons, $\Gamma$, has the largest
systematic variation, where the X-ray SED and surface brightness radial
profiles show more sensitivity than the gamma-ray emission. This is a result
of the smaller region used to extract the X-ray and SED (see Figure 3 and
Figure 4).

The modelled X-ray surface brightness radial profiles for the $21$ and
$40\,$\mathrm{kyr}$$ models are steeper than observations, indicating that the
model over-predicts the synchrotron emission closer to PSR J1826-1334. This
may be corrected by decreasing the rate at which the magnetic field drops off
with distance from the pulsar ($\beta$), allowing electrons to escape the PWN
at a faster rate. The outer edges of the PWN experiences greater synchrotron
losses at the cost of $\mathrm{TeV}$ gamma-ray emission from IC interactions,
flattening out the gamma-ray surface brightness radial profile. This is
demonstrated in the $10$ and $20\%$ variation of $\beta$ shown in the top row
of Figure 12 and Figure 13. Alternatively, decreasing the overall magnetic
field strength, $B_{0}$, decreases synchrotron losses towards HESS J1825-137
at the cost of increasing the gamma-ray to X-ray flux ratio. With flux being
dependent on the observational area, any changes to the gamma-ray and X-ray
ratio will be more prominent in the X-ray SED as shown in Figure 12 and Figure
13.

To better fit the X-ray surface brightness radial profile, the diffusion
suppression coefficient, $\chi$, towards HESS J1825-137 could be increased to
allow electrons to escape further from the pulsar before losing their energy
to synchrotron radiation. High-energy electrons rapidly lose their energy to
radiative losses and remain close to the pulsar, resulting in a shallower
gamma-ray surface brightness radial profile as shown in Figure 12 and Figure
13. As the region used to extract the X-ray data is small ($<2pc$) compared to
the HESS region ($\approx 49\,$\mathrm{pc}$$), electrons quickly escape the
X-ray region while remaining in the HESS region. Thus the X-ray SED far more
sensitive to the value of $\chi$ than the gamma-ray SED.

Both the surface brightness radial profiles and SED are very sensitive to the
injected electron spectral index, $\Gamma$, as seen in Figure 12 and Figure
13. If $\beta$, $B_{0}$ or $\chi$ was altered to fit the observed Suzaku X-ray
surface brightness radial profile, the predicted SED from the model will no
longer fit to the data. In turn, the spectral index can be modified to refit
the modelled SED. Consequently, the X-ray surface brightness radial profile
will no longer match the Suzaku observations.

The cutoff energies for the $21\,$\mathrm{kyr}$$ and $40\,$\mathrm{kyr}$$
models are $45\,$\mathrm{TeV}$$ and $500\,$\mathrm{TeV}$$ respectively. As the
cutoff energy for an exponential cutoff power-law increases, the energy
spectra starts to follow a power-law. Hence, the systematic variation of
$E_{C}$ is less apparent for $40\,$\mathrm{kyr}$$ than $21\,$\mathrm{kyr}$$ as
seen in the fourth row of Figure 12 and Figure 13 respectively.

##  Appendix E Additional Figures

Figure 14: SED towards HESS J1825-137 for Model 1 ($36\,$\mathrm{kyr}$$,
green) shown vs Model 1 ($40\,$\mathrm{kyr}$$). A ‘bump’ is present in the SED
above $10\,$\mathrm{TeV}$$ for Model 1 ($40\,$\mathrm{kyr}$$) where radiative
losses are balanced by the electron injection luminosity. The
$36\,$\mathrm{kyr}$$ model has the same parameters as the $40\,$\mathrm{kyr}$$
model (see Table 2).

Figure 15: SED towards HESS J1825-137 for Model 1 ($21\,$\mathrm{kyr}$$,
$n=2$, black, see Table 2) vs Model 1 ($21\,$\mathrm{kyr}$$, $n=2.5$, green)
and Model 1 ($21\,$\mathrm{kyr}$$, $n=3$, purple).

Figure 16: The distance that electrons are transported before losing their
energy through radiative cooling (synchrotron, inverse Compton and
Bremsstrahlung) assuming purely diffusive (solid line, $\chi=0.1$) or
advective transport (dashed line, $v=0.02c$).

Figure 17: The slice profiles (top & top-middle left), surface brightness
radial profiles (top-middle right & bottom-middle left) and SED towards HESS
J1825-137 (bottom-middle right) and HESS J1826-130 (bottom) for Model 2
$0.001c$ (black), $0.002c$ (green) and $0.003c$ (purple). See Table 2 for
model parameters.