UNDERSTANDING COMPACT OBJECT FORMATION AND NATAL KICKS. III. THE CASE OF CYGNUS X-1

In recent years, the number of observed black hole (BH) X-ray binaries (XRBs) has grown significantly. For these binaries, there exists a wealth of observational information about their current physical state: BH and donor masses, orbital period, donor's position on the H-R diagram and surface chemical composition, transient or persistent X-ray emission, and Roche lobe overflow (RLO) or wind-driven character of the mass transfer (MT) process. Furthermore, proper motions have been measured for a handful of these binaries (e.g., Mirabel et al. 2001, 2002; Mirabel & Rodrigues 2003). Together with the earlier measurements of center-of-mass radial velocities and distances, we can obtain information about the three-dimensional kinematic properties of these binaries. Given this plethora of observational results, the current observed sample of BH XRBs provides us with a unique opportunity to understand the formation and evolution of BHs in binaries. This paper is the third in a series where we investigate in detail the BH formation in XRBs, especially focusing on the mass relationship between BHs and their immediate progenitors and the possible BH natal kick magnitude imparted during the core-collapse event.

In the first paper of this series, Willems et al. (2005, hereafter Paper I) showed how using the currently available constraints one could uncover the evolution history of an XRB from the present state back to the time just prior to the core-collapse event. They applied their analysis to the BH XRB GRO J1655-40. In the second paper, Fragos et al. (2009, hereafter Paper II) performed the same analysis for the case of the BH XRB XTE J1118+480. In this paper, we focus on the case of the BH XRB Cygnus X-1. The MT mechanism in Cygnus X-1 is different from the XRBs studied in previous cases. Both donors in GRO J1655-40 and XTE J1118+480 are transferring mass to the BH under RLO, whereas the BH in Cygnus X-1 is accreting mass from the stellar wind of its companion.

The plan of the paper is as follows. In Section 2, we review Cygnus X-1's currently available observational constraints. A general outline of the analysis used to reconstruct the system's evolutionary history is presented in Section 3, while individual steps of the analysis are discussed in more detail in Sections 4–7. In Section 8, we derive constraints on the formation of the BH. The final section is devoted to a summary of our results and discussion of some of the assumptions introduced in our analysis.

Cygnus X-1 was first detected in Aerobee surveys in 1964 by Bowyer et al. (1965). Soon after the discovery, it was identified as an XRB, which consisted of a compact object and a visible star HDE 226868 (Murdin & Webster 1971; Webster & Murdin 1972; Bolton 1972a). Spectroscopic observations led Walborn (1973) to classify HDE 226868 as an O9.7 Iab supergiant. Bregman et al. (1973) estimated the distance to be 2.5 kpc and set a lower limit of 1 kpc, based on the colors of field stars in the vicinity of the supergiant. Using a combination of data from David Dunlap Observatory and the Royal Greenwich Observatory, Bolton (1972b) derived the orbital period, eccentricity, and systemic radial velocity (V₀) to be 5.5995 ± 0.0009 days, 0.09 ± 0.02, and −6.0 ± 0.1 km s⁻¹, respectively. Based on the absence of X-ray and optical eclipses, the author gave a lower limit of 7.4 M_☉ on the mass of the compact object. This implied that the compact object was too massive to be a white dwarf or a neutron star. Thus, the author confirmed Webster & Murdin's (1972) finding that the compact object is a very strong BH candidate.

Using the orbital period obtained from spectrometry and a range in the assumed degree of Roche filling of the supergiant, Gies & Bolton (1982, 1986a) found a lower mass limit of 7 M_☉ for the compact object. This confirmed that the compact object observed in Cygnus X-1 was a BH. The same authors also refined the orbital period and eccentricity to be 5.59974 ± 0.00008 and 0.021 ± 0.013, respectively, and measured V₀ to be −2.0 ± 0.7 km s⁻¹. Ninkov et al. (1987) used the relationship between the equivalent width of the Hγ spectral line and the absolute magnitude of early-type supergiants to estimate the distance as 2.5 ± 0.3 kpc.

Herrero et al. (1995) performed a detailed spectroscopic analysis on the supergiant, and derived the masses to be 10.1 M_☉ and 17.8 M_☉ for the BH (M_BH) and the supergiant (M₂), respectively, if an orbital inclination angle of 35° was assumed. Using the Isaac Newton telescope, LaSala et al. (1998) measured the orbital period as 5.5997 ± 0.0001 days and V₀ as −5.4 ± 0.1 km s⁻¹. With all the accumulated radial velocity measurements and their own spectroscopy of the supergiant, Brocksopp et al. (1999) refined the orbital period to 5.599829 ± 0.000024. The proper motion of Cygnus X-1 was observed with the Very Large Baseline Interferometry (VLBI) between 1988 and 2001 (Lestrade et al. 1999; Stirling et al. 2001; Mirabel & Rodrigues 2003). During this period, the system's position shifted at a rate of −4.2 ± 0.2 mas yr⁻¹ in right ascension (R.A.) and −7.6 ± 0.2 mas yr⁻¹ in declination (decl.). Meanwhile, a trigonometric parallax of 0.73 ± 0.30 mas was also measured with VLBI, which gave a distance of 1.4^+0.9_{− 0.4} kpc (Lestrade et al. 1999).

By studying the spectra obtained with the 0.9 m coudé feed telescope of Kitt Peak National Observatory, the 2.1 m telescope of University of Texas McDonald Observatory, and the 1.9 m telescope of University of Toronto David Dunlap Observatory between 1998 and 2002, Gies et al. (2003) derived V₀ as −7.0 ± 0.5 km s⁻¹ and estimated M_BH/M₂ ≈ 0.36 ± 0.05. Caballero-Nieves et al. (2009) examined the supergiant's ultraviolet spectra from the Hubble Space Telescope. Their results gave masses of 23⁺⁸_{− 6} and 11⁺⁵_{− 3} M_☉ for the supergiant and the BH, respectively. On the other hand, Shaposhnikov & Titarchuk (2007) used the X-ray quasi-periodic oscillation and spectral index relationship and deduced M_BH to be 8.7 ± 0.8 M_☉, which overlapped with the lower end of the M_BH range derived by Caballero-Nieves et al. (2009).

