Broadband X-ray Spectroscopy of the Pulsar Wind Nebula in HESS J1640-465

TL;DR

This study analyzes the X-ray properties of the pulsar wind nebula in HESS J1640-465 using Chandra and NuSTAR data, revealing spectral softening, high magnetic fields, and potential PeVatron status.

Contribution

It provides updated X-ray measurements, accounts for dust scattering effects, and employs a model to suggest the nebula's potential as a PeVatron candidate.

Findings

Spectral softening observed at higher energies

High magnetic field and maximum electron energy up to PeV

HESS J1640-465 may be a PeVatron candidate

Abstract

We present updated measurements of the X-ray properties of the pulsar wind nebula associated with the TeV $\gamma$-ray source HESS J1640-465 derived from Chandra & NuSTAR data. We report a high $N_{H}$ value along line of sight, consistent with previous work, which led us to incorporate effects of dust scattering in our spectral analysis. Due to uncertainties in the dust scattering, we report a range of values for the PWN properties (photon index and un-absorbed flux). In addition, we fit the broadband spectrum of this source and found evidence for spectral softening and decreasing unabsorbed flux as we go to higher photon energies. We then used a one zone time dependent evolutionary model to reproduce the dynamical and multi-wavelength spectral properties of our source. Our model suggests a short spin-down time scale, a relatively higher than average magnetized pulsar wind, a strong PWN magnetic field and maximum electron energy up to PeV, suggesting HESS J1640-465 could be a PeVatron candidate.

Tables

Table 1. Table 6: Chandra and NuSTAR’s PWN joint fit results using different pulsar properties

Parameter	Low Boundary		Best Fit		High Boundary
Joint fit between Chandra and NuSTAR’s PWN (3$-$10 keV)
	(1)	(2)	(3)	(4)	(5)	(6)
${\Gamma }_{PWN}$	$1.9\pm 0.1$	$1.7\pm 0.1$	$2.0\pm 0.1$	$1.6\pm 0.1$	$2.1\pm 0.1$	$1.6\pm 0.1$
$F_{PWN}$	$\left(8.5\pm 0.2\right)\times {10}^{-13}$	$\left(8.2\pm 0.2\right)\times {10}^{-13}$	$\left(9.5\pm 0.2\right)\times {10}^{-13}$	$\left(8.2\pm 0.2\right)\times {10}^{-13}$	$\left(10.9\pm 0.3\right)\times {10}^{-13}$	$\left(9.0\pm 0.3\right)\times {10}^{-13}$
${\chi }_{\nu }^2$(DOF)	1.14 (554)	1.14 (554)	1.14 (554)	1.12 (554)	1.16 (554)	1.14 (554)
NuSTAR’s PWN (10$-$30 keV)
	(7)	(8)	(9)	(10)	(11)	(12)
${\Gamma }_{PWN}$	${8.8}_{-2.0}^{+0.0}$	$3.2\pm 0.5$	$3.5\pm 0.5$	$1.8\pm 0.2$	$2.1\pm 0.3$	$1.6\pm 0.2$
$F_{PWN}$	$\left(3.14\pm 0.2\right)\times {10}^{-14}$	$\left(2.0\pm 0.4\right)\times {10}^{-13}$	$\left(2.0\pm 0.4\right)\times {10}^{-13}$	$\left(5.8\pm 0.5\right)\times {10}^{-13}$	$\left(4.9\pm 0.5\right)\times {10}^{-13}$	$\left(7.9\pm 0.5\right)\times {10}^{-13}$
${\chi }_{\nu }^2$(DOF)	1.4 (765)	1.09 (765)	1.09 (765)	1.07 (765)	1.07 (765)	1.08 (765)

Equations

$$t_{age}$$

$${\dot{E}}_0$$

$${\chi }^2=\sum _{i=1}^{i=35}{\left(\frac{{\mathcal{D}}_i-{\mathcal{M}}_i}{{\sigma }_i}\right)}^2$$

$${\Gamma }_{PSR}=2.08-0.029{\dot{E}}_{40}^{-1/2},$$

$$P_0=P{\left(1+\frac{t_{\text{age }}}{{\tau }_{\text{sd }}}\right)}^{-\frac{1}{p-1}}$$

$$\dot{P}(t)$$

$${\dot{P}}_0$$

$$B_{s0}\equiv 3.2\times 10^{19}\sqrt{P_0{\dot{P}}_0}\approx 1.8\times 10^{13}\mathrm{G}$$

$$B_{\mathrm{LC}}={\left(\frac{3I8{\pi }^4\dot{P}}{c^3{P}^5}\right)}^{1/2}\approx 2.94\times 10^8{\left(\dot{P}{P}^{-5}\right)}^{1/2}\mathrm{G}$$

Full text

