CONSTRAINTS ON THE MASS AND RADIUS OF THE NEUTRON STAR XTE J1807-294

One of the primary goals of neutron star astrophysics is the observational determination of the mass–radius relation of neutron stars, thus allowing the determination of the equation of state (EOS) of cold supernuclear density material (Bhattacharyya 2010). Modeling the pulsed light emitted from the surface of a neutron star has the potential to attain this goal (e.g., Leahy 2004a for mass–radius constraints for Her X-1). The pulse shape encodes the emission region properties and the paths of photons traveling through the gravitational field of the spinning neutron star to the observer. For a description of the Schwarzschild metric and of geodesics of light rays, see, e.g., Misner et al. (1973). For an introduction to emitting hot spots on neutron stars, see, e.g., Pechenick et al. (1983).

Excellent candidates for a pulse-shape analysis are the accreting millisecond period pulsars, where the energy released from accretion leads to X-ray emission from the surfaces of the neutron stars. SAX J1808.4-3658 (hereafter SAX J1808) was the first discovered pulsar of this class (Wijnands & van der Klis 1998). The pulsations are thought to be produced from the accretion of plasma funneled onto a hot spot at the neutron star's magnetic poles (see, e.g., Figure 12 of Gierliński et al. 2002). Magnetohydrodynamic simulations for millisecond pulsars with weak magnetic fields (e.g., see Kulkarni & Romanova 2005) show that the material flows from the inner edge of the accretion disk in some complicated way down to the surface to create a hot spot at some position angle offset from the rotation axis. A further complication is that pulse shapes of slowly rotating neutron stars give evidence that their magnetic fields are non-dipolar (e.g., Leahy & Li 1995; Leahy 2004b). Thus, we do not assume any fixed relation between the spot latitude and inner accretion disk radius. Spectral models (Gilfanov et al. 1998; Gierliński et al. 2002) provide strong evidence that the X-rays correspond to blackbody emission from a spot on the star which is then Compton scattered by electrons in a thin layer adjacent to the hot spot. Millisecond pulsars are low-luminosity systems (XTE J1807-294 in outburst was 10³⁶–10³⁷ erg s⁻¹, Falanga et al. 2005). SAX J1808-3658 rises from <10³² erg s⁻¹ in quiescence to 2–4 × 10³⁶ erg s⁻¹ in outburst (Hartman et al. 2008), so the accretion column is expected to be of small height above the surface. Other studies of pulse shapes of millisecond pulsars did not need to add emission above the surface to adequately model the observed pulse shapes (e.g., Poutanen & Gierliński 2003; Leahy et al. 2008). We note also that emission from the sides of an accretion column produces a distinctive fan beam pattern in the observed pulse shape (e.g., Leahy 2004b), which is not seen in the pulse shapes from millisecond X-ray pulsars. We make the assumption that the accretion column is of negligible height for simplicity.

A pulse-shape analysis requires tracing rays from the star to the observer. This is complicated by the rapid rotation of the star, since the metric and geodesics of a rotating neutron star must be computed numerically. However, the Schwarzschild metric, describing the exterior of a non-rotating star is simple and geodesics are simply found by quadrature. This leads to an often-used approximation, the Schwarzschild plus Doppler (S+D; Miller & Lamb 1998; Poutanen & Gierliński 2003). In the S+D approximation, the gravitational field of the star is approximated as purely Schwarzschild, and all relativistic Doppler effects are added to the final expressions for flux as though the star's gravitational field were negligible. The validity of the S+D approximation was investigated (Cadeau et al. 2005, 2007) by computing the geodesics of rapidly rotating neutron stars, using the exact numerical metric to describe the gravitational field and the shape of the star. Cadeau et al. (2005, 2007) showed that it is necessary to include travel times of photons from different parts of the star, but that it was sufficient to use the Schwarzschild geometry to calculate geodesics as the effect of frame dragging is small. They also showed that the neutron star's oblate shape has a significant effect on the light curves due to the "shadowing" effects caused by the non-spherical shape of the star's surface. The shape of the star needs to be included when modeling the light curves of rapidly rotating neutron stars. In Morsink et al. (2007), the oblate Schwarzschild (OS) approximation was introduced as a simple but sufficiently accurate calculation for modeling pulse shapes of accreting millisecond pulsars.

