Long-term Study of the Double Pulsar J0737-3039 with XMM-Newton: Spectral Analysis

14 years after its discovery, the double pulsar PSR J0737-3039 is still studied extensively. This system is composed of two neutron stars: an old, fast, mildly recycled 22.7 ms pulsar (hereafter referred to as Pulsar A; Burgay et al. 2003) and a younger and much slower pulsar with a period of 2.77 s (hereafter Pulsar B; Lyne et al. 2004). The compact objects orbit each other in a very tight orbit in only 2.4 hr, with mean orbital velocities of about 1 million km hr⁻¹ (see, e.g., Kramer & Stairs 2008, for a review).

PSR J0737-3039 represents the most compact relativistic system, and the only binary system in which both neutron stars have been detected as radio pulsars. The Keplerian and post-Keplerian parameters offer a unique test for theories of strong field gravity (Lyne et al. 2004; Kramer et al. 2006; Kramer & Stairs 2008), and constrain the masses of both neutron stars with high accuracy. The physical riches of the double pulsar have been highlighted from radio observations (Kramer & Stairs 2008; Kramer & Wex 2009): whereas the pulse profile of Pulsar A has been stable, the pulsed emission of Pulsar B showed strong orbital flux and profile variations, before disappearing entirely in 2008 (Perera & McLaughlin et al. 2010). Moreover, as a direct consequence of the inclination of a system observed nearly edge-on (i $\sim \,89^\circ$), radio eclipses were detected when Pulsar A passed behind Pulsar B (Kaspi et al. 2004; McLaughlin et al. 2004; Breton et al. 2008). The study of the light curves testified to a new type of interaction between Pulsar A's relativistic particle wind and Pulsar B's magnetosphere in the previously unexplored close environment between the two neutron stars. This aroused the interest to carry out high-energy observations in the X-ray band to investigate further the intra-binary environment.

Early observations performed by the Chandra and XMM-Newton satellites (Campana et al. 2004; McLaughlin et al. 2004; Pellizzoni et al. 2004, 2008; Chatterjee et al. 2007; Possenti et al. 2008) pointed out the difficulty in constraining the origin of the multifold nature of the double pulsar's X-ray emission. The non-thermal pulsed emission from Pulsar A, although very soft (power-law slope ${\rm{\Gamma }}\sim 3.3$), is clearly predominant in the X-ray flux (Chatterjee et al. 2007; Pellizzoni et al. 2008; Possenti et al. 2008).

In the frame of a first XMM-Newton "Large Program" of 235 ks exposure in 2006, pulsed X-ray emission from Pulsar B was also detected in part of the orbit. Due to its own low rotational energy loss, X-ray emission from Pulsar B can be only powered externally through the spin-down energy of Pulsar A (Pellizzoni et al. 2008). This emission, consistent with thermal radiation of temperature kT ∼ 30 eV and bolometric luminosity of $\sim {10}^{32}$ erg s⁻¹, was ascribed to the heating of Pulsar B's surface by Pulsar A's wind. A hotter (∼130 eV) and fainter ($\sim 5\times {10}^{29}$ erg s⁻¹) thermal component, possibly originating from backfalling material heating the polar caps of either Pulsar A or Pulsar B, was also suggested from the spectral analysis by Pellizzoni et al. (2008).

In 2011, a second deep XMM-Newton observation of 370 ks was carried out. Comprehensive timing analysis over such a large time span has confirmed X-ray pulsed emission from Pulsar B even after its radio disappearance (Iacolina et al. 2016). The unusual phenomenology of Pulsar B's X-ray emission includes orbital pulsed flux and profile variations, as well as a loss of pulsar phase coherence on the timescale of years, suggesting orbital-dependent penetration of Pulsar A's wind plasma onto Pulsar B's closed field lines. Furthermore, timing analysis of the full XMM-Newton data set provided the first evidence of orbital flux variability (7 ± 1%), possibly involving a bow shock between pulsar structures (Iacolina et al. 2016). An additional, possibly hard, spectral component associated with this intra-binary environment is then expected.

