UNDERSTANDING COMPACT OBJECT FORMATION AND NATAL KICKS. IV. THE CASE OF IC 10 X-1

Over the past few decades, it has become clear that neutron stars (NS) receive recoil kicks at birth (also known as natal kicks) during the core collapse event. This conclusion is based on proper motion studies of pulsars (see, e.g., Gunn & Ostriker 1970; Lyne et al. 1982; Lyne & Lorimer 1994; Brisken et al. 2003; Hobbs et al. 2005; Chatterjee et al. 2009) and evolutionary studies of NS-hosting binaries (see, e.g., Brandt & Podsiadlowski 1995; Pfahl et al. 2002; Thorsett et al. 2005; Willems et al. 2006; Martin et al. 2009; Wong et al. 2010). However, whether black holes (BHs) receive similar natal kicks during the core collapse event is still uncertain (see, e.g., Brandt et al. 1995; Nelemans et al. 1999; Gualandris et al. 2005; Dhawan et al. 2007; Repetto et al. 2012). If BH kicks are ubiquitous, then BH formation must be closely associated to that of NS before the formation of the event horizon. Otherwise, the formation of BHs through more than one physical process will be favored.

To shed light on questions related to BH formation, we perform detailed evolutionary modeling of the observed BH X-ray binaries (XRBs), which enables us to derive robust constraints on the mass of the BH immediate progenitor and the magnitude of the natal kick imparted to BH during the core collapse event. In the previous three papers of this series, we studied the XRB GRO J1655–40 (Willems et al. 2005), XTE J1118+480 (Fragos et al. 2009), and Cygnus X-1 (Wong et al. 2012). For GRO J1655–40, we constrained the mass of the BH immediate progenitor to ≃ 5.5–11.0 M_☉ and the magnitude of its natal kick to be ≲ 210 km s⁻¹. In the study of J1118+480, we found that a natal kick of magnitude ≃ 80–310 km s⁻¹ is required to explain the formation of the system, and also derived a lower limit of ≃ 6.5 M_☉ on the mass of the BH immediate progenitor. Finally, for Cygnus X-1 we constrained the mass of the BH immediate progenitor and the magnitude of its natal kick to be ≳ 15 M_☉ and ≲ 77 km s⁻¹, respectively. Similarly, Valsecchi et al. (2010) studied the formation of M33 X-7 using binary modeling and derived the mass of the BH immediate progenitor and the magnitude of the natal kick to be in the range of ≃ 15.0–16.1 M_☉ and ≃ 10–850 km s⁻¹, respectively. In this paper, we study the formation of the BH in the extragalactic XRB IC 10 X-1.

IC 10 X-1 is one of the four observed BH XRBs that are known to host a Wolf–Rayet (W-R) star as the mass donor (Clark & Crowther 2004). The other three systems are Cygnus X-3 (see, e.g., Zdziarski et al. 2013), NGC 300 X-1 (see, e.g., Crowther et al. 2010), and M 101 ULX-1 (see, e.g., Liu et al. 2013). The X-ray emission is powered by the accretion of stellar wind material onto the BH. At the present time, the BH in IC 10 X-1 is the most massive known stellar-mass BH (≃ 23–34 M_☉; Prestwich et al. 2007; Silverman & Filippenko 2008). Because the supergiant progenitor of the observed W-R star cannot fit into the tight orbit at the present time (orbital period = 34.93 hr; Silverman & Filippenko 2008), it is natural to consider the system's evolution via a common envelope (CE) evolution phase, which involves the BH and the massive progenitor of the observed W-R star. However, such a binary is likely to merge at the end of the CE phase, as the envelopes of massive stars are tightly bound (Podsiadlowski et al. 2003).

De Mink et al. (2009) suggested an alternative formation scenario for IC 10 X-1 that does not invoke any CE phase, which they called "Case M" evolution. They considered two massive stars in a tight orbit, such that tidal interactions always kept them spinning rapidly. This led to an efficient mixing of elements throughout their stellar interiors via rotational effects, and thus both stars went through chemical homogenous evolutions. They stayed compact throughout their main sequence evolution and turned into abnormally massive helium (He) stars without going through any CE evolution. Although Case M helps to explain the high masses of the BH and the W-R star in IC 10 X-1, the short orbital period gives rise to a difficulty: the intense mass loss suffered by the W-R star and its progenitor will widen the orbit. It is very hard to explain the current tight orbit without a CE phase. As a result, the evolution history of IC 10 X-1 has remained uncertain.

Instead of studying its past evolution, Bulik et al. (2011) performed binary modeling to predict the fate of IC 10 X-1. They estimated that the observed W-R star will go through core collapse in ≲ 0.3 Myr, leading to the formation of a close double BH binary with a short coalescence time (∼3 Gyr).

The plan of the paper is as follows. In Section 2, we review the current available observational constraints of IC 10 X-1. A general outline of our analysis methodology is presented in Section 3, while detailed discussions of each individual step are in Sections 4–6. Our derived constraints related to the formation of the BH and the past evolution of IC 10 X-1 are discussed in Section 7, and we offer our conclusions in the final section.

IC 10 X-1 is a persistent X-ray source in the local starburst galaxy IC 10. It was discovered by Brandt et al. (1997) from X-ray observations of IC 10 made with ROSAT. Bauer & Brandt (2004) derived the 0.5–8.0 keV unabsorbed X-ray luminosity of IC 10 X-1 to be 1.50 × 10³⁸ erg s⁻¹ from their Chandra observations. The optical counterpart of IC 10 X-1 was identified as the W-R star [MAC92] 17A by Clark & Crowther (2004). Based on the Schaerer & Maeder (1992) mass-luminosity relationship, they estimated the mass of this W-R star to be 35 M_☉. Using the data from Chandra and Swift observations, Prestwich et al. (2007) determined the orbital period of IC 10 X-1 to be 34.40 ± 0.83 hr and obtained a mass function of 7.8 M_☉. Using their mass function and considering an inclination angle of 90 degrees, they estimated the mass of the BH to be 23–34 M_☉. Their results indicated that IC 10 X-1 hosted the most massive stellar mass BH known at the present time. The observed orbital period also implied that the BH is currently accreting mass from the intense stellar wind of its W-R companion. Silverman & Filippenko (2008) precisely measured the radial velocity amplitude of the W-R star and refined the orbital period and mass function to 34.93 ± 0.04 hr and 7.64 ± 1.26 M_☉, respectively. For convenience, our adopted observational constraints are summarized in Table 1.

Among the evolutionary scenarios that could potentially explain the formation of IC 10 X-1, we adopt the simplest possible evolutionary history that provides a consistent explanation of all observational constraints. First, we assume that the BH progenitor and its companion were born at the same time. Toward the end of the BH progenitor's life it lost its hydrogen (H) rich envelope because of mass loss via a stellar wind or binary interactions. Hence, the BH immediate progenitor is a helium (He) rich star. Soon after the birth of the BH, the companion star evolved off the main sequence and became a supergiant. Eventually, the star overfilled its Roche lobe and underwent a phase of dynamically unstable mass transfer, which inevitably led to CE evolution. During the CE phase, the BH was engulfed into the H-rich envelope of its companion. The CE phase ended when the envelope was ejected, leaving behind a binary consisting of a BH and an He star in a tight, circular orbit. The intense mass loss via stellar wind from the He star continuously drove the binary components further away from each other. A very small fraction of this stellar wind material was accreted onto the BH, resulting in the X-ray emission from IC 10 X-1 that we see today.

This is the simplest possible evolutionary channel one can envision given the current properties of the system. Therefore, in this paper, we restrict ourselves to the formation of IC 10 X-1 through the above evolutionary channel. In order to derive constraints on the formation of the BH in IC 10 X-1, we track the evolutionary history of the system back to the instant just before the core collapse event. Our analysis incorporates a number of calculations that can be grouped in three main steps. Hereafter, we add the prefix "pre-" and "post-" to the name of any event that occurred in the evolutionary history of IC 10 X-1, in order to indicate the instant just before and right after that event, respectively.

In the first analysis step, our goal is to derive the post-CE binary and stellar properties of IC 10 X-1. We start by constructing He star models that satisfy the current mass and luminosity constraints given in Table 1. Using the properties of our He star models and different BH masses and post-CE binary semi-major axes, we evolve the post-CE binary's orbit forward in time until the current epoch. This calculation accounts for tides, wind mass loss, wind accretion onto the BH, and orbital angular momentum loss via gravitational radiation. We examine the evolutionary sequence of every binary model to find out if, at the present time, the BH mass and the orbital period simultaneously match the observational data, in which case we classify that evolutionary sequence as successful. Furthermore, we consider 1σ of observational uncertainties when matching the properties of our models with the observational data. The post-CE binary properties can then be obtained from our successful evolutionary sequences.

