MULTI-EPOCH ANALYSIS OF PULSE SHAPES FROM THE NEUTRON STAR SAX J1808.4−3658

In this paper, we restrict the models to only one visible spot. Prior to doing this we carried out fits to the data from all seven time periods with a two-spot model. The best-fit two-spot models in all cases had zero amplitude for the second spot, indicating that no second spot is required by the observed pulse shapes. This is consistent with the model introduced by Ibragimov & Poutanen (2009) where the accretion disk hides the antipodal spot on the other side of the star during the early stages of the outburst (stages 3 and 4 in the terminology of Hartman et al. 2008; stages P and SD in the terminology of Ibragimov & Poutanen 2009). In the later periods of the burst (stages 5 and 6) the disk recedes away from the star in the Ibragimov & Poutanen (2009) model, exposing the antipodal spot and creating a more complicated pulse shape. The low flux and signal-to-noise ratio in the pulse shapes corresponding to stages 5 and 6 of the 2002 and 2005 outbursts mean that it is not feasible to use these pulse shapes in the present analysis.

In an earlier paper (Leahy et al. 2009), the light curve of a different neutron star, XTE 1814, was analyzed. In that case, a one-spot model led to a residual (or bump) in the pulse shape that was seen in all wavelengths. In that case, the residual provided the motivation for a two-spot model for the data which did indeed fit the data significantly better. In the case of SAX J1808, the one-spot models (shown in Section 4) do not lead to a significant excess feature at all wavelengths that would motivate a two-spot model.

Poutanen & Gierliński (2003) have shown using analytical approximations that the theoretical pulse shape resulting from a hot spot is independent of the spot size, as long as the spot radius is small compared to the size of the neutron star radius. Previous computations of pulse shapes (e.g., Leahy et al. 2008, 2009) have confirmed that the assumed spot size does not strongly affect the resulting pulse shape. For this reason, we use a simplified computational procedure that treats the spot size as infinitesimally small. This speeds up the computations considerably compared to the use of spots that cover a larger area. Since the distance to the star and the overall intensity normalization are kept arbitrary, the area of the spot is not needed in our models.

The best-fit models using one infinitesimal spot have unphysically low masses. For example, in row 1 of Table 1 we show the best-fit model for the 2002B3 epoch with 2M/R fixed at 0.5. (Results for other values of the gravitational radius are similar.) Since there are 2 energy bands with 32 time bins each and a total of 8 free parameters, there are 56 degrees of freedom. The best-fit model has a chi-squared value of 69.7 for 56 degrees of freedom, and a mass of only 0.3 M_☉ which is unreasonably low for neutron stars, or for quark stars. Increasing the radius of the spot so that it has a radius of 0.1 R or 0.2 R does not change the mass or χ² significantly, as can be seen in rows 2 and 3 of Table 1.

We also tested the hypothesis that a non-circular spot could affect the results. We tested a long, thin uniform intensity spot that had a length of 0.4 R azimuthally. The result for this extreme shaped spot is shown in row 4 of Table 1 where it can be seen that there is again no significant change from the infinitesimal spot. A similar test of a long, thin spot with only latitudinal extent gives similar results. An alternate type of non-circular spot can be approximated by a two-spot model where both spots are artificially forced to have the same intensity, but the relative angular locations are allowed to be free. This would allow for two spots to be close together, giving an effective irregular shape. The result for this two-spot model is shown in row 5 of Table 1. This model marginally improves the χ², but the best-fit mass is still only 0.354 M_☉. We conclude that changes in the spot size and shape do not help the models to converge to physical values for neutron stars.

The hot-spot model described in Section 3 was used to fit two different data sets simultaneously. The resulting fits were found to be very poor. For example, data from 1998B4 and 2002B3 were fit using this model. The case for 2M/R = 0.5 is shown in the first row of Table 2 (labeled "no scatter"). The resulting poor χ² of 200 for 113 degrees of freedom is typical for fits done for other values of 2M/R and for other pairs of data.

