ORBITAL VARIATION OF THE X-RAY EMISSION FROM THE DOUBLE NEUTRON STAR BINARY J1537+1155

Double neutron star binaries (DNSBs) are rare systems comprised of two neutron stars (NSs) in a tight orbit. In addition to the gravitational interaction, in which general relativistic effects often play a significant role, the two NSs may also interact in other ways (e.g., via electromagnetic torque or pulsar wind). The manifestations of such interaction(s) can be found via precise radio timing, provided that one of the NSs is seen as a radio pulsar. In the exceptional case of DNSB J0737−3039 (hereafter J0737; see Kramer & Stairs 2008 for recent review), both the NSs are pulsars, which has enabled a suite of high-precision measurements of post-Keplerian parameters (Kramer et al. 2006), determination of the component masses (Lyne et al. 2004), and analysis of the interaction between the wind of pulsar A and the magnetosphere of pulsar B (Lyutikov 2004). The evolutionary processes in which DNSBs may form have been discussed by, e.g., van den Heuvel (2005).

J0737 had been the only DNSB detected in X-rays (McLaughlin et al. 2004), until Kargaltsev et al. (2006, hereafter KPG06) discovered X-ray emission from DNSB J1537+1155 (hereafter J1537). KPG06 observed J1537 with the Advanced CCD Imaging Spectrometer (ACIS) on board Chandra. In the 36 ks exposure, they detected 16 photons, whose energies suggested a steep spectrum (photon index Γ ≈ 3). Intriguingly, KPG06 found that the arrival times did not appear to be evenly distributed in orbital phase, showing a deficit of photons near apastron. However, the statistical significance of the deficit in the first observation was marginal (in fact, as we show below, it was overestimated by KPG06). Here we present analysis of the second Chandra ACIS observation (with a re-analysis of the first one), which provides additional support for the dependence of its emission on orbital phase.

In Section 2, we describe the data analysis, with details about the timing (Section 2.1) and the spectrum (Section 2.2). We present conclusions and a brief summary in Section 3.

In both observations of J1537, the ACIS was operated in the very faint (VF) mode. For each of the data sets, we applied VF-mode background filtering (Garmire et al. 2003), removed positional randomization applied by the pipeline, and applied the sub-pixel resolution tool to split-pixel events (Tsunemi et al. 2001; Mori et al. 2001).

The first observation (ObsID 5591) was carried out on 2005 April 10, with a total exposure of 36080 s (after filtering for good time intervals). The source was imaged on-axis on the most sensitive S3 chip, operated in timed exposure mode with a frame time of 3.24104 s (full-chip readout, giving 457 s dead time and 35623 s net live time). We reprocessed this data set and generated a new Level 2 event file,⁴ in order to ensure consistency with a more recent pipeline version (DS 8.2) used to process the second data set. As a result of slightly different processing, the extraction aperture of r = 123 (used in KPG06) now includes 17 photons (one extra photon compared to KPG06) in the energy range of 0.3–10 keV. For spectral analysis, photons were extracted using the CIAO task psextract (CIAO version 4.2 and CALDB version 4.2.0). Counting the number of background counts in an annulus (492 < r < 4428) centered on J1537, we expect to find (on average) 0.37 background photons within the source extraction aperture.

The second observation (ObsID 20224) was carried out on 2009 September 27 and had a total exposure time of 38181 s. There were no bad time intervals. Again, the ACIS S3 chip was used, this time in 128 row sub-array mode (i.e., 1/8 of the chip), providing 441.04 ms time resolution. The dead time for such short frame time is 9.3%, so the total live time after the dead-time correction is only 34628 s. The less efficient high time resolution mode was used in order to make it possible to detect pulsations in the data, should slow radio pulsations from NS B be detected in the radio. As above, we extracted photons from an r = 123 aperture centered on the source. The total number of photons is 16, in the 0.3–10 keV energy range. Of these, on average, 0.69 photons are expected to come from the background. The image from the second observation is shown in Figure 1. There appear to be no detector artifacts in the vicinity of the source, and the background beyond the extraction circle is quite uniform.