supernova remnants, neutron stars, pulsars, pulsar wind nebula

††journal: ApJ

\savesymbol

tablenum \restoresymbolSIXtablenum

1 Introduction

Since neutron stars are rapidly rotating and highly energetic, they power a highly relativistic electron positron wind which in turn interacts with the surrounding environment as it expands outwards creating a pulsar wind nebula (PWN). PWNe emit electromagnetic radiations through synchrotron emission in a wide range of energy going from radio up to X-rays due to the presence of highly relativistic electrons. In addition, the relativistic electrons can cause highly energetic emissions in the $\gamma$-ray part of the spectrum through inverse Compton scattering of ambient photons.

HESS J1640-465 was first discovered by the HESS collaboration in 2006 (Aharonian et al., 2006). It was found to be spatially coincident with G338.3-0.0, a shell like supernova remnant (SNR) previously identified through radio observations to have an 8*′* diameter (Whiteoak & Green, 1996). The source was analyzed by Landi et al. (2006) using data obtained by Swift X-ray telescope and a possible connection with pulsar wind nebula was suggested. Later, Funk et al. (2007) reported the detection of a highly absorbed X-ray point source with extended emission in the field of view of HESS J1640-465 using XMM-Newton observations. However, due to low angular resolution, neither a pulsar nor the associated PWN were spatially resolved. Lemiere et al. (2009) then managed to spatially resolve the pulsar and its associated PWN using observations from Chandra. No X-ray emission has been detected from the SNR, likely due to its low temperature and the high value of interstellar absorption along the line of sight ($N_{H}\sim 1.4\times 10^{23}{\rm cm}^{-2}$) (Gotthelf et al., 2014). Using the 21 cm HI absorption line, a distance of 8$-$13 kpc was derived for the source (Lemiere et al., 2009) which makes HESS J1640-465 the most luminious TeV source in our galaxy. Castelletti et al. (2011) derived upper limits on the radio flux with no detection of a radio counterpart and as of now, there is no $\gamma$-ray or radio pulsation from the source. Subsequently, Gotthelf et al. (2014) reported a discovery of X-ray pulsations originating from a 206 ms pulsar (PSR J1640-4631) using data from NuSTAR as part of a survey covering the Norma Arm region of the galactic plane

(Harrison et al., 2013). The pulsar was found to be spatially coincident with HESS J1640-465. It has a rapid spin-down rate of $\dot{P}=9.8\times 10^{-13}$, a spin-down luminosity of $\dot{E}=4.4\times 10^{36}$ erg $s^{-1}$ and a characteristic age of $\tau_{c}$ = 3350 yr. Archibald et al. (2016) then performed a 2.3 year phase-connected timing analysis on this source and calculated its braking index to be p = $3.15\pm{0.03}$ making it to be the first pulsar with a braking index $>$ 3.

Slane et al. (2010) reported the discovery of the Fermi source 1FGL J1640.8-4631 that is spatially coincident with HESS J1640-465. They concluded that the observed TeV emission arises from the PWN through inverse Compton Scattering. On the contrary, Abramowski et al. (2014) proposed that the observed emission would arise due to proton-proton interactions in the northern part of the spatially coincident SNR G338.4+0.1. Mares et al. (2021) conducted morphological analysis of the source using 8 years of Fermi-LAT data. This showed that HESS J1640-465 is largely extended in the GeV energy band. In addition to that, they modelled the broadband emission coming out of the source assuming two different scenarios, one where the emission is originating from the PWN linked to the pulsar PSR J1640-4631 and the second where the accelerated protons or electrons inside the SNR shock are responsible from the observed emission. For the first scenario, they managed to set an upper limit on the PWN’s magnetic field as well as put some constrains on the spin-down power of the pulsar. As for the second scenario, the ratio of electrons to protons in the SNR shock must be high. And though the broadband emission could be reproduced under this scenario, they make the remark that it is less likely because of the large ratio between the radio upper limit and the $\gamma$-ray emission.

In §2, we outline the Chandra and NuSTAR observations used in this paper and describe our analysis steps and methods. In §3, we report our X-ray spectral analysis results. And then in §4, we apply an evolutionary model to multi-wavelength observations of HESS J1640-465 and then discuss the implications of our spectral energy distribution modelling in §5. Finally, we summarize our findings and conclusions in §6.

2 Observations and Data Reduction

2.1 NuSTAR Observations and Data Reduction

The Nuclear Spectroscopic Telescope Array (NuSTAR) is an X-ray facility operating between 3$-$79 keV. It has two co-aligned telescopes with detectors on their focal plane modules named FPMA and FPMB. Each telescope provides a $12^{\prime}\times 12^{\prime}$ field of view with a point spread function (PSF) that has half power diameter (HPD) of $58^{\prime\prime}$ and full width half maximum (FWHM) of $18^{\prime\prime}$. HESS J1640-465 was observed by NuSTAR over the years 2013, 2014 and 2015 for the purpose of following the timing evolution of the pulsar. This amounted to a total exposure time of $\sim$ 759 ks accumulated in 24 separate pointings. A log of these observations is presented in Table 1.