The poor fits that result from the application of the hot-spot model motivate us to introduce an extra component to the model that includes light that travels from the hot spot on the star, scatters off of material near the star, and then arrives at the observer. The scattering material could be either an accretion column or an accretion disk. If the light is scattered by the disk, then there will be a bright pattern on the disk that will appear to rotate around the star at the same frequency as the star, as studied by Sazonov & Sunyaev (2001). In this section, we will outline a simple model for scattering off of an optically thin cloud of electrons, such as the ionized surface layer of an accretion disk.

In this simple model, we do not include the effects of gravitational light bending or Doppler boosting for the scattered light in order to include only the most important features. The reason why we can neglect these normally important relativistic effects is that the amplitude of the scattered light is a very small fraction of the amplitude of the light that travels directly from the spot to the observer. We have tested that the inclusion of Doppler boosting to the scattered-light model does not make any significant difference to the outcome of the fits.

Consider a spot at co-latitude θ, azimuthal angle ϕ on the surface of a star of radius R. The unit vector pointing from the center of the star to the spot is

$$\begin{equation} \hat{n} = \sin \theta (\cos \phi \hat{x} + \sin \phi \hat{y}) + \cos \theta \hat{z}. \end{equation} \tag{ 2 }$$

Photons emitted by the spot are scattered by material at co-latitude ϑ and azimuthal angle φ at a distance r. The unit vector pointing to the location of the scatterer is

$$\begin{equation} \hat{r} = \sin \vartheta (\cos \varphi \hat{x} + \sin \varphi \hat{y}) + \cos \vartheta \hat{z}. \end{equation} \tag{ 3 }$$

When light bending is neglected, the photon travels on a straight line defined by the vector $\vec{\ell }= r \hat{r} - R \hat{n}$ which has magnitude

$$\begin{equation} \ell = r (\sin ^2\sigma + (\cos \sigma - R/r)^2)^{1/2}, \end{equation} \tag{ 4 }$$

where σ is the angle between the vectors $\hat{r}$ and $\hat{n}$ defined by

$$\begin{equation} \cos \sigma = \cos \theta \cos \vartheta + \sin \theta \sin \vartheta \cos (\phi -\varphi). \end{equation} \tag{ 5 }$$

The vector pointing toward the observer is $\hat{k} = \sin i \hat{x} + \cos i \hat{z}$, so the photon must be scattered through an angle ρ in order to be detected, where ρ is defined by

$$\begin{eqnarray} \cos \rho = \hat{l} \cdot \hat{k} &=& \frac{1}{\ell } ({-} R(\cos \theta \cos i + \sin \theta \sin i \cos \phi)\nonumber\\ && + r(\cos \vartheta \cos i + \sin \vartheta \sin i \cos \varphi)). \end{eqnarray} \tag{ 6 }$$

The differential cross section for Thompson scattering is $\frac{d\sigma }{d\Omega } \propto (1 + \cos ^2\rho),$ where ρ is the angle that the photon is scattered through. If the flux of light originating from the hot spot at angle ϕ that impacts the scatterer at φ is F(ϕ), the flux of scattered light reaching the observer is then F_sc ∝ F(ϕ)(1 + cos²ρ) where the constant of proportionality depends on the optical depth of the scattering material.

In order to model scattering of light by an optically thin cloud of electrons close to the accretion disk, the co-latitude of the scattering material is set to ϑ = π/2. Photons emitted by the spot at angle ϕ = 2πt/P + ϕ₀ will impact the disk over a range of azimuthal angles φ and radii r. The pattern appearing on the disk will be centered around an angle φ equal to ϕ. However, the photons emitted by the spot and then scattered by the disk reach the observer at a later time than the photons that travel directly to the observer without scattering. A simple way to treat the time lag is to define an effective phase lag Δϕ so that φ = ϕ + Δϕ, where Δϕ includes the delays due to the light propagation time from the star to the disk as well as the light-crossing times for the photons to cross the disk. The lowest order expression for the phase lag (ignoring relativistic corrections) is

$$\begin{equation} \Delta \phi = 2 \pi \left(\frac{r (\cos ^2 \theta + (\sin \theta - R/r)^2)^{1/2}}{c P} - \frac{r}{c P} \sin i \cos \phi\right), \end{equation} \tag{ 7 }$$

