NuSTAR DISCOVERY OF A YOUNG, ENERGETIC PULSAR ASSOCIATED WITH THE LUMINOUS GAMMA-RAY SOURCE HESS J1640−465

To search for pulsations from CXOU J164043.5−463135, we examined NuSTAR data obtained during the two Norma Arm survey pointings, adjacent in time. The day-long data set yields a total of 913 photons from the two focal plane modules, extracted in the optimal 3–25 keV energy bandpass using a 30'' radius aperture centered on the source location. Photon arrival times were corrected for clock drift and converted to barycentric dynamical time (TDB) using the JPL DE200 ephemeris and the Chandra coordinates of the point source. The arrival times were binned at 2 ms and searched for coherent pulsations up to the Nyquist frequency using a 2²⁶ bin fast Fourier transform (FFT). We found a 206 ms coherent signal with a power of P_FFT = 58, significant at the 99.997% (≃ 4σ) confidence level. This motivated a second-epoch observation from which we were able to confirm the detection (P_FFT = 77) and measure the spin-down rate of the pulsar as described next.

We computed the spin-down rate from the difference in frequency over the 100 day time span between the June and September observations. For each data set, we generated a refined frequency measurement by oversampling the signal using the $Z^2_1$ test statistic around the known frequency (Table 1). We then re-fitted for frequency by including the frequency derivative in the light-curve folds, as each observation is sufficiently long that spin-down smears the pulsed signal (Δϕ = 0.11 and 0.32 cycles for the first and second epoch, respectively). The final period measurement used the iterated period derivative reported in Table 2. The uncertainties are estimated from a Monte Carlo simulation of the light curve using the method described by Gotthelf et al. (1999). The final power spectra are shown in Figure 2 (left). The derived physical parameters in Table 2 are from the relations $\dot{E}=4\pi ^2I\dot{P} P^{-3}$, where I = 10⁴⁵ g cm², and the dipole spin-down relations $\tau _c\equiv P/2\dot{P}$ and $B_s=3.2\times 10^{19}(P\dot{P})^{1/2}$ G.

Zoom In Zoom Out Reset image size

Table 2. Timing Parameters of PSR J1640−4631

Parameter	Value
R.A. (J2000.0)a	16h40m4352
Decl. (J2000.0)a	−46°31'354
Epoch (MJD TDB)	56,466.0
Frequencyb, f	4.843950957(40) Hz
Frequency derivativeb, $\dot{f}$	−2.290(10) × 10−11 Hz s−1
Periodb, P	0.206443048(33) s
Period derivativeb, $\dot{P}$	9.758(44) × 10−13
Spin-down luminosity, $\dot{E}$	4.4 × 1036 erg s−1
Characteristic age, τc	3350 yr
Surface dipole magnetic field, Bs	1.4 × 1013 G

Notes. ^aChandra position from Lemiere et al. (2009). ^b1σ uncertainties given in parentheses.

Download table as:  ASCIITypeset image

Figure 2 (right) shows the highly modulated pulse profile characterized by a relatively sharp peak and broad trough compared to a pure sinusoid signal. The pulsed fraction in the 3–25 keV band is f_p > 48 ± 10% after allowing for the background level in the source aperture, estimated using data from a concentric annulus. This increases to f_p ≈ 82% after taking into account contamination from the PWN in the source aperture using the spectral results presented below. Here, we define the pulsed fraction as the ratio of the pulsed emission to the net source (background subtracted) flux. To determine the unpulsed level we take the average of the two lowest bins in the 16-bin folded light curve of Figure 2. The pulse shape shows no variation with energy, within statistics.