In the second analysis step, we study the CE event to determine whether our derived properties of the post-CE binary can be achieved. To do so we employ the standard α-prescription (Webbink 1984), also examining several alternative formulations of the energy budget prescription, to compute the CE efficiency, α_CE. If this parameter is ⩽1, we conclude that our derived properties of the post-CE binary can be explained by the current understanding of CE evolution.

In the final analysis step, we utilize the results from the first two steps and perform a Monte Carlo simulation of the orbital dynamics involved in the core collapse event. Our goal is to derive limits and construct probability distribution functions (PDF) for the BH immediate progenitor mass and the potential natal kick magnitude imparted to the BH. We start with randomized properties of the pre-core collapse (pre-supernova (SN)) binary, which include the BH immediate progenitor mass, orbital semi-major axis, eccentricity, and phase. The core collapse event is approximated as occurring instantaneously, with the mass ejection and possible asymmetries in the explosion imparting to the BH a potential kick that has random magnitude and direction. Using the equations of orbital energy and angular momentum, we determine whether the binary can survive the core collapse event. If it does, we map the pre-SN binary properties to the post-core collapse (post-SN) phase space. Then, we apply the constraint related to CE evolution and retain the data points with α_CE ⩽ 1. Finally, our constraints related to the BH formation in IC 10 X-1 are derived from the retained data points.

To build the stellar models necessary for our analysis, we use a version of the STARS code originally developed by Peter Eggleton (Eggleton 1971, 1973; Pols et al. 1995; Eggleton & Kiseleva-Eggleton 2002; Glebbeek et al. 2008). The adopted opacity tables combine the OPAL opacities from Iglesias & Rogers (1996), the low temperature molecular opacities from Ferguson et al. (2005), the electron-conduction opacities from Cassisi et al. (2007), and the Compton scattering opacities from Buchler & Yueh (1976). The assumed heavy element composition in each stellar model is scaled according to the solar abundances described by Anders & Grevesse (1989). Convective boundaries are determined using the Schwarzschild criterion. The mixing of elements due to convection and semi-convection are taken into account (Eggleton 1973; Langer 1991), but the effects of rotational mixing and meridional circulation are excluded. We notice that these effects can be very important. However, their efficiencies are poorly constrained and need to be tuned using the observational data. We choose not to add extra free parameters because our models are intended to be the simplest ones able to satisfy the observational constraints. For mass loss during the main sequence evolution, we adopt the prescriptions of Vink et al. (2001). If the surface H abundance drops below 0.4, we switch to the W-R prescription developed by Nugis & Lamers (2002) with the metallicity scaling law determined by Vink & de Koter (2005). When the effective temperature falls below the working regime (i.e., <10 kK) of the mass loss prescriptions developed by Vink et al. (2001), we use the mass loss rates determined by de Jager et al. (1988) instead. This stellar-wind mass loss recipe is extrapolated into the post-main sequence evolution, including the regime of luminous blue variables (LBV). Indeed, all the above mass loss prescriptions give a high mass loss rate in the LBV regime, with a range of 10⁻⁵–10⁻² M_☉ yr⁻¹ for the massive stars involved in this analysis. This range of rates is consistent with the measured mass loss rates of LBVs (see, e.g., Stahl et al. 1990, 2001; Leitherer et al. 1994; Hillier et al. 2001; Vink & de Koter 2002; Umana et al. 2005; Mehner et al. 2013). In addition, there is no prescription included to simulate the mass loss that is occurred in giant eruptions of LBV stars, nor to account for the Humphreys & Davidson (1979) limit phenomenon. Stellar evolutionary codes generally approximate this episodic mass loss with the averaged continuous mass loss. We note that the mass loss rates of massive stars, especially LBVs, are poorly constrained at the present time. The effect of this uncertainty on our study of the required CE evolution will be discussed in Section 5.3.

To model the post-CE phase of IC 10 X-1, we start by modeling the observed W-R star. We construct He star models with different initial masses at the measured metallicity for IC 10 (Z = 0.004, see Table 1) using the stellar evolutionary code described in Section 3. Then, we evolve the models and retain those that, at a certain time during their evolution, have a luminosity and mass in agreement with the observations listed in Table 1 (see Figure 1). Hereafter, we refer to these He star models as successful. At this stage, we use the observational constraint on the W-R mass, rather than effective temperature. This is because the effective temperature shown in Table 1 is the spectroscopic temperature, which is generally unreliable for W-R stars (see Crowther 2007).

Zoom In Zoom Out Reset image size

Next, we use the properties of the successful He star models to follow the evolution of the post-CE binary orbit and the He star's spin until the present time. Our goal is to constrain the post-CE binary properties. For each successful He star model, we vary the post-CE orbital semi-major axis in steps of 0.1 R_☉ and consider different BH masses according to the measured mass function (f(M_BH)) and inclination angle (i) listed in Table 1. Specifically, we take the mass of each successful He star model at the current epoch and compute the maximum and minimum BH mass at present (M_{BH, max} and M_{BH, min}, respectively) using the measured f(M_BH) and i with 1σ of observational uncertainties. Then, we define

$$\begin{equation} M_{{\rm BH, mean}} \equiv \frac{M_{{\rm BH, max}}+ M_{{\rm BH, min}}}{2}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \Delta M_{{\rm BH}} \equiv \frac{M_{{\rm BH, max}} - M_{{\rm BH, mean}}}{2}, \end{equation} \tag{ 2 }$$

and vary the BH mass at present from M_{BH, min} to M_{BH, max}, in steps of 0.2 · ΔM_BH. The amount of mass accreted onto the BH since its birth is negligible, as the Bondi & Hoyle (1944) accretion of the stellar-wind material leads to a very small capturing fraction. Hence, the post-CE BH mass is only slightly different from the one at present. For the post-CE configuration of the binary orbit and the He star's spin, we assume that CE evolution circularizes and synchronizes the binary orbit. We note that our assumption that the post-CE He star's spin is synchronized with the orbital frequency is completely arbitrary. However, this arbitrary choice only affects the strength of the tide exerted on the He star by the BH. This tide is weak in general, as the orbit of the BH–He star binary in this study is not very tight. From post-CE to current epoch, the He star is always far from filling its Roche lobe, with the Roche lobe radius being roughly three to nine times larger than the He star's radius. Thus our arbitrary assumption of the post-CE He star's spin does not have any significant influence on our analysis. Under these assumptions, we follow the evolution of the binary orbital semi-major axis and eccentricity, and the He star's spin, accounting for tides exerted on the He star by the BH, stellar wind mass loss, orbital angular momentum loss due to gravitational radiation, and accretion of the He star's wind material onto the BH. In order to determine how much wind material is accreted onto the BH, we compute the accretion rate $\dot{M}_{{\rm acc}}$. Specifically, we adopt the Bondi & Hoyle (1944) accretion model and follow the formalism in Sections 4.1 and 4.2 of Belczynski et al. (2008) to obtain $\dot{M}_{{\rm acc}}$.

The relevant ordinary differential equations (ODEs) governing the orbital evolution are integrated forward in time. To do this, we use the orbital evolution code described in the supplementary information of Valsecchi et al. (2010) with the following two modifications.

• 1.  
  For the second-order tidal coefficient E₂, we take a stellar model from Claret (2005) with an initial mass of 25.2 M_☉ and metallicity of 0.004 (both in agreement with the observations shown in Table 1), and derive

  $$\begin{eqnarray} {\rm log}(E_{2}) &= &-5.49491 \nonumber \\ &&- 1.94284 t_{{\rm MS}}^{51.277} -1.99707 t_{{\rm MS}}^{2.41139}, \end{eqnarray} \tag{ 3 }$$

  where t_MS is the star's evolutionary time expressed in units of the main sequence lifetime. As we are dealing with an He star, t_MS is taken to be the time the star has spent burning He at its center.
• 2.  
  For the effect of wind mass loss on the spin of the He star, we approximate that the wind carries away its angular momentum from a thin shell at the star's surface. Hence, we set

  $$\begin{equation} \dot{J}_{{\rm spin}} = \frac{2}{3}\dot{M_{2}}R_{2}^2\omega, \end{equation} \tag{ 4 }$$

  where J_spin is the spin angular momentum, R₂ is the stellar radius, and ω is the spin angular frequency of the He star.