Zoom In Zoom Out Reset image size

2.1. Timing

For the timing analysis, we used the photon arrival times from the Level 2 event files, after correcting them to the solar system barycenter with the axbary tool from CIAO 4.2. For each photon, we calculated the orbital phase using the radio ephemeris (courtesy of Ingrid Stairs). The phase of each detected photon is plotted in Figure 2. There appears to be a portion of orbital phase in which few/no photons are detected from the source. We have checked that the fluctuations in the background rate were consistent with Poissonian statistics and no decrease in the background rate was seen in the phase intervals with no counts detected in the source aperture.

Zoom In Zoom Out Reset image size

As an alternative representation of the same data, Figure 3 shows the phase-folded light curve, with three bins, where we choose one of the bins to be centered at the phase of apastron, ϕ_A = 0.5. We also show confidence intervals for each bin.

Zoom In Zoom Out Reset image size

We will now assess the probability that the observed deficit of photons is caused by random fluctuations, assuming that the number of counts within a certain fixed time interval obeys Poissonian statistics.

Assuming the null hypothesis that the photon arrivals are uniformly distributed in time, the probability of detecting n photons within a time interval t is given by the Poissonian distribution P_Poiss(n, λt), where λ is the mean count rate (see Equation (A3) in the Appendix). In our case, it can be estimated as λ = (N₁ + N₂)/(T₁ + T₂) = 16.154 counts per binary period P_orb, where N₁ = 17 and N₂ = 16 are the numbers of counts detected in the time intervals T₁ = 0.9925P_orb and T₂ = 1.0503P_orb, respectively, and P_orb = 36, 352 s.

The probability of finding an interval between two neighboring events larger than (or equal to) τ at a given event rate λ is given by Equations (A4) and (A1). For the first observation (in which the largest empty interval is Δt₆ = 0.352P_orb), these equations give the probability⁵ P₁(τ₁ ⩾ 0.352P_orb) = 3.81%.

For the second observation, the largest intervals, Δt₁ = 0.140P_orb and Δt₁₇ = 0.134P_orb, are between the beginning of observation and the first count and between the last count and the end of observation. In this case, the exposure time in the gap is the sum of the two intervals at the start and end of the observations (see Figure 2), Δt₁ + Δt₁₇ = 0.274P_orb. The probability of obtaining such an interval by chance for uniformly distributed counts is estimated in the Appendix: P₂(τ₂ ⩾ 0.274P_orb) = 19.0%.

We must assess the probability of finding this set of two gaps together, i.e., how likely is it that we will find by chance gaps as large or larger than observed in two subsequent observations if the photons' arrivals are uniformly distributed in time? We have calculated the probabilities of measuring a gap of a given (or larger) size within a single observation (P₁ and P₂) above. An intuitive choice of statistics for calculating the combined probability is then Z = P₁P₂ (e.g., Fisher 1970). (Note that P₁P₂ is not itself the probability.⁶) Since P₁ and P₂ are both uniformly distributed, the cumulative distribution function is given by P(Z ⩽ z) = z − zln z (Springer 1979, p 102). In our case, z = 0.0072 and P(Z ⩽ 0.007) = 0.043. Exactly the same probability follows from the recipe by Fisher (1970): P(Z < z) = P(χ² < χ²₀, ν = 2N), where the right-hand side is the probability of observing a value of chi-squared larger than χ²₀ for ν degrees of freedom, χ²₀ = −2ln z = 2ln ∏_{i = 1}^NP_i, and N is the number of observations (N = 2 in our case).

Finally, we account for the fact that the largest interval in the first observation and the largest interval in the second observation happen to overlap in orbital phase space by a significant margin. For two randomly chosen phase intervals Δϕ' and Δϕ'', the probability of overlap by at least an amount ϕ_o is

$$\begin{equation} P_o=\Delta \phi ^{\prime }+\Delta \phi ^{\prime \prime }-2\phi _o{,} \end{equation} \tag{ 1 }$$