In this paper, we present a comprehensive spectral analysis of the full XMM-Newton data set from both large programs, which allows us, for the first time, to constrain high-energy components above 4 keV. Because of the weakness of the source and the consequent need to scrutinize the low count statistics, an in-depth revision of the background subtraction procedure is provided. The standard spectral analysis was coupled with an independent source detection method based on likelihood spatial analysis. The complementarity of both approaches allows us to speculate, for the first time, on the possible presence of an iron line in the double pulsar.

Both XMM-Newton Large Programs were carried out with similar instrumental configurations, suitable for simultaneous spectral and high-resolution timing analysis of the double pulsar. For these purposes, the EPIC-pn (Strüder et al. 2001) was used in Small Window mode whose time resolution of 5.67 ms compares with Pulsar A's period of 22.7 ms. The EPIC-MOS cameras (Turner et al. 2001) were also operated in Small Window mode with a time resolution of 0.3 s only suitable for Pulsar B timing analysis. Medium and thin optical filters were applied to the EPIC instruments in 2006 and 2011, respectively. The thin filter nearly doubles the instrument effective area in the soft band (0.3–0.5 keV) at the expense of higher background contamination: early X-ray observations of the double pulsar showed a very soft source spectrum motivating the use of the thin filter for subsequent observations in order to improve the overall counting statistics.

The XMM-Newton data were processed using Science Analysis Software (SAS) version 12. The calibrated and concatenated EPIC event lists were obtained by running the meta-tasks epproc and emproc, standard pipeline tasks for EPIC-MOS and EPIC-pn observations, respectively. We performed a barycentric correction on the event files with the task barycen, using the JPL DE405 ephemeris to convert the time of arrival of photons from the satellite position to the solar-system barycenter.

Because of the faintness of the double pulsar in X-rays, careful attention was provided in the assessment and mitigation of background contamination and in the selection of source photons. Strong flares, most likely associated with cosmic soft proton events, significantly affected different observations. In order to discard these episodes, we produced the light curves between 10 and 12 keV to determine the flaring times. We followed a background rejection procedure according to the general prescriptions by de Luca & Molendi (2004) and further investigated and tested differing filtering options in analogy with the espfilt SAS task, and reviewed the procedure adopted by Iacolina et al. (2016) and Pellizzoni et al. (2008) for the same data. Two different photon pattern selections were considered for the distribution of incident photons over the pixels of the EPIC-pn (carrying out the analysis for both options): single and double events (pattern $\leqslant 4$) or single events only (pattern = 0), slightly reducing the effective area. As for the EPIC-MOS, single, double, triple, and quadruple events were selected (pattern $\leqslant 12$).

In parallel, we cross-checked our analysis by adopting another method, which consists in analyzing the average quiescent count rates free from flaring particle background. We extracted a 10 s binned light curve in the 0.15–10 keV energy range from the whole field of view of each camera, omitting the source region, selecting events below a threshold of $5\sigma$ of the quiescent rate for the EPIC-pn for pattern = 0 selection, and $4\sigma$ for pattern $\leqslant 4$. As for the EPIC-MOS, the threshold was defined at $5\sigma$ from the quiescent count rate. The resulting effective observing time after this operation is reported in Table 1 together with the original total exposure time for each EPIC instrument (no significant dead-time differences are related to the two photon pattern options considered for the EPIC-pn). Alternative photon pattern selections and threshold levels provide consistent results, although the choice of pattern $\leqslant 4$ for the EPIC-pn is best suited for hard X-ray spectral analysis providing a higher effective area at the expenses of a background increase only in the softer energy channels ($\lt 0.4$ keV).

Once the EPIC event lists were cleaned from strong background contamination, we extracted the source and background spectra. The associated source's radius of extraction was set at $18^{\prime\prime}$ for the EPIC-pn and $15^{\prime\prime}$ for the EPIC-MOS, providing maximum source detection significance (pulsed emission from Pulsar A), and we carefully checked that residual flare contamination was not present within the selected region. The background spectra were extracted from different rectangular regions for consistency checks (see Figure 1). The task backscale was used to calculate the area of source and background regions. Standard energy response matrices (RMFs; rmfgen) and effective area files (ARFs; arfgen) were created for each spectrum via the standard SAS tools. The background filtering and source extraction procedures provided an overall harvest of ∼15,000 source X-ray photons considering the EPIC-pn, MOS1, and MOS2 data from both large programs.

