Spatio-spectral-temporal Modelling of Two Young Pulsar Wind Nebulae

Kes 75 (Kesteven, 1968) is a prototypical composite SNR that contains a PWN powered by pulsar PSR J1846$-$0258 (hereafter J1846). It is also one of the youngest SNRs in the Galaxy. Becker & Helfand (1984) first argued that this is a “composite SNR” containing both a shell and a PWN radio component, the latter having a radio index $\alpha\sim 0.0$ ($S_{\nu}\propto\nu^{\alpha}$) and a highly polarised, centrally-peaked radio distribution, upon imaging this source with the Very Large Array (VLA) at 2, 6, and 20 cm. Salter et al. (1989) added data at 84 GHz and concluded there may be a spectral steepening between $15-84$ GHz ($\alpha\sim-0.3$). Its distance was estimated to be $d=5.8\pm 0.5$ kpc (Verbiest et al., 2012), but recently Ranasinghe & Leahy (2018) updated the estimate to $d=5.6\pm 0.3$ kpc.

Kes 75 has been detected at VHEs by the High Energy Stereoscopic System (H.E.S.S.), with a flux of 2.4 $\times$ 10^-12 erg cm^-2 s^-1 and a photon spectral index of $\Gamma\sim 2.3$ (Terrier et al., 2008). At these energies, the source manifests as a point source, and it is not clear what the relative flux contribution of the SNR shell versus that of the PWN might be (a similar observational situation is found in the detection by INTEGRAL; cf. McBride et al. 2008). However, the size of this point-like TeV source is $0.010\pm 0.013^{\circ}$, which is much smaller than the diameter of the SNR shell of $3^{\prime}$ but comparable to the PWN’s diameter of $40^{\prime\prime}$ in the radio band (H.E.S.S. Collaboration et al., 2018; Ng et al., 2008). Furthermore, assuming a distance of 6 kpc for this source, the estimated PWN spectral energy distribution (SED) and the derived conversion efficiency of the spin-down luminosity of the pulsar to the $\gamma$-ray emissions are consistent with the origin of a PWN (see, e.g., Torres et al., 2014; Tanaka & Takahara, 2011). Moreover, this TeV source is classified as a PWN in the TeVCat catalogue (Djannati-Ataï et al., 2008a). Therefore, we also consider this VHE source to be the TeV counterpart of PWN Kes 75.

In the GeV band, Straal et al. (2023) detected flux above 100 MeV which they considered to be a combination of both pulsar magnetospheric and PWN emission. They concluded that the emission in the range 100 MeV $-$ 2 GeV should originate from the pulsar’s magnetosphere, and they model the emission above 10 GeV as originating from the PWN. They regarded the point source with a log-parabola spectrum, 4FGL J1846.9$-$0247c, to be the counterpart of Kes 75. This has been the case up to and including the Data Release 3 of the Fourth Fermi Point Source Catalog (4FGL-DR3; Abdollahi et al., 2020, 2022). However, with the release of 4FGL-DR4 (Abdollahi et al., 2022; Ballet et al., 2023) and subsequent GeV PWN studies, a new faint point-like source, 4FGL J1846.4$-$0258, is now added as the counterpart to Kes 75 instead (private communication, J. Eagle). Also, since the GeV emission from Kes 75 is faint, and 4FGL J1846.9$-$0247c is at a distance of $>0.2^{\circ}$ from Kes 75, it can further be argued that the latter source is not associated with the Kes 75 PWN. In what follows, we indicate the published Fermi Large Area Telescope (LAT) data from Straal et al. (2023) and compare the SED prediction from our PWN model to these older data. Since our model does not include pulsed emission from the pulsar, we show their data below 2 GeV (interpreted to be from the pulsar) for representation purposes, but we do not model those data.

In the X-ray regime, the Chandra X-ray Space Telescope was, however, able to resolve the morphology of the source as well as infer its expansion from 2006 to 2016, tentatively detecting no discernible proper motion of the pulsar (Reynolds et al., 2018). Ng et al. (2008) found a jet-torus PWN structure in the X-rays, with the absence of an observable counterjet, yielding a lower limit on the flow speed of $0.4c$. On the other hand, Reynolds et al. (2018) found a northern and southern jet in the keV band, and observed great changes in their morphologies on a timescale of a decade, attributable to the short-term pulsar activity (See their Fig. 3). Absence of shell emission in the east points to a density gradient of the medium into which Kes 75 is expanding. Recently, Hu et al. (2022) analysed archival Chandra data for various sources (spanning $2000$ to $2016$ for Kes 75). They fitted spatial features of several sources using a pure diffusion model, suggesting the dominance of diffusion in the particle transport in the majority of these sources. To explain the discrepancy in consistency of parameters for a couple of their targets, they proposed the need for a spatially-dependent magnetic field, which is something we consider in this paper.