For spectral study we extracted photons from the observations listed in Table 1 using an elliptical region of diameter 34 × 26 whose major-axis is oriented at position angle −40°, centered on the apparent enhanced PWN emission at (J2000) 16^h40^m4205, −46°31'478, which is offset by 20'' from the pulsar (see Figure 1). We estimate the background using nearby 50'' radius circular regions carefully chosen to account for stray light in the focal planes in two of the observations. Response matrices were generated for each spectral file using the NuSTAR analysis software. All eight spectra were combined, as were their matching response matrices, and fitted using the XSPEC software package (Arnaud 1996) with χ² minimization. The fitting is limited to the 3–20 keV range due to low signal-to-noise at higher energies. In this range, the combined spectrum yields a total of 14,100 source plus background photons, and is grouped to obtain >5σ signal-to-noise per spectral bin.

The NuSTAR spectrum in the large aperture is dominated by PWN emission and is well fit by an absorbed power-law model as expected for non-thermal emission from the nebula. We used the tbabs absorption model, with Wilms et al. (2000) abundances and Verner et al. (1996) cross sections, and obtain a best-fit N_H = (1.7 ± 0.9) × 10²³ cm⁻² and Γ = 1.9 ± 0.4. Errors are 90% confidence level (Δχ² = 4.61 for two interesting parameters) determined from the error ellipse contour extrema. This provides an acceptable fit with a reduced $\chi ^{2}_{\nu }=0.87$ for 36 degrees of freedom (dof). The absorbed 2–10 keV flux is $(8.0^{+0.4}_{-2.0})\times 10^{-13}$ erg cm⁻² s⁻¹ (90% confidence). This is consistent with the values reported by Funk et al. (2007) using XMM-Newton data extracted from a 25 diameter aperture centered on the PWN.

To further constrain the PWN emission and to isolation the pulsar contribution to the NuSTAR spectrum, we include the Chandra spectrum of the pulsar and PWN in the fitting process. We supplemented the Chandra data set (ObsID 7591) previously analyzed as part of the study of HESS J1640−465 by Lemiere et al. (2009) with data obtained on 2011 June 6 (ObsID 12508). These data were acquired with the Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2003) and reprocessed using the chandra_repro package of the Chandra Interactive Analysis of Observations (CIAO) software suite. Spectra and their response matrices from the two observations were produced using specextract. The pulsar and PWN spectra are grouped with a minimum of 10 and 40 counts per spectral bin, respectively, and each fitted with the absorbed power-law model, using a common column density.

We used an extraction radius of 2'' for the pulsar; the background is negligible in this small region. For the PWN we used the NuSTAR elliptical region and similar background region. We combined spectra from the individual observations, and their response matrices, as above. Although the ACIS-I detector is sensitive in the 0.5–10 keV energy range, we restricted the spectral fits to the 3–7 keV range due to the high absorption and limited statistics. A total of 142 and 1369 photons were fitted for the pulsar and PWN, respectively, during the 47 ks exposure. The second column of Table 3 presents the resulting spectral fits for the pulsar and PWN spectra using the Chandra data alone.¹³

Table 3. X-Ray Spectrum of PSR J1640−4631 and its Wind Nebula

Parameter	Chandra Only	Chandra +NuSTAR
NH (cm−2)	(1.2 ± 0.6) × 1023	(1.8 ± 0.6) × 1023
ΓPSR	$1.2^{+0.9}_{-0.8}$	$1.3^{+0.9}_{-0.5}$
FPSR (2–10 keV)	$1.9^{+0.2}_{-1.4} \times 10^{-13}$	(1.8 ± 0.4) × 10−13
ΓPWN	$2.3^{+1.2}_{-1.0}$	$2.2^{+0.7}_{-0.4}$
FPWN (2–10 keV)	$5.4^{+0.6}_{-2.3} \times 10^{-13}$	(5.5 ± 0.8) × 10−13
$\chi ^2_{\nu }$ (dof)	1.0 (56)	0.82 (79)

Notes. Absorbed power-law model fits to the pulsar and the PWN spectra with their column densities linked. The simultaneous fit to NuSTAR and Chandra data is described in the text. The uncertainties are 90% confidence limits determined from the error ellipse extrema. The given fluxes are absorbed, in units of erg cm⁻² s⁻¹.

Download table as:  ASCIITypeset image

