© The Authors 2024

1. Introduction

Our exploration of the sky at the highest photon energies has benefited from major advances over the past 10 − 15 yr, including: (i) the extension of the accessible energy range beyond a few tens of TeV to a PeV and above (Amenomori et al. 2019; Abeysekara et al. 2020; Cao et al. 2021), mostly thanks to facilities like the Tibet ASγ experiment, the High-Altitude Water Cherenkov Observatory (HAWC), and the Large High-Altitude Air Shower Observatory (LHAASO); (ii) the investigation of gamma-ray emission on larger angular scales (Abramowski et al. 2014; Amenomori et al. 2021; Cao et al. 2024), thanks to the rise of specific detection techniques and to progress in data analysis approaches (Knödlseder et al. 2019; Mohrmann et al. 2019; Abdalla et al. 2021; iii) the routine release of consistent high-level data sets covering significant portions of the sky, especially the Galactic Plane (HESS Collaboration 2018a; Albert et al. 2020; Cao et al. 2024). These developments provide a broader view on the cosmic-ray phenomenon and should help us in connecting the acceleration of particles in very localized astrophysical sites to their release in the vicinity of sources and the subsequent merging into a large-scale galactic population.

A number of recent observations seem to reveal that pulsars hold a predominant role in shaping the very-high-energy (VHE) and ultra-high-energy (UHE) appearance of our Galaxy. A significant number of gamma-ray sources above 50 − 100 GeV are indeed either counterparts to known pulsar-powered objects, or positionally coincident with relatively strong pulsars (HESS Collaboration 2018b; Albert et al. 2020, 2021; Cao et al. 2021). The observed emission is fully compatible with a pulsar origin in terms of spectral shape, maximum detected energy, and flux level (Sudoh et al. 2021; Breuhaus et al. 2021, 2022; de Oña Wilhelmi et al. 2022), and population synthesis efforts confirm that the majority of currently known TeV sources are powered by pulsars (Fiori et al. 2022; Martin et al. 2022a). A reliable exploitation of the growing body of gamma-ray observations therefore necessitates a solid understanding of pulsar-related emission, all the more so that deep and extensive surveys with next generation instruments like the Cherenkov Telescope Array (CTA) will be undertaken soon (Acharyya et al. 2023; The CTA Consortium 2023) and are expected to significantly increase the number of known pulsar-powered gamma-ray sources.

Pulsar-powered objects in the VHE/UHE range are primarily pulsar wind nebulae (PWNe), sources typically 3 − 30 pc in physical size involving pulsars with characteristic ages 1 − 100 kyr (HESS Collaboration 2018b). The recent discovery of TeV halos or pulsar halos (Linden et al. 2017; Abeysekara et al. 2017; Linden & Buckman 2018) seems to extend the contribution of pulsars in the VHE/UHE sky to larger physical extents (possibly > 30 pc) and older systems (> 100 kyr pulsars). It is not clear how this new kind of object relate to PWNe in terms of evolutionary path (see the reviews by Fang 2022; López-Coto et al. 2022), and how frequently the phenomenon occurs (Giacinti et al. 2020; Martin et al. 2022b). Part of the problem stems from the difficulty to characterize observationally the physical boundaries of the system and the medium that radiating particles are located in; depending on whether this is undisturbed interstellar medium (ISM), stellar ejecta in the supernova remnant (SNR), or shocked pulsar wind, the physical interpretation differs. Pulsar halos are a window on late stages of (some) pulsar-related systems and non-thermal particle transport in the vicinity of accelerators, and our inability to apprehend them reflects our limited understanding of both these topics.

As we gain access to lower surface brightness emission over larger angular scales, there is growing evidence that particle escape or leakage across the various components of the system is a key aspect in accounting for the properties of both dynamically young and evolved pulsar-related systems. This applies for instance to HESS J1825−137 (Abdalla et al. 2019; Principe et al. 2020), HESS J1809−193 (HESS Collaboration 2023), HESS 1813−178 (HESS Collaboration 2024), Vela X (Hinton et al. 2011), or Geminga and pulsar B0656+14 (Abeysekara et al. 2017).

The goal of this paper is therefore to propose a model for the gamma-ray emission from the whole system harbouring a PWN, taking into account the diffusion of particles out of the PWN into the SNR and subsequently to the surrounding ISM. Based on previous modelling efforts (Gelfand et al. 2009, hereafter GSZ09), we implemented a one-zone model for the dynamics and radiation from the PWN including a prescription for the properties of magnetic turbulence and the corresponding diffusive escape. The originality of our work is to provide a description for the subsequent fate and emission of escaping particles, namely their trapping in the remnant for some time followed by their eventual release in the vicinity of the source. The model framework is introduced in Sect. 2, we examine its typical output in Sect. 3, and apply it to the case of HESS J1809−193 and HESS J1825−137 in Sect. 4. We summarise our work in Sect. 5.

2. Modelling framework

The model framework presented here extends that introduced in GSZ09, of which we provide a detailed description including the most relevant quantities and formulae in Appendix A. A first additional feature is the possibility of diffusive escape from the nebula, computed under the assumption of particle scattering in a fully-developed Alfvenic turbulence. A second additional feature is the description of the temporary trapping of escaped particles in the remnant, followed by their eventual escape into the interstellar medium. In this section, we first provide a brief qualitative description of the main ideas behind the original framework of GSZ09, and then we introduce the assumptions underlying our model.

2.1. Main features of the original model

The model initially presented in GSZ09 describes a fast relativistic outflow emanating from the pulsar, whose size is neglected, and its interaction with the expanding stellar ejecta resulting from the parent supernova explosion. The model assumes spherical symmetry of the system and we will henceforth use r to denote the position from the centre of the system and t the time since supernova explosion and simultaneous appearance of the pulsar. We neglect any motion of the pulsar with respect to the system (but the original model includes the possibility of a non-zero velocity of the pulsar and its escape from the nebula).

The stellar ejecta interact with a cold circumstellar medium and rapidly take a characteristic structure, in which the expanding remnant drives a forward shock (FS) that propagates in the ambient medium and has position R_SNR(t), while a reverse shock (RS) propagates back in the ejecta and has position R_RS(t). The description of this double-shock dynamics is based on analytical formulae that hold only in the non-radiative stage, when energy losses from thermal radiation of the system are negligible.

Meanwhile, the central pulsar is assumed to spin down as a result of magnetic dipole radiation, which powers a relativistic outflow, the pulsar wind. At sufficiently large distance from the pulsar, the wind is mostly composed of a toroidal magnetic field and electron-positron pairs. Upon interaction with the surrounding medium, the wind is halted at a termination shock (TS), where particles are efficiently accelerated up to very high energies. The acceleration mechanism is not yet fully elucidated and may involve several sites and processes (diffusive shock acceleration at the TS, turbulent acceleration downstream of it, etc; see Bucciantini et al. 2011; Amato 2020), but observations indicate that it drains a large fraction of the wind kinetic energy, of the order of several tens of percent (Zhang et al. 2008; Bucciantini et al. 2011; Torres et al. 2014).

The PWN is the volume excavated by the shocked pulsar wind as it expands outwards, treated as a single zone in which quantities are assumed to be uniform (including energy injection). It is bounded on the inner side at a radius R_TS(t) by the TS, and on the outer side at a radius R_PWN(t) by a thin shell of swept-up ejecta material. Energy is injected in the nebula in the form of non-thermal pairs and magnetic field. The original model assumes magnetic flux conservation to determine the evolution of magnetic energy as the volume of the nebula changes, while the internal energy of non-thermal particles evolves as a result of adiabatic and radiation losses.

The bounding shell mass is made of ejecta material collected when the shell velocity exceeds the ejecta velocity immediately ahead of it. Its dynamics is controlled by its inertia and the pressure difference between the nebula on the inner side, and the stellar ejecta on the outer side. At early times, the innermost unshocked stellar ejecta are assumed to have a negligible pressure because of their fast expansion and rapid cooling. Later on, the reverse shock propagating back through the ejecta eventually reaches the nebula, heats and compresses the material ahead of the bounding shell such that the pressure immediately outside the nebula becomes non-negligible and possibly well in excess of that in the nebula. This may lead to a strong compression of the nebula, initiating several cycles of compression and re-expansion referred to as reverberation.

The spectral distribution of particles at any given time depends on the assumed injection properties (pulsar spin-down history and initial particle spectrum), the magnetic and radiation fields in which they are immersed (which sets the magnitude of radiative losses), and the dynamics of the nebula (which drives adiabatic energy losses or gains). It is used to compute the broadband radiation from the nebula, consisting of synchrotron emission in the nebular magnetic field and inverse-Compton scattering in the ambient photon fields.

2.2. Additional developments

2.2.1. Particle escape from the nebula

We modified the energy injection equation to include a term for turbulence, not part of the original model of GSZ09:

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{inj}}(t)& ={\dot{E}}_{\mathrm{inj},B}(t)+{\dot{E}}_{\mathrm{inj},\mathrm{e}}(t)+{\dot{E}}_{\mathrm{inj},\mathrm{T}}(t)={\dot{E}}_{\mathrm{PSR}}(t),\hfill \end{array}$$(1)

$$\begin{array}{cc}& ={\eta }_B(t){\dot{E}}_{\mathrm{PSR}}(t)+{\eta }_{\mathrm{e}}(t){\dot{E}}_{\mathrm{PSR}}(t)+{\eta }_{\mathrm{T}}(t){\dot{E}}_{\mathrm{PSR}}(t).\hfill \end{array}$$(2)

Term Ė_inj,e (and corresponding η_e) describes the injection of relativistic pairs in the nebula. Term Ė_inj,B (and corresponding η_B) describes a large-scale regular magnetic field component, similar to the original model, while term Ė_inj,T (and corresponding η_T) denotes a turbulent magnetic field component.

Turbulence is assumed to be composed of Alfvén waves so the injected energy consists for one half of kinetic energy of the turbulent fluid, and for the other half of its magnetic energy. For simplicity and consistency, the evolution with time and volume of turbulent energy and pressure is taken to be similar to what was assumed in GSZ09 for the large-scale magnetic field term.

The turbulence spectrum is ideally obtained by solving a complete wave transport equation for a given phenomenology (e.g. Kolmogorov or Kraichnan). In our context, turbulence is expected to be injected at a spatial scale which is a fraction of the TS radius. Turbulence cascading as a result of non-linear wave-wave interactions, which can be described as diffusion in wavenumber space (Miller & Roberts 1995), transports turbulent energy to larger wavenumbers (smaller physical scales), while volume expansion pushes turbulence towards smaller wavenumbers (larger physical scales).

For typical parameters of the problem, the Alfvén speed is quite high, up to a few tens of percent of the speed of light, such that the development of the cascade occurs on time scales that are much shorter than the dynamical time scale. We therefore considered a simplified description of turbulence, avoiding the solving of a complete transport equation, and assumed that the magnetic turbulence spectrum is fully developed at each time, from a maximum scale that is a fraction κ_T ≤ 1 of the radius of the nebula down to an arbitrarily small scale. In practice:

