THE MASS OF THE BLACK HOLE IN CYGNUS X-1

There is no shortage of published observational data for Cygnus X-1 in the literature. Brocksopp et al. (1999b) provide U, B, and V light curves containing nearly 27 years worth of observations (1971–1997) from the Crimean Laboratory of the Sternberg Astronomical Institute. The binned light curves contain 20 points each and are phased on the following ephemeris:

$$\begin{equation} {\rm Min\ {I}} = {\rm JD }2441163.529(\pm 0.009) + 5.599829(\pm 0.000016)E, \end{equation} \tag{ 1 }$$

where Min I is the time of the inferior conjunction of the O-star and E is the cycle number. In addition, Brocksopp et al. (1999b) provide 421 radial velocity measurements, which are phased to the same ephemeris. The light and velocity curves were kindly sent to us by C. Brocksopp.

In addition to the velocity data of Broscksopp et al., we also made use of the radial velocities published by Gies et al. (2003). We fitted a sine curve to their data and those of Brocksopp et al. (1999a). Fixing the uncertainties on the individual velocity measurements to 7.47 km s⁻¹ and 5.06 km s⁻¹, respectively (these values give χ² ≈ N), we found K = 74.46 ± 0.51 km s⁻¹ for the Brocksopp et al. data and K = 75.57 ± 0.70 km s⁻¹ for the Gies et al. data. The 1σ intervals of the respective K-velocities overlap, the residuals of the fits show no obvious structure, and there is no notable difference between the residuals of the two data sets, apart from the slightly greater scatter in the Brocksopp et al. data. We therefore combined the two sets of radial velocities while removing the respective systemic velocities from the individual sine curves. We phased the Gies et al. data on the above ephemeris before merging these data with the velocity data of Broscksopp et al. The combined data set has 529 points.

All of these velocity and light curve data are discussed further and analyzed in Section 2.4, where they are shown folded on the orbital period and seen to exhibit minimal scatter about the model fits. This small scatter in these data sets, each spanning a few decades, attests to the strongly dominant orbital component of variability. Meanwhile, Cygnus X-1 is well known to be variable in the radio and X-ray bands, including major transitions between hard and soft X-ray states (see Figure 1 in Gou et al. 2011). This raises the question of whether any non-orbital variability in the light curve data could significantly affect the component masses and other parameters determined by our model. (We focus on the light curve data, which are more susceptible to being affected by variability.)

We believe that our results are robust to such non-orbital variability for several reasons, including the following. (1) While Brocksopp et al. (1999a) note that the U, B, and V light curves were correlated with each other, they found no correlation with the radio and X-ray light curves, which is reasonable given that the time-averaged bolometric X-ray luminosity is only ≲ 0.3% of the bolometric luminosity of the O-star (Section 2.3). (2) Although Brocksopp et al. (1999b) did not specifically discuss the photometric variability seen in their 27 year data set, Brocksopp et al. (1999a) do discuss the multiwavelength variability of the source over a 2.5 year period. They noted that there is very little variability in the optical light curves apart from the dominant ellipsoidal/orbital modulations (see Section 2.4) and one weaker modulation with about half the amplitude on a 142 day period (which is thought to be related to the precession of the accretion disk). (3) The light curve and velocity data give remarkably consistent results for the small (but statistically very significant) measured values of eccentricity and argument of periastron (Section 2.4.2). (4) As highlighted above, the small scatter in the folded light curves for a data set spanning 27 years is strong evidence against a significant component of variability on timescales other than the orbital period. In summary, we conclude that non-orbital variability is unlikely to significantly affect our results.