Using NuSTARDAS version 2.0.0 and HEASOFT version 6.28 available from the NASA High Energy Astrophysical Science Archive Research Center (HEASARC), we followed standard pipeline processing to generate cleaned and calibrated “Level 2 Data products” using the “nupipeline” script. Once this stage was completed, we used the “nuproducts” script with the option extended=yes to generate the associated response matrix (RMF) and ancillary response files (ARF) for an extended source.

We extracted individual spectra from the PWN using an elliptical region of diameter $3.4^{\prime}\times 2.6^{\prime}$ to ensure consistency with Gotthelf et al. (2014). We also used an offset circular region of radius $76.6^{\prime\prime}$ to account for background emissions. Figure 1 shows the source and background regions used for extraction. Spectra from each focal module were combined as well as their respective response files/matrices using the addascaspec routine in XSPEC version 12.11.1 (Arnaud, 1996) and subsequent fitting was also performed in XSPEC.

2.2 Chandra Observations and Data Reduction

HESS J1640-465 has been observed using Chandra’s Advanced CCD Imaging Spectrometer (ACIS)(Garmire et al., 2003). Each ACIS chip detects photons in the energy range of 0.2–10 keV with a field of view 8*′.3 $\times$ 8′.3 and imaging resolution of $\sim$ 1′′*. For our analysis, we analyze ACIS-I observations taken in 2015 with a total observation time of $\sim$ 150 ks, as detailed in Table LABEL:tab:Chandra_log. All Observations were in timed exposure and very faint mode. Obs.Id 17582 was removed from our analysis due to its short exposure time $(\sim$9 ks). This data-set comprises $\sim$ 3x more exposure time than that analyzed in previous work by Gotthelf et al. (2014).

Each observation was reprocessed individually using the $chandra\_repro$ package of the Chandra Interactive Analysis of Observations (CIAO) software suite version 4.13 (Fruscione et al., 2006), which resulted in a reprocessed $event\_2$ file for each observation. The reprocessed event files were then merged using the $reproject\_obs$ routine to create a deeper combined image of the source. Figure 2 shows the resulting image of the source after combining the four data sets. This is a zoomed in, smoothed image of the source highlighting the region of the pulsar and its extended diffuse emission associated with the PWN. We extracted spectra from the individual observations logged in Table 2 along with their response matrices using CIAO’s $specextract$ tool. The spectra were then co-added using the $combine\_spectra$ routine in Sherpa version 4.13.0 (Freeman et al., 2001).

For the pulsar, we used a circular region of radius 2.8*′′* to extract the spectrum and its associated response files for a Point-like source using CIAO version 4.13 (Fruscione et al., 2006) with CALDB version 4.9.4. For the whole source, we used the same elliptical region we used earlier for NuSTAR to account for scattering along line of sight as discussed in §3.3. By using this large source region, we ensure all the scattered emission is taken into consideration. We extracted the spectrum and response files for an extended object also using CIAO version 4.13 (Fruscione et al., 2006). We used a rectangular background region chosen to be far from the source region to minimize possible contamination due to extended emission coming from the PWN. Figure 3 shows the chosen extraction region for the whole source, pulsar only and their associated background region. The background subtracted total of 405 photons across all observations were detected from the pulsar. We investigated the possibility of pileup using the CIAO routine $pileup\_map$, but we found no evidence of pileup in the data.

3 Spectral Analysis

In NuSTAR observations, the pulsar and PWN are not spatially resolved. Only Chandra has sufficient angular resolution to spatially resolve the emission coming out of the pulsar from that of the PWN. Therefore, to constrain the PWN spectrum in the data obtained by NuSTAR, we will do the following:

• •

Quantify the flux systematic offsets (Cross calibration constant) between Chandra and NuSTAR as discussed in §3.1

• •

Analyze Chandra’s pulsar data alone to account for scattering along line of sight as discussed in §3.3 and §3.3.1

• •

Use the un-scattered pulsar flux along with the cross normalization factor to jointly fit Chandra and NuSTAR spectra to obtain accurate measurement of the PWN spectrum as discussed in §3.4

For spectral fitting, we used XSPEC version 12.11.1 (Arnaud, 1996) and we accounted for the column density using the Tuebingen-Bolder ISM absorption model, implemented in Xspec as the TBabs with solar abundances set to wilms (Wilms et al., 2000). The fit statistic was set to $\chi^{2}$ in all of our fits and XSPEC’s Markov Chain Monte Carlo (MCMC) method was used to obtain the uncertainties in our model parameters and all the reported uncertainties in this paper are reported at the $90\%$ confidence level.