$$\begin{array}{cc}& w(k,t)=\frac{1}{3}\frac{E_{\mathrm{PWN},\mathrm{T}}(t)}{V_{\mathrm{PWN}}(t)}\frac{1}{(k_{\mathrm{inj}}^{-2/3}-k_{\mathrm{cut}}^{-2/3})}{k}^{-5/3},\hfill \end{array}$$(3)

$$\begin{array}{cc}& k_{\mathrm{inj}}(t)=\frac{2\pi }{{\kappa }_{\mathrm{T}}{R}_{\mathrm{PWN}}(t)}\mathrm{and}{k}_{\mathrm{cut}}(t)\gg k_{\mathrm{inj}}(t).\hfill \end{array}$$(4)

A justification for the assumption of fully-developed turbulence at all times is provided in Appendix B. We do not describe the energization of particles by magnetic reconnection in the turbulent flow, or stochastic scattering off turbulent fluctuations (Luo et al. 2020; Lu et al. 2023). These processes were considered in the context of the origin of the radio-emitting particle population, and in relation to the question of the confinement of the nebula. An implicit hypothesis in our model framework is that the conversion of ordered toroidal field into magnetic turbulence and the dissipation of magnetic energy occur on sufficiently small temporal and spatial scales downstream of the TS, and that the assumed turbulent state and particle injection spectrum in our model may be the result of these processes.

Electron-positron pairs are assumed to experience advection in the nebula, from the termination shock to the bounding shell. Advection is assumed to occur following the prescription used in Zhu et al. (2023), which involves a 1/r velocity profile down to a minimal value that is the velocity of the thin bounding shell, while diffusion results from resonant interactions with magnetic turbulence on the scale of the Larmor radius of the particle, which is described by the coefficient:

$$\begin{array}{cc}& D_{\mathrm{PWN}}(p,t)=\frac{1}{3}\frac{c{r}_{\mathrm{L}}(p,t)}{k_{\mathrm{res}}W(k_{\mathrm{res}},t)},\hfill \end{array}$$(5)

$$\begin{array}{cc}& W(k,t)=\frac{w(k,t)}{\frac{E_{\mathrm{PWN},B}(t)+E_{\mathrm{PWN},\mathrm{T}}(t)/2}{V_{\mathrm{PWN}}(t)}},\hfill \end{array}$$(6)

where p is the momentum of the particle. We used the approximation introduced in Ptuskin & Zirakashvili (2003, their Eq. 5), valid for an arbitrary turbulence level, in which the Larmor radius r_L is taken in the total ordered+turbulent magnetic field, k_res = 2π/r_L is the corresponding wavenumber at the resonant scale, and W(k) is the turbulence spectrum normalized to the ordered+turbulent magnetic field energy density.

The spectral distribution of the escaping flux at each time, Q_esc, PWN(E, t), is computed from the spectral distribution of particles in the nebula N_e, PWN(E, t) and the escape time scale τ_esc, PWN(E, t):

$$\begin{array}{cc}& Q_{\mathrm{esc},\mathrm{PWN}}(E,t)=-\frac{N_{\mathrm{e},\mathrm{PWN}}(E,t)}{{\tau }_{\mathrm{esc},\mathrm{PWN}}(E,t)},\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\tau }_{\mathrm{esc},\mathrm{PWN}}(E,t)={(\frac{1}{{\tau }_{\mathrm{diff},\mathrm{PWN}}}+\frac{1}{{\tau }_{\mathrm{adv},\mathrm{PWN}}})}^{-1},\hfill \end{array}$$(8)

where the escape time scale combines the diffusion and advection crossing times τ_diff, PWN and τ_adv, PWN:

$$\begin{array}{cc}& {\tau }_{\mathrm{diff},\mathrm{PWN}}(E,t)=\frac{R_{\mathrm{PWN}}^2(t)}{6D_{\mathrm{PWN}}(E,t)},\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\tau }_{\mathrm{adv},\mathrm{PWN}}(t)=\frac{(R_{\mathrm{PWN}}^2(t)-R_{\mathrm{TS}}^2(t))}{2R_{\mathrm{PWN}}(t)v_{\mathrm{PWN}}(t)},\hfill \end{array}$$(10)

where the diffusion coefficient is expressed as a function of kinetic energy E rather than momentum p as in Eq. (5). We assume here three-dimensional diffusion, although we did not make any explicit assumption on the topology of the ordered field. The escape flux term is included on the right-hand side of Eq. (A.14).

This implementation essentially describes the nebula as a magnetized turbulent bubble that holds particles only until they reach its outskirts, and implies that the surrounding bounding shell plays no role in confining them. This assumption can be justified from radio polarimetric observations of the Crab showing that the magnetic field in the nebula evolves from mostly toroidal close to the termination shock to predominantly radial at larger distances (Bietenholz & Kronberg 1990). Particles can then slide outwards along the field lines, which could explain the limited spectral index variations measured in X-rays along filaments in the outer regions (Seward et al. 2006).

2.2.2. Particle trapping in the remnant

Particles escaping from the nebula enter the remnant. The transport of pairs across the remnant will depend on the actual magnetic conditions in the volume and may well be pretty complicated. There are frontiers that particles may not easily cross before reaching the outer FS, first and foremost the PWN bounding shell but also the ejecta-ISM contact discontinuity, both of which are expected to be turbulent because of hydrodynamical instabilities. This complexity is beyond the scope of our model and we adopted a minimalist description of particle confinement in the SNR volume, focussing on what happens at its outer edge.

Because of particle acceleration and magnetic field amplification at the FS, that frontier acts as a magnetic barrier confining energetic particles up to a certain maximum momentum downstream of the shock. That maximum momentum evolves with time and we adopted the prescription used in Celli et al. (2019):