Recently, Reid et al. (2011) measured the trigonometric parallax of Cygnus X-1 with the National Radio Astronomy Observatory's Very Long Baseline Array and found a distance of 1.86^+0.12_{− 0.11} kpc. The authors also reported proper-motion measurements of Cygnus X-1, which were −3.78 ± 0.06 mas yr⁻¹ in R.A. and −6.40 ± 0.12 mas yr⁻¹ in decl. Meanwhile, Xiang et al. (2011) studied the X-ray dust scattering halo of Cygnus X-1 and determined the distance to be 1.81 ± 0.09 kpc, after considering the compatibility with the parallax result. Building on the trigonometric parallax distance measurement of Reid et al. (2011), Orosz et al. (2011) performed optical data modeling of Cygnus X-1 and found the mass of the supergiant to be 19.2 ± 1.9 M_☉ and the BH to be 14.8 ± 1.0 M_☉. Using the results of Reid et al. (2011) and Orosz et al. (2011), Gou et al. (2011) determined that Cygnus X-1 hosts a near-extreme Kerr BH, with a spin parameter a_* > 0.95.

Unlike most of the XRBs known to host a BH, Cygnus X-1 is a persistent X-ray source. Since the supergiant is currently not overfilling its Roche lobe (Gies & Bolton 1986a), the observed X-rays are mainly powered by the accretion of stellar wind. The X-ray luminosity of Cygnus X-1 varies between two discrete levels, namely, the "hard (low) state" and the "soft (high) state." As the system spends most of its time (∼90%, see Cadolle Bel et al. 2006) in the hard state, we focus on the hard state X-ray luminosity (L_X). Frontera et al. (2001) observed Cygnus X-1 with the Narrow Field Instruments of the BeppoSAX satellite at different epochs in 1996. The authors obtained the L_X (0.5–200 keV) and the extrapolated bolometric luminosity (L_bol) as 2.0 × 10³⁷ and 2.4 × 10³⁷ erg s⁻¹, respectively, assuming a distance of 2 kpc. Using observational data obtained by the Compton Gamma Ray Observatory between 1991 and 2000, McConnell et al. (2002) derived L_bol to be (1.62–1.70) × 10³⁷ erg s⁻¹, with the distance to the source fixed at 2 kpc. Cadolle Bel et al. (2006) observed Cygnus X-1 with the International Gamma-Ray Astrophysics Laboratory between 2002 and 2004 and measured L_X (20–100 keV) as 6.5 × 10³⁶ erg s⁻¹, assuming a distance of 2.4 kpc. The authors also gave L_bol as 2.2 × 10³⁷ erg s⁻¹.

For the systemic parameters relevant to our analysis, we adopt the most recent observational constraints, with the exception of L_bol. We consider all the L_bol values mentioned above, assuming they represent the typical X-ray variability range for this system. After rescaling their values to the parallax distance measurement by Reid et al. (2011) and considering the uncertainty in that distance, we adopt L_bol to be (1.17–2.35) × 10³⁷ erg s⁻¹. For ease of reference, our adopted observational constraints are summarized in Table 1.

Table 1. Properties of Cygnus X-1

Parameter	Notation	Value	References
Distance (kpc)	d	1.86+0.12− 0.11	(9)
Galactic longitude (deg)	l	71.3	(2)
Galactic latitude (deg)	b	+3.1	(2)
Proper motion in R.A. (mas yr−1)	μR.A.	−3.78 ± 0.06	(9)
Proper motion in decl. (mas yr−1)	μdecl.	−6.40 ± 0.12	(9)
Systemic velocity (km s−1)	V0	−7.0 ± 0.5	(5)
Orbital period (days)	Porb	5.599829 ± 0.000016	(1)
Orbital eccentricity	eorb	0.018 ± 0.003	(8)
Inclination angle	i	27.06 ± 0.76	(8)
Black hole mass (M☉)	MBH	14.81 ± 0.98	(8)
Black hole spin	a*	>0.95	(10)
Companion mass (M☉)	M2	19.16 ± 1.90	(8)
Companion Radius (R☉)	R2	16.50 ± 0.84	(8)
Companion Luminosity (L☉)	L2	(1.91–2.75) × 105	(8)
Companion Effective temperature (K)	Teff	30,000–32,000	(8)
Companion surface rotation speed (km s−1)	Vrotsin i	95 ± 6	(7)
Bolometric luminosity of the X-ray source (erg s−1)	Lbol	$(1.3\hbox{--}2.1) (\frac{{\rm d}}{1.86\,{\rm kpc}} )^2 \times 10^{37}$	(3), (4), (6)

References. (1) Brocksopp et al. 1999; (2) Lestrade et al. 1999; (3) Frontera et al. 2001; (4) McConnell et al. 2002; (5) Gies et al. 2003; (6) Cadolle Bel et al. 2006; (7) Caballero-Nieves et al. 2009; (8) Orosz et al. 2011; (9) Reid et al. 2011; (10) Gou et al. 2011.

Download table as:  ASCIITypeset image

In our analysis, we assume that Cygnus X-1 formed in the Galactic disk, from the evolution of an isolated primordial binary at solar metallicity. In fact, Mirabel & Rodrigues (2003) suggest that Cygnus X-1 belongs to Cygnus OB3 (Cyg OB3), which is an OB association located close to the Galactic plane. We also assume that no MT via RLO occurred in the evolutionary history of this binary.

According to our current understanding, in order to form a ∼15 M_☉ stellar BH at solar metallicity, the BH progenitor in the primordial binary needs to be more massive than 120 M_☉ (Belczynski et al. 2010). Such a massive star loses its hydrogen rich envelope via stellar wind and exposes its naked helium core. At the end of nuclear evolution, it collapses into a BH. During the core-collapse event, the orbit is altered by the asymmetric mass loss from the system and a possible recoil kick imparted to the BH. If the binary survives through the core-collapse event, angular momentum loss via gravitational radiation and tidal effects causes the orbit to shrink, although wind mass loss leads to orbital expansion. In the meantime, the more evolved BH companion is losing mass via its own stellar wind at a higher rate. The system becomes a BH XRB when the BH captures a non-negligible amount of mass from its companion's stellar wind.

In this paper, we restrict ourselves to the formation of BH XRBs through the above evolutionary channel. Like the first two papers, our goal is to track the evolutionary history of Cygnus X-1 back to the time just prior to the core-collapse event. Our analysis incorporates a number of calculations which can be summarized in four steps.

