The Complete Evolution of a Neutron-Star Binary through a Common Envelope Phase Using 1D Hydrodynamic Simulations

TL;DR

This study models the entire common envelope phase of a neutron star and red supergiant using 1D hydrodynamic simulations, revealing how the envelope is ejected and predicting the formation of a close double neutron-star system.

Contribution

It provides a detailed, self-consistent simulation of the common envelope evolution with a neutron star using MESA, including drag forces and envelope dynamics, which is a novel approach.

Findings

Most of the hydrogen envelope is ejected, forming a large, dusty, optically thick envelope.

The final binary consists of a 2.6 solar mass helium star and a neutron star with a separation of 3.3-5.7 solar radii.

The estimated efficiency of envelope ejection is approximately 5, indicating highly efficient common envelope evolution.

Abstract

Over forty years of research suggests that the common envelope phase, in which an evolved star engulfs its companion upon expansion, is the critical evolutionary stage forming short-period, compact-object binary systems, such as coalescing double compact objects, X-ray binaries, and cataclysmic variables. In this work, we adapt the one-dimensional hydrodynamic stellar evolution code, MESA, to model the inspiral of a 1.4M$_{\odot}$ neutron star (NS) inside the envelope of a 12$M_{\odot}$ red supergiant star. We self-consistently calculate the drag force experienced by the NS as well as the back-reaction onto the expanding envelope as the NS spirals in. Nearly all of the hydrogen envelope escapes, expanding to large radii ($\sim$10$^2$ AU) where it forms an optically thick envelope with temperatures low enough that dust formation occurs. We simulate the NS orbit until only 0.8M$_{\odot}$ of the hydrogen envelope remains around the giant star's core. Our results suggest that the inspiral will continue until another $\approx$0.3M$_{\odot}$ are removed, at which point the remaining envelope will retract. Upon separation, a phase of dynamically stable mass transfer onto the NS accretor is likely to ensue, which may be observable as an ultraluminous X-ray source. The resulting binary, comprised of a detached 2.6M$_{\odot}$ helium-star and a NS with a separation of 3.3-5.7R$_{\odot}$, is expected to evolve into a merging double neutron-star, analogous to those recently detected by LIGO/Virgo. For our chosen combination of binary parameters, our estimated final separation (including the phase of stable mass transfer) suggests a very high $\alpha_{\rm CE}$-equivalent efficiency of $\approx$5.

Equations

$$\Delta E_{grav}={\left(\frac{d{E}_{grav}}{d{E}_{orbit}}\right)}_f\Delta E_{orbit},$$

$$E_{grav}(m_{c,0})-E_{grav}(m_c)=$$

$${\left(\frac{d{E}_{grav}}{d{E}_{orbit}}\right)}_f\left(\frac{G{m}_c{m}_{NS}}{2r_c}-E_{orbit,0}\right),$$

$$r_c=\frac{1}{2}\frac{G{m}_c{m}_{NS}}{{\left(\frac{d{E}_{orbit}}{d{E}_{grav}}\right)}_f\left[E_{grav}(m_{c,0})-E_{grav}(m_c)\right]+E_{orbit,0}}.$$

Full text

stars: binaries: close, stars: evolution, X-rays: binaries

††journal: AAS Journals

1 Introduction

The recent detection of gravitational waves from coalescing binary black holes and most recently of a binary neutron star (NS) merger (The LIGO Scientific Collaboration et al., 2018, and references therein) and the associated gamma-ray burst and kilo-nova explosions (Abbott et al., 2017), sparked a renewed interest in the formation of compact-object binaries. For two NSs to merge within a Hubble time, they must have an orbital separation of \lesssim 10$${R_{\odot}}, yet on their way to NS formation, all massive stars go through a supergiant phase in which they expand to radii of \sim 500-1000$${R_{\odot}}, two orders of magnitude larger than their ultimate orbital size. The common envelope (CE) phase, originally discussed in the context of cataclysmic variable formation (Paczyński, 1976), has been widely adopted by the binary evolution community as the mechanism responsible for forming short-period, compact-object binaries (e.g., Han et al., 1995; Belczynski et al., 2016). For a modern description (and history) of the CE, we refer the reader to the thorough review by Ivanova et al. (2013).

