Fast orbital shrinkage of black hole X-ray binaries driven by circumbinary disks

Recently, the black hole X-ray binary (BHXB) Nova Muscae 1991 has been reported to be experiencing an extremely rapid orbital decay. So far, three BHXBs have anomalously high orbital period derivatives, which can not be interpreted by the standard stellar evolution theory. In this work, we investigate whether the resonant interaction between the binary and a surrounding circumbinary (CB) disk could produce the observed orbital period derivatives. Analytical calculations indicate that the observed orbital period derivatives of XTE J1118+480 and A0620-00 can originate from the tidal torque between the binary and a CB disk with a mass of $10^{-9}~\rm M_{\odot}$, which is approximately in agreement with the dust disk mass detected in these two sources. However, Nova Muscae 1991 was probably surrounded by a heavy CB disk with a mass of $10^{-7}~\rm M_{\odot}$. Based on the CB disk model and the anomalous magnetic braking theory, we simulate the evolution of the three BHXBs with intermediate-mass donor stars by using the MESA code. Our simulated results are approximately consistent with the observed donor star masses, orbital periods, and orbital-period derivatives. However, the calculated effective temperatures of the donor stars are higher than indicated by the observed spectral types of two sources.


INTRODUCTION
Stellar mass black holes (BHs) are products of collapsing massive stars after they exhausted all nuclear fuel. Due to the ultra-strong gravitational field, anything including particles and electromagnetic radiation can not escape from inside of BHs. Therefore, the best objects detecting BHs are X-ray binaries where the dynamical masses of BHs can be estimated. At present, there exist two dozen BH candidates that have been identified in X-ray binaries (Remillard & McClintock 2006;Casares & Jonker 2014, for a review). Most of them (19 sources) have been defined as BH low-mass X-ray binaries (BHLMXBs) because their donor star masses are less than 1 M ⊙ . Study of BHLMXBs will be of importance in understanding astrophysical process associated with ultra-strong gravitational fields, stellar and binary evolution, and common envelope (CE) evolution (Li 2015, for a review).
In a standard CE model, it is difficult for low-mass donor stars to eject the massive envelope of BH progenitors during the CE phase (Portegies Zwart et al. 1997;Podsiadlowski et al. 2003). As a result, the population synthesis predicted a birth rate to be two orders of magnitude lower than that derived from observations (Li 2015). This difference can be solved by adopting an anomalously high CE efficiency parameter (α CE ) (Kalogera 1999;Yungelson & Lasota 2008;Kiel & Hurley 2006). As an alternative evolutionary channel, BHLMXBs may have evolved from BH intermediate-mass X-ray binaries driven by the anomalous magnetic braking of Ap/Bp stars (Justham et al. 2006) or surrounding circumbinary disks (Chen & Li 2006).
Recently, Wang et al. (2016) found that BHLMXBs can be formed if most BHs are produced through a failed supernovae mechanism, in which the BH mass is equal to that of the He or CO core mass of the progenitor.
In the standard theory forming BHLMXBs, the angular-momentum-loss mechanisms usually include three cases as follows: gravitational radiation, magnetic braking (Verbunt & Zwaan 1981), and mass loss (Rappaport et al. 1982).