In order to constrain the dynamical model, it is crucial to have a good estimate of the radius of the companion star. However, customary methods of determining this radius fail because the Cygnus X-1 system does not exhibit eclipses nor does the companion star fill its Roche equipotential lobe. We obtain the required estimate of the stellar radius as we have done previously in our study of LMC X-1 (Orosz et al. 2009). The radius, which critically depends on distance, additionally depends on the apparent magnitude of the O-type star and interstellar extinction, and also on the effective stellar temperature and corresponding bolometric correction. The absolute magnitude of the star is M_abs = K  +  BC_K(T_eff, g)  −  (5log D  −  5)  −  0.11A_V, where K is the apparent K-band magnitude, BC_K is the bolometric correction for the K band, D is the distance, and A_V is the extinction in the V band. The luminosity and radius of the star in solar units are $L=10^{-0.4(M_{\rm abs}-4.71)}$ and $R=\sqrt{L(5770/T_{\rm eff})^4}$, respectively. In computing these quantities, we use D = 1.86^+0.12_{− 0.11} kpc (Reid et al. 2011) and a K-band apparent magnitude of K = 6.50 ± 0.02 (Skrutskie et al. 2006), which minimizes the effects of interstellar extinction. For the K-band extinction, we adopt E(B − V) = 1.11  ±  0.03 and R_V = 3.02  ±  0.03 (e.g., A_V = 3.35, Caballero-Nieves et al. 2009) and use the standard extinction law (Cardelli et al. 1989). The bolometric corrections for the K band were computed using the OSTAR2002 grid of models with solar metallicity (Lanz & Hubeny 2003). We note that the K-band bolometric corrections for the solar metallicity models and the models for half-solar metallicity differ only by 0.02 dex (T. Lanz 2007, private communication), so our results are not sensitive to the metallicity.

Figure 1 shows the derived radius and luminosity of the star as a function of its assumed temperature in the range 28,000 K ⩽ T_eff ⩽ 34,000 K. For T_eff = 28,000, the radius is R_dist = 19.26 ± 0.98 R_☉.⁵ As the temperature increases, the radius decreases rapidly at first, and then it plateaus midway through the range, attaining a value of R_dist = 16.34 ± 0.84 R_☉ at T_eff = 34,000 K. Meanwhile, the luminosity increases with temperature, rising from L = 2.1 × 10⁵ L_☉ at T_eff = 28,000 K to L = 3.2 × 10⁵ L_☉ at T_eff = 34,000 K.

Zoom In Zoom Out Reset image size

The effective temperature of the companion star can be determined from a detailed analysis of UV and optical line spectra (Herrero et al. 1995; Karitskaya et al. 2005; Caballero-Nieves et al. 2009). However, it is often difficult to determine a precise temperature for O-type stars owing to a correlation between the effective temperature T_eff and the surface gravity parameter log g. Model atmospheres with slightly smaller values of T_eff and log g give spectra that are very similar to those obtained for slightly higher values of these parameters. Fortunately, log g for the companion star in Cygnus X-1 is tightly constrained for a wide range of assumptions about the temperature, mass ratio, and other parameters because of a peculiarity of Roche-lobe geometry (Eggleton 1983). We show below that our dynamical model constrains the value of the surface gravity to lie in the range log g = 3.30–3.45 (where g is in cm s⁻²).

Based on an analysis of both optical and UV spectra, Caballero-Nieves et al. (2009) obtained for their favored model T_eff = 28, 000 ± 2500 K and a surface gravity log g ≳ 3.00 ± 0.25, which is outside the range of values implied by our model. Based on the plots and tables in Caballero-Nieves et al. (2009), we estimate a best-fitting temperature of T_eff = 30, 000  ±  2500 K when the surface gravity is forced to lie in the range determined by our dynamical model. Using optical spectra, Karitskaya et al. (2005) found T_eff = 30, 400 ± 500 K and log g = 3.31 ± 0.07, which is consistent with our dynamically determined value. Likewise, Herrero et al. (1995) report T_eff = 32, 000 and log g = 3.21 (no uncertainties are given in their Table 1). In the following, we adopt a temperature range of 30,000 K ⩽ T_eff ⩽ 32,000 K.

After considering several previous determinations of the projected rotational velocity of the O-type star and corrections for macroturbulent broadening, Caballero-Nieves et al. (2009) adopt V_rotsin i = 95  ±  6 km s⁻¹. We use this value as a constraint on our dynamical model, which we now discuss.