$$\begin{array}{c}\hfill p_{\max }(t)=\{\begin{array}{cc}{p}_{\max }^{\mathrm{ST}}(t/t_{\mathrm{ST}})\hfill & t\le t_{\mathrm{ST}},\hfill \\ p_{\max }^{\mathrm{ST}}{(t/t_{\mathrm{ST}})}^{-\xi }\hfill & t>t_{\mathrm{ST}}.\hfill \end{array}\end{array}$$(11)

Equation (11) implicitly encapsulates the intricate physics governing the evolution of the amplified magnetic field at the FS. We handle the absolute maximum momentum $p_{\max }^{\mathrm{ST}}$ and the time dependence index ξ as free parameters of the model, with typical values in the range ∼0.1 − 1 PeV and ∼2.0 − 4.0, respectively. As an example, the interpretation of gamma-ray observations of the γ-Cygni SNR in the framework of an escape model based on that prescription resulted in estimated values of $p_{\max }^{\mathrm{ST}}=78$ TeV/c and ξ = 2.55 (MAGIC Collaboration 2023).

The SNR is treated as a one-zone region hosting a specific population of energetic particles, fed by a continuous inflow of particles escaping from the PWN and progressively depleted of particles with energies above p_max(t). As the SNR evolves, the confined particles lose energy as a result of adiabatic expansion, inverse-Compton scattering, and synchrotron radiation. Particles trapped in the SNR, until their complete release into the ISM when the remnant is deep into the Sedov-Taylor stage, power an extended emission component around the PWN, with a significant flux comparable to or well in excess of that of the nebula in some cases, as we will see below.

The magnetic field in the remnant can have several origins: swept-up interstellar magnetic field, possibly including some amplification as a result of particle acceleration at the FS, magnetic energy frozen in the stellar ejecta, magnetic field amplification by particles streaming out of the nebula, ... These different contributions are expected to have specific spatial distributions, and in some cases to be restricted to compressed layers or filaments. For simplicity, and because of the one-zone assumption for this volume, we considered a magnetic field in the remnant similar to that in the surrounding ISM. This is obviously a lower limit case, mostly relevant to evolved systems with ages ≳2 − 3 kyr and size ≳10 − 20 pc.

3. Model predictions

3.1. Model applicability

We present below predictions for a reference pulsar-PWN-SNR system. We caution beforehand that the application of the model framework defined in Sect. 2 should be restricted to early evolutionary stages, which we define as the time period from pulsar birth through free expansion of the nebula and up to beyond reverse-shock interaction until modest levels of compression of the PWN are reached (less than a factor of ten). Indeed, it was demonstrated that analytical prescriptions for the hydrodynamical properties of the remnant such as those used in our model are inappropriate to accurately describe the evolution of the system past reverse-shock crushing (Bandiera et al. 2023a,b). Second, reverse-shock crushing can reasonably be expected to be asymmetric, thereby reflecting gradients in the surrounding ISM density, which combined with a non-zero pulsar kick not necessarily aligned with the density gradient can give rise to a myriad of situations that our spherical model can definitely not describe properly (see examples in Kolb et al. 2017). In its current version, the model cannot be used to describe very evolved systems, for instance those with two detached nebulae (a relic nebula and a younger bow-shock nebula in the SNR or ISM) or TeV halos.

We also emphasise that we do not describe the spatial transport of particles that manage to escape into the ISM. As it is, the model cannot self-consistently predict the transport conditions in that medium, so we only compute the spectral evolution of the escaped particles in time and the corresponding integrated non-thermal emission. Our purpose is to show that the ISM contribution can be significant in terms of total signal (see Sects. 3.4 and 3.5 below), such that observations should tell us something about diffusion/transport conditions close to the source.

3.2. Model setups

In our reference system, the SNR component is parametrized as an ejecta of mass M_ej = 10 M_⊙ and initial kinetic energy E_ej = 10⁵¹ erg, expanding into a uniform medium with hydrogen density n_ISM = 0.1 H cm⁻³ and mean molecular mass μ = 1.4. The initial ejecta density profile consists of a uniform core and a power-law envelope with index n_ej = 9, and the initial velocity profile is a linear growth starting from zero at the centre. Particle confinement within the remnant is defined by default by $p_{\max }^{\mathrm{ST}}={10}^5$ GeV/c and ξ = 2.5, but we explored variations of these parameters.

The pulsar component has an initial spin-down power L₀ = 5 × 10³⁸ erg s⁻¹ and spin-down time scale τ₀ = 2000 yr. The pulsar spins down over time with a constant braking index n_PSR = 3. For simplicity, we set the natal kick velocity of the pulsar to zero. At the wind termination shock, a fraction η_e = 90% of the wind power is converted into non-thermal electron-positron pairs, and we consider various possibilities for the splitting of the remaining 10% into ordered magnetic field and Alfvenic turbulence (via parameters η_B and η_T). Downstream of the termination shock, non-thermal particles are injected in the nebula with an energy spectrum consisting of a broken power-law with indices α₁ = 1.5 and α₂ = 2.3 respectively below and above a break energy E_brk = 500 GeV, while an exponential cutoff terminates the spectrum at E_cut = 1 PeV.

The whole system is located at a distance of 3 kpc from us, in the direction of the Galactic centre. The interstellar radiation field density at this position is taken from the model of Popescu et al. (2017), and the interstellar magnetic field is assumed to have a strength of B_ISM = 5 μG.

We explore the effects of different prescriptions for the turbulence, hence different particle escape profiles, on the resulting non-thermal emission. We considered two families of setups: low level of magnetic turbulence (η_B = 0.09 and η_T = 0.01), high level of magnetic turbulence (η_B = 0.01 and η_T = 0.09), and for each case we tested three spatial scales for turbulence injection (κ_T = 0.01, 0.1, 1.0). We emphasise that, because turbulent energy is half kinetic and half magnetic in the adopted Alfvenic turbulence, the different families of setups correspond to different levels of magnetic energy injection: 9.5% for the L setups, and 5.5% for the H setups. In addition to the impact on particle escape, there will thus be some impact on synchrotron radiation losses, hence on particle spectra at the high energy end.

The different η_B, η_T, and κ_T parameter combinations have dynamical implications, in the sense that they have consequences on the amount and form of the energy contained in the PWN, hence on its expansion rate or ability to withstand compression. Conversely, parameters $p_{\max }^{\mathrm{ST}}$ and ξ controlling particle escape from the SNR into the ISM have no effect on the dynamics of the system, and they only affect the radiation from the ISM relative to that from the SNR.

3.3. Energetics and dynamics

Figure 1 displays the dynamics of the forward and reverse shocks in the remnant, and of the thin bounding shell forming the outer frontier of the nebula. The reverse shock hits the nebula at t ≃ 7000 yr, which triggers its compression only about one thousand years later owing to the inertia of the bounding shell.

Fig. 1. Evolution of the size of the PWN for the two extreme model setups. Also shown as reference are the radii for the forward and reverse shocks of the SNR. The curves for the other setups lie in between the red and purple curves. In this and subsequent figures, the gray shaded area denotes the strong compression phase that the model is not suited to describe.

Figure 2 displays the time evolution of particle escape compared to that of radiative losses and pulsar spin-down. One can first note the similarity of the curves in the two cases with highest turbulence (red and orange), which reflects the fact that advection is then the dominant escape process. Our model suggests a maximum in escape losses somewhat around the decrease of the pulsar spin-down power at t ≃ τ₀. In high-turbulence cases, there is a second peak at the time of compression, following reverse-shock interaction. The actual dependence of the escape term as implemented in the model can be translated into an overall dependence as $B_{\mathrm{PWN}}^{-1/3}{R}_{\mathrm{PWN}}^{-4/3}$. Because of this scaling and the actual evolution of B_PWN and R_PWN, conditions are less and less favourable to particle escape as time goes by, at least until compression. The first peak in the escape luminosity results from a trade-off between the increasing escape time scale (mostly because the nebula grows) and the rising lepton energy content in the nebula (because of continued particle injection at a high rate while radiative losses diminish). A second peak can occur in cases with sufficient turbulence and results from the strong decrease of the PWN size and the simultaneous energization of the particle population upon compression. In all cases, particle escape is so efficient that the nebula is emptied of most of its leptons upon compression. Interestingly, this prevents the appearance of a superefficiency period, during which radiative losses exceed the instantaneous pulsar spin-down. The difference in particle escape luminosities between our low and high turbulence scenarios reaches one to two orders of magnitude. As illustrated in Fig. 3, this leads to very different lepton energy content when the nebula gets close to its maximum extent before compression, which will translate into distinctive radiation signatures as we will see below.

Fig. 2. Evolution of the radiation and escape luminosities in the PWN volume for the different model setups. Also shown as reference is the pulsar spin-down power.

Fig. 3. Evolution of the energy content in the PWN volume for the different model setups. Also shown as reference is the cumulative energy injected by the pulsar and the cumulative energy lost to radiation in the nebula.

Conversely, the impact of the turbulence prescription on the extent of the PWN is very modest as shown in Fig. 1. Low turbulence setups, and the higher escape they allow, result in a slower growth of the nebula and a later interaction with the reverse shock. The main reason for such a moderate impact can be understood from Figs. 2 and 3. In the low turbulence case, particle escape dominates over radiation losses after a few decades to a few centuries, and the particle energy content of the nebula remains below the ordered+turbulent magnetic energy at all times. The latter component therefore determines the pressure in the nebula that drives the expansion of the bounding shell. In the high turbulence case, particle escape becomes a significant loss term only after a few centuries to a few millenia. As a consequence, the particle energy content is much higher, by a factor of a few up to a few tens. Yet, the corresponding pressure only rivals that of the ordered+turbulent magnetic component at late times, when the pulsar has released most of its rotational power, and it never exceeds it by much such that the impact on the PWN dynamics remains limited.

We can discuss the impact of neglecting any kick velocity of the pulsar. In the reference model runs, the PWN has reached a radius of ∼10 pc in ∼8000 yr at the onset of compression, with a radius growing faster than linearly in time over that period (R_PWN ∝ t^1.2; see van der Swaluw et al. 2001). Consequently, all pulsars with a kick velocity below about 1200 km s⁻¹ would still be contained in the nebula by then, which accounts for the vast majority of the population (Verbunt et al. 2017). So the kick velocity would not have an influence. Model setups exist for which the pulsar leaves the nebula during compression, implying that particles would be injected almost directly in the remnant after some time, thereby enhancing even more the relative contribution of the SNR compared to the PWN (see Sect. 3.4).

Figure 4 displays the time evolution of the spatial diffusion coefficient for the different model setups, at a particle energy of 100 TeV. For reference, the effective value of the diffusion coefficient at 100 TeV for large-scale transport in the Galactic disk is on the order of 10³⁰ cm² s⁻¹. Turbulence in the nebula as prescribed here results in a diffusion coefficient that is initially four to eight orders of magnitude smaller than in the ISM. With the decrease of the magnetic field strength and increase of the largest scale of the turbulence, at least up to reverse-shock interaction, the diffusion coefficient progressively rises to about the interstellar value at best, for the model setup with the strongest escape, or to a value that is three orders of magnitude smaller, for the model setup confining particles the most efficiently.

Fig. 4. Evolution of the spatial diffusion coefficient at 100 TeV in the different model setups. For reference, the effective value of the diffusion coefficient at 100 TeV for large-scale transport in the Galactic disk is on the order of 1030 cm2 s−1.

3.4. Gamma-ray luminosity

Figure 5 shows the time evolution of the GeV and TeV luminosities of all emission components of the system for four different model setups. In the low turbulence cases, our model predicts that GeV–TeV radiation from the SNR significantly exceeds that from the PWN, by at least about one order of magnitude and up to nearly three orders of magnitude in scenarios with significant escape. For old enough systems, integrated emission from the surrounding ISM is at least comparable to or much higher than that of the SNR in the TeV range, while it is largely subdominant in the GeV range past a few decades. The latter statement is however dependent on the parameters $p_{\max }^{\mathrm{ST}}$ and ξ controlling the time evolution of particle confinement within the remnant. In the high turbulence cases, GeV–TeV emission from the SNR is unsurprisingly weaker than in the low turbulence case but still predicted to be a significant contribution to the overall signal, with a total flux within a factor of a few of the emission from the PWN over most of the valid time range. TeV radiation from the ISM is also much reduced and becomes comparable to that of the SNR past a few centuries, while the corresponding GeV radiation remains subdominant. The behaviour observed during compression of the nebula also differs between the low and high turbulence cases: in the low turbulence cases, compression results in an abrupt drop of the PWN emission, both in the GeV and TeV ranges, while the high turbulence cases exhibit a modest brightening of the nebula.

Fig. 5. Evolution of the gamma-ray luminosity in different bands for the PWN, SNR, and ISM components in four model setups.

The bulk of the detectable population can be expected to be made up of objects maximizing the age and luminosity product (because the former increases the size of the population, while the latter increases the distance up to which it can be detected). For the reference setups studied here, this would be objects with ages ∼5 − 10 kyr, for which our model predicts a non-negligible if not dominant contribution from the SNR. In the TeV range, an additional contribution from the surrounding medium can also exist, with a total flux possibly exceeding that of the other components (depending on the parameters of the system), but likely spread over a much larger volume such that its detection may be challenging. Emission from the PWN is the dominant contribution only for younger and less evolved objects, with ages ranging from a few centuries to a few millenia, under the assumption of high turbulence (see the red and green curves in Fig. 5). Yet, even in this case, the contribution from the SNR is at a comparable level at both GeV and TeV energies. This seems consistent with the results obtained for source MSH 15−52 that, at an age of order of 1000 yr, already shows signs of particle escape from the nebula to the remnant (Tsirou et al. 2017). It may also be consistent with the latest GeV and TeV analysis of the Crab nebula, which reveals a small extension of the source, hardly above the angular resolution limit of the instruments (Aharonian et al. 2024): the gamma-ray emission is more extended than the X-ray synchrotron emission, and possibly more than the radio and optical emission especially at GeV energies, which could be indicative of escape out of the nebula.

The bottom panel of Fig. 5 shows the predicted luminosities in the 25 − 100 TeV range. In most model setups considered here, emission in that band is dominated by the ISM contribution, which may spread over spatial scales much larger than the SNR+PWN system size. Observatories like Tibet ASγ, HAWC, and LHAASO may therefore be appropriate to search for this emission component and thereby probe particle escape and transport in the vicinity of the source. Interestingly, a large fraction of the sources in the first LHAASO catalogue seem to be positionally coincident with pulsars and very extended, with a 39% containment radius above 0.5° and up to 2°, even for the KM2A subsystem (Cao et al. 2024). The exact ISM flux level relative to the other components however depends on the details of particle confinement within the remnant. The ISM luminosity has a characteristic shape resulting from the assumed time evolution of the maximum energy of particles contained in the SNR. First, the flux in a given band increases until p_max(t) rises to a value high enough to confine the emitting particles within the remnant. The supply of fresh particles to the ISM is then suspended and the emission exhibits a plateau. At a later time, past the Sedov-Taylor time, the forward shock progressively weakens and p_max(t) eventually drops down below the value for particles confinement within the SNR and the ISM is fed again with new particles, which produces a rise in the emission while at the same time the emission from the remnant decreases as it gets emptied. The duration of the plateau is set by parameters $p_{\max }^{\mathrm{ST}}$ and ξ and increases with decreasing energy, as can be seen in Fig. 5, where the GeV plateau is much longer than the TeV one.

3.5. Gamma-ray spectra

Figure 6 shows the gamma-ray spectra of each emission component for four different model setups and at three times: close to the maximum lepton content of the nebula (1 kyr), shortly before reverse-shock interaction and after the first peak of particle escape (5 kyr), and at the beginning of compression (10 kyr). Importantly, these spectra are integrated over quite different spatial regions. From 1 to 10 kyr, the PWN component arises from a volume increasing from 1 to 10 pc, while the SNR component is produced in a sphere growing from 5 to 20 pc. The spatial extension of the ISM component is unclear and will depend on the specifics of particle transport in the vicinity of the remnant: it could be nearly as small as the SNR in case of strong confinement, or several hundreds of pc across in case of diffusion with the average large-scale properties for the Galactic disk (typically ∼10³⁰ cm² s⁻¹ at 10 TeV).

Fig. 6. Gamma-ray spectra of the PWN, SNR, and ISM components in four model setups at different ages. The gray dotted line corresponds to an order-of-magnitude estimate of the pion-decay emission from cosmic rays in the remnant (see text for details).

In all panels, we provide for comparison a spectrum for pion-decay emission from cosmic rays in the remnant. The latter is computed assuming that, at all times, a constant fraction 5 × 10⁻⁷ of the material swept-up by the forward shock has undergone diffusive shock acceleration, which produced a power-law distribution in momentum having a constant index 2.3 and extending up to the maximum momentum p_max(t) defined by Eq. (11). The injection fraction of 5 × 10⁻⁷ results in about 8% of the ejecta kinetic energy being converted into accelerated particles at an age of 10 kyr. These particles interact with gas in the SNR, swept-up interstellar material and ejecta, assumed to be uniformly distributed within the volume for simplicity. The resulting gamma-ray emission is scaled up by a nuclear enhancement factor of 1.845 (Mori 2009) to account for nuclei in the cosmic-ray population and ejecta. This calculation provides an order-of-magnitude estimate of the emission from the SNR, for comparison to that powered by particles escaping the PWN. It neglects the loss of particles escaping upstream of the forward shock, higher-than-solar metallicity in the ejecta or ISM, and any radial structure of thermal and non-thermal gas in the SNR.

The low turbulence spectra in Fig. 6 show the dominance of the SNR and ISM components over the emission from the PWN over 1 − 10 kyr. The SNR component is characterized by hard emission up to 0.1 − 1 TeV, at which energy even harder emission from the ISM takes over. This crossover energy depends on the parameters controlling confinement in the remnant, and it decreases with time as more and more particles escape from the system. In the high turbulence spectra of Fig. 6, the PWN component weighs more in the total signal, and is probably the dominant one in terms of surface brightness as it arises from the smallest volume. With time, the contribution from the SNR (at ≲1 TeV) and the ISM (at ≳1 TeV) gains in intensity such that at ages ≳5 − 10 kyr, it constitutes a non-negligible, if not utterly dominant, fraction of the total integrated signal.

In all model setups, it seems that inverse-Compton emission powered by particles escaping the PWN dominates pion-decay emission from cosmic rays accelerated at the forward shock and advected downstream in the SNR. The pion-decay spectrum provided for comparison is admittedly a simple estimate, and enhancements by a factor of a few are conceivable assuming a higher acceleration efficiency or hadronic interactions in a denser gas layer within the remnant.

Neglecting the ISM component, assuming for instance rapid particle diffusion after decoupling from the SNR, the predicted spectra are characterized by a very distinctive pattern that may allow one to indirectly infer the turbulence level: in the low turbulence case, the emission from the system drops off rather abruptly beyond ∼1 − 3 TeV, and the PWN takes over beyond ∼10 − 30 TeV with an intensity that is much lower, about two orders of magnitude or more; conversely, in the high turbulence case, the total PWN+SNR emission is more smooth over the full spectral range. We will see in Sects. 4.2 and 4.3 that the latter case seems to be more representative of real systems.

In Fig. 7, we illustrate the effect on the spectra at 10 kyr of the maximum particle momentum $p_{\max }^{\mathrm{ST}}$ allowed in the remnant, for three values 10⁴, 10⁵, and 10⁶ GeV/c. With increasing maximum momentum, the cutoff in the SNR emission is pushed to higher and higher energies, while the ISM emission is restricted to a smaller and smaller energy range at the high end of the spectrum. One should note that the sum of the two components is not constant because allowing particles to escape more rapidly and easily from the remnant preserves them from adiabatic losses.

Fig. 7. Gamma-ray spectra of the PWN, SNR, and ISM components at 10 kyr in a high-turbulence case, as a function of the maximum particle momentum $p_{\max }^{\mathrm{ST}}$ allowed in the remnant. The PWN emission is the same in all three runs since it is not affected by this parameter. The gray dotted line corresponds to an order-of-magnitude estimate of the pion-decay emission from cosmic rays in the remnant (see text for details).

Although the main goal of this paper is to investigate the effect of escape in PWNe on their gamma-ray emission, the model self-consistently predicts synchrotron emission. The corresponding spectra for both low and high turbulence model setups are presented in Appendix C, with emphasis on radio emission where our model is more reliable.

3.6. Gamma-ray morphology

The results discussed above showed that emission from the SNR may be a non-negligible, if not dominant, contribution to the emission from the system (leaving aside emission from the surrounding ISM). The exact share depends on the turbulence parameters, and on the energy at which the system is observed. As such, the gamma-ray emission from the whole system may have a specific energy-dependent morphology indirectly conveying information on turbulence. By construction, because of the one-zone assumption for each component of the system, the morphological information that can be extracted from the model is limited. Nevertheless, it provides interesting trends that are reminiscent of some observed patterns.

We computed the intensity distribution of the PWN+SNR emission as a function of the angular distance from the centre. Assuming uniform emission properties within each spherical volume, PWN or SNR, we integrated the volumetric emissivity along each line of sight from the inside out (or from the central pulsar to the forward shock), and then assessed the 68% containment radius of the emission at each gamma-ray energy. The results are displayed in Fig. 8, for the six model setups considered and at a system age of 10 kyr.

Fig. 8. The 68% containment radius of the PWN+SNR emission as a function of energy, for the different model setups and at a system age of 10 kyr. The gray and black dotted lines show the angular extent of the PWN and SNR, respectively, for reference.

In the low turbulence cases, the SNR contribution is so strong that the morphology is dominated by the remnant, with a 68% containment radius at about two thirds of the forward shock extent. This holds until a few TeV, the maximum energy radiated by the highest-energy leptons still contained in the SNR at this age. Beyond that limit, the SNR contribution collapses, as illustrated in the bottom panel of Fig. 6, and the morphology becomes dominated by the PWN, which leads to a drop by a factor ∼2 of the typical size of the emission. In the high turbulence cases, this trend is mitigated by the fact that the SNR contribution is more comparable to that of the PWN.

Adding the possibility of a contribution from the surrounding ISM, especially in the case of some confinement in the vicinity the system, would lead to an even richer set of likely morphologies. This is particularly relevant for instruments like HAWC or LHAASO, with their capabilities to detect and image very extended sources above 10 TeV.

4. Application to known sources

We applied our model on sources HESS J1809−193 and HESS J1825−137, that we selected for the following reasons: (i) there is in each case an identified pulsar with measured properties (spin period and derivative); (ii) they have extended gamma-ray coverage from GeV to nearly 100 TeV; (iii) they are intermediate in terms of dynamical evolution, with ages of order 10−20 kyr, which falls into the applicability domain of our model framework; (iv) they were discussed in past literature as sources for which particle escape is an important aspect, notably because of the large physical size of their gamma-ray emission.

About point (iii), there are indications from observations that both systems have undergone reverse-shock interactions (strongly asymmetric development of the nebula with respect to the pulsar; see the review of observations in Sects. 4.2.1 and 4.3.1). Yet, we caution that it is hard from observations alone to quantify for how long the systems have been in such a stage, and therefore to guarantee that they perfectly fit within the domain of applicability of our model.

About point (iv), this is what justifies the development and use of the three-zone model presented here. One-zone models in their original form, like that introduced in GSZ09, have difficulties reproducing the physical extent of such sources without stretching the remnant or pulsar parameters to extreme values (e.g. the discussion on HESS J1825−137 in de Jager & Djannati-Ataï 2009). Conversely, our model, which allows for the propagation of emitting particles across the whole system out to the surrounding medium, can naturally produce very extended sources (at least qualitatively, since the current version of the model does not compute the radial distribution of emitting particles, especially in the interstellar medium, as emphasised in Sect. 3.1).

4.1. Fitting strategy

In our application of the model to HESS J1809−193 and HESS J1825−137, we will focus on gamma-ray observations only at this stage. The purpose is to obtain a quantitative description of their spectra, and a semi-quantitative description of their morphological properties (because this three-zone model does not yet include a detailed computation of the spatial transport).

The model has a large number of free parameters in total, and some degree of non-linearity, which altogether offers some freedom in fitting data and some dependency on initial values. In order to get meaningful results and a better comprehension of the physical effects, we restricted the fitting procedure to the main subset of the following parameters:

1. Pulsar spin-down time scale τ₀.

2. Ejecta mass M_ej.

3. Ejecta energy E_ej.

4. Interstellar medium density n₀.

5. Turbulence maximum scale κ_T.

6. Turbulent energy injection efficiency η_T.

7. Particle injection spectrum low-energy index α₁.

8. Particle injection spectrum high-energy index α₂.

9. Particle injection spectrum break energy E_b.

10. Particle escape maximum momentum $p_{\max }^{\mathrm{ST}}$.

We emphasise that there is some degeneracy between parameters κ_T and η_T, in that the particle escape flux term scales roughly as ${\kappa }_{\mathrm{T}}^a/{\eta }_{\mathrm{T}}^b$, where the exponents a and b depend on the turbulence regime. But because there is no strict anti-correlation over the entire parameter space, we kept both quantities in the fitted set.

We aimed at solutions involving pulsar and supernova remnant properties in line with the statistical distributions inferred from population studies. For pulsars, we expected to obtain an initial spin-down power and spin-down time scale yielding an initial spin period and a magnetic field strength in agreement with population studies of radio pulsars (Faucher-Giguère & Kaspi 2006) and gamma-ray pulsars (Watters & Romani 2011). For supernova remnants, we expected E_ej and n₀ values in agreement with those inferred by Leahy et al. (2020) from a sample of Galactic X-ray SNRs, and we aimed at values for M_ej in line with those obtained in the explosion models of Sukhbold et al. (2016). Similarly, we aimed at getting particle injection properties consistent with the spread of values inferred from spectral studies of a number of established PWNe (Torres et al. 2014).

We emphasise that model complexity and computation time prevent a fine sampling of the large parameter space and a reliable estimation of the uncertainties on the model parameters. For that reason, the best-fit model setups presented below should be understood as sets of values yielding a satisfactory description of the experimental data, rather than measurements of the physical quantities involved. Also, we caution that the domain of applicability of the model (i.e., not too deep into reverse shock crushing of the nebula) tends to favour relatively large ejecta masses, as this delays the RS-PWN interaction, which in turn implies larger t_ST times and smaller $p_{\max }^{\mathrm{ST}}$ values. In practice, this drove the best-fit M_ej values to the maximum acceptable value according to Sukhbold et al. (2016), 15 M_⊙. Some of these issues may be solved in future developments of the model.

4.2. HESS J1809−193

4.2.1. Summary of observations

We first applied the model to HESS J1809−193, an unassociated gamma-ray source discovered in 2007 and possibly connected to PSR J1809−1917, a pulsar with period P = 0.083 s, characteristic age τ_c = 5.1 × 10⁴ yr and present-day spin-down power of Ė_PSR = 1.8 × 10³⁶ erg s⁻¹. The estimated distance to the pulsar is d = 3.3 kpc (Manchester et al. 2005). Alternative explanations or complementary contributions such as an SNR interacting with its molecular environment were put forwards (Castelletti et al. 2016; Araya 2018; Boxi & Gupta 2024; Li et al. 2023; HESS Collaboration 2023). A compilation of the multi-wavelength context and detection history can be found in HESS Collaboration (2023).

In the most recent analysis of H.E.S.S. observations by HESS Collaboration (2023), HESS J1809−193 can be resolved into two morphologically and spectrally distinct components: (i) an extended component A, described spatially by an elliptical 2D Gaussian intensity distribution (with a 1σ size of 0.6° and 0.3° along the major and minor axis respectively), and spectrally by a curved spectrum with a cutoff at 13 TeV; (ii) a more compact component B, described spatially as a symmetric 2D Gaussian intensity distribution (with a 1σ size of 0.1°) and spectrally by a flat power-law spectrum seemingly extending beyond a few tens of TeV without cutoff. PSR J1809−1917 and its X-ray nebula are offset from the centre of the compact component, at a distance of about its characteristic extent, which may result from a combination of asymmetric reverse-shock crushing and pulsar natal kick. HESS J1809−193 was detected up to 100 TeV and beyond from HAWC observations (Abeysekara et al. 2019). It is possibly extended at gamma-ray energies > 56 TeV, with a reported 1σ size of 0.3° for a symmetric 2D Gaussian morphology.

In the GeV range, an extended source known as 4FGL 1810.3−1925e was reported and is positionally very close to PSR J1809−1917 (Abdollahi et al. 2020; Araya 2018). It has an extension that is intermediate between that of the compact and extended components A and B observed at TeV energies (a 1σ size of 0.3° for a symmetric 2D Gaussian intensity distribution), and a relatively steep spectrum modelled as a simple power law above 1 GeV. It seems possible to describe the corresponding GeV emission assuming the morphology of TeV component A, which results in a modest increase in the flux suggesting that most of the emission is indeed captured with the intermediate size morphology (HESS Collaboration 2023).

4.2.2. Best-fit model setup

In the framework of our model, it seems natural to ascribe TeV component B to the PWN, while TeV component A could arise either from the SNR only, or from a combination of the SNR and the ISM components. The origin of the GeV emission is less obvious. The source has an extension comparable to that of component A, and a morphology apparently consistent with it, which suggests 4FGL 1810.3−1925e could be the GeV counterpart of the total emission from the system (PWN+SNR+ISM), but the relatively steep spectrum below 10 GeV does not seem a priori compatible with that interpretation (see below). Because of this uncertainty on the origin of the GeV emission, we searched for the parameter set best accounting for the TeV spectra only.

The best-fit parameters are given in Table 1. The best-fit model involves a system age t_age ≃ 16500 yr, at which time the PWN and SNR have a radius of 9 pc and 23 pc, respectively. The predicted size for the PWN is consistent with the 86% containment radius of TeV component B (a radius twice the sigma of the fitted 2D Gaussian intensity distribution). The model setup features M_ej = 15 M_⊙, which leads to a reverse-shock crushing occurring at 15 kyr, such that the PWN has been undergoing compression and/or disruption for about two thousand years. This is qualitatively consistent with X-ray observations showing a hard-spectrum nebula extending away from the pulsar in the direction of the TeV source peak (Klingler et al. 2020).

Enforcing t_age = τ_c − τ₀ yields an initial spin-down time scale of 34500 yr. Assuming a constant braking index of 3, applicable for magnetic dipole radiation, this implies an initial spin-down power of 3.9 × 10³⁶ erg s⁻¹, an initial spin period of 68 ms, and an initial magnetic field of 1.7 × 10¹² G. These initial properties agree with those inferred from the known pulsar population (Faucher-Giguère & Kaspi 2006; Watters & Romani 2011).

The particle injection spectrum for the PWN is a broken power law with index 1.2 below a break energy of 800 GeV, and index 2.4 above it, plus a fixed cutoff at 1 PeV. A relatively high cutoff is required by the flat spectrum of TeV component A. These values are typical of known young PWNe (Bucciantini et al. 2011; Torres et al. 2014), and so are the efficiencies with 90% of the pulsar power transferred to relativistic pairs and the remaining 10% injected predominantly into turbulence. The latter has a maximum spatial scale of 2% of the PWN radius, which is about 1/4 of the TS radius at the current age of the system.

After diffusive escape from the nebula, following Eq. (11), particles are confined in the remnant up to a maximum momentum of about $p_{\max }^{\mathrm{ST}}=150$ TeV/c at around t_ST ≃ 4500 yr and subsequently decreasing in time with a fixed power-law index ξ = 2.5. Such a value is consistent with the range inferred for γ-Cygni from a modeling of the gamma-ray emission in and near the remnant (MAGIC Collaboration 2023). Yet, the model fitting was performed on TeV data only, so $p_{\max }^{\mathrm{ST}}$ is not strongly constrained. Lower values by up to one order of magnitude can be considered before the fit to the TeV spectra starts to degrade.

The fitted spectrum in the upper panel of Fig. 9 shows that HESS J1809−193 is dominated by emission from particles spreading out in the ISM down to an energy of about 1 TeV. In that respect, the elongated shape of TeV component A, with a low inclination with respect to the Galactic plane, might indicate anisotropic propagation with faster diffusion along large-scale magnetic field lines that are preferentially aligned with the plane of the Galaxy. The model provides an explanation for the 10 − 100 GeV emission detected with the Fermi Large Area Telescope (LAT): in that range, the radiation has a flat or hard spectrum and comes primarily from particles that escaped into the SNR and are still trapped within it. The model cannot account for emission below 10 GeV, which therefore has to be ascribed to another component (for instance the radiation of cosmic rays accelerated at the forward shock, which would account for the steeper spectrum).

Fig. 9. Model fits to the spectrum of HESS J1809−193 and HESS J1825−137. For HESS J1809−193 (top panel), the model was fitted to the H.E.S.S. data points only. The emission below 10 GeV is attributed to another component not accounted for in our model framework (see text). For HESS J1825−137 (bottom panel), the model was fitted to Fermi-LAT and H.E.S.S. data points. Upper limits are not shown for clarity but model predictions are consistent with them. In each panel, PWN refers to the emission from particles still trapped in the shocked pulsar wind, SNR refers to the emission from particles that escaped into the remnant and are trapped downstream of the forward shock, and ISM refers to the emission from the halo of particles released in the surrounding medium.

An alternative interpretation exists for HESS J1809−193, in which the compact component B actually originates from hadronic cosmic rays interacting with nearby molecular clouds (HESS Collaboration 2023). A solution for such a mixed scenario can easily be found in the present model framework, for instance by using a higher κ_T parameter, thereby allowing more particle escape and depressing the spectrum of the PWN component, while other parameters are slightly adjusted so that component A still is described as an SNR+ISM component.

4.3. HESS J1825–137

4.3.1. Summary xnof observations

We further applied the model to HESS J1825−137. This gamma-ray source is classified as PWN, presumably powered by PSR J1826−1334 (also known as PSR B1823−13), a young pulsar with period P = 0.1015 s, a characteristic age τ_c = 2.14 × 10⁴ yr and a present-day spin-down power of Ė_PSR = 2.8 × 10³⁶ erg s⁻¹. The estimated distance to the pulsar is d = 4 kpc (Manchester et al. 2005). A compilation of the multi-wavelength context of the source can be found in Abdalla et al. (2019), and we will focus here on its gamma-ray properties.

In the TeV range, HESS J1825−137 is among the brightest sources known, with a flux of about 60% of the Crab nebula, and one of the most extended, with an emission detected over about 1.5° (Abdalla et al. 2019). The source is highly asymmetric and actually extends predominantly to the southwest of the pulsar. This is interpreted as resulting from an asymmetric reverse-shock interaction, in which the nebula was first hit and compressed on the north-eastern side. This possibility is supported by the discovery of molecular gas north of HESS J1825−137 and at a compatible distance (Lemiere et al. 2005): the evolution of the SNR in that direction would have occurred more rapidly, with an earlier formation of the reverse shock that pushed parts of the nebula to the southwest and constrained its subsequent evolution to occur in that direction. The TeV emission is markedly energy-dependent, with a spectrum softening with distance from the pulsar, which was interpreted as the progressive cooling of electrons and positrons as they are transported further and further away from the pulsar (Aharonian et al. 2006). The observed pattern is a rich data set that allows us to constrain the respective contributions of diffusion and advection, and the evolution of the magnetic field within the volume (Van Etten & Romani 2011; Collins et al. 2024). A direct consequence is that source morphology is also strongly energy-dependent, with an emission that shrinks around the pulsar with increasing gamma-ray energy (Abdalla et al. 2019). The source was also observed with HAWC and detected as extended (Albert et al. 2021), with a measured spectrum that is above the one inferred from the H.E.S.S. data, by a factor of about two at energies above 5 TeV.

In the GeV range, a very significant source is detected and found to be even more extended than in the TeV range (Principe et al. 2020). The Fermi-LAT observations extend the trend observed with H.E.S.S. and the source size further increases as photon energy decreases, up to nearly 2° in radius at 1 GeV. The centroid of the emission is also observed to evolve with energy, but the measured pattern is hard to interpret. Spectrally, emission in the Fermi-LAT band is hard and curved. The full SED from 1 GeV to 100 TeV is well described by a log-parabola or broken power-law shape, with a peak at 100 GeV.

4.3.2. Best-fit model setup

The very large size of HESS J1825−137 has been a challenge for most modelling attempts (Van Etten & Romani 2011; Khangulyan et al. 2018; Liu & Yan 2020; Collins et al. 2024). Assuming that the observed emission arises from shocked pulsar wind within the nebula implies that the PWN has a size of 80 pc at least and the SNR an even larger size probably exceeding 100 pc (for the southern part of the system). Both requirements can be achieved for a high explosion energy E_ej ≳ 3 × 10⁵¹ erg, in a very tenuous medium n₀ ≲ 10⁻³ H cm⁻³ (de Jager & Djannati-Ataï 2009; Van Etten & Romani 2011). While such values cannot be excluded, they are hardly compatible with the average properties inferred for Galactic SNRs (Leahy et al. 2020). In this work, we take a different approach and aim at a solution based on more typical parameters of the SNR. As we will see, given the properties of PSR J1826−1334, this yields a PWN much smaller than 80 pc. Unsurprisingly, the emission observed at large distances from the pulsar is therefore attributed to particles that leaked out of the nebula into the surrounding medium.

Because of the difference between the H.E.S.S. and HAWC spectrum for HESS J1825−137, we performed two different model fits: one on Fermi-LAT and H.E.S.S. data, and another one on Fermi-LAT and HAWC data. They yielded very similar parameters sets, the main difference being on pulsar properties. The fit to the Fermi-LAT and H.E.S.S. data is better, in terms of χ² given the number of degrees of freedom, and we discuss only this one in the following. The corresponding best-fit parameter set is summarised in Table 1 and the fitted spectrum is displayed in Fig. 9.

The best-fit model setup features a system age t_age ≃ 21000 yr, at which time the PWN and SNR have a radius of 26 pc and 33 pc, respectively. As illustrated in Fig. 9, the gamma-ray emission is dominated by particles spreading out in the ISM at energies below 30 TeV. As discussed further below, the observed extent of the GeV–TeV emission can therefore be used to constrain particle transport in the medium surrounding HESS J1825−137. As we approach 100 TeV, the contribution from the PWN becomes increasingly significant, which is consistent with the H.E.S.S. observation that the emission from a core region of 0.4° radius converges towards that of a larger region of 0.8° radius above ∼50 − 60 TeV (Abdalla et al. 2019).

The model setup features M_ej = 15 M_⊙, which leads to a reverse-shock crushing occurring at 12 kyr, such that the PWN has been undergoing compression and/or disruption for the past few millenia. We emphasise however that this may hold primarily for the south-western part of the system, while the north-eastern side may have experienced earlier compression due to the presence of molecular material there. We are reaching here the limit of the spherically symmetric model to describe an object that is obviously asymmetric. An initial spin-down time scale of 470 yr is found, which implies an initial spin-down power of 5.8 × 10³⁹ erg s⁻¹, an initial spin period of 15 ms, and an initial magnetic field of 3.3 × 10¹² G. These initial properties agree with those inferred from the known pulsar population (Faucher-Giguère & Kaspi 2006; Watters & Romani 2011). The short initial spin period implies that PSR J1826−1334 was originally very energetic, with a total rotational energy of 8 × 10⁴⁹ erg, which may be the reason why HESS J1825−137 is such a bright source now. Although this is comparable to the 5 × 10⁵⁰ erg kinetic energy of the stellar ejecta, it is unlikely that this have a significant impact on the dynamics of the SNR because more than half of the rotational energy has been radiated away by the nebula after about 1 kyr.

The particle injection spectrum for the PWN is a broken power law with index 1.0 below a break energy of 210 GeV, and index 2.5 above it, all values that are typical of known young PWNe (Bucciantini et al. 2011; Torres et al. 2014). The injection efficiencies are similar to those inferred for HESS J1809−193, with about 80% of the pulsar power transferred to relativistic pairs and the remaining 20% going preferentially into magnetic turbulence. The latter is injected with a maximum spatial scale of 1% of the PWN radius, which is about one fourth the radius of the TS at the current age of the system.

The fit to the spectrum, and the large extent of HESS J1825−137 over a large GeV–TeV range, imply a largely subdominant contribution from particles trapped in the remnant. In practice, we had to freeze the $p_{\max }^{\mathrm{ST}}$ parameter, the maximum particle energy reached at t_ST ≃ 14 kyr, to a low value of 1 TeV because it could not be well constrained in the fit. With the fixed parameter ξ = 2.5, this means that particles with energies above ≃360 GeV and above should have decoupled from the remnant at the current age of the system, otherwise emission below 10 GeV would be dominated by the SNR component. Overall, in the specific framework of the model, a low $p_{\max }^{\mathrm{ST}}$ makes it possible for particles to quickly escape into the ISM, which could be the reason for the puzzling very large extent of HESS J1825−137 over such a broad energy range. A more generic interpretation of this requirement is that the parent remnant of HESS J1825−137 was relatively inefficient at confining particles for some reason: partial disruption of the FS upon interaction with the surrounding material, uneven magnetic conditions over the FS surface, maybe in relation to the ambient large-scale magnetic field structure.

Since our model is not spatially resolved and consists in three concentric and uniform zones, it cannot be quantitatively tested against the measured evolution with distance of the 0.3 − 5 TeV photon index reported in Abdalla et al. (2019). At small angular distances from the pulsar < 0.3 − 0.4°, the 0.3 − 5 TeV surface brightness is most likely dominated by the PWN contribution that has a relatively flat average spectrum in agreement with observations. As we move to larger angular distances, the emission becomes dominated by the ISM contribution, and the slightly steeper spectrum predicted by the model is an average of the emission from particles increasingly affected by radiative cooling as they are transported away from the PWN, which eventually may be consistent with the observed trend of photon index increasing from about 2 to more than 3 between 0.3° and 1.5°. This needs to be verified from further development of the model.

Computing the spatial distribution from these particles that escaped the SNR-PWN system and are now experiencing transport in the surrounding ISM is beyond the scope of this study and will be addressed in a subsequent paper. Several works like those of Van Etten & Romani (2011), Liu & Yan (2020) and Collins et al. (2024) have shown the potential of HESS J1825−137 for probing the transport conditions in the medium around the source. We make below some qualitative remarks about this, taking into account the GeV results of Principe et al. (2020).

Under the assumption of purely diffusive three-dimensional transport, the typical range of propagating particles can be approximated by:

$$\begin{array}{c}\hfill r_{\mathrm{diff}}(E)=\sqrt{6D_{\mathrm{ISM}}(E)\times \min (t_{\mathrm{age}},t_{\mathrm{cool}}(E))},\end{array}$$(12)

where t_cool is the cooling time for particles in the ISM. For a diffusion coefficient with an energy dependence of the form D_ISM ∝ E^δ, and considering leptonic energy losses only in the Thomson limit, this yields r_diff ∝ E^δ/2 for t_age < t_cool (diffusion-limited regime) and r_diff ∝ E^{(δ − 1)/2} for t_age > t_cool (loss-limited regime). The characteristic extent of the gamma-ray emission decreases from about 105 to 15 pc as gamma-ray energy increases from 10 GeV to 10 TeV. Considering the above scalings, such a monotonic decrease with energy can be explained either by loss-limited diffusion across the whole energy range, or by a peculiar energy dependence of the diffusion coefficient (or a combination of both). The former case can be dismissed because 100 TeV particles (radiating typically at 1 − 10 TeV) are indeed in the loss-limited regime, with t_cool ≃ 4 × 10³ yr, but 100 GeV particles (radiating typically at 1 − 10 GeV) are not, with t_cool ≃ 10⁶ yr. The observed extent of the gamma-ray emission therefore suggests a diffusion coefficient somehow decreasing with energy. Interestingly, such a trend is predicted in calculations of the self-excitation of Alfvenic turbulence by cosmic rays escaping an SNR and streaming in the hot phase of the ISM (see Nava et al. 2019, especially their Fig. 4 at an age of 2 × 10⁴ yr). The levels of diffusion suppression do not quite match, though. At 100 GeV, observations indicate a characteristic propagation length of 105 pc over the system age of 2 × 10⁴ yr, which translates into a diffusion coefficient of 2.6 × 10²⁸ cm² s⁻¹. This is about the value assumed as large-scale average coefficient for the Galaxy and consistent with the small levels of diffusion suppression found at these energies in Nava et al. (2019). At 100 TeV, the characteristic propagation length is 15 pc over a cooling time of 4 × 10³ yr, which translates into a diffusion coefficient of 2.8 × 10²⁷ cm² s⁻¹. This is about a factor 300 below the large-scale average value for the Galaxy, and an order of magnitude below the level of diffusion suppression obtained in Nava et al. (2019).

5. Conclusions

We presented a multi-zone model for the gamma-ray emission from PWNe, taking into account the escape of particles out of the nebula into the parent remnant and subsequently to the surrounding ISM. In its current version, the application of the model framework should be restricted to early evolutionary stages, up to modest levels of compression of the PWN upon reverse-shock interaction. The extension to later evolutionary stages will be addressed in a subsequent work, involving numerical hydrodynamical simulations to determine the properties of the remnant.

Relativistic electron-positron pairs injected at the pulsar wind termination shock can escape the nebula as a result of advection and diffusion in the volume. The level and maximum spatial scale of magnetic turbulence in the nebula are the parameters controlling the amount of escape. Particles exiting the nebula are subsequently confined within the remnant because particle acceleration and magnetic field amplification at the forward shock transforms it into a magnetic barrier that keeps energetic particles below a certain maximum momentum downstream of the shock. This confinement is described through a simple time-dependent prescription for this maximum momentum, such that the escape of particles from the SNR into the surrounding ISM is controlled via the highest momentum reached in particle acceleration at the forward shock and the power-law index of its subsequent decay in time (see Eq. 11).

For a typical reference system and a wide range of turbulence properties in the nebula, particle escape losses exceed radiative losses after a few centuries. The GeV–TeV emission from the remnant is then comparable to that of the nebula in the case of high turbulence, or completely outshines it in the case of low turbulence. In our reference model setup, this emission dominates the pion-decay radiation from cosmic rays accelerated at the forward shock and advected downstream in the SNR. In the TeV and especially PeV range, the contribution from particles escaped into the surrounding ISM has a total flux possibly exceeding by far that of the SNR+PWN components. This is relevant for instruments like Tibet ASγ, HAWC, and LHAASO, whose capabilities to detect very extended sources above 10 TeV make them useful probes of particle escape and transport in the vicinity of the source.

We applied the model to sources HESS J1809−193 and HESS J1825−137 with the objective of accounting for the spectral and morphological properties of their gamma-ray emission. The solutions found involve supernova remnants and pulsars with parameters that are typical of the currently known populations, and a nebula that has been undergoing reverse-shock interaction over the past millenia. In the case of HESS J1809−193, the compact TeV component is ascribed to the PWN, while the more extended TeV component is mostly explained from particles escaped into the ISM. Particles that escaped the PWN but are still trapped in the SNR can account for most of the hard-spectrum signal detected in the 10 − 100 GeV range. In the case of HESS J1825−137, the broadband GeV–TeV emission is accounted for mostly from particles escaped into the ISM, while hard-spectrum emission from the PWN dominates the signal above 30 TeV. In both cases, high turbulence is required to account for compact hard-spectrum > 10 TeV emission. Despite that, the total emission is completely dominated by particles that escaped into the ISM. In the framework of our model, this implies that the parent remnants were relatively inefficient particle accelerators, reaching cosmic-ray energies of about 100 − 150 TeV at most, and maybe far less in the case of HESS J1825−137. A more general conclusion, independent of the practical implementation of particle transport in the SNR, is that the parent remnants were inefficient at confining pairs escaped from the PWN, for some reason possibly related to the forward shock structure or integrity. Eventually, the prediction of significant particle escape into the ISM provides a convenient explanation for the relatively large sizes of the gamma-ray emission in both systems (more than 30 pc in HESS J1809−193, and 100 pc in HESS J1825−137). In the case of HESS J1825−137, our best-fit model setup involves an initially very energetic pulsar, which may be the reason why the source is so bright now.

The model has some potential for the description of other sources exhibiting both compact and extended emission components, for instance MSH 15−52 (Tsirou et al. 2017), HESS J1813−178 (HESS Collaboration 2024), or HESS J1834−087 (Abramowski et al. 2015). It also provides in a consistent way the temporal and spectral properties of the flux of particles originally energized by the pulsar wind and ultimately released in the ISM. It can therefore be used to constrain the escape and transport of particles in the vicinity of pulsar-PWN-SNR systems from broadband gamma-ray observations, or in studies of the contribution of pulsar-related systems to the local electron and positron flux.

Acknowledgments

The authors acknowledge financial support by ANR through the GAMALO project under reference ANR-19-CE31-0014. This work has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Bibliographic Services. The preparation of the figures has made use of the following open-access software tools: Astropy (Astropy Collaboration 2013), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and SciPy (Virtanen et al. 2020). Pierrick Martin wishes to thank Sarah Recchia for useful discussions on particle transport, and Francisco Salesa Greus for providing us with the HAWC data points.

References

Appendix A: Detailed summary of original model

We provide here a description of the original model from GSZ09, including an introduction of all important quantities and formulae, repeated here for convenience, and clarifications for the few minor variations we implemented. The appendix is organized into the main components of the system.

A.1. Supernova remnant

The stellar ejecta interact with a cold circumstellar medium here assumed to have uniform mass density ρ₀ and negligible pressure P₀. This interaction first goes through a non-radiative stage, during which energy losses from thermal radiation of the system are negligible, and then moves towards a radiative stage, where these losses become dynamically important.

The whole evolution in an idealized setup with spherical symmetry has been the focus of numerous studies over the past decades. In the non-radiative stage, the system rapidly takes the following characteristic structure: the expanding ejecta drive a FS propagating in the ambient medium, while a RS propagates back in the ejecta; neglecting instabilities, the shocked interstellar medium and shocked ejecta, downstream of the forward and reverse shock respectively, are adjacent and separated by a contact discontinuity (CD). The non-radiative stage is extensively described in Truelove & McKee (1999, hereafter TM99), which provides useful analytical formulae for the dynamics of the FS and RS as a function of the main parameter of the system¹: the ejecta mass M_ej and kinetic energy E_ej, an ejecta density structure parameter n_ej, and the ambient medium density ρ₀.

The ejecta are initially in homologous expansion with velocity being proportional to radius, v_ej(r, t) = r/t (Eq. 19 in TM99), and they remain so until they interact with the RS. Their density structure ρ_ej(r, t) is defined as a uniform core and power-law envelope with power-law index n_ej (Eqs. 20 and 24 in TM99). The position of the core-envelope transition is set so that the total mass and kinetic energy of the ejecta, given the assumed density and velocity profiles, integrate to M_ej and E_ej. The outer envelope of the ejecta is assumed to be quite steep with n_ej > 5 (see the discussions and references in Truelove & McKee 1999; Bandiera et al. 2021), our adopted value being n_ej = 9.

The evolution of the FS in the initial so-called ejecta-dominated stage is described by Eqs. 75-76 with the parameters of Table 6 in TM99. As the mass swept-up by the FS becomes significant compared to the ejecta mass, and a significant fraction of the ejecta energy has been transferred to the ambient gas, the dynamics of the FS evolves towards a different scaling with time, described as an offset power-law form of the self-similar Sedov-Taylor solution, defined by Eq. 56 in TM99. The transition occurs at time t_ST given by Eq. 81, while values for different density structure parameters n_ej are listed in Table 6 in TM99. The dynamics of the RS is initially obtained from that of the FS via the parameter ℓ_ED, following Eq. 77 in TM99. Once the RS enters the core part of the ejecta, its position is described by Eq. 83 in TM99.

A.2. Pulsar wind

Meanwhile, the central pulsar is assumed to spin down as a result of magnetic dipole radiation. The spin-down luminosity Ė_PSR evolves in time as:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{PSR}}={\dot{E}}_{\mathrm{PSR},0}{(1+\frac{t}{{\tau }_0})}^{-\frac{n+1}{n-1}}\mathrm{with}n=3.\end{array}$$(A.1)

