© ESO 2021

1. Introduction

In the last five years, the discoveries from ground-based gravitational wave (GW) detectors have driven a large effort to understand the origin of compact binary black holes (BHs). With 50 detections with a high significance (Abbott et al. 2019, 2020a), the majority of the observed GW sources correspond to binary BH mergers covering a large range of masses ∼5−100 M_⊙. Important discoveries contained in the first Gravitational Wave Transient Catalog (Abbott et al. 2019) include the discovery of BHs more massive than those detected through electromagnetic observations in the Galaxy (Abbott et al. 2016) and the direct association of short gamma-ray bursts with binary neutron star mergers (Abbott et al. 2017). The increased size of the sample on the second Gravitational Wave Transient Catalog (Abbott et al. 2020a) points to additional features in the properties of merging binary BHs, indicative of the contribution of multiple formation channels (Abbott et al. 2021; Zevin et al. 2021). Accurate predictions from different formation channels are then necessary to understand their relative contributions.

A large number of formation scenarios have been proposed to form the merging binary BHs observed by ground-based detectors. Scenarios that involve binary systems include evolution through a common-envelope (CE) phase (e.g., Paczynski 1976; van den Heuvel 1976; Tutukov & Yungelson 1993; Belczynski et al. 2002; Dominik et al. 2012; Stevenson et al. 2017; Giacobbo & Mapelli 2018), chemically homogeneous evolution (Mandel & de Mink 2016; Marchant et al. 2016; de Mink & Mandel 2016; du Buisson et al. 2020; Riley et al. 2021), stable mass transfer (MT; van den Heuvel et al. 2017; Neijssel et al. 2019; Bavera et al. 2021) and Population III stars (Belczynski et al. 2004; Kinugawa et al. 2014; Inayoshi et al. 2017). Dynamical processes are also predicted to contribute to the observed sample, including isolated triple systems (Thompson 2011; Antonini et al. 2017; Vigna-Gómez et al. 2021) and interactions in globular (e.g., Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Rodriguez et al. 2015; Di Carlo et al. 2019) and nuclear clusters (Antonini & Perets 2012). Additional formation scenarios include the pairing and growth of stellar mass BHs in the disks of active galactic nuclei (McKernan et al. 2014; Stone et al. 2017; Bartos et al. 2017), and primordial BH formation (Bird et al. 2016; Sasaki et al. 2018).

Our focus is on the formation of merging binary BHs through both CE evolution and stable MT. The majority of predicted binary BHs that formed through CE evolution or stable MT in an isolated binary involve an initial phase of stable MT. This removes the hydrogen envelope of the more massive star and leads to the formation of a wide single degenerate binary (Dominik et al. 2012; Langer et al. 2020). When the secondary evolves and thus fills its Roche lobe, depending on its response to mass loss and the mass ratio of the system, it can undergo a CE phase which ejects the envelope of the secondary at the cost of hardening the binary (cf. Belczynski et al. 2016). If the second phase of MT is stable, depending on the mass ratio of the system, a phase of nonconservative mass transfer can harden the binary without the need for CE evolution, as van den Heuvel et al. (2017) argue, thus potentially allowing for the formation of compact object binaries that can merge within a Hubble time. Recent studies have claimed that the contribution from stable mass transfer to the observed sample of merging binary BHs can be comparable or even larger than those that formed through CE evolution (Neijssel et al. 2019; Bavera et al. 2021).

However, binary evolution models have large uncertainties associated with the CE phase (see Ivanova et al. 2013 for a recent review). Important open issues include the unknown efficiency of CE evolution (Zorotovic et al. 2010; Davis et al. 2012; Toonen & Nelemans 2013), the mass boundary at which CE evolution terminates (the core-envelope boundary; Han et al. 1994; Dewi & Tauris 2000; Ivanova 2011), and the actual conditions for the onset of CE evolution (Hjellming & Webbink 1987; Soberman et al. 1997; Pavlovskii et al. 2017). Observational and theoretical constraints for the outcome of CE evolution in massive stars are also limited as the vast majority of observed post-CE systems as well as multidimensional hydrodynamical simulations of CE evolution correspond to low and intermediate-mass stars (see Iaconi & De Marco 2019 for a recent compilation of CE simulations and observed post-CE systems).

One uncertainty of CE evolution that is studied in population synthesis calculations is the impact of the evolutionary stage of the star at the onset of CE evolution. Belczynski et al. (2010) proposed that in low metallicity environments, a larger number of massive stars are expected to survive CE evolution, thus leading to an enhanced formation rate of merging binary BHs. This is argued as a consequence of halted expansion after core-hydrogen burning, with stellar models undergoing core-helium burning as blue supergiants and covering a larger range of radii after core-helium ignition. Belczynski et al. (2010) argued that interaction after core-helium ignition favors the ejection of the envelope as, at that point, a well defined core-envelope structure is formed. In various population synthesis codes, this argument is encoded in “pessimistic” and “optimistic” options for CE survival. In the pessimistic approach, all systems undergoing CE evolution before core-helium ignition are assumed to merge during CE evolution, while the optimistic approach allows for their survival. Recently, Klencki et al. (2020, 2021) have argued that this approach likely overestimates the number of merging binary BHs that formed as most of the systems predicted to survive CE in the optimistic scenario would initiate CE evolution with a radiative envelope that has a binding energy too large to eject through an inspiral.

The aim of this work is to study, in detail, how the evolutionary stage of the stellar progenitor of a BH affects MT stability and the outcomes of CE evolution. We do this by computing detailed models with the MESA code (Paxton et al. 2011, 2013, 2015, 2018, 2019) of a low metallicity massive star with a BH companion, meant to reproduce the evolution of a binary system after the formation of the first BH. Recent work on the formation of binary BHs using detailed binary evolution models (in contrast to rapid population synthesis calculations) has been presented by Eldridge & Stanway (2016) and Pavlovskii et al. (2017), covering a broad range of donor masses and metallicities. Our work, in contrast, focuses on a single donor mass of 30 M_⊙ at a metallicity of Z_⊙/10, but with a large resolution in the initial mass ratio and period. This allows us to study, in detail, how the stability of MT and the outcomes of CE evolution are affected by the evolutionary stage of the donor. For these calculations, we have made improvements to the commonly used MT prescription of Kolb & Ritter (1990, hereafter KR90) and implemented a method to model CE evolution which self-consistently determines the core-envelope boundary. Our methods are described in Sect. 2 and we present our results in Sect. 3. We conclude by discussing our results in Sect. 4. All input files necessary to reproduce our simulations as well as associated data products are available for download online¹.

2. Methods

We performed our simulations using version 15140 of the MESA stellar evolution code. Our models consist of binary systems with a 30 M_⊙ donor star and a BH companion. The objective of these simulations is to model a MT event happening after the formation of the first BH in the system, which is argued to be the most important MT phase where a CE can lead to the formation of a merging binary BH (cf. Dominik et al. 2012). Since we did not model the previous evolution of the system, we approximated the starting state of the donor star in the system as a zero-age main-sequence (ZAMS) star of its mass. Evolution was computed until either core carbon was depleted or until a merger between the donor star and the BH occurred.

Mass transfer rates were computed using the method described in Sect. 2.1, while the actual accretion rate into the BH was limited to its Eddington rate as described in Marchant et al. (2017). Our model for CE evolution is described in Sect. 2.2. Mass loss, either due to stellar winds from the donor or from mass ejected from the vicinity of the BH by MT above the Eddington limit, was assumed to take an amount of angular momentum corresponding to the specific orbital angular momentum of each component. We also accounted for the possibility of a star growing to the point that not only it overfilled its Roche lobe, but also overflowed its outer Lagrangian point. Whenever this happened, we accounted for the loss of mass and angular momentum as described in Sect. 2.1.5. Angular momentum loss from GW radiation was also included (Peters 1964), but it does not play a role in the timescales of evolution of our models which are of the order of millions of years. Our simulations do not take stellar rotation or spin-orbit coupling into account.

Our stellar models were computed at a low metallicity of Z = Z_⊙/10, where we took the relative metal mass fractions from Grevesse & Sauval (1998) and the solar abundance to be Z_⊙ = 0.0142 (Asplund et al. 2009). Opacities were computed using opacity tables from the OPAL project (Iglesias & Rogers 1996), together with the low temperature opacity tables of Ferguson et al. (2005). The equation of state used by MESA consists of a combination of OPAL (Rogers & Nayfonov 2002), HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), and SCVH (Saumon et al. 1995). Nuclear reaction rates were taken in order of preference from Cyburt et al. (2010) and Angulo et al. (1999).

Convection was modeled using the mixing length theory of Böhm-Vitense (1958) as described by Cox & Giuli (1968), using a mixing length parameter α_MLT = 2. Convective regions were determined using the Ledoux criterion (Ledoux 1947). We included overshooting of the hydrogen burning convective core using a step-overshooting scheme, extending the size of the convective core by α_ov = 0.335 pressure scale heights following the calibration of Brott et al. (2011). Overshooting from other convective regions is not well understood, so for convective cores after the main-sequence, we included a small amount of exponential overshooting (Herwig 2000) given by a length scale for exponential decay of the mixing coefficient of f = 0.01². We also included semi-convective mixing (Langer et al. 1983) and thermohaline mixing (Kippenhahn et al. 1980). Following Schootemeijer et al. (2019), we adopted a large value for the efficiency parameter for semiconvection α_sc = 100, while for thermohaline mixing we adopted an efficiency parameter of unity.

Stellar winds were accounted for using a combination of different mass loss rate prescriptions as described by Brott et al. (2011). This includes the line-driven mass loss rates derived by Vink et al. (2001) for stars with a hydrogen surface mass fraction of X > 0.7 and the Wolf–Rayet mass loss rate of Hamann et al. (1995) for X < 0.4, scaled by a factor of 10 to account for wind clumping (Yoon et al. 2010). For surface hydrogen mass fractions between X = 0.7 and X = 0.4, we interpolated between the rate of Vink et al. (2001) and Hamann et al. (1995) to provide a continuous transition. At temperatures below that of the bi-stability jump (as derived by Vink et al. 2001), we took the maximum between the rate of Nieuwenhuijzen & de Jager (1990) and the one resulting from the combination of the Vink et al. (2001) and Hamann et al. (1995) rates we just described. Additionally, we scaled the Nieuwenhuijzen & de Jager (1990) rate by the same factor (Z/Z_⊙)^0.85 predicted by Vink et al. (2001) for line-driven winds. This significantly lowers the mass loss rates of red-supergiants in low metallicity environments, although observations suggest there is only a weak dependence on metallicity (van Loon et al. 2005; Goldman et al. 2017). Nevertheless, we kept this scaling for the sake of comparison, as it is commonly used in population synthesis calculations.

2.1. Mass transfer

We modeled MT using an extension of the method developed by KR90, which accounts for MT from an extended atmosphere as described in Ritter (1988, hereafter R88) as well as the case where optically thick regions overflow the Roche lobe of the donor. In particular, this does not assume that the photosphere of the star operates as a hard rim, potentially leading to significant amounts of overflow during MT. Here we describe the method of KR90 in detail and our modifications to it, which include the possibility of outflows from the outer Lagrangian point of the donor. For the case of overflow from optically thick regions, we also discuss similarities between our method and that of Pavlovskii & Ivanova (2015).

2.1.1. The Roche potential