For each combination of the post-CE binary component masses and orbital semi-major axis, the integration of the relevant ODEs proceeds forward in time only if the He star is spinning slower than the break-up frequency $\Omega _c \approx \sqrt{G M/R^3}$ and is not filling its Roche lobe.

At the end of the calculation, we retain the evolutionary sequences that at the present time match the observed orbital period. Hereafter, we refer to them as successful evolutionary sequences. The X-ray luminosity of each evolutionary sequence at the present time is computed by following the formalism in Section 9.1 of Belczynski et al. (2008). Using the relation between the radius of the accretor and the unknown BH spin, the X-ray luminosity is given by

$$\begin{equation} L_{X} = \eta _{{\rm bol}} \epsilon \frac{\dot{M}_{{\rm BH}}c^{2}}{r_{{\rm isco}}}, \end{equation} \tag{ 5 }$$

where c is the speed of light, and η_bol and are adopted to be 0.8 and 0.5, respectively. The variable r_isco is the radius of the inner most stable circular orbit around the BH expressed in units of GM_BH/c². It equals six for a non-spinning BH and one for a maximally spinning BH. Figure 2 shows that for every successful evolutionary sequence the observed X-ray luminosity always falls between the upper (r_isco = 1) and lower (r_isco = 6) limits of our predicted L_X at the present time. Thus, we are not formally imposing the observed L_X as a constraint, but our predicted L_X at the present time for all successful evolutionary sequences are naturally consistent with the observed values.

Zoom In Zoom Out Reset image size

The post-CE binary component masses and the orbital semi-major axis of IC 10 X-1 given by all successful sequences are illustrated in Figure 3. The figure shows that the post-CE binary consists of a 25–39 M_☉ BH and a 32–42 M_☉ He star and has an orbital separation of 17–22 R_☉. The equivalent orbital period is 25–35 hr. Due to the intense stellar wind suffered by the He star, the orbital separation increases continuously. Such a binary evolves to become the observed XRB IC 10 X-1 in ⩽0.4 Myr.

Zoom In Zoom Out Reset image size

Prior to the Monte Carlo simulation of the orbital dynamics involved in the core collapse event, we study whether our derived post-CE binary properties of IC 10 X-1 are achievable based on the current understanding of CE evolution. We construct stellar models of the pre-CE He star progenitor candidates, which have core masses covering the entire range of the post-CE He star's mass given by the successful evolutionary sequences. These models are built based on the stellar evolution of isolated stars. We note that the BH progenitor could have potentially transferred mass to its companion before the core collapse event. However, we assume the He star progenitor will quickly adjust itself after the end of the mass transfer. Thus, its structure right before the CE phase is not expected to be significantly different from an isolated star of the same mass. Furthermore, in order to initiate CE evolution, the immediate progenitor of the He star needs to fill its Roche lobe and start an unstable mass transfer before its radius reaches its maximum value at the instant t_Rmax. If this is not the case, the star will shrink rapidly afterward due to intense stellar wind mass loss. Meanwhile, most of this wind material will leave the binary system, resulting in the widening of the binary orbit. Because of the increase in orbital separation, the Roche lobe of the star expands. Since the star is shrinking and its Roche lobe is expanding, it cannot overfill its Roche lobe and initiate CE evolution after t_Rmax. Thus, we only consider the properties of our He star progenitor models before t_Rmax when studying the CE event.

To examine whether CE evolution can explain our derived post-CE properties for IC 10 X-1, we use the energy formalism (Webbink 1984) and compute the corresponding CE efficiency α_CE,

$$\begin{equation} \alpha _{{\rm CE}} \cdot \Delta E_{{\rm orb}} = E_{{\rm bind}}, \end{equation} \tag{ 6 }$$

where ΔE_orb is the change in orbital energy and E_bind is the energy required for dispersing the envelope to infinity. Here, ΔE_orb is simply

$$\begin{equation} \Delta E_{{\rm orb}} = -\frac{GM_{{\rm BH}}M_{{\rm WRpro}}}{2A_{{\rm preCE}}} + \frac{GM_{{\rm BH}}M_{{\rm WR}}}{2A_{{\rm postCE}}}, \end{equation} \tag{ 7 }$$

where M_WRpro and M_WR are the masses of the pre-CE He star progenitor and the post-CE He star, and A_preCE and A_postCE are the pre- and post-CE orbital semi-major axis, respectively. Since the outcome of the CE phase is given by the successful evolutionary sequences, the second term in Equation (7) is predetermined. Thus, it is obvious that ΔE_orb reaches its maximum when A_preCE approaches ∞. Depending on the post-CE parameters given by each successful evolutionary sequence, maximum ΔE_orb ranges from 7.3 × 10⁴⁹ to 1.6 × 10⁵⁰ erg (see Figure 4).

Zoom In Zoom Out Reset image size

When computing E_bind of our He immediate progenitor models, we set the core-envelope boundary at the base of the convective thick hydrogen burning shell. Since this is generally located close to where the hydrogen abundance (X_H) is 10%, we simply define the core boundary at X_H = 0.1. Dewi & Tauris (2000) also considered the same definition for core boundary in their studies of CE evolution. Our choice of core boundary definition is justified in Section 5.2.

Based on our adopted energy formalism, shown as Equation (6), α_CE ⩽ 1 means that the loss in the orbital energy is sufficient to unbind the envelope of our He star's immediate progenitor model. In order words, this means that CE evolution is capable of producing our derived post-CE properties of IC 10 X-1. On the other hand, α_CE > 1 means that there is not enough energy to unbind the envelope and the pre-CE binary will eventually merge instead. We examine five different treatments of CE evolution: the original Webbink (1984) prescription (Section 5.1), the Webbink (1984) prescription with enhanced convective overshooting (Section 5.2) or increased mass-loss rates (Section 5.3), allowing hyper-critical accretion onto the BH during CE evolution (Section 5.4), and the "enthalpy" formalism (Ivanova & Chaichenets 2011; Section 5.5). We find that only the "enthalpy" formalism is capable of explaining the post-CE binary properties given by all successful evolutionary sequences with α_CE ⩽ 1. In the following subsections, we will discuss each of the aforementioned treatments of CE evolution.

5.1. Original Webbink (1984) Prescription

First, we use the standard definition of E_bind, which can be written as

$$\begin{equation} E_{{\rm bind}} = -\int ^{{\rm surface}}_{{\rm core\;boundary}} (\Phi (m) + \epsilon (m))\; dm. \end{equation} \tag{ 8 }$$

Here Φ(m) = −Gm/r is the gravitational potential and (m) is the specific internal energy, which only includes the thermal energy of the plasma gas and does not include the recombination energy of H and He, or the association energy of H₂ (Han et al. 1994, 1995). For our progenitor models, the sum of the recombination and association energies is less than a thousandth of the thermal energy of the plasma gas. Hence, E_bind will not change significantly by including these energy sources or not. In Figure 5, we illustrate the typical variation of E_bind throughout the evolution of an He star progenitor. After the termination of the main-sequence at t_tms, E_bind increases as the He star progenitor shrinks. It reaches a maximum shortly after the beginning of convective thick H shell burning at t_Hsb. Then, it starts to decrease because of the envelope expansion. Figure 5 also shows that throughout its evolution prior to t_Rmax, the E_bind of an He star progenitor is always the smallest at t_Rmax.

Zoom In Zoom Out Reset image size

Using the above definition of E_bind and Equation (6), we search for He star progenitor models that will give us α_CE ⩽ 1 at a certain time during its evolution before t_Rmax. To achieve this goal, we create a grid of stellar models at the observed metallicity of IC 10 X-1 by varying the initial mass from 65 to 105 M_☉ in steps of 0.1 M_☉. This grid covers the entire initial mass range of massive stars, which can have core masses fall within the mass range of the post-CE He stars constrained by all successful evolutionary sequences (i.e., 32–42 M_☉).

For each model in our grid, we match the core mass at any time before t_Rmax with the post-CE He star masses of the successful evolutionary sequences. Then, we compute ΔE_orb with Equation (7) using the relevant post-CE binary properties given by those matched sequences. Using this ΔE_orb and the E_bind defined by Equation (8), we can obtain α_CE with Equation (6). We find that for all models α_CE > 1, as E_bind is always several times larger than ΔE_orb. To illustrate this, we select six representative models of He star progenitors from our grid and list them in Table 2 as Model I–VI. Here, we consider the CE evolution involving these models at t_Rmax, which is the time when these models have the lowest E_bind. The mass range of the cores in these models at t_Rmax is the same as that of the post-CE He stars given by all successful evolutionary sequences. In Column 8, we list the maximum ΔE_orb during the CE evolution involving these models, under the assumption that A_preCE approaches ∞ (see Figure 4). It is obvious that for these models E_bind is at least 3.5 times larger than the corresponding maximum ΔE_orb. Thus, α_CE is always >1.

