REFINED NEUTRON STAR MASS DETERMINATIONS FOR SIX ECLIPSING X-RAY PULSAR BINARIES^*

In this section, we review the analytic method introduced by Rappaport & Joss (1983) and Joss & Rappaport (1984) that is widely used to measure the mass of the neutron star and its optical companion (e.g., see, van Kerkwijk et al. 1995b; van der Meer et al. 2007). To begin, one can write the masses of the optical companion and the X-ray source (M_opt and M_X, respectively) in terms of the mass functions:

$$\begin{equation} M_{\mathrm{opt}} = \frac{K_{\mathrm{X}}^3 P(1 - e^2)^{3/2}}{2\pi G \sin ^3 i}(1+q)^2 \end{equation} \tag{ 1 }$$

and

$$\begin{equation} M_{\mathrm{X}} = \frac{K_{\mathrm{opt}}^3 P(1 - e^2)^{3/2}}{2\pi G \sin ^3 i} \left(1+\frac{1}{q}\right)^2, \end{equation} \tag{ 2 }$$

where K_X and K_opt are the semiamplitudes of the respective radial velocity curves, P is the period of the orbit, i is the inclination of the orbital plane to the line of sight, e is the eccentricity of the orbit, and q is the mass ratio defined as

$$\begin{equation} q \equiv \frac{M_{\mathrm{X}}}{M_{\mathrm{opt}}} = \frac{K_{\mathrm{opt}}}{K_{\mathrm{X}}}. \end{equation} \tag{ 3 }$$