BHLMXBS
The orbital-angular momentum of a BHLMXB is J = Ωa 2 M bh M d /(M bh + M d ), where a is the orbital separation, Ω the orbital angular velocity of the binary, M bh , and M d are the BH mass, and the donor star mass, respectively. Differentiating this equation, the change rate of the orbital period iṡ where β = −Ṁ bh /Ṁ d is the BH accreting efficiency, q = M d /M bh is the mass ratio of the binary. According to the first and the third term on the right hand side of Equation (1), the orbital-angular-momentum loss and the mass loss of the system can cause the orbit to shrink. However, the second term would produce a positive orbital-period derivative if material transferred from the less massive donor star to the more massive BH. In general, the angular-momentum-loss rate of BHLMXBs isJ =J gr +J mb +J ml +J ot , whereJ gr ,J mb ,J ml ,J ot represent the angular-momentum-loss rate caused by gravitational radiation, magnetic braking, mass loss, and other mechanisms, respectively. Table 1 lists the relevant observed parameters of the three BHLMXBs. The orbital-period-change rate originating from gravitational radiation iṡ where G is the gravitational constant, c the light velocity in vacuo. According to Equation (2), the orbital-period derivatives produced by gravitational radiation for 1118, 0620, 1991 are respectively ∼ 3.0, 2.0, 4.0 × 10 −13 s s −1 , which are obviously 2−3 orders of magnitude lower than the observed results. Based on the standard magnetic braking prescription given by Rappaport et al. (1983), the corresponding orbital-period derivative can be estimated to bė where R d is the donor-star radius. Adopting γ = 1, the orbital-period derivatives given by magnetic braking are ∼ 7.8, 3.8, 2.2 × 10 −12 s s −1 for 1118, 0620, 1991, respectively. These estimations are still one order of magnitude lower than these observed. Actually, 1118 should have a fully convective donor star, which is not generally thought to produce magnetic braking (Rappaport et al. 1983;Spruit & Ritter 1983) To account for the formation of compact BHLMXBs, Justham et al. (2006) proposed an anomalous magnetic braking (AMB) mechanism, which is caused by the coupling between the strong magnetic field of Ap/Bp stars and an irradiation-driven wind induced by the X-ray flux. The orbital-period derivative predicted by the AMB model is given bẏ where B s is the surface magnetic field of the donor star, f is the wind-driving efficiency. According to the equation given by King et al. (1996), the accretion rate of BHs can be estimated to beṀ bh ∼ 0.1, 0.5, 2.0 × 10 −9 M ⊙ yr −1 for 1118, 0620, and 1991, respectively (González Hernández et al. 2017). Assuming a winddriving efficiency of f = 0.001 and a surface magnetic field of B s = 5000 G, the resulting orbital-period derivatives areṖ ∼ 5.3, 6.1, 3.5×10 −11 ss −1 for 1118, 0620, and 1991, respectively. Even if taking such an ultra-strong field of 5000 G, the orbital-period derivative induced by AMB mechanism it still one order of magnitude lower than that of 1991. Therefore, it seems that there are other efficient angular-momentum-loss mechanisms to cause the rapid orbital decay of the three BHLMXBs. Dramatically, Muno & Mauerhan (2006) have detected that the excess mid-infrared-emission area are obviously larger than the binary-orbit areas these systems, and they suggested that it probably arise from a contribution of circumbinary (CB) disks. Recently, observations performed by the Wide-Field Infrared Survey Explorer have confirmed that these two sources should be surrounded by CB disks (Wang & Wang 2014). In this work, we attempt to explore whether a CB disk around these three sources could be responsible for their observed orbital period derivatives. In section 3, we describe the CB disk model, and constrain the CB disk masses. In Section 4, we use the MESA code to simulate the formation of the three BHLMXBs. Finally, we summarize the results with a brief conclusion and discussion in Section 5.