The eclipsing light curve (ELC) model (Orosz & Hauschildt 2000) has parameters related to the system geometry and parameters related to the radiative properties of the star. For the Cygnus X-1 models, the orbital period is fixed at P = 5.599829 days (Brocksopp et al. 1999b). Once the values of P, the K-velocity of the O-star, and the O-star's mass M_opt are known, the scale size of the binary (e.g., the semimajor axis a) and the mass of the black hole M are uniquely determined. With the scale of the binary set, the radius of the star R_opt determines the Roche-lobe filling fraction ρ. Not all values of R_opt are allowed (for a given P, M_opt, and K): if R_opt exceeds the effective radius of the O-star's Roche lobe, we then set the value of ρ to unity.

The main parameters that control the radiative properties of the O-star are its effective temperature T_eff, its gravity darkening exponent β, and its bolometric albedo A. Following standard practice for a star with a radiative outer envelope, we set β = 0.25 and A = 1.

The ELC model can also include optical light from a flared accretion disk. In the case of Cygnus X-1, the O-star dominates the optical and UV flux (Caballero-Nieves et al. 2009), where the ratio of stellar flux to accretion disk flux at 5000 Å is about 10,000. Consequently, we do not include any optical light from an accretion disk.

We turn to the question of the X-ray heating of the supergiant star and its effect on the binary model. The X-ray heating is computed using the technique outlined in Wilson (1990). The X-ray source geometry is assumed to be a thin disk in the orbital plane with a radius vanishingly small compared to the semimajor axis (this structure should not be confused with the much larger accretion disk that potentially could be a source of optical flux). Points on the stellar surface "see" the X-ray source at inclined angles and the proper foreshortening is accounted for.

Measurement of the broadband X-ray luminosity of Cygnus X-1 (L_xbol; hereafter in units of 10³⁷ erg s⁻¹, adjusted to the revised distance of 1.86 kpc) requires special instrumentation and considerations. There are soft and hard states of Cygnus X-1 (e.g., Gou et al. 2011). The hard state is especially challenging because the effective temperature of the accretion disk is relatively low (T < 0.5 keV), while the hard power-law component (with photon index ∼1.7) must be integrated past the cutoff energy ∼150 keV (Gierlinski et al. 1997; Cadolle Bel et al. 2006). Since the ground-based observations (i.e., photometric data and radial velocity measurements) are distributed over many years, both the range and the long-term average of the X-ray luminosity must be estimated.

The archive of the Rossi X-ray Timing Explorer (RXTE) contains several thousand observations of Cygnus X-1 collected in numerous monitoring campaigns conducted over the life of the mission. We processed and analyzed 2343 exposure intervals (1996 January to 2011 February; mean exposure 2.2 ks) with the Proportional Counter Array (PCA) instrument and we then computed normalized light curves in four energy bands (2–18 keV) in the manner described by Remillard & McClintock (2006). In the hardness–intensity diagram, the soft and hard states of Cyg X-1 can be separated by the value of hard color (HC; i.e., the ratio of the normalized PCA count rates at 8.6–18.0 versus 5.0–8.6 keV), using a simple discrimination line at HC = 0.7. Based on this, we determine that Cygnus X-1 is found to be in the hard state 73% of the time.

Zhang et al. (1997) studied the broadband spectra of Cygnus X-1 with instruments of RXTE and the Compton Gamma Ray Observatory. During the hard and soft states of 1996, they found broadband X-ray luminosities in the range 1.6–2.2 for the hard state and 2.2–3.3 for the soft state. Samples of the hard and intermediate states during 2002–2004 with the International Gamma-Ray Astrophysics Laboratory yielded L_xbol in the range 1.2–2.0 (Cadolle Bel et al. 2006). Additional measurements of the soft state with BeppoSAX (Frontera et al. 2001) found L_xbol in the range 1.7–2.1, while Gou et al. (2011) analyzed bright soft state observations by ASCA and RXTE (1996) or Chandra and RXTE (2010) to determine L_xbol in the range 3.2–4.0.