Kes 75’s embedded pulsar, J1846, was discovered as a rotation-powered pulsar in 2000 (Gotthelf et al., 2000) in the X-ray band, but no pulsations were detected in the radio band (Archibald et al., 2008) and also not significantly (4.2$\sigma$) in the $30-100$ MeV gamma-ray band (Kuiper et al., 2018) nor in the $>100$ MeV band (Straal et al., 2023), due to a lack of statistics. It is one of the most energetic pulsars known, with a period of $\sim 326$ ms and a spin-down luminosity of 8.1$\times 10^{36}$ erg s^-1. For a period derivative of 7.1$\times 10^{-12}$ s s^-1 (Gotthelf et al., 2000), it has a small characteristic age of $\tau_{\rm c}\equiv P/2\dot{P}\sim 728$ yr (for the canonical value of the braking index, $n=3$). Its inferred magnetic field of $5\times 10^{13}$ G places it on the verge of the magnetar population. Indeed, in 2006, this pulsar exhibited magnetar-like short X-ray outbursts (Gavriil et al., 2008; Kumar & Safi-Harb, 2008), with a subsequent softening of the X-ray spectrum (Gavriil et al., 2008) and a change in the pre-burst braking index value from $n=2.65\pm 0.01$ (Livingstone & Kaspi, 2006) to $n=2.16\pm 0.13$ (Livingstone et al., 2011). Archibald et al. (2015) found an index $n=2.19\pm 0.03$ over a seven-year period following the outburst (i.e., from 2006 to 2013), with evidence for a glitch that occurred between 2005 and 2008. After 14 years of quiescence a pulsar reactivation (outburst in August 2020) was reported by Krimm et al. (2020); Blumer et al. (2021). Hu et al. (2023) detected a large spin-up glitch in NICER (Neutron star Interior Composition ExploreR) monitoring data and an increase in pulsed X-ray flux, similar to the behaviour of the 2006 outburst. They reported a braking index of $n=2.7\pm 0.2$ before this outburst, indicating that between 2006 and 2020, the spin-down relaxed to the original evolution before the 2006 outburst. The braking index value post-2020 outburst has not been estimated yet due to a lack of data.

G21.5 has been known as a radio source in our Galaxy from the 1970s. It was later observed as an extended X-ray and radio source by the Einstein Observatory and radio instruments (Altenhoff et al., 1970; Becker & Szymkowiak, 1981). It is a composite SNR, where the central region is plasma-filled and powered by a pulsar (PSR J1833$-$1034), and it also has an expanding surrounding shell of material powered by the supernova blast wave. The period of the pulsar is $P=61.9$ ms and its derivative is $\dot{P}=2\times 10^{-13}$s s^-1 (Camilo et al., 2006). Its current spin-down power is $\simeq 3.4\times 10^{37}$ erg s^-1, ranking it as the fifth highest $\dot{E}_{\rm rot}$ pulsar in our Galaxy (Manchester et al., 2005). The characteristic age of the pulsar is $\tau_{\rm c}=4.9$ kyr (Gupta et al., 2005). The revised source distance is $d=4.4\pm 0.2$ (Ranasinghe & Leahy, 2018) compared to earlier estimations of $\sim 4.7$ kpc (Becker & Szymkowiak, 1981; Tian & Leahy, 2008).

The radio spectrum of G21.5 is similar to that of two other young PWNe: Crab and 3C 58. It is uniform over the circular projected image of the PWN of $\sim 1^{\prime}$ in extent (Bietenholz & Bartel, 2008). Using the expansion speed of the nebula, its age is estimated¹¹1If one assumes $\tau_{\rm exp}=P/(n-1)\dot{P}$, one finds $n\sim 1.4$. This is similar to the estimation made using the Giant Metrewave Radio Telescope of $n=1.857$ (Roy et al., 2012). to be $\tau_{\rm exp}=870^{+200}_{-150}$ years (Bietenholz & Bartel, 2008).

Infrared (IR) observations of this object were performed using the Very Large Telescope (VLT), the Canada-France-Hawaii Telescope, and the Spitzer Space Telescope by Zajczyk et al. (2012). Significant linear polarisation was found, indicating a non-thermal origin of the IR emission, which helps to constrain the synchrotron part of the spectrum. However, their observations are from a 2^′′ region around the central pulsar and also the flux levels were lower compared to the IR observations by Gallant & Tuffs (1999). In our modelling, we have considered the IR observations based on Gallant & Tuffs (1999), that also indicate that G21.5’s PWN has a morphology very similar to that of the Crab.

In the $0.1-10$ keV band, Chandra X-ray observations revealed a non-thermal emitting region of size of $40^{\prime\prime}$ in which the X-ray spectrum becomes steeper with radial distance from the pulsar (Matheson & Safi-Harb, 2005). NuSTAR also observed this object within a radius of $165^{\prime\prime}$, and reported that the spectra in the energy range of $3-45$ keV are better explained by a broken power-law with a break at $9.7\pm 1.3$ keV (Nynka et al., 2014). The spectral indices below and above the break energy are $1.996\pm 0.013$ and $2.093\pm 0.013$, respectively. However, Hitomi analysis selected a spatial extent of $140^{\prime\prime}$. Their observations in the energy range $0.8-80$ keV revealed a break at $7.1\pm 0.3$ keV, characterised by a lower spectral index value of $1.74\pm 0.02$ before the break energy. Beyond the break, the index value increased to $2.14\pm 0.01$, with the emission predominantly originating from the PWN. The exact cause of the lower value of the break energy remains unknown and warrants further investigation (Hitomi Collaboration et al., 2018). Recently, Hu et al. (2022) used $0.5-8$ keV Chandra X-ray archival observations (spanning $1999-2014$) to derive a spectral index of $\sim 1.4$ in the innermost region, steepening to $\sim 2.2$ at a radius of $35^{\prime\prime}$. This steepening might be connected with the injection of electrons from the central pulsar and SR cooling as they traverse the PWN magnetic field (Slane et al., 2000). This source is useful to cross-calibrate instruments, for example X-ray satellites, due to its lower brightness compared to the Crab PWN (Tsujimoto et al., 2011).

