NEW EVIDENCE FOR A BLACK HOLE IN THE COMPACT BINARY CYGNUS X-3

We extracted data from the HEASARC archives covering an approximately nine year span starting at modified Julian day 50,319 (where MJD=JD–2,400,000.5). RXTE was our primary database, but we also examined several epochs of INTEGRAL and BeppoSAX coverage as well, but ultimately chose to use the RXTE data almost exclusively in our final mass-estimation analysis. Some 35 spectra were examined, spanning a large range of hardness and intensity variations. Figure 1 depicts the distribution of these data in hardness–intensity space. As is well known, Cyg X-3 undergoes extreme spectral evolution, the full range of which cannot be represented by a single model; see, e.g., Hjalmarsdotter et al. (2009). As we describe in detail in the next section, we found satisfactory fits to our model for the subset of the data which were not at the extreme range of Figure 1, i.e., the very hard and very soft states of the system.

Zoom In Zoom Out Reset image size

The data reduction and analysis was performed using current HEADAS software release (ver. 6.7) and the RXTE instrumental response file release in 2009 August. For our spectral analysis, we use the Proportional Counter Array (PCA) Standard-2 data mode and the Standard Archive HEXTE Mode. Standard dead-time correction was applied to all spectra. Spectra were modeled using the XSPEC analysis package. We used approximate energy intervals of 3.0–30.0 keV and 20.0–200.0 keV for the PCA and HEXTE data correspondingly, although depending on the spectral and intensity state in question, additional channel selection cuts were made as warranted. A systematic error of 1% was assumed to account for uncertainty of the absolute instrumental calibration. The stated uncertainties on spectral model parameters were calculated using 1σ confidence interval.

We examined each individual observation applying our spectral fitting model which we describe below. We note that in many instances, we were simply unable to achieve a satisfactory fit to the data. The complex variations of the local absorption are likely the underlying issue here, as has been pointed out by others; see, e.g., Zdziarski et al. (2010). This is consistent at least with Figure 1, where the subset of the data we ultimately selected for our mass-determination analysis is depicted. The hardness–intensity space is roughly divided there into extreme-soft, moderate, and extreme-hard spectral regimes. In the latter, the effects of Comptonization from the ambient plasma or from electron populations associated with a collimated outflow or with coronal electrons lead to spectral hardening. For the extreme-soft case, the hard power law is very poorly constrained (or undetected). An example is depicted in Figure 2, which is based on an INTEGRAL observation. We note that although the soft-X-ray detector on INTEGRAL, JEM-X, is much less sensitive than the RXTE PCA, the IBIS hard-X-ray detector is superior in the hard-X-ray regime. Clearly, the hard-power law is not reliably determined. Furthermore, with the bolometric luminosity strongly skewed to the keV regime, uncertainties in the flux and spectral shape due to absorption are substantial.

Zoom In Zoom Out Reset image size

We do, however, obtain satisfactory fits to our model for a subset of about half of the observations studied. In these cases, the reduced χ² statistics are of order unity resulted and meaningful constraints on the key parameters of our model were achieved. These data are represented in the hardness–intensity space by the blue points in Figure 1, which span a range of about 0.1–0.35 in hardness defined by the ratio of PCA counts in the (9–20)keV/(2–9)keV bands. The rationale for this selection is as follows. The extreme-soft-state spectra cannot be used for index–mass accretion rate correlation. As noted, this is because these emergent spectra are the result of reprocessing of the emitted spectra in the optically thick cloud covering the central source. In this case, the entirety of the information on the nature of the central source is lost. In the extreme-hard-state case, we speculate that an additional Comptonization component associated with plasma ambient to the central accretion flow contributes to a hardening and curvature of the continuum above 10 keV. This is counter to the basic converging inflow scenario we propose. The specific data sets included in our mass-determination analysis, along with the inferred spectral parameters, are enumerated in Table 1.

Our spectral model consists of the so-called bulk-motion Comptonization (BMC) model (see Titarchuk et al. 1997, hereafter TMK97), convolved with model corrections for (Galactic plus local) absorption, Fe Kα emission and absorption edge, and an exponential cutoff. The BMC model is a convolution of an input thermal spectrum having normalization N_bmc and color temperature kT_c with Green's function encoding the Comptonization process, which is a broken power law. This approach is applicable to any type of Comptonization process. The generic shape of a Comptonization Green function is a broken power law which is followed by exponential cutoff at energies higher than the kinetic energy of matter and it is independent of the type of the Comptonization process, e.g., bulk or thermal. An example of our model fitted to the data is depicted in Figure 3.