The computation of the MT rate depends on the Roche potential Φ, which is taken to be that of two point masses in a circular orbit:

(1)

where M_d is the mass of the donor and the mass ratio is defined as q ≡ M_a/M_d with M_a being the mass of the accretor. We used a coordinate system centered on the donor star (see Fig. 1), denoting distances normalized by the orbital separation a as ; and are the normalized distances of a point to the center of each star, and indicates the x-coordinate of the center of mass of the system. The x-coordinate of the first Lagrangian point is given by X_L1. The x-coordinate of the outer Lagrangian point of the donor (which can be either L₂ or L₃ depending on the mass ratio) is given by X_L out, which is defined as being negative. The Roche lobe radius of the system R_RL is defined as the volume equivalent radius of the region below the equipotential at the first Lagrangian point Φ_L1, and we used the fitting formula from Eggleton (1983) to approximate it:

(2)

Fig. 1. Definition of the coordinate system and variables used in the calculation of MT. The coordinate system is centered in the donor star, with the x-axis defined as the line joining the donor and the accretor, with the accretor located at x = a where a is the orbital separation. The z-axis is oriented along the direction of the orbital angular momentum of the system, while the y-axis completes a standard right-handed oriented Cartesian set of coordinates. The Lagrangian point located behind the donor is either L3 or L2, depending on whether the donor is more or less massive than the accretor, respectively.

Since we considered stars with significant overflow, potentially reaching the outer Lagrangian point of the donor, we also defined volume-equivalent radii beyond the equipotential of L₁ by dividing space with a plane crossing L₁, which is perpendicular to the line joining both stars, as depicted in Fig. 1. We refer to this plane as the L₁ plane. The volume equivalent radius for an arbitrary value of Φ is then defined as V(Φ) = 4πR(Φ)³/3, where V(Φ) is the volume contained within the equipotential Φ and on the side of the L₁ plane where the donor is located. In particular, we care about the radius associated with the equipotential of the outer Lagrangian point of the donor, for which we find that the following equation provides a fit with an error < 0.15% in the range of −10 < log₁₀q < 10³:

(3)

Figure 2 shows the value of R_L out/R_RL for different values of the mass ratio. In the limit of log₁₀q → ±∞, we have that R_L out/R_RL → 1, while in the range of −1 < log₁₀q < 1, which is typical of interacting binary stars, we find that R_L out/R_RL is between ∼1.2 and 1.3. Thus, a donor star growing beyond ∼20%−30% of its Roche lobe radius is expected to also fill its outer Lagrangian point.

Fig. 2. Ratio between the volume equivalent radius associated with the outer Lagrangian point of the donor and its Roche lobe radius. Bottom panel: result from numerical integration of the Roche potential together with the fit given in Eq. (3). Top panel: relative error between the fit and the data, with the small variations being due to the finite precision of the volume calculation.

2.1.2. Mass transfer through the L₁ point

Restricting the discussion first to flows through the L₁ point, the MT rate can be computed as an integral over the L₁ plane,

(4)

where ρ and v are the density and velocity of the fluid in the L₁ plane. The flow is assumed to be steady in which case the Bernoulli equation is satisfied,

(5)

where the integral is done along a streamline of the fluid with i and f denoting an initial and final point along a streamline. Following the work of Lubow & Shu (1975), the velocity near the Lagrangian point is expected to be equal to the sound speed. Depending on whether the photosphere of the star is inside the Roche lobe or has expanded beyond it, the flow is assumed to be isothermal (R88) or adiabatic (KR90), respectively. Combining either the adiabatic or isothermal approximation with Eq. (5) allows for the computation of ρ and v in the vicinity of the Lagrangian point as a function of the Roche potential Φ.

Either in the isothermal or adiabatic approximation, the surface integral in Eq. (4) can be expressed as an integral over the potential,

(6)

where A(Φ) is the area in the L₁ plane below the equipotential Φ. Figure 1 illustrates the area A(Φ_L out). The dA/dΦ term describes how the area enclosed by an equipotential in the L₁ plane changes with Φ and encodes the properties of the Roche geometry. Along the L₁ plane, given the symmetry of Φ and the L₁ point being a local minimum, the Roche potential can be approximated as

(7)

where ΔΦ ≡ Φ − Φ_L1. As shown in Appendix A, C₃ = C₅ = C₄/2 and the value of dA/dΦ can be computed up to first order in ΔΦ as

(8)

Both R88 and KR90 approximate dA/dΦ as a constant, but as we want to consider the case where the star can significantly overfill its Roche lobe, we included the term of order ΔΦ when computing MT in the optically thick case. R88 and KR90 also rely on fits to the coefficients C₁ and C₂ as a function of the mass ratio. To avoid relying on fitting formulae with a limited range of validity, we directly computed the location of the L₁ point, X_L1, and computed all C_i coefficients explicitly from the Roche potential⁴. The variation of dA/dΦ with increasing overflow has also been taken into account for optically thick Roche lobe overflow (RLOF) by Pavlovskii & Ivanova (2015), but using tables with numerical integrations of the area in the L₁ plane instead of a higher-order analytical approximation.

2.1.3. Optically thin mass transfer (R_ph < R_RL)

For the case of MT from a star with a photosphere radius R_ph below its Roche lobe radius R_RL, R88 considers the flow to have a constant temperature equal to T_eff, in which case the relevant sound speed is the isothermal one,

(9)

where k is the Boltzmann constant, μ is the mean molecular weight, and m_h is one atomic mass unit. This sound speed corresponds to the equation of an ideal gas plus radiation,

(10)

where a_rad is the radiation constant. We considered streamlines that start with negligible velocity from the photosphere of the star by taking v_i = 0 and Φ_i = Φ_ph in Eq. (5), where Φ_ph is the value of the Roche potential at the photosphere. Combining Eqs. (5), (9), and (10), one can solve for the density of the flow at the L₁ point,

(11)

where the ph and L1 subscripts indicate properties at the photosphere and the Lagrangian point.

Equation (11) can be evaluated in the vicinity of the L₁ point to express the variation in density across the L₁ plane as a function of the variation in Φ and v,

(12)

As for a steady flow, the L₁ point operates as a nozzle, and acceleration perpendicular to the flow is expected to vanish on the L₁ plane,

(13)

which combined with Eq. (12) results in dv = 0 (and thus v = v_th) in the vicinity of L₁.

These results can be used to evaluate the MT rate using Eq. (6). Transforming the integral into one over the density with Eq. (13) results in

(14)

As the density decays exponentially with the potential, see Eq. (11), only the regions in the very vicinity of L₁ are expected to contribute to the MT rate, so we approximated dA/dΦ by its value at L₁ to obtain

(15)

where M_0, thin is the MT rate for a star just filling its Roche lobe,

(16)

The factor that is still required to compute the MT rate in an evolutionary code is the difference Φ_ph − Φ_L1. R88 opted to use the derivative of the potential with respect to volume equivalent radii r_V covered by different equipotentials and approximated

(17)

Instead, we followed the method of Jackson et al. (2017) who developed an expansion of the Roche potential around the donor as a function of the equivalent radius r_V,

(18)

with a₁ and a₂ given by

(19)

Using this, we computed Φ_L1 with the Roche lobe radii from Eq. (2) and computed Φ_ph using the photospheric radius of our stellar model.

In all our simulations, only a small fraction of the total MT comes from the contribution of this thin MT rate. However, it still operates as a physically motivated mechanism to smoothly turn on MT which helps prevent numerical instabilities in the calculations when a star initiates RLOF.

2.1.4. Optically thick mass transfer (R_ph  >  R_RL)

Following KR90, if layers of the star below the photosphere are overflowing the Roche lobe, then the MT rate is computed as

(20)

where Ṁ_thick is the contribution to the integral in Eq. (4) from the overflowing, optically thick regions,

(21)

To calculate the density and velocity of the flow at L₁, we again made use of the Bernoulli equation, Eq. (5), by assuming the flow moves parallel to the equipotential surfaces and is adiabatic such that the pressure and sound speed are given along a streamline by

(22)

The constant k can be computed in terms of a reference pressure and density for each streamline, . This assumes that Γ₁ is constant along the streamline; Pavlovskii & Ivanova (2015) considered the case of a variable Γ₁ and found that the resulting MT rates are only modified by a few percent.

We assumed that along each equipotential surface, far from the Lagrangian point, the fluid is near hydrostatic equilibrium (dP₀ ≃ −ρdΦ and v ≪ v_s) and took the value of P₀ and ρ₀ there to compute k. Combining Eq. (22) with Eq. (5), and taking Φ_i = Φ_f, v_i = 0 and v_f = v_s, one can compute the density of the flow in the L₁ plane for a given value of ρ₀,

(23)

Combining this density with Eq. (21) gives the contribution to the MT rate from the overflowing optically thick layers as

(24)

where

(25)

In our calculations, we considered P₀ and ρ₀ to be the density in the overflowing layers of our hydrostatic model, from which we also obtained Γ₁. In this case, P(R_RL) and P(R_ph) are the pressure of the stellar model in the layer were R = R_RL and R = R_ph, respectively.

Equation (24) is exactly as derived by KR90, but in our calculations there are two important differences as to how we computed it. First, KR90 made the assumption of an ideal gas, replacing the P₀/ρ₀ ratio in Eq. (24) by kT/μm_h. In layers dominated by radiation pressure, this underestimates the mass loss rate significantly, so we did not make this assumption and computed P₀/ρ₀ from our stellar model. The second difference is that we did not assume dA/dΦ to be constant, but instead evaluated it using Eq. (8), with ΔΦ being computed in the overflowing layers of the model as

(26)

2.1.5. Mass loss from the outer Lagrangian point of the donor

Models of interactive massive stars using the prescription of KR90 can exhibit radii much larger than their Roche lobe radii. In particular, in the simulations presented in this work involving a 30 M_⊙ donor, we find cases where the outer Lagrangian point would overflow (R > R_L out). For these cases, we need a model for the overflow from this outer Lagrangian point, which is L₃ if M_d > M_a and L₂ otherwise.

The entire method described in Sects. 2.1.3 and 2.1.4 can be applied for an outflow through the outer Lagrangian point by replacing the Roche lobe radius with R_L out using Eq. (3) and evaluating the dA/dΦ factor at the outer Lagrangian point instead of L₁. However, an assumption needs to be made as to what happens to this outflow. We assume the material is ejected from the system carrying the specific angular momentum that corresponds to the outer Lagrangian point,

(27)

In order for this to be satisfied, the outflow needs to have a velocity much larger than the orbital velocity, otherwise tidal forces from the binary can modify the angular momentum content of ejected material. This has been studied in the case of overcontact binary systems that overflow the L₂ point (Shu et al. 1979; Pejcha 2014; Pejcha et al. 2016) or have their stellar winds torqued as they accelerate (MacLeod & Loeb 2020). If the energy of the outflowing material is not sufficient to unbind it from the binary, it can also accumulate in a circumbinary disk.

Owing to these uncertainties, we also experimented with models where the specific angular momentum in Eq. (27) is increased by factors of 4 and 9, corresponding to points corotating with the orbit at 2 and 3 times the distance of L_out from the center of mass. However, any of our simulations that undergo stable MT while overflowing L₂ or L₃ should be considered with care, since this is an extreme regime where three-dimensional hydrodynamic effects become important.

