CONSTRAINTS ON THE PROPERTIES OF THE NEUTRON STAR XTE J1814−338 FROM PULSE-SHAPE MODELS

Pulse shapes for the accretion-powered pulsations are constructed using the ephemeris given by Chung (2007). Data are limited to the first 23 days of the 2003 outburst, 2003 June 5–27 (MJD 2452795–2452817), in order to avoid the later period of the outburst when the flux and pulse shape became more erratic (see, for example, Watts et al. 2005). X-ray bursts were cut out of the data during the interval between 100 s before and after the start of each burst.

Although the Rossi X-Ray Timing Explorer (RXTE) observations include data in the range 2–50 keV, we have chosen to concentrate on the lower energy range from 2–10 keV for two reasons. First, the data are noisier at energies above 10 keV. Second, the Chandra observations by Krauss et al. (2005) constrain the spectrum in the 2–10 keV range. It is also useful to separate the data into narrow energy bands in order to separate the different spectral components. We have chosen two narrow bands, the 2–3 keV band and the 7–9 keV band based on the Chandra spectrum. The narrowband pulse shapes constructed using data from the full observation period (June 5–27) are shown in Figure 1.

Zoom In Zoom Out Reset image size

We have also investigated the variability of the pulse shape with time. In order to do this, the data were separated into one-day segments and separate light curves constructed for each day. It is not computationally feasible to model all days simultaneously, so we instead focus on two separate days. The days were chosen by comparing light curves in the 2–10 keV range for different days using a χ² test and selecting two days which differ the most from each other. This also has the effect of selecting days with intrinsically smaller error bars. The days resulting from this selection process correspond to June 20 and 27. Light curves for these two days, in the two narrow energy bands are shown in Figure 2.

Zoom In Zoom Out Reset image size

Krauss et al. (2005) observed XTE J1814 with Chandra on 2003 June 20. They modeled the spectrum in the 0.5–10 keV range and found that the best-fit solution corresponds to a combination of a 0.95 keV blackbody and power-law emission with a photon spectral index of Γ = 1.4. The ratio of flux from the blackbody to the power law in their model is about 10%. We use the Krauss et al. (2005) spectral model and assume that it holds for the other days covered by the RXTE data. This assumption is motivated by the fact that the relative normalization of different energy bands is approximately constant from day to day, although the overall flux at all wavelengths changes with time.

The spectral model of Krauss et al. (2005) motivates the use of two narrow bands in our pulse-shape models. A low-energy band is necessary in order to capture the blackbody component, so we choose the lowest possible X-Ray Timing Explorer (XTE) energy band at 2–3 keV. The spectrum in this band is dominated by the power-law component, but the blackbody contribution is still important. We also choose the 7–9 keV band as the highest energy band covered by the Chandra observation. In this high-energy band the blackbody flux is negligible.

Our method for modeling the observed emission is very similar to the method presented in Leahy et al. (2008). The spectral model has three components: (1) Comptonized flux in the high-energy band (7–9 keV); (2) Comptonized flux in the low-energy band (2–3 keV); and (3) blackbody flux in the low-energy band. The observed flux for the ith component, F_i, integrated over the appropriate observed energy band is given by

$$\begin{equation} F_i(E) = I_i \eta ^{3+\Gamma _i} (1 - a_i \mu). \end{equation} \tag{ 1 }$$

In Equation (1), I_i is a constant amplitude, η is the Doppler boost factor, Γ_i is the photon spectral index in the star's rest frame, μ is the cosine of the angle between the normal to the star's surface and the initial photon direction, and the constant a_i describes the anisotropy of the emitted light. For a definition of η as well as a more complete description of the modeling method, please see Leahy et al. (2008).