First, we identify the current evolutionary stage of the BH companion, so that all the observational constraints are satisfied. Under the assumption that the BH companion mass has not been altered by MT in the past, we model it as an isolated star. Using a stellar evolution code, we calculate a grid of evolutionary sequences of isolated stars at different zero-age main-sequence (ZAMS) masses. We examine each sequence to find whether there exists a point in time that the calculated stellar properties, i.e., mass, radius, luminosity, and effective temperature, are all simultaneously in agreement with the currently observed properties of the BH companion. If such a period of time exists, we classify that sequence as "successful." The current age of the BH companion can be estimated from these successful sequences, and the time expired since the BH formation can then be derived by subtracting the approximate lifetime of the BH progenitor.

Next, we consider the kinematic evolutionary history of the XRB in the Galactic potential. Starting from the current location, we follow the methodology of Gualandris et al. (2005) and use the observed three-dimensional velocity to trace the Galactic motion of Cygnus X-1 backward in time. Together with the constraints on the current age of the system derived in the first step, this allows us to determine the location and velocity of the binary at the time of BH formation (we denote these as "birth" location and velocity). By subtracting the local Galactic rotational velocity at the "birth" location from the systems's center-of-mass velocity, we derive constraints on the peculiar velocity of the binary right after the formation of the BH.

In the third step, we analyze the orbital dynamics of the core-collapse event due to mass loss and possible natal kicks imparted to the BH. In this paper, we refer to the instants right before and after the formation of the BH by the terms "pre-SN" and "post-SN," respectively. We start with the constrained parameter space of (M_BH and M₂) derived in the first step and perform a Monte Carlo simulation scanning over the parameter space of the pre-SN binary properties. This parameter space is limited by requirements of orbital angular momentum and energy conservations, and by the post-SN binary peculiar velocity constraint derived in the second step. This calculation yields a population of simulated post-SN binaries for each successful sequence.

Finally, we follow the orbital evolution of these simulated binaries to the current epoch. Our calculation accounts for tides, wind mass loss, wind accretion onto the BH, and orbital angular momentum loss via gravitational radiation. At the end of the calculations, we require agreement between the observed and calculated orbital period and eccentricity.

Under the assumption that the companion mass has not been altered by MT in its past, we model the companion as an isolated star using a modified version of the stellar evolution code EZ (originally developed by Paxton 2004).

We calculate the evolution of our stellar models at solar metallicity, which is the same metallicity that Orosz et al. (2011) used in deriving the properties of the BH companion. When we place the companion's observational constraints on an H-R diagram, we find that the current location of the companion does not seem to be consistent with any evolutionary tracks calculated by the stellar evolution code. As shown in Figure 1, the companion is overluminous for a star of its mass. This cannot be explained by earlier MT from the BH progenitor to the companion. Braun & Langer (1995) studied the effects of mass accretion onto massive main-sequence stars and found that the accreting stars would not appear overluminous for their new masses during the rest of their main-sequence lifetime. If mass accretion leads to a so-called rejuvenation of the accreting star, which means its central hydrogen abundance substantially increases, it would have the same luminosity as a star of its new mass. If rejuvenation does not occur, the accreting star would appear underluminous for its new mass during the rest of its main-sequence lifetime.

Zoom In Zoom Out Reset image size

One possible solution for matching the observed companion's luminosity is increasing the core overshooting parameter α_ov to ∼0.45. Although this value is relatively high, it is not unphysical. Claret (2007) compared the data from 13 double-line eclipsing binary systems with theoretical predictions of stellar modeling and found α_ov could be as high as 0.6 for massive stars. We vary α_ov from 0.35 to 0.5, in steps of 0.01. We note that the need for such higher values of α_ov in the modeling of massive stars may very well be connected to the significant presence of internal rotation and associated rotational mixing. Effectively increasing α_ov leads to stronger internal mixing and in a way allows the stellar model to behave more like a rotating model.

Besides the observational constraints on the companion's properties, there are three additional constraints. The first one comes from the fact that the companion is currently not overfilling its Roche lobe (Gies & Bolton 1986a). Thus, we require the stellar radius R₂ in our models to be

$$\begin{equation} R_2 \le A_{\rm orb}r_{\rm Egg} + \Delta R, \end{equation} \tag{ 1 }$$

where r_Egg is the effective Roche lobe radius given by Eggleton (1983). Here, we make an approximation that the orbit is circular and synchronized. The parameter ΔR is a constant accounting for the difference in the calculated stellar radii among stellar evolution codes (Valsecchi et al. 2010). We set ΔR to 2.5 R_☉.

Another constraint is that the calculated bolometric luminosity (L_bol) resulting from the stellar wind accretion process needs to fall within the observational range, which is (1.17–2.35) × 10³⁷ erg s⁻¹. By adopting the Bondi & Hoyle (1944) accretion model and following Belczynski et al. (2008), the orbital-averaged accretion rate is given by

$$\begin{equation} \dot{M}_{{\rm acc}} = - \frac{F_{{\rm wind}}}{\sqrt{1-e_{{\rm orb}}^2}} \left(\frac{GM_{{\rm BH}}}{V_{{\rm wind}}^2} \right)^2 \frac{\alpha _{{\rm wind}}}{2A_{{\rm orb}}^2} \frac{\dot{M}_2}{(1+V^2)^{3/2}}. \end{equation} \tag{ 2 }$$

Here, M_BH is the BH mass in our models, which varies within the 1σ range of the observational constraint, in steps of 0.098 M_☉. Since the total mass that the BH could have accreted from its companion stellar wind is negligible, M_BH in each evolutionary sequence is fixed throughout our analysis. $\dot{M}_2$ is the wind mass-loss rate of the companion in our models. F_wind is a parameter such that $\dot{M}_{{\rm acc}}$ never exceeds $-0.8 \dot{M}_2$, and α_wind is the accretion efficiency, which varies between 1.5 and 2.0 (Boffin & Jorissen 1988). A_orb and e_orb are the orbital semimajor axis and eccentricity, respectively. A_orb is derived from the mean measured orbital period P_orb, which is

$$\begin{equation} A_{{\rm orb}} = \left[ {\frac{G \left(M_{{\rm BH}} + M_2 \right) P_{{\rm orb}}^2}{4\pi ^2}} \right] ^{\frac{1}{3}}, \end{equation} \tag{ 3 }$$