2.1.6. Outflows from L₂ when M_a/M_d < 1

An alternative approach to the problem of large overflow in simulations was taken by Misra et al. (2020). Considering systems with neutron star accretors, they computed the equivalent radii associated with the full L₂ equipotential (including both the donor and the accretor), and they assumed that whenever a stellar model exceeded this radius, then MT becomes unstable owing to the formation of an outflow with a high specific angular momentum. However, since for a compact object accretor there is not a hydrostatic structure within the Roche lobe of the accretor, the volume associated with the full L₂ equipotential cannot be directly compared to the radii of a hydrostatic stellar model in order to assess overflow of the outer Lagrangian points.

Nevertheless, in cases where M_d > M_a, some of the material streaming through L₁ could potentially be ejected from the vicinity of L₂ before the donor overflows L₃. In this case, our simulations would underestimate angular momentum loss from the system, since the specific angular momentum associated with L₂ can be up to 47% larger than that of L₃, as is shown in Fig. 3. Figure 3 also indicates the volume equivalent radii for the donor associated with the L₃ and L₂ equipotentials when M_a/M_d < 1, where for the donor’s L₂ equivalent radius, we used the fit of Marchant et al. (2016),

(28)

Fig. 3. Top: specific angular momentum in units of a2Ωorb for the donor, the accretor, and the first three Lagrangian points as a function of the mass ratio. Bottom: volume equivalent radii for the donor corresponding to the L2 and L3 equipotentials. Volume equivalent radius for L2 is taken from the fit of Marchant et al. (2016), which due to small errors in the fit results in a slightly larger value than the radius for the L3 equipotential at q = 1.

Although for a broad range of mass ratios, the donor can expand significantly before filling out the L₃ equipotential, much less expansion beyond the Roche lobe is required to fill out L₂.

In order to assess potential uncertainties during phases of large overflow, where material streaming through L₁ could be ejected from the vicinity of L₂ when M_a/M_d < 1, we also considered an alternate model. In this model, instead of accounting for overflows from the donor’s outer Lagrangian point, we assumed part of the mass streaming through L₁ is ejected from L₂. To do this, we computed an additional optically thin rate from Sect. 2.1.3 by replacing R_RL with R_L2 and evaluating dA/dΦ at L₂ rather than L₁. Whenever R_ph > R_L2, we considered the contribution to the optically thick mass transfer rate transferred through the L₁ plane, but above the equipotential of L₂, which is given by Eq. (24) with a modified lower interval for the integral:

(29)

Since this expression describes a flow through L₁, dA/dΦ has to be evaluated at L₁. We assumed both the optically thin component through L₂ and the optically thick component going through L₁ above the L₂ equipotential are ejected from the system, taking the specific angular momentum that corresponds to the L₂ point,

(30)

Similar to what was discussed in the previous section, material could be torqued or form a circumbinary disk before being ejected, so we also consider extreme cases where we increased the specific angular momentum removed by factors of 4 and 9.

The default setup of our simulations is the one described in Sect. 2.1.5, but we also performed simulations with this modified model with L₂ outflows to test the robustness of our results. For simplicity, in the remainder of the manuscript we do not use the subscript 0 to refer to properties from our hydrostatic stellar models and refer to the photospheric radius of a star as R.

2.2. Common envelope evolution

Whenever our MT prescription gives a MT rate exceeding a given threshold Ṁ_high, we assume evolution proceeds through a CE phase. We treated CE evolution following the standard energy prescription (Tutukov & Yungelson 1979; Webbink 1984) which equates the binding energy of ejected layers to the difference in orbital energy product of an inspiral,

(31)

where α_CE is a free parameter that represents the efficiency with which the orbital energy ejects the envelope. Using subscripts i and f to represent the pre- and post-CE properties of the system, the difference in orbital energy is

(32)

and we considered M_d, f to be the core mass of the donor M_core, while ignoring accretion into the accretor such that M_a, f = M_a, i. The binding energy depends on the value of M_core and was computed by adding up the internal and gravitational potential energy of the removed layers at the onset of CE,

(33)

where u is the specific internal energy of the gas and we included an additional free parameter, α_th, which was introduced by Han et al. (1995) to represent the efficiency with which thermal energy can be used to eject the envelope. In this work, we include the contribution of recombination energy from hydrogen and helium on u. For simplicity, we also assume α_th = 1 throughout.

Computing the binding energy requires one to know M_core, which is a nontrivial problem that can lead to large variations in the predicted binding energy (see for instance the discussion in Ivanova et al. 2013). Our goal is to determine M_core self-consistently in our simulations by modeling the stripping of the stellar envelope. For this, at the onset of CE, we computed the value of E_bind for all choices of M_core in the stellar model. The CE phase was then modeled by artificially removing mass from the star and computing, at each step, the binding energy one would have obtained from the pre-CE model if M_core were assumed to be the current mass of the star. In this way, Eq. (32) can be used to determine the final orbital separation a_f as a function of M_core, and the mass losing stellar model can be used to determine the point at which the star would contract inside its Roche lobe.

We did not simply force a rapid mass loss rate on our simulations until we found R < R_RL, but instead softly turned off the mass loss rate as the star went inside its Roche lobe to determine the mass coordinate at which the star would be undergoing contraction in the absence of mass loss. In particular, we forced a mass loss rate from CE evolution of

(34)

where we interpolated between a high MT rate Ṁ_high, meant to reproduce a near adiabatic response, and a small MT rate, Ṁ_low, which should be well below the thermal timescale of the envelope. In particular in our simulations, we chose Ṁ_high = 1 M_⊙ yr⁻¹, Ṁ_low = 10⁻⁵ M_⊙ yr⁻¹, and δ = 0.02. Our choice of Ṁ_high was based on the thermal timescale of the envelope estimated as τ_th = E_bind/L and the associated MT rate Ṁ_th = M/τ_th. For the 30 M_⊙ donor considered in this work evolved as a single star, and assuming a core-envelope boundary at the innermost point where X > 0.1 (Dewi & Tauris 2000), we find that Ṁ_th does not exceed 0.1 M_⊙ yr⁻¹. Similarly, our choice of Ṁ_low is such that it is comparable to MT happening on a nuclear timescale. The dependence of our results on the free parameters Ṁ_high, Ṁ_low, and δ is discussed in Appendix B, where we show how the outcome of some of our CE models are affected by these choices.

If the star contracts below R/R_RL = 0.98, then we assumed the CE phase finished and proceeded with the evolution of the model in the standard way, applying our MT formalism described in the previous subsections. As pointed out by Ivanova (2011), a star can significantly contract during a phase of rapid mass loss, but re-expand by orders of magnitude on a few thermal timescales if mass loss is suddenly shut off. Because of this, they determined the core mass of a stellar model by performing rapid mass loss simulations and finding at which mass coordinate the resulting star would not re-expand after mass loss is shut off. In our simulations, we did find re-expansion after CE, which, according to Ivanova (2011), would indicate that the core boundary is at a deeper mass coordinate. However, but we find that the ensuing MT phase is stable when using our MT prescription.

This method is similar to the one presented by Deloye & Taam (2010) and Ge et al. (2010), who computed the adiabatic response of the stellar radii to rapid mass loss and compared it to the evolution of the Roche lobe radii in a binary to determine at which point a binary would detach. The main difference is that our method determines the mass coordinate at which detachment would happen without mass loss, capturing the nonadiabatic response of the star near the end of CE. Another similar method is the one presented by Eldridge et al. (2017), which was integrated into detailed binary evolution calculations but does not consider the contribution of thermal or recombination energy.

3. Results

In this work, we focus on the evolution of a 30 M_⊙ star at a metallicity of Z_⊙/10 in order to assess how some of the differences expected from evolution in low metallicity environments can affect the outcome of either MT or CE evolution. As we did not perform a broad analysis for different donor masses and metallicities, we cannot perform population synthesis to compute, for example, the resulting rates of binary BHs observable by ground-based detectors. However, this serves to assess if some of the assumptions commonly used in rapid population synthesis codes are justified. In particular, Belczynski et al. (2010) argued that for BH progenitors, ejection of the envelope during CE evolution is more likely to happen in low metallicity environments, as these systems are predicted to halt their expansion in the Hertzsprung gap (HG) before rising on the Hayashi line as red supergiants, spending a significant fraction of their core helium burning lifetime as blue supergiants. After core helium ignition, it is argued that a well defined core-envelope structure is formed, which allows for a successful envelope ejection. A halted expansion after the main sequence is a common feature in stellar evolution calculations at low metallicity (see Appendix C).

Before discussing grids of binary simulations, it is instructive to study how our 30 M_⊙ donor evolves as a single star. The top panel of Fig. 4 shows the evolution in the Hertzsprung-Russell (HR) diagram for this star under the physical assumptions described in Sect. 2, for which we find that expansion is halted before becoming a red supergiant. Most of the helium-burning phase of this star would be spent as a blue supergiant with a convective envelope developing only near core helium depletion. This type of evolution is sensitive to physical assumptions in the model. In particular, Fig. 4 shows the effect of increasing overshooting during core hydrogen burning, as well as increasing wind strengths. With increased overshooting, the star undergoes rapid expansion after the main sequence and becomes a red supergiant before core helium ignition. The model with a higher mass loss rate slows down its expansion before becoming a red supergiant, but it expands to a larger radius during the HG than our default model. These variations are not captured in rapid population synthesis codes, in particular those reliant on the fits to single star evolution of Hurley et al. (2000) which are based on detailed stellar models without mass loss and with a fixed amount of overshooting. Stellar winds may be overestimated in our calculations (Björklund et al. 2021), while observations of massive stars suggest that our choice for overshooting from a convective hydrogen-burning core is underestimated for stars with masses > 15 M_⊙ (Castro et al. 2014). Keeping in mind that due to uncertainties in stellar evolution, a halted expansion in the HG is uncertain, we can still study what would be the consequences in the evolution of a star if this effect is real.

Fig. 4. Evolution of a 30 M⊙ single star with a metallicity of Z⊙/10, separated into different evolution phases. Dots indicate equal time intervals of 105 years. The Hertzsprung gap is defined as the evolutionary stage between the terminal-age main sequence and the moment where the ratio of the nuclear burning luminosity and the surface luminosity satisfies |log10(Lnuc/L)| < 0.001 for at least a thermal timescale estimated as GM2/RL, or the star develops a convective envelope > 1 M⊙. Top: evolution under our default set of physical inputs. Middle: evolution with core overshooting increased from 0.335 to 0.535 pressure scale heights. Bottom: evolution with wind mass loss rates tripled.

Focusing only on our model that uses our default set of assumptions, we assessed if core helium ignition has a significant impact on the core-envelope structure of the star which could favor CE ejection. Figure 5 shows a Kippenhahn diagram for the 30 M_⊙ model, indicating how the stellar radius and binding energy of the envelope evolve. Assuming that the core-envelope boundary is located at the points where the hydrogen mass fraction drops below either X = 0.1 or X = 0.3 (similar to Dewi & Tauris 2000), we find that the binding energy of the star remains at ∼10⁵⁰ erg from terminal-age main sequence (TAMS) to the point where a convective envelope is developed. This happens despite a ∼30 fold increase in the stellar radius. This model experiences the formation and merger of multiple intermediate convection zones during the HG phase, so we performed convergence tests to verify that our results are not significantly affected by changes in temporal and spatial resolution during this phase (see Appendix D).