In our modeling, the amplitudes I₁ and I₂ are free parameters while the third amplitude I₃ is defined through the constant $b = \bar{F_3}/\bar{F_2}$, the ratio of the phase-averaged blackbody to Comptonized flux in the low-energy band. In the spectral model by Krauss et al. (2005), b = 0.1, but we include this parameter as a fitting parameter with 1σ limits from their spectral model. The photon spectral indices for the Comptonized components are fixed at Γ₁ = Γ₂ = 1.4 as given by their model. In the narrow range of the low-energy band the blackbody component of 0.95 keV can be modeled by a power law with photon spectral index Γ₃ = 0.85. The anisotropy parameters for the Comptonized components a₁ = a₂ = a are assumed to be equal, and the parameter a is kept as a free parameter.

In the modeling of the nonaccreting ms pulsars, it was found (Zavlin & Pavlov 1998; Bogdanov et al. 2007, 2008) that a limb-darkened hydrogen atmosphere spectral model is required by the data. It is reasonable to expect that the blackbody component of the spectrum should also be limb-darkened. We tested this hypothesis by multiplying the blackbody flux by a limb-darkening function of the form e^−τ/μ. We then computed the best-fit neutron star models for two types of models: (1) models with nonzero optical depth τ and (2) models with zero optical depth. The best-fit models for these two cases are almost identical: the mass and radius of the best-fit model changes by less than 0.5% when a nonzero optical depth is added, and the value of δχ² = 0.1 when the limb-darkening is added. Since the change in χ² and the physical parameters are negligible we conclude that adding an extra parameter to model limb-darkening is not warranted by the data. The reason for this is due to the Chandra model which restricts the blackbody contribution in the 2–3 keV band to only 10% of the Comptonized contribution, and effectively sets the blackbody component to zero in the high-energy band. Since the Comptonized flux is dominant and has fan-beaming included, small changes to the anisotropy of the blackbody component do not affect the final models. For this reason we have set the anisotropy parameter for the blackbody component to zero (a₃ = 0). This is consistent with the results found for SAX J1808.4−3658 (Leahy et al. 2008) which also did not require any limb-darkening. The final set of free parameters describing the spectrum are I₁, I₂, b, and a.