where M₂ is the companion mass in our models. e_orb is set equal to the mean measured orbital eccentricity. V_wind denotes the wind velocity. V² equals to V²_BH/V_wind², where V²_BH is the orbital velocity square of the BH and is approximated as G(M_BH + M₂)/A_orb. We adopt the spherically symmetric wind velocity law given in Lamers & Cassinelli (1999),

$$\begin{equation} V_{\rm wind}(r) \ = \ V_{\rm esc} \ + \ (V_\infty \ -\ V_{\rm esc})\left(1 \ - \ \frac{R_2}{r} \right)^\beta, \end{equation} \tag{ 4 }$$

where r is the distance from the companion to the BH and is set equal to A_orb. β is a free parameter varying from 0.6 to 1.6 (Gies & Bolton 1986b; Lamers & Leitherer 1993), in steps of 0.1. V_∞ is the wind velocity at infinity, while V_esc is the effective escape velocity at the surface of the companion. Within the typical range of O star surface temperature, V_∞ is scaled as 2.65V_esc (Kudritzki & Puls 2000). Following Lamers & Cassinelli (1999),

$$\begin{equation} V_{\rm esc} = \sqrt{2(1-\Gamma _e)GM_2/R_2}, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} \Gamma _e \ = \ \frac{\sigma _eL_2}{4\pi cGM_2} \end{equation} \tag{ 6 }$$

is the mass correcting factor for the radiative force due to electron scattering, and c is the speed of light in vacuum. Lamers & Leitherer (1993) scaled the electron scattering coefficient per unit mass σ_e as

$$\begin{equation} \sigma _e \ = \ 0.401\left(\frac{1\ +\ q\epsilon }{1\ + \ 3\epsilon } \right), \end{equation} \tag{ 7 }$$

where q is the fraction of He⁺⁺ and (1 − q) is the fraction of He⁺, with q = 1 if T_eff ⩾ 35, 000 K, q = 0.5 if 30,000 K ⩽ T_eff < 35,000 K, and q = 0 if T_eff < 30,000 K. The abundance ratio = He/(H + He) is fixed at 0.15, which is appropriate for an O star with a spectral type of Class I. Using $\dot{M}_{\rm acc}$ from Equation (2), we follow Belczynski et al. (2008) and calculate the bolometric luminosity resulting from the companion's stellar wind being accreted onto the BH as

$$\begin{equation} L_{\rm bol} = \frac{1}{2} \frac{GM_{\rm BH}\dot{M}_{\rm acc}}{R_{{\rm acc}}}, \end{equation} \tag{ 8 }$$

where R_acc denotes the radius of the accretor. For the case of BH, R_acc is the radius of the inner most stable circular orbit, which we calculate with Equation (2.21) in Bardeen et al. (1972). Given the observationally inferred spin a_* > 0.95 (Gou et al. 2011), we adopt the median a_* = 0.97 and find

$$\begin{equation} R_{\rm acc} = 2.57 \left(\frac{M_{\rm BH}}{M_\odot } \right) \,{\rm km}. \end{equation} \tag{ 9 }$$

This calculated L_bol needs to fall within the observational range.

Zoom In Zoom Out Reset image size

The last additional constraint is that the observational constraints on L₂ and L_bol have to be evaluated at the same distant estimation. To examine this, we calculate the ratio

$$\begin{equation} f_L \ = \ \left(\frac{L_2}{10^5 L_\odot } \right) \left(\frac{L_{{\rm bol}}}{10^{37} {\rm \,erg \,s^{-1}}} \right)^{-1}\,, \end{equation} \tag{ 10 }$$

which is independent of distance. From Figure 1 in Orosz et al. (2011), L₂ is 2.09 × 10⁵ L_☉ at T_eff = 30, 000 K, and is 2.51 × 10⁵ L_☉ at T_eff = 32, 000 K, assuming a distance of 1.86 kpc. Together with the measured range of L_bol rescaled at the same distance estimation, the upper and lower limits of f_L are 1.01 and 1.90, respectively. We can assure that both luminosity constraints are evaluated at the same distance estimation if f_L falls within that range.

In order to find the current evolutionary stage of the BH companion, we apply these constraints to a set of evolutionary sequences, which cover the parameter space of the companion's ZAMS mass (M_{2, zams}), α_ov, M_BH, α_wind, and β. For each sequence, we find whether there exists a point in time that the calculated properties simultaneously satisfy all observational constraints: the BH companion's mass, luminosity, temperature, and radius, L_bol, f_L, and not overfilling the Roche lobe of the BH companion. Similar to the Roche lobe constraint, we also consider an uncertainty of ±2.5 R_☉ in the calculated stellar radii when we apply the observational constraint of the BH companion's radius. If such a period of time exists, we classify that evolutionary sequence as "successful." The behavior of some relevant parameters is illustrated in Figure 2 for two selected successful sequences that are chosen mainly to provide a clear and instructive picture. The displayed sequences therefore do not represent our best possible matches to the observed properties of Cygnus X-1. Figure 3 shows the parameter space of M_{2, zams}, α_ov, and M_BH covered by all successful sequences. For α_wind and β, the successful sequences covered the entire allowed parameter space, which are 1.5 ⩽ α_wind ⩽ 2.0 and 0.6 ⩽ β ⩽ 1.6.

Zoom In Zoom Out Reset image size

The current age of the BH companion could be derived from the time interval at which all observational constraints are satisfied. Assuming that the BH progenitor and its companion formed at the same time, we could compute the time since the BH formation (t_sys) by

$$\begin{equation} t_{\rm sys} = t_2 - t_{\rm BH}, \end{equation} \tag{ 11 }$$

where t_BH is the approximate lifetime of the BH progenitor. We follow Belczynski et al. (2010) to calculate M_BH and t_BH for different progenitors using the stellar evolution code SSE (Hurley et al. 2000), and adopting the mass-loss prescriptions which were classified as "Vink et al. winds." The calculated t_BH are fit as a function of M_BH,

$$\begin{equation} \frac{t_{\rm BH}}{10^6 \,{\rm yr}} = \frac{\frac{M_{\rm BH}}{M_\odot }}{19.26 - 4.902 \left(\frac{M_{\rm BH}}{M_\odot } \right) + 0.3841 \left(\frac{M_{\rm BH}}{M_\odot } \right)^2} + 3.341, \end{equation} \tag{ 12 }$$