Fig. 5. Kippenhahn diagram for our default 30 M⊙ model, separating different evolutionary stages as well as showing the evolution of the stellar radii. “He dep.” and “C dep.” stand for core helium and carbon depletion respectively. In addition, we show the result of computing Ebind under the assumption that the core–envelope boundary Mcore at which CE would finish is at the innermost mass coordinates where the hydrogen mass fraction is X = 0.1 or X = 0.3. We only show Ebind after the formation of the helium core, so no information is given on the main sequence. To provide a comparison point to rapid population synthesis codes, we also show the resulting Ebind from the fits of Claeys et al. (2014) and Xu & Li (2010).

If one computes the binding energy using, for example, the fit provided by Claeys et al. (2014) rather than the structure of the model, the resulting binding energy is always underestimated, with more that an order of magnitude difference in late stages⁵. In particular, this fit to the binding energy is used in population synthesis calculations of GW sources, including the MOBSE (Giacobbo & Mapelli 2018) and COSMIC (Breivik et al. 2020) codes. Another parametrization of the binding energy that is used in published rapid population synthesis results is that of Xu & Li (2010), which is an option in the STARTRACK (Dominik et al. 2012) and COMPAS (Neijssel et al. 2019) codes⁶. In this case, we also find a consistent underestimation of the binding energy except for late stages near core-carbon burning where the envelope becomes almost fully convective. Both Claeys et al. (2014) and Xu & Li (2010) include internal energy in their fits, but not recombination energy. A systematic underestimation of the binding energy in population synthesis calculations implies an overestimation of post-CE orbital separations, or a prediction of successful envelope ejection in systems that would merge during CE. We also compared the binding energies of our models with different overshooting and wind strengths to the results obtained from the Claeys et al. (2014) and Xu & Li (2010) fits, and find that the fits consistently underestimate the binding energy of the envelope as well (see Fig. 6). Binding energies computed using the Xu & Li (2010) fit experience a jump when the star transitions out of the HG, which for the track with extra overshooting does not happen until the star begins to rise in the Hayashi line.

Fig. 6. Binding energy of the envelope of a 30 M⊙ post main-sequence star as a function of radius. Tracks are shown for our default set of physical assumptions, as well as for a model with increased overshooting and one with increased mass loss rates. For each model, we also show the binding energy resulting from the fits of Xu & Li (2010) and Claeys et al. (2014).

The constancy of the binding energy while the star expands with a radiative envelope is also illustrated by Fig. 7, which shows the full binding energy profiles at different evolutionary phases. In particular, between the end of the HG and the point where a deep convective envelope starts to develop, the radius of the star increases by more than an order of magnitude, but this expansion only happens in the outermost mass layers of the star. It is only after the formation of the convective envelope that a large fraction of the envelope mass is pushed to large radii with correspondingly low binding energies that would make envelope ejection during CE feasible. All of this points to the ignition of helium at the stellar core having no impact on the binding energy of the envelope, and also since the ignition of helium does not modify the mode of energy transport in the outer layers, the stability of MT should not be impacted in a fundamental way either.

Fig. 7. Radius (top) and binding energy (bottom) profiles for a 30 M⊙ single star at different evolutionary stages. A profile near the end of core helium-burning (Yc = 0.08), but before the formation of a convective envelope, is included. Also, a profile is shown corresponding to the beginning of the formation of a deep convective envelope, taken as the first point in the evolution where the mass of the convective envelope Mconv env is larger than 1 M⊙.

The near constancy and high value of the binding energy before the formation of a convective envelope has also been studied by Kruckow et al. (2016) and Klencki et al. (2021) comprehensively, covering a large range of masses and metallicities. However, the binding energies computed by Kruckow et al. (2016) and Klencki et al. (2021) as well as those shown in Fig. 5 rely on an assumed boundary for the core mass. By performing binary simulations using the method described in Sect. 2.2 to determine the core-envelope boundary self-consistently, we can verify if these computed binding energies with an assumed boundary are a good approximation.

3.1. Grid of simulations

Using our prescription for MT and our method to calculate the outcome of CE evolution, we computed a grid of models consisting of a 30 M_⊙ donor with a BH accretor in a circular orbit. Our method for CE evolution still leaves the efficiency α_CE as a free parameter, and we computed simulations with both α_CE = 1, 0.4, 0.2, and 0.1. The grid spans a period range between −0.3 < log₁₀(P_i/d) < 3.5 and initial mass ratios q = M_BH/M_1, i between 0.01 and 0.59. This includes BHs with masses as low as 0.3 M_⊙ that are not expected to form through stellar evolution, but they were included to test the limits of MT stability. This allowed us to cover the entire parameter space where CE happens, as well as to cover systems from the boundary of RLOF at ZAMS, to effectively single star evolution at large periods. Owing to the small amount of radial expansion occurring in the phases where the star has a convective envelope, we adopted two different resolutions in the period for our grid. For log₁₀(P_i/d) < 3.2, we used a spacing of Δlog₁₀(P_i/d) = 0.1, while at higher initial orbital periods we increased the resolution in the period to Δlog₁₀(P_i/d) = 0.01. The resolution in mass ratio of the grid is Δq = 0.02.

The final outcomes of our grid with α_CE = 1 are shown in Fig. 8. As a simple approximation to determine if a simulation can produce a merging binary BH, for stars that reach core carbon depletion we assumed the donor star collapses directly to form a BH with a mass equal to its baryonic mass. Our 30 M_⊙ donor produces final helium core masses ∼14 M_⊙, for which the formation of a BH through direct collapse is indeed expected (Fryer 1999). From our simulations, we did not find any cases of successful envelope ejection from CE evolution for initial orbital periods below log₁₀(P_i/d) < 3.2. Below this period range, the boundary between systems undergoing stable MT and those merging during CE evolution varies continuously with the period, going from q ∼ 0.5 for systems interacting close to ZAMS, down to q ∼ 0.15 for cases where interaction happens right before the formation of a deep convective envelope. The boundary is not sensitive to whether an interaction happens before or after the TAMS, or before or after the end of the HG. This is to be expected as the radial response to rapid mass loss is dependent on the mode of energy transport in the stellar envelope (Hjellming & Webbink 1987), and the envelope of the star across all these phases is radiative. Pavlovskii et al. (2017) also find stable MT between a star and a BH for such extreme mass ratios for different donor masses, but owing to their coarse mass ratio sampling (no more than three different mass ratios for each donor mass), they cannot describe in detail how the boundary for stability shifts with the initial orbital period.

Fig. 8. Summary of final outcomes for simulations of circular binaries consisting of a BH and a 30 M⊙ star at a metallicity of Z⊙/10, and a CE efficiency parameter αCE = 1. Horizontal dotted lines indicate boundaries for an interaction before different evolutionary stages of the star. For systems that undergo stable MT, black horizontal lines indicate regions where the donor exceeded the L2 equipotential, while cross hatched regions mark regions where also the L3 equipotential is exceeded. Black rectangles indicate systems that do not interact or would result in a Roche lobe filling system at ZAMS. Systems that undergo stable MT or eject their envelope during CE evolution are denoted by “st. MT” and “CE ej.”, respectively, and they are separated into systems forming binary BHs that would merge in less or more than a Hubble time. Systems marked as “CE merger” merge during the CE phase.

The lack of systems surviving CE evolution between the end of the HG and the formation of a deep convective envelope can be understood in terms of the binding energy of the envelope. As was discussed for the single 30 M_⊙ model, the binding energy during this phase is ∼10⁵⁰ erg. Considering evolution between the end of the HG and the formation of a convective envelope, we find that CE evolution only happens for BH masses ≲5 M_⊙, and that stripping the hydrogen envelope of the 30 M_⊙ star would result in a ∼14 M_⊙ helium core. For these masses and assuming a fully efficient use of the orbital energy to eject the envelope, in order to release 10⁵⁰ erg, the orbital separation needs to be reduced below 1.3 R_⊙, an orbit for which a helium star of 14 M_⊙ would overfill its Roche lobe. Population synthesis calculations that underestimate the binding energy during this phase and find successful CE ejection, for example, by using fits such as those of Claeys et al. (2014) or by assuming a fixed value of the λ parameter, potentially overestimate the rate of formation of merging binary BHs through CE evolution (Gallegos-Garcia et al., in prep.).

The boundary for unstable MT is also right at the place where it allows for a narrow band of systems undergoing stable MT that produce a binary BH merging in less than a Hubble time. In this case, the shrinkage of the orbit from stable MT is sufficient to produce a compact BH binary (van den Heuvel et al. 2017), and we find this outcome for a wide range of initial orbital periods ranging from 1 to 1000 days. Multiple systems undergoing stable MT fill their outer Lagrangian point during that process, but remain stable despite the increased angular momentum loss. For wider orbital periods, where MT occurs after the development of a deep convective envelope, the critical mass ratio for stability shifts upwards (up to q ∼ 0.5 for systems interacting after core helium depletion). The majority of our models with an assumed efficiency of α_CE = 1 that manage to eject their envelopes during CE evolution result in wide BH binaries, for which GW radiation is insufficient to produce a merger in less than a Hubble time. Figure 9 shows the outcome of our simulations with CE efficiencies of α_CE = 0.4, 0.2, and 0.1, in which case we find that most systems that produce wide binary BHs at α_CE = 1 would merge during CE evolution rather than produce merging binary BHs at lower efficiencies.

Fig. 9. Same as Fig. 8, but for CE efficiency parameters of αCE = 0.4, 0.2, and 0.1.

One important concern in our simulations that form merging binary BHs through stable MT is that an important fraction of these undergo overflow of L₂ or even L₃. The methods we have developed to account for this are necessarily an approximation, as they rely on hydrostatic one-dimensional stellar models. Moreover, the simulations shown in Fig. 8 do not account for potential outflows through L₂ when M_a < M_d, which could remove significant mass and angular momentum from the system. Using the alternate method described in Sect. 2.1.6 for L₂ outflows, we can test if the systems that form merging binary BHs through stable MT are sensitive to the potential presence of such outflows. Figure 10 shows simulations done using this method that cover the region where we find that stable MT leads to merging binary BHs, including cases where we assume that the specific angular momentum removed by material outflowing L₂ is up to nine times the specific angular momentum of the L₂ point. Compared to our simulations that only account for overflow of the outer Lagrangian point of the donor, we find that the boundary for instability is pushed to higher BH masses. However, we still find systems across a large range of periods that form merging binary BHs, as additional angular momentum losses lead to the formation of more compact binaries post-MT. Given this, the conclusion that stable MT can lead to the formation of merging binary BHs across a large range of periods appears to be robust.

Fig. 10. Same as Fig. 8, but for a subset of models computed using the method described in Sect. 2.1.6 that considers outflows from L2 rather than L3 when Ma < Md. Each panel indicates simulations performed with different specific angular momentum assumed for the material outflowing through L2, which is considered to be jL2, 4jL2, or 9jL2. The purple dotted line indicates the boundary between stable MT and CE evolution as determined by our models that account only for outflows from the outer Lagrangian point of the donor.