The hot spot is assumed to be infinitesimal in size. We have performed many fits using an extended hot spot of various sizes and shapes (Morsink & Leahy 2011). Within reasonable bounds (less than 3 km spot diameter for a 10 km neutron star radius) the spot shape and size has a negligible effect on the derived mass and radius of the neutron star. The same conclusion was reached by Lamb et al. (2009a, 2009b) who tested variation of spot shapes for spots with angular sizes up to 45°.

The pulsations from two accreting millisecond pulsars have been analyzed in order to constrain the neutron star EOS. Poutanen & Gierliński (2003) modeled the pulse shape of SAX J1808 using the S+D approximation and gave constraints on the mass-to-radius ratio. Leahy et al. (2008) applied the OS model, which includes oblateness and time delays, to give improved constraints on SAX J1808. Both studies only included data from the 1998 outburst which skewed the analyses to very small stars. Subsequently, this neutron star has gone through further outbursts and the resulting pulse profiles were observed to vary over time (Hartman et al. 2008). This variation over time was included in the study of Morsink & Leahy (2011), which separated data with different pulse shapes and modeled these with common mass, radius, and inclination angle parameters, but time-variable hot-spot parameters. The resulting mass–radius constraints (their Figure 7) gave a moderate-mass compact star consistent with a number of mass–radius relations predicted by some common soft EOSs. A similar multi-epoch analysis was carried out for the millisecond pulsar XTE J1814-338 (hereafter XTE J1814; Leahy et al. 2009) which gave a high-mass large-radius star, consistent with a stiff EOS. Studies of more accreting millisecond pulsars are needed to further explore possible EOSs, and to determine the overlap in the mass–radius plane for different cases. Here we present a multi-epoch pulse-shape study of XTE J1807-294 (hereafter XTE J1807).

The accreting millisecond-period pulsar XTE J1807 was discovered (Markwardt et al. 2003) in 2003 February when it went into outburst. The Rossi X-Ray Timing Explorer Proportional Counter Array (RXTE/PCA) observed the pulsar during most of its outburst, which lasted approximately 120 days. Zhang et al. (2006) noticed four events that they call "puny" flares, when the flux increases significantly, although not as much as in a type I X-ray burst. Chou et al. (2008) identified two more "puny" flares and showed that the pulsar phase jumps discretely during the flare events (see Figure 3 of Chou et al. 2008). Pulse profiles can be constructed using data from periods between the flare events.

This paper is a summary of the analysis of the RXTE/PCA pulse shapes of XTE J1807. The data used in the models are described in Section 2. Our model of the pulse profiles includes emission from a hot spot on the surface of the neutron star plus scattered emission from the surface of the accretion disk. The model is described in detail in Section 3. We use this model to put limits on the parameters of the pulse-shape model, including constraints on mass and radius of the pulsar in Section 4.

The observed pulse shapes of XTE J1807 are from RXTE/PCA observations of the 2003 outburst. Since the pulse phase evolution is not continuous during the flare periods, we consider data from times when the pulsar is not undergoing a flare. We group the data into three non-flare epochs: Epoch 1 corresponds to 12697 ⩽ TJD < 12703, Epoch 2 is 12716 < TJD < 12720, and Epoch 3 is 12730 < TJD < 12759. These three non-flare epochs are illustrated in Figure 2(c) of Patruno et al. (2010) using circles; their figure shows a fourth non-flare epoch with very low flux that we do not use since it is too noisy.

All event times of non-flare epochs were first corrected to the barycenter of the solar system and then folded with the ephemeris, consisting of orbital Doppler effect and a fourth-order polynomial pulse phase drift, proposed by Chou et al. (2008). A time averaged pulse shape for each epoch was created by binning the event phases into 32 bins. All the pulse shapes were made separately for eight narrow energy bands (bands 1–8, given by 1: 2–3.7 keV, 2: 3.7–4.9 keV, 3: 4.9–5.7 keV, 4: 5.7–7.3 keV, 5: 7.3–9.0 keV, 6: 9.0–11.5 keV, 7: 11.5–14.8 keV, and 8: 14.8–19.4 keV).

The χ² test was applied to assess the difference in pulse shapes between different epochs and different energy bands. To avoid differences in normalization we compared pulse shapes normalized to unit amplitude. Comparing the pulse shapes of a given energy band between epochs showed that, in all cases, the pulse shapes were significantly different. As a result, we analyze all three epochs of data.