3.1 Cross Normalization/Calibration

Since we jointly fit data from two different observatories (Chandra and NuSTAR), cross calibration is important. To find the cross normalization between the two missions for our source, we extracted the Chandra spectra for the whole source using the same elliptical region used to extract the NuSTAR spectra. In doing that, emission from both the pulsar and PWN is included in our selected region. Figure 3 shows the extraction region used for the whole source used to find the cross normalization factor.

For NuSTAR, we initially fitted spectra from FPMA and FPMB independently and got consistent results. Therefore throughout this paper, we will be fitting them simultaneously.

We first jointly fit the NuSTAR and Chandra spectra in their overlapping 3$-$10 keV energy range using an absorbed power-law model, $tbabs\times pegpwrlw\times constant$ (in XSPEC) with photon index $(\Gamma)$ and hydrogen column density ($N_{H}$) linked across the two data-sets. This constant was fixed to 1 for Chandra and allowed to vary for NuSTAR that resulted in a value of 0.75 for the NuSTAR cross calibration constant, which is comparable to the value obtained by Madsen et al. (2017). This value will be subsequently used for all our joint fits in this paper.

3.2 Chandra Data

For the pulsar, we used the circular extraction region shown in Figure 3. We restricted our fit to be between 0.2$-$10 keV and grouped the pulsar spectrum to have 25 counts per bin and used an absorbed power law ($tbabs\times pegpwrlw$) to fit the data (Yang et al., 2016).

3.2.1 Results

For the pulsar, we obtained a best fit $N_{H}=(1.4^{+0.5}_{-0.4})\times 10^{23}$ ${\rm cm}^{-2}$, $\Gamma=0.6^{+0.6}_{-0.5}$ and an absorbed flux at 2$-$10 keV, F = $2.5^{+0.4}_{-0.3}\times 10^{-13}$ erg ${\rm cm}^{-2}$ $s^{-1}$. For this fit, the reduced $\chi^{2}_{\nu}$ = 0.98 for 18 degrees of freedom and the absorbed flux at 2$-$10 keV is reported for comparison with previous work done by Gotthelf et al. (2014). The resulting spectral fit is shown in Figure 4 and the best fit parameters are tabulated in Table 3 in which we compare our results with earlier work on this source done by Gotthelf et al. (2014). While our results generally agree with Gotthelf et al. (2014), we set a tighter constraint on the pulsar properties due to the larger exposure time analyzed in this paper.

3.3 Scattering Effects

As X-ray photons travel through the interstellar medium (ISM), they can be absorbed or scattered, which affects the measured spectrum within a particular source region around astrophysical sources. The Scattering of X-ray photons by dust along the line of sight forms x-ray halo around these objects. As a result, the spectrum of photons detected in a small region around the object is both lower and harder than the intrinsic spectrum of the source (Smith et al., 2016). Since our source is highly absorbed ($N_{H}$ = $1.4\times 10^{23}$ ${\rm cm}^{-2}$), and consistent with earlier studies (e.g, Slane et al. 2010 & Gotthelf et al. 2014 ), we need to take scattering into consideration in our spectral analysis. We incorporate these effects using the $xscat$ model in XSPEC (Smith et al., 2016). This model calculates the un-scattered intrinsic pulsar spectrum taking into account the fitted photon energy range, model for dust, extraction region and distance between the dust and the source. Since we have no prior information on either the type of dust or its location along the line of sight, we decided to explore the derived intrinsic values of pulsar photon indices and un-absorbed fluxes for different dust models over assumed range of $N_{H}$ values derived from our analysis of the pulsar data.

For dust model, we used two of $xscat^{\prime}s$ models; first is the MRN model based on the work of Mathis et al. (1977) and the second is the ZDABAS model based on the work of Zubko et al. (2004), and we allow the relative location of dust grains along line of sight from the source to vary between 0 (the dust grains are at the source) and 1 (the dust grains are at the observer).

3.3.1 Results of Scattering Effects

Table 4 shows our results when we include the $xscat$ model in our spectral fits. For the same value of $N_{H}$, the MRN and ZDABAS dust models resulted in significant effects on the pulsar spectrum depending on how far the dust is from our source. It is worth mentioning that the four italicized values without error bars in Table 4 in the MRN dust model resulted in a bad fit with $\chi_{\nu}^{2}$ $>$ 2, so we don’t take them into consideration. In addition, at $N_{H}=1.9\times 10^{23}{\rm cm}$, our fits result in higher $\chi_{\nu}^{2}$ values than at lower $N_{H}$ as detailed in Table 4, which suggests that a lower column density value is favoured. However, we decided to include all results in our PWN analysis for completeness.