where P is the neutron star's spin period. The time lags for the scattered light should be much less than the star's spin period. Since cP = 750 km for SAX J1808.4−3658, physical solutions should have a scattering radius r ≪ cP. In this model, we are approximating the extended pattern on the disk as an effective point source. We tested more complicated models with a smeared out pattern on disk and found that the changes from the point-source model were insignificant.

The minimum radius of impact on the disk for a photon emitted at co-latitude θ is r_min = R/sin θ. Similarly, parts of the disk at radii r>r_shadow can be observed at all values of φ, where the shadow radius is r_shadow = R/cos i. As a result, the pattern on the disk is always visible to the observer if θ ⩽ π/2 − i. If the spot is at a high latitude (i.e., θ ≪ 1) as suggested by Lamb et al. (2009), then the illuminated part of the disk will be at r ≫ R and the illuminated part of the disk will not be eclipsed by the star.

As the spot moves around the star, the resulting pattern on the disk appears to move around the star at the same rate. In the frame rotating with the star, the pattern will not appear to move, so the flux F(ϕ) impacting the disk will be a constant. As a result, the scattered flux observed will be proportional to just the differential cross section. In addition to the 14 parameters required for the hot-spot model introduced in the previous section, the scattered-light model requires that we add two more free parameters: I_sc and r. The scattered light is added to the model by adding the term

$$\begin{equation} F_{{\rm sc}} = I_{{\rm sc}} (1 + \cos ^2\rho) \end{equation} \tag{ 8 }$$

to each energy band. Since we are implementing a simple model for scattering from a disk (ϑ = π/2), the scattered flux is

$$\begin{eqnarray} F_{{\rm sc}} &=& I_{{\rm sc}} \bigg(1 + \bigg(\frac{r}{\ell }\bigg)^2 \bigg(\sin i \cos (\phi +\Delta \phi)\nonumber\\ && - \frac{R}{r}(\cos \theta \cos i + \sin \theta \sin i \cos \phi) \bigg)^2\bigg), \end{eqnarray} \tag{ 9 }$$

where Δϕ depends on the free parameters r, i,  and θ through Equations (4) and (7). The normalization for the hot-spot model is defined so that the amplitudes I₁ and I₂ are of order unity, so a physical scattered-light model should have an amplitude I_sc ≪ 1.

In row 2 of Table 2, the results of a joint fit to the 1998B4 and 2002B3 data sets using the hot-spot model along with the simple scattering model given by Equation (9) are shown (labeled "Model 1"). Three new parameters were introduced: the flux amplitudes I_s1andI_s2 for the scattered light in 1998B4 and 2002B3, respectively, and r, the effective radius at which the scattering mainly occurs. The addition of these three parameters to the model brings the χ² value down to 120, which is an acceptable fit for 110 degrees of freedom. The resulting best-fit stellar model has a mass of 1.44 M_☉. The best-fit spot latitudes are close to the north pole in both time periods, while the inclination angle is large, placing the observer close to the equatorial plane of the star. The best-fit scattering radius of 116 km satisfies both inequalities r>R/sin θ_1,2 and r>R/cos i so that this location will be illuminated by the spot and will always be visible to the observer. The relative amplitude of the scattered radiation in 1998B4 is 0.2% of the flux coming directly from the spot. Similarly, the relative amplitude for the 2002B3 data is 0.8%. Both amplitudes are very small fractions of the direct flux, as would be expected for scattered light. The minimum in χ² is fairly broad in the parameters corresponding the mass of the star and inclination, with the 1σ allowed region including a range of masses between 1.2 and 1.6 M_☉. Only one model (for 2M/R = 0.5) is shown in Table 2, but the results are similar for other values of 2M/R.

The best-fit solutions using the scattering model described above all have the property that the location of the scatterer is large compared to the radius of the star, r ≫ R. This motivates the approximation of Equation (9) to lowest order in R/r:

$$\begin{equation} F_{{\rm sc}} = I_{{\rm sc}} (1 + \sin ^2 i \cos ^2(\phi +\Delta \phi)). \end{equation} \tag{ 10 }$$

Furthermore, the inclination angles are fairly close to 90°, so sin i ≃ 1. Using this approximation leads to