INTEGRAL observations suggest that there is a relatively small contribution ($<20\%$) from the pulsar in the $0-200$ keV energy range, and most of it is due to a PWN of size $40^{\prime\prime}$ (de Rosa et al., 2009). On the other hand, Fermi LAT observations reported no detection of GeV radiation from this object (Ackermann et al., 2011). TeV gamma-ray observations of this region identified the source HESS J1833$-$105 (Djannati-Ataï et al., 2008b). This object was also observed by H.E.S.S. during their Galactic Plane Survey and the emission properties were consistent with the previous observations (H.E.S.S. Collaboration et al., 2018). The H.E.S.S. source spectral shape $(\Gamma=2.1\pm 0.2)$ is very similar to the INTEGRAL hard X-ray spectral index of $\Gamma_{X}=2.2\pm 0.1$ (de Rosa et al., 2009).

The leptonic PWN model we use in this paper is basically the same as described in van Rensburg et al. (2018a), van Rensburg et al. (2018b), and van Rensburg (2020); the only update is adding an increase in the bulk flow speed $V_{0}$ after some time $t_{0}$ in the case of Kes 75 that may be due to the magnetar-like flares that deposit additional energy into the system. Martin et al. (2020) suggest that such bursts may lead to increased pair formation and hence an injection of additional energetic leptons into the PWN, or boost the mean strength of the nebular magnetic field. Such perturbations may yield increases in the flux in some radiative bands, depending on model parameters. Here, we investigate the scenario where the bulk flow is suddenly increased by a few percent, given the energy deposition via the bursts (See Appendix A). We furthermore exploit new model outputs and fit them to new data (expansion rate, time-dependent X-ray flux; Section 3.4), in addition to previous model outputs. In this section, we therefore summarise only the main elements of the spatio-temporal model. More details may be found in (van Rensburg et al., 2018b).

As a first step, we solve the following Fokker-Planck-type transport equation:

where $N_{\rm{e}}(\mathbf{r},E_{\rm{e}},t)$ is the number of leptons per unit energy and volume, V is the bulk velocity of leptons, $\kappa$ is the spatially-independent diffusion coefficient, $E_{\rm{e}}$ is the particle energy, $\dot{E}_{\rm{e,rad}}$ is the total, i.e., the sum of the SR and inverse Compton (IC) scattering, energy loss rates, and $Q$ is the particle injection spectrum. Here $r$ is the radial dimension (assuming spherical symmetry) and $t$ is the time since the PWN’s birth. We solve this transport equation in time and energy, including one spatial dimension.

The particle injection spectrum is assumed to be a broken power-law:

where $Q_{0}(t)$ is the time-dependent normalisation constant that is determined by equating the first moment of the injection spectrum to a constant fraction $\eta$ of the time-dependent pulsar spin-down luminosity $L(t)$, $E_{\rm{b}}$ is the injection spectral break energy, $E_{\rm{e,min}}$ ($E_{\rm{e,max}}$) is the minimum (maximum) energy of the injected leptons, and $\alpha_{1}$ and $\alpha_{2}$ are the spectral indices. To limit the number of free parameters in this model, we assume that $\alpha_{1}$ and $\alpha_{2}$ are time-independent. The choice of the injection spectrum as a broken power-law is motivated observationally and has been invoked in phenomenological models, where the value of $\alpha_{1}$ typically varies in the range $1.0-1.8$, and after the break the value of $\alpha_{2}$ is found to lie in the range $2.0-3.1$ (Tanaka & Takahara, 2011; Torres et al., 2013; Vorster et al., 2013b; Zhu et al., 2018). Since we are injecting the particle spectrum at the inner boundary of our radial grid, we can formally write $Q(\mathbf{r},E_{\rm{e}},t)=Q(E_{\rm{e}},t)\delta(r-r_{\rm min})$, where $r_{\rm min}$ is taken to be the PWN shock radius.

For the low-energy part of the broken power-law distribution, Tanaka & Asano (2017) found that the low-energy particles can originate from stochastic acceleration of relic particles or of particles injected at the interface of the supernova ejecta (Tanaka & Kashiyama, 2023; Tanaka & Ishizaki, 2024). Although the injection point of the low-energy particles is still a matter of contention (e.g., Aharonian & Atoyan 1996; Amato et al. 2000), the assumption that the low-energy particles are also injected from the PWN shock radius does not affect our results, since the low-energy particles are responsible for the radio emission of the spatially-integrated SED.

The energy losses in the model are due to radiative and adiabatic energy losses. The calculations for radiative losses (SR and IC scattering) are similar to the globular cluster modelling of Kopp et al. (2013) and the expressions can be found in Section $2.3$ of van Rensburg et al. (2018b).

The adiabatic losses caused by the bulk motion of the particles in the expanding PWN are given by the $\frac{1}{3}(\nabla\cdot\mathbf{V})$ term in Equation (1). Since this is a radiative model, not a magnetohydrodynamic (MHD) one, we parameterise the bulk flow speed profile of the leptons as well as the magnetic field’s spatial and temporal profile. The first is parameterised as $V(r)=V_{0}\left(r/r_{0}\right)^{\alpha_{\rm{V}}}$, where $V_{0}$ is the bulk flow velocity²²2As in Kennel & Coroniti (1984), we assume that $V_{0}$ is a constant. This means that the velocity profile only changes with space, and not with time. It may be interesting to consider the impact of $V_{0}(t)$ in future, but this will introduce more free parameters to the model, that would need to be constrained by, e.g., MHD modelling. at the termination shock at $r_{0}=r_{\rm min}$, and $\alpha_{\rm{V}}$ is the bulk-flow index parameter. We parameterise the latter as $B(r,t)=B_{\rm{age}}\left(r/r_{0}\right)^{\alpha_{\rm{B}}}\left(t/t_{\rm{age}}\right)^{\beta_{\rm{B}}},$ where $B_{\rm{age}}$ is the present-day magnetic field at $r=r_{0}$ and $t=t_{\rm{age}}$ ($t_{\rm{age}}$ is the age of PWN). With Kennel & Coroniti (1984), we assume that the magnetic field is toroidal (azimuthal) and the bulk flow is purely radial (further discussion of this assumption is provided in van Rensburg et al. 2018b). To limit the number of free parameters for this work, we assume a fixed $\beta_{B}=-1$ (the variation for $\beta_{B}$ will be explored in our follow-up works). In the free expansion stage of the PWNe, a range of values of $\beta_{B}$ ($-1.15$ to $-1.45$) have been used (Vorster et al., 2013b), while in other works, this parameter was fixed to $-1.3$ (Reynolds & Chevalier, 1984; Tanaka & Takahara, 2010). Varying $\beta_{B}$ in our model between this range did not substantially change the model outputs.