Table 2. Properties of Selected He Star Progenitor Models

Model	Initial Mass	αova	MWRprob	Mcorec	RWRprod	Ebinde	Max. ΔEorbf
	(M☉)		(M☉)	(M☉)	(R☉)	(1050 erg)	(1050 erg)
I	80.4	0.12	65.9	32.2	1406	4.05	1.14
II	84.1	0.12	67.3	34.1	1337	4.33	1.23
III	88.6	0.12	69.9	36.1	1268	4.72	1.32
IV	92.6	0.12	72.9	38.2	1214	5.13	1.42
V	95.9	0.12	75.6	40.0	1175	5.48	1.50
VI	99.9	0.12	78.2	41.6	1136	5.85	1.58
VII	75.0	0.20	60.0	32.2	1429	3.84	1.14
VIII	78.5	0.20	62.4	34.1	1343	4.14	1.23
IX	83.0	0.20	66.3	36.1	1262	4.63	1.32
X	86.8	0.20	68.2	38.2	1205	4.85	1.42
XI	89.9	0.20	70.9	40.0	1171	5.16	1.50
XII	93.6	0.20	73.0	41.6	1140	5.46	1.58
XIII	67.0	0.30	53.8	32.2	1403	3.49	1.14
XIV	70.9	0.30	56.9	34.1	1317	3.82	1.23
XV	74.8	0.30	59.7	36.1	1239	4.13	1.32
XVI	78.6	0.30	62.4	38.2	1188	4.40	1.42
XVII	81.8	0.30	65.1	40.0	1149	4.69	1.50
XVIII	84.8	0.30	67.0	41.6	1118	4.92	1.58

Notes. The models with α_ov = 0.12 have normal strength of convective core overshooting. The parameters listed in Columns 4–7 are the properties of the He star progenitor models when their radii reach their maximum values at the instant t_Rmax. ^aConvective overshooting parameter. ^bMass of the He star immediate progenitor. ^cCore mass of the He star immediate progenitor (same as the mass of the descendent He star right after CE evolution). ^dRadius of the He star immediate progenitor. ^eEnergy required to disperse the envelope of the He star immediate progenitor. ^fMaximum change in orbital energy during CE evolution involving the He star progenitor at t_Rmax (see Figure 4).

Download table as:  ASCIITypeset image

5.2. Enhanced Convective Overshooting

A possible reason for the negative result in the previous trial could be that the envelopes of our He star progenitor models are too massive. Tauris & Dewi (2001) showed that the envelope mass and E_bind vary significantly with different definitions of core boundary. To check whether there is a definition of core boundary that can decrease the E_bind of our He star progenitor models, we follow Tauris & Dewi (2001) and consider these definitions of core boundary: energy production rate (max _nuc), binding energy profile (Han et al. 1994), mass-density gradient (∂²log ρ/∂m² = 0; Bisscheroux 1998), and specific entropy profile. We find that the E_bind of our He-star progenitor models resulted from these core boundary definitions are similar to those obtained from our canonical choice (X_H = 0.1), with a difference of <3%. In Figure 6, we use Model I listed in Table 2 as an example to illustrate how similar the locations of the core boundary given by different definitions are.

Zoom In Zoom Out Reset image size

Alternatively, we can reduce the envelope mass by increasing the strength of convective core overshooting. In this trial, we increase the convective core overshooting (α_ov) parameter from the canonical value of 0.12 to 0.2 and 0.3. Our canonical α_ov = 0.12 leads to an overshooting length of 0.32 H_p, where H_p is the pressure scale height (see Schroder et al. 1997). We again compute α_CE using Equation (6) with E_bind defining as Equation (8).

As before, we make grids of He star progenitor models at the observed metallicity of IC 10 X-1 with α_ov = 0.2 and 0.3, by varying the initial mass between 60–95 and 55–88 M_☉ in steps of 0.1 M_☉, respectively. Each grid covers the entire initial mass range of massive stars, which can have core masses fall within the mass range of the post-CE He stars constrained by all successful evolutionary sequences (i.e., 32–42 M_☉).

After computing α_CE in the same way as our previous trial (see Section 5.1) for all the He star progenitor models in our grids of α_ov = 0.2 or 0.3, we do not find any cases of α_CE ⩽ 1. To help explain this result, we select a group of representative models from each grid and list them as Model VII–XII (α_ov = 0.2) and Model XIII–XVIII (α_ov = 0.3) in Table 2. Here, we consider the CE evolution involving these models at t_Rmax, which is the time when these models have the lowest E_bind (see Figure 5). For each group of models, their core masses cover the entire range of post-CE He star's mass given by the successful evolutionary sequences. The maximum ΔE_orb during the CE evolution involving these models (see Figure 4) are listed in Column 8. When comparing the E_bind of these models to that of the models with the normal strength of convective overshooting (i.e., Model I–VI), we find that for the models with the same core mass E_bind only decreases slightly with the increasing α_ov. Considering the difference between E_bind and maximum ΔE_orb of these models, that amount of decrease in E_bind is not enough to make α_CE ⩽ 1.

5.3. Enhanced Mass Loss Rates

As we have mentioned in Section 3, the mass loss rates for massive stars are not well constrained. The uncertainty in these rates not only can affect the envelope mass of the He star progenitor that needs to be ejected in CE evolution, but also the density structure of the envelope and hence E_bind (see, e.g., Podsiadlowski et al. 2003). To study the variance in E_bind due to this uncertainty, we compute He star progenitor models with all considered mass loss rates being enhanced by a factor of two. When comparing these models to our canonical models (α_ov = 0.12) according to the core mass at t_Rmax, we find that the E_bind of these models is 3%–14% smaller. However, this decrease in E_bind is not sufficient to make α_CE ⩽ 1.

On the other hand, there is also uncertainty in the mass loss rate for W-R stars. If it is higher than what we adopted when modeling the post-CE binary (see Section 4), the post-CE orbit could be tighter than what we found. This is because the post-CE orbital evolution is dominated by wind mass loss from the system, which always increases the orbital separation. To get an idea of how much this uncertainty can change the post-CE orbital separation, and hence the available orbital energy for ejecting the envelope in CE evolution, we let the mass loss rate for W-R stars to be a factor of two higher than what we considered. Then, the new mass of the post-CE He star ($M^{\prime }_{{\rm He,postCE}}$) can be written as

$$\begin{equation} M^{\prime }_{{\rm WR,postCE}} = 2M_{{\rm WR,postCE}} - M_{{\rm WR,now}}, \end{equation} \tag{ 9 }$$

where M_{WR, now} is the observed mass of the He star in IC 10 X-1 and M_{WR, postCE} is the mass of the post-CE He star derived from our post-CE binary modeling. Using the equation of Jeans-mode mass loss (see, e.g., Belczynski et al. 2008), we can write the new post-CE orbital separation ($A^{\prime }_{{\rm postCE}}$) as

$$\begin{equation} \frac{A^{\prime }_{{\rm postCE}}}{A_{{\rm postCE}}} = \frac{M_{{\rm BH}} + M_{{\rm WR,postCE}}}{M_{{\rm BH}}+M^{\prime }_{{\rm WR,postCE}}}. \end{equation} \tag{ 10 }$$

With Equation (7) and taking the limit that A_preCE tends to ∞, we can express the new maximum change of orbital energy during CE evolution (max. $\Delta E^{\prime }_{{\rm orb}}$) as

$$\begin{equation} max.\;\Delta E^{\prime }_{{\rm orb}} = \frac{M^{\prime }_{{\rm WR,postCE}}}{M_{{\rm WR,postCE}}} \frac{A_{{\rm postCE}}}{A^{\prime }_{{\rm postCE}}} \cdot {\rm max.}\;\Delta E_{{\rm orb}}. \end{equation} \tag{ 11 }$$

Using Equations (9)–(11), we find that doubling our adopted mass loss rates for W-R stars will lead to a <45% increase in the maximum change of orbital energy during CE evolution. This increase is not enough to make α_CE ⩽ 1, as the E_bind of our He star progenitor models are at least three times larger than the original max. ΔE_orb (see Table 2).

We show that the required α_CE will still be >1 even if we increase our adopted mass loss rates by a factor of two. However, we note that the uncertainties in these mass loss rates are larger than what we consider, especially for the stars that evolved off the main sequence. For instance, when studying the mass transfer in massive binaries, Petrovic et al. (2005) allowed a factor of six uncertainty in their mass loss rates for W-R stars. Given these large uncertainties, we can potentially obtain α_CE < 1 by adopting mass loss rates at even higher values. Since it is numerically challenging to evolve very massive stars with extraordinarily high mass loss rates beyond the main sequence, we choose to seek an alternative CE treatment that will naturally give α_CE < 1.