$$\begin{equation} F_{{\rm sc}} = I_{{\rm sc}} (1 + \cos ^2(\phi +\Delta \phi)). \end{equation} \tag{ 11 }$$

If the best-fit value of r/R ≫ 1 and sin i ∼ 1, Equation (11) should be a good approximation to Equation (9).

As an alternative to the scattered-light model given by Equation (9), we now test the simpler model given by Equation (11), with Δϕ as a free parameter, instead of r. The best-fit stellar model that results when Equation (11) is added to the basic hot-spot model and fit to the 1998B4 and 2002B3 data sets is shown in row 3 of Table 2 (labeled "Model 2"). The best-fit star (with 2M/R = 0.5) using Model 2 has χ² = 120 for 110 degrees of freedom, so it is as good a fit as Model 1. The best-fit values of the various parameters are different from the values for Model 1, however, the differences are smaller than the 1σ differences allowed by either model. For instance, the mass using Model 2 is 1.57 M_☉, which lies in the 1σ region allowed by Model 1. We conclude that the two models produce results that are statistically the same. For this reason, we use the simpler model using Equation (11) to model the scattered light for the rest of this study.

Since the addition of three parameters (I_s1,2 and Δϕ) brings about a significant reduction in χ², and the masses are raised to realistic values, we include this scattered-light model in subsequent fitting of the pulse shapes. In summary, this component is required in order to obtain acceptable fits and to obtain physically reasonable masses.

A small constant (or "dc") flux component was required in a similar analysis of XTE J1814 (see Section 3.2 of Leahy et al. 2009). In the case of XTE J1814, two free parameters (one for each energy band) were included to improve the quality of the fits. The amplitudes of the dc components were less than 2% of the pulsed component, similar to what is seen in the present analysis of SAX J1808.

The joint analysis of the 2002B3 data with the 1998B4 data results in best-fit neutron star models with larger masses compared to the best-fit models that only use the 1998B4 data. Table 3 shows the values of the best-fit parameters for a range of neutron star compactness values (2M/R). Each row of Table 3 is generated by keeping the compactness ratio 2M/R (at the location of the spot) at a fixed value and varying all other parameters until the lowest value of χ² is found. For each value of 2M/R the following variables are displayed: the mass of the neutron star M, the equatorial radius of the neutron star R, the angular location of the spot θ₁ (as measured from the north pole) in 1998, the angular location of the spot θ₂ in 2002, the inclination angle i, and the anisotropy parameters for 1998 (a₁) and 2002 (a₂). We assume that the orbital and spin axes are aligned, so i is also the binary's inclination angle. In addition, there are four normalizations, two phases, two blackbody-to-Compton ratios, and four scattered-light parameters which are not displayed in the table, since they are of no interest here. A χ² penalty was used to keep the blackbody-to-Compton ratios b₁ and b₂ within 1σ of the spectral models.

The best-fit models in Table 3 should be compared with Table 5 of Leahy et al. (2008) where only data for 1998 were used. Consider the best-fit models with 2M/R = 0.4 in both tables. When only 1998 data were used, the best-fit model had small mass and radius (M = 0.90 M_☉ and R = 6.7 km) and both the spot's angular location and the inclination angle were close to 30°. The addition of data from 2002B3 allows for a larger mass and radius (M = 1.43 M_☉ and R = 10.8 km) for the same value of compactness. The larger radius is compensated for by moving the spot closer to the spin axis with θ₁ close to 10° in 1998 and 5° in 2002B3.

The best-fit values of the scattered-light amplitudes A₁ are 2 × 10⁻³ for 1998B4 and 8 × 10⁻³ for 2002B3, where the flux is normalized to 1.

The best-fit model with 2M/R = 0.6 in Table 3 is shown in Figures 1 and 2 with dashed curves for the low-energy bands and solid curves for the high-energy bands. All four light curves are fit simultaneously, but for clarity, the data and best-fit models are plotted in two separate figures. The relative phases between the high- and low-energy bands for each data set are included in the models.