CB DISK MODEL
In this section, we investigate whether the rapid orbital-decay observed in the three BHLMXBs could be interpreted by CB disks around these sources. The resonant theory between a binary and its CB disk is based on a standard thin disk (Goldreich & Tremaine 1979;Artymowicz & Lubow 1994), in which H/R = 0.01−0.1 (H, and R are the thickness and the half angular momentum radius of the CB disk, respectively). This resonant torque can be estimated using the viscous torque of the CB disk, which can be written as the following relation (Lubow & Artymowicz 1996;Dermine et al. 2013) where M cb is the CB-disk mass, ν = R 2 (H/R) 2 αΩ d is the disk viscosity (α, and Ω d are the viscous parameter and the angular velocity of the CB disk, respectively). Therefore, the orbital separation derivative of BHLMXBs is given by (Lubow & Artymowicz 1996;Dermine et al. 2013) where l, and m are the time-harmonic number and the azimuthal number (Artymowicz & Lubow 1994), and µ is the reduced mass of the binary. Differentiating the Keplerian third law G(M bh + M d )/a 3 = 4π 2 /P 2 , we can obtain the orbital-period  derivative of BHLMXBs as follows, Assuming that the mass-loss rate of BHLMXBs during the mass transfer isṀ bh +Ṁ d = −1.0 × 10 −7 M ⊙ yr −1 (an ultra-high mass-loss rate), and M bh + M d = 10 M ⊙ , we can estimate the second term on the right hand side of Equation (7) to be −5 × 10 −9 yr −1 . For a binary with an orbital period of 0.5 d, the contribution of this term iṡ P ∼ −0.2 ms yr −1 , which is obviously lower than those of the three BHLMXBs. Therefore, in this section we ignore the effect of the mass loss on the orbital-period derivative. Combining equations (6) and (7), and considering the resonances are very weak (m = l) when the eccentricity e ≤ 0.1 √ α (Dermine et al. 2013), the orbital-period derivative predicted by the CB disk model iṡ The inner radius of the CB disk should locate a distance to the mass center of the BHLMXBs as r in = 1.7a, at which the disk would be tidally truncated (Taam et al. 2003;Dubus et al. 2004). In addition, the lack of excess flux at 24µm in the observation for 1118 and 0620 imply that the outer radius of the disk is near r out = 3a (Muno & Mauerhan 2006). Therefore, we can obtain a half angular momentum radius to be R = (r in + r out )/4 + √ r in r out /2 = 2.3a. Assuming that the CB disks in these three sources have the same relation between R and a, Equation (8)  ) and a binary parameter (1/µ). Based on the observed masses of two components, and taking H/R = 0.1, α = 0.1, we can constrain the CB disk mass for the three BHLMXBs.
In Figure 1, we compare the orbital-period derivatives predicted by the CB disk scenario with observations in theṖ − 1/µ diagram. Muno & Mauerhan (2006) provided an estimation (∼ 10 −9 M ⊙ ) for CB disk masses surrounding 1118 and 0620. Since the CB disk masses of different BHLMXBs should have a dispersion, a two orders of magnitude mass range is considered. In Figure  1, the solid, dashed, and dotted curves represent the pre-dictedṖ derived by equation (8)  fitted by the theoretical line of 10 −9 M ⊙ , which is the estimated CB disk masses of these two sources. For 1991, a relatively heavy CB disk (∼ 10 −7 M ⊙ ) would be expected in order to account for the observed orbital-period derivative.

Input physics
In this section, we use a MESAbinary update version (8118) in MESA module (Paxton et al. 2015) to simulate the formation of the three BHLMXBs. The evolutionary beginning is assumed to be a binary system containing an intermediate-mass donor star (with a mass of M d ) and a BH (with a mass of M bh ). For the donor-star compositions, we adopt a solar compositions (X = 0.70, Y = 0.28, and Z = 0.02). Meanwhile, the two components are thought to be circularized at all times.
Once the donor star overflows its Roche lobe by a longterm nuclear evolution, the material would be transferred from the donor star to the BH through the inner Lagrangian point at a rate ofṀ tr . The accretion rate of the BH is limited to the Eddington rate as followṡ where X is the hydrogen abundance in the accreting material, and is the energy conversion efficiency of the BH, where M bh,0 is the initial BH mass (see also Bardeen 1970;King & Kolb 1999). Therefore, the accretion rate of the BH isṀ bh = min[Ṁ Edd , −Ṁ tr ]. If the accretion process is super-Eddington, we assume that a constant fraction δ of the lost mass feeds into the CB disk surrounding the BHXB, i. e. the mass increasing rate of the CB disk iṡ Similar to Chen & Podsiadlowski (2016), we consider the wind loss from the donor star is driven by X-ray irradiation. The irradiation-driving wind loss rate is given byṀ where f ir is the irradiation efficiency (in this work, we take f ir = 10 −3 ). We calculate the X-ray luminosity by L X = ηṀ bh c 2 . Therefore, the mass loss rate of the donor star isṀ d =Ṁ tr +Ṁ w . Assuming that 1118 originated from the Galactic disk and the donor has solar metallicity, Fragos et al. (2009) found that this system includes a ∼ 6.0 − 10.0 M ⊙ BH and a ∼ 1.0 − 1.6 M ⊙ donor star. However, some clues indicate that an intermediate-mass ( 2.0 M ⊙ ) should be a plausible range for the progenitor mass of the donor stars in BHLMXBs. First, it still remains controversial whether a donor star with a mass less than 1.5 M ⊙ can provide sufficient orbital energy to eject the envelope of the black-hole progenitor (Podsiadlowski et al. 1995;Portegies Zwart et al. 1997;Kalogera 1999;Podsiadlowski et al. 2003). Second, CNO-processed elements were observed on the surface of 1118 (Haswell et al. 2002), which implies its progenitor should be an intermediate-mass star.
In the input physics calculating the evolution of binary stars, orbital angular-momentum losses are key issue. In the MESA code, we consider four types of orbital angular-momentum loss during the evolution of BHXBs: (1) gravitational-wave radiation; (2) anomalous magnetic braking: we adopt the same magnetic braking prescription given by Justham et al. (2006) and Chen & Podsiadlowski (2016); (3) mass loss: the mass loss from the vicinity of the BH is assumed to be ejected in the form of isotropic winds and to carry away the specific orbital angular of the BH, while the donor star winds carry away that of the donor star; (4) tidal torque produced by the interaction between the CB disk and the BHXB. According to Equation (5), the angular momentum loss rate extracting by the CB disk can be written asJ