The values for K_X and P can be obtained very accurately from X-ray pulse timing measurements (the projected semimajor axis of the pulsar's orbit in light-seconds, a_Xsin i, is usually quoted in publications, from which one finds K_X = 2πc a_Xsin i/P) and optical and/or UV spectra can provide a value for K_opt.

Assuming a spherical companion star, the inclination of the system is related to the eclipse half-angle⁵ θ_e, the stellar radius R, and the orbital separation a by

$$\begin{equation} \sin i = \frac{\sqrt{1-(R/a)^2}}{\cos \theta _e}. \end{equation} \tag{ 4 }$$

Following the approach in Rappaport & Joss (1983), the radius of the companion star is some fraction of the effective Roche lobe radius

$$\begin{equation} R = \beta R_{\rm L}, \end{equation} \tag{ 5 }$$

where R_L is the sphere-equivalent Roche lobe radius. We will refer to the fraction β as the "Roche lobe filling factor." Combining Equations (4) and (5) yields

$$\begin{equation} \sin i = \frac{\sqrt{1 - \beta ^2 \left({R_{\rm L}}/{a} \right)^2}}{\cos \theta _e}. \end{equation} \tag{ 6 }$$

Rappaport & Joss (1983) provide an approximate expression for R_L/a, the ratio of the effective Roche lobe radius and the orbital separation:

$$\begin{equation} \frac{R_{\rm L}}{a} \approx A + B \log q + C \log ^2 q, \end{equation} \tag{ 7 }$$

where the constants A, B, and C are

$$\begin{equation} A = 0.398 - 0.026 \Omega ^2 + 0.004 \Omega ^3 \end{equation} \tag{ 8 }$$

$$\begin{equation} B = -0.264 + 0.052 \Omega ^2 - 0.015 \Omega ^3 \end{equation} \tag{ 9 }$$

$$\begin{equation} C = -0.023 - 0.005 \Omega ^2. \end{equation} \tag{ 10 }$$

Here, Ω is the ratio of the rotational frequency of the optical companion to the orbital frequency of the system. In other words, it is a measure of the degree of synchronous rotation, where Ω = 1 is defined to be synchronous. These four expressions give the value of R_L to an accuracy of about 2% over the ranges of 0 ⩽ Ω ⩽ 2 and 0.02 ⩽ q ⩽ 1 (Joss & Rappaport 1984). If the orbit is eccentric, the value of β is defined at periastron. The star is assumed to have the same volume over its entire orbit (Avni 1976), and at any given phase outside periastron the value of β is adjusted accordingly. With a given eccentricity e, the orbital phase of the X-ray eclipse is determined by the argument of periastron ω.

For a given system, one can calculate the neutron star mass using Equations (1) through (10) when given values of P, a_Xsin i, θ_e, K_opt, Ω, and β (and if the orbit is eccentric, e and ω). It is possible to estimate Ω by measuring the projected rotational velocity, v_rotsin i, of the optical companion star. This process is described in van der Meer et al. (2007) and is employed there for the systems SMC X-1, LMC X-4, and Cen X-3. Finally, one must assume some value for the Roche lobe filling factor, β. Many of the wind-fed systems discussed here are thought to be close to filling their Roche lobes, and a range of β between 0.9 and 1.0 is usually adopted (Rappaport & Joss 1983; van der Meer et al. 2007). Since all the measured input quantities are not known exactly, a simple Monte Carlo technique may be used to derive the most likely values of M_X, M_opt, and i, and the corresponding 1σ confidence limits.

To determine neutron star masses in this manner for the six eclipsing systems where all of the necessary quantities are known or estimated (see Tables 1 and 2 and references therein), we assume that any value within the range 0.9 ⩽ β ⩽ 1.0 is equally likely. We then use a Monte Carlo technique to generate a distribution of resulting neutron star masses, as shown in Figure 1. This provides a clear visual representation of uncertainties through the shapes of each distribution and it provides an estimate of the "most likely" mass for each system. The neutron star masses and system inclinations found in this manner are presented in Table 3. The masses derived for SMC X-1, LMC X-4, and Cen X-3 agree very well with those cited in van der Meer et al. (2007).

Zoom In Zoom Out Reset image size

We note that there is no physical solution for either Vela X-1 or 4U 1538-52 when using this technique with an eccentric orbit, because the quantity in Equation (6) is larger than unity. Adjusting the argument of periastron ω and/or the eclipse duration θ_e can force a solution. However, for 4U 1538-52, no solution exists within the 1σ uncertainties of ω and θ_e. This discrepancy arises due a high inclination and is discussed further in Sections 4.1 and 4.2. As a workaround, when employing this analytic technique we use a larger eclipse width (θ_e = 33° ± 3°; from van Kerkwijk et al. 1995a) for Vela X-1 and adopt a circular orbit (e = 0) as in Clark (2000) for 4U 1538-52.

The analytic method presented in Section 2.1 is straightforward and is easy to implement on a computer. However, this method relies on the following two approximations.

• 1.  
  The computation of the effective Roche lobe radius, R_L/a, from Equations (7)–(10).
• 2.  
  The computation of the X-ray eclipse duration, 2θ_e, from Equation (6).

We use the Eclipsing Light Curve (ELC) code of Orosz & Hauschildt (2000), which is based on Roche geometry, to test these two approximations, and we discuss each one in turn.

Strictly speaking, "Roche geometry" applies only to binary systems with circular orbits and co-rotating stars. Numerous authors have presented generalizations of the Roche potential to account for situations in which one or both of these assumptions are not met (e.g., Avni & Bahcall 1975; Avni 1976; Wilson 1979; Avni & Schiller 1982). These generalizations do not fully describe complete dynamics of the star, and as a result some small approximations are involved (e.g., Wilson 1979). The main assumption is that the timescale for the internal motions of the star that are required for the star to adjust to the varying potential is considerably shorter than the orbital period. Given this, one can compute an effective potential locally at each orbital phase without significant inconsistency (Wilson 1979). Although these generalized potentials are widely used, it is not known how well they work in practice. For most of the systems discussed here, the orbits are circular and the stars rotate close to the synchronous rate, so the modified potential we use (the ELC code uses the potential given in Wilson 1979) should be fairly accurate. Finally, we note that Rappaport & Joss (1983) also adopted a modified "Roche potential" in their analysis—their fitting functions are approximations to numerical integrations of the critical potential surface. Therefore, any systematic error that ELC would have owing to improper generalizations of the Roche model would also be present in the work of Rappaport & Joss (1983) and others that used these fitting functions such as van der Meer et al. (2007).

The shape and size of the critical potential surface (hereafter the "Roche lobe") depend only on the mass ratio q and the parameter Ω. When given q and Ω, it is straightforward to define the equipotential surface (from the value of the gravitational potential at the inner Lagrangian point) and to numerically integrate its volume. The sphere-equivalent radius R_L then follows. We define a large grid of points in the q–Ω plane and compute values of R_L, and compare them with the values of R_L found from Equations (7) through (10). The results are shown in Figure 2. Rappaport & Joss (1983) claim an accuracy of their fitting functions of about 2% over the stated range (0 ⩽ Ω ⩽ 2 and 0.02 ⩽ q ⩽ 1), and our results confirm this.

Zoom In Zoom Out Reset image size

From Equation (6), one can see that the duration of the X-ray eclipse depends on the inclination i, the Roche lobe filling factor β, and R_L/a, which is a function of the mass ratio q and the parameter Ω. The ELC code can be used to compute the duration of the X-ray eclipse for a given geometry. Rather than using ray tracing to determine whether a point is eclipsed by the companion star (e.g., see, Chanan et al. 1976), ELC locates the limb of the star to high accuracy by testing the viewing angles of each surface element. ELC then uses bisection to find, at each latitude row, the longitude of the point that has a viewing angle of μ = 0. Once found, the points on the limb define a polygon in sky coordinates. At that same phase, the location of the X-ray source (assumed to be a point source) in sky coordinates is determined. A simple test is used to determine if the sky coordinate of the X-ray source is inside or outside the polygon defined by the horizon of the star. ELC uses another bisection routine to find the orbital phase of the X-ray eclipse ingress to high accuracy. If the orbit is circular, the eclipse half-angle θ_e is equal to the ingress phase. If the orbit is eccentric, the X-ray eclipse egress phase is also computed, and the eclipse half-angle is computed from both the ingress and egress phases.

Setting Ω = 1, we use ELC to compute the full duration of the X-ray eclipse for a wide range of values in the q–i plane, using β = 1.0 and β = 0.9 and assuming a circular orbit. The full eclipse duration was also computed from Equation (6), and the difference between the numerically computed value of 2θ_e and the analytically computed value of 2θ_e was determined. Figure 3 shows the differences for β = 1 and Figure 4 shows the differences for β = 0.9. The differences can be quite extreme: they are in excess of 10° for small mass ratios and large inclinations and are less than −10° for grazing eclipses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The approximate locations in the q–i plane for three systems (SMC X-1, LMC X-4, and Cen X-3) are also shown in Figures 3 and 4. Since Her X-1 is a low-mass X-ray binary, its mass ratio does not appear within the limits of the two figures. Vela X-1 and 4U 1538-52 are excluded because of their eccentric orbits. The instantaneous Roche lobe filling factor of an eccentric system during X-ray eclipse will be less than the value of β as calculated for non-eccentric systems. When β = 1, all three systems shown in Figures 3 and 4 are near the contour denoting zero difference. When β = 0.9, all three systems are above the zero-difference contour, which indicates that the eclipse durations computed numerically are longer than those computed analytically.

To illustrate the differences between the analytic and ELC results, Figures 5 and 6 show a system resembling Cen X-3 in sky coordinates, where the Roche lobe filling factor is β = 0.9. Figure 5 demonstrates that the companion star is not spherical. The parts of the star near its equator extend beyond the circle denoting the volume-equivalent sphere. Hence, one would expect that the duration of an X-ray eclipse for very high inclinations would be longer than what one would compute from the analytic approximations, and a glance at Figures 3 and 4 confirms this. In a similar manner, Figures 3 and 4 show that for grazing eclipses (i.e., lower inclinations), the numerically computed durations are shorter than the analytic approximations. One can see from Figure 5 that the polar regions of the companion star are inside the circle denoting the volume-equivalent sphere, and as a result there would be inclinations at which the analytic approximations indicate X-ray eclipses when in fact none occur.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows a magnified view of Figure 5 near the neutron star eclipse. In this example, the egress phase of the X-ray eclipse is very close to 33°, since the neutron star is just crossing the limb of the companion star. The limb of the analytically computed volume-equivalent sphere is well inside the limb of the star. The egress phase of X-ray eclipse would be near 32° when the analytic expressions are used, since the neutron star is just crossing the limb of the sphere at that phase. Thus, in this example, the full duration of the X-ray eclipse computed numerically is a full 2° longer than the duration computed analytically.

Having tested both approximations, we conclude that the ELC code does offer a significantly more accurate representation of the physical system than the analytic method presented in Section 2.1. The size of the effective Roche lobe radius as a function of q and Ω is relatively well represented by the analytic formulae and ELC offers only a modest improvement. However, the difference in X-ray eclipse duration can be quite extreme (as much as ±10°) and has a direct effect on the calculation of the neutron star mass.

Vela X-1 has an 8.96 day orbital period and an eccentric orbit with e = 0.0898 ± 0.0012 (Barziv et al. 2001). The argument of periastron for the companion star is ω = 33259 ± 092 (Bildsten et al. 1997), where ω = ω_X + 180. The optical V light curve shown in the upper left panel of Figure 8 for Vela X-1 is binned data from the All Sky Automated Survey (Pojmansky 2002). In our preliminary analysis of Vela X-1 (e.g., Figure 7), we adopt a semiduration of the X-ray eclipse θ_e = 33° ± 3° (van Kerkwijk et al. 1995a). We later use the more precise value from Kreykenbohm et al. (2008), θ_e = 34135 ± 05. Unfortunately, the analytic technique from Section 2 does not arrive at a physical solution with this longer eclipse duration due to the system's high inclination. Specifically, solving for the inclination i is not possible using Equation (6) as sin  i>1 for all the Monte Carlo simulations. The numerical ELC code does not have this same limitation. As a workaround, we keep θ_e = 33° ± 3° for all instances of the analytic case instead. Our final derived mass for Vela X-1 is 1.77 ± 0.08 M_☉ with a system inclination of 778 ± 12.

4U 1538-52 has a 3.73 day orbital period and most likely an eccentric orbit. Clark (2000) and Mukherjee et al. (2006) give e = 0.174 ± 0.015 and 0.18 ± 0.01, respectively, which are in agreement. However, Makishima et al. (1987) give a much lower e = 0.08 ± 0.05 and van Kerkwijk et al. (1995a) adopt e = 0 in their analysis; even Clark (2000) provides an alternate set of fit parameters for a circular orbit. As with Vela X-1, we find no physical analytic solution for 4U 1538-52 with an eccentric orbit and its reported eclipse duration (θ_e = 285 ± 15) due to the system's high inclination. Therefore, we follow the approach of Clark (2000) and present numerical results for both an eccentric orbit and a circular orbit.

In addition, Clark (2000) and Mukherjee et al. (2006) give significantly different values for ω, the argument of periastron: 244° ± 9° and 220° ± 12° for the optical source, respectively. Following the approach of Zhang et al. (2005), we have extrapolated a linear decrease of ω based on the time difference between the observations made by Clark (2000) and Mukherjee et al. (2006). We estimate ω is decreasing by 00010 ± 00006 per day and that our subsequent observations should therefore have ω = 198° ± 14°. This is the value we adopt in the numerical model. Even though using a lower ω nudges the parameter space closer to a single physical analytic solution, the reported eclipse duration is still too short to allow sin i < 1 (as in Equation (6)), and we are forced to retain a model where e = 0 for the analytic case.

Due to the above discrepancies, the reportedly low neutron star mass (∼1 M_☉; van Kerkwijk et al. 1995a), and the lack of a published light curve, we found this system especially worthy of our attention.

Optical light curves in BVI for 4U 1538-52 were obtained at the Cerro Tololo Inter-American Observatory on the 1.3 m SMARTS telescope with the ANDICAM in 2009 June–September. There are a total of 39 images in each filter taken on different nights, each with a 60 s exposure time. Standard pipeline reductions were done on all images, and differential photometry was performed in IRAF for the target and several comparison stars. Light curves are shown in Figure 9, and observations in all three filters were incorporated into our numerical analysis.

Twenty-one high-resolution spectra of 4U 1538-52 were also taken on several nights in 2009 July and August at Las Campanas Observatory on the 6.5 m Clay Magellan telescope with the MIKE spectrograph. Standard pipeline reductions were done on all frames, including heliocentric corrections. The spectroscopic analysis was performed in IRAF using cross-correlation of the blue half of the spectrum (4750–4950 Å) with a model B0 star. Radial velocities were then computed and incorporated into our numerical analysis. The radial velocity curve is plotted in the bottom panels of Figure 9 with the adopted ELC model solution.

We ran two sets of model fits. We assumed an eccentric orbit with eccentricities in a narrow range around e = 0.18 and also a circular orbit. We find K_opt = 14.1 ± 1.1 km s⁻¹ for an eccentric orbit and K_opt = 21.8 ± 3.8 km s⁻¹ for a circular orbit. The difference between these two measurements is due in part to the noisy data and incomplete phase coverage near phase 0.25. For comparison, Reynolds et al. (1992) found K_opt = 19.2 ± 1.2 km s⁻¹ (uncorrected) and K_opt = 19.8 ± 1.1 km s⁻¹ (corrected for tidal distortions) assuming a circular orbit. Our final mass for 4U 1538-52 is quite low: 0.87 ± 0.07 M_☉ for an eccentric orbit and 1.00 ± 0.10 M_☉ for a circular orbit. If the orbit is indeed eccentric, this is an extremely low neutron star mass.

SMC X-1 has a 3.89 day orbital period. It also exhibits a superorbital X-ray cycle that is probably caused by precession of either the accretion disk or neutron star (e.g., Priedhorsky & Holt 1987). The optical V light curve for SMC X-1 shown in the lower left panel of Figure 8 is from van Paradijs & Kuiper (1984). Our final derived mass for SMC X-1 is 1.04 ± 0.09 M_☉ with a system inclination of 685 ± 52. We discuss the implications of this low-mass result and that of 4U 1538-52 in Section 5.

LMC X-4 has a 1.41 day orbital period and, like SMC X-1, also exhibits a superorbital X-ray cycle. Heemskerk & van Paradijs (1989) performed an extensive study of the long-term variations in X-ray flux of LMC X-4 and concluded that it contained a warped, precessing accretion disk. The optical B light curve for LMC X-4 shown in the upper right panel of Figure 8 is from Ilovaisky et al. (1984). Since no data table was available, the points were extracted from their "X-ray OFF states" plot using DEXTER available via the SAO/NASA Astrophysics Data System. Our final derived mass for LMC X-4 is 1.29 ± 0.05 M_☉ with a system inclination of 670 ± 19.

Cen X-3 has a 2.09 day orbital period. The optical V light curve for Cen X-3 shown in the lower right panel of Figure 8 is from van Paradijs et al. (1983). Our final derived mass for Cen X-3 is 1.49 ± 0.08 M_☉ with a system inclination of 667 ± 24.

Her X-1 has a 1.7 day orbital period and a much lower companion star mass than the other five systems of interest (∼2 M_☉). Like SMC X-1 and LMC X-4, Her X-1 also exhibits superorbital X-ray cycles.

No optical light curve was used for Her X-1 because the optical K-velocity has a very large uncertainty (90 ± 20 km s⁻¹) due to uneven X-ray heating of the companion star. As a result, any spectral lines measured do not on average emanate from the star's center of mass. An optical light curve would therefore be ineffective in constraining the fit. Further, the companion's Roche lobe filling factor β is likely very close to 1 (between 0.95 and 1) as adopted by and discussed in van Kerkwijk et al. (1995a) and Bahcall & Chester (1977). We therefore set β = 1 in our analysis. Our final derived mass for Her X-1 is 1.07 ± 0.36 M_☉ with a system inclination greater than 859.