For both dust models, we observe a similar trend, but systematic offsets in the two parameters we are observing. To be more specific, for the same value of $N_{H}$ at the same distance, the MRN model always results in higher flux value and softer (larger) photon index than the ZDABAS model.

For different values of $N_{H}$ at the same location of dust along line of sight, increase in $N_{H}$ results in an increase in flux and photon index. This trend is seen in both models as shown in Figure 5.

To accommodate as large of an uncertainty range as we can and propagate that to the NuSTAR data-set to obtain a more accurate PWN spectrum that will be used in the modelling later, we therefore selected two points for each $N_{H}$ value; one representing the lowest flux and hardest spectrum (lowest $\Gamma_{PSR}$) and the other representing the highest flux and softest spectrum (highest $\Gamma_{PSR}$). Table 5 shows our final derived values for the pulsar properties after taking dust scattering into consideration. In total, we have 6 different combinations of parameters that will be used in fitting the NuSTAR data in the next section.

3.4 NuSTAR Data

For our source, we used the extraction region shown in Figure 1 and grouped the spectrum to have at least 5 $\sigma$ per spectral bin. Though NuSTAR can detect photons with energies up to 79 keV, emission from our source is only brighter than the background below 30 keV, and therefore we only fit the measured spectrum below this energy as shown in the right panel of Figure 6.

As previously mentioned, NuSTAR lacks sufficient angular resolution to spatially resolve the emission coming from the pulsar and the PWN. To overcome this, we will fit the composite spectrum with two absorbed power law components, one to account for emission from the pulsar and the other to account for emission from the PWN. The properties of the pulsar component will be fixed to the values in Table 5.

3.4.1 Results

Since Chandra and NuSTAR overlap in the 3$-$10 keV energy range, we jointly fit them together. For that, we used the entirety of the Chandra selection region shown previously in Figure 1 and used a $tbabs\times(pegpwrlw+pegpwrlw)\times constant$ model with the constant value set to 1 for the Chandra data and 0.75 for both FPMA and FPMB to account for the cross-normalization issues between the two telescopes as explained in §3.1. The parameters of the first power law $-$representing the derived pulsar properties$-$ are shown in Table 5 and the joint fit results representing the PWN emission are shown in Table 6 and two representative joint fits are shown in Figure 6. All the fits yielded good results with an average $\chi^{2}_{\nu}$ of 1.14. It is worth mentioning that fit number (7) resulted in a much worse $\chi_{\nu}$ than all the other, therefore we decided to exclude it from our consideration. This is the result of having a very hard pulsar photon index which means most of the emission is dominated by the pulsar with little to almost no contribution from the PWN.

In the 3$-$10 keV energy range and for the same value of $N_{H}$, the PWN gets harder (decrease in $\Gamma_{PWN}$) and emits less flux as the pulsar gets softer (increase in $\Gamma_{PSR}$). For the six different combinations of the pulsar parameters we used, we get a value of $\Gamma_{PWN}$ ranging between 1.6 to 2.1 and a value of un-absorbed flux ranging from $8.2\times 10^{-13}$ to $10.9\times 10^{-13}$ erg ${\rm cm}^{-2}$ $s^{-1}$.

In the 10$-$30 keV energy range, we jointly fit FPMA and FPMB only, excluding the Chandra data since it doesn’t extend beyond 10 keV. The assumed pulsar emission between 10$-$30 keV is a power-law extrapolation of the pulsar emission between 3$-$10 keV, as measured by Chandra after accounting for cross calibration offset between these two telescopes in their overlapping energy range as explained in §3.1

In this higher energy range, and for the same value of $N_{H}$, the PWN gets harder (decrease in $\Gamma_{PWN}$) and emits more flux as the pulsar gets softer (increase in $\Gamma_{PSR}$). For the six different combinations of the pulsar parameters we used, we get a value of $\Gamma_{PWN}$ ranging between 1.6 to 3.5 and a value of un-absorbed flux ranging from $2.0\times 10^{-13}$ to $7.9\times 10^{-13}$ erg ${\rm cm}^{-2}$ $s^{-1}$.

In addition, for 2 out of the 3 derived pulsar parameters, we generally observe a decrease in the PWN flux and spectral softening (increase in $\Gamma_{PWN}$) as we go from the 3$-$10 keV part of the spectrum to the 10$-$30 keV. However, at an $N_{H}$ value of $1.9\times 10^{23}$ ${\rm cm}^{-2}$, the $\Gamma_{PWN}$ remains unchanged across energy bands while the flux goes down, similar to what we observed at the other two $N_{H}$ values. This decrease in flux across energy bands was observed before by Hattori et al. (2020) in their analysis of PWN G21.5-0.9.