Figure 7 shows 2σ and 3σ confidence contours on the mass–radius plane for joint fits with 1998B4 and 2002B3. This figure is generated by specifying a mass and then varying all parameters until the lowest value of χ² for the mass is found. Two constraints were added to the variation of parameters. First, the ratio 2M/R was only sampled in the range 0.2 ⩽ 2M/R ⩽ 0.6 in order to stay within the bounds of reasonable EOSs. Second, an unconstrained variation of the inclination angle leads to acceptable fits for solutions with i>90°, which is very unlikely. Hence, we constrained the inclination angle to values with i ⩽ 90°. The resulting confidence contours in Figure 7 should be compared with the corresponding confidence contours for the 1998 data shown with bold curves in Figure 3 of Leahy et al. (2008) (labeled "oblate and time delays"). The 3σ confidence region only allows very small stars with R < 7 km if only the 1998 epoch is used. The additional data from 2002B3 and allowance for scattered light used in the present analysis allow stars with larger radii, extending to 13 km at the 3σ confidence level. The largest mass stars allowed (at 3σ) have M = 1.7 M_☉. At 3σ, the data allow for a wide variety of modern EOSs.

Zoom In Zoom Out Reset image size

The best-fit models displayed in Table 3 have inclinations ranging from 70° to 90°. However, for each best-fit set of values for mass and radius shown in Table 3, there is a range of inclination angles allowed by the data. The range of inclination angles allowed at the 3σ level is found by keeping the mass and radius fixed at the values shown in the table, and allowing all other parameters to vary. In Figure 8, the locations on the mass–radius plane for the five models shown in Table 3 are marked with triangles. The allowed 3σ range of inclination angles (in degrees) for each of these points is shown in Figure 8. For example, the model with M = 1.43 M_☉ and R = 10.8 km is allowed inclination angles in the range of 73° < i < 81°, although the lowest χ² results for i = 755.

Zoom In Zoom Out Reset image size

In Figure 8, the same 3σ contours on the mass–radius plane that were shown in Figure 7 are redisplayed. For selected points on the contours (shown with circles), the one value of inclination that is allowed is shown. For example, Figure 8 shows that the data allow at the 3σ confidence level a star with M = 1.44 M_☉, R = 12.5 km, and i = 79°. In general, the lowest mass stars require smaller inclination angles (near 40°), while the highest mass stars require inclinations that are as high as is allowed.

There are a number of observations that suggest that the inclination angle may be smaller than 90°. The strongest evidence is the absence of eclipses, which constrains i < 82° assuming a Roche lobe filling companion (Chakrabarty & Morgan 1998). Wang et al. (2001) observed V4580 Sgr, the optical counterpart of SAX J1808, during the 1998 outburst. They modeled the optical and IR data using an X-ray-heated accretion disk and found limits on the inclination of 42° < i < 88° if the distance is d = 2 kpc and 20° < i < 65° if d = 3 kpc. Further modeling of the optical data observed during quiescence by Wang et al. (2009) suggests that i < 70° using the distance of d = 3.5 kpc derived by Galloway & Cumming (2006). Deloye et al. (2008) modeled the optical data from quiescence using a model that includes emission from both the disk and the companion. Using a distance of 3.5 kpc, Deloye et al. (2008) constrain the inclination to 32° < i < 74° if 10% uncertainty in the distance is assumed.

A feature interpreted as an iron line relativistically broadened through the orbital motion of the disk was detected in X-ray observations during the 2008 outburst by Cackett et al. (2009) and Papitto et al. (2009). Both groups modeled the emission feature and constrained the inclination angle. Cackett et al. (2009) found i = 55⁺⁸₋₄ deg, while Papitto et al. (2009) provide a less precise constraint of approximately i>60°. A different type of constraint on the inclination arises from the analysis of the 2002 outburst data by Ibragimov & Poutanen (2009). In their analysis, Ibragimov & Poutanen (2009) assume that the X-ray emission is produced by two antipodal spots, one of which is sometimes hidden by the accretion disk, due to movement of the inner edge of the accretion disk. In this model, they find that the data can be explained by inclinations in the range i < 73°.