We fitted the Chandra and NuSTAR spectra together using two power laws as before, with their column densities tied. To constrain the pulsar and PWN contributions to the NuSTAR spectrum, we fixed the power-law indices to the Chandra models and tied their flux together in the overlapping 2–10 keV band. To allow for flux calibration differences between the two missions we introduce an overall normalization constant to the NuSTAR model, with best fit value of 1.11. The resulting spectrum is shown in Figure 3, and the parameters are reported in Table 3. The NuSTAR spectra of the pulsar and PWN are evidently successfully modeled, as the fitted values are in good agreement with the Chandra results for each component. To estimate the PWN contribution to the pulse profile shown in Figure 2, we repeated our spectral analysis using a 30' radius aperture; the PWN is found to account for ≈53% of the background-subtracted source flux in pulse profile.

Zoom In Zoom Out Reset image size

The discovery of PSR J1640−4631 supports the conjecture that CXOU J164043.5−463135, HESS J1640−465, and 1FHL J1640.5−4634 are all manifestations of a middle-aged pulsar/PWN system. For the measured distance of 12 kpc to G338.3−0.0 (see below), the 0.2–10 TeV luminosity of HESS J1640−465 is 6% of the pulsar's present $\dot{E}$, while the ratio F(0.2–10 TeV)/F_PWN(2–10 keV) ≈ 13, indicating that, in a leptonic model where the TeV emission results from up scattering of ambient photons by relativistic electrons, inverse Compton losses now dominate over synchrotron emission.

A result that was perhaps unexpected is the young characteristic age of the pulsar. The pulsar's characteristic age, $\tau _c \equiv P/2\dot{P} =$ 3350 yr, is an approximate measure that applies when a simple vacuum dipole spin-down model is assumed, and when the spin period is significantly larger than that at birth. Although pulsars for which the model has been tested do not rotate like simple vacuum dipoles (Kaspi et al. 2001), it is notable that, excluding magnetars, there are only seven pulsars with τ_c < 3000 yr in the Australia Telescope National Facility (ATNF) catalog¹⁴ (Manchester et al. 2005), followed by PSR J1640−4631 with τ_c = 3350 yr. This is younger than most authors have assumed in modeling the evolution of HESS J1640−465/G338.3−0.0 (T = 20 kyr, Funk et al. 2007; T = 15 kyr, Lemiere et al. 2009; T = 10 kyr, Slane et al. 2010). Fang & Zhang (2010) presented two models similar to that of Slane et al. (2010), one with $\dot{E}=1\times 10^{38}$ erg s⁻¹ and T = 4500 yr, and the other with $\dot{E}=1.65\times 10^{37}$ erg s⁻¹ and T = 8200 yr. While these ages are closer to the likely true age, the actual spin-down power of PSR J1640−4631 is considerably smaller than the assumed values.

Prior models for G338.3−0.0 estimate the age of the SNR from its radius assuming that it is in the Sedov phase with $r_s=(2.02\,{\cal E}/\rho)^{1/5}\,t^{2/5}$ (Spitzer 1978). These estimates are imprecise due to the unknown energy $\cal E$ of the explosion and the ambient mass density ρ, and the strong dependence of age t on radius, which magnifies any uncertainty in distance. Also, the pulsar is offset from the center of the SNR in a direction that, given the morphology of the PWN, cannot be explained by high kick velocity. Lemiere et al. (2009) conclude that the geometrical center of the remnant is not the explosion site. Instead, the structure of the remnant may be affected by complicated interactions with the local ISM, which could affect the estimation of its age.