The CE phase typically occurs when the companion of a giant star is engulfed by the giant’s envelope, which now fills the combined binary potential, surrounding both the accretor and the giant star’s core. The frictional torque of the companion orbiting within the fluid of the giant’s envelope transforms orbital energy and angular momentum into heat and spin angular momentum of the CE, dramatically shrinking the orbit. The frictional heat dumped in the envelope causes its expansion and – if there is enough orbital energy available – its eventual expulsion from the system, leaving behind the core of the giant star and its companion in a close orbit.

Quickly after its conception, the CE was studied using numerical hydrodynamics in 1D by Taam et al. (1978); Meyer & Meyer-Hofmeister (1979); Taam (1979). These authors showed some promising initial results; however it was quickly realized that the resulting evolution was not spherically symmetric (Bodenheimer & Taam, 1984), which led to the abandonment of 1D simulations for over 15 years in favor of multi-dimensional simulations to account for the non-axisymmetric geometry of the problem and the turbulent processes involved. Modern simulations employ adaptive mesh refinement techniques (Ricker & Taam, 2008; Taam & Ricker, 2010; Ricker & Taam, 2012; Passy et al., 2012), moving mesh simulations (Ohlmann et al., 2016), or smooth particle hydrodynamics (Nandez et al., 2015; Pejcha et al., 2016a, b). These studies highlight the difficulties of modeling the CE evolution in 3D, as they typically require several million CPU-hours per simulation and still may not include all the necessary physics such as energy and radiation transport and a realistic equation of state. Even without all the necessary physics required, efforts to model this process are hampered by the wide dynamic range in both temporal and spatial scales (Terman et al., 1994; Rasio & Livio, 1996; Taam & Sandquist, 2000).

Building off previous models by Taam et al. (1978), Meyer & Meyer-Hofmeister (1979), Podsiadlowski (2001), and Ivanova (2002) (and more recent work by Ivanova & Nandez, 2016; Clayton et al., 2017), we present one-dimensional simulations of the evolution of a binary through the CE phase, from the onset of the dynamically unstable mass-transfer until the successful ejection of the CE. In this first work, we present the methodology of our approach and describe in detail one test case simulation, relevant to the formation of a coalescing binary NS. This simulation is performed with a modified version of the stellar evolution code MESA (Paxton et al., 2011, 2013, 2015, 2018), taking advantage of the accurate energy transport and detailed microphysics of a stellar structure code as well as the hydrodynamic capabilities implemented in recent code releases. In section 2 we describe our adaptations to MESA. We provide the results and discussion of our simulations in Section 3. Finally, we provide some conclusions in Section 4.

2 Method

In what follows we describe the setup of our simulation and the modifications we have made to the standard MESA code. We use MESA version 9793 and the MESA software development kit version 20170802.

2.1 Model initialization and onset of the dynamical instability

We initialize our binary system, consisting of a 12${M_{\odot}}$ zero-age main sequence star and a 1.4${M_{\odot}}$ point-like companion (representing a NS) at an orbital separation of 972 R*⊙*, using the binary module within MESA. We assume a circular orbit and neglect both wind mass loss and tidal interactions. Therefore, the orbital separation is constant until the star expands to $\approx 550\,\rm R_{\odot}$ as a giant star, overfilling its Roche lobe. The rate of mass transfer due to Roche lobe overflow is calculated using the implicit numerical scheme formulated by Kolb & Ritter (1990), which allows for a potentially significant overflow of the Roche lobe by the donor star. Mass transfer is assumed to be conservative up to the Eddington limit, with the excess material being ejected with the specific orbital angular momentum of the NS.

We continue our simulation until the donor star overfills either of the outer Lagrangian points $L_{\rm 2}$ or $L_{\rm 3}$ (using the fitting formulae for the equivalent $R_{L_{\rm 2}}$ and $R_{L_{\rm 3}}$ radii by Misra et al., in preparation), at which point we consider that the binary enters dynamically unstable mass transfer. The donor’s mass loss rate at the time of $L_{\rm 2}$/$L_{\rm 3}$ overflow onset is $\approx 10^{-2}\rm\,M_{\odot}\,yr^{-1}$, and a few percent of the donor’s mass has already been lost.