5.4. Hyper-critical Accretion

It has been suggested that a compact object might accrete a significant amount of mass after being engulfed into the envelope of its companion, due to hyper-critical accretion (Blondin 1986; Chevalier 1989, 1993; Brown 1995). Because only part of the companion's envelope will be dispersed to infinity, the E_bind in the energy formalism shown as Equation (6) needs to be adjusted accordingly. Hence, we write the equation of energy balance as

$$\begin{equation} \alpha _{{\rm CE}} \cdot \Delta E_{{\rm orb}} = \Delta E_{{\rm bind}} = f_{{\rm ej}} \cdot E_{{\rm bind}}, \end{equation} \tag{ 12 }$$

where

$$\begin{equation} f_{{\rm ej}} \equiv \frac{M_{{\rm env}} - \Delta M_{{\rm BH}}}{M_{{\rm env}}} \end{equation} \tag{ 13 }$$

is the fraction of the envelope mass ejected to infinity. Here, M_env is the mass of the envelope and ΔM_BH is the amount of mass accreted onto the BH. Based on Equations (12) and (13), we follow Belczynski et al. (2002) and derive the rates of change in the BH mass and the binary semi-major axis during the phase of hyper-critical accretion. The detailed derivation of these rates with respect to the mass of the He star progenitor (M_com) during the accretion phase can be found in the Appendix. Our rates are different from those derived by Belczynski et al. (2002), since we consider a fraction instead of the total E_bind when balancing the energy budget of the envelope ejection. In other words, we use Equation (12) instead of Equation (6) to incorporate α_CE.

Using the He star progenitor models in our constructed grids (see Sections 5.1 and 5.2) and considering their properties at any time prior to t_Rmax, we numerically integrate Equations (A10) and (A11) from M_com = M_core to M_WRpro with different α_CE between 0 and 1. Here, we assume that α_CE is constant throughout the whole CE evolution, including the hyper-critical accretion phase. From this integration, we can obtain the required radius of the He star progenitor and the BH mass at the onset of the hyper-critical accretion phase (i.e., when the BH touches its surface). Then, using Equation (A16), we can compute the corresponding pre-CE binary semi-major axis and justify whether it is physically possible.

We find that the BH could have accreted ∼10 M_☉ during the hyper-critical accretion phase. However, the radius of the He star progenitor at the onset of the hyper-critical accretion phase is required to be at least 1.3 times larger than that at the beginning of CE evolution (i.e., when the He star progenitor fills its Roche lobe). According to Equation (A16), a pre-CE binary semi-major axis leading to such condition does not exist. Therefore, we conclude that the hyper-critical accretion formalism cannot explain the post-CE binary properties given by the successful evolutionary sequences.

5.5. "Enthalpy" Formalism

Ivanova & Chaichenets (2011) argued that enthalpy should be considered when calculating the binding energy of the envelope and introduced

$$\begin{equation} E^{\prime }_{{\rm bind}} = -\int ^{{\rm surface}}_{{\rm core\;boundary}} \left(\Phi (m) + \epsilon (m) + \frac{P(m)}{\rho (m)} \right)\; dm. \end{equation} \tag{ 14 }$$

Using Model III, listed in Table 2, as an example, we illustrate the behavior of its envelope binding energy computed by Equations (8) and (14) on Figure 7. It is clear that $E^{\prime }_{{\rm bind}}$ is always smaller than E_bind during the evolution of Model III. This is because the term P(m)/ρ(m) is always positive within a star. Including it as an energy source lowers the binding energy of the envelope, and thus the required CE efficiency. Figure 7 also shows that $E^{\prime }_{{\rm bind}}$ at t_Rmax is smaller than the corresponding maximum ΔE_orb during the CE phase (see Figure 4 and Table 2). According to Equation (6), α_CE could be ⩽1 if ΔE_orb during the CE evolution involving this progenitor model is close to that maximum value.

Zoom In Zoom Out Reset image size

To examine whether the "enthalpy" formalism can explain our derived post-CE properties of IC 10 X-1, we consider the models of He star progenitors in our constructed grid with the normal strength of convective overshooting (α_ov = 0.12, see Section 5.1). Throughout the evolution of each model, we match its core mass with the masses of the post-CE He stars given by the successful evolutionary sequences. Then we compute α_CE with Equations (6) and (14), using the post-CE binary properties of the matched sequences. As we expected, we find many cases with α_CE ⩽ 1.

To summarize, the post-CE binary properties given by the successful evolutionary sequences can be explained by our understanding of CE evolution, using the standard α_CE prescription (Webbink 1984) with the "enthalpy" formalism (Ivanova & Chaichenets 2011). Since the "enthalpy" formalism can account for the formation of IC 10 X-1, we suggest that this strengthens the evidence in favor of the "enthalpy" formalism being an appropriate description of the energy budget for envelope ejection during CE evolution, at least for the massive stars involved in the formation of IC 10 X-1. We will discuss alternative CE treatments in Section  8, which are not considered in this study but might also be able to explain the existence of IC 10 X-1 without the need to invoke α_CE > 1. In addition, we find that the core boundaries of our He star progenitor models are roughly unchanged when considering alternative boundary definitions that existed in the literature.

Using the progenitor properties obtained from the previous steps (see Sections 4 and 5) as constraints, we perform Monte Carlo simulations on the binary orbital dynamics involved in the core collapse event. Our goal is to derive the constraints on the properties of the BH immediate progenitor and the magnitude of the natal kick imparted to the BH.

Just before the core collapse event, the binary consists of the BH immediate progenitor and its companion star, with an orbital semi-major axis A_preSN and eccentricity e_preSN. Note that since it is not necessary for the pre-SN progenitor to have experienced any mass-transfer phase, we cannot assume that the pre-SN orbit is circular. Instead, we consider the full range of possibilities with eccentric pre-SN orbits. The masses of the BH immediate progenitor and its companion are M_BHpro and M₂, respectively. During the core collapse event, the mass loss from the BH immediate progenitor and the potential natal kick imparted to the BH alter the binary orbital properties. Hence, the post-SN orbital semi-major axis and eccentricity become A_postSN and e_postSN, respectively. As we assume the companion star is not affected by the instantaneous core collapse event, its properties remain unchanged.

Using the equation of binary orbital energy and angular momentum, the pre- and post-SN binary properties are related as (Hills 1983; Wong et al. 2012):

$$\begin{eqnarray} V^2_k &+& V^2_{{\rm BHpro}} + 2 V_k V_{{\rm BHpro}} \cos \theta _k \nonumber \\ && = G(M_{{\rm BH}}+M_2)\left(\frac{2}{r} - \frac{1}{A_{{\rm postSN}}} \right) \end{eqnarray} \tag{ 15 }$$

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

Here, V_k is the magnitude of the natal kick imparted to the BH, while θ_k and ϕ_k describe its direction in the frame of the BH immediate progenitor. Specifically, θ_k is the polar angle of the natal kick with respect to the relative orbital velocity of the BH immediate progenitor, and ϕ_k is the corresponding azimuthal angle (see Figure 1 in Kalogera 2000, for a graphic representation). The variable r, which is the separation between the BH immediate progenitor and its companion at the moment of core collapse, can be expressed as

$$\begin{equation} r = A_{{\rm preSN}} (1-e_{{\rm preSN}}\cos \mathcal {E}_{{\rm preSN}}), \end{equation} \tag{ 17 }$$

where the pre-SN eccentric anomaly $\mathcal {E}_{{\rm preSN}}$ is related to the pre-SN mean anomaly $\mathcal {M}_{{\rm preSN}}$ as

$$\begin{equation} \mathcal {M}_{{\rm preSN}} = \mathcal {E}_{{\rm preSN}} - e_{{\rm preSN}} \sin \mathcal {E}_{{\rm preSN}}. \end{equation} \tag{ 18 }$$

The relative orbital speed V_BHpro of the BH immediate progenitor in the pre-SN binary can be written as

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

Finally, the angle ψ is the polar angle of the position vector of the BH immediate progenitor with respect to its relative orbital velocity in its companion's frame, which is related to the pre-SN orbital parameters as

$$\begin{equation} \sin \psi = \left[ \frac{A^2_{{\rm preSN}} (1-e^2_{{\rm preSN}})}{r (2A_{{\rm preSN}} - r)} \right] ^{1/2}. \end{equation} \tag{ 20 }$$