Zoom In Zoom Out Reset image size

Up to five significant spectral components could in principle be expected from the double pulsar in X-rays according to theoretical models: surface thermal emission and non-thermal magnetospheric emission from both neutron stars, and orbital phase dependent bow-shock emission due to the interaction between Pulsar A's wind and Pulsar B's magnetosphere. This latter non-thermal component, if present, was proved to be much weaker than expected: no clear evidence of shock emission was reported in the spectral analysis of 2006 data. While the spectrum related to that data set was correctly fitted by a two-component model (a power law plus a blackbody or two blackbodies), the phase-resolved analysis indicated that the spectrum was clearly dominated by pulsed non-thermal emission from Pulsar A and the off-pulse required two-components thermal emission models (Pellizzoni et al. 2008). Furthermore, only very weak orbital flux modulation (7%) was detected as a possible signature of a bow-shock in the overall timing analysis of both large programs (Iacolina et al. 2016). Thus, a two- (pl+bb) or three-component (pl+bb+bb) pulsar emission model could satisfactorily account for most of the X-ray flux from the double pulsar.

We used the XSPEC version 12.8 (Arnaud 1996). All uncertainties are given at the 90% confidence level (${\rm{\Delta }}{\chi }^{2}=2.706$). The EPIC-pn, MOS1, and MOS2 data were rebinned to have at least 25 counts per energy channel. The bulk of the X-ray emission from PSR J0737-3039 is at low energy, below ∼1 keV. However, in order to better constrain the models, spectral analysis was performed in the 0.15–11 keV energy range. The probed models were all modified at low energies to account for interstellar photoelectric absorption through the XSPEC Tbabs model (Wilms et al. 2000). Tbabs calculates the cross section for X-ray absorption by the ISM as the sum of the cross sections for X-ray absorption due to the ISM in its gas, grain, and molecular phases.

We first verified that the spectral analysis of the observations performed in 2011 provides results comparable with those related to 2006 data alone (i.e., no significant spectral variability is evident over years). For the sake of completeness, we tested both two-component and three-component fits as suggested by the phase-resolved studies on the 2006 data (Pellizzoni et al. 2008). The models were applied on the EPIC-pn data, then EPIC-pn plus MOS1 and MOS2 data. The parameters associated with the different models are reported in Table 2. The results obtained separately for the two large programs are perfectly in agreement (see Pellizzoni et al. 2008). The observed flux related to the 2011 data in the 0.2–3 keV energy range is $\sim 4.3\times {10}^{-14}$ erg cm⁻² s⁻¹ (using the power law plus blackbody model). By comparison, the observed flux was estimated at $\sim 4.0\times {10}^{-14}$ erg cm⁻² s⁻¹ during the observations in 2006. We therefore analyzed both data sets together. A summary of model parameters is reported in Table 3.

Table 2.  Comparison of Three Models on the Second Large Program Data

Detector	Modela	NH	Γ	Normpob	kTbb1	Normbb1b	kTbb2	Normbb2b	${\chi }^{2}/$dof	${\chi }_{\mathrm{red}}^{2}$
		(1020 cm−2)			(eV)		(eV)
pn	1	${0.8}_{-0.7}^{+1.2}$	2.9 ± 0.3	5.3 ± 1.2	${150}_{-10}^{+20}$	2.0 ± 0.6			239/206	1.16
	2	-c			100 ± 8	3.8 ± 0.3	250 ± 30	${2.0}_{-0.3}^{+0.4}$	254/206	1.23
	3	-	2.6 ± 0.4	${3.4}_{-1.4}^{+1.2}$	110 ± 30	${2.1}_{-1.9}^{+0.8}$	${210}_{-100}^{+160}$	${1.1}_{-0.9}^{+1.1}$	237/204	1.16