To summarize our spectral analysis, in the 3$-$10 keV energy band, we see a photon index ($\Gamma_{PWN}$) ranging between ($1.6\pm{0.1}$) and ($2.1\pm{0.1}$) and flux values ranging from $(8.2\pm{0.2})\times 10^{-13}$ to $(10.9\pm{0.3})\times 10^{-13}$ erg ${\rm cm}^{-2}$ $s^{-1}$. For the 10$-$30 keV energy range, we see a photon index ($\Gamma_{PWN}$) ranging between ($1.6\pm{0.2}$) and ($3.5\pm{0.5}$) and flux values ranging from $(2.0\pm{0.4})\times 10^{-13}$ to $(7.9\pm{0.5})\times 10^{-13}$ erg ${\rm cm}^{-2}$ $s^{-1}$.

Therefore, the X-ray part of our multi-wavelength modelling will be divided into two energy ranges; soft X-ray range (3$-$10 keV) from the Chandra & NuSTAR joint fit and hard X-ray range which will incorporate the the 10$-$30 keV results from NuSTAR. The range of X-ray spectral parameters used when modelling the data is provided in Table 7.

4 PWN Modelling

PWNe are often studied to infer the properties of their progenitor supernova, associated neutron star along with its pulsar wind since these properties are hard to be measured directly. Currently, this is best achieved by using an evolutionary model of a spherical PWN inside a SNR to fit the observed dynamical and radiative properties of these PWNe (see Gelfand 2017 for a recent review of such models).

Our work is based on the model described by Gelfand et al. (2009) and was previously implemented to study other PWNe systems such as the Eel (Burgess et al., 2022), Kes 75 (Straal et al., 2023; Gotthelf et al., 2021; Gelfand et al., 2014), G21.5-0.9 (Hattori et al., 2020), G54.1+0.3 (Gelfand et al., 2015) and HESS J1640-465 (Gotthelf et al., 2014; Mares et al., 2021).

A Metropolis Markov Chain Monte Carlo (MCMC) algorithm was implemented to identify the best fit values of our 15 model input parameters $\Theta$ outlined in Table 8 that best reproduce the observed properties of our source (see §3.2 of Gelfand et al. 2015 for a full description). In addition to requiring the X-ray spectrum predicted by the model to fall within the range of X-ray spectral parameters listed in Table 7, this includes the radius of the associated SNR G338.3-0.0, $\theta_{snr}$ = $4.45^{\prime}\pm 0.5^{\prime}$ (Shaver & Goss, 1970), the radius of the PWN, $\theta_{pwn}$ = $3.3^{\prime}\pm 0.2^{\prime}$ (Lemiere et al., 2009), the upper limit on the 660 MHz radio flux density of the PWN, $S_{660}<690\pm 300$ mJy (Castelletti et al., 2011) and updated $\gamma$-ray fluxes listed in Table 6 of Mares et al. (2021).

The true age $t_{age}$ of the system, initial spin-down luminosity $\dot{E}_{0}$ and and spin-down timescale $\tau_{sd}$ of the associated pulsar J1640-4631 are calculated by the model using its current spin-down luminosity of $\dot{E}=4.4\times 10^{36}$ ergs/s and its characteristic age $t_{ch}$ = 3350 years for an assumed constant braking index of p = 3.15 as measured by Archibald et al. (2016) :

[TABLE]

The MCMC algorithm makes use of the maximum likelihood estimation method described in detail in Hattori et al. 2020. Table 8 outlines the combination of our 15 model input parameters $\Theta$ that yielded the biggest likelihood of our MCMC runs or equivalently, the lowest $\chi^{2}$, and Table 9 shows the output dynamical properties of our model for this source. In addition, we also include previous modelling results for this source that were done by Mares et al. (2021) using the same model in both tables for comparison purposes and to show how our updated measurement of the X-ray spectrum can lead to different predictions by our evolutionary model. Figure 7 shows our well fitted SED for the observed properties of this source. Statistically, $\chi^{2}=29.9$ for 22 degrees of freedom resulting in $\chi_{\nu}^{2}=1.36$ where we define $\chi^{2}$ as:

[TABLE]

where $\mathcal{D}_{i}$ is the observed property, $\mathcal{M}_{i}$ is our model-predicated value for the observed property and $\sigma_{i}$ is the error on observed quantity $\mathcal{D}_{i}$.

5 Discussion and Conclusions

For the supernova properties ($E_{sn},M_{ej}$), our values are in agreement with the values reported by Mares et al. (2021), which in turn are considered to be typical values for supernova explosions (e.g Baade & Zwicky 1934).

For PSR J1640$-$4631, the pulsar associated with HESS J1640$-$465, we got a spin-down timescale of $\tau_{sd}\approx$ 5 years. While this value is in agreement with the value reported by Mares et al. (2021), this is an extremely short timescale in comparison with a typical value of $\tau_{sd}\sim 1000$ years reported for other pulsars as shown in Table 10. (also see Tanaka & Takahara 2011, Torres et al. 2014, Gelfand et al. 2015 & Hattori et al. 2020).