It is not computationally feasible to use all eight energy bands in our modeling. It has been shown through spectral modeling (Poutanen & Gierliński 2003) of SAX J1808 that it is sufficient to use only two energy bands when modeling an accreting millisecond pulsar's pulse profile. In particular, a low-energy band is required to capture the low-energy blackbody emission of the spot which typically peaks near 1 keV. The blackbody emission falls off quickly at higher energies, leaving the higher energy bands almost completely dominated by the emission that has been Compton up-scattered by the electron plasma lying overhead of the hot spot. For this reason we need to choose two energy bands with pulse shapes that are as different as possible in our analysis.

In order to identify the energy bands to be used in our analysis, we compared pulse shapes for the same epoch between energy bands using the χ² test. For Epoch 1, the 2.0–3.7 keV pulse shape is significantly different from all other energy bands (reduced chi-squared: χ²_r > 1.9). Most pairwise comparisons of pulse shapes showed significant differences. The exceptions for the 28 pairs for Epoch 1 were for comparisons between band n and n  +  1 for 1 < n < 7, which had χ²_r < 1.32 (the value that corresponds to less than a 2σ significance). This demonstrates the gradual change of pulse shape with energy and that band 1 is more different than the other bands.

Epoch 2 had similar results: the same pattern of low χ²_r was obtained for comparisons of neighboring energy bands, except with additions of band 3 versus band 5 and band 4 versus band 6 included in the low χ²_r set of pairs. Epoch 3 also gave a similar pattern to Epoch 1, except that band 2 versus band 3 (χ²_r = 2.0) and band 3 versus band 4 (χ²_r = 1.35) pairs were not in the low χ²_r set. The differences in the χ²_r between the three epochs are in large part due to the differing amount of data averaged, and thus the size of the error bars (largest for Epoch 2 and smallest for Epoch 3).

In summary, the 2.0–3.7 keV band is the most different from the other energy bands, and the other energy bands only exhibit small changes in shape. This is consistent with the spectral analysis of data from Epoch 3 by Falanga et al. (2005), which finds a Comptonized component in all energy bands, and an additional blackbody component in the low-energy band. Our pulse-shape model was developed for this case of two spectral components in the low-energy band and one spectral component in a high-energy band. In order to include a significant amount of blackbody emission in our low-energy band, we chose the 2.0–3.7 keV band. For the high-energy band we chose to combine bands 5 and 6, i.e., 7.3–11.5 keV. This gave pulse shapes with smaller error bars and also avoided any possible contamination by iron emission lines or the blackbody component. Figures 1–3 show the observed pulse shapes in these two energy bands for Epochs 1–3, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pulse-shape models used in this analysis are similar to those used by Leahy et al. (2008) and Morsink & Leahy (2011) to model SAX J1808 and by Leahy et al. (2009) to model XTE J1814. These models are similar to those of Poutanen & Gierliński (2003) for SAX J1808, but include time delays and oblateness, which are important, e.g., see Morsink et al. (2007). Our model computes pulse shapes in two energy bands: a high-energy band dominated by Comptonized flux from the hot spot and a low-energy band which includes both Comptonized and blackbody emission from the hot spot. As explained above, we make the assumption that the Comptonized emission from the accretion column is of negligible height above the surface for simplicity.

The hot-spot model is described in detail in our previous work (e.g., Morsink & Leahy 2011; Leahy et al. 2008, 2009). The previous work shows that the pulse shape depends very weakly on the hot-spot size or shape (Morsink & Leahy 2011), so for the current work we use an infinitesimal spot. It is assumed that the Comptonized and blackbody emissions both originate from the same region on the surface of the star. The emissivity of the blackbody component is isotropic and the Comptonized component emissivity is described by an analytic approximation, I∝1 − acos α, where α is the angle between the normal to the surface and the initial direction of the photon and a is a free parameter. The low-energy band includes both blackbody and Compton components, and the phase averaged ratio of the low-energy band blackbody to Compton flux is represented by the parameter b. The blackbody and Comptonized components have effective power-law indices Γ over each observed energy band, which are determined by the observed X-ray energy spectrum. We made use of the spectral model of Falanga et al. (2005) to determine the values of Γ in these two energy bands.