2.2 Hydrodynamic modelling of the CE phase

After we resolve the onset of dynamical instability, we switch to the CE MESA module that we developed specifically for this work. Throughout this stage of simulation, we use an energy-conserving implicit hydrodynamic solver which has improved performance when simulating rapid variations (such as shocks) outside the hydrostatic regime (for details see section 4 of Paxton et al., 2015). When using this numerical scheme, the outer boundary conditions are not set by standard stellar atmosphere models, which in any case are inadequate for the problem at hand due to supersonic outflow velocities and very low gas densities developed in the outer part of the CE during the simulation. Instead, to close the stellar structure equations at the outer boundary of the star, we ensure that the compression vanishes ($d\rho/dm=0$) and the temperature is equal to that of a black body for the given luminosity and radius of the shell, and we arbitrarily set the optical depth of the outer shell to $10^{-4}$. In our calculations, we also include corrections to the energy equation due to rotation, although the current implementation of the implicit scheme does not allow analogous rotation corrections to the momentum equation. During our code testing, we encountered numerical instabilities when simulating stellar envelopes of sufficiently low density and temperature which resulted from a transition between the OPAL and the HELM equations of state. To handle this, we use numerical partial derivatives of the internal energy, entropy and gas pressure, instead of the interpolated values provided by the equation of state tables. In the version of MESA used in our simulations this feature was in an experimental stage, but since then, it has been further refined and is now a standard feature in the most recent MESA versions.

Using the above adaptations, we initialize the second phase of our simulation by placing the NS just below the surface of the envelope of the donor star, at a radius equal to 99% of the total donor’s radius. We further set the donor star to have a uniform rotation equal to 95% of the orbital angular frequency. This amounts to the assumption that during the expansion phase of the donor star and the onset of the mass transfer, there was enough time to, almost, synchronize the rotation of the donor’s envelope with the orbit.

We consider the NS to be a point mass moving in a circular orbit with a Keplerian velocity within the envelope of the giant. The NS feels a gravitational drag force causing it to spiral deeper within the giant; it also feels a hydrodynamic drag force and tidal forces, however both 3D hydrodynamic simulations and order of magnitude estimates indicate the gravitational drag force dominates (e.g. Taam et al., 1978; Passy et al., 2012). We add the dissipated energy, determined by the gravitational drag, as an extra heating term within the star’s envelope which is deposited within roughly one accretion radius, $R_{A}$, of the compact object. Each shell within the accretion radius is weighed by a factor $\exp{\left[-(\Delta r/R_{A})^{2}\right]}$, where $\Delta r$ is the radial distance of each mass shell from the NS orbit. The drag force is calculated using the fitting formula by MacLeod, & Ramirez-Ruiz (2015) in the supersonic regime, which provides corrections to the standard expression due to the density gradient in the structure of the star, and by analytic estimates from Ostriker (1999) in the subsonic regime, with the two being smoothly blended for Mach numbers between 0.9 and 1.1. We additionally include accretion luminosity as an energy source, calculated from the Bondi rate but limited to the Eddington rate. This provides a reasonable lower limit to the accretion luminosity as the flow is unable to get to the neutrino cooling regime (MacLeod, & Ramirez-Ruiz, 2015). These simulations find accretion rates somewhat larger than the Eddington-limited rate, but nevertheless only total to $\lesssim$0.1 ${M_{\odot}}$.

At every time step we update the orbital energy based on the gravitational drag, from which we calculate the NS’s position, assuming a circular orbit. Although we conserve energy, since we enforce orbital circularity, we cannot simultaneously conserve angular momentum. In future revisions of our simulations we plan to relax the assumption of a circular binary orbit, which will allow us to also self-consistently spin up the envelope during the inspiral. Here, we still follow the angular momentum transport within the envelope as it restructures, but we do not apply any torque to it due to the NS’s inspiral.