for 0 < ϕ_o < min(Δϕ', Δϕ''). In our case, Δϕ' = 0.352, Δϕ'' = 0.223, and ϕ_o = 0.177 (see Figure 2), leading to the overlap probability P_o = 22.1%.

The combined probability of the two gaps being at least as large as observed and overlapping is

$$\begin{equation} P= P_{\rm comb}P_o = 0.22\times 0.043 = 0.0094. \end{equation} \tag{ 2 }$$

Thus, with a 99.1% probability (2.6σ significance), the observed two large, overlapping gaps suggest a physical dependence of the count rate on the orbital phase rather than a random fluctuation in the count rate of a steady source. We confirmed the validity of the above calculations by extensive Monte Carlo simulations.

Finally, we used an alternative method based on the binned light curve (shown in Figure 3) to estimate the probability of the observed count deficit. We produced the phase-folded light curve by binning source counts from both observations. Note that we chose the bin with the deficit to be centered precisely on the physically significant phase of apastron, ϕ_A = 0.5, rather than picking a special phase which minimizes the number of counts in the bin. The Poissonian probability of a bin containing only three or less photons when expecting 11 is 0.49%. We can also evaluate the likelihood of such a light curve being a realization of a constant-rate model, with the expectation number of counts in each bin e_i = 33/N for N phase bins. Since the numbers of counts per bin are small, the commonly used χ² statistic (appropriate for large numbers of counts) is biased. Instead we can use the modified Cash (1979) statistic

$$\begin{equation} C = 2\sum _{i=1}^N n_i \times \left[ \frac{e_i - n_i}{n_i} - \ln \left(1 + \frac{e_i - n_i}{n_i}\right) \right], \end{equation} \tag{ 3 }$$

where n_i is the measured number of counts in the ith bin. For three bins, we find C = 10.95. To test how likely it is to find a higher value of C by chance, we generate 10 000 sets of photon arrival phases, drawn from a uniform distribution in the 0–1 phase range. For each realization, we generate binned light curve as above and calculate C. We find that C > 10.95 in 0.4% of cases. Thus, the null hypothesis of a count rate independent of phase is ruled out with 99.6% confidence. The apparently higher significance reflects our choice of the number of bins in the light curve.

To visualize the placement of the observed deficit in relation to the NS orbits and the observer's line of sight, we plotted in Figure 4 the configurations of the orbits and marked the true anomaly values corresponding to the photon arrival phases. The relativistic advance of periastron by $\delta \omega =\dot{\omega }T=1\mbox{$.\!\!^\circ $}756$ yr⁻¹ × 4.465 yr =784 is accounted for.

Zoom In Zoom Out Reset image size

2.2. Spectrum

We extracted the spectra and produced response files following standard ACIS data reduction procedure. Spectral analysis was performed with Sherpa⁷ in the 0.3–8 keV range. Due to the low numbers of photons, we fitted the models by minimizing the C-statistic (Cash 1979) without binning. We fitted a power-law (PL) model, modified by interstellar absorption (equivalent to wabs*powerlaw in XSPEC), and fixed the absorbing column to N_H = 3.6 × 10²⁰ cm⁻². This value corresponds to the pulsar's dispersion measure for 10% ionization and it is close to the H i column density through the Galaxy in this direction. Note that such a low extinction has virtually no effect on the spectrum above 0.3 keV. We find that, between the two data sets, the best-fit parameters are consistent within their 1σ uncertainties (see Table 1). The combined fit results are shown in Table 1, and the confidence contours corresponding to these parameters are plotted in Figure 5. The unabsorbed luminosity is L_{0.3–8 keV} = (4.0 ± 1.0) × 10²⁹ erg s⁻¹, for the distance of 1.02 kpc.

Zoom In Zoom Out Reset image size

Table 1. Spectral Fits

Parameter	Olda	New	Joint
Γ	2.9 ± 0.5	3.3±0.6	3.07+0.37− 0.36
${\cal N}$b	6.4+1.8− 1.5	5.7+1.8− 1.5	6.1+1.2− 1.1
NH (1020 cm−2)	3.6	3.6	3.6
F0.3–8 keVc	3.2 ± 0.9	3.3 ± 1.1	3.2 ± 0.8
Credd	0.14	0.10	0.12

Notes. All fits are performed in the 0.3–8 keV energy range, the background ignored, and C-statistics are used without binning. N_H is held fixed. ^aFit slightly different from KPG06, due to the newer calibration used and an extra photon included. ^bNormalization in units of 10⁻⁷ photons keV⁻¹ cm⁻²s⁻¹ at 1 keV. ^cFlux in units of 10⁻¹⁵ erg s⁻¹ cm⁻², in 0.3–8 keV range. ^dC-statistic (Cash 1979) per detector energy channel.

Download table as:  ASCIITypeset image

We also jointly fitted an absorbed blackbody (BB) model and obtained a slightly worse fit (C = 154 versus C = 127 for the PL). The best-fit BB parameters, kT = 262⁺³⁰_{− 18} eV and R = 20⁺¹⁵_{− 9} m, are at the hotter/smaller end of the confidence regions calculated by KPG06.

We cannot exclude the possibility that the spectrum is described by a multicomponent model, but the statistics are too poor to analyze such models.

KPG06 suggested three possible scenarios to explain the X-ray emission from J1537: (1) the emission is generated in the magnetosphere or heated polar caps of the recycled pulsar A; (2) the emission results from interaction of the pulsar A's wind with the NS companion, such that the strength of the interaction depends on the varying distance between the two NSs in the eccentric orbits; and (3) emission only happens when the companion NS B passes through the equatorial outflow of the A's wind. Model (1) predicts no dependence on orbital phase, (2) predicts fainter emission at apastron and brighter emission at periastron, and (3) predicts two segments with faint/no emission at opposite phase intervals, around apastron and periastron. Note that modulation would be possible for scenario (1) via absorption in the magnetosphere of NS B; however with the small magnetosphere of NS B and the inclination of the orbit, it is impossible for this system. We find a deficit of photons around apastron and no evidence for enhanced emission at periastron. Thus, our data appear to support model (3). As the pulsar's equatorial plane and the orbital plane are misaligned (by δ = 25° ± 4°; Stairs et al. 2004), the NS B will intersect the wind for only part of its orbit, consisting of two segments (see Figure 4). One expects more emission as the NS B plunges through the pulsar wind and less when it is outside the wind. Following the calculation in KPG06, the current results support a limit on the equatorial outflow opening angle of about Δ ⩾ 20°. X-rays can be generated in the intrabinary shock or can originate from the J1537B's magnetosphere or surface, bombarded by the particles of J1537A's wind.

We should note that the phases of superior conjunction of pulsar A (slightly different in the two observations because of the relativistic advance of periastron) are also within the Δϕ ≈ 0.18 phase interval in which no photons were detected in the two observations (see Figures 2 and 4). However, since an excess of flux around superior conjunction would be more natural than a deficit, due to the Doppler boost, we believe that the observed deficit is connected with NS B being outside of the pulsar A's wind around apastron and is not related to superior conjunction.

In the light of this interpretation, it is illustrative to compare the X-ray properties of J1537 and J0737, the only other DNSB detected in X-rays. For J0737, the pulsed flux from the recycled pulsar J0737A accounts for ≈76% of the total DNSB flux in the energy range 0.15–4 keV, the pulsed fraction growing with energy (Pellizzoni et al. 2008). There is also some evidence for pulsations from the slower, non-recycled pulsar J0737B, but only for small (Δϕ = 0.25) part of the orbit, and the pulsed fraction is rather uncertain. Measuring pulsations with the J1537A's period in the J1537 X-ray emission would provide a lower limit on the J1537A's contribution to the total observed X-ray flux, but the low time resolution of the Chandra ACIS observations does not allow such a measurement.

In contrast with J1537, no orbital phase dependence is seen in the X-ray emission of J0737, down to the limit of ≲ 15% on pulsed fraction (Pellizzoni et al. 2008). It itself does not exclude a sizeable fraction of the J0737 luminosity to be powered by the J0737A's wind. Indeed, both the eccentricity (e = 0.0877; Possenti et al. 2008) and the angle between the orbital angular momentum and the A's spin (<14° at 95% confidence; Ferdman et al. 2008) are much smaller than those in J1537. As a result, the distance between the pulsars changes very little, and pulsar J0737B is always near the equatorial plane of J0737A, i.e., within the wind of J0737A, even if the opening angle of the equatorial outflow is small. One would not, therefore, expect the emission caused by the interaction of the J0737A's wind with J0737B to vary significantly with orbital phase (except, possibly, for variations caused by the changing view of the shock front).

The spectra of both J1537 and J0737 are rather soft in comparison with a typical pulsar wind nebula (PWN; Kargaltsev & Pavlov 2008), and their slopes, Γ ≈ 3.1 and Γ ≈ 3.4 (for the phase-averaged spectrum⁸; Pellizzoni et al. 2008), respectively, are at the upper end of the Γ ≃ 2–3 range typical for recycled and non-recycled old pulsars whose spectra can often be approximated by an absorbed PL if the number of counts is small⁹ (e.g., Bogdanov et al. 2010; Pavlov et al. 2009, and references therein). The similarity of the J1537 and J0737 spectra seems to contradict the interpretation of the J1537 X-ray emission as powered by the J1537A's wind because the X-ray luminosities of the two DNSBs are dominated by physically different mechanisms. However, even though most spectra of PWNe around isolated pulsars fit harder PLs (Kargaltsev & Pavlov 2008), we know very little about emission from pulsar winds shocked in the vicinity of their parent pulsar due to the presence of a binary companion. It is possible that spectra of intrabinary PWNe are softer than those of ordinary PWN (e.g., due to a smaller amount of the magnetic field energy being transferred to the particles, which also implies a higher magnetization, σ; see Pellizzoni et al. 2008). We also note that J1537's spectrum is marginally harder than that of J0737, and the large uncertainties leave room for even more appreciable differences. Deeper X-ray observations of J1537 should better constrain its spectrum and orbital dependence.

A significant difference between the two DNSBs is that the X-ray radiative efficiency, $\eta _{\rm 0.3\hbox{--}8\,keV} = L_{0.3\hbox{--}8\,{\rm keV}}/\dot{E}_A \approx 2.2\times 10^{-4}$ for J1537, is a factor of seven lower than for J0737 (with account for the revised distance to J0737, d ≃ 1150 pc; Deller et al. 2009). The X-ray efficiency of J0737, η_{0.3–8 keV} ≈ 1.6 × 10⁻³, is typical for millisecond pulsars detected in X-rays (see Figure 4 in Pavlov et al. 2007), in accordance with the conclusion that the J0737's emission is dominated by pulsar J0737A (Pellizzoni et al. 2008). On the contrary, the J1537's efficiency is lower than those of nearly all other recycled pulsars. If most of the X-ray luminosity in J1537 comes from the interaction of the pulsar A's wind with pulsar B, as follows from our interpretation, the efficiency of J1537A is even lower, making it one of the least efficient recycled pulsars known.

An upper limit on the wind-powered J0737's luminosity, ∼1.7 × 10³⁰ erg s⁻¹ in the 0.3–8 keV band (which includes the J0737B's luminosity), can be estimated as the difference between the total observed luminosity and the luminosity pulsed with the A's period. This limit is a factor of four larger than the total luminosity of J1537, but this fact does not exclude the possibility that the unpulsed luminosity of J0737 is wind-powered because the J0737 orbit is tighter.