We start our calculation at the instant just before the core collapse event. The properties of the companion star are taken from a stellar model in a grid constructed by varying the initial mass between 65 and 105 M_☉, in steps of 0.1 M_☉. It is indeed the same grid of models used in finding the correct formalism describing CE evolution, which is capable of explaining the post-CE binary properties given by all successful evolutionary sequences (see Sections 5.1 and 5.5). To obtain the properties of these companion star models at the core collapse event of the BH progenitor, we need to approximate when the BH is formed in the evolutionary timeframe of its companion (t_BH). Under the assumption that the BH progenitor and its companion are born at the same time, t_BH simply equals the lifetime of the BH progenitor. We adopt that to be the lifetime of a stellar model with an initial mass of 150 M_☉, which is approximately 2.9 Myr. Other pre-SN binary properties—namely M_BHpro, A_preSN, e_preSN, and $\mathcal {M}_{{\rm preSN}}$—are drawn randomly from uniform distributions. Supplemented with a natal kick magnitude V_k and direction angles (θ_k, ϕ_k) drawn randomly from uniform and isotropic distributions, we can obtain A_postSN and e_postSN from the pre-SN binary properties using Equations (15) and (16). In this calculation, M_BH is taken directly from a successful evolutionary sequence, because the BH in our adopted formation scenario of IC 10 X-1 has only accreted a negligible amount of mass since its birth.

For each combination of companion star model and successful evolutionary sequence, we perform 2000 Monte Carlo trials. We only retain the data points that satisfy all of the following constraints and classify those data points as successful.

• 1.  
  M_BHpro is set to be less than 60 M_☉, which is a conservative upper limit guided by the study of Belczynski et al. (2010) on the maximum mass of stellar BHs. We will discuss the impact of this limit on our derived constraints related to the BH formation in Section 7.1.
• 2.  
  The binary must survive through the core collapse event. This means A_postSN and e_postSN obtained from Equations (15) and (16) needs to have realistic values: A_postSN > 0 and 0 ⩽ e_postSN < 1.
• 3.  
  Since the core collapse event of the BH immediate progenitor is instantaneous, the separation between the pre-SN binary components is the same as that of the post-SN binary. This gives a constraint as

  $$\begin{eqnarray} A_{{\rm preSN}}&& (1-e_{{\rm preSN}}\cos \mathcal {E}_{{\rm preSN}}) \nonumber \\ &&= A_{{\rm postSN}} (1-e_{{\rm postSN}}\cos \mathcal {E}_{{\rm postSN}}), \end{eqnarray} \tag{ 21 }$$

  which has to be satisfied with a realistic post-SN eccentric anomaly: $|\cos (\mathcal {E}_{{\rm postSN}}) | \le 1$.
• 4.  
  Both pre-SN binary components cannot spin faster than the breakup angular velocity $\Omega _c \approx \sqrt{G M/R^3}$. Here, we assume their spins are pseudo-synchronized to the pre-SN orbital angular velocity. Since the BH immediate progenitor is expected to be He rich due to potential binary interaction or intense mass loss via stellar wind, we approximate its radius using Equation (3) in Fryer & Kalogera (1997).
• 5.  
  According to our adopted formation scenario of IC 10 X-1, both components of the pre-SN binary need to fit within their Roche lobes at periapsis (Sepinsky et al. 2007).
• 6.  
  The BH companion in the post-SN binary must later fill its Roche lobe at periapsis before its radius reaches its maximum value (see Section 5), which leads to CE evolution. As massive stars evolve at roughly the same nuclear time scale, the time difference between the formation of the BH and the onset of the CE event is small. Hence, we assume that the binary semi-major axis and eccentricity remain unchanged within this period of time. The outcome of the CE phase is indeed constrained by the post-CE binary properties of the corresponding successful evolutionary sequence. When the BH companion fills its Roche lobe at periapsis, its core mass needs to match the corresponding mass of the post-CE He star, with a tolerance of 1 M_☉. This tolerance value is chosen according to the initial mass resolution in our grid of post-CE He star models, whose properties are used in constructing the successful evolutionary sequences (see Section 4). Furthermore, the CE efficiency α_CE is determined by the standard α prescription (Webbink 1984) with the "enthalpy" formalism (see Section 5.5) must be ⩽1.

The elements presented in the previous sections can now be combined to establish a complete picture of how we track the evolution of IC 10 X-1 backward in time, and derive constraints related to the BH formation in this system. We first use the modeling of binary evolution and observational constraints to determine the post-CE binary properties. Specifically, our successful evolutionary sequences at the present time simultaneously match the measured component masses, He star luminosity, and binary orbital period of IC 10 X-1. Then, we search for the correct formalism and treatment of CE evolution leading to the formation of IC 10 X-1. We find that the standard α prescription (Webbink 1984) with the "enthalpy" formalism (Ivanova & Chaichenets 2011) is capable of explaining the post-CE binary properties given by our successful evolutionary sequences. Last, we use our findings in the two previous steps as part of the constraints applied on a Monte Carlo simulation of the binary orbital dynamics involved in the core collapse event. Each data point in this simulation contains seven random parameters: the BH immediate progenitor mass, the pre-SN orbital mean anomaly, semi-major axis and eccentricity, the magnitude of the natal kick velocity imparted to the BH, and two angles describing the kick direction. These random parameters are drawn from uniform or isotropic prior distributions. If a data point satisfies all constraints mentioned in Section 6, we classify it as a successful data point. Our results and derived constraints (at 95.4% of confidence) presented in the following sections are all obtained from successful data points, as well as the marginalized PDFs illustrated in Figure 8.

Zoom In Zoom Out Reset image size

7.1. Core Collapse Constraints

Just before the core collapse event, we find the BH immediate progenitor mass (M_BHpro) to be 46 ± 14 M_☉ and its highly evolved main-sequence companion mass (M₂) to be 82 ± 7 M_☉. At this time, the companion star has ∼4–18% of H left in its core. The orbital separation (r) between the BH immediate progenitor and its companion is $8100_{-7100}^{+71000} \,R_\odot$. During the core collapse event, the BH immediate progenitor loses ⩽50% of its mass, which is ⩽20% of the total mass in the pre-SN binary. Possible asymmetries developed in the core collapse event can lead to a natal kick (V_k) of ⩽130 km s⁻¹ imparted to the BH. The PDFs of V_k and M_BHpro are presented in Figure 8.

We notice that the M_BHpro PDF contains a plateau between 38 and 60 M_☉, which indicates there is no clear upper limit on M_BHpro aside from the expected maximum value based on the stellar evolutionary theory. Fortunately, our derived upper limit on V_k depends very weakly on our adopted value of maximum M_BHpro. Figure 9 illustrate the two-dimensional joint V_k–M_BHpro confidence levels. For M_BHpro above 38 M_☉, the boundaries of the confidence levels are almost perpendicular to the V_k axis. This means the V_k PDF remains roughly the same even if we omit the data points in a certain mass range in this regime. Hence, our derived upper limit on V_k will not change significantly if we adopt a lower M_BHpro upper limit.

Zoom In Zoom Out Reset image size

7.2. Common Envelope Evolution Constraints

In this analysis, we track the evolution of IC 10 X-1 backward in time up to the instant just before the core-collapse event and study the formation of the BH in this system. This covers the following evolutionary phases: binary orbital dynamics at core collapse, CE evolution, and evolution of the BH–He star binary progenitor of the observed system. We first focus on the latter, using the modeling of binary evolution to construct successful evolutionary sequences to determine the post-CE binary Properties. These sequences are referred as successful because their properties at present match the following observational constraints of IC 10 X-1: binary orbital period, component masses, and luminosity of the W-R star. Our predicted X-ray luminosity at the present time, resulting from the stellar-wind accretion onto the BH, is consistent with the observed values. We then analyze the evolution through the necessary CE phase. We employ different CE treatments, as the standard treatment leads to unphysical results. We find that only the "enthalpy" formalism (Ivanova & Chaichenets 2011), along with an energy-based CE efficiency (Webbink 1984) can explain physically the post-CE binary properties of the IC 10 X-1 progenitor. Finally, we perform a Monte Carlo simulation on the orbital dynamics involved in the core collapse event. Each data point contains seven free parameters drawn from uniform and isotropic distributions, which describe the properties of the pre-SN binary and the natal kick imparted to the BH. Aside from the constraints related to the core collapse event, we also use what we learned about the CE event involved in the formation of IC 10 X-1 as additional constraints to reject the data points. If a data point satisfies all the constraints mentioned in Section 6, such as the survival of the binary through the core collapse event and the CE efficiency α_CE being ⩽1, we classify it as a successful data point. Our constraints (at 95.4% of confidence) related to the BH formation and the CE event occurred in the past of IC 10 X-1 are derived from all successful data points. We find that the BH immediate (He rich) progenitor has a mass of 46 ± 14 M_☉ and constrain the magnitude of the natal kick imparted to the BH to be ⩽130 km s⁻¹.