As an additional caveat, we note that we ignore the back-reaction of the NS’s gravitational pull on the giant star’s envelope. Properly accounting for this effect requires multi-dimensional simulations.

Our models do not formally include mass loss, and we track the evolution of the outer envelope even as it expands to radii beyond 104 ${R_{\odot}}$. This allows us to analyze the long-term behavior of the envelope and determine whether it remains truly unbound or whether it goes through oscillation cycles as recently studied by Clayton et al. (2017). Our simulation is terminated due to numerical reasons when the outer layers of the expanding envelope reach temperatures and densities near the edges of our tabulated equation of state.

3 Results and Discussion

3.1 Phenomenological Description

The results of our simulation of a 1.4 ${M_{\odot}}$ NS inspiraling within the envelope of a 12 ${M_{\odot}}$ giant star are shown in Figure 1. The inspiral shows the three distinct phases of CE evolution described in Podsiadlowski (2001) and Ivanova et al. (2013). The initial phase when the NS’s position is roughly constant shows the “loss of co-rotation”. Throughout this phase, the envelope is nearly co-rotating with the engulfed NS, and therefore the weak gravitational drag on the NS leads to only a gradual decay in the NS orbit. Since the thermal timescale from the position of the NS to the surface of the CE is extremely short ($10^{-4}-10^{-2}\,\rm yr$), much of the injected energy quickly radiates away, leaving the envelope relatively unaffected by the presence of the NS.

After $\sim$1.5 years of slow inspiral, when co-rotation is broken, the system enters the “plunge-in” phase in which the NS’s orbit rapidly shrinks from several hundreds of ${R_{\odot}}$ to $\sim$20 ${R_{\odot}}$. Concurrent with the shrinking of the NS’s orbit, the giant star’s envelope shows rapid expansion to $\sim$2000 ${R_{\odot}}$. The expansion of the whole envelope can be seen by following the time evolution of the thin, black, dotted lines in Figure 1 which denote surfaces of constant mass. After the dynamical “plunge-in” phase, the inspiral decelerates and almost stalls during the “self-regulated regime”, which we model for $\simeq$18 years. For clarity, Figure 1 shows only the first 10 years of the simulation.

The top left panel of Figure 1 shows the evolution of the density profile of the donor star. The effect of the NS’s inspiral on the envelope structure is profound; at a radius of 100 ${R_{\odot}}$, the density drops by roughly four orders of magnitude, from $\sim$10*-6* g cm*-3* at the start of the simulation to $\sim$10*-10* g cm*-3* after 10 years. By the end of the simulation, the majority of the envelope’s mass is found in an extended, low density envelope at very large radii ($\sim$104${R_{\odot}}$), with temperatures low enough ($\sim 600\,\rm K$) that dust formation can occur. This is the cause of the increased opacity at the outer layers and the rise of optical depth 10 line closer to the surface after $\simeq$8 yr, as seen in the bottom right panel of Figure 1. In addition, it is worth noting that throughout the self-regulated regime, the NS is found at an almost constant density of $\sim$10*-8* g cm*-3*, just outside a region with a very steep density gradient.

The top right panel of Figure 1 shows the radial velocity profile of the envelope, scaled to the local escape velocity. After the first $\simeq$4 yr, the envelope expansion occurs at a nearly constant velocity of $\simeq$20 km/s. White indicates material leaving at the escape velocity, with red showing material that escapes faster. By the end of the simulation, all the material at large radii ($>$8.5 ${M_{\odot}}$ of the envelope) has a velocity above the escape velocity, indicating the envelope has been successfully ejected.