In summary, we have confirmed the deficit of X-ray emission in the J1537 DNSB around apastron. This implies that a significant fraction of the X-ray emission is powered by the equatorial wind of the recycled pulsar J1537A, and the emission probably occurs when the companion NS crosses this wind. Our interpretation also implies that the apparent X-ray efficiency of pulsar J1537A is very low compared to other millisecond pulsars.

We are deeply indebted to Ingrid Stairs for providing us with the up-to-date radio ephemeris. This work was supported by Chandra award GO9-0073X and NASA grant NNX09AC84G. The research was also partly supported by the National Science Foundation under grants 0908733 and 0908611. The work by G.G.P. was partly supported by the Ministry of Education and Science of Russian Federation (contract no. 11.G34.31.0001).

Facilities: CXO (ACIS) - Chandra X-ray Observatory satellite, XMM (EPIC) - Newton X-Ray Multimirror Mission satellite

The probability that, for an observation with N counts detected, the largest of the N + 1 time intervals between neighboring events (which include the beginning and the end of the observation) is as large as or larger than τ can be calculated as the difference between unity and the probability of all the intervals being smaller than τ:

$$\begin{eqnarray} &&P(\Delta t\ge \tau \,|\,N;T)=1-N!T^{-N} \int _0^\tau d(\Delta t_1)\ldots \int _0^\tau d(\Delta t_N)\nonumber\\ &&\qquad\quad\times\Theta \left(T-\sum _{i=1}^{N}\Delta t_i\right)\,\Theta \left(\sum _{i=1}^{N}\Delta t_i+\tau -T\right),\nonumber\\ \end{eqnarray} \tag{ A1 }$$