As our models do not include the evolution prior to the formation of the first BH, we cannot assess the relative likelihood of each of our individual simulations to occur in nature. However, if we make the simplistic assumption that after BH formation the orbital period is distributed flat in log₁₀P and the mass ratio distribution is flat in q, our grid with α_CE = 1 gives a ratio between the number of merging binary BHs that formed through CE evolution to those that formed via stable MT of 0.017. The actual outcome of our simulations for systems surviving CE evolution should be considered with care, since we approximated a complex three-dimensional phenomenon with one-dimensional hydrostatic models. However, even under the extreme case that all our systems that survive CE evolution would produce a binary BH that merges in less than a Hubble time, we find that the ratio between merging binary BHs that formed through CE evolution to those that formed via stable MT is 0.35. This resembles the results from population synthesis calculations of Neijssel et al. (2019), who find that depending on the metallicity distribution of the star formation history, most merging binary BHs that will be observed by current generation GW detectors at design sensitivity would have been formed through stable MT. Our results suggest that the yield of merging binary BHs from CE evolution is overestimated at low metallicities, which would further increase the relative contribution of stable MT relative to CE evolution in population synthesis calculations such as those of Neijssel et al. (2019).

3.2. Mass transfer and outer Lagrangian point overflow

We consider here two cases from our simulation grid that lead to the formation of a merging binary BH: a system interacting before the end of HG (M_BH = 7.5 M_⊙, P_i = 31.6 days) and one interacting after it (M_BH = 4.5 M_⊙, P_i = 1000 days). The short period system evolves without overflowing its outer Lagrangian point, while the long period one is expected to overflow it, so we can use it to study how the resulting outcome is modified when accounting for this on our MT prescription. To evaluate the impact of the changes we have made to the MT prescription of (KR90), we used four different setups for our simulations:

1. A full implementation of our MT prescription including overflows from the outer Lagrangian point as described in Sect. 2.1.5.

2. Same as 1, but without the assumption P/ρ = kT/μm_H made in KR90 (see the discussion in Sect. 2.1.4)

3. Same as 1, but considering the area term dA/dΦ to be constant and ignoring the term linear in ΔΦ of Eq. (8).

4. Using the prescription of KR90 without modification and without outflows from the outer Lagrangian point.

Figure 11 shows the evolution in the HR diagram for both the P_i = 31.6 day and P_i = 1000 day system using our updated MT prescription. In both cases, MT first proceeds rapidly, stripping the majority of the hydrogen envelope in a few thousand years, followed by a longer phase (∼10⁵ yr) driven by the evolution of the core during core-helium burning. In both cases, MT lowers the mass of the donor to ∼13 M_⊙ and carbon is depleted while still undergoing RLOF. The 31.6 day system depletes carbon while still retaining hydrogen in its surface (surface mass fraction of X_s = 0.40), with a radius of 8.1 R_⊙ and an orbital period of 1.8 days. After it first interacts, the 1000 day system always remains near or at RLOF, reaching carbon depletion with a surface hydrogen mass fraction X_s = 0.36, a radius of 7.5 R_⊙, and an orbital period of 1.5 days. The resulting binary BHs have large delay times for a merger due to GW radiation, with the 31.6 day system being expected to merge after 6.4 Gyr and the 1000 day systems after 6.0 Gyr.

Fig. 11. Evolution in the HR diagram for a 30 M⊙ donor star in two systems leading to the formation of a merging binary BH through stable MT. Dots in the tracks indicate intervals of 105 years, while thicker lines indicate periods of RLOF.

Focusing on the 31.6 day system, which we do not predict to overflow its outer Lagrangian point, Fig. 12 shows the evolution of the MT rate and the Roche filling factor with the variations to the MT prescription described before. For our updated prescription, and in the absence of outflows from the outer Lagrangian point, the computed MT rate is expected to be higher than that of KR90. This is because dropping the assumption of P/ρ = kT/μm_H of KR90 and accounting for the linear term in dA/dΦ increases the value of the integrand in Eq. (24). As shown in Fig. 12, we do find that both of these assumptions result in a lower MT rate and a higher amount of RLOF compared to our full MT scheme. This is very noticeable in the slower second MT phase, where our new MT prescription predicts a maximum Roche filling factor of R/R_RL = 1.16, while the standard prescription from KR90 reaches R/R_RL = 1.37 during this phase. For this example system, we find that the largest difference from the KR90 prescription comes from dropping the assumption that P/ρ = kT/μm_H, while the contribution of the linear term in dA/dΦ has a much smaller impact. Without our updated MT prescription, the KR90 prescription predicts that the system overflows its outer Lagrangian point.

Fig. 12. Evolution of the MT rate and overflow for a system with a 30 M⊙ donor star, a 7.5 M⊙ BH, and an initial orbital period of 31.6 days. Top: ratio of the computed MT rate under different assumptions, compared to the result from our full MT prescription. Middle: ratio of the stellar radius to the Roche lobe radius and to the volume equivalent radius of the L2 and L3 equipotentials (for the L2 equipotential, we only include the results of our default model for ease of visualization). Bottom: MT rate.

The case for the 1000 day system, which undergoes overflow of its outer Lagrangian point in our model, is shown in Fig. 13. In this case, we also find that our updated prescription keeps the star significantly more compact, with the KR90 prescription reaching a Roche filling factor of R/R_RL = 1.48, while our model reaches R/R_RL = 1.34. With the updated MT scheme, the donor reaches an overflow of the outer Lagrangian point of R/R_L out = 1.04, resulting in an outflow from L₃ that at its maximum accounts for less than 9% of the total mass loss rate of the star. Since most mass is still transferred through L₁, the additional angular momentum loss from mass loss through L₃ only produces a small change in the final orbital separation, with the model using the prescription of KR90 resulting in a final orbital period of 1.6 days compared to the 1.5 day result when using our updated MT prescription. For this binary system, we find that during overflow of the outer Lagrangian point including the linear term in dA/dΦ has an effect comparable to that of assuming P/ρ = kT/μm_H. As discussed in Appendix A, it is indeed expected that the linear term in dA/dΦ has a larger impact on mass outflows through L₃ than from L₁.

Fig. 13. Same as Fig. 12, but for a system with a 30 M⊙ star with a 4.5 M⊙ BH companion at an initial orbital period of 1000 days.

As the specific angular momentum that would be carried by material being ejected through the outer Lagrangian point is uncertain (we assume by default that it corresponds to the location of the Lagrangian point), we also repeated the 1000 day simulation but increased the angular momentum loss rate associated with this effect by factors of 4 and 9. In both cases, we still find the system to be stable, but they result in much shorter final periods of 1.2 day and 0.83 day. This corresponds to delay times of 3.6 Gyr and 1.2 Gyr, compared to the 6.0 Gyr delay time for our default model. This additional uncertainty could then allow for stable MT to produce systems with short delay times, where tidal spin-up can result in a high spin second formed BH (cf. Bavera et al. 2020). Another potential source of mass and angular momentum loss is the occurrence of L₂ outflows when M_a < M_d for which our default model does not account. For the 100 day system, we find that the donor slightly overfills the L₂ equipotential, exceeding the L₂ equivalent radius by 0.2%. This small amount of overflow is not expected to play a significant role, and indeed in our experiments with L₂ outflows shown in Fig. 10, we do not see an important difference. For the 1000 day simulation, however, our default model overflows the L₂ equivalent radius by 21%. Due to this, in the simulations with L₂ outflows shown in Fig. 10, we find that the outcome is sensitive to potential outflows through L₂, resulting in CE evolution and a merger if the specific angular momentum of material ejected through L₂ is large.

In summary, we find that our modifications to the MT prescriptions of R88 and KR90 do not have a significant effect on the MT rate or the final masses after MT concludes. The changes, however, do affect the size of the donor through the MT phase, resulting in more compact stars with less overflow during MT. Systems where we predict overflow of the outer Lagrangian point result in smaller orbital periods than when using the prescriptions of R88 and KR90, but we do not find that the additional source of angular momentum loss has a significant impact on the stability of these systems.

3.3. Common envelope in donors with convective envelopes

To illustrate the result of our simulations that survive CE evolution, we consider a particular model that results in successful ejection in both our α_CE = 1 and α_CE = 0.1 grids, which has a 30 M_⊙ star with a 14.1 M_⊙ BH in an initial orbit of 2344 days. This model undergoes RLOF after core helium depletion and the evolution of their stellar and Roche lobe radii is shown in Fig. 14 for both CE efficiencies considered. Before the onset of RLOF, wind mass loss reduces the mass of the donor to 28.2 M_⊙ and before hitting the threshold value of Ṁ_high = 1 M_⊙ yr⁻¹ at which we started modeling the system as a CE, the mass of the donor is reduced to 24.7 M_⊙. The onset of CE happens after the star slightly overflows its outer Lagrangian point (R/R_L out = 1.08), at which point our MT prescription predicts 53% of the mass to be lost by the star going through L₃. During the CE calculation, we artificially removed the mass and adjusted the separation according to the binding energy of the removed layers, resulting in a sharper decrease in mass of the Roche lobe radius for the lower CE efficiency. The stellar radius remained larger than 10³ R_⊙ until the mass of the star was lowered below 15 M_⊙, at which point the radius began to decrease sharply with the mass. The point at which we found detachment without any enforced mass loss from CE, see Eq. (34), varies with efficiency, with the α_CE = 1 simulation detaching at a mass of 14.3 M_⊙ and the α_CE = 0.1 model detaching at 13.7 M_⊙.

Fig. 14. Evolution of the radius and Roche lobe radius for two models undergoing successful envelope ejection from CE evolution. Both models correspond to a simulation of a 30 M⊙ star with a 14.1 M⊙ BH with an initial orbital period of 2344 days, with the only difference between the simulations being the efficiency of CE, αCE = 1 or 0.1. The shaded blue region indicates pre-CE evolution which is identical for both cases. For each case, empty circles denote the moment our algorithm determines that CE evolution would end as the star contracts within its Roche lobe in the absence of mass loss from CE evolution.

Figure 15 shows the binding energy profiles of the star at different phases of evolution. Between the beginning of RLOF and the onset of CE evolution, the donor loses ∼3 M_⊙ of mass which modifies the binding energy of the convective envelope but not of the deeper layers. The binding energy of the ejected envelope changes significantly between the α_CE = 1 and α_CE = 0.1 simulations with E_bind = −3.9 × 10⁴⁷ erg and E_bind = −4.3 × 10⁴⁸ erg, respectively. For the fully efficient simulation, all ejected mass formed part of the convective envelope at the onset of CE. This factor of about 10 difference on the binding energy, coupled with the factor of 10 difference in CE efficiency leads to very different outcomes post-CE. The α_CE = 1 simulation terminates CE evolution with an orbital separation of 570 R_⊙, while the α_CE = 0.1 simulation finishes CE with an orbital separation of 8.4 R_⊙. As discussed by Deloye & Taam (2010), the core-envelope boundary is not uniquely defined, resulting in different values depending on the mass ratio and efficiency of the process. One caveat of our simulations is that the amount of mass loss previous to CE evolution is sensitive to our choice of Ṁ_high, which modifies the binding energy of the outer layers of the star (see Appendix B). This modifies the post-CE separation for systems we predict to remain wide, such as the α = 1 system discussed here. Accurate predictions for these wide post-CE models requires a better description of how the star transitions from stable MT to CE evolution.