There is no near/far distance ambiguity for G338.3−0.0 because 21 cm absorption is seen up to the maximum negative value at the tangent point in both it and adjacent H ii regions. Lemiere et al. (2009) concluded that the distance to G338.3−0.0 is in the range 12–13.5 kpc, while Kothes & Dougherty (2007) had already found a value of ≈11–12 kpc from the same 21 cm data, but assuming a Galactic center distance of 7.6 kpc instead of 8.5 kpc. We adopt d = 12 kpc here. While the distance appears to be constrained well, there is no direct measurement of the ambient density; the above referenced models assume values of n_ISM between 0.1 cm⁻³ and 10 cm⁻³. Castelletti et al. (2011) estimated an electron density of 100–165 cm⁻³ for the western and northern part of the SNR shell by assuming that the low-energy turnover of its radio spectrum is due to local free-free absorption. They also searched the NANTEN CO survey data for molecular gas that could be target material for hadronic production of γ-rays from protons accelerated in the SNR shell, but did not find any that could be associated with G338.3−0.0.

HESS Collaboration (2014) present a model in which the TeV emission is produced by hadronic interactions in the north and west part of the SNR shell. It is motivated largely by the continuity of the Slane et al. (2010) Fermi spectrum with the HESS spectrum, which they argue is difficult to fit in a PWN model, and also by the overlap of the TeV source with the side of the SNR shell that is adjacent to H ii regions. They require a high ambient density of 150 cm⁻³ to reproduce the GeV/TeV spectrum, but a small density to account for the size of the SNR shell. These are achieved by assuming that the explosion occurred inside a wind-blown bubble of the progenitor star, of average density 0.1 cm⁻³, which is then consistent with their assumed age of 2500 yr. As they remark, the model requires a large fraction of the SN energy to be channeled into cosmic rays, especially since only one side of the SNR contributes to this process.

As we argue below, the Slane et al. (2010) Fermi spectrum used by HESS Collaboration (2014) is not likely to be correct, and there may not even be a detectable source in Fermi at energies <10 GeV. This is why we use different Fermi results in our own model, which contributes to the diverging conclusions of the two papers.

Slane et al. (2010) modeled the γ-ray emission from HESS and the associated broad-band spectrum with an evolving, one-zone PWN model, where the γ-rays are ambient photons inverse Compton scattered by the relativistic electrons that also produce the synchrotron emission. Crucial to these models are the age and spin-down power of the pulsar. Slane et al. (2010) assumed a remnant age of 10 kyr and a pulsar with spin-down power of $\dot{E}=4\times 10^{36}$ erg s⁻¹. The latter came from the $\dot{E}/L_x$ relation of Possenti et al. (2002) and, given the observed scatter, it was fortuitously a very accurate prediction. It was necessary to add a Maxwellian distribution of electrons around 0.1 TeV to a power-law tail to account for the strong Fermi GeV γ-rays that they found from this source relative to its TeV flux.

The pulsar age assumed by Slane et al. (2010) was likely too large. The true age of a pulsar is given by

$$\begin{equation*} T = {P \over (n-1)\dot{P}}\left[1-\left(P_0 \over P\right)^{n-1}\right], \end{equation*}$$

where $n\equiv f\ddot{f}/\dot{f}^2$ is the braking index, assumed by Slane et al. (2010) to be 3, consistent with a vacuum dipole, P₀ is the spin period at birth, and P is the currently measured period. For n = 3 the pulsar's true age is less than its characteristic age, and refined models will likely favor a younger system and/or a smaller braking index.

Assuming again that the PWN of PSR J1640−4631 is primarily responsible for the γ-ray emission, we fit an evolutionary model of a PWN inside an SNR to the broad-band spectrum. This can constrain the properties of the central neutron star, the pulsar wind, progenitor supernova, and the surrounding ISM. We include the PWN emission in the X-ray band and the radio upper-limit reported in the literature (see Table 4). However, in examining more recent Fermi publications, we note that the source 2FGL J1640.5−4633 (Nolan et al. 2012) is fainter than the flux (based on ∼1 yr of data) that was extracted by Slane et al. (2010). More recent analyses of this Fermi source have been restricted to energies >10 GeV (Acero et al. 2013; Ackermann et al. 2013), and we use only the latter in our revised modeling, in particular the spectrum of 1FHL J1640.5−4634 (Ackermann et al. 2013) derived from three years of data. Indeed, it is not obvious from inspecting Fermi counts maps that there is a significant source at energies <10 GeV, which is why we do not attempt to model this part of the Fermi band. This is an important difference from HESS Collaboration (2014), who accept and model the Fermi results of Slane et al. (2010).