Since this implies that $B(r,t)$ becomes infinite at $t=0$ (i.e., for negative values of $\beta_{\rm B}$) in this approximate parameterised expression, we limit the maximum magnetic field to $10~{}B_{\rm age}$. To test the robustness of this prescription, we also tested the value $B_{\rm max}=30~{}B_{\rm age}$ and did not find significant changes to our model outputs, given our assumed value for current age, $t_{\rm age}$ (because the $B$-field is only very large for a very short time and for a small region in space). Raising this maximum value increases the computational time significantly while affecting the results only marginally, hence we consider this to be a reasonable limit, since we do not expect assumptions for $t_{\rm age}$ to change too much. The magnetic field variation for different radial distances over the first $50$ years of Kes 75’s lifespan is shown in Fig. 1.

We furthermore assume that, since the nebular plasma is a good conductor, we can apply the ideal MHD equations to describe the PWN wind such that the magnetic field is frozen into the outflowing plasma. In this case, Ohm’s law becomes

By combining this with Faraday’s law, we find (e.g., Ferreira & de Jager, 2008)

To link the radial profiles of the magnetic field and bulk motion of the particles, we assume with Schöck et al. (2010); Holler et al. (2012); Lu et al. (2017, 2019) that the temporal change in the magnetic field is slow enough that we can set $\partial\mathbf{B}/\partial t\simeq 0$ in the above equation, even for a time-dependent prescription of $B$. (This assumption holds exactly for steady-state models such as those of Kennel & Coroniti 1984; Vorster et al. 2013a; it may not be fully applicable for the first decade of the evolution, see Fig. 1, although this represents only a small fraction of the total age.) From this follows

which for our simplified parametric specifications of the magnetic field and bulk flow implies that

This relation is used to reduce the number of free parameters by one and to simplify our exploration.

The diffusion is assumed to be either Bohm-type,

with $\kappa_{\rm B}=c/3e$ and $e$ denoting the elementary charge, or a power-law type (e.g., Kolmogorov-type),

where $\kappa_{\rm X}$ and $\delta$ are the diffusion coefficient normalisation and exponent respectively, similar to what has been used in Galactic cosmic-ray studies. We present a motivation of the choices for $\delta$ below in Section 3.2, whereas $\kappa_{\rm X}$ is left as a free parameter to be constrained by the spatial and spectral data.

Although a considerable amount of theoretical work (see, e.g., Schlickeiser, 2002; Shalchi, 2009, 2020; Engelbrecht et al., 2022) has been done in the context of the heliosphere, very few studies of particle diffusion in PWNe have been performed (see, however, Tang & Chevalier, 2012; Vorster et al., 2013a; Porth et al., 2016; Olmi et al., 2016; Olmi & Bucciantini, 2019; Zhu et al., 2023; Lu et al., 2023). This renders studies of the relative influence of diffusion as a transport mechanism in these environments difficult. These difficulties are further exacerbated by the lack of observational constraints pertaining to both the large-scale behaviour of PWN plasma quantities as well as the nature of the plasma turbulence. The latter information is essential when modelling the diffusion coefficients of energetic charged particles.

However, some information as to the energy dependencies of these coefficients can be gleaned from existing theories, by making assumptions as to the nature of the turbulence, as well as the large-scale PWN magnetic field. Assuming the pulsar wind magnetic field (PWMF) is essentially frozen into the outflow of charged particles when $\sigma$ (the ratio of electromagnetic energy density to that of particles) is low (i.e., a high plasma-$\beta$ scenario), and a radial outflow of said particles, the PWMF can be assumed to spatially scale as a divergence-free Parker (1958) magnetic field, given by

where the Parker field is assumed to emanate from a source surface located at a distance $r=r_{0}$ (with a value there of $B_{0s}$), nominally where $\sigma$ drops below 1. Furthermore, $V_{0}$ denotes the pulsar wind speed, $\Omega$ is the rotation rate of the pulsar wind at $r_{0}$, and $\Psi$ is the winding angle between the magnetic field and the radial direction. As a 1D approach is being taken in this paper, as a first approach we can assume that the co-latitude $\theta=90^{\circ}$ so that, given the high rotation rate and pulsar wind speed, the above reduces to the familiar $r^{-1}$ radial dependence for the magnetic field magnitude (beyond the pulsar light cylinder, where the corotation speed equals the speed of light). Given the self-consistent constraint of Equation (6), this would mean that $\alpha_{B}=-1$ and $\alpha_{\rm V}=0$. This, of course, is only valid prior to any later shock in the plasma. For instance, in the heliosphere the stellar wind drops below supersonic speeds at $\sim 85$ AU generating the heliospheric termination shock, after which the Parker field would not adequately describe the heliospheric magnetic field (see, e.g., Kleimann et al., 2022, and references therein). Note that the presence of analogous termination shocks in PWNe are expected from theory and observation (see, e.g., Lyubarsky, 2003; Kirk et al., 2009; Slane, 2016; Cerutti & Giacinti, 2020). Assumptions about the PWN field structure leads to a constraint for parameter studies, particularly for $\delta$, which governs the energy dependence of the diffusion coefficient in Equation (8). For a pulsar, such a Parkerian field would be essentially azimuthal, with the implication that the radial diffusion coefficient, which is what is technically employed in our 1D transport model, would be dominated by the diffusion coefficient perpendicular to the magnetic field.