Fig. 15. Binding energy profiles of the donor in a system with an initial mass of 30 M⊙ and a BH companion of 14.1 M⊙ at an orbital period of 2344 days. Profiles are shown at RLOF, at the onset of CE evolution, and after detachment from CE for αCE = 1 and 0.1. Vertical lines indicate the final mass after CE evolution at the two efficiencies considered.

After detachment, for both efficiencies we find that the donor undergoes a thermal pulse and re-expands, filling its Roche lobe but undergoing stable MT rather than developing a CE again. Figure 16 shows the evolution in the HR diagram both before and after CE evolution. For the case of fully efficient CE evolution, during this additional phase of MT the star would appear as an extended blue or yellow supergiant with a surface hydrogen mass fraction of X_s ∼ 0.4. The duration of this MT phase for the case of α_CE = 1 is ∼2000 years, during which mass is transferred at a rate above the Eddington limit of the 14.1 M_⊙ BH companion. For the α_CE = 0.1 case, the donor star has a radius of just 2 R_⊙ with a surface hydrogen mass fraction of X_s ∼ 0.3, and it transfers mass at super-Eddington rates for ∼800 years. Given that the Eddington luminosity of a 14.1 M_⊙ BH is ∼2 × 10³⁹ erg s⁻¹, during this final MT phase the system can potentially appear as an ultra-luminous X-ray source (see Kaaret et al. 2017 for a recent review), but owing to the short time before core carbon depletion, the number of observable sources could be negligible compared to those that formed by other evolutionary pathways.

Fig. 16. HR diagram showing the evolution of a system undergoing CE evolution composed of a 30 M⊙ donor with a 14.1 M⊙ BH companion at an initial period of 2344 days. Results are shown for two values of the efficiency of CE evolution, αCE = 1 and αCE = 0.1, with the evolution prior to CE being identical.

Not all our simulations that survive CE evolution have short lifetimes as X-ray sources. For the system we have discussed so far, this is the case because CE evolution happens after core helium depletion. If we consider a system with a 7.5 M_⊙ BH at an initial orbital period of 2000 days and α_CE = 1, we find that the onset of CE happens when the donor has a central helium mass fraction of Y_c = 0.02, and for the entirety of its post-CE lifetime until core carbon depletion (26 000 years), it is a Roche-lobe filling binary. This matches the results of Quast et al. (2019) who argue that post-CE systems can undergo MT on the nuclear timescale of the donor. We expect that in cases where the star does not halt its expansion in the HG, CE evolution is initiated with a much larger value of Y_c, leading to a longer lifetime as an X-ray source post-CE.

The single model we considered with different values of α_CE illustrates how the binding energy is dependent on the efficiency of CE evolution, but using the models from our grid we can see how the resulting binding energy of systems surviving CE evolution depends on the mass ratio as well. Figure 17 illustrates the binding energy as a function of the radius at RLOF resulting from models at three different initial mass ratios, α_CE = 1, and varying initial orbital periods. For most of these simulations, we find that CE evolution terminates near or above the mass coordinate of the bottom of the convective envelope at RLOF, resulting in binding energies < 10⁴⁸ erg. These low binding energies are below those predicted by the fits of Xu & Li (2010) and Claeys et al. (2014), although for a mass ratio of q = 0.17 we find binary models with binding energies between the value of these two fits. The q = 0.17 simulations with binding energies > 10⁴⁸ erg behave similarly to the simulation with α_CE = 0.1 shown in Fig. 15, where the binary only detaches after removing layers with a high binding energy that were radiative at the onset of CE evolution. This results in an order of magnitude variation of the binding energy with a varying mass ratio.

Fig. 17. Binding energy of the envelope of a star with an initial mass of 30 M⊙ as a function of the radius at RLOF. Solid lines indicate the results from our αCE = 1 grid with three different mass ratios q = 0.17, 0.29, and 0.45. Dotted lines show the binding energy at a given radius deduced from a single star model, including the results obtained with the fit of Xu & Li (2010), the fit of Claeys et al. (2014), and assuming the bottom of the envelope to be at X = 0.1, X = 0.3, or to correspond to the bottom of the convective envelope. Symbols indicate whether the binary undergoes RLOF before or after core helium depletion.

3.4. Common envelope in donors with radiative envelopes

As all of our simulations where CE evolution is triggered for a star with a radiative envelope are expected to merge during CE evolution, we cannot properly define a core-envelope boundary to compute the binding energy. However, we can artificially produce systems that eject their outer envelope by using artificially high values of α_CE. Here, we consider systems consisting of a 30 M_⊙ donor and a 3 M_⊙ BH with initial orbital periods between 10¹ days and 10³ days, and values of α_CE of 10 and 20. The low mass of the BH, which would fall in the putative lower BH mass gap (Özel et al. 2010; Farr et al. 2011; Kreidberg et al. 2012; Abbott et al. 2020b), was chosen such that MT would be unstable across the range of chosen orbital periods.

All of these simulations managed to eject their envelopes during CE evolution, and Fig. 18 shows the resulting masses post-CE as well as the binding energy of the ejected layers. We find that CE terminates when the surface hydrogen mass fraction is relatively high (X ∼ 0.4−0.5) and that the binding energy of the envelope varies between 1.3 × 10⁵⁰ erg and 6.6 × 10⁴⁹ erg, a factor of 2 variation despite a change by a factor of 20 in radius at RLOF. Figure 18 also illustrates that whether RLOF occurs before or after the end of the HG has no impact on the binding energy. Estimating the binding energy using either the Claeys et al. (2014) or Xu & Li (2010) fits would underestimate the binding energy of the envelope for all these simulations. Since the binding energy just changes by a factor of 2 despite over an order of magnitude change in the radius of the star, population synthesis calculations that approximate E_bind using a fixed value of the dimensionless λ parameter cannot reproduce the variation of E_bind with radius, independent of the choice of λ.

Fig. 18. Properties after CE evolution for the case of donor stars with a radiative envelope, with artificially high values of αCE = 10, 20 in order to allow for a successful CE ejection. Simulations correspond to binary systems composed of a star with a ZAMS mass of 30 M⊙ and a BH companion of 3 M⊙, with periods ranging from 101 days to 103 days. Results are shown as a function of the radius of the star when it fills its Roche lobe. Top panel: central helium mass fraction Yc at the onset of CE, showing RLOF before and after the end of the HG as a large drop in Yc. Second panel: surface hydrogen mass fraction after envelope ejection. Third panel: stellar mass after envelope ejection, also showing the mass of the helium core defined as the innermost mass coordinate with Y > 0.01. Bottom panel: binding energy of the ejected envelope, as well as the binding energy that would be predicted by the Claeys et al. (2014) and Xu & Li (2010) fits, and the one corresponding to λ = 1.

4. Discussion

We have studied the outcomes from the interaction of low-metallicity massive stars with a companion BH using detailed binary evolution calculations performed with the MESA code. Our simulations, which consist of a 30 M_⊙ low-metallicity donor with a companion BH, show several features relevant to understand the formation of merging binary BHs through both stable MT and CE evolution. For these simulations, we have also updated the MT prescription developed by KR90, including the possibility of overflow from an outer Lagrangian point, and we also implemented a method to self consistently determine the core-envelope boundary in cases of CE evolution. The main conclusions we derive from our simulations are as follows:

1. Stable MT after the formation of the first BH in a binary can produce merging binary BHs for a broad range of orbital periods. The boundary between stable and unstable MT is located at a point that allows for orbital hardening just from stable MT, without requiring a CE phase (van den Heuvel et al. 2017). This includes cases where MT between the star and the first formed BH is initiated during the main sequence. We also find cases where the low density envelope of the donor would overflow its outer Lagrangian point. Using our model to quantify this outflow, we find that binaries can remain stable despite the high specific angular momentum of this outflow, resulting in more compact binaries after MT than if the process is ignored.

2. We do not find any case of CE evolution initiated with a star with a radiative envelope where the envelope is successfully ejected. Between the TAMS and the start of the rise in the Hayashi line, the binding energy of the envelope remains at ∼10⁵⁰ erg, varying by just a factor of two despite an increase of more than an order of magnitude in radius. Before the formation of a deep convective envelope, we find that the prescriptions of Xu & Li (2010), as developed by Dominik et al. (2012) for M > 20 M_⊙, and Claeys et al. (2014) underestimate the binding energy of the envelope by up to an order of magnitude.

3. The parameter space where systems can survive CE is limited to a small range in the initial orbital period (∼0.2 dex) where RLOF happens after the formation of a deep convective envelope. Assuming that all internal and recombination energy as well as all the energy from inspiral during CE is used to eject the envelope, we find that the majority of these systems surviving CE would form binary BHs too wide to merge within a Hubble time.

4. As in Deloye & Taam (2010), we find that the core-envelope boundary is not uniquely defined, with different CE efficiencies (or similarly, different mass ratios) resulting in different binding energies for the envelope. In particular for the case discussed in Sect. 2.2, we computed binding energies that differ by over an order of magnitude depending on the efficiency of CE evolution. This is in contrast to the way binding energies are determined in population synthesis codes, where a star has a uniquely defined binding energy at each point of its evolution.

5. Under the simple assumption of a flat distribution in the mass ratio and a flat distribution in log P after the formation of the first BH, we find that stable MT would dominate the formation rate of merging binary BHs. If we consider CE efficiency to be 100%, the ratio of binaries surviving CE evolution (independent of whether they form a compact or wide binary BH) to those that form merging binary BHs from stable MT is 0.35. This is further reduced to 0.017 if we take only the CE simulations that produce merging binary BHs. This is in significant contradiction with rapid population synthesis calculations that predict the vast majority (> 90%) of BHs are formed via an initial phase of stable MT before the formation of the first BH, followed by a CE phase when the secondary fills its Roche lobe (cf. Table 4 of Dominik et al. 2012).

These results broadly agree with those of Klencki et al. (2020, 2021), who have used single star models with a fixed definition of the core-envelope boundary to propose that current population synthesis calculations may overestimate the production rate of merging binary BHs from CE evolution. Since our simulations are restricted to a 30 M_⊙ donor at Z = Z_⊙/10, it is necessary to further explore different masses and metallicities to understand how generally our conclusions apply (Gallegos-Garcia et al., in prep.). Recent investigations using rapid population synthesis codes have indicated that given current uncertainties in binary evolution such as the cosmic star formation rate (Neijssel et al. 2019) or the critical mass ratio for stable MT (Andrews et al. 2020), stable MT could be an important production channel for the merging binary BHs observed with gravitational waves. Our results indicate these rapid population synthesis calculations still significantly overestimate the number of systems surviving CE evolution at low metallicities, and accounting for this would further reduce the contribution from the CE evolution channel.

As we find that the boundary between stable and unstable MT is located right at a point that allows for the formation of merging binary BHs via stable MT, any uncertainties that can shift this boundary can have a large effect on the number and properties of expected GW sources. Among such uncertainties, we can point out that our simulations do not include stellar rotation or tidal coupling to the orbit. This additional sink of orbital angular momentum can potentially affect the stability of MT in our calculations. We have explicitly checked in our stable MT simulations that the moment of inertia of the star is too low to undergo the Darwin instability (Darwin 1879), which is the limit were tidal coupling leads to a runaway shrinkage of the orbit. Another important simplification is that we did not account for the effect of the companion in the momentum equation, which would make the adiabatic response of the star to mass loss dependent on the mass ratio of the system and its separation. It has also been argued that hydrodynamics need to be considered to assess the stability of MT for giant donors (Pavlovskii & Ivanova 2015), although a consistent treatment of hydrodynamics needs to consider the perturbation to the gravitational potential from the companion star as well. Additionally, a significant number of the systems we predict that form merging binary BHs also undergo overflow of the L₂ and even the L₃ equipotentials. Although for these large overflow cases, we have verified that merging binary BHs are formed even under extreme assumptions of mass and angular momentum outflows, three-dimensional hydrodynamic studies are required to accurately model outflows in such systems. All these uncertainties need to be studied in detail to determine if stable MT can indeed produce merging binary BHs.