To calculate a model pulse shape, the flux in each phase-bin is computed using the oblate Schwarzschild approximation (Morsink et al. 2007). This accounts for the light-bending caused by gravity, the time delays caused by the finite travel time for photons across the star, and the oblate shape of the rotating star. Since the light-bending depends only on the ratio of the neutron star's mass, M, to radius, R, at the spot location, the star's gravitational radius, 2GM/Rc², is kept fixed for a series of fits. (For simplicity, when referring to the gravitational radius, we use units with G = c = 1.) The lowest χ² for a given 2M/R is found. There are eight free parameters for a fit to a single pair of two-energy-band pulse shapes: R, θ, i, I₁, I₂, b, a, and ϕ. The angle between the spot's center and the spin axis is θ, and the observer's inclination angle to the rotation axis is i. I₁ and I₂ are the normalizations for the high-energy and low-energy pulse shapes and the angle ϕ is an arbitrary rotation phase offset. It is very unlikely that the spot and observer are in opposite hemispheres, so we restrict both θ and i to be less than 90°. The anisotropy is restricted to 0 < a < 1.

In this paper we restrict the models to only one visible spot. Prior to doing this we carried out fits to the data for the three epochs 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 presence of an accretion disk which hides the antipodal spot.

In the current work, our goal is to simultaneously fit the two-energy-band pulse shapes from all three different epochs. In the simultaneous fits, the parameters M, R, and i are taken to have the same values. All other parameters are allowed to change their values. This brings the number of free parameters to 20. Of these, six are normalizations and three are phase offsets which are not of physical interest.

First we carried out test fits for a single epoch (Epoch 1), which had a small number of free parameters (8). We found that the best fits had unacceptably large χ² (205 for 56 degrees of freedom) and had radius values that coincided with our imposed lower limit (4 km). This is similar behavior to what we found for SAX J1808 (Morsink & Leahy 2011), thus we added a scattered light component (described below) and obtained much better fits (χ² = 72 for 54 degrees-of-freedom). Single-epoch fits of Epoch 2 and Epoch 3 were acceptable without the need to add the scattered light component. Joint fits of Epoch 1 with Epoch 2 or with Epoch 3 also had the same problem of high χ² which was cured when a scattered light component was added for Epoch 1. A joint fit to Epoch 2 and Epoch 3 gave acceptable χ² without the scattered light component.

The scattered light model is described in Morsink & Leahy (2011). It accounts for light that travels from the hot spot to an optically thin plasma near the disk and then scatters into the line-of-sight to the observer. Light-travel time delays and the Compton scattering angular dependence are included. The formulae for the scattered light component are given in Morsink & Leahy (2011), specifically we use Equation (9) with an amplitude of I_sc and distance of scattering region from the center of the star of r_s. This equation is the form for scattering from a disk in the equatorial plane of the neutron star. Since the hot spot is at high latitude on the neutron star, the light from the hot spot that scatters to the observer is mainly from parts of the disk that are far enough out that they are not moving relativistically. As shown in Morsink & Leahy (2011), because the amplitude of the scattered light component is small (of the order of a few %), the model used for scattered light is adequate.

With the scattered light component added to the Epoch 1 pulse-shape model, the fits produce acceptable χ² and reasonable parameter values. Since the simultaneous three-epoch fits produce tighter constraints on parameters than single-epoch fits (or simultaneous two-epoch fits), we present those results here.