To get a rough idea of how these special observations relate to typical conditions, we used the available contemporaneous RXTE observations to scale the measured L_xbol values against the PCA count rates, considering hard and soft states separately. We then estimated that the average hard and soft states would correspond to L_xbol values of 1.9 and 3.3, respectively. Finally, the time-averaged luminosity would then be roughly 2.1 × 10³⁷ erg s⁻¹, which corresponds to 0.01 of the Eddington limit for Cygnus X-1. This value is much smaller than the bolometric luminosity of the O-star (L_Bol ≈ 7.94 × 10³⁹ erg s⁻¹), and so we expect that X-ray heating in Cygnus X-1 will not be a significant source of systematic error. We ran some simple tests in which we increased the bolometric X-ray luminosity by up to an order of magnitude and found the light curves to be essentially identical.

When computing the light curve, ELC integrates the various intensities of the visible surface elements on the star. At a given phase and viewing angle, each surface element on the star has a temperature T, a gravity log g, and a viewing angle μ = cos θ, where θ is the angle between the surface normal and the line of sight. ELC has a pre-computed table of specific intensities for a grid of models in the T_eff–log g plane. Consequently, no parameterized limb darkening law is needed. For the specific case of Cygnus X-1, the range of temperature–gravity pairs on the star contained points that are outside the OSTAR2002+BSTAR2006 model grids (Lanz & Hubeny 2003, 2007). We therefore computed a new grid of models, assuming solar metallicity. The range in effective temperature is from 12,000 K to 40,000 K in steps of 1000 K, and the range in the logarithm of the surface gravity log g (cm s⁻²) is from 2.00 to 3.75 in steps of 0.25 dex. The wavelength-dependent radiation fields I_λ(μ), as a function of the cosine of the emergent angle μ, were computed from 232 spherical, line-blanketed, LTE models using the generalized model stellar atmosphere code PHOENIX, version 15.04.00E (Aufdenberg et al. 1998). These models take into consideration between 1,176,932 and 2,296,243 atomic lines in the computation of the line opacity. The radiation fields were computed at 26,639 wavelengths between 0.1 nm and 900,000 nm and at 228 angle points at the outer boundary of each model structure. The models include 100 depth points, an outer pressure boundary of 10⁻⁴ dyn cm⁻², an optical depth range (in the continuum at 500 nm) from 10⁻¹⁰ to 10², a microturbulence of 2.0 km s⁻¹, and they have a radius of 18 R_☉. (At the gravities of interest, the models are very insensitive to the precise value of the stellar radius.) A table of filter-integrated specific intensities for the Johnson U, B, and V filters for use in ELC was prepared in the manner described in Orosz & Hauschildt (2000).

As noted above, the observational data we model include the U, B, and V light curves from Brocksopp et al. (1999b), and the radial velocities from Brocksopp et al. (1999b) and Gies et al. (2003), which are combined into one set. We use the standard χ² statistic to measure the goodness of fit of the models:

$$\begin{eqnarray} &&\chi ^2_{\rm data} = \sum _{i=1}^{20}\left({U_i-U_{\rm mod}(\vec{a})\over \sigma _i}\right)^2 + \sum _{i=1}^{20}\left({B_i-B_{\rm mod}(\vec{a})\over \sigma _i}\right)^2 \nonumber \\ &&\quad+ \sum _{i=1}^{20}\left({V_i-V_{\rm mod}(\vec{a})\over \sigma _i}\right)^2 + \sum _{i=1}^{529}\left({RV_i-RV_{\rm mod}(\vec{a})\over \sigma _i}\right)^2.\quad \end{eqnarray} \tag{ 2 }$$

Here, the notation U_i refers to the observed U magnitude at a given phase i; σ_i is the uncertainty and $U_{\rm mod}(\vec{a})$ is the model magnitude at that same phase, given a set of fitting parameters $\vec{a}$. A similar notation is used for the B and V bands, and for the radial velocities.