Our simulations for the outcome of CE evolution also have important uncertainties. Although we do not have the core-envelope mass boundary as a free parameter, we still have the CE efficiency as an unknown variable. The efficiency can be reduced by radiative processes or also by a surplus of energy being removed as kinetic energy from the system (cf. Nandez & Ivanova 2016). Within the context of one-dimensional stellar models, it might be possible to constrain the efficiency of the CE process by means of simplified hydrodynamical simulations of the inspiral (Taam et al. 1978; Clayton et al. 2017; Fragos et al. 2019) or, ideally, with full three-dimensional simulations of CE evolution in massive stars (Taam & Sandquist 2000; Ricker et al. 2019; Law-Smith et al. 2020). Despite these uncertainties that significantly affect the outcome of systems that survive CE evolution, we consider our conclusion that massive stars would merge when undergoing CE prior to the formation of a convective envelope as much more certain. We find this result for the case of full efficiency for the CE process, including recombination energy. Our method to determine the core-envelope boundary can also be considered to produce a lower bound on the binding energies. If we were to consider the boundary deeper inside down at the point where thermal pulses would not happen post-CE (Ivanova 2011), we would find larger binding energies and less systems surviving CE evolution.

Our simulations do not exclude CE evolution as an efficient mechanism to form merging binary neutron stars. At lower masses, a larger fraction of the stellar expansion of a star during its lifetime happens after the formation of a deep convective envelope (Fig. 19). This results in a broader range of orbital periods where we expect the ejection of a CE to be possible, leading to binaries compact enough to merge from GW radiation. Following the work of Misra et al. (2020), the boundaries for stable MT between a star and a neutron star companion do not allow for stars massive enough to form a second neutron star in the system. This excludes stable MT as a formation channel for merging neutron stars. The case of MT and CE evolution in the case of neutron star companions will be studied in future work.

Fig. 19. Radial evolution of single stars of various masses at a metallicity of Z⊙/10 modeled until core-carbon depletion. Numbers above each track indicate the mass of the model in solar masses, and symbols indicate different stages in the evolution of the star. The color of the tracks indicates the fraction of the hydrogen envelope that is part of an outer convective zone.

Acknowledgments

We thank the referee for detailed feedback on our work. PM acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. M.G.-G. is grateful for the support from the Ford Foundation Predoctoral Fellowship. CPLB is supported by the CIERA Board of Visitors Professorship. VK is supported by a Senior CIERA Fellowship and by Northwestern University. PM thanks Alina Istrate for useful discussion.

References

Appendix A: Computation of d A/dΦ

The coefficients of the quadratic terms in the Taylor expansion of Eq. (7) are given by

and the coefficients of the quartic terms are

When we consider outflows from the outer Lagrangian point of the donor, is switched by the normalized x-coordinate of the outer Lagrangian point, , which, with our definition, has a negative value.

In order to estimate the area on the plane crossing L₁ for a given ΔΦ, it is useful to use cylindrical coordinates y = ϖ sin θ, z = ϖ cos θ, in which case

(A.1)

where γ ≡ C₁sin²θ + C₂cos²θ. This can be integrated to obtain the area A(Φ),

(A.2)

which after differentiation by Φ gives Eq. (8).

The relative importance of the constant and linear terms for dA/dΦ shown in Eq. (8) can be assessed by using the ratio of dA/dΦ(ΔΦ = 0) and d²A/dΦ²(ΔΦ = 0). This essentially quantifies how large ΔΦ needs to be to produce a significant change in dA/dΦ, with smaller values indicating a higher sensitivity to ΔΦ, and thus a larger importance of the higher-order term. Figure A.1 shows this value at different mass ratios both for the L₁ and L_out points of a star. This value is always positive, indicating that the inclusion of the higher-order term results in an increased value of the computed MT rate. We also see that the value of the ratio is much smaller for the outer Lagrangian point, and in the case where the donor is the more massive component in the binary, a variation of Φ on the order of 0.1 GM_d/a can modify dA/dΦ significantly. This points out that the inclusion of the higher-order term has a larger impact on the computation of outflows from the outer Lagrangian point than from L₁.

Fig. A.1. Ratio of dA/dΦ to d2A/dΦ2 at L1 and at the outer Lagrangian point of a star expressed in units of GMd/a.

Appendix B: Dependence on numerical parameters of CE algorithm

As our algorithm to model CE evolution has various numerical parameters, it is important to evaluate in detail how our results depend upon these. Here we modify the three parameters described in Sect. 2.2: Ṁ_high, Ṁ_low, and δ, for which in our default model we used 1 M_⊙ yr⁻¹, 10⁻⁵ M_⊙ yr⁻¹, and 0.02, respectively. Here, we consider the difference in the outcome of the CE simulations shown in Sect. 3.3 for variations in these three parameters, with the properties of the system post-CE being shown in Table B.1. For all parameters, we find only ∼5% variations in the post-CE properties, except for variations of Ṁ_high when using an efficiency α_CE = 1. Detachment in these systems happens in the proximity of the mass coordinate from the bottom of the convective envelope at the onset of CE. Variations in Ṁ_high modifies the mass of the star at the onset of CE, as is shown in Table B.1. This modifies the binding energy in the outermost layers, which is illustrated in Fig. B.1, making evolution sensitive to the exact point at which we switched between MT and CE evolution. This is not the case for the system with α_CE = 0.1, as the core-envelope boundary is computed to be deeper in the star, in layers that are not affected by the exact point at which we switched between MT and CE evolution.

Fig. B.1. Binding energy profiles of the donor in a system with an initial mass of 30 M⊙ and a BH companion of 14.1 M⊙ at an orbital period of 2344 days. Profiles are shown at RLOF and at the onset of CE evolution for different values of Ṁhigh. Symbols indicate the mass coordinate at which CE evolution terminates for each choice of Ṁhigh and for CE efficiencies of αCE = 1 and αCE = 0.1.

Appendix C: Comparison to other stellar evolution tracks

A comparison of our single metal-poor 30 M_⊙ model to similar stellar evolution tracks that are publicly available is shown in Fig. C.1. These include tracks by Meynet et al. (1994), Pols et al. (1998), Tang et al. (2014), and Klencki et al. (2020). Out of these, the only ones that use the same stellar evolution code are our own models and those of Klencki et al. (2020). From the models of Meynet et al. (1994), we took the 25 M_⊙ stellar track which is the closest one to 30 M_⊙. The track taken from Pols et al. (1998) corresponds to their models with overshooting, and the one from Klencki et al. (2020) is from their set of nonrotating stellar tracks. All tracks experience a long-lived phase as blue supergiants, halting their rapid expansion in the HG phase at log₁₀(T_eff/K) ∼ 4.2 (except the lower mass one from Meynet et al. 1994 that does so at a higher T_eff).

Fig. C.1. Evolution of stars with a mass of ∼30 M⊙ and a metallicity of Z⊙/10 computed by different research groups and codes. Dots are plotted every 105 years.

Appendix D: Single star resolution test

Given the formation of multiple intermediate convection zones in our models during the HG phase, we performed a resolution test to see if our results are significantly modified by increased spatial and temporal resolution. We made use of the mesh_delta_coeff and time_delta_coeff options that allow for an approximate scaling of all spatial and temporal resolution controls, and used models with approximately two, four, and eight times the resolution of our default setup. The resulting tracks in the HR diagram are shown in Fig. D.1. We see small variations with changes in resolution on the evolution during the HG phase, producing a variation of about 10% on the effective temperature at which the HG phase ends.

Fig. D.1. Evolution in the HR diagram of a single 30 M⊙ model with a metallicity of Z⊙/10 computed with increasing spatial and temporal resolution. For each track, we indicate the average number of zones in the model as well as the steps taken in the simulation. The lowest resolution track shown corresponds to our default setup.

In particular, for this work we are interested in how the binding energy of the envelope of a star compares to prescriptions present in the literature. Figure D.2 shows the binding energies obtained for models at increasing temporal and spatial resolution, as well as the results we obtain with the Xu & Li (2010) and Claeys et al. (2014) fits. Binding energies were computed from the stellar models taking the bottom of the envelope at the innermost mass coordinate where X > 0.1 and X > 0.3. Between the different models, we see variations in the computed binding energies at a given radius smaller than 12%. The larger variations happen at the end of the HG and are due to the slightly different radii achieved at the end of the HG. Ignoring the values between log₁₀(R/R_⊙) = 1.8−1.95, we find that the differences are smaller than 6%.

Fig. D.2. Binding energy of the envelope of a post-main-sequence 30 M⊙ model with a metallicity of Z⊙/10 computed with increasing spatial and temporal resolution. For each track we indicate the average number of zones in the model as well as the steps taken in the simulation. The lowest resolution track shown corresponds to our default setup. For comparison, we include the binding energies obtained using our default resolution together with the fits of Xu & Li (2010) and Claeys et al. (2014). Before the formation of a convective envelope, the binding energies computed at X = 0.1 and X = 0.3 almost overlap.

Appendix E: Application of λ fits to stellar models

The fits from Claeys et al. (2014) and Xu & Li (2010) depend on the mass of the star, its radius, and its evolutionary stage in terms of the star-types defined by Hurley et al. (2002). A model analogous to our 30 M_⊙ model is assigned types “HG” during the HG, “CHeB” during core helium burning but before ascending the Hayashi line, and “EAGB” while on the Hayashi line. To distinguish between the CHeB and EAGB types in our models, we considered models with convective envelopes > 1 M_⊙ to be of type EAGB, as we found this is close to the point where the star rises on the Hayashi line. As we are only concerned with these three assigned types, we only describe the fits associated with them.

E.1. Claeys et al. (2014)

The fits of Claeys et al. (2014) were first presented by Izzard (2004), who took them directly from the source code of the BSE code (Hurley et al. 2002). Claeys et al. (2014) do not specify the range of masses used to compute these fits. The fit for λ including internal but not recombination energy is

(E.1)

where

(E.2)

with R_ZAMS being the radius of the stellar model at the ZAMS. The value of λ₁ is given by

(E.3)

with

(E.4)

This fixes a typo (0.9 instead of −0.9 in the minimum function) from Claeys et al. (2014) and Izzard (2004) that is inconsistent with the implementation in the BSE code from which the fit was taken. Claeys et al. (2014) also describe variations in the λ₁ and λ₃ values for stellar types that are not relevant for our 30 M_⊙ model, so we do not describe those.

E.2. Xu & Li (2010)

The λ fits of Xu & Li (2010) cover up to a mass of 20 M_⊙, while fits for higher masses were developed by Dominik et al. (2012). These fits are not published, but can be found in the source code of the COMPAS code. Fits are provided for the binding energy with and without thermal energy included, that is α_th = 0 or 1 in Eq. (33). As we worked with α_th = 1, we compared our results to the λ fit including thermal energy. Fits are provided individually for each of the HG, CHeB, and EAGB phases, and the ones that correspond to our 30 M_⊙ model at Z_⊙/10 are the Population II fits for masses 18 M_⊙ < M < 35 M_⊙. For HG stars λ is computed as