The initial spin-down luminosity and spin-down time scale are Ė₀ and τ₀, respectively, both of which are related to the initial values of the spin period and magnetic field P₀ and B₀:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{PSR},0}=\frac{2{\pi }^2{I}_{\mathrm{NS}}}{{\tau }_0{P}_0^2}=\frac{8{\pi }^4{B}_0^2{R}_{\mathrm{NS}}^6}{3c^3{P}_0^4}.\end{array}$$(A.2)

I_NS and R_NS are the moment of inertia and radius of the neutron star, respectively, for which we take typical values of 10⁴⁵ g cm² and 12 km.

Under the assumption that spin-down results from magnetic dipole radiation and that the perpendicular component of the magnetic dipole moment does not change significantly with time, the spin-down time scale is related to the pulsar true age t_age and its characteristic age τ_c as:

$$\begin{array}{c}\hfill t_{\mathrm{age}}={\tau }_{\mathrm{c}}-{\tau }_0=\frac{P}{2\dot{P}}-\frac{3c^3{P}_0^2{I}_{\mathrm{NS}}}{4{\pi }^2{B}_0^2{R}_{\mathrm{NS}}^6}.\end{array}$$(A.3)

In the limit that P ≫ P₀, the age of the pulsar can be approximated from the measurable quantity τ_c = P/2Ṗ.