For young, highly energetic ( $\dot{E}$ >  $4\times 10^{36}$ erg $s^{-1}$ ) pulsars, Gotthelf (2003) observed a correlation between the X-ray spectra of pulsars and their $\dot{E}$ following this equation:

[TABLE]

where $\dot{E}_{40}$ is the spin-down energy in units of $10^{40}$ erg $s^{-1}$. By plugging in our pulsar’s $\dot{E}$ = $4.4\times 10^{36}$ erg $s^{-1}$, we get a $\Gamma_{PSR}=0.7$. Therefore, looking back at Figure 5, this observed correlation favors a low $N_{H}$ value of $1.0\times 10^{23}{cm}^{-2}$ and a nearby location for the dust grain to our source.

Using these values along with the current pulsar period of P $\approx$ 206 ms (Gotthelf et al., 2014), we can calculate the initial period $P_{0}$ and initial period derivative $\dot{P}_{0}$ using the following equation (e.g., Pacini & Salvati 1973 & Gaensler & Slane 2006 and the references therein):

[TABLE]

[TABLE]

Plugging our values in the above two equations, we get an initial period value of $P_{0}\approx 10.1$ ms, and an initial period derivative value of $\dot{P}_{0}\approx 3.13\times 10^{-11}\frac{\mathrm{s}}{\mathrm{s}}$. Using these values for $P_{0}$ & $\dot{P}_{0}$, we get an initial spin-down inferred surface dipole magnetic field strength of:

[TABLE]

only $\sim 30\%$ larger than the current value of $B_{s}\approx 1.4\times 10^{13}\mathrm{G}$, as reported by Gotthelf et al. (2014). It is also worth mentioning that while the predicted initial spin down luminosity $\dot{E}_{0}=1.2\times 10^{42}$ erg $s^{-1}$ of this pulsar is significantly larger than the values derived for other systems (e.g., Tanaka & Takahara 2011; Torres et al. 2014; Gelfand et al. 2015; Hattori et al. 2020).

As for the magnetization of the pulsar wind, we get a value of $\eta_{B}\approx 0.1$ which is extremely higher than a typical value of $\eta_{B}\approx 2.2\times(10^{-2}$– $10^{-3})$ inferred from other pulsars (e.g., Tanaka & Takahara 2011; Torres et al. 2014; Gelfand et al. 2015; Hattori et al. 2020). This value, however, is comparable with the also high wind magnetization value predicted for Kes 75 (Gotthelf et al., 2021) as shown in Table 10. It is noteworthy that both of them have an inferred dipole surface magnetic field strength $B_{ns}\gtrsim 10^{13}$ G associated with their pulsars. To put this into perspective, Table 10 shows a comparison between the pulsar parameters of PSR J1640-4631 and other systems. There seems to be a possible relation between the dipole magnetic field strength ($B_{ns}$) and the resulting wind magnetization ($\eta_{B}$); the higher the $B_{ns}$, the higher the resulting $\eta_{B}$. However, the sample size isn’t big enough to draw rigorous conclusions from this, and therefore more systems need to be modelled in order to investigate how real this correlation is. It is important to mention that all these properties were inferred using the same model described by Gelfand et al. (2009) to minimize systematic errors.

A closely related quantity is the magnetic field strength at the light cylinder $B_{LC}$, which using the current P and $\dot{P}$ listed in Table 10, can be calculated using the following equation (check Goldreich & Julian 1969; Abdo et al. 2010):

[TABLE]

with the radius defined as ($R_{\mathrm{LC}}=cP/2\pi$). This will result in a value of $B_{LC}\approx 0.42\mathrm{G}$ for our source. The reason why $B_{LC}$ is a quantity of interest is because the pulsar wind is equal to the plasma that leaves the magnetosphere at the light cylinder. Therefore, the physical properties of the wind could depend on the physical properties at the light cylinder. To investigate this, we included the $B_{LC}$ value of other systems in Table 10 for comparison purposes. Similar to $B_{NS}$, there seems to be a possible connection between the resulting value of wind magnetization $\eta_{B}$ we inferred from our model and the value of ${B_{LC}}$. However, our sample size needs to be bigger to be able to draw rigorous conclusions and check if this possible trend is real or not.

As for the injected particle spectrum, we get: $p_{1}$ = 1.44 and $p_{2}=2.68$, both are harder than what was previously found by Gotthelf et al. (2014) & Mares et al. (2021). In addition, our model predicts having particles with $E_{min}$ of 1 GeV, 270 times lower than what Mares et al. (2021) found, $E_{break}$ of 1 TeV, 16 times lower than what Mares et al. (2021) found. And quite importantly, particles injected in the PWNe can be accelerated up to $E_{max}$ = 1.24 PeV, which makes HESS J1640-465 to be an electronic PeVatron candidate.