where T is the duration of the observation, Θ(x) is the Heaviside step function (Θ(x > 0) = 1, Θ(x < 0) = 0), and Δt_i is the interval between the (i − 1)th and ith events (i = 0 and i = N + 1 correspond to the beginning and the end of the observation, respectively). Equation (A1) takes into account that, for N fixed, each possible configuration of intervals ought to have the same probability of occurrence. The step functions in Equation (A1) correspond to the following constraint:

$$\begin{equation} 0 < T -\sum _{i=1}^N \Delta t_i = \Delta t_{N+1} < \tau {.} \end{equation} \tag{ A2 }$$

The analytical calculation of the N-dimensional integral in (A1) leads to very cumbersome expressions; in particular, for each given N, they are different for the separate intervals (k + 1)⁻¹ < τ/T < k⁻¹, where k = 1, 2, ..., N. However, the integral can be evaluated numerically.

The number of counts is itself a random quantity that obeys some statistical distribution. For the null hypothesis we are investigating, the probability of detecting n photons in a time interval t is given by the Poissonian distribution

$$\begin{equation} P_{\rm Poiss}(n,\lambda t)=\frac{(\lambda t)^n\, e^{-\lambda t}}{n!}\,, \end{equation} \tag{ A3 }$$

where λ is the mean count rate. For a given λ, the probability that the largest interval between neighboring events is as large as or larger than τ is the expectation value of P(Δt ⩾ τ | N; T):