The bottom left panel of Figure 1 shows the available ionization energy of the envelope. At a given time, the sharp transition in radius from red to white/blue denotes the hydrogen recombination front; the majority of the envelope has recombined by the end of the simulation. At the same time, a helium recombination front exists, which can be seen by the transition between two shades of red, at radii of $\sim 10^{2.5}$ ${R_{\odot}}$ initially and $\sim 10^{2}$ ${R_{\odot}}$ after $\simeq 5\,\rm yr$. The latter occurs at high enough opacities (shown in the bottom right panel of Figure 1) that photons can be reabsorbed by escaping material before they expand to large enough radii, cool and recombine hydrogen as well. There has been much discussion in recent literature addressing the question of whether the radiation released by recombination will be absorbed by the envelope and aid in its ejection (e.g., Ivanova, 2018) or whether the radiation will be quickly transported to the outer layers of the star where it can be radiated away (e.g., Grichener et al., 2018). However, here we note that, as Kruckow et al. (2016) have also demonstrated, in massive stars, the total amount of energy released by recombination adds up to only a small fraction ($\lesssim 10\%$ in our case) of the energy required to eject the envelope.

3.2 Energetics of the inspiral

To evaluate the importance of the various physical processes, we calculate each of their energies throughout the simulation. Figure 2 shows the cumulative energy budget over the entire star as a function of time for each relevant source (solid lines) and sink (dashed lines). Specifically we show the change in potential/gravitational, thermal (which includes energy stored in radiation), ionization and kinetic energy of the star. Furthermore, we track the cumulative energy injected in the envelope due to the gravitational drag onto the NS that converts orbital energy into heat ($\Delta E_{orbit}$) as well as accretion onto the NS ($E_{accretion}$), limited to the Eddington rate. Note that while Figure 2 shows that accretion is not currently an important energy source, if the rate were super-Eddington by a factor of 10 or more, accretion could play a significant role (the reader is referred to MacLeod, & Ramirez-Ruiz, 2015, for a detailed discussion of accretion onto a NS within a common envelope). We plan to explore the effect of super-Eddington accretion in future work. Finally, we track the total energy that is radiated away from the surface of the star and the nuclear energy released. We have checked that the sum of all energy sources and sinks equals approximately zero throughout the evolution, ensuring energy conservation to $\sim 1\%$.

During the first 1.5 yr of the CE evolution, i.e. phase of “loss of co-rotation”, there is little evolution in the energies shown in Figure 2 since the envelope structure remains relatively unchanged. However, during the “plunge-in” phase (between $\simeq$1.5 yr and $\simeq$3 yr) the envelope expands quickly, and its gravitational binding energy increases (becomes less negative) by $\simeq 1.38\times 10^{48}\,\rm erg$, a factor of 2.7 larger than the orbital energy released from the inspiraling NS ($\simeq 0.51\times 10^{48}\,\rm erg$). The excess energy resides in the initial thermal energy content stored in the star, which as the envelope expands, does work and, as shown in Figure 2, releases an additional $\approx 0.84\times 10^{48}\,\rm erg$. Other energy sources and sinks can be ignored, as they add up to $\lesssim 0.05\times 10^{48}\,\rm erg$. The second dotted line in Figure 2 shows the end of the “dynamical plunge-in” phase. The position of the NS within the giant’s envelope where this transition occurs can be roughly estimated, directly from the initial stellar model, by finding the point where the gravitational binding energy of the envelope and its internal energy (integrated from the surface to that point) balances the change in the NS’s orbital energy.

This implies that for every erg released by the decay of the NS orbit, an additional 1.7 ergs are tapped from the thermal energy content of the envelope. According to the virial theorem, the thermal energy of the envelope relates to its gravitational energy as $E_{thermal}=-[1/(3\gamma-3)]E_{grav}$, where $\gamma$ is the adiabatic index of the envelope that is associated to the equation of state. During the inspiral, if we inject a given amount of energy, in our case $\Delta E_{orbit}$, to the total energy of the system ($E_{thermal}+E_{grav}$), then its thermal energy will decrease by $\Delta E_{thermal}=-[1/(4-3\gamma)]\Delta E_{orbit}$ and its gravitational energy will increase by $\Delta E_{grav}=[(3-3\gamma)/(4-3\gamma)]\Delta E_{orbit}$.

For envelopes of massive stars, $\gamma$ lies between the values of 4/3 and 5/3, corresponding to a purely radiative gas and a monoatomic, ideal gas, respectively. In our specific model, the outer 9 M⊙ of the star initially has $\gamma$ ranging from 1.4 to 1.5, which would translate, based on the virial theorem, to $\Delta E_{grav}/\Delta E_{orbit}\simeq 3$, close to the value of 2.7 that we find (Figure 2).