Zoom In Zoom Out Reset image size

As a side note, this Green function technique is well suited to the analysis of a wide variety of problems and has been applied, for example, to a description of the income distribution in financial analysis by Titarchuk et al. (2009). The normalization N_bmc is a ratio of the source luminosity to the square of the distance. The spectral index parameter of the model has been found to correlate approximately linearly with the normalization up to a saturation level, typically about 2.2–2.8, after which it flattens. For various sources the same effect is seen in photon index Γ versus quasi-periodic oscillation frequency (see ST09). Furthermore, the self-similar nature of an ensemble of observations can be related by scaling laws.

The resulting spectra are characterized by the inferred parameters log(A) related to the Comptonized fraction f as f = A/(1 + A) and to a spectral index α = Γ − 1. There are several advantages of using the BMC model with respect to the more common approach, i.e., the application of an additive combination of a blackbody/multi-color-disk plus power-law/thermal Comptonization. First, by its nature, it is applicable to the general case of photon energy gain through not only thermal Comptonization but also via dynamic (bulk) motion Comptonization (see Shaposhnikov & Titarchuk 2006 for details). Second, the BMC spectral shape has an appropriate low energy curvature, which is not the case with a simple additive power law. This is essential for a correct representation of the lower energy spectrum. This can lead, for example, to inaccurate N_H column values and thus produce a non-physical component "conspiracy" involving the high-energy spectral cutoff and spectral-index parameterizations. Specifically, when a multiplicative component comprised of an exponential cutoff is combined with our model, the e-fold energies are in the expected range of greater than 20 keV. When it is applied to an additive model including a power law, the cutoff can often extend below 10 keV, resulting in unreasonably low values for the inferred photon index. As a result, the implementation of such phenomenological models makes it much harder to correctly identify the spectral state of the source, which is a central task of our study. Additionally, a more important property of the model is that it self-consistently calculates the normalization of the "seed" spectral component, which is expected to be a robust mass accretion rate indicator. We further note that the Comptonized fraction is properly evaluated by the BMC model.

We considered the Cyg X-3 spectra approximately resembling the galactic BH low-hard state and intermediate state as the basis of our study; see, e.g., McClintock & Remillard (2006). We note that only a small part of the disk emission component is seen directly. The energy spectrum is dominated by a Comptonization component very often approximated by a power law at energies above ∼5 keV. To calculate the total normalization of the "seed" blackbody disk component, we model the spectrum using the model described in the previous section. As detailed in Section 2, we argue that the disk emission normalization calculated using this approach produces a more robust correlation than that one obtains using additive models.

The presence of index saturation in various BH-ray binaries—Cyg X–1, GRO J1655–40, GRS 1915+105, GX339–4, H 1743–322, 4U 1543–47, and XTE J1650–500—has been demonstrated previously (ST07, ST09, TS09). The self-similar evolutionary tracks of those analyses have led to BH mass estimates that are generally in good agreement with dynamical determinations. The basic idea involves a Compton-scattering cloud. Its scale, temperature, and optical depth regulate the nature of the emergent spectrum. We postulate that such a cloud is a natural consequence of an adjustment from a Keplerian flow to a turbulent, innermost sub-Keplerian boundary near the central compact object. This barrier ultimately leads to the formation of a transition layer between these regions. This interpretation predicts two phases or states dictated by the photon upscattering produced in the transition layer: (1) an optically thin and very hot (kT_e∼50 keV) medium producing photon upscattering via thermal Comptonization; the photon spectrum index Γ ∼ 1.5–1.7 for this state is dictated by gravitational energy release and Compton cooling in an optically thin shock near the adjustment radius; and (2) an optically thick and relatively cool medium (kT_e ∼ 5–10 keV); the index for this state, Γ≳ 2.1 is determined by soft-photon upscattering and photon trapping in converging flow into the BH. Details are presented in ST09 and references therein.

The scaling technique is based on the empirical correlation patterns observed during spectral transitions which effectively form a one-parameter space; the mass of the BH defines a specific correlation curve (ST07, ST09). After scalable state transition episodes are identified for two sources, the correlation pattern Γ–N_bmc for a reference transition is parameterized in terms of the analytical function in Equation (1). By fitting this functional form to the correlation pattern, we find a set of parameters A, B, D, N_tr, and β that represent a best-fit form of the function F for a particular correlation curve:

$$\begin{equation} F(x)=A-D B\ln \lbrace \exp [(1.0-(N/N_{\rm tr})^\beta))/D]+1\rbrace. \end{equation} \tag{ 1 }$$

For N ≫ N_tr, the correlation function F(x) converges to a constant value A. Thus, A is the value of the index saturation level, β is the power-law index of the left-hand portion of the curve, and N_tr is the value at which the index transitions. Parameter D controls how smoothly the fitted function saturates to A. We scale the data to a template by applying a transform N → s_N × N until the best fit is found. We note that the scaling fit to determine s_N has to be performed with all the parameters describing the shape of the pattern fixed in order to comply with the pattern scalability requirement. We then use the Γ − N_bmc data points of GRS 1915+105, which has a reasonably good dynamical mass determination of 14 ±  4 M_☉ (Greiner et al. 2001), which is consistent with the value of 15.6 ± 1.5 M_☉ of ST07. We adopt the latter value in our subsequent analysis. The reference data set is well represented by the function defined in Equation (1). We then fit F(x) varying the parameters A, N_tr, β to best represent the GRS 1915+105 measurements. To test for self-consistency we have used both the least-squares and χ² statistical minimization methods to determine the best-fit parameter values. The least-squares method uses the values of index and normalization with equal weighting, while in the χ² minimization the data are weighted by the inverse error square values. Because data points have different error bars due to different exposure, power law fraction, etc., and because the points can be clustered, the least-squares method is less biased. Moreover, the χ² method results in a very high reduced χ² statistics due to large scatter in the data, which makes it less robust during the error calculation. We therefore select the least-squares method as our primary minimization technique. Parameters D and B are not well constrained by the data. We fix those parameters at 0.3 and 1.7 respectively based on our previous experience with correlation parameterization (ST09). For free parameters we obtain the following best-fit values with the least-squares method: A = 3.02 ± 0.04 (saturation level), N_tr = 0.138 ± 0.004, and β = 2.03 ± 0.14. We then obtain the scale factor, s_N = 0.85 ± 0.07 (see Figure 4). We use this value for our BH mass calculations in Cyg X-3. We note that the scale factor obtained by the χ² minimization approach, s_N = 0.75 ± 0.08, leads to higher BH mass estimates, are consistent within uncertainties. In this sense, the least-squares method we chose represents the more conservative of the two approaches. According to the scaling laws (see Equations (8) and (9) in ST09), the expression for the BH mass in Cyg X-3 is

$$\begin{equation} m_{\rm Cyg x-3}=m_{1915}\frac{f_G}{s_N}\left(\frac{d_{\rm Cyg{} X-3}}{d_{1915}}\right)^2. \end{equation} \tag{ 2 }$$

Using published mass and distance determinations for GRS 1915+105 of d₁₉₁₅ = 12.1 ± 0.8 kpc (Greiner et al. 2001), M₁₉₁₅/M_☉ = 15.6 ± 1.5 (ST07) along with the value of s_N determined in our least-squares analysis leads to a BH mass in Cyg X-3 of 10.8 ± 3.4 M_☉ and 6.2 ±  2.0 M_☉ for Cyg X-3 distance of 9.3 and 7.2 kpc, respectively. In the above calculation, we used f_G = 1.0 as discussed below. Thus, over the full range of parameters and uncertainties considered we cannot justify the mass of the X-ray star being less than 4 solar masses. This seems to rule out the possibility that a neutron star is the compact object in Cyg X-3. Accepting the recent distance estimate of 9.0 kpc (Predehl et al. 2000) we estimate a compact object mass of 10.1 ± 3.2.

Zoom In Zoom Out Reset image size

We note that the index-normalization method is subject to several limitations that do not affect the index-quasi-periodic oscillation method. First, there is a dependence on the relative inclination angles of the binary system and that of the template system. Second, it is not always the case that self-similar patterns are seen in the data. For example, while our Cyg X-3 results were easily reconcilable with GRS 1915+105, the data could not be fit to an analogous template derived from XTE J1550–564. Finally, the relative source distances and geometry come into play, so that our mass determination is dependent on a factor f_G. In the case of a spherical coronal geometry, f_G is unity. This is a reasonable approximation for the case of low-hard and intermediate spectral-state BHs. The ST09 results confirm that for most cases, when the system inclinations are unknown, the f_G = 1.0 works very well. In the unlikely case of a flat disk geometry, f_G is the ratio of the inclination angle cosines between a reference source and the target source.