Initially we assumed a circular orbit and synchronous rotation, which left us with five free parameters: the orbital inclination angle i, the amplitude K of the O-star's radial velocity curve (the "K-velocity"), the mass of the star M_opt, the radius of the star R_opt, and a small phase shift parameter ϕ to account for slight errors in the ephemerides. We quickly found, upon an examination of the post-fit residuals, that circular/synchronous models are inadequate and more complex models are required. In the end, we computed four classes of models. (1) Model A has a circular orbit and synchronous rotation, and it uses the five free parameters discussed above: $\vec{a}=(i,k,M_{\rm opt},R_{\rm opt},\phi)$. (2) Model B has a circular orbit and nonsynchronous rotation, requiring one additional free parameter Ω, which is the ratio of the rotational frequency of the O-star to the orbital frequency: $\vec{a}=(i,k,M_{\rm opt},R_{\rm opt},\phi,\Omega)$. (3) Model C has an eccentric orbit and synchronous rotation at periastron. For this model there are seven free parameters, namely the five parameters for Model A plus the eccentricity e and the argument of periastron ω: $\vec{a}=(i,k,M_{\rm opt},R_{\rm opt},\phi,e,\omega)$. Finally, (4) Model D incorporates the possibility of nonsynchronous rotation and an eccentric orbit and has a total of eight free parameters: $\vec{a}=(i,k,M_{\rm opt},R_{\rm opt},\phi,e,\omega,\Omega)$.

We also have additional observed constraints that do not apply to a given orbital phase, but rather to the system as a whole. These include the stellar radius (computed from the parallax distance R_dist for a given temperature) and the rotational velocity of V_rotsin i = 95 ± 6 km s⁻¹. For Cygnus X-1, X-ray eclipses are not observed and hence there is no eclipse constraint. As discussed in Orosz et al. (2002), the radius and velocity constraints are imposed by adding additional terms to the value of χ²:

$$\begin{equation} \chi ^2_{\rm con}=\left({R_{\rm opt}-R_{\rm dist}\over \sigma _{R_{\rm dist}}}\right)^2+ \left({V_{\rm rot}\sin i-95\over 6}\right)^2. \end{equation} \tag{ 3 }$$

Finally, for a given model, we have as our ultimate measure of the goodness of fit

$$\begin{equation} \chi ^2_{\rm tot}=\chi ^2_{\rm data}+\chi ^2_{\rm con}. \end{equation} \tag{ 4 }$$

As with any fitting procedure, some care must be taken when considering relative weights of various data sets. In the case of Cygnus X-1, the quality of the light curves is similar in all three bands. We therefore scaled the uncertainties by small factors in order to get χ² ≈ N for each set separately. This resulted in mean uncertainties of 0.00341 mag for U, 0.00155 mag for B, and 0.00226 for V. The uncertainties for the radial velocity measurements are given in Section 2.1.

2.4.1. Principal Results

The fits to the light curve and radial velocity data for Model A (left panels) and Model D (right panels) are shown in Figure 2, and the best-fit values of the parameters are given for all four models in Table 1. These results are for T_eff = 31,000 K; similar figures and tables for other temperatures are presented in the Appendix. For a given temperature, the tabular data reveal a key trend: the total χ² of the fit rapidly decreases as one goes from Model A to Model D, indicating that eccentricity and nonsynchronous rotation are important elements, and we therefore adopt Model D. For final values of the fitted and derived parameters given in Table 2, we take the weighted average of the values derived for each of the temperatures in the range 30,000 K ⩽ T_eff ⩽ 32,000 K. A schematic diagram of the binary based on our best-fitting model is shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Cygnus X-1 Model Parameters