In Figure 3 we show the ratio of the total change, since the beginning of the simulation, of the envelope’s gravitational energy to the NS’s orbital energy ($\Delta E_{grav}/\Delta E_{orbit}$; orange line), as well as the instantaneous ratio of the changes of the two energies as a function of time ($dE_{grav}/dE_{orbit}$; blue line). Supporting our expectation from applying the virial theorem to the unperturbed model of our star, during the dynamical plunge-in and up to $\approx 5\,\rm yr$ of the simulation, when radiative losses are still negligible, Figure 3 shows that both $\Delta E_{grav}/\Delta E_{orbit}$ and $dE_{grav}/dE_{orbit}$ have average values of $\simeq$3. We note that Ohlmann et al. (2016), find a similar ratio of $\Delta E_{grav}/\Delta E_{orbit}\simeq 2.5$ in their 3D hydrodynamic simulation of a CE, albeit using a very different mass range and simulation set-up.

After 6.5 yr, most of the giant star’s envelope has recombined, releasing $0.23\times 10^{48}\,\rm erg$ of energy into the CE. This recombination reduces the envelope’s opacity by 1-2 orders of magnitude, which in turn reduces the energy transport timescale from the NS to the surface of the CE by two orders of magnitude to $\approx 10^{7}\,\rm s$. At the same time, the inspiral timescale slows to $\approx 10^{9}\,\rm s$. Radiative losses become considerable, and the evolution is no longer adiabatic. Figure 2 shows that by the end of the simulation $0.56\times 10^{48}\,\rm erg$ have been lost from the system via radiation.

While the injection of orbital energy still causes a commensurate release of thermal energy, radiative losses restrict that energy from being entirely used to expel the envelope. Furthermore, as the envelope expands and cools, it becomes less radiation-dominated, leading to an increase in its adiabatic index, $\gamma$. The result of both effects is a reduction in the instantaneous $dE_{grav}/dE_{orbit}$ - and therefore a gradual decrease in the time averaged $\Delta E_{grav}/\Delta E_{orb}$ - seen after $\approx$5 years in Figure 5. One should also keep in mind that at these longer timescales, the energy added to the envelope due to accretion onto the NS and nuclear reactions from the hydrogen-burning shell underneath become non-negligible.

At the end of our hydrodynamic simulation a non-negligible hydrogen envelope still remains around the helium core. This remaining hydrogen-rich layer between the NS’s orbit and the helium core, although small in mass, dominates the overall gravitational binding energy of the envelope due to its small radius. Therefore, exactly how much hydrogen remains around the helium core when the binary exits the CE phase greatly affects its final orbital separation (Ivanova et al., 2013). In the following section, we discuss how we extrapolate our simulation and estimate the final orbital separation and the resulting value of the common envelope efficiency parameter, $\alpha_{CE}$.

3.3 Exiting the Common Envelope and Final Outcome

Figure 3 shows that for the last $\sim$10 years of the simulation, $dE_{grav}/dE_{orbit}$ has converged to a constant value of $\approx$1.4. Furthermore, one can see from the top, right panel of Figure 1 that nearly the entire envelope above the orbit of the NS is unbound ($v/v_{esc}>1$). Therefore, the NS’s inspiral rate in mass coordinate of CE during this phase can be interpreted as the rate at which the envelope material becomes unbound. This allows us to calculate the subsequent inspiral of the NS, after our simulations have stopped, using the equation:

[TABLE]

where the subscript $f$ on the $dE_{grav}/dE_{orbit}$ term refers to the value from Figure 3 at the end of the simulation, which we assume remains constant throughout future evolution. Replacing $\Delta E_{grav}$ and $\Delta E_{orb}$ one finds:

[TABLE]