for M_BH ⩾ 9.5 M_☉. Figure 4 shows the variations of t₂ and t_BH against M_{2, zams}. We find that t_sys is between 4.8 and 7.6 Myr.

Zoom In Zoom Out Reset image size

Here, we assume that Cygnus X-1 formed in the Galactic disk. The consideration of Cyg OB3 being the parent association of Cygnus X-1 is discussed in Section 9.2. Given the observed position and measured proper motion of Cygnus X-1, we derive the post-SN peculiar velocity of the binary's center of mass by tracing its orbit in the Galaxy back to the time of BH formation. We describe the motion of the binary with respect to a right-hand Cartesian reference frame, whose origin coincides with the Galactic center. The Z-axis points to the northern Galactic pole, while the X-axis points in the direction from the projected position of the Sun onto the Galactic plane to the Galactic center. In this reference frame, the Sun is located at (X_☉, Y_☉, Z_☉) = (−8.5, 0, 0.03) kpc (Joshi 2007; Ghez et al. 2008; Gillessen et al. 2009; Reid et al. 2009) and has a peculiar motion (U_☉, V_☉, W_☉) = (11.1, 12.24, 7.25) km s⁻¹(Schönrich et al. 2010). Cygnus X-1 is currently located at a distance of 1.86^+0.12_{− 0.11} kpc from the Sun, with a Galactic longitude l = 713, and a Galactic latitude b = 31 (Lestrade et al. 1999; Reid et al. 2011). This means Cygnus X-1 is currently ∼130 kpc above the Galactic plane.

To model the Galaxy, we adopt the Galactic potential of Carlberg & Innanen (1987) with updated model parameters of Kuijken & Gilmore (1989). The equations governing the system's motion in the Galaxy are integrated backward in time, up to the time corresponding to the current system's age t_sys given by the successful sequences. We follow the methodology of Gualandris et al. (2005) to initialize the parameters for the integration, which accounts for the uncertainties in the estimated distance and measured velocity components. We generate the initial system's position by the Galactic coordinates (l, b) and a random distance drawing from a Gaussian distribution. We generate initial system's velocity by drawing randomly the proper motions (μ_R.A., μ_decl.) and heliocentric radial velocity (V₀) from Gaussian distributions. The current system's age is uniformly distributed between 4.8 and 7.6 Myr (see Figure 4).

Figure 5 shows the possible positions of Cygnus X-1 at the time of BH formation, obtained from integrating 3000 trajectories backward in time. As there is no trajectory crossing the Galactic plane and the end points of all trajectories fall within 110 pc from the Galactic plane, we consider each end point as a possible birth site of the BH. The post-SN peculiar velocity V_{pec, postSN} of the binary is obtained by subtracting the local Galactic rotational velocity from the center-of-mass velocity of the binary at the birth sites. We find V_{pec, postSN} ranges from 22 to 32 km s⁻¹ and is time independent. The distribution of V_{pec, postSN} against the time expired since the formation of the BH are displayed in Figure 5.

Zoom In Zoom Out Reset image size

For each successful sequence, we perform a Monte Carlo simulation which consists of 20 million pre-SN binaries. The properties of the BH progenitor's companion are taken from the stellar model of that sequence, at the time when the age of the star is equal to t_BH. During a supernova (SN) explosion, the mass loss from the system and possibly the kick imparted to the BH change the binary's orbital parameters. The pre- and post-SN component masses, orbital semimajor axis, and orbital eccentricity are related by the conservation laws of the orbital energy and angular momentum. In the following, we add the subscripts "preSN" and "postSN" to the notation of the orbital elements to distinguish between their values just prior to, and right after, the SN explosion that formed the BH.

We start with seven free parameters: the BH immediate (He-rich) progenitor mass (M_He), pre-SN orbital semimajor axis (A_preSN) and eccentricity (e_preSN), the mean anomaly (m), the magnitude (V_k), and direction (θ, ϕ) of the kick velocity imparted to the BH. θ is the polar angle of the kick with respect to the relative orbital velocity of the BH progenitor just prior to the SN explosion, and ϕ is the corresponding azimuthal angle (see Figure 1 in Kalogera 2000 for a graphical representation). The first five parameters are drawn from uniform distributions, while the last two are drawn from isotropic distributions. It is obvious that the progenitor must of course be more massive than the BH, but there is no absolute upper limit for the progenitor mass. We adopt M_He ⩽ 20 M_☉, and provide a discussion on this upper limit in Section 9.1.

The relations between pre- and post-SN parameters have been derived by Hills (1983):

$$\begin{eqnarray} V^2_k + V^2_{\rm He,preSN} + 2V_kV_{\rm He,preSN}\cos \theta &=& G(M_{\rm BH}+M_2)\nonumber\\ &&\times\,\left(\frac{2}{r} - \frac{1}{A_{\rm postSN}} \right),\nonumber\\ \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&G (M_{\rm BH} + M_2) A_{\rm postSN}\big(1-e^2_{\rm postSN}\big)\nonumber\\ &&\quad =\, r^2 \big(V^2_k \sin ^2 \theta \cos ^2 \phi + [ \sin \psi (V_{\rm He,preSN} + V_k \cos \theta) \nonumber \\ &&\qquad -\, V_k \cos \psi \sin \theta \sin \phi]^2 \big). \end{eqnarray} \tag{ 14 }$$

Here, r is the orbital separation between the BH progenitor and its companion at the time of SN explosion,

$$\begin{equation} r = A_{\rm preSN}(1-e_{\rm preSN}\cos E_{\rm preSN}), \end{equation} \tag{ 15 }$$

where E is the eccentric anomaly and is related to m as

$$\begin{equation} m=E-e \sin E. \pagebreak \end{equation} \tag{ 16 }$$

V_{He, preSN} is the relative pre-SN orbital velocity of the BH progenitor,

$$\begin{equation} V_{\rm He,preSN} = \left[ G(M_{\rm He}+M_2)\Bigl (\frac{2}{r} - \frac{1}{A_{\rm preSN}}\Bigr) \right]^{1/2}. \end{equation} \tag{ 17 }$$

The angle ψ is the polar angle of the position vector of the BH with respect to its pre-SN orbital velocity in the companion's frame. It is related to the pre-SN parameters as

$$\begin{equation} r^2V^2_{\rm He,preSN} \sin ^2 \psi = G(M_{\rm He} + M_2) A_{\rm preSN} \big(1 - e^2_{\rm preSN}\big). \end{equation} \tag{ 18 }$$