It should be noted, however, that the perpendicular diffusion coefficient would be expected to be considerably smaller than the diffusion coefficient parallel to the PWMF (see, e.g., Engelbrecht et al., 2022; Herbst et al., 2022). The strong background PWMF would imply the presence of strongly anisotropic transverse 2D turbulence, based on MHD simulations performed by Shebalin et al. (1983), which could in principle lead to effective cross-field pitch-angle scattering. In such a scenario, a reasonable model for the perpendicular diffusion coefficient would be the nonlinear guiding centre (NLGC) theory of Matthaeus et al. (2003). Cosmic-ray modulation studies employing diffusion coefficients derived from this theory have successfully reproduced spacecraft observations of Galactic cosmic-ray protons and antiprotons (e.g. Engelbrecht & Moloto, 2021), and said diffusion coefficients are in good agreement with numerical test-particle simulations in synthetic turbulence (e.g., Minnie et al., 2007). Although more complex expressions for electron diffusion coefficients have been proposed and employed in modulation studies (see, e.g., Engelbrecht, 2019; Dempers & Engelbrecht, 2020), the NLGC expression has the benefits of being relatively tractable as well as successful in reproducing heliospheric observations.

The NLGC theory, however, requires at the very least an estimate for the form of the 2D magnetic turbulence power spectrum, which is obviously currently impossible to ascertain from observations for PWN contexts. MHD simulations can provide useful insights, such as those of Porth et al. (2016), who analyse turbulence generated via their MHD model and report a spectrum with a wavenumber-independent energy-containing range, and a Kolmogorov inertial range. The NLGC diffusion coefficient also needs, at the very least, an estimate of the energy/rigidity dependence of the parallel mean free path. A reasonable approach, based on the assumption that highly energetic particles would only be resonantly scattered by turbulent fluctuations with characteristic length scales comparable to the particle Larmor radii, and that such fluctuations would be found within the wavenumber-independent energy-containing range of the turbulence power spectrum, would be to assume that the parallel mean free path scales as $R^{2}$, where $R$ denotes the particle’s rigidity, following what is expected from magnetostatic quasilinear theory (Jokipii, 1966; Teufel & Schlickeiser, 2003). The NLGC theory predicts a perpendicular mean free path that scales as $\lambda_{\perp}\sim\lambda_{\parallel}^{1/3}$, with $\lambda_{\parallel}$ the parallel mean free path, leading to a perpendicular diffusion coefficient that scales, at very high energies, as $R^{2/3}$, implying that $\delta=2/3$. Should, however, the particles in question be resonantly scattered by fluctuations in the inertial range, whether the spectral index of this be the Kolmogorov value ($-5/3$) or the Iroshnikov-Kraichnan value ($-5/2$), the parallel mean free path would scale as $R^{1/3}$ or $R^{1/2}$, respectively (see, e.g., Caballero-Lopez et al., 2019). This would then imply an NLGC perpendicular mean free path that scales as $R^{1/9}$ or $R^{1/6}$, implying that $\delta=1/9$ or $\delta=1/6$. In what follows, we investigate the effect of assuming different values of $\delta$, and also consider $\kappa_{\rm X}$ to be free, to find optimal by-eye fits to the data. See Table 2 and 3 for the best qualitative fit values of these parameters.

Our model is spherically-symmetric, spatially multi-zonal, and time-dependent, and is implemented in a code that calculates the transport and radiation of particles as they traverse a PWN. We use a spatial grid of concentric spherical shells of increasing radius, centred on the pulsar, which we call zones. The code spans three dimensions: space/radius, time, and energy. The radial binning is linear and is determined by $\delta_{\rm r}=\left(r_{\rm max}-r_{\rm min}\right)/\left(N-1\right)$, where $\delta_{\rm r}$ is the radial step, $r_{\rm min}$ and $r_{\rm max}$ are the chosen minimum and maximum values of $0.02$ and $10.0$ pc, respectively, and $N$ is the total number of radial bins. The results are a bit sensitive to the choice of $r_{\rm min}$. If the inner boundary is chosen to be too distant, we cannot reproduce some of the spatial data at lower radii. Therefore, we tested different values of $r_{\rm min}$ and found that $0.02$ is most suitable. We choose $r_{\rm max}=10.0$ pc, which is much larger than the observed radii of both our sources ($\sim 0.8$ pc), to ensure that the solution is not prematurely cut off due to an imposed boundary condition³³3While we assume that the PWN is in a freely-expanding stage, we effectively ignore the outer SNR conditions; there is therefore no hydrodynamic boundary in our model.. We tested that our spectral results are not too dependent on this choice by checking for convergence within acceptable limits for larger values of this parameter. The energy bins are calculated on a logarithmic scale, using $\delta_{\rm E}=\ln\left(E_{\rm max}/E_{\rm min}\right)/\left(M-1\right)$, where $\delta_{\rm E}$ is the logarithmic energy step, and $E_{\rm max}$ and $E_{\rm min}$ are the maximum and minimum lepton energy values (different for particles and also different for the various radiation mechanisms), and $M$ is the total number of energy bins. For a complete description of the binning procedure, refer to Appendix A.2 in van Rensburg (2020). We tested for convergence of our outputs for both the spatial and energy dimensions for a different number of bins. We found that beyond a combination of $N=180$ and $M=200$, there are not any appreciable differences in the SED and hence, we fix our number of bins at these values.