In order to fit a set of light curves we also need to introduce a set of parameters describing the star and the emission geometry. These parameters are the mass M and equatorial radius R of the star, the co-latitude of the spot θ, the inclination angle i as well as a free phase ϕ. The radius of the circular spot (in the star's rest frame) is kept fixed at 1.5 km, as given by the Chandra spectral model (Krauss et al. 2005).

Our models make use of light curves for two different days' data, which requires a separate set of parameters for each day. However, on the two different days the parameters M, R, and i do not change. In order to simplify the analysis, we also assume that the photon spectral indices and the parameters a and b are also fixed. The full set of free parameters are {I₁, I₂, θ, ϕ} for each day plus M, R, i, a, and b, for a total of 13 free parameters. However, for each of our fits, the ratio M/R is kept fixed, so for any one value of M/R, there are only 12 free parameters.

We use the oblate Schwarzschild approximation (Morsink et al. 2007) to connect photons emitted at the star's surface with those detected by the observer. In previous studies (Cadeau et al. 2005, 2007) we have shown that, to the accuracy required for extracting the parameters of a rapidly rotating neutron star, it is sufficient to use the Schwarzschild metric to compute the bending of light rays and the relative time delays of photons emitted at different locations on the star. The extra time delays and light bending caused by frame-dragging or higher-order rotational corrections in the metric are negligible. However, the rotation of the star causes a deformation of the star into an oblate shape, which changes (relative to a sphere) the directions that photons can be emitted into. We have developed a simple approximation (Morsink et al. 2007) that allows an empirical fit to the oblate shape of a rotating star to be embedded in the Schwarzschild metric and make use of it in this analysis.

The pulse profiles in Figure 1 show a feature in the phase interval between 0.24 and 0.4. This feature is seen in all of the other energy bands as well. In order to investigate the nature of this feature, we restrict the analysis here to just the 7–9 keV light curve. Since the blackbody contribution in this energy band is negligible, we use a simplified model which only includes the Comptonized component of the radiation.

The simplest model for the emission is a single spot. We fitted the 7–9 keV light curve shown in Figure 3 with a single-spot model by first fixing 2M/R = 0.4 and varying the parameters M, i, θ, a, I, and ϕ (similar results are obtained for other values of 2M/R). Since there are 32 data points this corresponds to 25 degrees of freedom (dof). The best-fit solution for a single-spot model has χ² = 50.7, which is not a very good fit. This best-fit solution is shown as a solid curve in Figure 3.

Zoom In Zoom Out Reset image size

We now turn to a two-spot model, where the second spot is allowed to have an arbitrary location relative to the first spot, but the spectrum is assumed to be the same. The introduction of a second spot introduces three new parameters to the model: an intensity and two angles. The resulting fit for 2M/R = 0.4 has χ² = 28.3 for 23 dof, which is a significant improvement. The light curve for this model is shown with a dotted curve in Figure 3. The mass for this model is 2.04 M_☉. In this model, the second spot's location is situated so that the second spot is almost never seen and its light only contributes during the phase interval between 0.24 and 0.4. Outside of this interval the light received only originates from the primary spot. This suggests a simpler one-spot model where bins in the phase interval 0.24–0.4 (marked with a circle in Figure 3) are removed from the data set. This reduces the dof to 19 (32 − 6 data points and 5 parameters for a one-spot model). The resulting best-fit model for 2M/R = 0.4 has χ² = 19.2, which is also a significant improvement from the one-spot model that uses all of the data. The mass for this model is 2.08 M_☉ and the light curve for this model is shown as a dashed line in Figure 3.

Comparing the two-spot model with the one-spot model with the circle-bins excluded, we see that the difference in χ² is not significant, and there is little change between the best-fit values of mass and radius. This leads us to the conclusion that there is good evidence for a second spot (or a feature that mimics a second spot), but that the amount of data encoded in those bins affected by the second spot is not sufficient to allow us to model the details of the second spot with any confidence. Since the inclusion or exclusion of the second spot does not change the best-fit values of the main physical parameters of the star (M, R, i) it is more appropriate to choose the simpler one-spot model. The qualitative results do not change when we look at different energy bands. As a result, for the remaining modeling reported in this paper we only use one-spot models where data in the 0.24–0.4 phase range are removed from the analysis.

Our procedure for modeling the data from the two energy bands shown in Figure 1 is to assume a one-spot model that includes both blackbody and Comptonized emission. Data in the phase period between 0.24 and 0.4 are omitted, as described in Section 3.1. We do a number of fits, each with a fixed value of 2M/R (fixing the ratio of M/R simplifies the fitting procedure since the light-bending and time delays depend on this ratio). Once 2M/R is fixed, all other parameters are allowed to vary and the minimum value of χ² is found. In addition to the parameters described in Section 2.2, we also added two parameters corresponding to DC offsets for the two energy bands. This allows for small errors in the background subtraction. Once 2M/R is fixed, we have a total of 10 free parameters: M, θ, i, a, b, 2 amplitudes, 2 DC offsets, and 1 overall phase. The parameter b is restricted to have a value that is within 1σ of the value found by Krauss et al. (2005). Since each energy band has 32 points, but we exclude 6 of these points we have 64 − 12 − 10 = 42 dof.

For a fixed value of 2M/R, we find that there are two local minima. These two minima are shown in Table 1 and we label these two best-fit solutions as the high- and low-mass solutions. The lowest value of χ² corresponds to the high-mass solution and the lower-mass solution has a higher value of χ². Although the high-mass solution is a better fit, we exclude this solution on physical grounds. First of all, it requires a neutron star radius of 28 km, which is not allowed by any known EOS. Secondly, once the neutron star mass and the inclination angle are known, the companion's mass can be calculated (shown in the column labeled M_c in Table 1). In the high-mass case the companion's mass is 1.7 M_☉. Due to the dim nature of the companion, Krauss et al. (2005) have shown that the companion (if a main-sequence star) would have to have a mass that is no bigger than 0.5 M_☉. Clearly this excludes the high-mass solution but allows the lower-mass solution. For these reasons, we exclude the high-mass solutions.

In Table 2, the best-fit solution for each value of 2M/R is shown. In each case only the lower-mass solution is shown. The best-fit solution shown as a solid curve in Figure 1 corresponds to the 2M/R = 0.4 solution shown in Table 2. Although we call this these solutions "lower mass," clearly they still correspond to high-mass neutron stars. In Table 2, only solutions for 2M/R = 0.3–0.6 are shown. In the case of less compact stars, such as 2M/R = 0.2, the only solution corresponded to the "high-mass" branch of unphysical solutions. We did not test solutions that were more compact than 2M/R = 0.6 since these solutions would allow spots to have multiple images, and our program is unable to handle multiple images. Technically, the solutions with 2M/R = 0.6 could have spots with multiple images, but we have checked that the relative values of θ and i do not lead to this problem for our solutions for the most compact stars.

In Figure 4, we show a number of mass–radius curves for stars spinning at 314 Hz as well as the 2 and 3σ confidence regions for the "low-mass" solutions (in this figure, the radius R refers to the equatorial radius). These regions are found by fixing the value of 2M/R, and varying all other parameters and finding the solutions that have χ² larger than the global minimum of χ²_min = 61.0 by δχ² = 4 (2σ) or 9 (3σ). The region allowed with 99.7% confidence (3σ) only includes large stars with high mass. Of the EOS displayed in Figure 4, the only one lying in the 3σ allowed region is L, corresponding to pure neutron matter computed in a mean-field approximation (Pandharipande et al. 1976). A pure neutron core is unlikely to be the correct description of supra-nuclear density matter. However, it is possible for an EOS that includes some softening due to the presence of other species to be allowed by these data.

Zoom In Zoom Out Reset image size

It is possible for variability of the data to affect the fit results, as appears to be the case for SAX J1808 (Leahy et al. 2008). For this reason we have rebinned the data into one-day segments in order to see if there is any significant change in the pulse shape on a day-to-day basis. We performed a χ² test to see how closely each day's data matched the other days' data. Comparisons of one day with an adjacent day gave values of χ²/dof ranging from 0.8 to 1.6, indicating that the day-to-day changes are small. The largest change is between day 810 (June 20) and day 817 (June 27) with χ²/dof = 4.8. Light curves in the two narrow energy bands for these two days are shown in Figure 2.

The data corresponding to these two days are binned into 32 bins per period. Since we continue to remove the six bins corresponding to the second "spot," this corresponds to a total of 104 data points. We fit the data for these two days by assuming that the parameters M, R, i, a are the same for both days. The spot's latitude is allowed to vary, as are the amplitudes of the energy bands and the value of phase. Since the DC contributions found in our previous fits are very small, we do not include terms for DC offsets. In addition, we keep the value of b fixed at the Krauss et al. (2005) value in order to simplify the fits. This corresponds to 12 parameters, and a total of 92 dof.

The best-fit solutions for this two-day joint fit are shown in Table 3. For each fixed value of 2M/R we only found one minimum, unlike the case with all days included. The angular locations of the spot on the two days are labeled θ₁ and θ₂. The change in angular location of the spot between the two days is less than 2° in all cases. The solutions continue to have large masses and radii, as in the case of fits using all of the data. In the case of 2M/R = 0.2 a low-mass (1.2 M_☉) solution is allowed, but it has a very large radius. The 2σ and 3σ confidence regions for the two-day joint fits are shown in Figure 5. The 3σ confidence region is somewhat larger than the same region computed using all of the data, but the two methods have a significant overlap. The two-day joint fit also only allows the stiffest EOS L.

Zoom In Zoom Out Reset image size

The models in this paper depend on the results of the spectral models of Krauss et al. (2005). We now consider the effect of allowing the parameters in the spectral model to vary within the error bars. As mentioned in Section 2.2 we already allow the ratio of the blackbody to power-law components (b) to vary within the 1σ limits given by the Krauss et al. (2005) spectral model. Another spectral parameter that could affect the fits is the photon spectral index Γ for the power-law component of the spectrum. Krauss et al. (2005) found a value of Γ = 1.41 ± 0.06 and in all of our models presented in the previous section we kept the photon spectral index fixed at a value of Γ = 1.40. We would expect that a larger value of Γ would allow for smaller stars. This is because the flux in Equation (1) is proportional to the Doppler boost factor η raised to the power Γ +  3. The Doppler boost factor is mainly responsible for introducing higher harmonics into the signal, so a larger value of Γ creates a more asymmetric pulse shape. In order to compensate, the best-fit solution will require a smaller value of stellar radius R in order to decrease the value of η. In order to test the dependence of the best-fit values of the parameters on Γ, we chose a value of Γ = 1.50 which is somewhat larger than the range allowed by Krauss et al. (2005) and fitted the data using the same method described earlier in this paper. The results of the best-fit parameters for the two values of Γ for the case of 2M/R = 0.4 are shown in Table 4. As expected, increasing Γ allowed for a smaller star, but the decrease is only by 3%. Similar results occur for other values of M/R. Clearly the dependence on the photon spectral index is not sensitive enough to affect the resulting large size of the best-fit stars.

In our models, we keep the spot size (as measured on the star's surface) fixed at a diameter of 3 km. In previous modeling (Leahy et al. 2008) of SAX J1808.4−3658, we found that the final values of the best-fit parameters were not sensitive to changes in the spot size, assuming that the spot is small compared to the size of the star. For this reason we have kept the spot size fixed at 3 km for all models in this analysis. We have also assumed that the hot spot is circular, although recent MHD models (Kulkarni & Romanova 2005) have more complicated spot shapes. In this paper, we have only attempted to use the simplest possible model that is still consistent with the data. Adding extra parameters to our models in order to describe more complicated spot shapes is not yet warranted by the quality of the data.

For all models computed so far, we have made use of an empirical formula for the oblate shape of the star, and have included the relative time delays for photons emitted from different parts of the star. In Table 4 the effects of oblateness and time delays on the fits are shown. For the model labeled "sphere," a spherical initial surface was assumed, but relative time delays were included in the computation. The resulting best-fit solution is about 10% larger than the corresponding oblate model (labeled "oblate" and Γ = 1.4 in Table 4). This shrinkage of the star's radius when oblateness is included has been observed in the modeling of SAX J1808.4−3658 (Leahy et al. 2008). For the model labeled "no td" time delays were omitted from the calculation and a spherical surface was used. The comparison of the two spherical models in Table 4 shows that the model that includes time delays is about 3% smaller than the model that omits time delays.

In our analysis of XTE J1814 we have only included data from the accretion-powered pulsations and have omitted any data corresponding to an X-ray burst. Bhattacharyya et al. (2005) analyzed the light curves constructed from the X-ray bursts for this neutron star. In their analysis they assumed a spherical surface for the star and traced the paths of the X-rays using the Kerr metric. They also made use of a limb-darkened blackbody emission (2 keV) spectral model appropriate for X-ray bursts. Due to the method that they adopted, it was necessary for them to assume one of two different EOS. The stiffer EOS used by Bhattacharyya et al. (2005) is the same as the EOS that we label "APR" and is the A18+ δv +  UIX computed by Akmal et al. (1998). The analysis of the X-ray bursts by Bhattacharyya et al. (2005) allows the APR EOS, while our analysis of the accretion-powered pulsations only allows stiffer EOS. From their analysis it is difficult to determine whether or not their analysis of the X-ray burst data is consistent with a very stiff EOS, as indicated by our analysis.

Watts et al. (2005) provide a detailed analysis of many aspects of both the X-ray bursts and the nonburst emission. One of the quantities that they measured was the fractional amplitude of the pulsations at the fundamental frequency and the first harmonic for both the burst and nonburst emission. They found that the nonburst emission (modeled in this paper) has a larger harmonic content than the burst emission studied by Bhattacharyya et al. (2005). Since Doppler boosting is partially responsible for the harmonic content, it is perhaps not surprising that our analysis of the nonburst pulsations implies a larger Doppler factor for the star, which in turn implies a larger radius than a similar analysis for the X-ray burst oscillations.