$$\begin{equation} P(\Delta t\ge \tau \, |\, \lambda T)=\sum _{N=0}^{\infty } P(\Delta t\ge \tau \, |\, N;T)\, P_{\rm {Poiss}}(N,\lambda T).\, \end{equation} \tag{ A4 }$$

For the first observation (τ = 0.352, T = 0.9925, and λ = 16.15; the unit of time is the orbital period) this gives P(Δt > 0.352 | 16.03) ≡ P₁ = 0.0381.

For the second observation, the largest interval in the phase space is not the one between the consecutive events but the sum of the waiting times before the first count and after the last one (see Figure 2). The probability of the sum of two intervals to be larger than τ when the total number of events N in the observation time T is given can be calculated as

$$\begin{eqnarray} &&P(t_1+t_{N+1}\ge \tau \, |\, N;T)=1 - N! T^{-N} \int _0^\tau d(\Delta t_1)\ldots \nonumber \\ &&\qquad\times\int _0^\tau d(\Delta t_N)\Theta\left(T-\sum _{i=1}^{N}\Delta t_i\right)\Theta \left(T-\sum _{i=2}^{N}t_i-\tau \right).\nonumber\\ \end{eqnarray} \tag{ A5 }$$

Since the expectation value of P(t₁ + t_{N + 1}  ⩾  τ | N; T) is P(t₁ + t_{N + 1}  ⩾  τ | λT) = ∑^∞_{N = 0}P(t₁ + t_{N + 1} ⩾ τ | N; T)P_Poiss(N, λT), for the second observation we obtain P(t₁ + t₁₇ ⩾ 0.274 | 16.70) ≡ P₂ = 19.0%.