The spin-down powers a relativistic outflow, the pulsar wind, that, at sufficiently large distance from the pulsar, is mostly composed of a toroidal magnetic field and electron-positron pairs. Upon interaction with the surrounding medium, the wind is halted at a TS, where particles are efficiently accelerated up to very high energies. The acceleration mechanism is not yet fully elucidated and may involve several sites and processes (diffusive shock acceleration at the TS, turbulent acceleration downstream of it, etc; see Bucciantini et al. 2011; Amato 2020), but observations indicate it drains a large fraction of the wind kinetic energy, of the order of several tens of percent (Zhang et al. 2008; Bucciantini et al. 2011; Torres et al. 2014).

We neglect the possibly complex radial and latitudinal dependence of the pulsar wind and assume an isotropic outflow. The TS forms where the wind ram pressure equals the pressure of the surrounding medium (Gaensler & Slane 2006):

$$\begin{array}{c}\hfill R_{\mathrm{TS}}(t)=\sqrt{\frac{{\dot{E}}_{\mathrm{PSR}}(t)}{4\pi c{P}_{\mathrm{PWN}}(t)}}.\end{array}$$(A.4)

The surrounding medium is initially the SNR, but very rapidly becomes the PWN itself, that is downstream of the shock.

Energy is injected in the nebula in the form of non-thermal pairs and magnetic field. In the original model:

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{inj}}(t)& ={\dot{E}}_{\mathrm{inj},\mathrm{B}}(t)+{\dot{E}}_{\mathrm{inj},\mathrm{e}}(t)={\dot{E}}_{\mathrm{PSR}}(t),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill & ={\eta }_{\mathrm{B}}(t){\dot{E}}_{\mathrm{PSR}}(t)+{\eta }_{\mathrm{e}}(t){\dot{E}}_{\mathrm{PSR}}(t).\hfill \end{array}$$(A.6)