where $r_{c}$ and $m_{c}$ are the radius and mass coordinate of the NS inside the CE, $m_{NS}$ is the mass of the NS and $G$ the gravitational constant. $E_{grav}(m_{c})$ is the gravitational binding energy of the envelope, integrated from the surface down to mass coordinate $m_{c}$, calculated based on structure of the envelope at the last timestep of our hydrodynamic simulation. Finally, $E_{grav}(m_{c},0)$ and $E_{orb,0}$ are the gravitational binding energy of the envelope down to the position of the NS and the orbital energy of the NS at the end of the hydrodynamic simulation. After solving for $r_{c}$, we find:

[TABLE]

The solid orange line in the left panel of Figure 4 shows the inspiral of the NS from our hydrodynamic simulation. In the right panel, we zoom in on the last stages of the CE evolution resolved by our simulation. There, the dashed orange line shows the extrapolation as predicted by Equation 3.

At the end of our hydrodynamic simulation, the NS is inspiraling into the CE at a steady rate of $\sim 0.003\,\rm M{{}_{\odot}}/yr$ and, equivalently, the donor is ejecting its remaining envelope at the same rate. Calculations of adiabatic mass-loss from giant stars (Deloye & Taam, 2010; Ge et al., 2010, 2015) show that when a large part of the envelope has been removed and the surface hydrogen abundance drops below a critical limit, the reaction of the envelope’s radius to mass-loss (i.e. $\zeta_{ad}=\left(d\ln R/d\ln M\right)_{adiabatic}$) changes suddenly and the envelope contracts quickly, leading to the detachment of the binary.

We repeat such a calculation for our initial model of the 12 M⊙ and 550 R*⊙* giant star, rapidly removing mass by hand at four different rates. The cyan, light blue, purple, and magenta lines in Figure 4 show the radius evolution, from right to left, of the giant star as the star’s mass decreases at each of these mass loss rates. All simulations show that after an initial phase of gradual expansion, once the star’s mass becomes small enough (with the exact mass dependent on the specific mass-loss rate chosen), the radius sharply contracts. Given that at the end of the hydrodynamic simulation the inferred envelope mass-loss rate is $\sim 0.003\,\rm M{{}_{\odot}}/yr$, we expect that the remaining envelope will recede and the binary will detach somewhere between the crossings of the orange, dashed line with the mass-radius relations for mass-loss rates of $0.001\,\rm M{{}_{\odot}}/yr$ (magenta line) and $0.005\,\rm M{{}_{\odot}}/yr$ (blue line). Therefore we conclude that the NS will have inspiraled for another $\sim$0.075-0.225 M⊙, to an orbital separation of $\sim$4.5-8.0 R⊙. The duration of the remaining inspiral is expected to be $\sim 25-75\,\rm yr$.

As soon as rapid mass-loss stops, the envelope is expected to re-expand on a thermal timescale to a giant star structure, albeit with a reduced thermal equilibrium radius, since the surface hydrogen abundance has now dropped to $\approx 0.3$. This implies that upon exiting a CE, the remaining hydrogen-rich envelope will overfill its Roche lobe in a semi-detached configuration and proceed to dynamically stable mass transfer on a nuclear timescale, as shown in recent work by Quast et al. (2019). In their calculations, they find that these type of binaries drive mass-transfer rates of $\sim 10^{-5}\,\rm M_{\odot}\,yr^{-1}$, which in our case would translate to a mass-transfer duration of $\sim 10^{5}\,\rm yr$. The super-Eddington mass-transfer rate is expected to lead to a highly non-conservative mass transfer which will further shrink the orbit. Assuming that the envelope will be removed down to a surface hydrogen abundance of $\approx$10% ($\approx$1%) before the remaining core permanently contracts to radii characteristic of naked helium stars ($\approx$1 R⊙), one can use Equation 8 in Sørensen et al. (2017) to calculate the final orbital separation of the binary. Assuming fully non-conservative mass-transfer, this leads to final orbital separations of $\approx$3.5-5.7 R⊙ ($\approx$3.3-5.3 R⊙).