Table 4. Properties of HESS J1640−465

Parameter	Observed	Modeled	Reference
SNR radius	45 ± 05	451	Shaver & Goss (1970)
Sν [660 MHz] (mJy)	⩽690	337	Giacani et al. (2008)
FX [2–25 keV] (erg cm−2 s−1)	(1.68 ± 0.4) × 10−12	1.0 × 10−12	This work
ΓX	$2.2^{+0.7}_{-0.4}$	2.4	"
Fγ [10–500 GeV] (erg cm−2 s−1)	(3.15 ± 1.00) × 10−11	1.6 × 10−11	Ackermann et al. (2013)
Γγ	1.92 ± 0.24	1.96	"

Note. The quoted 2–10 keV flux is corrected for interstellar absorption.

Download table as:  ASCIITypeset image

Our model is based on the work of Gelfand et al. (2009), modified to include background photon fields in addition to the cosmic microwave background (CMB), and a broken power-law spectrum of particles injected into the PWN at the termination shock, as favored by recent studies of these objects (e.g., Bucciantini et al. 2011). We assumed two additional background photon fields: one with temperature T₁ = 15K and energy density u₁ = 4u_CMB, where u_CMB = 4.17 × 10⁻¹³erg cm⁻³ is the energy density of the CMB, and the other with T₂ = 5000 K and u₂ = 1.15u_CMB, the same as used by Slane et al. (2010). We fixed the distance at 12 kpc. To obtain the best-fit values of the free parameters and their errors, which are listed in Table 5, we employed a Markoff Chain Monte Carlo algorithm similar to the one used to fit the properties of the PWNe in Kes 75 (Gelfand et al. 2014) and G54.1+0.3 (J. D. Gelfand et al. 2014, in preparation).

Table 5. PWN Evolutionary Model Results for HESS J1640−465

Parameter	Value
SN explosion energy, log (${\cal E}/10^{51}$ erg)	−0.78
Ejecta mass, log (Mej/M☉)	−0.34
ISM density, log (nISM/cm−3)	−1.53
Braking index, n	+1.87
Spin-down timescale, log (τ0/yr)	+2.89
Wind magnetization, log (ηB)	−1.76
Low-energy particle index, p1	+0.28
High-energy particle index, p2	+2.59
Minimum energy, log (Emin/GeV)	−3.86
Break energy, log (Ebr/TeV)	−0.71
Maximum energy, log (Emax/PeV)	+0.12

Download table as:  ASCIITypeset image

For a constant braking index, the spin-down power of the pulsar evolves as

$$\begin{equation*} \dot{E}(t)=\dot{E}_0\left(1+{t \over \tau _0}\right)^{-({n+1 \over n-1})}, \end{equation*}$$

where $\dot{E}_0$ is the initial spin-down power and τ₀ is the initial spin-down timescale. The present P and $\dot{P}$ are known, which fixes the values of $\dot{E}_0$ and the present age T corresponding to each trial value of n and τ₀.

As shown in Figure 4, the model is able to reproduce the observed broad-band spectrum with χ² = 24.7 for 15 degrees of freedom, using the best-fit parameters given in Table 5. In general, the values of these parameters are similar to what has been inferred for other systems using similar models (e.g., Bucciantini et al. 2011). However, our best-fit model requires a small braking index n ≈ 1.9 and a short initial spin-down timescale τ₀ ∼ 800 yr (see Figure 5). This implies that the pulsar was born spinning rapidly, with P₀ ≈ 15 ms and a high initial spin-down power $\dot{E}_0\approx 10^{40}$ erg s⁻¹. Its present age would then be 6800 yr. We note that, in this case, the model is not entirely self-consistent, as the initial rotational energy of the pulsar, 9 × 10⁴⁹ erg, is a significant fraction of the fitted SN explosion energy, and could affect the dynamics of the explosion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size