Since the core collapse is instantaneous, r remains unchanged. This gives a constraint

$$\begin{eqnarray} r & = A_{\rm preSN}(1 - e_{\rm preSN}\cos E_{\rm preSN}) \nonumber \\ & =A_{\rm postSN}(1 - e_{\rm postSN}\cos E_{\rm postSN}), \end{eqnarray} \tag{ 19 }$$

which needs to be satisfied with |cos E_postSN| ⩽ 1.

The mass loss from the system and a natal kick imparted to the BH can induce a post-SN peculiar velocity (V_{pec, postSN}) at the binary's center of mass. Its magnitude is determined by Equations (28)–(32) in Paper I and is required to fall within the range derived in Section 5, which is 22–32 km s⁻¹.

In addition, there are two more restrictions on the properties of pre- and post-SN binary components. First, we require that both components have to fit within their pre- and post-SN Roche lobe at periapsis. We impose this condition to avoid complications arising from MT-induced changes in the stellar structure of the MS companion that later becomes the BH companion of the XRB. To calculate the Roche lobe radius of each component in eccentric pre- and post-SN orbits, we adopt the fitting formulae of Sepinsky et al. (2007). When calculating the pre-SN Roche lobe radii, we assume that the pre-SN orbit is pseudo-synchronized. Again, due to the difference in calculated stellar radii among stellar evolution codes, we consider an uncertainty of ±2.5 R_☉ on the companion radius (Valsecchi et al. 2010). The radius of the BH immediate progenitor can be approximated by Equations (3) in Fryer & Kalogera (1997), since we assume that it is a Helium star. Second, the pre-SN spin of the BH immediate progenitor and its companion need to be less than the breakup angular velocity Ω_c ≈ (GM/R³)^1/2. As the calculated stellar radius R₂ associates with an uncertainty ΔR = 2.5 R_☉ (Valsecchi et al. 2010),

$$\begin{equation} \Omega _c = \sqrt{\frac{GM_2}{R_2^3} } \left(1 + \frac{3}{2} \frac{\Delta R}{R_2}\right) \end{equation} \tag{ 20 }$$

for the BH companion.

The orbital evolution of the simulated binaries, which are generated from the Monte Carlo simulations described in Section 6, is calculated up to the current epoch. After the formation of the BH, the orbital parameters of the binary are subject to secular changes due to the tidal torque exerted by the BH on its companion, and due to the loss of orbital angular momentum via gravitational radiation and stellar wind. Since the tidal interactions depend on both the orbital and rotational properties of the MS companion, the star's rotational angular velocity (Ω) right after SN explosion that formed the BH enters the problem as an additional unknown quantity. Here, we assume the rotational angular velocity of the BH companion is unaffected by the SN explosion and is pseudo-synchronized to the pre-SN orbital frequency. The system of equations governing the tidal evolution of the orbital semimajor axis A, eccentricity e, and the BH companion's rotational angular velocity Ω has been derived by Hut (1981):

$$\begin{eqnarray} \left(\frac{dA}{dt} \right)_{\rm tides} &= & {-} 6 \frac{k_2}{T} \frac{M_{\rm BH}}{M_2} \frac{M_{\rm BH} + M_2}{M_2} \left(\frac{R_2}{A}\right)^8 \frac{A}{(1-e^2)^{15/2}}\nonumber \\ && \times\, \left[ f_1(e^2) - (1-e^2)^{3/2}f_2(e^2) \frac{\Omega }{n}\right], \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \left(\frac{de}{dt} \right)_{\rm tides} &=& {-} 27 \frac{k_2}{T} \frac{M_{\rm BH}}{M_2} \frac{M_{\rm BH} + M_2}{M_2} \left(\frac{R_2}{A}\right)^8 \frac{e}{(1-e^2)^{13/2}}\nonumber \\ && \times\, \left[ f_3(e^2) - \frac{11}{18}(1-e^2)^{3/2}f_4(e^2) \frac{\Omega }{n}\right], \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \left(\frac{d\Omega }{dt} \right)_{\rm tides} &= & 3 \frac{k_2}{T} \left(\frac{M_{\rm BH}}{M_2} \right)^2 \frac{M_2R_2^2}{I_2} \left(\frac{R_2}{A} \right)^6 \frac{n}{(1-e^2)^6}\nonumber \\ && \times\, \left[ f_2(e^2) - (1-e^2)^{3/2}f_5 (e^2) \frac{\Omega }{n} \right]. \end{eqnarray} \tag{ 23 }$$

Here, k₂ and I₂ are the apsidal-motion constant and moment of inertia of the MS companion, respectively. T is a characteristics timescale for the orbital evolution due to tides and n = 2π/P_orb is the mean orbital angular velocity. The coefficient functions f_i(e²) for i = i, 2, ⋅⋅⋅, 5 are given in Equations (11) in Hut (1981). As the BH companion in Cygnus X-1 is a massive MS star with a radiative envelope, the factor k₂/T can be approximated as

$$\begin{eqnarray} \left(\frac{k_2}{T} \right)_{\rm rad} &=& 1.9782 \times 10^4 \left(\frac{R_2}{R_\odot } \right) \left(\frac{R_\odot }{A} \right)^{5/2} \left(\frac{M_2}{M_\odot } \right)^{1/2}\nonumber \\ && \times\, \left(\frac{M_{\rm BH}+M_2}{M_2} \right)^{5/6} E_2 \;{\rm yr^{-1}}. \end{eqnarray} \tag{ 24 }$$

The constant E₂ comes from a fit to the tables in Claret (2004),

$$\begin{equation} {\rm log} E_2 = - \frac{t / t_{{\rm ms}}}{2.20489 - 1.89579 (t/t_{{\rm ms}})} - 5.51039, \end{equation} \tag{ 25 }$$

for 15.85 ⩽M_{2, zams} ⩽ 25.12 M_☉. Here, t_ms is the main-sequence lifetime. We define the end of the main sequence as the hydrogen abundance at the core being less than 0.01.

To follow the secular changes of the orbital parameters associated with emissions of gravitational waves, we adopt Equations (35) and (36) in Junker & Schäfer (1992), which are derived up to 3.5 post-Newtonian order.

The rates of change in A and e due to wind mass loss and wind accretion onto the BH are determined by following Equations (15) and (16) in Hurley et al. (2002),