pn+MOS	1	1.5 ± 0.9	3.1 ± 0.2	6.0 ± 0.9	${160}_{-10}^{+20}$	1.8 ± 0.4			355/318	1.12
	2	-			105 ± 6	3.7 ± 0.2	260 ± 20	1.9 ± 0.2	367/318	1.15
	3	≤1.0	${2.7}_{-0.3}^{+0.4}$	2.6 ± 0.1	110 ± 20	${2.4}_{-0.8}^{+0.5}$	${230}_{-40}^{+50}$	1.3 ± 0.5	346/316	1.10

Notes.

^aModels: (1) tbabs*(po+bb), (2) tbabs*(bb+bb), (3) tbabs*(po+bb+bb). ^bThe power-law and blackbody normalizations are in units of 10⁻⁶ photons cm⁻² s⁻¹ keV⁻¹ and 10⁻⁷, respectively. ^cThe "-" in the N_H column indicates that the value is found to be very low and not constrained.

Download table as:  ASCIITypeset image

Table 3.  Comparison of Three Models Using the Data of Both Large Programs

Detector	Modela	NHc	Γ	Normpob	kTbb1	Normbb1b	kTbb2	Normbb2b	${\chi }^{2}/$dof	${\chi }_{\mathrm{red}}^{2}$
		(1020 cm−2)			(eV)		(eV)
pn	1	1.7 ± 0.9	3.0 ± 0.2	6.1 ± 0.9	${147}_{-13}^{+20}$	1.6 ± 0.5			347/333	1.04
	2	-c			105 ± 5	3.9 ± 0.2	${260}_{-20}^{+25}$	${1.8}_{-0.2}^{+0.3}$	360/333	1.08
	3	≤1.7	${2.5}_{-0.6}^{+0.7}$	${2.7}_{-1.3}^{+2.8}$	110 ± 20	${2.7}_{-1.5}^{+0.7}$	${230}_{-55}^{+90}$	${1.2}_{-0.6}^{+0.7}$	343/331	1.04
pn+MOS	1	2.2 ± 0.7	3.2 ± 0.2	6.1 ± 0.7	160 ± 15	${1.6}_{-0.3}^{+0.4}$			533/520	1.02
	2	-			105 ± 5	3.9 ± 0.2	270 ± 20	1.8 ± 0.2	551/520	1.06
	3	≤1.6	${2.5}_{-0.5}^{+0.7}$	${2.3}_{-0.9}^{+3.2}$	110 ± 10	${2.8}_{-1}^{+0.5}$	${230}_{-30}^{+40}$	${1.4}_{-0.7}^{+0.5}$	524/518	1.01

Notes.

^aModels: (1) tbabs*(po+bb), (2) tbabs*(bb+bb), (3) tbabs*(po+bb+bb). ^bThe power-law and blackbody normalizations are in units of 10⁻⁶ photons cm⁻² s⁻¹ keV⁻¹ and 10⁻⁷, respectively. ^cThe "-" in the N_H column indicates that the value is found to be very low and not constrained.

Download table as:  ASCIITypeset image

Since the hypothetical presence of significant bow-shock emission should be associated with a hard photon index (Sturner et al. 1997), we tried to fix the photon index of the power law to a lower value, such as ${\rm{\Gamma }}\sim 2$. The resulting fit leads to a ${\chi }_{\mathrm{red}}^{2}\sim 1.34$. These results are also consistent with Pellizzoni et al. (2008).

The reported weak orbital flux variability of 7% from timing analysis (Iacolina et al. 2016) could in principle be constrained by orbital phase-resolved analysis, in spite of the relatively low count statistics involved. In fact, splitting the orbital phase interval in different sections, we did not find any significant variation of the above spectral components. In particular, we cannot claim any significant spectral changes corresponding to peculiar orbital phases (around neutron star conjunctions and quadratures, and periastron/apoastron passages) where putative X-ray emission could in principle be enhanced due, for example, to a hypothetical anisotropic bow-shock emission (Granot & Mészáros 2004; Lyutikov 2004).

Though perfectly adequate for the soft X-ray band, none of these two- or three-component models properly fits the data above 4 keV. By applying a two- or three-component model from Table 3, an excess is clearly observed at 6–7 keV (see Figure 2), while error bars become very important above 8 keV where the background dominates.

Zoom In Zoom Out Reset image size