The simulation begins at $t=0$ representing the time of the birth of the PWN. The simulation is allowed to go up to the $t_{\rm age}$. The transport of the population of particles inside the PWNe is solved numerically using Equation (1) over the entire age of the PWN (van Rensburg et al., 2018a). The pulsar is allowed to spin down according to the usual braking formula (e.g., Pacini & Salvati, 1973) $L(t)=L_{0}(1+t/\tau_{0})^{-\gamma},$ with $\tau_{0}\equiv P_{0}/(n-1)\dot{P}_{0}=P/(n-1)\dot{P}-t_{\rm age}$ the characteristic age at birth. The canonical value of $n=3$ (implying $\gamma=(n+1)/(n-1)=2$) represents a dipolar magnetic field, and measurements yield typical values between $n\sim 2-3$ (see, for e.g., Gaensler & Slane 2006), although this value maybe be greater, or vastly greater when measured on small time scales (Archibald et al. 2016; Ekşi et al. 2016; Li & Gao 2023; Tong & Kou 2017 or Parthasarathy et al. 2020; Lower et al. 2021; Johnston & Karastergiou 2017, respectively). Such large braking indices may be connected, among other things, to changes in the pulsar moment of inertia, magnetic field decay, glitches, or the alignment of the magnetic and spin axes with time.

We lastly perform a line-of-sight (LOS) calculation (van Rensburg et al., 2018b) to project radiation from multiple spherical shells onto the plane of the sky, which enables us to estimate the SB and flux as a function of radial distance.

By solving Equation (1) numerically, we are able to predict the following PWN properties in a spherically-symmetric scenario:

• •

  Broadband SED: The SED shows the full radiation spectrum for the PWN, integrated over the spherical radial coordinate (or equally, over the cylindrical LOS-projected radius), spanning from radio, through X-rays to the TeV energy range, and including both an SR and an IC component, but neglecting the pulsed GeV component expected from the central pulsar.

• •

  Integral flux for different epochs: The integrated flux over a particular energy band is then calculated at different epochs (snapshots in time during the PWN evolution), but again integrated over all radii. In the case of Kes 75, this will be compared with the results in Table 3 of Reynolds et al. (2018), where the flux was spatially integrated over the whole PWN and the pulsar was excluded. It should be noted that the observed decrease in X-ray flux was not spatially uniform, with most of the decline associated with the dimming of the northern knot, thus violating our spherically symmetric model that we apply as a first approximation here. We therefore emphasise that we do not attempt to model a significant flux decrease limited to a particular region, as our model cannot resolve such a small spatial scale. But we can compare changes in the total flux from our model distributed over the whole PWN with that given in their paper.

• •

  SB profiles: LOS calculations are then performed by finding the intersection volume of radiating spheres and cylindrical⁴⁴4These should actually trace out cones, but given the vast distance to the PWN compared to its size, a cylindrical approximation is well motivated. lines of sight, thereby projecting the 3D PWN as a 2D image on the plane of the sky. This yields the predicted SB profile (as a function of projected radius $\rho$) at any snapshot in time, for a particular energy band. Details may be found in Section 2.9 of van Rensburg et al. (2018b). A significant update from previous works using this model is that instead of calculating the SB profile (and some other features) for a single energy value (the midpoint value of the relevant energy band), we integrate quantities over the energy band as specified by the observations (provided in Table 1).

• •

  PWN expansion for different epochs: Reynolds et al. (2018) calculated an average expansion of Kes 75 by taking account of changes in the physical image (i.e., spatial) and SB (i.e., flux) over time, with the expansion centred on the pulsar. They focused on the relatively faint northwestern rim that is far removed from the jets where significant morphological and flux variations were expected. Combining three independent expansion rate measurements over different epochs, they finally found an average expansion rate of $\sim 0.249\%\pm 0.023\%$ yr^-1. This corresponds to a PWN expansion velocity of $\sim 1000$ km s^-1, depending on the distance. To match these observations, we calculate the cumulative flux contained within a particular projected radius (specifically, the integral of $dN_{\gamma}/dE_{\gamma}$ over $0.7-8$ keV, normalised to the total flux found by integrating over all radii) and then use interpolation to find the $\rho_{68}$, which is the radius containing $68\%$ of this integrated flux⁵⁵5We fit the integrated flux profile using a Stretched Exponential Decay function of the form $a{\rm e}^{(-r/b)^{c}}$ to ensure smoothness of the profile for interpolation.. We use $\rho_{68}$ as the effective PWN size at a particular time, and compute the relative expansion of this radius over several different epochs. We compared our outputs to the expansion rates at different epochs (spanning $\sim 7$, $10$, and $16$ years), as listed in Table 2 of Reynolds et al. (2018).

• •

  X-ray photon index versus distance: The SR component of the emitted SED for each projected annulus on the sky (i.e., after the LOS calculation) is used to calculate the X-ray photon index as a function of $\rho$ by fitting a power-law function to the model SED in the X-ray energy band for several annuli centred on the pulsar. This is done over different epochs, and for different $\rho$ (lines of sight).

• •

  X-ray photon index versus epoch: We calculate the photon index as before, but for the SED integrated over all radii, taking a snapshot at a particular epoch.

• •

  X-ray photon index versus energy: We calculate the photon index as before but it is averaged over energy for specific energy intervals to be able to compare the X-ray part of our SED to the X-ray energy data from Hattori et al. (2020) in the case of G21.5.