Results
In our calculation, the donor stars in BHXBs are assumed to be Ap/Bp star with an initial mass of 3.0 M ⊙ and a surface magnetic field of 500 G, and the initial masses of BHs are 6.0 M ⊙ (for 1118 and 0620) and 10.0 M ⊙ (for 1991). To fit the CB disk mass inferred in Section 3, a faction δ = 5.0 × 10 −9 and 5.0 × 10 −7 of the mass loss during the super-Eddington accretion is thought to feed into the CB disk. By changing the initial orbital periods, we can diagnose whether the relevant BHXBs can evolve into the three observed sources by comparing the donor star masses, orbital periods, and orbital-period derivatives.
Our calculation show that the CB disk masses are approximately consistent with the inferred mass in Section 3 when the initial orbital-periods are 1.21 d and 1.71 d for 1118 and 1991, respectively. In Figure 2, we plot the evolution of BHXBs in the P orb − M d diagram. It is clear that 1118 and 0620 can evolved from a BHXB with an initial orbital-period P i = 1.21 days, while the progenitor of 1991 should have an initial orbital-period P i = 1.71 d. Figure 3 presents the evolution of orbitalperiod derivatives with the orbital periods. Both cases are approximately in agreement with the observed values of 1118 or 0620, and 1991. Figure 4 summarizes the evolution of the mass transfer rate and the CB disk mass. When P i = 1.21 days, the CB disk mass reaches a maximum of ∼ 1.4 × 10 −9 M ⊙ at the donor-star mass of ∼ 0.4 M ⊙ (at this moment P orb ≈ 0.8 d) due to the mass transfer of super-Eddington. Because the CB-disk mass increases very slowly, the orbitalperiod first increases, and then sharply deceases due to a relatively high disk mass (see also equation 13, results in an efficient angular-momentum loss) when P orb ≈ 2.0 d. Subsequently, two short reversals of the orbital period also correspond to the two new CB-disk masses. For the case P i = 1.71 days, the CB disk mass rapidly increases to a maximum of M cb ∼ 1.0 × 10 −7 M ⊙ when the orbital period is ∼ 1.1 d and the donor-star mass is 1.9−2.0 M ⊙ . At the current stage, 1118 and 1991 have an mass transfer rate of ∼ 2 × 10 −9 and ∼ 2 × 10 −7 M ⊙ yr −1 , respectively.
In figure 5, we also compare the simulated results with the effective temperatures indicated by the observed donor-star spectral types and orbital periods in the T eff − P orb diagram. The effective temperature of the donor star in 0620 is approximately consistent with the simulated result. However, the donor stars in 1118 and 1991 were detected as cool spectral types. A similar problem had already been noticed by the previous works performed by Justham et al. (2006) and Chen & Li (2006).