From the formation studies of low mass BH XRBs, envelope ejection of massive stars during CE evolution has long been known to be energetically problematic (e.g., Kalogera 1999; Podsiadlowski et al. 2003; Wiktorowicz et al. 2013). In this study, we adopt the energetic formalism of Webbink (1984) to calculate the CE efficiency (α_CE). Using the original Webbink (1984) prescription, we find that the binding energies (E_bind) of our He star progenitor models are at least three times larger than the available orbital energy, leading to α_CE > 1 (see Section 5.1). As E_bind of our He star progenitors can potentially be smaller by adopting a different core definition, we examine four other core definitions (see Section 5.2). Surprisingly we find that, for the very massive supergiants relevant here, these definitions all give roughly the same E_bind. This result is in contrast to what Tauris & Dewi (2001) found in their study of E_bind for red-giant branch (RGB) and asymptotic-giant branch (AGB), in which they argued that E_bind can be varied significantly by adopting a different core definition. However, the RGB and AGB stars considered by Tauris & Dewi (2001) are relatively much-lower mass stars, with initial masses up to 20 M_☉.

We also look at alternative CE treatments, but find that neither enhancing the convective core overshooting (Section 5.2) nor doubling the mass loss rates (Section 5.3) of our He star progenitor models can sufficiently decrease E_bind to make α_CE ⩽ 1. We also consider the formalism of hyper-critical accretion onto the BH during CE evolution (see Section 5.4). In order to have α_CE ⩽ 1, we find that the radius of the He star progenitor at the onset of the hyper-critical accretion phase needs to be at least 1.3 times larger than that at the onset of the CE phase, but there cannot be any pre-CE binary configuration satisfying this requirement. The last CE treatment considered in our analysis is the "enthalpy" formalism (Ivanova & Chaichenets 2011). We find that it naturally provides α_CE ⩽ 1 with realistic pre-CE binary configurations, because including the term P(m)/ρ(m) as an energy source lowers E_bind sufficiently. By adopting this CE treatment and considering all constraints on our derived evolutionary history of IC 10 X-1, we find α_CE to be in a range of 0.6–1 (at 95.4% of confidence).

Whether this "enthalpy" formalism is physically justified is not yet a settled issue (see Ivanova et al. 2013). There are other potential energy sources discussed in the literature that were not explored in our analysis. Ivanova (2002) and Voss & Tauris (2003) suggested that the energy released from accretion onto the compact object during CE evolution can contribute to envelope ejection. Soker (2004) argued that this type of accretion will produce jets that help disperse the envelope. A detailed discussion of all potential energy sources can be found in a recent review on CE evolution by Ivanova et al. (2013).

Alternatively, there are uncertainties in building our He star progenitor models that can lead to a decrease in E_bind. As mentioned in Section 3, our models did not account for the effects of rotation. Rotational mixing (Maeder & Meynet 2000; Maeder 2009) can increase the mass of the He core and change the internal structure of the H rich envelope. Also, we did not include rotationally enhanced mass loss rates, which can decrease the mass of the envelope. Furthermore, although we find that doubling the mass loss rates does not significantly lower the E_bind of our models (see Section 5.3), the uncertainties in the mass loss rates of massive stars beyond the main-sequence phase are larger than a factor of two. For instance, Petrovic et al. (2005) allowed a factor of six uncertainty in their mass loss rates for W-R stars. Given the uncertainties in modeling the evolution of massive stars, we cannot rule out the possibility that by fine-tuning our He star progenitor models we can decrease E_bind sufficiently and obtain α_CE ⩽ 1 without invoking the contribution of extra energy sources during CE evolution.

In this study, we uncover the evolution of IC 10 X-1 back to the point just prior to the formation of the BH. Although we did not extend our detailed binary modeling further backward in time, we can illustrate one scenario of how a primordial binary evolves to the current state of IC 10 X-1, which is based on the results of one successful data point in our analysis. We start with a binary consisting of two zero age main sequence stars, which are ∼ 150 and ∼87 M_☉. The binary orbit is initially circular, with an orbital radius of ∼8300 R_☉ (equivalent to a period of ∼16 yr). This orbit is so wide that throughout the evolution of the primary, its radius is at least four times smaller than its Roche lobe. Hence the orbital evolution is dominated by mass loss from the system via stellar wind. Soon after the more massive primary evolves off the main sequence, it loses its H rich envelope due to its massive stellar wind and becomes an He rich star. This He star also suffers from intensive mass loss. Just before collapsing to a BH, the mass of this He star is ∼44 M_☉. At the same time, the secondary approaches the end of its main sequence evolution, with a core H abundance of ∼9%. Its mass and radius are ∼82 M_☉ and ∼40 R_☉, respectively. Due to huge mass loss from the system, the orbital period increases to ∼55 yr. During the core-collapse event, the BH (He rich) immediate progenitor loses ∼11 M_☉ and forms a ∼33 M_☉ BH. Meanwhile, a natal kick of ∼70 km s⁻¹ is imparted to that BH. Right after the core-collapse event, the binary orbit is very eccentric (e ≈ 0.85) and the orbital period increases to ∼64 yr. In ≲ 0.3 Myr, the secondary becomes a supergiant and fills its Roche lobe at periastron, leading to CE evolution. At the end of CE evolution, the binary consists of a ∼33 M_☉ BH and a ∼35 M_☉ He star in a circular orbit, with an orbital period of ∼34 hr. This binary continues to evolve to the current state of IC 10 X-1. We emphasize that we do not necessarily consider this scenario as unique; instead we use it as the simplest possible scenario that can be consistent with the evolution of a primordial binary.

Based on our derived evolutionary history of IC 10 X-1, the spin angular momentum of the BH immediate progenitor is likely to be low. This is because once it loses its H rich envelope the BH progenitor will suffer from the high mass loss rates of W-R stars. This intense mass loss via stellar wind will take most of the angular momentum away from the BH progenitor and spin it down quickly. Tidal effects could have kept the BH progenitor from spinning down. However, the pre-SN binary orbit is relatively wide, with an orbital period of ⩾0.5 yr. This means that the tides exerted on the BH immediate progenitor by its companion star are expected to be weak. Even if the tidal interactions are much stronger than expected, they can at best synchronize the spin of the BH immediate progenitor with the orbital frequency. Hence, under the assumption that spin angular momentum is conserved during the core collapse event, we argue that the natal spin of the BH in IC 10 X-1 is likely to be small. Furthermore, since the BH has accreted a negligible amount of mass from the stellar wind of its companion, it cannot be significantly spun up after its formation. Therefore, we expect the current spin of the BH in IC 10 X-1 to be small as well.

The BH spin of IC 10 X-1 is expected to be measured in the coming years. If the spin of this BH turns out to be fairly high, it may imply that the BH was spun up by accreting a significant amount of mass during the CE evolution at a hyper-critical accretion rate. Another possible explanation could be that the BH was spun up during the core collapse event (see Blondin & Mezzacappa 2007).

This work has been primarily supported by the NSF Grant AST-0908930 awarded to V.K.; T.W. also acknowledges partial support through the CIC/Smithsonian pre-doctoral fellowship program and the hospitality of the CfA; T.F. acknowledges support from the CfA and the ITC prize fellowship programs. E.G. acknowledges support from the NWO under grant 639.041.129. The calculations for this work were performed on the Northwestern University Fugu cluster, which was partially funded by NSF MRI grant PHY-0619274, and on the Northwestern University Quest High Performance Computing cluster.

Let us denote the mass of the black hole (BH) by M_BH, the mass of the BH companion and its core mass by M_com and M_core, and the binary semi-major axis by A. Accretion onto the BH will be initiated when the binary semi-major axis equals to the BH companion's radius. As the time interval between the onset of the common envelope (CE) evolution and that of the accretion phase is relatively short, the masses of the BH and its companion are expected to be unchanged. Because of the very short circularization timescale due to the Roche lobe filling BH companion, the binary orbit is assumed to be circular at the onset of accretion. CE evolution and accretion onto the BH will end when the envelope of the BH companion is ejected to infinity.

Following Belczynski et al. (2002), the rate of change in M_BH and A with respect to M_com due to accretion onto the BH is given by:

$$\begin{eqnarray} [ c_d (M_{{\rm com}} &+& M_{{\rm BH}}) - M_{{\rm com}} ] \frac{dM_{{\rm BH}}}{dM_{{\rm com}}} \nonumber\\ &&= -\frac{M_{{\rm BH}} M_{{\rm com}}}{A}\frac{dA}{dM_{{\rm com}}} + M_{{\rm BH}}, \end{eqnarray} \tag{ A1 }$$

where c_d is the drag coefficient of the BH traveling in its companion's envelope. We adopt c_d to be 6 (Shima et al. 1985; Bethe & Brown 1998). In order to express these rates in ordinary parameters, we use the equation of the energy balance during CE evolution.

During the phase of accretion, the binding energy of the BH companion's envelope can be expressed as

$$\begin{eqnarray} E_{{\rm bind}} & = f_{{\rm ej}} \cdot \frac{G M_{{\rm com}} (M_{{\rm com}} - M_{{\rm core}})}{\lambda A} \nonumber \\ & = \frac{G M_{{\rm com}} (M_{{\rm com}} - M_{{\rm core}} - \Delta M_{{\rm BH}})}{\lambda A}, \end{eqnarray} \tag{ A2 }$$

where f_ej is the mass fraction of the ejected envelope as defined in Equation (13) and ΔM_BH is the amount of mass accreted onto the BH. Since the outcome of the CE event is known, as given by the successful evolutionary sequences, we express ΔM_BH using the known post-CE BH mass (M_{BH, f})

$$\begin{equation} \Delta M_{{\rm BH}} = M_{{\rm BH},f} - M_{{\rm BH}}, \end{equation} \tag{ A3 }$$

By substituting Equation (A3) into Equation (A2), we obtain

$$\begin{equation} E_{{\rm bind}} = \frac{G M_{{\rm com}} (M_{{\rm com}} - M_{{\rm core}} - M_{{\rm BH},f} + M_{{\rm BH}})}{\lambda A}. \end{equation} \tag{ A4 }$$

The parameter λ is a numerical factor scaling the binding energy of the BH companion's envelope during the phase of accretion, and is defined as

$$\begin{equation} \lambda \equiv \frac{G M_{{\rm com},i} (M_{{\rm com},i} - M_{{\rm core}})}{E_{{\rm bind},i} R_{{\rm com},i}}. \end{equation} \tag{ A5 }$$

Here, M_{com, i} and R_{com, i} are the pre-CE mass and radius of the BH companion, respectively. We notice that R_{com, i} is different from the BH companion's radius at the onset of the accretion phase (R_{com, acc}). The envelope binding energy E_{bind, i} of the pre-CE BH companion is defined in Equation (8).

The envelope of the BH companion is ejected at the expense of the binary orbital energy. According to Equation (12), we can write

$$\begin{equation} \alpha _{{\rm CE}} \cdot \frac{dE_{{\rm orb}}}{dM_{{\rm com}}} = \frac{dE_{{\rm bind}}}{dM_{{\rm com}}}. \end{equation} \tag{ A6 }$$

Here, the rate of change in orbital energy with respect to M_com is given by

$$\begin{equation} \frac{dE_{{\rm orb}}}{dM_{{\rm com}}} = -\frac{G}{2A} \left(M_{{\rm com}}\frac{dM_{{\rm BH}}}{dM_{{\rm com}}} - \frac{M_{{\rm BH}}M_{{\rm com}}}{A}\frac{dA}{dM_{{\rm com}}} + M_{{\rm BH}} \right). \end{equation} \tag{ A7 }$$

From Equation (A4), we derive the rate of change in E_bind with respect to M_com as

$$\begin{eqnarray} \frac{dE_{{\rm bind}}}{dM_{{\rm com}}} = & \frac{G}{\lambda A} \left [ M_{{\rm com}}\frac{dM_{{\rm BH}}}{dM_{{\rm com}}}\right. \nonumber\\ &\left. + \frac{ M_{{\rm com}}(M_{{\rm core}}-M_{{\rm com}} + M_{{\rm bh},f} - M_{{\rm BH}})}{A}\frac{dA}{dM_{{\rm com}}}\right. \nonumber \\ &\left. + 2M_{{\rm com}}-M_{{\rm core}} - M_{{\rm bh},f} + M_{{\rm bh}} \right ]. \end{eqnarray} \tag{ A8 }$$

Using Equations (A6)–(A8), we can write

$$\begin{eqnarray} \left(1 + \frac{2}{\alpha _{{\rm CE}} \lambda } \right) & M_{{\rm com}} \frac{dM_{{\rm BH}}}{dM_{{\rm com}}} \nonumber \\ = \frac{M_{{\rm com}}}{A} & \left[ M_{{\rm BH}} - \frac{2}{\alpha _{{\rm CE}} \lambda } (M_{{\rm core}}-M_{{\rm com}}+M_{{\rm BH},f}-M_{{\rm BH}}) \right] \nonumber\\ \times\, &\frac{dA}{dM_{{\rm com}}}- M_{{\rm BH}} \nonumber \\ -\, & \frac{2}{\alpha _{{\rm CE}} \lambda }(2M_{{\rm com}} - M_{{\rm core}} - M_{{\rm BH},f} + M_{{\rm BH}}). \end{eqnarray} \tag{ A9 }$$

Then, we can derive a system of two ordinary differential equations using Equations (A1) and (A9):

$$\begin{equation} \frac{dA}{dM_{{\rm com}}} = \frac{A}{M_{{\rm com}}} \left(1 + \frac{2 M_{{\rm com}}}{\alpha _{{\rm CE}} \lambda h_2} \right) \end{equation} \tag{ A10 }$$

$$\begin{equation} \frac{dM_{{\rm BH}}}{dM_{{\rm com}}} = -\frac{2 M_{{\rm BH}} M_{{\rm com}}}{\alpha _{{\rm CE}} \lambda h_1 h_2}, \end{equation} \tag{ A11 }$$

where

$$\begin{equation} h_1 = c_d (M_{{\rm BH}} + M_{{\rm com}}) - M_{{\rm com}} \end{equation} \tag{ A12 }$$

$$\begin{eqnarray} h_2 & =& \frac{M_{{\rm BH}} M_{{\rm com}}}{h_1} \left(1 + \frac{2}{\alpha _{{\rm CE}} \lambda } \right) + M_{{\rm BH}} \nonumber\\ &&+ \frac{2}{\alpha _{{\rm CE}} \lambda } (M_{{\rm com}} - M_{{\rm core}} - M_{{\rm BH},f} + M_{{\rm BH}}). \end{eqnarray} \tag{ A13 }$$

Using the known post-CE parameters as initial conditions, we can obtain the BH mass (M_{BH, i}) and semi-major axis (A_acc) at the onset of the accretion phase by numerically integrating Equations (A10) and (A11) from M_core to M_{com, i}. In order to calculate the pre-CE semi-major axis (A_i), we again use the energy balance equation α_CEΔE_orb = ΔE_bind. Here

$$\begin{equation} \Delta E_{{\rm orb}} = \frac{1}{2}G M_{{\rm BH},i} M_{{\rm com},i} \left(\frac{1}{A_{{\rm acc}}} - \frac{1}{A_i} \right) \end{equation} \tag{ A14 }$$

and

$$\begin{equation} \Delta E_{{\rm bind}} = \frac{G M_{{\rm com},i} (M_{{\rm com},i}-M_{{\rm core}})}{\lambda } \left(\frac{1}{R_{{\rm com},i}} - \frac{1}{R_{{\rm com,acc}}} \right). \end{equation} \tag{ A15 }$$

Since the radius of the BH companion equals the binary semi-major axis at the onset of the accretion phase, we set R_{com, acc} = A_acc and obtain

$$\begin{equation} A_i = \left[ 1 + \frac{2 (M_{{\rm com},i} - M_{{\rm core}})}{\alpha _{{\rm CE}} \lambda M_{{\rm BH},i}} \left(1 - \frac{A_{{\rm acc}}}{R_{{\rm com},i}}\right) \right]^{-1} A_{{\rm acc}}. \end{equation} \tag{ A16 }$$

As the binary orbit continues to shrink during CE evolution, A_i > A_acc. Also, the BH companion must completely fill its Roche lobe in the pre-CE binary orbit. Hence, A_i needs to satisfy the constraint

$$\begin{equation} R_{{\rm com},i} = A_i (1-e_i) r_L \end{equation} \tag{ A17 }$$

with a pre-CE orbital eccentricity (e_i) <1. Here, r_L is the approximated Roche lobe radius (Eggleton 1983).