We neglect a possible ionic component and assume constant injection efficiencies.

Pairs are assumed to have a spectrum at injection Q_inj, e, as a function of particle kinetic energy E, consisting of a broken power-law distribution with exponential cutoff:

$$\begin{array}{c}\hfill Q_{\mathrm{inj},\mathrm{e}}(E,t)=Q_0(t)\times e^{-E/E_{\mathrm{c}}}\times \{\begin{array}{cc}{(E/E_{\mathrm{b}})}^{-{\alpha }_1}\hfill & \mathrm{if}E\le E_{\mathrm{b}},\hfill \\ {(E/E_{\mathrm{b}})}^{-{\alpha }_2}\hfill & \mathrm{if}E>E_{\mathrm{b}}.\hfill \end{array}\end{array}$$(A.7)

The modelling of the broadband emission from observed PWNe suggest α₁ values in the range 1.0 − 2.0, α₂ values in the range 2.0 − 2.8, and E_b in the range 100 − 500 GeV (Bucciantini et al. 2011; Torres et al. 2014). We assume a constant cutoff energy E_c with typical values in the range 0.1 − 1 PeV. For relatively old systems, with ages of a few 100 kyr, the time evolution of the potential drop can be expected to shift the maximum attainable energy to below 100 TeV (de Oña Wilhelmi et al. 2022).

In practice, we consider a kinetic energy grid running from 0.1 GeV to 1 PeV, and the injection spectrum integrated over that range yields:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{inj},\mathrm{e}}(t)={\displaystyle {\int }_{0.1\mathrm{GeV}}^{1\mathrm{PeV}}}E{Q}_{\mathrm{inj},\mathrm{e}}(E,t)dE.\end{array}$$(A.8)