As for the surrounding medium, the $\gamma$-ray emission we observed from this source is due to inverse Compton scattering of the particles inside the PWN off two photon fields on top of the Cosmic Microwave Background (CMB). The first has a temperature of 355 K, which is a typical temperature of the dust inside supernova remnants (Temim et al., 2006) and the other is a very hot photon field with a temperature of T $\approx 16000$ K, which is a typical temperature of massive stars. The change in the particle injection spectrum described in the previous paragraph affects the properties of the colder background photon field inferred from $\gamma$-rays causing it to become hotter as we can see from Table 8.

As previously mentioned in §1, PWNe are composed of highly relativistic electrons which cause them to emit electromagnetic radiation through mainly two dominant mechanisms; synchrotron radiation through the interaction of leptons with the PWN’s magnetic field and inverse Compton radiation through the interaction of these leptons with ambient low-energy photons (e.g, Gaensler & Slane 2006). In order to fit the updated X-ray spectrum presented in this paper, a different particle injection spectrum than that reported by Mares et al. (2021) was needed. This invariably affected the background photon fields needed to reproduce the $\gamma$-ray measurements. As a result, the aforementioned highlighted differences between our model parameters and that of Mares et al. (2021), namely the injection particle spectrum and background photon fields, can be directly attributed to our updated measurements of the X-ray spectrum of the PWN.

Lastly, for the PWN’s evolution phase, our model shows that the PWN has started to interact with the SNR’s reverse shock as indicated by their radii values listed in Table 9 ($R_{rs}<R_{PWN}$). This result is in line with the asymmetrical morphology of the source that we can clearly see in Figure 2.

6 Summary and Conclusions

We have analyzed archival X-ray (Chandra and NuSTAR) data for HESS J1640-465. The high value of $N_{H}$ along the line of sight suggests that scattering will significantly affect measurements of the pulsar spectrum. Consequently, we used the $xscat$ model incorporated in XSPEC to quantify these effects in the puslar spectrum extracted from the Chandra data so we can propagate them to the NuSTAR data.

We analyzed the X-ray spectrum of the PWN in two bands; soft X-ray in the 3$-$10 keV energy range, and hard X-ray in the 10$-$30 keV energy range. Our spectral analysis showed evidence for spectral softening across bands which ultimately resulted in having range of values for the PWN properties in both energy bands. The values were reported in Table 7.

We then used a one-zone time-dependent model for the evolution of a PWN inside a SNR to reproduce the observed dynamical and multi-wavelength spectral properties of HESS J1640-465. Our model shows that PWNe can be PeVatron candidates as there are particles that can be accelerated up to $E_{max}$ = 1.24 PeV for our system. Also, Our model requires the associated pulsar to have a high initial spin down luminosity of $\dot{E}_{0}\approx 1.2\times 10^{42}$ erg $s^{-1}$ and a short spin-down timescale of $\tau_{sd}\approx 5$ years compared to other pulsars as shown in Table 10. In addition, our value for the wind magnetization ( $\eta_{B}$ = 0.07 ) is higher than average and only comparable to the value obtained for Kes 75 (Gotthelf et al., 2021) as detailed in Table 10. Lastly, there seems to be a possible connection between the inferred wind magnetization ($\eta_{B}$) values and the values of the surface dipole magnetic field ($B_{NS}$) and magnetic field strength at the light cylinder ($B_{LC}$). However, more systems need to be modelled in order to draw rigorous conclusions about this.

Acknowledgement

The contributions of JDG and SMS was supported by the National Aeronautics and Space Administration (NASA) under grant number NNX17AL74G issued through the NNH16ZDA001N Astrophysics Data Analysis Program (ADAP). JDG and SH are also supported by the NYU Abu Dhabi Research Enhancement Fund (REF) under grant RE022. The research of JDG is also supported by NYU Abu Dhabi Grant AD022. The work of J.D.G. and S.M.S. is also supported by Tamkeen under the NYU Abu Dhabi Research Institute grant $CAP^{3}$.

This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. It also made use of data from the NuSTAR mission, a project led by the California institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. This research made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). We also made use of data and asoftware provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. We also made use of NASA’s Astrophysics Data System Bibliographic Services.

Facilities: CXO, NuSTAR

Software: Xspec (Arnaud, 1996), CIAO (Fruscione et al., 2006), HEASoft (Nasa High Energy Astrophysics Science Archive Research Center (2014), Heasarc), Sherpa (Freeman et al. 2001; Doe et al. 2007) and SAOImageDS9 (Joye 2006; Smithsonian Astrophysical Observatory 2000)