$$\begin{eqnarray} \left(\frac{dA}{dt} \right)_{\rm wind} &=& {-}A \left[ \frac{\dot{M}_2}{M_{\rm BH}+M_2} + \left(\frac{2-e^2}{M_{\rm BH}} + \frac{1+e^2}{M_{\rm BH}+M_2}\right)\right.\nonumber\\ &&\times\, \left.\frac{\dot{M}_{\rm acc}}{1-e^2} \right] \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} \left(\frac{de}{dt} \right)_{\rm wind} = & {-}e\dot{M}_{\rm acc} \left(\frac{1}{M_{\rm BH}+M_2} + \frac{1}{2M_{\rm BH}}\right). \end{eqnarray} \tag{ 27 }$$

The mass loss via stellar wind also induces a loss in the spin angular momentum of the BH companion. Hurley et al. (2000) showed that if all the mass is lost uniformly from a thin shell at the surface of the MS star,

$$\begin{equation} \dot{J}_{\rm 2,spin} = \frac{d}{dt} \left(I_2 \Omega \right) = \frac{2}{3}\dot{M}_2 R^2_2 \Omega, \end{equation} \tag{ 28 }$$

where J_{2, spin} is the spin angular momentum of the BH companion.

For each of simulated binaries, we follow the secular changes of its orbital properties due to all the mechanisms mentioned in this section. The properties of binary components are adopted from the corresponding successful sequence. Unlike finding V_{pec, postSN} in Section 5, the orbital evolution of the binary goes forward in time, from t_BH to t₂. Within this period of time, the BH companion has to always fit within its Roche lobe at periapsis. In other words, its calculated radius is constrained to be less than the Roche lobe radius at periapsis given by Sepinsky et al. (2007). Again, we allow an uncertainty of ±2.5 R_☉, due to the difference in calculated stellar radii among stellar evolution codes (Valsecchi et al. 2010). Furthermore, the rotational angular velocity of the BH companion has to be smaller than the breakup angular velocity Ω_c. If the orbital period and eccentricity of the simulated binary at t₂ match the measured values of Cygnus X-1, we classify that binary as a "winning binary." Figure 6 illustrates the time evolution of orbital parameters for one selected winning binary and shows that the change in the semimajor axis is mainly determined by the stellar wind mass loss, while the change in the eccentricity is overwhelmingly dominated by the tidal effects.

Zoom In Zoom Out Reset image size

The elements presented in the previous sessions can now be combined to establish a complete picture of the evolution of Cygnus X-1 and the dynamics involved in the core-collapse event that formed the BH. After finding the successful evolutionary sequences that satisfy all the observed properties of the BH companion and the bolometric X-ray luminosity, as discussed in Section 4, we trace the motion of the system in the Galaxy back in time to the formation of the BH. We adopt the methodology of Gualandris et al. (2005) to account for the uncertainties in the measured distance and velocity components of Cygnus X-1. The time of BH formation is different for each successful sequence. It is estimated by the BH mass of the sequence, which connects to an approximate lifetime of the corresponding BH progenitor. This procedure gives us a constraint on the system's peculiar velocity right after the BH formation. We then perform Monte Carlo simulations on the orbital dynamics at core collapse for each successful sequence. There are seven free parameters: the BH immediate progenitor mass, the pre-SN orbital semimajor axis and eccentricity, the mean anomaly, the magnitude of kick velocity imparted to the BH, and the two angles specifying the direction of the kick velocity. The Monte Carlo simulations produce a population of simulated binaries, which satisfy the post-SN system's peculiar velocity constraint derived already. Lastly, we evolve the orbits of these simulated binaries forward in time to the current epoch. If the orbital period and eccentricity of the simulated binary at current epoch match the measured values of Cygnus X-1, we classify that simulated binary as a "winning binary." The results presented in what follows are derived from the winning binaries of all successful sequences.

In Figure 7, we present the probability distribution functions (PDFs) of the BH immediate (He-rich) progenitor mass (M_He) and natal kick magnitude (V_k). We find M_He to be in a range of 15.0–20.0 M_☉, and V_k to be ⩽77 km s⁻¹, both at 95% confidence. Figure 8 illustrates the two-dimensional joint V_k–M_He confidence levels, which shows that if M_He is less than ∼17 M_☉, the BH might have received a non-zero natal kick at the core-collapse event. For small M_He, a minimum V_k of ∼55 km s⁻¹ is necessary for explaining the current observed properties of Cygnus X-1. Furthermore, both the M_He PDF and the two-dimensional joint V_k–M_He confidence levels show that the maximum M_He is constrained by our adopted upper limit of 20 M_☉. We impose this limit based on the physics involved in the evolution of massive stars. A discussion on this limit can be found in Section 9.1. Given our understanding of mass loss from Helium stars, it seems that the BH has potentially received a small natal kick velocity of ⩽77 km s⁻¹ (95% confidence) during the core-collapse event.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on the dynamical model of Orosz et al. (2011), Gou et al. (2011) found that the BH in Cygnus X-1 has a spin parameter a_* > 0.95 at 3σ. To determine whether the BH was born with an extreme spin, we first estimate how much mass the BH could have accreted from its companion's stellar wind since the time of BH formation. The winning binaries of all successful sequences show that at maximum the BH has accreted ∼2 × 10⁻³ M_☉. Since it is impossible to spin the BH up to a_* > 0.95 by accreting that negligible amount of mass, the BH needs to have an extreme spin at birth. This high spin has implications about BH formation and the role of rotation in core collapse. Axelsson et al. (2011) also concluded that the observed spin connects to processes involved in core collapse and is not likely to originate from the synchronous rotation of the BH progenitor.

Besides the constraints on the BH formation, our results also shed light on the evolutionary picture of Cygnus X-1. We find that right after the formation of the BH, the BH companion has a mass of 19.8–22.6 M_☉, in an orbit with period of 4.7–5.2 days. Since then, the orbital separation of Cygnus X-1 has been increasing with time, as the rate of change in the semimajor axis is dominated by the influence of stellar wind mass loss from the system. On the other hand, the orbital eccentricity has decreased slightly since the BH formation. This is because the tides exerted on the companion by the BH, as the dominant mechanism of circularizing the orbit, are not strong enough to decrease the orbital eccentricity significantly within the time period of several million years since the time of BH formation. We find that e_postSN ranges from 0.015 to 0.022. However, this does not suggest that e_preSN has to be small. An eccentric pre-SN orbital could become fairly circular if there is a natal kick imparted to the BH at the right direction. As illustrated in Figure 9, there are winning binaries with e_preSN being as high as ∼0.53.