Depending on data availability for each PWN, many of these predictions are simultaneously explored. Whether a 1D model parameter set represents data well or not, is judged by eye (i.e., not via a statistical prescription, but qualitative fits only) and the preferred model parameters are thus determined. We summarise the energy ranges for the data and their references in Table 1.

In this section, we show the influence of several free parameters in our model, corresponding to the present-day magnetic field, bulk-flow motion, diffusion, and braking index on various model outputs for Kes 75, and explain the physical reasons for the corresponding variations. Apart from the parameter under investigation, all the model parameters are fixed as listed in Table 2.

In the subsections below, we focus on describing the impact of varying the diffusion and braking index parameters on the model outputs. We do not show here the magnetic field variation study for G21.5, since the conclusions were found to be same as for Kes 75, discussed in Section 4.1.2.

Our preferred models for Kes 75 are shown in Fig. 13. The corresponding model parameters are provided in Table 2.

For the SED (see Fig. 13(a)), our predicted radio emission is slightly shifted in energy (the index is slightly hard), but it explains the X-ray and TeV emission properly. Fermi data is shown only for representation purposes (see Section 2.1). Some of the previous works like those of Bucciantini et al. (2011); Gotthelf et al. (2021); Straal et al. (2023) report excellent fitting for the SED of Kes 75, and we note that our parameter values are similar to those of the latter two works. The high-energy particle index ($\alpha_{2}$) in our model is preferred to be $2.4$, which is closer to the reported values of $2.17$ and $2.3$ by Gelfand et al. (2014) and Torres et al. (2014), respectively. We do not require a softer spectrum of $\alpha_{2}\sim 3.0$ as in the works of Gotthelf et al. (2021); Straal et al. (2023) (note that the latter study is also different from this work because of the inclusion of GeV data in their modelling). The slight difference in the quality of the model fit comes from the fact that we have fine-tuned our parameters to produce simultaneous model solutions for multiple (spatial) data features, and not just spectral data. See Section 5.3 below for further comparison to the best-fit parameters of a 0D code.

The predicted integrated X-ray flux and expansion rate over several epochs as shown in Fig. 13(b) and 13(c), are too low by a factor $\sim 2$. Since repeated measurements are done to improve precision, not simply to measure temporal evolution, we also calculate the average expansion rate per year to compare with a similar observational quantity. We calculate this from our model by dividing our predicted expansion rate (in percent) by the spans of the corresponding epochs ($16,10,$ and $7$ yr, respectively) and taking the mean of the three values. The average expansion rate we obtain is $0.104\%$ yr^-1, which is similar to the yearly expansion rates over individual epochs, and within a factor $\sim 2$ of the average measured value of $0.249\%\pm 0.023\%$ yr^-1, as given in Reynolds et al. (2018). Similarly, for the integral flux vs time, our model reflects the trend of no substantial time variation, as also seen in data, while overshooting the X-ray flux by a factor of $\sim 2$.

We find good correspondence between the model and the data for the (lack of) temporal and spatial variation of the X-ray photon index, as shown in Fig. 13(d) and 13(e). Within the dynamic range of the $y$-axis, and given the large error bars, the X-ray photon index does not seem to vary significantly, which is also reflected by the model results. Although the SB profile is underestimated by a factor $\sim 5$ at some radii, we show in Fig. 6(c) that it is well explained with a Bohm diffusion scenario; however, this does not explain other data sets.

In Fig. 10, the SED for G21.5 is shown, for possible diffusion scenarios within the PWN, as discussed in subsection 3.2. On the basis solely of the SED fit, it is difficult to identify the most favourable diffusion coefficient inside the PWN. However, for the power-law diffusion model, the normalisation $\kappa_{\rm X}$ at 1 TeV could be constrained to $\sim 10^{26}~{}{\rm cm^{2}\,s^{-1}}$ for $\delta\leq 1/6$ and it lowers by a factor of $\sim 3$ for $\delta=1/3$ and by $\sim 10$ for $\delta=2/3$, respectively. The normalisation $D_{0}$ at 1 TeV, for a value of $\delta=0.3$ in the ISM is approximately $2.5\times 10^{29}{\rm cm}^{2}\,{\rm s}^{-1}$ (Vladimirov et al., 2012). Hence we consider a diffusion model with $\kappa_{\rm X}\sim 10^{26}{\rm cm}^{2}\,{\rm s}^{-1}$ as a suitable choice, as also discussed in other works (Di Mauro et al., 2020).

In addition to the SED, we also model the SB and X-ray photon index profile (Fig. 11(a) and 11(b)). The SB profile is better explained by the power-law type diffusion parameterisation compared to the Bohm diffusion one. However, no parameter combination can reproduce the X-ray photon index profile below a radius $\sim 20^{\prime\prime}$, given the very hard photon index of $\sim 1.4-2.0$ compared to that of Kes 75, which is well explained by our model. For $\delta=1/6$, we reach the observed spectral index asymptotically, as shown in Fig. 11(b). Further, the normalisation of the diffusion coefficient was also found to be $~{}\sim 10^{26}~{}{\rm cm^{2}\,s^{-1}}$ by utilising the spatial variation of the SB and X-ray photon index data, as shown in Fig. 12. Our range of inferred model parameters based on the above model is listed in Table 3. Based on these model parameters, we estimated the integrated flux and X-ray photon indices over different epochs, which are shown in Fig. 14. The model outputs are shown in comparison to the available data, where two points are picked corresponding to the epochs 2012 and 2013, with observation IDs 4002 and 6003, as given in Table 4 of Hattori et al. (2020). We note that while Hattori et al. (2020) do present the results of observations at different times, they do not claim that the differences between measurements in flux and photon index in the same energy ranges are due to temporal variations of these quantities, but rather the result of statistical and systematic uncertainties in these measurements that may mimic a temporal variation. As such, we simply show our model outputs vs time in order to compare with these data, without claiming temporal variations. We do not show the intermediate data available for epochs within this one year, as we have not investigated the implications of such finer temporal details in this work. Within error bars, the average value of the NuSTAR X-ray flux is higher in the FPMB telescope and slightly lower in FPMA telescope, and forms an envelope for the model predictions, with rather small deviations. However, the average measured X-ray photon index is lower than what is predicted by the model, although this deviation is also rather small. We also show the variation of the X-ray photon index with respect to energy in Fig. 15. The data points are taken from Hattori et al. (2020) without a pulsar blackbody component, based on Chandra, Hitomi, and NuSTAR observations. As seen in the figure, our model captures the trend of the softening of the photon index with energy.