The best-fit model that simultaneously fits the data from the three epochs is shown as solid and dashed curves in Figures 1–3. The best-fit model corresponds to a neutron star with mass M = 1.47 M_☉ and equatorial radius R = 14.7 km (gravitational radius 2M/R = 0.3) and χ²/dof = 197/172. Other parameter values for the best-fit model are shown in Table 1 in the row for 2M/R = 0.3. A series of fits keeping the gravitational radius (at the spot's location) fixed was performed, and as can be seen from Table 1, the fits for other values of 2M/R < 0.6 are effectively as good as the fit for 2M/R = 0.3 allowing for a wide range of allowed mass and radius values.

Zoom In Zoom Out Reset image size

The parameter θ is the angle between the spin axis and the spot. In our initial fits of these data, we allowed the angle θ to vary between Epochs 1–3. However, the best-fit solutions were always consistent with the value being the same during Epochs 2 and 3. For this reason, we simplified the analysis and introduced a single parameter θ₂₃ to represent the angle between the spin axis and the spot during Epochs 2 and 3. The best-fit values for θ in all epochs are small, consistent with the predictions made by Lamb et al. (2009b). The values of θ in Epoch 1 are smaller than those of Epochs 2 and 3 in Table 1.

Based on the lack of eclipses and the assumption that the companion fills its Roche Lobe, Falanga et al. (2005) show that the binary's inclination angle, i, is less than 83°. This restriction on the inclination angle was used in our fits. The resulting best-fit inclination angle is 38° for 2M/R = 0.3. However, as can be seen by the values shown in Table 1, a wide range of inclination angles is still allowed by the fits. The full range of allowed inclination angles is illustrated in Figure 4, where the value of the inclination is shown as a number next to selected points on the 2σ and 3σ curves. The smallest inclination consistent with the data is 28°, and all values up to a limit of 83° are allowed.

Zoom In Zoom Out Reset image size

The parameter a quantifies the anisotropy of the Comptonized emission, where a = 0 corresponds to completely isotropic emission and a = 1 corresponds to "fan" beaming where emission is most intense for light emitted tangent to the surface. As in the case of the angle θ, the anisotropy parameter was consistent with a constant value in Epochs 2 and 3, so we combined the anisotropy parameters for the last two epochs into one parameter a₂₃. In Table 1, the best-fit values for the anisotropy in Epoch 1 (a₁) are small 0.05 ⩽ a₁ ⩽ 0.08 while the anisotropy in the latter epochs is large 0.61 ⩽ a₂₃ ⩽ 0.76. These values can be understood by noting that in Epoch 1 there is very little difference between the low- and high-energy bands. The low-energy band includes both isotropic blackbody radiation and anisotropic Comptonized radiation, while the high-energy band only includes the Comptonized radiation. A difference in anisotropy will add relative phase shifts as well as pulse-shape changes between the two energy bands. Hence, in order to keep the pulse shapes very similar (and to avoid a phase lag) the best-fit solution minimizes the anisotropy of the Compton scattered component, resulting in a small value for a₁. On the other hand, the high- and low-energy pulse shapes are significantly different in both the latter epochs, allowing for large values of anisotropy.

The parameter b represents the ratio of blackbody radiation to Comptonized radiation in the soft energy band. The pulsar was observed by the XMM-Newton (Campana et al. 2003; Kirsch et al. 2004) and INTEGRAL (Falanga et al. 2005) satellites simultaneously for a period of time (2003 March 22) during Epoch 3. Falanga et al. (2005) fit the combined XMM-Newton, INTEGRAL, and RXTE data (labeled rev-52 in their paper) using a thermal Comptonization model along with blackbody emission from the disk. This spectral model constrains the ratio of the spot's blackbody emission to Comptonized emission to be near 30%, which is included in our models by adding a χ² penalty for values of b₃ that lie outside of the Falanga et al. (2005) error bars. In the first two epochs, the parameter b is allowed to freely vary. In all of the fits, the ratio of blackbody to Compton components in the low-energy band (5% for Epoch 1, 10% for Epoch 2, and 30% for Epoch 3) increased as the source luminosity decreased.

The disk scattering model introduced by Morsink & Leahy (2011) is a simple two-parameter model for the scattering of light emitted by the star by optically thin material near the disk. The two parameters are I_sc, the amplitude of the scattered flux (where the flux from the spot has an amplitude of 1), and r_s, the distance (in km) from the neutron star where the light is scattered. In these models, adding the scattered light component to the first epoch improved the fits, however, the addition of the scattered light to the latter two epochs did not improve the fits. For this reason, the scattered light was only included in Epoch 1. In the best-fit models in Table 1, the scattered light (in Epoch 1) corresponds to 2%–7% of the direct flux, and the scattering location is near 180 km.

Each best-fit solution provides an inclination angle and the neutron star's mass. Making use of the known mass function (Chou et al. 2008) the companion's mass can be calculated. Further, using the assumption that the companion is filling its Roche lobe, its radius can be determined. For the best-fit models shown in Table 1, the companion's mass is near 0.01 M_☉ and its radius is near 0.04 R_☉. According to the models of finite temperature low-mass white dwarf stars computed by Deloye & Bildsten (2003), these values of mass and radius for the companion are consistent with a hot C or O white dwarf, but not with a He white dwarf. The non-detection of the companion during the quiescent period as reported by D'Avanzo et al. (2009) is consistent with a very low mass white dwarf companion.

To explore limits on mass and radius for XTE J1807-294, we carried out an extensive grid of simultaneous three-epoch fits in the mass–radius plane. The resulting 2σ and 3σ limits on mass and radius are shown in Figure 4 with dashed and dot-dashed curves. These curves are the result of a full variation of all the free parameters in the model. In our calculations we only performed fits for gravitational radii in the range 0.2 ⩽ 2M/R ⩽ 0.6 which corresponds to the range of physical values expected for neutron stars. These two boundaries are shown as dotted lines in Figure 4. At selected points (filled circles) along the 2σ and 3σ boundaries, the value of the binary's inclination angle (in degrees) is shown.

We find that high-mass or large-radius stars are favored. As for example, at the 3σ level, the (M, R) values (in units of M_☉ and km) range from (1.67, 8.2) to (0.96, 14.3) along the lower-mass boundary and (2.55, 12.6) to (2.91, 17.3) to (1.56, 24.2) along the higher-mass boundary. The larger radii preferred by these models arise due to the pronounced harmonics visible by eye in Figures 1–3. In this type of pulse-shape model, harmonics are mainly due to the Doppler boosting effect, which requires a large velocity. A large velocity can come about from either a large radius, or a spot located near the equator. However, the small pulse-fraction tends to discriminate against larger values of θ as discussed by Lamb et al. (2009b).

For reference, mass–radius curves for stars spinning at a frequency of 191 Hz are shown for a few representative EOSs in Figure 4. All but one of the EOSs chosen for this figure have the property that a mass of 1.93 is allowed when the star spins at 317 Hz, consistent with the mass measurement of PSR J1614-2230 (Demorest et al. 2010). The exception is EOS Q160 (quark bag model, with a bag constant B^1/4 = 160 MeV Glendenning 2000), which is shown for illustrative purposes. The EOS used in Figure 4 includes hadronic EOSs APR (Akmal et al. 1998), BBB2 (Baldo et al. 1997), and L (Pandharipande et al. 1976); Hyperon EOS H4 (Lackey et al. 2006); mixed phase quark-hadron EOS ABPR1 (Alford et al. 2005); and a quark bag model EOS with a low value for the bag constant, Q140. The data for XTE J1807-294 are only marginally consistent with the Q160 EOS if the star is at the maximum allowed mass for this EOS. The data are consistent with all other EOSs shown in Figure 4. Clearly, any independent observation that constrains one of the free parameters (such as the inclination angle) would provide a stronger constraint on the EOS.

A pulse-shape fitting analysis using the same procedure described in this paper has been performed for two other accreting millisecond-period X-ray pulsars: SAX J1808 (Morsink & Leahy 2011) and XTE J1814 Leahy et al. (2009). Figure 5 shows the 3σ limits on mass and radius for all three pulsars. The fits for SAX J1808 allow for only small masses and radii, while XTE J1814 allows for large masses and radii. The 3σ region found in the present analysis for XTE J1807 encompasses almost all of the region for XTE J1814, and overlaps with a small part of the region for SAX J1808. The mass–radius curves shown in Figure 5 correspond to stars spinning at a frequency of 400 Hz (the spin frequency of SAX J1808). Although the mass–radius curves change when the spin frequency changes, the differences are small enough that for the purpose of this figure they do not alter the results significantly.

Zoom In Zoom Out Reset image size

As long as an EOS mass–radius curve has sections that pass through each of the star's 3σ allowed regions, it will be allowed by all three stars' data. None of the EOS curves in Figure 5 possess this property. However, an EOS that is somewhat stiffer than EOS APR and much softer than EOS L would pass through all three regions. For instance, an EOS mass–radius curve that predicts a radius of 12 km for all masses smaller than about 2.4 M_☉ would satisfy all three constraints. Interestingly, Steiner et al. (2010) found results consistent with this from an analysis of three type I bursters with photospheric radius expansion bursts, thermal emission from three transient low-mass X-ray binaries, and the isolated cooling neutron star RXJ1856-3754. They used a parameterized EOS and a Monte-Carlo/Bayesian analysis to determine an EOS consistent with observations. They found that the EOS is soft near nuclear saturation density, giving radii of 11–12 km for 1.4 M_☉, but the EOS stiffened at higher density to give a maximum mass of 1.9–2.3 M_☉. Their Figure 9 (right) summarized the probability distribution for M, R curves in the mass–radius plane, and their allowed area coincides with the allowed region deduced above from the pulse-shape analyses of SAX J1808, XTE J1814, and XTE J1807.

In summary, the resulting allowed region for mass–radius curves is small but consistent with a mass–radius relation that has nearly constant radius (∼12 km) for masses between 1 and 2.5 solar masses. Observations of additional millisecond pulsars, that have pulse shapes sufficiently bright for pulse-shape modeling, are clearly needed to confirm and refine the results of the current work.

This research was supported by grants from the Natural Sciences and Engineering Research Council of Canada. Y. Chou acknowledges partial support from Taiwan National Science Council grant NSC 100-2119-M-008-025.