From this last estimate of the final orbital separation ($a_{f}\approx 3.3-5.7\,\rm R_{\odot}$) we calculate a value of $\alpha_{CE}$ that can be directly compared to the results of rapid binary population synthesis codes. In Figure 5 we consider binaries consisting of a 12 M⊙ giant star and a NS that initiate a CE phase at different initial orbital separations. We calculate the final orbital separation of the binary with the traditional $\alpha_{CE}$ prescription (where $\alpha_{CE}$ is defined as the ratio between the envelope’s gravitational binding energy and the change in the binary orbital energy), using different values of $\alpha_{CE}$ and two different definitions of the core-envelope boundary (1% and 10% hydrogen fraction). Notably, the estimated final separation from our hydrodynamic simulation and the arguments we presented earlier, correspond to an $\alpha_{CE}$ value of $\sim 5$, pointing to an extremely efficient CE ejection.

4 Conclusions

Motivated by the challenges faced by three dimensional simulations, we model the evolution of a high-mass binary system through a common envelope using the one dimensional hydrodynamic stellar evolution code MESA. We are able to evolve a 12 ${M_{\odot}}$ donor star, as a 1.4 ${M_{\odot}}$ neutron star spirals into it, following the evolution through the three phases of CE evolution: “loss of co-rotation”, “plunge-in”, and “self-regulated inspiral”. While the initial phase of a CE may occur on a dynamical timescale, we demonstrate that the self-regulated regime occurs on a thermal timescale, where radiative transport becomes essential.

By the end of our simulation, the majority of the envelope ($\simeq 8{M_{\odot}}$) has expanded to radii larger than 10^{4}$${R_{\odot}}, and is formally unbound from the binary. At the same time, the NS has inspiraled from an initial orbital separation of 550 ${R_{\odot}}$ down to $\simeq$10 ${R_{\odot}}$. By analyzing the relevant energies in the system, we find the envelope is predominantly expelled by tapping into the original thermal energy of the envelope.

We estimate the final state of the binary by extrapolating the inspiral until the point where the donor star contracts. Subsequent re-expansion will cause the remaining envelope to be removed by stable mass-transfer. Such mass transfer could drive super-Eddington accretion onto a NS, potentially forming a NS ultra-luminous X-ray binary that would last for $\gtrsim 10^{5}$ yr. Although the binary would still be enshrouded in an optically thick envelope at a radius of $\sim$100 AU, such a phase may still be visible in X-rays. The recently observed NS ultraluminous X-ray source NGC 7793 P13 (Fürst et al., 2016; Israel et al., 2017) is potentially such a post-CE X-ray binary (Quast et al., 2019).

Using a quantitative prescription for the final inspiral phase and the potential subsequent stable mass transfer, we find that the final binary consists of a $\simeq 2.6\,\rm M_{\odot}$ helium star orbiting a NS at an orbital separation of $\simeq$3.3-5.7 ${R_{\odot}}$. According to detailed calculations by Ivanova et al. (2003) and Tauris et al. (2015), this binary configuration will lead to the formation of a binary NS that will merge within a Hubble time. Most interestingly, the final post CE separation translates to a very efficient CE ejection: $\alpha_{CE}\sim 5$ for the specific initial binary configuration. Our results suggest a higher efficiency than most previous studies since our simulations find that the most bound hydrogen layers surrounding the helium core are removed non-conservatively, after the binary has detached. Although there is no reason to believe that this estimated effective $\alpha_{CE}$ should have a universal value, it is tempting to think about the implications of such an efficient CE ejection. For example, it could resolve inconsistencies between rates predicted by the CE channels for coalescing binary black holes and binary NSs (e.g. Mapelli & Giacobbo, 2018) as well as the formation of black hole low-mass X-ray binaries (e.g. Podsiadlowski et al., 2003).

We thank Bill Paxton, Pablo Marchant, Fred Rasio, and Matteo Cantiello for useful conversations. The authors acknowledge support from the Swiss National Science Foundation grants (project numbers PP00P2 17686 and 200020-172505), the Danish National Research Foundation (DNRF132), the Marie Sklodowska-Curie RISE grant “ASTROSTAT” (project number 691164), the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 617001, the National Science Foundation grant AST-1517753 and a CIFAR grant from the Gravity and Extreme Universe Program and a Simons Foundation grant. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.