Since the pair injection efficiency is assumed constant, the pair injection spectrum is normalized at each time by the pulsar spin-down power.

A.3. Pulsar wind nebula

The PWN is the volume V_PWN excavated by the shocked pulsar wind as it expands. It is bounded on the inner side at a radius R_TS by the TS, and on the outer side at a radius R_PWN by a thin shell of swept-up ejecta material. The model treats the nebula as a single zone, in which quantities are assumed to be uniform.

The dynamics of this shell is controlled by its inertia and the pressure difference between the nebula on the inner side, and the stellar ejecta on the outer side. At early times, the innermost unshocked stellar ejecta are assumed to have a negligible pressure because of their fast expansion and rapid cooling. Later on, the reverse shock propagating back through the ejecta will eventually reach the nebula. The material ahead of the bounding shell is heated and compressed and the pressure immediately outside the nebula becomes non-negligible, and possibly well in excess of that in the nebula, thereby leading to a strong compression of the nebula. We assume that past this reverse-shock crushing occurring at time t_crush, the remnant rapidly settles into a state described by the Sedov-Taylor solution (but see Bandiera et al. 2023a). The conditions ahead of the shell are therefore:

$$\begin{array}{cc}& P_{\mathrm{ej}}(R_{\mathrm{PWN}},t)=\{\begin{array}{cc}0.0\hfill & \mathrm{if}t\le t_{\mathrm{crush}},\hfill \\ {\overline{P}}_{\mathrm{ST}}(R_{\mathrm{PWN}}/R_{\mathrm{FS}})P_{\mathrm{FS}}(t)\hfill & \mathrm{if}t>t_{\mathrm{crush}}.\hfill \end{array}\hfill \end{array}$$(A.9)

$$\begin{array}{cc}& {\rho }_{\mathrm{ej}}(R_{\mathrm{PWN}},t)=\{\begin{array}{cc}{\rho }_{\mathrm{ej}}(R_{\mathrm{PWN}},t)\hfill & \mathrm{if}t\le t_{\mathrm{crush}},\hfill \\ {\overline{\rho }}_{\mathrm{ST}}(R_{\mathrm{PWN}}/R_{\mathrm{FS}}).{\rho }_{\mathrm{FS}}(t)\hfill & \mathrm{if}t>t_{\mathrm{crush}}.\hfill \end{array}\hfill \end{array}$$(A.10)

Quantities P_FS(t) and ρ_FS(t) are the pressure and density immediately downstream of the FS and the terms with the ST subscripts are the self-similar Sedov-Taylor profiles, computed following Bandiera (1984). In practice, for typical values of the system’s parameters such as those used in Sect. 3, the reverse-shock crushing occurs at t_crush > t_ST > t_core, which justifies the above assumption.

The bounding shell mass M_PWN(t) is contributed to by ejecta material collected when shell velocity v_PWN(t) exceeds the ejecta velocity immediately ahead of the shell, v_ej(R_PWN, t) (we use here a corrected version of the original formula following footnote 1 in Bandiera et al. 2023a). The shell width is neglected, an assumption justified by one-dimensional hydrodynamical simulations that show a shell width of the order of two percent of the shell radius past a few centuries (Jun 1998). Two-dimensional calculations show that the shell is susceptible to a series of flow instabilities that lead to some mixing and thickening of the shell (Jun 1998; Porth et al. 2014). Past reverse-shock interaction, the thin-shell assumption becomes highly questionable (Bandiera et al. 2023a).

$$\begin{array}{c}\hfill \frac{d}{\mathit{dt}}\left(M_{\mathrm{PWN}}{v}_{\mathrm{PWN}}\right)=4\pi R_{\mathrm{PWN}}^2(P_{\mathrm{PWN}}-P_{\mathrm{ej}}(R_{\mathrm{PWN}})+C_{\mathrm{PWN}}).\end{array}$$(A.11)

The last term on the right-hand side is the contribution from the momentum of swept-up mass and is included only if v_PWN(t) > v_ej(R_PWN, t). The pressure inside the nebula P_PWN is contributed to by magnetic field and relativistic pairs:

$$\begin{array}{c}\hfill P_{\mathrm{PWN}}=P_{\mathrm{PWN},\mathrm{B}}+P_{\mathrm{PWN},\mathrm{e}}+P_{\mathrm{PWN},\mathrm{T}}.\end{array}$$(A.12)

The evolution of the magnetic pressure is computed following GSZ09, while the pressure from relativistic pairs derives from their total energy in the volume, assuming an adiabatic index of 4/3:

$$\begin{array}{c}\hfill P_{\mathrm{PWN},\mathrm{e}}(t)=\frac{E_{\mathrm{PWN},\mathrm{e}}(t)}{3V_{\mathrm{PWN}}(t)}=\frac{1}{3V_{\mathrm{PWN}}(t)}{\displaystyle {\int }_{0.1\mathrm{GeV}}^{1\mathrm{PeV}}}E{N}_{\mathrm{e}}(E,t)dE.\end{array}$$(A.13)

The spectral distribution of pairs at a given time, N_e(E, t), is obtained by solving a transport equation involving injection and continuous energy losses:

$$\begin{array}{c}\hfill \frac{\partial N_{\mathrm{e}}}{\partial t}+\frac{\partial {\dot{E}}_{\mathrm{loss},\mathrm{e}}{N}_{\mathrm{e}}}{\partial E}=Q_{\mathrm{inj},\mathrm{e}}.\end{array}$$(A.14)

The particle energy loss term Ė_loss,e includes radiative losses from inverse-Compton scattering and synchrotron radiation as well as losses from adiabatic expansion of the nebula. The expression for energy losses from radiative processes including inverse-Compton scattering in the Klein-Nishina regime was approximated following Moderski et al. (2005):

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{rad},\mathrm{e}}& ={\dot{E}}_{\mathit{IC}}+{\dot{E}}_{\mathit{syn}},\hfill \end{array}$$(A.15)

$$\begin{array}{cc}\hfill & =-\frac{4}{3}c{\sigma }_T{\gamma }^2[{\displaystyle \sum _i}\frac{U_i}{{(1+4\gamma {ϵ}_i)}^{3/2}}+\frac{B_{\mathrm{PWN}}^2}{8\pi }],\hfill \end{array}$$(A.16)

where σ_T is the Thomson cross-section and γ the Lorentz factor of the particle. The sum is performed over all components of the radiation field, described as graybodies with energy densities U_i and temperatures T_i, yielding normalized photon energies for the fields ϵ_i = 2.8k_BT_i/m_ec², where k_B is the Boltzmann constant, m_e the electron mass, and c the speed of light.

A.4. Numerical aspects

We implemented the model following Sect. 2.2 in GSZ09, using a logarithmic grid in both time and particle energy, and starting from the initial conditions exposed in their Appendix B. Given the rather crude numerical scheme, small enough time step is required, of the order of the few thousands steps per decade at least, in order to properly capture the compressions and re-expansions of the nebula after reverse-shock interaction (this has however a limited impact here since our results and discussions focus on the early stage, until modest levels of compression).

Appendix B: Investigation of time scales involved

We assessed some of the hypotheses subtending our model by computing the relevant time scales. As a reference, we use the dynamical time scale for the evolution of the nebula, which is:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dyna}}=\frac{R_{\mathrm{PWN}}}{v_{\mathrm{PWN}}}.\end{array}$$(B.1)