The EPIC-pn data obtained during the two Large Programs were studied separately in order to assess the spectral parameters of this newly-recognized hard component, including possible variability. The spectrum obtained in 2006 indicates an excess at 6–8 keV in the two XMM-Newton orbits. A narrow Gaussian line adequately fits this excess (see Figure 3). Its centroid energy is found at $6.2\,\pm \,0.5\,\,\mathrm{keV}$ with the associated line width constrained to $\sigma \lt \,0.8\,\mathrm{keV}$. The line flux $(3\pm 2)\times {10}^{-7}$ photons cm⁻² s⁻¹ keV⁻¹ represents 7% of the total absorption-corrected luminosity ${L}_{0.3-11\mathrm{keV}}$ = $6.1\,\pm \,0.3\times {10}^{30}$ erg s⁻¹ (assuming a distance of 1.1 kpc; Verbiest et al. 2012). No other simple model, such as a power law or a blackbody, is able to fit the high-energy excess. Since only a few source photons are detected at high energy ($\sim 1 \%$ in the 4–8 keV range), it is not possible to assess the improvement of the fit for the whole energy range including the Gaussian line. A spectral fit comparison with or without the line was performed on the restricted energy range 4–8 keV. A power law plus a blackbody component were applied to the data, with the parameters fixed to the best values found in the whole energy range. The associated ${\chi }^{2}$ is 11.4 for 7 dof (${\chi }_{\mathrm{red}}^{2}\sim 1.6$), indicating that an additional spectral contribution at high-energy is required. Indeed, the addition of a Gaussian line at ∼6.2 keV significantly improves the fit, reaching ${\chi }^{2}=3.4$ for 4 dof (${\chi }_{\mathrm{red}}^{2}\sim 0.85$). Comparable results are obtained assuming other models in Table 3.

Zoom In Zoom Out Reset image size

The use of the C-statistic (Cash 1979) instead of the ${\chi }^{2}$ is in principle more appropriate for the fitting procedure with only a few photons per energy bin in the spectra, although its exploitation in the frame of XSPEC is not straightforward. Since the background dominates at energy $\gt 8\,\mathrm{keV}$, we restricted the analysis in the energy range 0.15–8 keV when considering pattern = 0, and 0.3–8 keV in the case of pattern = 0–4. The previously discussed model was applied to the data, i.e., a power law plus a blackbody component. Adopting the C-statistic and pattern selection $\leqslant 4$, the Gaussian line centroid is found at $6.1\,\pm \,0.3\,\mathrm{keV}$, with line width $\sigma \lt \,0.8\,\mathrm{keV}$ and normalization $\sim 3\times {10}^{-7}$ photons cm⁻² s⁻¹ keV⁻¹, in agreement with results obtained through the use of ${\chi }^{2}$ statistic. As for the case of the ${\chi }^{2}$ statistic, in order to quantify the improvement of the fit provided by the introduction of a Gaussian component, we compared the values of the C-statistic using continuum models in the restricted 4–8 keV band, with or without the inclusion of the line. By applying only the continuum to the data, a C-statistic value of 195 for 182 dof indicates a bad fit. The addition of the Gaussian at 6.2 keV definitely improves the C-statistic providing 179.4 for 179 dof. Comparable results are obtained adopting the strictest criterion photon pattern = 0. Since our procedure can be assimilated to an F-test which is not strictly adequate to assess the statistical significance of the line (Protassov et al. 2002), we also provided spectral simulations of fake data based on power law plus blackbody models described above. These simulations show that (1) background counts dominate at $E\gt 8\,\mathrm{keV}$ and (2) continuum emission counts drops above 5 keV. Chance probability of a fluctuation of source continuum emission mimicking a line at 6–7 keV appears negligible. Even after 400 trials, fake continuum data never required the inclusion of a line feature in the fit, confirming line detection above 3-sigma level.

In the case of the second long-observation performed in 2011, the high-energy component appears to be different from that detected in 2006, with weaker evidence of a high-energy spectral line (see Figure 4). However, it is worth noting that the quality of the data is not as good as during the observation performed in 2006, with a stronger background emission partly due to the use of the EPIC thin filter. We simulated fake data including a continuum (power law plus blackbody) model and a line feature with the same parameters derived from 2006 data, assuming the background model obtained from 2011 data. An emission line with the same energy and flux as in the 2006 observation would appear at the limit of detectability in 2011 data.