DISCUSSION AND SUMMARY
Recently, the three BHLMXBs including 1118, 0620, and 1991 were reported to be experiencing an extremely fast orbital decay. The detected orbital-period derivatives are 1 − 3 orders of magnitude higher than those given by gravitational radiation, and standard magnetic braking. For the AMB, the estimatedṖ for 1991 is still one order of magnitude lower than observations even if the donor star has an ultra-strong magnetic field of 5000 G.
In this work, we attempt to explore whether the observed orbital decay can be interpreted by the existence of CB disks around BHLMXBs. Adopting some typical CB disk parameters H/R = 0.1, and α = 0.1, the observedṖ in the three sources could be explained by a surrounding CB disk with a mass of 10 −9 M ⊙ (for 1118 and 0620) or 10 −7 M ⊙ (for 1991). Dramatically, the inferred CB disk masses are approximately consistent with the observed results in mid-infrared emission for 1118 and 0620 (Muno & Mauerhan 2006). The CB disks surrounding BHLMXBs may originate from three following channels: (1) CB disks are the remnants of the common envelope. In principle, compact binary systems should experience a common envelope evolutionary phase (Ivanova et al. 2013). If the common envelope can not be fully ejected, the remaining material may collapse into a CB disk surrounding the binary sys-tem . (2) CB disks are the products of the mass transfer. A fraction of the mass loss probably form a disk structure surrounding the binary rather than leave it (van den Heuvel & de Loore 1974;van den Heuvel 1994). (3) CB disks could be fed by mass loss during single outburst or successive outburst in BHLMXBs (Xu & Li 2018).
In this work, we also simulate the formation of the three BHLMXBs 1118, 0620, and 1991 by using the MESA code. In the calculation, we assume that a fraction δ of the mass loss during super-Eddington accretion of BHXBs forms a CB disk surrounding the binary. To fit the CB disk mass inferred in Section 3, δ should be 5 × 10 −9 and 5 × 10 −7 for 1118 (or 0620), and 1991, respectively. Our simulations indicate that, the progenitor of 1118 (or 0620) may be a BH intermediate-mass Xray binary consisting a 6.0 M ⊙ BH and a 3.0 M ⊙ donor star, and with an initial orbital period of 1.21 d; while the progenitor of 1991 should have a heavy BH (10.0 M ⊙ ), and a relatively wide orbit (initial orbital period is 1.71 d). Our simulated donor-star masses, the donor-star radii, the orbital periods, and the orbital-period derivatives are approximately in agreement with the observed results. For 1118, the calculated mass-transfer rate is 2 × 10 −9 M ⊙ yr −1 at the current stage. However, the observed peak luminosity of 1118 is ∼ 10 −3 L Edd (L Edd is the Eddington luminosity) (Wu et al. 2010), which implies the accretion rate of the BH is ∼ 10 −10 M ⊙ yr −1 . Narayan & Yi (1995) proposed that the critical ratė M crit ∼ α 2Ṁ Edd (α is the viscous parameter of the accretion disk) for the advection-dominated accretion flow. If we take α = 0.1, thenṀ crit ∼ 10 −9 M ⊙ yr −1 , which is the same order of magnitude with our simulated masstransfer rate. Therefore, the advected energy are probably lost into the BH, and the radiation efficiency of the accretion disk in 1118 is relatively low.
Recently, Xu & Li (2018) also employed the CB disk model to account for the fast orbital decay in these three sources. However, their initial donor star masses are 1.0 M ⊙ , which is difficult to result in CNO-processed elements detected on the surface of 1118 (Haswell et al. 2002). In addition, they assumed that the CB disk is formed due to single outburst or successive outburst at current time. Therefore, their model only produced a high orbital period derivative in a relatively short timescale.
Certainly, our simulation present a relatively high effective temperature of the donor stars. The main reasons could be as follows. First, Torres et al. (2004) found that the donor star of 1118 was only detected ∼ 55 percent light during quiescence, hence the determination for the spectral types of the donor stars in such systems are controversial. Second, the irradiation process of Xray could alter the effective surface boundary condition of the donor stars, especially change the ionization degree of the hydrogen at the bottom of the irradiate layer (Podsiadlowski 1991).
We thank the referee for his/her very careful reading and comments that have led to the improvement of the manuscript. This work was partly supported by the National Natural Science Foundation of China (under grant number 11573016, and 11733009), the Program for Inno- vative Research Team (in Science and Technology) at the University of Henan Province, and the China Scholarship Council. This work has also been supported by a Humboldt Research Award to PhP at the University of Bonn.