For comparison to the above results, we have used TIDE, a one-zone leptonic, time-dependent 0D model to describe the SED and dynamical evolution of a PWN (for a detailed description refer to Martín et al. 2012; Torres et al. 2014; Martín et al. 2016; Martin & Torres 2022). We use here the latest and more complete incarnation of the model as described in the latter reference. This allows us to produce a fit of several parameters at once, and in particular, to note whether there are alternative sets of parameters producing good fits, which does not appear to be the case. One should note, however, that there is a difference in the magnetic field description compared to what is used in the 1D model. Differently to what has been commented above, TIDE defines the $B$-field (independent of position) by solving a time-dependent equation that balances energy injection into the PWN and losses by expansion. This model difference may impact the comparison of the results of the 0D and 1D codes.

Considering the young ages of the pulsars, both PWNe we consider are believed to be in the free-expansion phase. Our 0D model could obtain a good fit to the SED data for these two PWNe. All model parameters, measured and assumed, along with fitted values and explored ranges, are presented in Table 4. Additionally, Fig. 16 shows the best-fit SEDs for PWN Kes 75 and G21.5. In the case of the latter PWN, we approximated the piecewise X-ray power law fits (Hattori et al., 2020) to the data by points with errors to facilitate the automated fitting process. The deduced radii of these two PWNe are in good agreement with the measured values, and their magnetic fields are also typical compared to what is found in other PWNe.

As outlined in Section 2.1, two X-ray outbursts were detected from PWN Kes 75 in 2006 and 2020 and they resulted in similar changes in the braking index values, decreasing from 2.65 (2.7) to 2.16 (2.19). The similarity between these changes seems to suggest the possibility of repeated outbursts. Typically, such variations in braking index values do not have a significant impact on SED and PWN properties when the true age of the pulsar, $t_{\rm age}$, is significantly smaller than the age that would lead to a negative initial spin-down age $\tau_{0}$ (e.g., Martín et al. 2016; De Sarkar et al. 2022). By adopting $n=2.16$, the fitted value of $t_{\rm age}$ for Kes 75 is 1.05 kyr, similar to what was used above (i.e., 728 yr). When assuming $n=2.65$, any value of $t_{\rm age}$ larger than 0.883 kyr will result in a negative $\tau_{0}$. For $n=2.65$, if we use 0.882 kyr as an upper limit when fitting for $t_{\rm age}$, we find that $t_{\rm age}$ is very close to this upper limit, resulting in an unusually small $\tau_{0}$ and an unusually large $L_{0}$. By imposing a fixed $t_{\rm age}$ value of 0.7 kyr, our 0D model would yield distinct fitted parameters from those in Table 4 and an $R_{\rm PWN}$ of approximately 0.6 pc, which is notably lower than the observed value of 0.814 pc. This mismatch may be due to the influence of plausible repeated outbursts, or to the deviation of 2.65 from the value of $n$ that dominates along the entire evolution. Both frameworks infer (or assume) an age of $t_{\rm age}\sim 1000$ yr for G21.5.

The preferred values of $\alpha_{1}$ and $\alpha_{2}$ for Kes 75 (see Table 2) are similar but outside the range suggested by the 0D code, and agree extremely well for G21.5 (see Table 3). The difference in the case of Kes 75 could be explained as due to the degeneracy of the model parameters - while it may possible to find a good SED fit within the limits suggested by TIDE, it might not describe the rest of the five spatio-temporal features we have considered in our preferred model fitting method. The inferred break energies are also quite similar between the two models. The preferred present-age magnetic field values for Kes 75 and G21.5 by the 1D model are $135$ and $140$ $\mu$G, while the 0D model suggests much smaller⁷⁷7The first value is the $B$-field at $r_{0}$; the volume-averaged value for the 1D model is of the order of $\sim 10\mu$G for $R_{\rm PWN}\sim 0.8$ pc, which is closer to the values preferred by the 0D model, although not exactly the same. The value varies depending on the choice of the $R_{\rm PWN}$. values of $\sim 30$ and $\sim 47$ $\mu$G. These values should be compared to the equipartition magnetic fields of $\sim 40$ $\mu$G for Kes 75 (Ng et al., 2008) and $\sim 300$ $\mu$G for G21.5 (de Jager et al., 2008). The difference between the preferred magnetic field values between the two models could be attributed to the difference in the magnetic field prescription between the two approaches. The considered magnetic energy fractions and soft-photon temperatures and energy densities are a bit different, indicating some freedom in these parameters. Finally, the 1D model can constrain model parameters connected with spatial aspects, such as the diffusion coefficient, bulk flow speed, and $B$-field profile, while the 0D model may infer dynamically-linked quantities such as the ejected mass and PWN radius that the 1D model does not consider. We thus see a broad concordance in the predicted SEDs of the two models, as well as in the values of the overlapping model parameters, despite there also being differences as discussed above.