Zoom In Zoom Out Reset image size

Unfortunately, orbital phase-resolved analysis did not provide any significant constraints at E $\gt \,4\mbox{--}5\,\mathrm{keV}$ because of the low count statistics, and because pulsar phase-resolved analysis was not suitable in hard X-rays.

In order to rule out that the reported hard X-ray features are associated with background emission, the standard spectral analysis was coupled to an independent source detection method based on likelihood spatial analysis. The most sensitive procedure to estimate the strength of sources in X-ray data in common with other photon-counting data entails constructing composite models of the spatial distribution of events combining source and background components (Pollock 1987). The background near the double pulsar is made up of cosmic and instrumental contributions and is assumed to be flat.

A point spread function spatial model was applied to EPIC images between 4 and 8 keV, where the discrepancy between the low-energy spectral models and the hard X-ray component is the most relevant. Table 4 reports the results separately and combined for each of the five long observations comprising the two Large Programs in 2006 and 2011. The exposure time reported, T(s), includes the dead-time correction of the Small Window mode. The likelihood detection statistic, ln L, shows that the double pulsar is a weak source at high energies that is difficult to detect in a whole spacecraft revolution of nearly two days, but consistent in strength with the count rate accumulated in half a million seconds of elapsed time.

The combined log-likelihood detection statistic of 15.1 is a secure $\sim 4\sigma$ detection that amounts to 160 ± 40 photons above 4 keV from the X-ray source. Unlike spectral analysis, unpredictable background contamination cannot significantly affect likelihood analysis results, since background counts are not distributed in the image according to the source point spread function. Moreover, the spatial analysis testifies that there is no detection of source photons above 8 keV. Therefore our analysis proves that hard X-ray photons originate from the double pulsar in the 4–8 keV energy range.

The X-ray emission from pulsars is typically attributed to a magnetospheric and/or a thermal origin. Because of the small separation between the two neutron stars, a strong interaction is expected between the wind of Pulsar A and the magnetosphere of Pulsar B, inducing possibly more complex mechanisms of the X-ray emission. In the following, we discuss the origin of the soft X-ray emission (up to 4 keV) and hard X-ray emission, including the detection of the spectral feature at 6–7 keV.

5.1. Soft X-Ray Emission

5.2. Hard X-Ray Emission

Two Large Programs of observations of the relativistic double pulsar were performed by XMM-Newton in 2006 and 2011. The total observing time of about 600 ks offers the opportunity to investigate the low luminosity source with a high number of ∼15,000 source photon events.

No significant spectral variations in soft X-rays can be claimed in the 5 years between 2006 and 2011, substantially confirming the previously reported results based on the analysis of the 2006 data alone. Two-component (blackbody plus power law) or three-component spectral models (double blackbody plus power law) can in principle fit the whole data set, although the latter models seem better suited according to phase-resolved spectral analysis (Pellizzoni et al. 2008). However, for any applicable spectral scenario in the 0.15–3 keV range, the data above 4 keV are not correctly described by the models.

We investigated, for the first time, the high-energy part of the spectrum of the double pulsar. An intriguing emission feature, possibly variable, is detected at about 6–7 keV. Spatial analysis confirms the emission of hard energy photons between 4 and 8 keV ascribed to the source. This feature is most likely attributed to iron line emission, testifying the presence of a relic disk that survived the supernova explosions that terminated the lives of the double pulsar's stellar progenitors.

We acknowledge the referee for his/her constructive suggestions that have helped to improve the content of this paper. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). The authors thank the members of the XMM-Newton User Support Group. E.E. and M.M. acknowledge financial support from the Autonomous Region of Sardinia through a research grant under the program CRP-25399 PO Sardegna FSE 2007–2013, L.R. 7/2007, Promoting scientific research and innovation technology in Sardinia. N.R.I. acknowledges support of the Russian Scientific Foundation under the grant No. 14-50-00043.