Fig. B.1. Time scales of the dynamical evolution of the nebula and the turbulence cascade in the different model setups.

We first check the assumption that turbulence is fully developed at all times. We computed the time scale for diffusion in wavenumber space, as a result of non-linear wave-wave interactions, which results in the so-called cascade of turbulence from small wavenumbers (large physical scales) to large wavenumbers (small physical scales). The corresponding coefficient for a Kolmogorov phenomenology is (Miller & Roberts 1995):

$$\begin{array}{c}\hfill D_{\mathrm{kk}}=0.052v_A{k}^{7/2}W{(k)}^{1/2}.\end{array}$$(B.2)

In steady-state, this diffusion coefficient yields a Kolmogorov spectrum, w(k)∝k^−5/3, with constant energy flux over the relevant wavenumber range. The corresponding time scale is:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{casc}}=\frac{k^2}{D_{\mathrm{kk}}}.\end{array}$$(B.3)

This time is a decreasing function of k for a Kolmogorov spectrum, τ_k ∝ k^−2/3, so we evaluated the time scale at the largest spatial scale. The results are plotted in Fig. B.1. Unsurprisingly, the cascade proceeds all the more rapidly that the level of turbulence is high and its largest spatial scale is small².

As can be seen from Fig. B.1, the turbulence cascade time scale is at least one order of magnitude smaller than the dynamical time scale in our model setups, which justifies our assumption of fully-developed turbulence at all times.

Appendix C: Radio synchrotron emission

Figures C.1 and C.2 show the radio synchrotron spectra of each emission component for the different model setups and at the same three times as before. The magnetic fields in the nebula at these times are in the range 120 − 180 μG, 70 − 110 μG, and 55 − 65 μG, respectively. In this band, the ISM component plays a very minor, if not completely negligible, role because the low-energy particles involved in the radiation are released from the remnant very late.

In the low turbulence cases, most of the leptons injected into the nebula have escaped it past 1 kyr, and they are henceforth trapped in the SNR. The corresponding synchrotron emission is therefore about the same in all three model setups (for the three values of κ_T), and will hardly change from 1 to 10 kyr. The residual lepton population still residing in the PWN depends significantly on the escape conditions (the values of κ_T), and so does the corresponding synchrotron emission. Despite this, the emission from the PWN is initially stronger than that of the SNR because of a much stronger magnetic field in the latter. As time goes by, continued escape from the nebula and decrease of the magnetic field cause the radio signal from the PWN to drop below that of the SNR. The situation is very similar in the high turbulence cases, except that the emission from the PWN is weakly dependent on the turbulence scale. The main reason for both behaviours is that advection then dominates the escape, such that the details of turbulence matter little.

The model also predicts synchrotron emission in the X-ray range but we do not present or discuss it here as the one-zone assumption for the nebula seems little appropriate for this. X-ray emission from PWNe is observed to be more concentrated around the pulsar than radio or gamma rays, in a region where the magnetic field is probably not representative of the average value in the nebula volume.

Fig. C.1. Radio synchrotron spectra of the PWN, SNR, and ISM components in low turbulence model setups at different ages.

Fig. C.2. Radio synchrotron spectra of the PWN, SNR, and ISM components in high turbulence model setups at different ages.

All Tables

All Figures

	Fig. 1. Evolution of the size of the PWN for the two extreme model setups. Also shown as reference are the radii for the forward and reverse shocks of the SNR. The curves for the other setups lie in between the red and purple curves. In this and subsequent figures, the gray shaded area denotes the strong compression phase that the model is not suited to describe.
In the text

	Fig. 2. Evolution of the radiation and escape luminosities in the PWN volume for the different model setups. Also shown as reference is the pulsar spin-down power.
In the text

	Fig. 3. Evolution of the energy content in the PWN volume for the different model setups. Also shown as reference is the cumulative energy injected by the pulsar and the cumulative energy lost to radiation in the nebula.
In the text

	Fig. 4. Evolution of the spatial diffusion coefficient at 100 TeV in the different model setups. For reference, the effective value of the diffusion coefficient at 100 TeV for large-scale transport in the Galactic disk is on the order of 1030 cm2 s−1.
In the text

	Fig. 5. Evolution of the gamma-ray luminosity in different bands for the PWN, SNR, and ISM components in four model setups.
In the text

	Fig. 6. Gamma-ray spectra of the PWN, SNR, and ISM components in four model setups at different ages. The gray dotted line corresponds to an order-of-magnitude estimate of the pion-decay emission from cosmic rays in the remnant (see text for details).
In the text

Fig. 7. Gamma-ray spectra of the PWN, SNR, and ISM components at 10 kyr in a high-turbulence case, as a function of the maximum particle momentum $p_{\max }^{\mathrm{ST}}$ allowed in the remnant. The PWN emission is the same in all three runs since it is not affected by this parameter. The gray dotted line corresponds to an order-of-magnitude estimate of the pion-decay emission from cosmic rays in the remnant (see text for details).

In the text

	Fig. 8. The 68% containment radius of the PWN+SNR emission as a function of energy, for the different model setups and at a system age of 10 kyr. The gray and black dotted lines show the angular extent of the PWN and SNR, respectively, for reference.
In the text

Fig. 9. Model fits to the spectrum of HESS J1809−193 and HESS J1825−137. For HESS J1809−193 (top panel), the model was fitted to the H.E.S.S. data points only. The emission below 10 GeV is attributed to another component not accounted for in our model framework (see text). For HESS J1825−137 (bottom panel), the model was fitted to Fermi-LAT and H.E.S.S. data points. Upper limits are not shown for clarity but model predictions are consistent with them. In each panel, PWN refers to the emission from particles still trapped in the shocked pulsar wind, SNR refers to the emission from particles that escaped into the remnant and are trapped downstream of the forward shock, and ISM refers to the emission from the halo of particles released in the surrounding medium.

In the text

	Fig. B.1. Time scales of the dynamical evolution of the nebula and the turbulence cascade in the different model setups.
In the text

	Fig. C.1. Radio synchrotron spectra of the PWN, SNR, and ISM components in low turbulence model setups at different ages.
In the text

	Fig. C.2. Radio synchrotron spectra of the PWN, SNR, and ISM components in high turbulence model setups at different ages.
In the text