Although all of the constraints on inclination have some model dependence, there is a general agreement that inclinations smaller than about 70° are suggested by the different data. With this in mind, we refit the data for 1998B4 and 2002B3 subject to the constraint i ⩽ 70°. The best-fit models for fixed values of 2M/R are shown in Table 4. Comparing Table 4 with Table 3, it can be seen the most significant change is that lower mass stars are allowed when the inclination angle is constrained to i ⩽ 70°. The values of χ² increase, but the magnitude of the increase is not significant.

Figure 9 shows the 2σ and 3σ confidence contour on the mass–radius plane using the constraint i ⩽ 70° for joint fits to the 1998B4 and 2002B3 data. The 3σ confidence region shrinks when the constraint on inclination is added; however, stars as massive as 1.5 M_☉ are still allowed.

Zoom In Zoom Out Reset image size

In all of the hot-spot models for the 1998B4 and 2002B3 data sets, the co-latitude of the spot θ is larger (further from the pole) in 1998 compared to the co-latitude in 2002. This trend is seen in all models fitting the other data sets as well. Since a hot spot located close to the spin axis produces a less-modulated pulse shape, this is a natural outcome of the fact that the 1998 data are highly modulated (∼10%) and the data for the later outbursts have smaller modulation (∼5%).

In the best-fit models for the 1998B4 and 2002B3 data shown in Table 3, the shift in co-latitude Δθ = θ₁ − θ₂ ranges from 5° to 8°. When all models that fit the data within 3σ confidence are considered, the co-latitude in 1998 lies in the range 9° < θ₁ < 32° and that in 2002B3 lies in the range 4° < θ₂ < 17°. The 3σ range of co-latitude difference is in the range 6° < Δθ < 15°. This wandering of the hot-spot location is consistent with the model proposed by Lamb et al. (2009) where a hot spot close to the spin axis moves due to movement of the star's magnetic field.

The 2002B4 and 1998B4 data have similar pulse shapes (although a statistical test shows that they are different at high significance—chance probability <10⁻²⁰). Of all the data sets, these two epochs have the smallest error bars which are roughly equal in magnitude. The most important difference is in the modulation of the pulse shapes.

The best-fit models for joint fits to the 1998B4 and 2002B4 data sets are shown in Table 5. In this table, the subscript "2" corresponds to values in the 2002B4 epoch, while the subscript "1" corresponds to values in the 1998B4 epoch. These best-fit models should be compared with those shown in Table 3. The range of masses and radii found when the 2002B4 data are fit is similar to those found when data from 2002B3 are used. The values of mass and radius shown in Table 5 are consistent with the 3σ limits arising from the joint fit between the 1998B4 and 2002B3 data, with the exception of the model with a mass of 0.66 M_☉. The best-fit values of the 1998 angular location of the spot θ₁ and the inclination angle i differ in the best-fit models displayed in Tables 3 and 5. This is not surprising, since as we have shown earlier, given values of mass and radius, a wide range of angles θ and i will give acceptable fits to the data.

The best-fit model's light curve for the 2002B4 data is shown in Figure 3. (The light curves for 1998B4 generated by the best-fit model are not distinguishable from the curves shown in Figure 1 and are not plotted.)

The 2005 outburst was only about 70% as luminous as the 2002 outburst (Hartman et al. 2008). In addition, the RXTE satellite has lost sensitivity since 2002. As a result, the 2005 data are much noisier and have larger error bars than the 1998 and 2002 data as can be seen in Figures 4–6. One-spot models were used to simultaneously fit the 2005 data with the 1998 data. The best-fit model parameters are shown in Tables 6–8 and the best-fit light curves are shown in Figures 4–6.

The best-fit models that make use of the 2005 outburst data shown in Tables 6–8 are very similar to the best-fit models computed using the 2002B4 data and displayed in Table 5. The models with 2M/R = 0.4, 0.5, and 0.6 (and with M>1.0 M_☉) are consistent with the 3σ limits arising from the joint 1998B4–2002B3 fits.

The 2008 outburst was almost as bright as the 2005 outburst; however, the degradation of the RXTE Proportional Counter Array over the intervening three years has the result that the pulse shapes from the 2008 outburst (Hartman et al. 2009), while similar to the 2005 pulse shapes, are even noisier than the 2005 pulse profiles. For this reason, we have not attempted to fit the 2008 pulse profiles.