Zoom In Zoom Out Reset image size

In this paper, we constrained the progenitor properties and the formation of the BH in the persistent XRB Cygnus X-1. Our analysis accounts for the orbital evolution and motion through the Galactic potential right after the BH formation, and the binary orbital dynamics at the time of core collapse. We find that the mass of the BH immediate progenitor falls within a range of 15.0–20.0 M_☉ at 95% confidence. We note that the maximum progenitor mass is constrained by our adopted upper limit, which is discussed in Section 9.1. The BH has potentially received a small natal kick velocity of ⩽77 km s⁻¹ at 95% confidence. In fact, if the progenitor mass is less than ∼17 M_☉, a non-zero natal kick velocity is necessary to explain the currently observed properties of Cygnus X-1. Since the BH has only accreted mass from its companion's stellar wind, the total amount of mass accreted since the BH formation is less than ∼2 × 10⁻³ M_☉. This indicates that the observationally inferred BH spin of a_* > 0.95 (Gou et al. 2011) cannot be explained by mass accretion and has to be natal. This high spin has implications for BH formation and the role of rotation in core collapse. Right after the BH formation, the BH companion has a mass of 19.8–22.6 M_☉, in an orbit with period of 4.7–5.2 days and eccentricity of 0.015–0.022. Although the post-SN orbital eccentricity is small, the pre-SN orbit can potentially be fairly eccentric. This is possible if the BH receives a natal kick velocity at the right magnitude and direction.

The formation of the BH in Cygnus X-1 has been previously studied by Nelemans et al. (1999) and Mirabel & Rodrigues (2003). Both studies assumed symmetric mass loss during the core-collapse event and considered only the binary orbital dynamics at the time of core collapse. Comparing with these two earlier studies, we consider the possible asymmetries developed during the core-collapse event and the evolution of the binary since the BH formation. It is important to note that these two earlier studies do not consider the multitude of the observational constraints taken into account here and hence the suggested progenitors are not complete solutions for the evolutionary history of Cygnus X-1.

Finally, we discuss some of the assumptions introduced in our analysis in the following subsections.

9.1. Maximum BH Progenitor Mass

9.2. Association with Cyg OB3

The center of Cyg OB3 locates at l = 728 and b = 20, and at a distance of 1.4–2.7 kpc away from the Sun (Massey et al. 1995; Dambis et al. 2001; Mel'Nik et al. 2001; Mel'Nik & Dambis 2009). When comparing that to the location of Cygnus X-1 (Table 1), it is clear that not only are their Galactic coordinates close to each other, but also their distance estimations overlap with each other. Furthermore, the measurements of proper motion and radial velocity show that Cygnus X-1 is moving as the members of Cyg OB3 (Dambis et al. 2001; Mirabel & Rodrigues 2003; Mel'Nik & Dambis 2009). Based on these observations, Mirabel & Rodrigues (2003) argue that Cyg OB3 is the parent association of Cygnus X-1. This infers that V_{pec, postSN} due to the core-collapse event has to be small. If we change the constraint on V_{pec, postSN} to ⩽10 km s⁻¹, we find M_He to be in a range of 13.9–16.9 M_☉ and V_k to be ⩽24 km s⁻¹, both at 95% confidence. Besides the change in 95% limits, non-zero BH natal kicks are not needed for progenitors of M_He ⩽ 17 M_☉ in order to explain the observed properties of Cygnus X-1, but become necessary for M_He > 17.5 M_☉ (see Figure 10). Also, we note that a relatively small change on the range of V_{pec, postSN} affects the derived constraint on V_k qualitatively.

Zoom In Zoom Out Reset image size

9.3. Super-synchronized Orbit

After considering several previous measurements of the BH companion's surface rotation speed (V_rotsin i), Caballero-Nieves et al. (2009) adopted V_rotsin i = 95 ± 6 km s⁻¹. Orosz et al. (2011) found that the ratio of the BH companion's spinning frequency to the orbital frequency (f_Ω) was 1.400 ± 0.084, which was derived based on their results of the inclination angle i = 2706 ± 076 and the companion radius R₂ = 16.5 ± 0.84. This indicates that the BH companion is super-synchronized. We note that with the analysis presented here, we find that none of our winning binaries have super-synchronized BH companions at the current epoch. They are all sub-synchronized with f_Ω reaching ∼0.87 at maximum.

In an effort to examine how our standard assumptions can be modified and to investigate whether super-synchronism is at all allowed by the models as indicated by the observations, we make two modifications to our analysis. We first remove the assumption that the pre-SN orbit is pseudo-synchronized, and randomly distribute the pre-SN spin of the BH companion between zero and its breakup angular frequency Ω_c. Next, we reduce the secular changes of the orbital parameters due to the influence of tides by multiplying the right-hand side of Equations (21)–(23) by a constant f_tide.

As shown in Figure 11, by allowing the pre-SN spin of companion to be greater than pseudo-synchronization and keeping the tidal strength unchanged (i.e., f_tide = 1.0), the maximum f_Ω of the winning binaries increases to ∼1.2. Although it is getting close, this value is still below the observationally inferred one. Together with a weakened tidal strength of f_tide = 0.2 and 0.5, we comfortably find winning binaries with f_Ω = 1.4. Furthermore, the minimum pre-SN surface rotation speed of the companion in those winning binaries are ∼500 and 700 km s⁻¹ for f_tide = 0.2 and 0.5, respectively. Given the uncertainties in the physics of fast rotating massive stars, it seems that super-synchronism in Cygnus X-1 is allowed by our presented models, if the companion is spinning faster than orbital pseudo-synchronization right before the core-collapse event and the tides exerted on the companion are weaker than the nominal theoretical values.

Zoom In Zoom Out Reset image size

9.4. Constraints on the BH Kick Direction

This work has been primarily supported by the NSF Grant AST-0908930; V.K. also acknowledges partial support through the NSF Grant PHY-1066293 and the hospitality of the Aspen Center for Physics; T.F. acknowledges Fellowship support by the Harvard-Smithsonian Center for Astrophysics and the Harvard Institute for Theory and Computation.