Parameter	Model A	Model B	Model C	Model D	Adopted
i (deg)	51.45 ± 7.46	67.74 ± 2.83	28.50 ± 2.21	27.03 ± 0.41	27.06 ± 0.76
Ω	1.000	0.674 ± 0.043	1.000	1.404 ± 0.099	1.400 ± 0.084
e	...	...	0.025 ± 0.003	0.018 ± 0.002	0.018 ± 0.003
ω (deg)	...	...	303.5 ± 5.1	308.1 ± 5.5	307.6 ± 5.3
Mopt (M☉)	20.53 ± 2.00	26.27 ± 2.41	25.17 ± 2.54	18.97 ± 2.15	19.16 ± 1.90
M (M☉)	7.37 ± 1.19	6.98 ± 0.39	15.83 ± 1.80	14.69 ± 0.70	14.81 ± 0.98
Ropt (R☉)	15.07 ± 1.26	16.41 ± 0.72	18.46 ± 0.77	16.09 ± 0.65	16.17 ± 0.68
Rdist (R☉)	16.42 ± 0.84	16.42 ± 0.84	16.42 ± 0.84	16.42 ± 0.84	16.50 ± 0.84
log g (cm s−2)	3.394 ± 0.016	3.427 ± 0.008	3.306 ± 0.018	3.302 ± 0.012	3.303 ± 0.018
Vrotsin i (km s−1)	106.52 ± 6.39	92.57 ± 5.59	79.91 ± 4.79	92.99 ± 5.89	92.99 ± 4.62
χ2U	26.11	24.74	21.14	17.76	...
χ2B	46.71	43.57	24.96	19.33	...
χ2V	34.81	32.49	24.58	23.99	...
χ2RV	554.13	551.57	531.83	536.32	...
χ2tot	668.03	652.54	614.769	597.67	...
d.o.f.	584	583	582	581	...

Notes. The assumed stellar temperature is T_eff = 31,000 K. R_opt is the stellar radius derived from the models. R_dist is the stellar radius derived from the measured parallax and assumed temperature. V_rotsin i = 95 ± 6 km s⁻¹ is the projected rotational velocity of the O-star derived from the models. Model A assumes a circular orbit and synchronous rotation. Model B assumes a circular orbit and nonsynchronous rotation. Model C assumes an eccentric orbit and pseudosynchronous rotation. Model D assumes an eccentric orbit and nonsynchronous rotation. The adopted parameters are the weighted averages of the values for Model D derived for each temperature in the range of 30,000 K ⩽ T_eff ⩽ 32,000 K.

Download table as:  ASCIITypeset image

Table 2. Cygnus X-1 Final Parameters

Parameter	Value
i (deg)	27.06 ± 0.76
Ω	1.400 ± 0.084
e	0.018 ± 0.003
ω (deg)	307.6 ± 5.3
Mopt (M☉)	19.16 ± 1.90
Ropt (R☉)	16.17 ± 0.68
log g (cgs)	3.303 ± 0.018
M (M☉)	14.81 ± 0.98

Download table as:  ASCIITypeset image

2.4.2. Need for Non-zero Eccentricity

2.4.3. Sensitivity of the Results on the Model Assumptions

We computed models for 14 values of T_eff between 28,000 and 34,000 K. We fit for each temperature separately since the derived radius from the parallax distance depends on the temperature and the derived radius is used as an external constraint. For each model at each temperature, ELC's genetic code was run twice using different initial populations and the Monte Carlo Markov Chain code was run once. The best solutions were then refined using a simple grid search. We computed uncertainties on the fitted parameters and on the derived parameters (e.g., the black hole mass, M, the gravity of the O-star, log g, etc.) using the procedure discussed in Orosz et al. (2002). The results for all 14 temperatures are shown in Table 3 and Figures 4 and 5. As discussed in the main text, the improvement in the χ² values as one goes from Model A to Model D is evident. For Model A, we consistently find R_opt < R_dist. Furthermore, the rotational velocity derived from the model is consistently larger than the observed value. By allowing nonsynchronous rotation for the circular orbit (as in Model B), the model-derived stellar radius agrees with the radius computed from the distance, and the model-derived rotational velocity agrees with the measured one. Generally speaking, the star rotates slower than its synchronous value. However, an inspection of the light curves (see Figure 2 in the main text) shows that the maximum near phase 0.25 is slightly higher than the maximum near 0.75. Because an ellipsoidal model predicts maxima of equal intensity, the fit to the data is not optimal. We accommodate the unequal maxima by adding eccentricity to the synchronous model (as in Model C). However, in this case, the model stellar radius R_opt is consistently larger than the radius derived from the distance R_dist and the model rotational velocity is smaller than the observed value. By allowing nonsynchronous rotation in the eccentric orbit (Model D), R_opt agrees with R_dist and the model rotational velocity agrees with the measured value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size