© ESO 2020

1. Introduction

Most of the spin-down luminosity of young and very energetic pulsars is carried away in a magnetised wind of charged particles. The confinement of this particle wind outflow, which is predominantly composed of electron-positron pairs, leads to the development of a pulsar wind nebula (PWN). PWNe form when the particle wind collides with its surroundings, especially the slowly-expanding ejecta of the progenitor supernova, and forms a termination shock (Gaensler & Slane 2006). Evolved PWNe are ideal candidates for investigating particle transport mechanisms and electron cooling inside celestial object, due to their large γ-ray extension and possible energy-dependent morphology, which can provide spatially resolved spectra under certain circumstances. PWNe have been observed to emit photons up to TeV energies and, with more than 35 detections, they dominate the population of TeV gamma-ray sources in the Galactic plane (see Fig. 1). Despite deep observations of several PWNe, many open questions still remain; in particular, the mechanism by which the particles are accelerated at the termination shock is not yet understood (Slane et al. 2017).

Fig. 1. Sky map, in Galactic coordinates and Mollweide projection, showing the sources in the TevCat catalogue (http://tevcat.uchicago.edu, version April 2019) classified by their most likely association. All the 3FHL sources (Ajello et al. 2017) are also plotted, with grey points, for a comparison.

Among such objects, HESS J1825-137 is the largest and one of the most TeV efficient γ-ray PWNe currently known (H.E.S.S. Collaboration 2018, 2019). It is powered by a young and very energetic pulsar PSR J1826-1334 (also known as PSR B1823-13), which was discovered by Clifton et al. (1992). The pulsar has characteristics very similar to the Vela pulsar: it has a spin period of 101.48 ms, a characteristic age of 21 kyr, a spin-down energy of 2.8 × 10³⁸ erg s⁻¹, and is at a distance of 3.9 ± 0.4 kpc (Manchester et al. 2005; Cordes & Lazio 2002).

The first detection of the extended nebula was made by Finley et al. (1996) using X-ray observations with ROSAT, which observed a compact nebula of ∼20″ radius around the pulsar. Subsequent X-ray observations made by Chandra and Suzaku revealed its asymmetric morphology and an extended emission up to 15′ (17 pc) (Uchiyama et al. 2009). The discovery of the energy-dependent morphology of HESS J1825-137 at TeV energies (Aharonian et al. 2006) provided important proof that the emission is dominated by ‘relic’ electrons from the earlier epochs of the nebula in which the pulsar was spinning down more rapidly, therefore releasing more energy into the system. More recently, Abeysekara et al. (2020) has shown that HESS J1825-137, together with eHWCJ1908+063 and eHWCJ2019+368, are the only three sources detected above 100 TeV by the HAWC experiment.

In the GeV regime the source was detected first in 2011 by the Large Area Telescope (LAT, Atwood et al. 2009), on board the Fermi Gamma-ray Space Telescope, whilst at MeV energies the source is not significantly detected (Principe et al. 2018). Previous LAT analyses of the PWN HESS J1825-137 have been performed using 20 months of Pass 6 data in the 1–100 GeV energy band (Grondin et al. 2011) and subsequently six years of Pass 8 data in the 10 GeV–1 TeV energy band (Ackermann et al. 2017a). Taking advantage of more than 11 years of Fermi-LAT data now available, we have performed an analysis of the energy-dependent morphology and of the spectral parameters of the source in the energy range between 1 GeV and 1 TeV.

2. Fermi-LAT data and analysis

The LAT is a γ-ray telescope that detects photons by conversion into electron-positron pairs and has an operational energy range from 20 MeV to 2 TeV. It is comprised of a high-resolution converter tracker (for direction measurement of the incident γ-rays), a CsI(Tl) crystal calorimeter (for energy measurement), and an anti-coincidence detector to identify the background of charged particles (Atwood et al. 2009).

2.1. Data selection

For the LAT analysis of HESS J1825-137 we used 11.6 years of Pass 8 (P8R3) Source class events (Atwood et al. 2013; Bruel et al. 2018) collected between August 4, 2008, and March 20, 2020 (Fermi Mission Elapsed Time 239587201 s – 606355205 s) in the energy range between 1 GeV and 1 TeV. The data were taken in a region of interest (ROI) of radius 15° and centred on the source position given in Ackermann et al. (2017a). In the following we often use the term ‘PWN’ to refer to the source HESS J1825-137. We created sky maps with a pixel size of 0.1°. In order to eliminate most of the contamination from secondary γ-rays from the Earth’s limb, we excluded γ-rays with zenith angle larger than 105° (Abdo et al. 2009). We used the P8R3_Source_V2 instrument response functions (IRFs).

2.2. Region modelling

The model used to describe the sky includes all point-like and extended LAT sources, within 20° degrees of the source position, listed in the fourth Fermi-LAT source catalogue (4FGL, Abdollahi et al. 2020), as well as the Galactic diffuse and isotropic emission. We modelled the Galactic diffuse emission using two different templates, repeating the analysis for each, in order to study the systematic uncertainty due to the choice of the diffuse model. The first template used (D1) was the Galactic and isotropic diffuse templates¹ (labelled in our analysis as “D1”) used in the 4FGL model. For a crosscheck, we used as a second the template the diffuse model² derived in Ackermann et al. (2017b, labelled in our analysis as “D2”) which was especially developed for analysis of extended emission near the Galactic centre. The residual background and extragalactic radiation were described by a single isotropic component with the spectral shape in the tabulated model iso_P8R3_SOURCE_V2_v01.txt. For the analysis results presented here, we used the first template for the diffuse model unless specified otherwise.

2.3. Analysis procedure

The analysis was performed with Fermipy³ (version 0.17.4, Wood et al. 2017), a python package that facilitates analysis of data from the LAT with the Fermi Science Tools, of which the version 11-07-00 was used. In our analysis we applied the correction for the energy dispersion, as implemented in Fermipy, disabling it for the isotropic model.

In addition to the study of the general characteristics of the source (localisation, averaged extension and spectra) using the complete energy range between 1 GeV and 1 TeV, in this work, we studied also the PWN’s morphology in smaller energy bands. The description of the model used to describe the data is presented in Sect. 2.2. We modelled the PWN using a 2D-Gaussian model for the spatial template and a LogParabola spectral model

$$\begin{array}{c}\hfill ({\displaystyle \frac{\mathit{dN}}{\mathit{dE}}}=N_0{\left(\frac{E}{E_0}\right)}^{-[\alpha +\beta log(E/E_0)]});\end{array}$$(1)

as used in the 4FGL, as well as in Ackermann et al. (2017a); H.E.S.S. Collaboration (2019).

After a preliminary optimisation (fermipy.optimize) and fit (fermipy.fit) of the parameters of sources included in the model, we investigated (using fermipy.find_src) the possible presence of additional faint sources, not in the 4FGL catalogue, and we found three new candidate sources that we added to our model. The best-fit positions of these new sources are RA, Dec = (18^h13^m38^s, −17°49′47″),(18^h18^m28^s, −9°56′23″), and (18^h29^m34^s, −16°14′23″), with 95% confidence-level uncertainty R₉₅ = 6′44″. We verified also the possible influence of PSR J1826-1256 (4FGL J1826.1-1256), a bright Fermi-LAT source with significant emission at the energies considered in this analysis and that lies at 1° from the PWN centre. The gamma-ray steady emission from a point source spatially coincident to the pulsar position is significantly detected in our analysis (TS = 5285) with a high flux (F_{E >  1 GeV} = (5.55 ± 0.06 × 10⁻⁸ ph cm⁻² s⁻¹). We verified the possible influence of the pulsar by performing an analysis using only off-phase data, characterised in the forthcoming third Fermi-LAT pulsar catalogue using gamma-ray pulsar timing as detailed in Kerr et al. (2015). We found that the morphology and spectral results are compatible within the errors to the results obtained using the full phase data. Furthermore we did not see in our analysis (1 GeV–1 TeV) any extended residual around the PSR location that would be possibly associated to the VHE extended source HESS J1826-130 seen by Angüner et al. (2017). The emission of HESS J1826-130 is, however, detected by H.E.S.S. only at energies above 2 TeV, which are not covered by Fermi-LAT.

Following this, in order to improve the spatial and spectral modelling of HESS J1825-137, we performed a localisation (fermipy.localize), extension (fermipy.extension) and spectral (fermipy.sed) analysis of the PWN for the entire energy range. For the spectra derivation using the entire energy range, we divided all the photons between 1 GeV and 1 TeV in 24 energy bins (8 per decade) logarithmically distributed. We left all spectral parameters of the diffuse background as well as those of the sources within a 3° radius from our target free to vary. For the sources in a radius between 3° and 6° away, only the normalisation was fitted, while we fixed the spectral parameters of all the sources within the ROI at larger angular distances from our target.

Subsequently, we studied the energy-dependent morphology of the PWN, for which we grouped the photons into five energy bands: four logarithmically spaced bands between 1 and 100 GeV and one band between 100 GeV and 1 TeV. In the analysis of the extension variation with energy (see Sect. 3.2.1), we fixed the spectra of HESS J1825-137 (except the normalisation), as well as those of all the other nearby sources, using the resulting model from the full energy range spectral analysis. Only the normalisation of the diffuse and isotropic backgrounds were re-fitted during the analysis of the energy dependent morphology. Finally, we generated the SED by doing a spaced spectral analysis in each energy band with their corresponding spatial model.

2.4. Systematic uncertainties due to the emission models used

For the specific goal of a further investigation of the bias introduced by the choice of the diffuse model, we also performed the analysis with other two different models: a template used in Acero et al. (2016a, z = 4, ts = 150) and the LAT diffuse emission ring-hybrid model gll_iem_v06.fits (Acero et al. 2016b), and we compared the results of the four different models used in order to have an estimate of the systematic error due to the diffuse model. Comparing the results of the morphology for the entire energy range 1 GeV–1 TeV, we found a deviation of the source localisation of ∼0.08°, which is similar to the LAT PSF size at 10 GeV (i.e. ∼0.1°), and we considered it as an estimate of the systematic error due to the diffuse model. This deviation is more pronounced, ∼0.17° at low energy (E <  10 GeV), where the diffuse emission is much brighter than at high energy. Similarly, during preliminary studies of this source (Principe 2019; Principe et al. 2019), we performed the analysis modelling the ROI using the sources contained in the FL8Y⁴ which has been derived using standard LAT diffuse emission ring-hybrid model, and we compared the results. The resulting extension and spectral parameters obtained with these different background models are compatible within the uncertainties.

3. Fermi-LAT results for the entire GeV domain

In this section we present the results of the spectral and morphological analysis of the PWN HESS J1825-137 using 11.6 years of Fermi-LAT data between 1 GeV and 1 TeV. We also compare the morphological results with previous LAT analyses and combine the spectra obtained in this work with that published in H.E.S.S. Collaboration (2019).

3.1. Localisation and extension analysis for whole energy range

We performed the localisation and extension analysis using the entire energy range 1 GeV–1 TeV. For the determination of the extension, we investigated the source radius starting from a value of 0.05° (similar to the best PSF reached), which corresponds to the case of a point source, up to a maximum radius of 2.5°, in 41 linearly separated steps. Figure 2 shows the test statistic⁵ (TS) map of the region around the PWN for the energy range between 1 GeV and 1 TeV. The TS was evaluated by placing a point source at the centre of each pixel, Galactic diffuse emission and nearby sources being included in the background model. In the map, the size obtained for the PWN in this work is compared with the radius obtained in the previous analysis with Fermi-LAT data.

Fig. 2. TS map (in sigma units), in celestial coordinates, of the region around the PWN HESS J1825-137 in the energy range between 1 GeV and 1 TeV. The red circle and star indicate the 2D-Gaussian extension and centroid fit obtained in this work. The red dashed circles mark the uncertainty on the 2D-Gaussian extension. The green and white circles correspond to the extension obtained in the FGES catalogue (Ackermann et al. 2017a) and in Grondin et al. (2011), respectively. The black point indicates the position of PSR J1826-1334, the pulsar which is believed to power the nebula. The 4FGL sources, as well as the three candidate sources added in the model, are represented with light grey points.

The best-fitting radius obtained in this work is larger than those obtained in previous works (FGES catalogue, Ackermann et al. 2017a; Grondin et al. 2011). This is probably connected to the larger amount of data, down to 1 GeV, which allowed the more extended emission below 10 GeV to be resolved (see Sect. 3.2.1). Although the pulsar PSR J1826-1334 is detected in radio and X-rays (Duvidovich et al. 2019), its emission is not significantly detected by Fermi-LAT yet. A search for the pulsar was performed looking at a possible steady emission (no pulsation search has been done) and no significant point source emission from this position was found between 1 GeV and 1 TeV, only extended emission from the nebula is detected.

Table 1 reports the results of the localisation and extension analysis performed in the energy range between 1 GeV and 1 TeV using a 2D-Gaussian model as spatial template for the sources. The extension result reported here corresponds to the 68% containment radius.

The resulting centre position of the PWN is shifted by about 0.48° from PSR J1826-1334, the pulsar associated with the nebula. This asymmetry of the nebula extension with respect to the PWN position is related to the pulsar proper motion (discussed also in Sect. 5.1) as well as to the presence of a dense molecular cloud towards the north of the nebula. EVLA observation at 1.4 GHz made by Castelletti et al. (2012) reveals, in fact, a nearby molecular cloud with a density of ∼400 cm⁻³. At all wavelengths the nebula emission is observed to extend towards the south of the pulsar. The reason could be that in the past the external part of the supernova shell interacted with the nearby molecular cloud, leading to a relatively fast formation of a reverse shock on the northern side of the nebula. The recoil of this reverse shock forced the nebula to expand more towards the southern side.

3.2. Energy dependent analysis of HESS J1825-137

3.2.1. Energy-dependent extension analysis with a 2D-Gaussian template

During the analysis of the energy-dependent extension, we fixed the spectra of HESS J1825-137 and of the nearby sources using the resulting model from the initial spectral analysis on the whole energy range (see Sect. 2.3). The normalisation of the isotropic plus diffuse components are instead left as free parameters of the fit. In this part, we performed the localisation and extension analysis separately in each energy band, using 2 bands per decade between 1 and 100 GeV, and a single band between 100 GeV and 1 TeV, keeping the internally 8 bins per decade in the science tools analysis. Before performing the extension analysis, for each energy band the localisation is again optimised. The extension is then estimated by fitting a 2D-Gaussian template in each energy band and simultaneously refitting the source position. We estimated the systematic error of the extension σ_R, sys, due to the choice of the diffuse model, as:

$$\begin{array}{c}\hfill {\sigma }_{R,\mathrm{sys}}=\sqrt{{(R_{\mathrm{D}1}-R_{\mathrm{D}2})}^2+{({\sigma }_{R_{\mathrm{D}1}}-{\sigma }_{R_{\mathrm{D}2}})}^2+{\mathrm{\Delta }}^2},\end{array}$$(2)

where Δ is the distance between the source centroids for the two different diffuse models that we considered in our analysis. Table 2 contains the results of the extent measurements performed with the 2D-Gaussian template in 2 logarithmic bands per decade between 1 and 100 GeV and in a single band between 100 GeV and 1 TeV.

The extension estimates (R) obtained with the two different diffuse models (D1 and D2), are compatible for each energy band. Similarly the position of the source centroid in the two models are compatible or, in any case, they deviate by a distance which is smaller than the Fermi-LAT PSF at the corresponding energy. The PWN position is observed to vary with energy between the different bands. Moving from the lowest energy band, 1–3 GeV, to the highest energy one, 100 GeV–1 TeV, the fitted 2D-Gaussian centroid of the PWN moves towards the current position of the pulsar (see Fig. 3).

Fig. 3. TS maps (in sigma units), in celestial coordinates, of the region around HESS J1825-137 for the energy bands: 1–3 GeV (top left), 3–10 GeV (top centre), 10–32 GeV (top right), 32–100 GeV (bottom left) and 100 GeV–1 TeV (bottom right). The TS maps are smoothed with a Gaussian of radius 0.1°. The white circles represent the extension (solid line) and its statistical uncertainty (dashed lines) determined in the respective energy band. For comparison, the resulting extension obtained for the entire energy range (1 GeV–1 TeV) is overlaid with a red line. All extensions shown correspond to the 68% containment radius. In the 100 GeV–1 TeV (bottom right) plot, H.E.S.S. significance contours at 5, 10, and 15σ, for energies below 1 TeV (H.E.S.S. Collaboration 2019), are shown with light-blue lines for comparison. The H.E.S.S. contour also includes the excess of the nearby LS 5039 source (small circular excess at the southern boarder of the PWN).

3.2.2. Energy-dependent extension analysis with the radial profile method

For comparison, we additionally perform measurements of the nebula extent using the same approach as that adopted in H.E.S.S. Collaboration (2019). We estimated the extent of the nebula as a function of energy as the radial distance at which the emission in the southern half of the nebula drops to a factor 1/e relative to the maximum, starting from the position of the pulsar PSR J1826-1334 which powers the system.

The emission is considered only in one hemisphere due to the strong asymmetry of the source. The orientation of the semi-circular region is the same as that used in the H.E.S.S. analysis, with the major axis of the emission orientated along an angle of 208° with respect to the direction of positive declination, as shown in Fig. 4.

Fig. 4. Left: TS map (in counts) of the PWN HESS J1825-137 in the energy range 1 GeV–1 TeV with Fermi-LAT, with the region used to extract the radial profile (as used by H.E.S.S.) overlaid in white. The preferred emission direction (“major axis”) as found by H.E.S.S., along which the extent is evaluated, is indicated by the black dashed line. The position of the pulsar (black) and best-fit Fermi-LAT centre of a 2D Gaussian (magenta) are also indicated. Right: radial profile of the excess counts fit with Eq. (3) beyond the peak emission. The characteristic R(1/e) size of the nebula is indicated by a white (black) dashed line in the left (right) hand plot.

We estimated the extent of the emission by fitting a polynomial to the radial profile from a minimum radius, out to 4°, using the formula:

$$\begin{array}{c}\hfill y(x)=\{\begin{array}{cc}a{(r-r_0)}^n+c,\hfill & (x<r_0)\hfill \\ c,\hfill & (x\ge r_0)\hfill \end{array}\end{array}$$(3)

such that with increasing r the emission decreases out to a distance r₀ at which it approaches the constant value c. The minimum radius of the fit was chosen to be beyond the emission peak, the position of which was determined by a moving average approach. The position of the maximum emission was found to be offset from the pulsar and to vary with energy. The parameter a provides the overall normalisation, whilst the fitted value r₀ defines the maximum extent (see Table A.1 for the values of the fit parameters).

Since the distance r₀ was found to be highly sensitive to the order of the polynomial n, to mitigate this effect, the extension was taken as the radius at which the fitted function dropped to a fraction (1/e) of the peak value, R_1/e, as a characteristic extent of the nebula. This parameter, or indeed any other fixed fraction of the peak value, was found to be stable to the arbitrary choice of the polynomial index n, as tested in the range n = 2 − 5. We chose a value of n = 4 in Eq. (3). We evaluated the errors on the extension by performing the fit procedure on 1000 Monte Carlo realisations of the radial profiles in each energy band. The radial profiles were generated according to a random number selected from assuming a symmetric Gaussian distribution for each point, with width corresponding to the error on the point. The nature of the fitted function naturally results in an asymmetric error with larger extensions being more difficult to arise from statistical fluctuations, errors are therefore represented by the 16th and 84th percentiles of the extension distribution.

3.2.3. Comparison of the 2D-Gaussian and radial profile extent estimates

The extension of the PWN in the LAT energy range obtained by Fermipy is given as the 68% containment radius (R_68%) of a 2D symmetric Gaussian. In order to compare the two methods it is possible to approximate the polynomial function used for the radial profile analysis as a simple Gaussian without introducing a large bias. In this case, for a single Gaussian, the radius at which the function drops to the ratio 1/e of the peak value corresponds to $R(1/e)=\sqrt{2}\sigma$. Consequently, we can compare the extent results obtained with the two different methods (2D-Gaussian and radial profile) using the relation:

$$\begin{array}{c}\hfill R(1/e)=\sqrt{2}\frac{R_{68\%}}{\sqrt{-2log(1-0.68)}}=0.937R_{68\%}.\end{array}$$(4)

Since the extent obtained with the radial profile method is estimated starting from the pulsar position and considering the semi-circular region oriented along the major axis, to directly compare the results of two methods, we found from the 2D-Gaussian results the distance from the pulsar along the major axis at which the R(1/e) of Eq. (4) intersects the major axis. This correction accounts for the differences in the 2D-Gauss Ext. between Tables 2 and 3.

Table 3 presents the extension results obtained with the two different methods: 2D-Gaussian template and radial profile analysis.

The results from the 2D-Gaussian, including the offset of the centroid from the pulsar position, are broadly compatible with the results obtained from the radial profile (see Fig. 5), except for the energy band between 10 and 30 GeV. The differences may be due to the different regions considered in the extension analysis: the radial profile method takes the asymmetry of the source into consideration and performs the analysis only in the southern hemisphere, whilst the 2D-Gaussian assumes a symmetric source emission. Another possible reason is the different treatment of the extension as R(1/e); in the case of the radial profiles, the peak value is obtained from the excess counts independent of the fitted function; however, for the 2D-Gaussian, R_68% and the derived R(1/e) relies on the normalisation of the fit. If the normalisation of the Gaussian under- or overestimates the true peak value in the excess counts, which may occur due to the smoother curvature of the Gaussian function, then the distance R(1/e) is shifted accordingly.

Fig. 5. Nebula extent as a function of energy showing results from H.E.S.S. Collaboration (2019) and this analysis.

Our Fermi-LAT results extend the information of the energy-dependent extension of the PWN down to an energy of 1 GeV. The Fermi-LAT and H.E.S.S. extents for the common energy band (100 GeV–1 TeV) are compatible. The Fermi-LAT results reveal a continuous increase of the nebula extent towards lower energies, with a trend similar to that seen by H.E.S.S. above 100 GeV. The possible turnover observed by H.E.S.S. around 300 GeV appears to be ruled out by the LAT results at lower energies, leaving to the possibility of a turnover around few GeV.

3.3. Fermi-LAT spectral energy distribution

For the analysis of the source spectrum, we divided the photon events into 24 logarithmically distributed energy bins between 1 GeV and 1 TeV. We generated the spectral energy distribution (SED), using a 2DGaussian template for the source, by doing a spectral analysis in the various energy bands with their respective spatial parameters (see Table 2 in Sect. 3.2.1). We corrected for the jumps in the spatial model at the boundaries of the larger energy bands, by slightly enlarging the energy bands and averaging the flux estimates for the overlapped energy bins, as well as by increasing the error bars on the original results to include the uncertainty on the spatial modelling in the interim energy bands. The diffuse background was were left free to vary.

The SED is fitted with a LogParabola model (see Eq. (1)) as well as with a Broken Power Law (Broken PL) function, which is not a very natural or smooth model, but provides a better estimate of the break energy, E_b. The function used for the Broken PL law is:

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathit{dN}}{\mathit{dE}}}=N_0\times \{\begin{array}{c}{\left(\frac{E}{E_b}\right)}^{-{\mathrm{\Gamma }}_1}\mathrm{if}E<E_b\hfill \\ {\left(\frac{E}{E_b}\right)}^{-{\mathrm{\Gamma }}_2}\mathrm{otherwise}.\hfill \end{array}\end{array}$$(5)

The fit results for the LogParabola and Broken PL models are reported in the “Fermi” column of Table 4. The “Fermi + H.E.S.S.” column of the Table 4 contains the resulting spectral information for the PWN HESS J1825-137 from a joint fit to the energy flux points derived independently from Fermi-LAT and H.E.S.S. data (H.E.S.S. Collaboration 2019). The LogParabola is found to be preferred (χ²/ndf ∼ 1, with ndf = 35) for both the Fermi only and the combined spectral fits, and it nicely describes the spectra (see Fig. 6 for the combined SED). The Broken PL model, which fails to describe completely (χ²/ndf >  2) the spectral results above 10 TeV due to sharp cutoff, returns a energy break estimate of about 115 GeV.

Fig. 6. Combined spectra of the PWN HESS J1825-137 with the spectral measurements obtained in this work (red points) using 11.6 years of Fermi-LAT data from 1 GeV to 1 TeV and the H.E.S.S. results for the 100 GeV–90 TeV energy range (black points). The Fermi-LAT flux points were obtained doing a spectral analysis in the various energy bands with their relative spatial model. The combined SED has been fitted with both a LogParabola (blue line) and a Broken PL (green line). The vertical line corresponds to the energy break Eb of the Broken PL model. The bottom panel shows the normalised residual between the data and the LogParabola model.

4. Modelling of the Nebula

We modelled the combined GeV-TeV spectral energy distribution of the nebula as Inverse Compton (IC) scattering from a leptonic particle population. Several packages exist for spectral modelling with complementary features. We use both the NAIMA package (Zabalza 2015) for a single zone model due to the statistical fitting methods available, and the modular GAMERA package (Hahn 2016) for a multi-zone model due to the enhanced flexibility offered. The two models were found to be consistent except for an offset in the flux and magnetic field constraint, which arise from the different ambient radiation fields used. Whereas a lookup table for the total radiation fields from the model of Popescu et al. (2017) could be used directly in GAMERA, a black-body approximation to this model was used with the NAIMA package.

4.1. SED modelling with NAIMA

To test if a simple analytic electron distribution can explain the combined spectra of H.E.S.S. and Fermi-LAT, we used the NAIMA python package. We considered a leptonic population of particles, from an energy of 1 GeV up to 510 TeV, producing γ-ray by IC scattering. The radiation fields at the galactic position of HESS J1825-137 are obtained using the parametrization of Popescu et al. (2017). They are modelled as five black body components, with different temperatures and energy densities ε that correspond to Cosmic Microwave Background (CMB, T = 2.72 K, ε = 0.26 eV cm⁻³), Far Infra-Red (IFR, T = 43.07 K, ε = 0.78 eV cm⁻³), Infra-Red (IF, T = 238.4 K, ε = 0.17 eV cm⁻³), VISible (VIS, T = 3493 K, ε = 1.8 eV cm⁻³), and Ultra-Violet (UV, T = 19840 K, ε = 0.17 eV cm⁻³). The dominant contributors to the IC γ-ray flux are the FIR for energy between 3 GeV and 15 TeV, and the CMB otherwise. The contribution of the UV is negligible compared to the other radiation fields. The particle population is assumed to follow a broken power-law, its parameters are fitted using the MCMC method implemented in NAIMA. The normalisation of the broken power-law is derived from the total energy of the electrons W_e assuming that the source is located at a distance of 4 kpc. The results are presented in Fig. 7 and the parameters of the model are given in Table 5. This single population model is able to explain the complete range of energy covered by H.E.S.S. and Fermi-LAT data. We used this distribution in order to model the synchrotron component assuming different values of the magnetic field which we compare with X-ray data obtained by Suzaku (Uchiyama et al. 2009). The X-ray flux point corresponds to the sum of the X-ray emission from a comparatively small region (≲15′) around the pulsar, re-scaled for a comparable opening angle to that of the γ-ray analyses, with the X-ray emission assumed negligible outside of this region. This can be treated as an upper limit on the total X-ray flux averaged over the much larger region probed by the γ-ray emission. We found that the maximum magnetic field allowed by X-ray observations using this model is ∼4 μG.

Fig. 7. Results of the fitted IC NAIMA SED to the Fermi-LAT data and H.E.S.S. data, compared with the measured Fermi-LAT (red circle) and H.E.S.S. (black diamond) data points. A synchrotron component is computed from the electron distribution supposing 4, 5, 6, 7, and 8 μG, and compared with Suzaku data (blue square). The electron distribution parameters are presented in Table 5.

4.2. Multi-zone modelling with GAMERA

4.2.1. SED modelling

To simultaneously describe the SED (Fig. 6) and the variation in extent with energy (Fig. 5), we attempted a multi-zone modelling approach using the GAMERA package (Hahn 2016). A leptonic particle population producing γ-rays by IC scattering was again assumed, with the radiation fields at the location of the PWN obtained from Popescu et al. (2017) as for the NAIMA modelling; however, the black-body approximation is not necessary with the GAMERA package enabling a more precise parameterisation to be used. The difference this introduces is small, yet allows slightly higher magnetic field strength values of 5 − 6 μG.

We used a series of particle zones, added with burst-like injection at different times and left these free to evolve in time until the system age is reached. For describing the total γ-ray emission from the nebula at the current time we used a summation of 20 zones of particles of different ages, evenly split in time. The model parameters used for the curves shown in

Fig. 8 are given in Table 6. Several parameters were constrained to match current values at the present day, with the pulsar characteristic age assumed to be the age of the nebula system. The evolution of the pulsar spin-down luminosity L with time t is described by:

$$\begin{array}{c}\hfill L(t)=(1-\eta )L_0{(1+\frac{t}{{\tau }_0})}^{-\frac{n+1}{n-1}},\end{array}$$(6)

Fig. 8. Total SED of the nebula, combining X-ray and γ-ray data, is described by the summation of 20 zones of particles of different ages using the GAMERA modelling package Hahn (2016). The total SED of each zone is shown by a coloured line, from yellow for the youngest and smallest zone, through to blue for the oldest and largest. Model parameters are given in Table 6.

where η accounts for the conversion efficiency, n is the braking index for which a value of 3 was assumed, corresponding to pure magnetic dipole radiation, and τ₀ is the initial spin-down timescale of the pulsar. The spin-period P of the pulsar evolves as:

$$\begin{array}{c}\hfill P=P_0{(1+\frac{t}{{\tau }_0})}^{\frac{1}{n-1}},\end{array}$$(7)

from which τ₀ can be determined for an assumed initial spin-down period P₀ (a free parameter of the model) using the constrained parameters at the present day age of the system. For the electron population we used a broken power law injection spectrum.

The average magnetic field of the nebula was assumed to be 5 μG at the present day, evolving in time as:

$$\begin{array}{c}\hfill B(t)\propto {(1+\frac{t}{{\tau }_0})}^{-1}.\end{array}$$(8)

As the average magnetic field found for the whole nebula from the models is at ∼4 − 5 μG close to the average ISM magnetic field strength of ∼3 μG, any spatial dependence is likely to be weak. In reality, however, the magnetic field is expected to exhibit both temporal and spatial dependence, reducing with increasing distance from the pulsar within the nebula, with the strongest evolution in the region nearest the pulsar. Some spatial dependence of the B-field is included in this multi-zone model as a consequence of the different ages and sizes of the emission zones. The resulting SED is shown in Fig. 8, which is consistent with the available data.

4.2.2. Radial extent modelling

In a second step we described the radial extent as a function of energy. We stored the spectra for each zone of the model separately such that the zones could be arbitrarily arranged for the spatial modelling. The zones were treated as expanding shells in space, initially spherically symmetric, with the particle spectra filling the shell volume. That is, at a given radial distance r from the pulsar, the line-of-sight depth d_z through a given zone z of radial size R_z is given by:

$$\begin{array}{c}\hfill d_z=\{\begin{array}{cc}2\sqrt{R_z^2-r^2}\hfill & (r<R_z)\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\end{array}$$(9)

such that the depth through the zone along the line-of-sight decreases towards the edge of the zone and is zero at r >  R_z. The contribution from the spectrum of a given zone to the total emission at a distance r from the pulsar was therefore weighted by d_z.

We used the projection of the emission along the line of sight to form an emission profile from the model similar to that of Fig. 4. This projection was made in multiple energy bands, summing the relevant parts of the zone spectra. For determining the radial extent in each energy band we applied the aforementioned procedure (see Sect. 3.2.2) of fitting the radial profile to evaluate the distance at which the emission drops to a fraction 1/e of the peak value from the model.

The energy-dependence of the radial extent could be described using a radially dependent velocity profile v(r, t) with index β in the range [0.5,0.75] and initial outflow velocity v₀ = 0.03c in an advection dominated scenario:

$$\begin{array}{c}\hfill v(r,t)=v_0{\left(\frac{r}{r_{\max }}\right)}^{\beta }{\left(\frac{t}{T}\right)}^{-\beta },\end{array}$$(10)

where r_max is the maximum extent of the nebula, for which a value of 150 pc was assumed, corresponding to a maximum angular extent of ∼2° based on the extent measurements presented here. The spatial evolution of each of the zones was calculated from a minimum radius of 0.01 pc, corresponding to typical termination shock radii. The consistency of this range of β values with the experimentally measured extents is shown in Fig. 9.

Fig. 9. Radial extent of the nebula as a function of energy; the results of the model using GAMERA is indicated by grey shaded region for a compatible range of β, the index of the radially and time dependency of the velocity profile, between 0.5 (upper edge) and 0.75 (lower edge).

The energy dependent extent of the nebula obtained from this multi-zone model is seen to be approximately constant below a γ-ray energy of ∼0.1 TeV. Assuming that the energy dependent morphology of the nebula is due to electron ageing and cooling, then for un-cooled electrons, the nebula size remains constant in a spherically symmetric scenario.

Whilst a range of values of β compatible with previous models of the nebula are consistent with our results, a single β value for the full nebula seems not to be able to describe the data, implying either a changing velocity profile over time or non-spherical symmetry.

5. Discussion

We analysed 11.6 years of Fermi-LAT data in the energy range between 1 GeV and 1 TeV, and we performed, for the first time with Fermi-LAT data, an energy dependent analysis of the morphology of the PWN HESS J1825-137 in the GeV range.

5.1. System evolution

The nebula HESS J1825-137 presents a strong energy dependent morphology. The spectral index was seen by Aharonian et al. (2006), at TeV energies, to soften with increasing distance from the PSR; this implies that the population of electrons in the nebula had travelled and cooled out to large distances. The fact that the low-energy electrons at large distances from the pulsar produce the softest spectrum was interpreted to mean that these are the oldest electrons in the system.

In this work, we extended the observations down to 1 GeV confirming the emission at large distances by the lowest energetic particles. Figure 10 shows the size of the PWN for different energy bands. The centroid of the PWN moves with energy and, at high energies (above 30 GeV), it lies on the H.E.S.S. major emission axis. The centre position of the PWN seems to move from low to high energy emission, in a direction similar to the pulsar proper motion, which transverse velocity is v_⊥ = 440 km s⁻¹ (Pavlov et al. 2008) (see arrows in Fig. 10). The larger distance between the centroid of the PWN above 100 GeV (youngest electron population) and below 10 GeV (oldest electron population) with respect to the expected change in pulsar position (∼0.14°) due to the proper motion of the pulsar over the estimated characteristic spin-down age, could indicate an higher age for the source or, alternatively, a different preferred direction for the nebula extension in the past. Variation in the preferred emission direction may occur due to complex shock interactions. The nebula may be crushed preferentially in certain regions by the reverse shock returning from the progenitor supernova towards the centre of the system at different times for different directions, such as from nearby molecular clouds.

Fig. 10. Evolution of the system varying the energy (age) of the emitting particle. Black point: the current PSR position; blue point: birth PSR position at a characteristic age of 21 kyr; red point: birth PSR position for an hypothetical characteristic age of 60 kyr. The white points (circles) corresponds to the centroids (extensions) of the PWN estimated using the 2D-Gaussian method in different energy bands (see Table 2). For the plot, we use the TS map of the energy range 30–100 GeV as background colour map. The current day major axis for the emission as determined from H.E.S.S. is indicated by a black dashed line.

Whilst in the spherically symmetric advection dominated model a constant nebula size is expected below a cooling break energy (Fig. 9), a contribution to the energy dependent morphology below 1 TeV may be attributable to the pulsar proper motion.

5.2. Particle transport mechanisms

There are several mechanisms which could account for the travelling and cooling of the electron population inside the nebula. Three are the possibilities previously discussed in Aharonian et al. (2006): radiative cooling of electrons while they move away from the pulsar causing energy loss; particle transport mechanisms such as diffusion or advection; and variation in the electron injection spectrum over time.

The energy dependent morphology of the nebula can be reproduced with a multi-zone model for electron injection and evolution. The electron population is injected with the same spectrum at the pulsar and evolved in both time (including cooling) and space, assuming dominant advection with a velocity profile of Eq. (10). This is consistent with both previous models of the PWN (Van Etten & Romani 2011) and recent H.E.S.S. results (H.E.S.S. Collaboration 2019) as well as this analysis.

However, a constant value for β, the index of the radial and time dependence of the velocity profile, appears not to be favoured; at low energies the data in Fig. 9 favours the upper edge (β = 0.5), tending towards the lower edge (β = 0.75) at mid-high energies; implying either a change in velocity profile with time and/or that the description is incomplete. Further effects, such as an additional diffusion component as used by Van Etten & Romani (2011) are not included here and may also need to be incorporated for a fuller treatment of the PWN. Additional diffusive effects on top of the bulk advective motion assumed here are not ruled out.

5.3. The transitional PWN – TeV halo nature of HESS J1825-137

The discovery of extended TeV emission around the Geminga pulsar (Abeysekara et al. 2017a), whose properties are consistent with free particle propagation in the interstellar medium (ISM), suggests the presence of such halo phenomena in other sources (Abeysekara et al. 2017b) as well as its presence also at GeV energies (Di Mauro et al. 2019). Following the approach of Giacinti et al. (2020), we estimated the energy densities for the PWN HESS J1825-137, and compared it with the typical energy density of the ISM, ϵ_ISM = 0.1 eV cm⁻³, in order to determine to a first approximation whether the electrons that are responsible for the γ-ray emission occupy the relatively unperturbed ISM (TeV halo like, ϵ_e ≲ ϵ_ISM ), or if they are still contained in a region energetically and dynamically dominated by the pulsar (PWN like, ϵ_e >  ϵ_ISM).

For the estimation of the electron energy density, we divided the total energy of the electron (W_e) derived with the NAIMA package (see Table 5) by the volume of the nebula. As a first approximation, the nebula has been assumed to be a sphere of a radius of ∼1.35° (see Table 1), which corresponds to a physical radius of 91 pc (2.8 × 10²⁰ cm) and volume V ∼ 0.92 × 10⁶² cm³. Using this basic assumption the mean energy density obtained is ϵ_e, 1 = 0.16 eV cm⁻³. Taking into account the variability of the extent measurement versus energy, where low energetic particles are living in a wider space than high energetic particles that are more concentrated around the pulsar, we computed the gamma-ray intensity weighted mean of the volume using the full nebula spectrum obtaining a volume of V ∼ 0.84 × 10⁶² cm³. For this case the resulting electron energy density is ϵ_e, 2 = 0.17 eV cm⁻³. Both the obtained values for the electron energy density (ϵ_e, 1, ϵ_e, 2) are compatible to the value derived in Giacinti et al. (2020), ϵ_e = 0.25 eV cm⁻³, supporting the transitional scenario, PWN – TeV halo like, for this source. At this stage, high-energy electrons start to escape from the PWN, and propagate into the surrounding supernova remnant, with further escape into the surrounding ISM becoming possible.

6. Conclusion

Thanks to the first energy-dependent analysis of the HESS J1825-137 in the GeV range, we found continued morphological changes and increasing size of this PWN towards lower energies. The PWN extent was measured using two complementary approaches; a 2D-Gaussian template fit and the radial profile method as adopted by H.E.S.S. Collaboration (2019), with compatible results. Not only does the PWN extension continue increasing towards lower energies, indicating a possible turnover only below few GeV, but also the best fit centroid of the emission shifts in a direction opposite to the pulsar proper motion. Furthermore, considering that the change in the centroid position is larger than the change in pulsar position (see Fig. 10), this may be an indication that the system age is somewhat older than that suggested by the 21 kyr characteristic age of PSR J1826-1334 or, alternatively, that the preferred direction for the particle transport and therefore nebula extension has varied over time. The combined Fermi-LAT and H.E.S.S. SED of the nebula could be described by a single electron population model with a break at few hundred GeV. To simultaneously describe the SED and energy dependent morphology of the nebula, we used a multi-zone modelling approach with burst-like injection and an advective velocity profile is found to be consistent with the data. The order β of the radial and temporal velocity dependence must, however, vary within the range [0.5,0.75]; a constant velocity profile is not compatible. The estimated values of the electron energy density support the transitional scenario, PWN - TeV halo like, for this source.

HESS J1825-137 is one of the most γ-ray luminous and TeV efficient PWN known (H.E.S.S. Collaboration 2018), enabling rich and detailed analyses to be performed. Future studies and modelling of this source as well as similar systems will enable further insights into pulsar wind nebula formation and evolution to be gained.

Acknowledgments

The effort of the LAT-team Calibration & Analysis Working Group to develop Pass 8, as well as the LAT Galactic Working Group for the support given to this project, are gratefully acknowledged. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Etudes Spatiales in France. This work is performed in part under DOE Contract DE-AC02-76SF00515.

References

Appendix A: Visual comparison of the 2D-Gaussian and radial profile extent estimates

In this section we report the excess maps with the comparison between the extent obtained using the 2D-Gaussian (see Sect. 3.2.1) and the radial profile method (see Sect. 3.2.2).

We also report in Table A.1 the fit parameters for the polynomial parameterisation in Eq. (3). In particular, r₀ is the distance at which it approaches the constant value c, and a provides the overall normalisation.

Fig. A.1. Left: excess maps (in sigma units), in celestial coordinates, of the region around HESS J1825-137. The magenta points and the white dashed circles represent the 2D-Gaussian centroid and extension respectively, while the dashed semicircumference shows the extension obtained with the radial profile method considering only the southern hemisphere. For comparison, the region used to extract the radial profile (as used by H.E.S.S.) is overlaid in white. The preferred emission direction (major axis) as found by H.E.S.S., along which the extent is evaluated, is indicated by the black dashed line. Right: radial profile of the excess counts fit with Eq. (3) beyond the peak emission. The characteristic R(1/e) size of the nebula is indicated by a white (black) dashed line in the left (right) hand plot. The plots are related to the energy bands: 1–3 GeV (top) and 3–10 GeV (bottom).

Fig. A.2. Same as Fig. A.1 for the energy bands: 10–32 GeV (top), 32–100 GeV (middle) and 100 GeV–1 TeV (bottom).

All Tables

All Figures

	Fig. 1. Sky map, in Galactic coordinates and Mollweide projection, showing the sources in the TevCat catalogue (http://tevcat.uchicago.edu, version April 2019) classified by their most likely association. All the 3FHL sources (Ajello et al. 2017) are also plotted, with grey points, for a comparison.
In the text

Fig. 2. TS map (in sigma units), in celestial coordinates, of the region around the PWN HESS J1825-137 in the energy range between 1 GeV and 1 TeV. The red circle and star indicate the 2D-Gaussian extension and centroid fit obtained in this work. The red dashed circles mark the uncertainty on the 2D-Gaussian extension. The green and white circles correspond to the extension obtained in the FGES catalogue (Ackermann et al. 2017a) and in Grondin et al. (2011), respectively. The black point indicates the position of PSR J1826-1334, the pulsar which is believed to power the nebula. The 4FGL sources, as well as the three candidate sources added in the model, are represented with light grey points.

In the text

Fig. 3. TS maps (in sigma units), in celestial coordinates, of the region around HESS J1825-137 for the energy bands: 1–3 GeV (top left), 3–10 GeV (top centre), 10–32 GeV (top right), 32–100 GeV (bottom left) and 100 GeV–1 TeV (bottom right). The TS maps are smoothed with a Gaussian of radius 0.1°. The white circles represent the extension (solid line) and its statistical uncertainty (dashed lines) determined in the respective energy band. For comparison, the resulting extension obtained for the entire energy range (1 GeV–1 TeV) is overlaid with a red line. All extensions shown correspond to the 68% containment radius. In the 100 GeV–1 TeV (bottom right) plot, H.E.S.S. significance contours at 5, 10, and 15σ, for energies below 1 TeV (H.E.S.S. Collaboration 2019), are shown with light-blue lines for comparison. The H.E.S.S. contour also includes the excess of the nearby LS 5039 source (small circular excess at the southern boarder of the PWN).

In the text

Fig. 4. Left: TS map (in counts) of the PWN HESS J1825-137 in the energy range 1 GeV–1 TeV with Fermi-LAT, with the region used to extract the radial profile (as used by H.E.S.S.) overlaid in white. The preferred emission direction (“major axis”) as found by H.E.S.S., along which the extent is evaluated, is indicated by the black dashed line. The position of the pulsar (black) and best-fit Fermi-LAT centre of a 2D Gaussian (magenta) are also indicated. Right: radial profile of the excess counts fit with Eq. (3) beyond the peak emission. The characteristic R(1/e) size of the nebula is indicated by a white (black) dashed line in the left (right) hand plot.

In the text

	Fig. 5. Nebula extent as a function of energy showing results from H.E.S.S. Collaboration (2019) and this analysis.
In the text

Fig. 6. Combined spectra of the PWN HESS J1825-137 with the spectral measurements obtained in this work (red points) using 11.6 years of Fermi-LAT data from 1 GeV to 1 TeV and the H.E.S.S. results for the 100 GeV–90 TeV energy range (black points). The Fermi-LAT flux points were obtained doing a spectral analysis in the various energy bands with their relative spatial model. The combined SED has been fitted with both a LogParabola (blue line) and a Broken PL (green line). The vertical line corresponds to the energy break Eb of the Broken PL model. The bottom panel shows the normalised residual between the data and the LogParabola model.

In the text

	Fig. 7. Results of the fitted IC NAIMA SED to the Fermi-LAT data and H.E.S.S. data, compared with the measured Fermi-LAT (red circle) and H.E.S.S. (black diamond) data points. A synchrotron component is computed from the electron distribution supposing 4, 5, 6, 7, and 8 μG, and compared with Suzaku data (blue square). The electron distribution parameters are presented in Table 5.
In the text

	Fig. 8. Total SED of the nebula, combining X-ray and γ-ray data, is described by the summation of 20 zones of particles of different ages using the GAMERA modelling package Hahn (2016). The total SED of each zone is shown by a coloured line, from yellow for the youngest and smallest zone, through to blue for the oldest and largest. Model parameters are given in Table 6.
In the text

	Fig. 9. Radial extent of the nebula as a function of energy; the results of the model using GAMERA is indicated by grey shaded region for a compatible range of β, the index of the radially and time dependency of the velocity profile, between 0.5 (upper edge) and 0.75 (lower edge).
In the text

Fig. 10. Evolution of the system varying the energy (age) of the emitting particle. Black point: the current PSR position; blue point: birth PSR position at a characteristic age of 21 kyr; red point: birth PSR position for an hypothetical characteristic age of 60 kyr. The white points (circles) corresponds to the centroids (extensions) of the PWN estimated using the 2D-Gaussian method in different energy bands (see Table 2). For the plot, we use the TS map of the energy range 30–100 GeV as background colour map. The current day major axis for the emission as determined from H.E.S.S. is indicated by a black dashed line.

In the text

Fig. A.1. Left: excess maps (in sigma units), in celestial coordinates, of the region around HESS J1825-137. The magenta points and the white dashed circles represent the 2D-Gaussian centroid and extension respectively, while the dashed semicircumference shows the extension obtained with the radial profile method considering only the southern hemisphere. For comparison, the region used to extract the radial profile (as used by H.E.S.S.) is overlaid in white. The preferred emission direction (major axis) as found by H.E.S.S., along which the extent is evaluated, is indicated by the black dashed line. Right: radial profile of the excess counts fit with Eq. (3) beyond the peak emission. The characteristic R(1/e) size of the nebula is indicated by a white (black) dashed line in the left (right) hand plot. The plots are related to the energy bands: 1–3 GeV (top) and 3–10 GeV (bottom).

In the text

	Fig. A.2. Same as Fig. A.1 for the energy bands: 10–32 GeV (top), 32–100 GeV (middle) and 100 GeV–1 TeV (bottom).
In the text