(E.5)

where ξ = R/R_⊙ and the value of the fitting coefficients is given in Table E.1. For CHeB stars we have

(E.6)

with the expression for the range 755 R_⊙ < R < 900 R_⊙ being a linear interpolation from 0.1 to 0.2 with b = (900 − ξ)/(900−755). For the case of EAGB stars the fit is given by

(E.7)

where b = (900 − ξ)/(900−850).

All resulting values are capped between 0.05 and 1.5. The interpolation in the range 755 R_⊙ − 900 R_⊙ for CHeB stars is missing in the implementation of Dominik et al. (2012), and it erroneously applies the fitting polynomial in this range resulting in negative values, producing a positive value only because the final result is capped to a minimum of 0.05. Not including this correction and having λ = 0.05 in this range of radii does not change our conclusions, as the binding energy is still underestimated compared to our own models.

All Tables

All Figures

Fig. 1. Definition of the coordinate system and variables used in the calculation of MT. The coordinate system is centered in the donor star, with the x-axis defined as the line joining the donor and the accretor, with the accretor located at x = a where a is the orbital separation. The z-axis is oriented along the direction of the orbital angular momentum of the system, while the y-axis completes a standard right-handed oriented Cartesian set of coordinates. The Lagrangian point located behind the donor is either L3 or L2, depending on whether the donor is more or less massive than the accretor, respectively.

In the text

	Fig. 2. Ratio between the volume equivalent radius associated with the outer Lagrangian point of the donor and its Roche lobe radius. Bottom panel: result from numerical integration of the Roche potential together with the fit given in Eq. (3). Top panel: relative error between the fit and the data, with the small variations being due to the finite precision of the volume calculation.
In the text

Fig. 3. Top: specific angular momentum in units of a2Ωorb for the donor, the accretor, and the first three Lagrangian points as a function of the mass ratio. Bottom: volume equivalent radii for the donor corresponding to the L2 and L3 equipotentials. Volume equivalent radius for L2 is taken from the fit of Marchant et al. (2016), which due to small errors in the fit results in a slightly larger value than the radius for the L3 equipotential at q = 1.

In the text

Fig. 4. Evolution of a 30 M⊙ single star with a metallicity of Z⊙/10, separated into different evolution phases. Dots indicate equal time intervals of 105 years. The Hertzsprung gap is defined as the evolutionary stage between the terminal-age main sequence and the moment where the ratio of the nuclear burning luminosity and the surface luminosity satisfies |log10(Lnuc/L)| < 0.001 for at least a thermal timescale estimated as GM2/RL, or the star develops a convective envelope > 1 M⊙. Top: evolution under our default set of physical inputs. Middle: evolution with core overshooting increased from 0.335 to 0.535 pressure scale heights. Bottom: evolution with wind mass loss rates tripled.

In the text

Fig. 5. Kippenhahn diagram for our default 30 M⊙ model, separating different evolutionary stages as well as showing the evolution of the stellar radii. “He dep.” and “C dep.” stand for core helium and carbon depletion respectively. In addition, we show the result of computing Ebind under the assumption that the core–envelope boundary Mcore at which CE would finish is at the innermost mass coordinates where the hydrogen mass fraction is X = 0.1 or X = 0.3. We only show Ebind after the formation of the helium core, so no information is given on the main sequence. To provide a comparison point to rapid population synthesis codes, we also show the resulting Ebind from the fits of Claeys et al. (2014) and Xu & Li (2010).

In the text

	Fig. 6. Binding energy of the envelope of a 30 M⊙ post main-sequence star as a function of radius. Tracks are shown for our default set of physical assumptions, as well as for a model with increased overshooting and one with increased mass loss rates. For each model, we also show the binding energy resulting from the fits of Xu & Li (2010) and Claeys et al. (2014).
In the text

Fig. 7. Radius (top) and binding energy (bottom) profiles for a 30 M⊙ single star at different evolutionary stages. A profile near the end of core helium-burning (Yc = 0.08), but before the formation of a convective envelope, is included. Also, a profile is shown corresponding to the beginning of the formation of a deep convective envelope, taken as the first point in the evolution where the mass of the convective envelope Mconv env is larger than 1 M⊙.

In the text

Fig. 8. Summary of final outcomes for simulations of circular binaries consisting of a BH and a 30 M⊙ star at a metallicity of Z⊙/10, and a CE efficiency parameter αCE = 1. Horizontal dotted lines indicate boundaries for an interaction before different evolutionary stages of the star. For systems that undergo stable MT, black horizontal lines indicate regions where the donor exceeded the L2 equipotential, while cross hatched regions mark regions where also the L3 equipotential is exceeded. Black rectangles indicate systems that do not interact or would result in a Roche lobe filling system at ZAMS. Systems that undergo stable MT or eject their envelope during CE evolution are denoted by “st. MT” and “CE ej.”, respectively, and they are separated into systems forming binary BHs that would merge in less or more than a Hubble time. Systems marked as “CE merger” merge during the CE phase.

In the text

	Fig. 9. Same as Fig. 8, but for CE efficiency parameters of αCE = 0.4, 0.2, and 0.1.
In the text

Fig. 10. Same as Fig. 8, but for a subset of models computed using the method described in Sect. 2.1.6 that considers outflows from L2 rather than L3 when Ma < Md. Each panel indicates simulations performed with different specific angular momentum assumed for the material outflowing through L2, which is considered to be jL2, 4jL2, or 9jL2. The purple dotted line indicates the boundary between stable MT and CE evolution as determined by our models that account only for outflows from the outer Lagrangian point of the donor.

In the text

	Fig. 11. Evolution in the HR diagram for a 30 M⊙ donor star in two systems leading to the formation of a merging binary BH through stable MT. Dots in the tracks indicate intervals of 105 years, while thicker lines indicate periods of RLOF.
In the text

Fig. 12. Evolution of the MT rate and overflow for a system with a 30 M⊙ donor star, a 7.5 M⊙ BH, and an initial orbital period of 31.6 days. Top: ratio of the computed MT rate under different assumptions, compared to the result from our full MT prescription. Middle: ratio of the stellar radius to the Roche lobe radius and to the volume equivalent radius of the L2 and L3 equipotentials (for the L2 equipotential, we only include the results of our default model for ease of visualization). Bottom: MT rate.

In the text

	Fig. 13. Same as Fig. 12, but for a system with a 30 M⊙ star with a 4.5 M⊙ BH companion at an initial orbital period of 1000 days.
In the text

Fig. 14. Evolution of the radius and Roche lobe radius for two models undergoing successful envelope ejection from CE evolution. Both models correspond to a simulation of a 30 M⊙ star with a 14.1 M⊙ BH with an initial orbital period of 2344 days, with the only difference between the simulations being the efficiency of CE, αCE = 1 or 0.1. The shaded blue region indicates pre-CE evolution which is identical for both cases. For each case, empty circles denote the moment our algorithm determines that CE evolution would end as the star contracts within its Roche lobe in the absence of mass loss from CE evolution.

In the text

	Fig. 15. Binding energy profiles of the donor in a system with an initial mass of 30 M⊙ and a BH companion of 14.1 M⊙ at an orbital period of 2344 days. Profiles are shown at RLOF, at the onset of CE evolution, and after detachment from CE for αCE = 1 and 0.1. Vertical lines indicate the final mass after CE evolution at the two efficiencies considered.
In the text

	Fig. 16. HR diagram showing the evolution of a system undergoing CE evolution composed of a 30 M⊙ donor with a 14.1 M⊙ BH companion at an initial period of 2344 days. Results are shown for two values of the efficiency of CE evolution, αCE = 1 and αCE = 0.1, with the evolution prior to CE being identical.
In the text

Fig. 17. Binding energy of the envelope of a star with an initial mass of 30 M⊙ as a function of the radius at RLOF. Solid lines indicate the results from our αCE = 1 grid with three different mass ratios q = 0.17, 0.29, and 0.45. Dotted lines show the binding energy at a given radius deduced from a single star model, including the results obtained with the fit of Xu & Li (2010), the fit of Claeys et al. (2014), and assuming the bottom of the envelope to be at X = 0.1, X = 0.3, or to correspond to the bottom of the convective envelope. Symbols indicate whether the binary undergoes RLOF before or after core helium depletion.

In the text

Fig. 18. Properties after CE evolution for the case of donor stars with a radiative envelope, with artificially high values of αCE = 10, 20 in order to allow for a successful CE ejection. Simulations correspond to binary systems composed of a star with a ZAMS mass of 30 M⊙ and a BH companion of 3 M⊙, with periods ranging from 101 days to 103 days. Results are shown as a function of the radius of the star when it fills its Roche lobe. Top panel: central helium mass fraction Yc at the onset of CE, showing RLOF before and after the end of the HG as a large drop in Yc. Second panel: surface hydrogen mass fraction after envelope ejection. Third panel: stellar mass after envelope ejection, also showing the mass of the helium core defined as the innermost mass coordinate with Y > 0.01. Bottom panel: binding energy of the ejected envelope, as well as the binding energy that would be predicted by the Claeys et al. (2014) and Xu & Li (2010) fits, and the one corresponding to λ = 1.

In the text

	Fig. 19. Radial evolution of single stars of various masses at a metallicity of Z⊙/10 modeled until core-carbon depletion. Numbers above each track indicate the mass of the model in solar masses, and symbols indicate different stages in the evolution of the star. The color of the tracks indicates the fraction of the hydrogen envelope that is part of an outer convective zone.
In the text

	Fig. A.1. Ratio of dA/dΦ to d2A/dΦ2 at L1 and at the outer Lagrangian point of a star expressed in units of GMd/a.
In the text

	Fig. B.1. Binding energy profiles of the donor in a system with an initial mass of 30 M⊙ and a BH companion of 14.1 M⊙ at an orbital period of 2344 days. Profiles are shown at RLOF and at the onset of CE evolution for different values of Ṁhigh. Symbols indicate the mass coordinate at which CE evolution terminates for each choice of Ṁhigh and for CE efficiencies of αCE = 1 and αCE = 0.1.
In the text

	Fig. C.1. Evolution of stars with a mass of ∼30 M⊙ and a metallicity of Z⊙/10 computed by different research groups and codes. Dots are plotted every 105 years.
In the text

	Fig. D.1. Evolution in the HR diagram of a single 30 M⊙ model with a metallicity of Z⊙/10 computed with increasing spatial and temporal resolution. For each track, we indicate the average number of zones in the model as well as the steps taken in the simulation. The lowest resolution track shown corresponds to our default setup.
In the text

Fig. D.2. Binding energy of the envelope of a post-main-sequence 30 M⊙ model with a metallicity of Z⊙/10 computed with increasing spatial and temporal resolution. For each track we indicate the average number of zones in the model as well as the steps taken in the simulation. The lowest resolution track shown corresponds to our default setup. For comparison, we include the binding energies obtained using our default resolution together with the fits of Xu & Li (2010) and Claeys et al. (2014). Before the formation of a convective envelope, the binding energies computed at X = 0.1 and X = 0.3 almost overlap.

In the text