© The Authors 2024

1. Introduction

Our understanding of the progenitor systems leading to gravitational-wave (GW) emitting mergers of double compact objects is challenged by several major knowledge gaps (e.g., Dominik et al. 2012; Ivanova et al. 2013; Belczynski et al. 2016; Tauris et al. 2017; Vigna-Gómez et al. 2018; Giacobbo & Mapelli 2018). Two of which are the reaction of the system to the supernova (SN) explosion forming the compact objects and the ability of a common-envelope (CE) phase to shrink the orbit sufficiently for GW emission, initiating a merger within a Hubble time.

In the CE interaction, the more compact secondary object is engulfed by the envelope of the giant primary star due to a preceding dynamically unstable mass transfer (MT, e.g., Eggleton 2011) phase or a tidal interaction such as the Darwin instability (Darwin 1879), for example. When orbiting inside the CE, the companion and the core of the primary star experience drag forces caused by the gravitational pull of accumulated envelope material behind them through self-gravity (Chandrasekhar 1943; Dokuchaev 1964; Ostriker 1999; Kim & Kim 2007; Kim 2010). These drag forces tighten the orbit of the binary system that formed from the core of the primary star and the companion (referred to as “core binary system” in the following). During this process, the core binary system transfers angular momentum and orbital energy to the envelope material, which eventually leads to a partial or full ejection of the CE.

The motion of the companion object and the core of the primary star within the envelope material triggers numerous (magneto-)hydrodynamical instabilities that not only affect the final orbital separation, but also shape the flows of the magnetized plasma. The implied dynamics create unique features such as magnetic outflows perpendicular to the orbital plane (Ondratschek et al. 2022) and the potential formation of a centrifugally supported structure interacting with the inner post-CE binary during and after the dynamic CE interaction (Röpke & De Marco 2023; Gagnier & Pejcha 2023; Tuna & Metzger 2023; Wei et al. 2024).

Given the vast range of spatial and temporal scales involved in the problem – spanned by the scales of the progenitor star and the compact companion – three-dimensional (3D) simulations of massive CE systems face significant difficulties in following the CE phase to full envelope ejection (e.g., Law-Smith et al. 2020; Moreno et al. 2022; Lau et al. 2022a,b). During the simulated time span, these models fail to tighten the central binary to distances at which the emission of GWs becomes important, as predicted by the analytical prescription (e.g., Webbink 1984) used in population synthesis models (e.g., Hurley et al. 2002; Dominik et al. 2012; Stevenson et al. 2017; Vigna-Gómez et al. 2018; Mapelli & Giacobbo 2018; Kruckow et al. 2018; Belczynski et al. 2016). This strong discrepancy between theoretical expectations and numerical findings has to be settled in order to make reliable rate predictions of GW-induced mergers of compact objects. Following the discussion of Moreno et al. (2022), their simulations appear to be converged in the final orbital separations. In this case, two causes for the discrepancy are conceivable: First, the current simulations fail to capture all relevant physical effects in the plunge-in phase of the CE interaction, and second, there are mechanisms after the rapid spiral-in sufficiently contracting the core binary until the system reaches full envelope ejection.

Here, we investigate the second case. We report on the extensions of the simulations conducted by Moreno et al. (2022) of two CE interactions of a 10 M_⊙ post core-helium burning red super-giant (RSG) primary star (at an age of 25.5 Myr) with a black hole (BH) and a neutron star (NS) companion. After a SN explosion of the core of the primary star (most likely forming a NS, Podsiadlowski et al. 2004; Schneider et al. 2021), the remaining compact binary system would be a promising progenitor for a NS-BH or NS-NS merger event (e.g., Podsiadlowski et al. 2004; Vigna-Gómez et al. 2018; Giacobbo & Mapelli 2018; Belczynski et al. 2016; Mandel & Broekgaarden 2022). However, Moreno et al. (2022) found a final orbital separation of a_f,  NS = 16 R_⊙ for the NS companion and a_f,  BH = 47 R_⊙ for the BH companion at the end of their simulation, which are both too wide for the systems to merge within a Hubble time by GW emission. They further showed that even taking into account a subsequent second MT phase as well as an isotropically distributed natal NS kick, only 8.7% of the NS companion systems and 0.6% of the BH companion systems are expected to form a NS-NS or NS-BH merger, respectively. We mitigated the tremendous computational cost encountered in the simulations conducted by Moreno et al. (2022) with a new refinement approach of gas cells deep in the gravitational potential of the primary core and companion.

We summarize our model and methods in Sect. 2 and show therein the pivotal aspects of our initial primary star (Sect. 2.1), CE model (Sect. 2.2), and introduce the new refinement approach (Sect. 2.3). The results of our work are shown in Sect. 3, where we describe the temporal evolution of the orbital separations and the envelope ejection (Sect. 3.1), the appearance of a geometrically thick toroidal structure and magnetically driven bipolar outflows (Sect. 3.2, where the latter is described in more detail in the upcoming publication Vetter et al. in prep.), the identification of the torus as a circumbinary disk (CBD) and its characteristic properties (Sect. 3.3), updated α_CE values (Sect. 3.4), and possible final fates of the systems (Sect. 3.5) based on the outcome of our simulations. Last, we discuss our results in Sect. 4 and conclude in Sect. 5.

2. Methods

We use the identical CE models as presented by Moreno et al. (2022), where they closely adhere to the procedures outlined in Ohlmann et al. (2016a), Ohlmann et al. (2017), and Sand et al. (2020) for conducting CE simulations. We refer to these publications for a more detailed discussion and additional information, and only summarize basic aspects of both the initial primary star and the CE models in Sect. 2.1 and Sect. 2.2.

The three-dimensional simulations shown in this work are conducted with the 3D moving-mesh magnetohydrodynamics code AREPO (Springel 2010; Pakmor et al. 2011; Pakmor & Springel 2013), which employs a second-order finite-volume method. The Powell scheme (Powell et al. 1999; Pakmor & Springel 2013) is utilized to control the divergence of the magnetic field. Newtonian self-gravity was calculated with a tree-based algorithm. To account for recombination energy, we applied the OPAL (Iglesias & Rogers 1996; Rogers et al. 1996; Rogers & Nayfonov 2002) equation of state (EoS) similar to previous CE simulations conducted by, for example, Sand et al. (2020), Kramer et al. (2020), and Moreno et al. (2022).

Here, we introduce a new refinement approach (Sect. 2.3) for the cells within the softened potential around the primary core and companion (Ohlmann et al. 2016a) to ease the computational costs encountered in the late stages of the simulations presented by Moreno et al. (2022), which prevented a continuation to later times in the evolution of the considered systems.

2.1. Establishing a stable primary star model

For our initial primary star model, we consider a M₁ = 10 M_⊙ zero age main sequence star with solar metallicity (Z = 0.02). The star is first evolved with the one dimensional MESA code (Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019) to a RSG with an age of about 25.5Myr. At the desired evolutionary stage, the star has exhausted its core-helium abundance and already developed a deep convective envelope (Moreno et al. 2022, Fig. 1). The considered stellar model has a mass of M₁ = 9.4 M_⊙, a radius of R₁ = 395 R_⊙, an effective temperature of T_eff = 3564K and a logarithmic luminosity of logL/L_⊙ = 4.35.

Subsequently, the one-dimensional stellar profile is mapped onto the 3D grid of AREPO closely following the procedure outlined in Ohlmann et al. (2016a). During this process, two main challenges arise: First, we have to represent the inner core region, which is very small compared with the overall spatial scales of the RSG. As discussed by Moreno et al. (2022), resolving this core also implies a severe timescale problem. Second, the unavoidable introduction of discretization errors perturbs the star from its hydrostatic equilibrium (HSE) and gives rise to spurious velocities, which have to be kept under control.

The first challenge is addressed by replacing the star’s core with a point mass that interacts with the other material only via gravity (Ohlmann et al. 2017, Section 3.2). In this context, the core of the star is defined based on the cut radius R_cut, which is ideally chosen such that (i) most of the convective envelope is resolved on the AREPO grid (the bottom of the convective envelope is at r = 17.9 R_⊙ with mass coordinate m = 3.07 M_⊙), (ii) it encompasses key points, such as the compression point (Ivanova 2011; r = 0.40 R_⊙, m = 2.78 M_⊙) and the location where the hydrogen mass fraction X drops below 0.1 (Dewi & Tauris 2000) (r = 0.33 R_⊙, m = 2.76 M_⊙) and ideally (iii) below or near the final orbital separation to ensure convergence (Ohlmann et al. 2017; Moreno et al. 2022). The points in (ii) collectively define the amount of envelope material which needs to be ejected to prevent immediate re-expansion of the star.

While a cut radius on the order of criterion (ii) is beyond the capabilities of current numerical approaches, conditions (i) and (iii) can be met, and Moreno et al. (2022) find a hydrostatically stable primary star and converging results in the final orbital separations in their relaxation and CE simulations, respectively. The considered primary model was set up with R_cut ≈ 20 R_⊙ and a (linear) resolution of the softened sphere of N_cps = 40 grid cells (“cps” stand for “cells per softening length”, i.e., the actual spatial resolution of the softened volume is approximately given by N_cps³). With this choice, 97% of the mass of the hydrogen-rich envelope is contained on the hydrodynamic grid. The gas cells within the core region are adapted to the solution of a modified Lane-Emden equation to ensure a smooth transition of the gradients over the cut radius (Ohlmann et al. 2017).

The simulated 3D domain extends far beyond the radius of the star and spans 3159 R_⊙ in each spatial dimension. The star is embedded in a “pseudo-vacuum” with a density of 10⁻¹⁵g cm⁻³ and a temperature of 4000 K. We utilize a Lagrangian refinement criterion for the gas cells outside the softened spheres of point particles, that is, the mass of all cells (on average, our simulation contains 4 569 947 grid cells) fluctuates by less than a factor of two around a target mass resolution of m_cell = 1.6 × 10⁻⁶ M_⊙. It is important to note that during the CE simulation, the simulated box expands far beyond the initially set box size, and the number of cells contained in our domain also varies given the applied mass refinement criterion.

2.2. Common-envelope phase

We aim to simulate two CE interactions with secondary masses of M₂ = 5 M_⊙ and M₂ = 1.4 M_⊙, implemented as gravity-only point particles (similar to the core of the primary star). The companion masses are chosen to represent a BH and a NS companion, respectively. With the implementation of the companions as point particles, we neglect potential accretion, neutrino cooling and jet formation through accretion streams onto the compact objects in our simulations. The companions are initially in full corotation with the envelope material¹.

A seed dipolar magnetic field with a surface field strength of 10⁻⁶ G is applied to the primary star in all conducted simulations (Ohlmann et al. 2016b; Moreno et al. 2022). To quantify the amount of ejected (i.e., unbound) envelope material during the CE event, we adhere to the commonly used energy criteria (e.g., Ohlmann et al. 2016a; Sand et al. 2020; Moreno et al. 2022; Lau et al. 2022a,b; Chamandy et al. 2024), which read,

$$\begin{array}{cc}\hfill e_{\mathrm{kin}}+e_{\mathrm{pot}}>0,& \text{kinetic-energy criterion},\hfill \end{array}$$(1) $$\begin{array}{cc}\hfill e_{\mathrm{kin}}+e_{\mathrm{therm}}+e_{\mathrm{pot}}>0,& \text{thermal-energy criterion},\hfill \end{array}$$(2) $$\begin{array}{cc}\hfill e_{\mathrm{kin}}+e_{\mathrm{int}}+e_{\mathrm{pot}}>0,& \text{internal-energy criterion},\hfill \end{array}$$(3)

where e_kin is the kinetic, e_therm is the thermal, e_int is the internal energy (including thermal and ionization energy), and e_pot is the potential energy, respectively. We express the amount of unbound material by the sum over all cell masses, which satisfy the individual criteria, and relative to the total envelope mass f_ej = M_unbound/M_env (here, M_env  ≈  6.648 M_⊙ is defined by the total gas mass contained on the AREPO grid in our simulations).

In contrast to Moreno et al. (2022), we additionally employ Lagrangian tracer particles in our simulations to post-process the gas flows, since we are still restricted by the finite-volume approach of AREPO. These particles are implemented in AREPO as passively advected Lagrangian particles. The trajectories are advected with the local velocity field interpolated to the position of the tracer particles (Genel et al. 2013). In all simulations, we add a total of N_tr = 799 212 tracer particles (i.e., approximately a sixth of the initial number of hydro-cells) that sample the initial mass distribution over all hydro-cells in the entire simulation box. Given the target mass resolution in our simulations, each tracer particle represents the same mass of m_tr = M_env/N_tr  ≈  8.32 × 10⁻⁶ M_⊙.

2.3. New refinement criterion for cells inside the softening length of the point particles

As already briefly mentioned in Sect. 2.1, the primary-core approximation (Ohlmann et al. 2016a) requires a sufficient number of cells per softening length N_CPS for the star to establish HSE in the relaxation run. However, during the spiral-in of the companion, the softening length h is reduced proportional to the orbital separation a between the point particles as soon as the separation decreases to a(t) = 2.5h(t) (enforcing h(t)  ≤  0.4a(t) for any given time)². This behavior ensures that the softened spheres of the core and secondary do not overlap. In previous works (e.g., Ohlmann et al. 2016a,b, 2017; Sand et al. 2020; Moreno et al. 2022), the number of cells per softening length was kept constant during the CE interaction. This decreases the cell sizes inside the softened gravitational potential, leading to an increased physical resolution within the approximated core region. According to the Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928), the gradual decline in size of the cells reduces the time step in the binary simulation.

This problem becomes even more evident in the final stages of the CE interaction, when the core is in a close orbit with the companion and detached from its former envelope material. The necessity of the forced decrease in cell size caused by the constant resolution N_CPS = const criterion seems to be questionable, since we neither have to ensure HSE anymore nor are we able to resolve the core in the first place. As a result of the decreased cell sizes, the simulations of Moreno et al. (2022) became prohibitively expensive after the spiral-in due to too small time steps and only f_ej,  BH ≈ 48% (f_ej,  NS ≈ 54%) of the envelope mass in the CE simulation with the BH (NS) companion was unbound according to the kinetic energy criterion (Eq. (1)) over the affordable simulation time.

To ease the computational costs in our simulations, we thus follow the idea of controlling the number of cells per softening length N_CPS by establishing a new refinement criterion, for which the number of cells per softening length is a power law of h,

$$\begin{array}{c}\hfill N_{\mathrm{CPS}}\propto h^{\kappa }\Rightarrow N_{\mathrm{CPS}}(t)=N_{}\mathrm{CPS},0{\left(\frac{h(t)}{h_0}\right)}^{\kappa },\end{array}$$(4)

where N_CPS,  0 is the number of initially used cells per softening length and h₀ the initial softening length, which we set to the cut radius R_cut of our relaxed primary star model in Sect. 2.1. We note that κ = 0 represents the standard refinement criterion with N_CPS = N_CPS,  0  =  const used in previous works, while we use κ = 1 in this work.

Ultimately, κ = 1 leads to an approximately constant cell size inside the softening length throughout the CE simulation and thus increases the lowest time step compared to κ = 0 (see Fig. A.1 and Appendix B in general for more information). Neither κ = 0 nor κ = 1 are ideal solutions and are born out of the necessity to approximate the primary core and companion and in future studies, a more sophisticated model is desirable.

3. Results

3.1. Global orbital evolution and envelope ejection

Once the companions are placed above the surface of the primary star with initial separations of a_i,  BH  =  501 R_⊙ and a_i,  NS  =  479 R_⊙ for the BH and NS companion, respectively, we observe a similar dynamic behavior during the plunge-in and rapid spiral-in phases as Moreno et al. (2022). Consequently, we find a good overall agreement between the simulations with κ  =  0 and the simulations with the new criterion with κ  =  1 in the orbital evolution (almost perfect correlation with only small deviations, see Appendix B, Fig. B.1) as well as in the fraction of ejected envelope material f_ej as shown in Fig. 1. However, in the late-time evolution, both systems (especially the NS scenario) develop small changes in eccentricity (and even orbital precession in the system involving a NS companion, that is, Fig. 2b). At present, there is no clear explanation for this behavior, and we refer to Sect. 4.1 and Appendix B for further discussion.

Fig. 1. Temporal evolution of the orbital separation a (solid dark) and the fraction of unbound material fej according to the kinetic (dashed green, Eq. (1)), thermal (dashed rust, Eq. (2)), and internal energy (orange, Eq. (3)) criterion, respectively. The simulation involving the BH is shown on the left, while the NS scenario is on the right. The thin lines in the plots represent the simulations that contained the new refinement approach (see Sect. 2.3, κ  =  1), while the thick, faint lines correspond to the results of Moreno et al. (2022). On the left-hand side, we additionally show the unbound mass fraction only for cells contained within a radius of R  =  1.6 × 104 R⊙ around the central binary due to the artificial “rebinding” (see Sect. 3.1) of material in the pseudo-vacuum for t  ≥  4400 d (see Sect. 3.1). Furthermore, the hourglass markers represent the individual starting points of the magnetically driven outflow in the simulations with κ  =  1. The minor ticks in both plots in between 1 − 10 orbits are adapted for readability, and only the 5 orbit mark is shown.

Fig. 2. Slices of spherical radial velocity vr, sph. in the x − z plane in (a)-(f) and in the x − y plane in (g)-(j) for the BH (NS) companion in the left (right) column. The iso-surfaces in (a) and (b) enclose the bound mass according to the kin. energy criterion (Eq. (1)) in dashed gray. The gray streamlines in (c)-(j) represent the velocity field in the corresponding plane. The primary core (companion) is marked with a black cross (plus sign).

We were able to follow the dynamical CE phase to almost full envelope ejection (more than 97% of the former envelope in both systems is formally ejected according to Eq. (1)) which was established after ≈13 708d for the BH companion and ≈7 330d for the NS companion, respectively (Fig. 1). In the BH companion case, we observe an increasing portion of material that becomes bound again in the pseudo-vacuum (in the following referred to as “rebinding of material”) at R >  3 × 10⁵ R_⊙. There, the already formally ejected envelope material is forced to dissipate its kinetic energy to the ambient low-density material and, hence, appears to be bound again according to the kinetic energy criterion (Fig. 1a, dashed green line). Notably, the immediate fallback of this material is prevented by a pressure support established through the later on ejected material. A similar observation has also been made in the simulations conducted by Chamandy et al. (2019, 2024). Here, we consider it an artifact caused by the need to fill the outer regions of our simulation domain with low-density pseudo-vacuum. We correct for this unexpected rebinding effect in the fraction of unbound envelope material f_ej by considering cells with distances larger or equal to 1.6 × 10⁴ R_⊙ with respect to the central core binary system as unbound independently of the result Eq. (1) (Fig. 1a, dash-dotted green line).

Upon almost full envelope ejection, the conducted CE simulations result in a final orbital separation of a_f,  BH  ≈  37 R_⊙ (after 1220 orbits) for the BH scenario, and a_f,  NS  ≈  14 R_⊙ for the NS companion establishes (after 1831 orbits). Compared with the final distances reported by Moreno et al. (2022) of a_f,  BH  ≈  47 R_⊙ and a_f,  NS  ≈  15 R_⊙, the core binary hardened by ∼21% in the BH companion case and its distance changes only marginally for the CE interaction involving the NS companion.

3.2. Reshaping of the envelope

Similar to the low-mass system studied by Ondratschek et al. (2022), the seed magnetic field rapidly amplifies during the pre-plunge and in-spiral phases of the CE interaction (Ohlmann et al. 2016b) and leads to the formation of a magnetically driven bipolar outflow (see Fig. 1, black marker) perpendicular to the orbital x − y plane. This feature is best visible in the spherical radial velocity v_r,  sph. (the subscript “sph” stands for the calculation in spherical coordinates) in Fig. 2a and 2b for the BH and NS companion, respectively. The accelerated material reaches peak radial velocities of almost ∼300 km s⁻¹. Moreover, the rapidly expelled material forms characteristic bow shocks (Fig. 2a and 2b) resulting from the interaction with the already ejected envelope (similar to the findings of Ondratschek et al. 2022).³

Prior to the onset of the bipolar outflow, the bound material in the simulations (approximately 68% [BH companion] and 76% [NS companion] of the total envelope mass according to Eq. (1)) redistributes to a thick toroidal structure in the orbital plane of the core binary (see gray dashed contour in Fig. 2a and b) and together with the bipolar outflows show morphological similarities to bipolar pre-planetary nebulae (PNe) observed in lower-mass systems.

Furthermore, radially in-falling regions onto the central binary develop (see blue areas throughout the entire Fig. 2), suggesting the formation of a thick circumbinary accretion disk (Sect. 3.3) refueling the magnetically driven outflows as already proposed but not confirmed in Ondratschek et al. (2022). The velocity fields in the orbital plane exhibit circularized motion (Fig. 2i and j) and even show indications for the presence of turbulence (Fig. 2e and f), possibly enabling angular momentum and mass transport along the orbital x − y plane.

The circular motion of the bound material transitions into a radially outward pointing motion of the gas for radii larger than R  ≈  1000 − 4000 R_⊙ with respect to the core binary system (e.g., in Fig. 2g and 2h, the transition is around R  ≈  2000 R_⊙ for the BH companion and R  ≈  1000 R_⊙ for the NS companion). This transition region appears to coincide with that of hydrogen recombination (see Movie M2 and Movie M4 in Table C.1 in Appendix C) and we find that the ejection of this material is mainly driven by heating from recombination of helium and hydrogen in the cool outer parts of the disk (as we subsequently explain in more detail in Sect. 3.3) and resembles a radial wind with radial velocities of ≈40 − 50km s⁻¹. In fact, the envelope ejection at that stage is completely dominated by the radially outflowing material, with progressively decreasing mass loss rates on the order of Ṁ ∝ 10⁻⁴ M_⊙ d⁻¹ down to 10⁻⁵ M_⊙ d⁻¹ as the simulations evolve toward full envelope ejection. The point in time where the radial winds become dominant in the CE ejection is visible by the steeply increasing fraction of ejected material f_ej according to Eq. (1) in Fig. 1a (1b) around ∼1400 d (∼700 d) for the BH (NS) case. The reshaped envelope is sufficiently adiabatically cooled at this point that recombination commences.

In summary, we observe a morphological transformation of the envelope due to orbital energy and angular momentum transferred during the spiral-in. The initial spherical primary star is reshaped into a toroidal structure surrounding the core binary (for example, see Movie M1-M4 in Table C.1). The foundation for a CBD might thus have already been set as early as the end of the plunge-in phase (i.e., within the first ≈800 d for the BH and ≈500 d for the NS companion) and we further describe this in Sect. 3.3. The CBD is affected by two primary mass-loss channels. First, the inward pointing mass flux indicated by the negative radial velocities of the gas and second, the radial winds driven by recombination. The first channel is responsible for refueling the magnetically propelled jet-like outflows perpendicular to the orbital plane, which, together with the bound torus, resemble the morphological structure of bipolar pre-PN.

3.3. Formation of a circumbinary disk and recombination winds

During the rapid in-spiral of the companion, the envelope of the primary star is lifted and transformed from its initially spherical shape into a torus-like structure. We further find the material to be rotationally supported (i.e., circularized, and Keplerian motion) and still to be gravitationally bound. Hence, the basic requirements for the remaining envelope material to be a disk are already fulfilled (e.g., Pringle 1981) and in the following we want to characterize the system in more detail. The main open questions to address in this regard are thus: is the specific angular momentum distribution such that the system can maintain a quasi-steady state, does the material develop HSE in the vertical direction, and does the system develop a (partially) evacuated cavity around the core binary.

Given the presence of the magnetized outflows in our simulations (Sect. 3.1), the ejection of material perpendicular to the orbital plane (Fig. 3a) can act like a “pressure valve” (Chamandy et al. 2018) and prevent the virialization of the gas in the bound material around the core-binary. Therefore, the system cannot provide pressure support and, thus, disable radially inward mass migration. This already has become evident in the form of radially inward propagating material in the vicinity of the binary (Sect. 3.2) and we later follow up on this idea and measure the mass flux toward the central binary.

Fig. 3. Specific entropy s and density ρ slices of the simulation with the BH companion at t  =  8244 d. The quantities are shown to the left (x ≤ 0 for the spec. entropy) and right (x ≥ 0 for the density) in an edge-on view (x − z) in (b), respectively. Furthermore, the distribution of all cells along the z-direction in a cylindrical shell (indicated by the dashed rectangle) in (a) and (c). The gray streamlines in (b) show the velocity field in the x − z plane.

Regarding the development of a quasi-steady state and vertical HSE, we face significant difficulties with the confirmation of the establishment of a CBD in our complex 3D simulations. Especially in the simulation involving the NS companion, the orbital precession further complicates matters (Sect. 3.1) and we restrict our analysis to the simulation with the BH companion in the following. There, the density as well as the specific entropy in the toroidal structure (Fig. 3b, or Movie M1 in Table C.1) follow a clear cylindrical symmetry (except for local fluctuations). The material is perturbed in both quantities by the spiral shocks generated through gravitational torque exerted from the core-binary onto the gas in its immediate vicinity. The vertical (z-direction) distribution of the density and specific entropy point (Fig. 3a and c) toward a convectivly stable, Gaussian-like profile as expected for disks in vertical HSE (e.g., Pringle & King 2007).

To check the temporal persistence of the disk structure, we show the time evolution of the radial distribution of quantities evaluated on cylindrical⁴ shells in Fig. 4. For all taken averages within a cylindrical shell, we use mass-weighted averages (i.e., ⟨Q⟩  =  ∫Qρ dV/∫ρ dV for any Q). Furthermore, we consider 2 × a(t) as the inner border (e.g., Muñoz et al. 2020; Siwek et al. 2023b,a), while the outer boundary for our averaging is kept constant at R_out  =  4000 R_⊙. The choice for the outer radius is motivated by substantial loss of material radially outward in the orbital plane at around R_out  ≳  4000 R_⊙ (see streamlines in Fig. 2g, h), indicating the presence of radial winds (Fig. 4h) that further strip the envelope material as briefly described in Sect. 3.2. In the first four panels of Fig. 4, the temporal evolution of the average density ⟨ρ⟩ (Fig. 4a), the cumulative mass M_cum (Fig. 4b), the square of the specific angular momentum J² (Fig. 4c) and the average angular velocity ⟨Ω⟩ are shown. Unsurprisingly, the density profile in Fig. 4a decreases over time as the bound material is progressively ejected (Fig. 4b). The squared specific angular momentum (Fig. 4c) strictly increases along the radial direction, hence, pointing toward a Rayleigh-stable configuration (∂J²/∂r >  0, e.g., Pringle & King 2007) and the angular velocity (Fig. 4d) reflects the behavior of a pressure-supported sub-Keplerian disk (e.g., Adachi et al. 1976; Weidenschilling 1977, we find Ω ∝ r^−1.86). Hence, the system obeys the dynamic properties of a CBD, but is experiencing mass and angular momentum loss over time.

Fig. 4. Time evolution of averaged density ⟨ρ⟩ in (a), cumulative mass Mcum. in (b), specific angular momentum J2 in (c), angular velocity ⟨Ω⟩ in (d), isothermal pressure scale height H in (e), the aspect ratio H/r in (f), the temperature T in (g), the column density Σ in (h), the ionization fraction of HII (χHII), HeII (χHeII), and HIII (χHeIII) in (i), (j), and (k). Last, the average optical depth of cells τcell is shown in (l). The averages are taken over the range of r  ∈  [2a(t), 4 × 103 R⊙] (Sect. 3.3).

We further try to explore characteristic properties of a CBD in panels (e) to (h) of Fig. 4, such as the isothermal pressure scale height H, the aspect ratio θ ≡ H/r, the average temperature ⟨T⟩ and the vertically integrated column density Σ  =  ∫ρ dz. The pressure scale height in Fig. 4e is measured by fitting the density distribution along the vertical direction with a Gaussian profile (e.g., Fig. 3c). We find that H remains approximately constant throughout the simulation. The vertical scale height is resolved with approximately 10(3) to 60(28) cells per H at t  =  3555 d (t  =  13 698 d) for the inner and outer boundary, respectively. The aspect ratio (Fig. 3f) is fluctuating around ≈0.8 and thus confirming the CBD to be thick. Similar to the density, the temperature (Fig. 4g) and column density (Fig. 4h) decrease with time.

With these quantities at hand, we can now further measure the mass transport rates through the different ejection channels. Given the cylindrical symmetry of the system, the mass flux can be approximated by the average transport of cylindrical shells $\dot{M}=2\pi r{\overline{v}}_r\mathrm{\Sigma }$ with ${\overline{v}}_r$ being the average cylindrical radial velocity. With this definition, negative (positive) values correspond to radial inward (outward) propagation of the corresponding shell. We note that with this definition, the outward mass transport rates may be underestimated because the recombination energy is locally thermalized, and the resulting velocity field tends toward spherical symmetry, depending on the local pressure gradient (see Fig. 2c and d). The resulting radial profile of the transport rates is shown in Fig. 5. The mass flux is dominantly inward for radii r ≲ 10³ R_⊙ and points outward for lager radii, reflecting our findings from Sect. 3.2, where we concluded that the bound torus is eroded by radial winds and accretion toward the core binary, which, before reaching it, is channeled into a jet-like outflow.

Fig. 5. Radial distribution of mass transport rates $\dot{M}=2\pi r{\overline{v}}_r\mathrm{\Sigma }$ similar to Fig. 4. The faint bands over the entire radial domain in the background correspond to the overall envelope ejection rate Ṁej of the system at the corresponding time with the matching color.

The negative mass transport rates in the inner part of the disk (out to about 600 R_⊙ in Fig. 5, which we will in the following refer to as inflow rates) appear to be relatively constant throughout the simulation. They are on the order of Ṁ ∼ −10⁻⁵ M_⊙ d⁻¹ and decrease significantly only for very late times (for t = 13 698 d, where we measure Ṁ ∼ −10⁻⁶  M_⊙ d⁻¹, see Fig. 5). The notion of inflow rates over the disk is supported by computing the α-viscosity given by the mass weighted average (Balbus & Hawley 1998),

$$\begin{array}{c}\hfill \alpha =\langle \frac{\delta v_r\delta v_{\varphi }}{c_s^2}-\frac{B_r{B}_{\varphi }}{4\pi \rho c_s^2}\rangle ,\end{array}$$(5)

where δv_r  =  v_r and $\delta v_{\varphi }=v_{\varphi }-\sqrt{G(M_{\mathrm{core}}+M_2)/r}$ are the radial and azimuthal velocity perturbations, c_s is the sound speed, B_r and B_ϕ are the radial and azimuthal magnetic field components. This yields average α-viscosity values in the inner (i.e., well within the hydrogen recombination front; R  ≲  600 R_⊙) parts of the bound material of ≈0.06 to 0.13. This is in good agreement with the expectations from local magneto-rotational instability (MRI) simulations with α-values on the order of 10⁻² − 10⁻¹ (e.g., Davis et al. 2010). Furthermore, we can reproduce the inflow rates measured in our simulation by an α-disc model (Shakura & Sunyaev 1973),

$$\begin{array}{cc}\hfill {\dot{M}}_{\mathrm{visc}}& =-3\pi \alpha r^2\mathrm{\Omega }{\theta }^2\mathrm{\Sigma }\approx -({10}^{-4}-{10}^{-6})M_{\odot }{\mathrm{d}}^{-1},\hfill \end{array}$$(6)

for α  ≈  10⁻¹ at the inner disk edge (i.e., r  ≈  80 R_⊙) with the local angular velocity Ω  ≈  10⁻⁶ s⁻¹, the aspect ratio θ  ≈  0.8 (Fig. 4f) and the integrated column density Σ  ≈  8 × 10⁴ g cm⁻² (Σ  ≈  8 × 10² g cm⁻²) for t  =  3555 d (t  =  13 698 d). Most (if not all) of this inward-transported material, however, is found to be launched into the magnetically driven outflow mentioned in Sect. 3.1 (Vetter et al. in prep.). It will not participate in altering the orbital separation via angular momentum accretion onto the engulfed core binary system (e.g., Gagnier & Pejcha 2023; Wei et al. 2024; Siwek et al. 2023b).

The mass transport in the outer regions (i.e., beyond the hydrogen recombination front at maximal radii in Fig. 5) exceeds the inflow rates in the inner parts by approximately one order of magnitude (Ṁ ∼ 10⁻⁴−10⁻⁵ M_⊙ d⁻¹) and thus is expected to dominate the ejection rates in the entire system. In fact, we find good agreement in the overall mass ejection rates in the system Ṁ_ej (Fig. 5).

The outer transport region coincides with the recombination front of hydrogen and helium (see Fig. 4i–k, where we show the ionization fraction of H_II, He_II and He_III). Assuming full conversion of released recombination energy to kinetic energy, the expected velocities of the ejecta escaping the gravitational pull of the core binary (∼47 km s⁻¹ for helium and hydrogen and ∼38 km s⁻¹ for only hydrogen recombination)⁵ approximately match the characteristic outflow velocities observed in our simulation of ∼42 km s⁻¹ (Fig. 2g). This supports the argument for a recombination-driven radial wind. As the material is progressively ejected, the remaining bound material cools adiabatically (Fig. 4g), leading to inwardly propagating recombination fronts (Fig. 4i, j, k) and finally to the ejection of the material according to Eq. (1).

Given the sharp decline in optical depth of the cells τ_cell  =  2κ_opacρR_cell (with κ_opac being the opacity) closely behind the hydrogen recombination front (see Fig. 4l), the ejection of the bound material in the radial winds could be influenced by cooling. The adiabatic assumption in our simulations (i.e., local thermalization of the released recombination energy) could thus overestimate the amount of ejected material in the radial winds through hydrogen recombination. In this case, more material would be expected to remain bound in the CBD and further evolve toward a thin disk structure.

3.4. Common-envelope ejection efficiency

In the following, we are trying to compare the expectations from the α-formalism to the final orbital separations found in our simulations. Applying the alpha-energy formalism (Webbink 1984),

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{CE}}=\frac{E_{\mathrm{bind}}}{\mathrm{\Delta }{E}_{\mathrm{orb}}},\end{array}$$(7)

where α_CE is the envelope ejection efficiency and E_bind the binding energy of the envelope (Webbink 1984; Moreno et al. 2022),

$$\begin{array}{cc}\hfill E_{\mathrm{bind}}& =-{\displaystyle {\int }_{M_{\mathrm{core}}}^{M_{\mathrm{surf}}}}\frac{\mathit{Gm}}{r}\mathrm{d}m+{\alpha }_{\mathrm{th}}{\displaystyle {\int }_{M_{\mathrm{core}}}^{M_{\mathrm{surf}}}}u\mathrm{d}m\hfill \end{array}$$(8) $$\begin{array}{cc}& \equiv -\frac{G{M}_1{M}_{\mathrm{env}}}{\lambda R_1},\hfill \end{array}$$(9)

and ΔE_orb the released orbital energy,

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{orb}}=-\frac{G{M}_{\mathrm{core}}{M}_2}{2a_{\mathrm{f}}}+\frac{G{M}_1{M}_2}{2a_{\mathrm{i}}},\end{array}$$(10)

we can calculate the final post-CE orbital separation for α_CE and α_th chosen to be unity. Here, M_core is the primary-core mass, u the internal energy (containing recombination energy), m the mass coordinate with corresponding radius r, and λ as well as α_th parameters accounting for the available binding energy of the envelope (de Kool 1990) and the fraction of total internal energy aiding to unbind the envelope. We note that the integral in Eq. (8) is performed from an a priori unknown mass coordinate of the stellar core to the surface of the star. This intrinsically raises the problem of defining the core and, hence, the envelope of the star and is often approximated in binary evolution calculations by, for instance, the maximum compression point or the point, where the hydrogen abundance drops below X = 0.1 (see Sect. 2.1, Moreno et al. 2022).

The binding energy of the envelope E_bind is obtained via integrating Eq. (8) over the one-dimensional MESA stellar profile Moreno et al. (2022), Fig. 1b) and assuming full accessibility of the internal energy reservoir (i.e., α_th  =  1). Here, the stellar core mass in Eq. (8) (and consequently the envelope mass in Eq. (9)) is chosen to be defined by the compression point of the star (Sect. 2.1). By rearranging Eq. (7) and Eq. (10) and further assuming α_CE  =  1, the expected final separations are a_f,  BH  ≈  19.2 R_⊙ and a_f,  NS  ≈  5.9 R_⊙ for the BH and NS companion, respectively (Moreno et al. 2022). These values are in contrast to the larger final separations of our CE simulations of a_f,  BH  ≈  37 R_⊙ and a_f,  NS  ≈  14 R_⊙. This discrepancy was also found in lower mass systems (e.g., Ohlmann et al. 2016a; Iaconi et al. 2017; Sand et al. 2020; Kramer et al. 2020).

This difference may be because the maximum compression point is not contained on our simulation domain, leading to a lower envelope binding energy. In the following, we use the binding energy of the envelope from the relaxed AREPO RSG progenitor model, which is up to a factor of 6 smaller (see Moreno et al. 2022, Fig. 1b). It is therefore plausible for the above estimate to yield final orbital separations reduced by the same factor. The integral is performed from the mass coordinate defined by the initial cut radius (M_core  =  2.97 M_⊙) to the surface of the star (M₁  =  9.4 M_⊙). For α_th  =  0.0,  0.5 and 1.0 we find through Eq. (9) λ  =  0.483,  0.763 and 1.810, respectively (Table 1, Moreno et al. 2022).⁶

Combining Eq. (9) and Eq. (10) yields

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{CE}}=\frac{G{M}_1{M}_{\mathrm{env}}}{\lambda R_1}{[\frac{G{M}_{\mathrm{core}}{M}_2}{2a_{\mathrm{f}}}-\frac{G{M}_1{M}_2}{2a_{\mathrm{i}}}]}^{-1},\end{array}$$(11)

which results for our measured final orbital separations (Sect. 3.1) in α_CE,  BH  =  0.5 − 1.86 and α_CE,  NS  =  0.57 − 2.12, depending on the choice of α_th (see Table 1).

3.5. Final fate of the systems

Given the final orbital separations of a_f,  BH  ≈  37 R_⊙ and a_f,  NS  ≈  14 R_⊙, and eccentricities of e_BH  =  0.12 and e_NS  =  0.055 for the BH and NS companion, respectively, we find merging times of T_{merge,  BH}  =  2320 Gyr and T_{merge,  NS}  =  325 Gyr by GW emission (Peters 1964). Hence, the systems are not expected to merge within a Hubble time of 13.787 Gyr (Planck Collaboration VI 2020).

The extended lifetime of the post-CE binaries set by the – compared with the CE phase – inefficient inspiralling process due to GW emission and the remaining lifetime of the RSG remnant of 6 × 10⁴ yr (Moreno et al. 2022) opens up several evolutionary pathways, that could potentially alter their fates tremendously. These include a second stable MT phase initiated by the re-expansion of the primary core (e.g., Tauris et al. 2017; Woosley 2019; Laplace et al. 2020; Jiang et al. 2021; Vigna-Gómez et al. 2022; Wei et al. 2024), the interaction of the engulfed post-CE binary with the left-over bound material in the orbital plane (e.g., Gagnier & Pejcha 2023; Tuna & Metzger 2023; Wei et al. 2024) and the SN event of the primary core forming a NS (Sect. 2.1) and a connected NS kick (e.g., Brandt & Podsiadlowski 1995; Wei et al. 2024). The results of Wei et al. (2024) underpin the inevitability of the MT as well as the SN event, narrowing down the possible pathways.

To account for the MT, we closely follow the analysis of Moreno et al. (2022) and consider the companions to be genuine compact objects with Eddington-limited accretion.⁷ This effectively leads to almost completely nonconservative MT. For quantifying the evolution of the orbital separations, we utilize Eq. (3) of Tauris (1996) and, for example, Eq. (8) of Podsiadlowski et al. (2002). The impact of an MT event on the final orbital separations and the merging time, given our assumptions, is listed in the second row of Table 2. Notably, the systems involving a NS show a decrease in orbital separations, whereas the systems with the BH companions expand.

The interaction of the binary with the left-over bound material – unaccounted for in the preceding work of Moreno et al. (2022, Sect. 4) – can drastically alter the orbital separation (e.g., Shi et al. 2012; Tiede et al. 2020; Siwek et al. 2023b; Gagnier & Pejcha 2023; Tuna & Metzger 2023; Wei et al. 2024) and therefore the final fate of the binary on timescales of 10³ − 10⁵ yr (when photo-evaporation is expected to become important) even before the MT event takes place. In fact, the presence of a thick CBD may even harden massive binaries down to distances, where the emission of GWs becomes important or lead to a merger before the SN of the remnant He core (Dittmann & Ryan 2022; Wei et al. 2024), depending on the mass and the lifetime of the disk. In the following, we suppose that recombination-driven ejection is not perfectly efficient due to nonadiabatic effects. In this scenario, the remaining material can form a long-lived CBD after the dynamic CE phase that interacts with the engulfed core binary on timescales much longer than the simulated duration (Sect. 3.1). The orbital change due to CBD interaction in our systems can then be described by parameterized torque prescription provided in Tuna & Metzger (2023),

$$\begin{array}{c}\hfill \frac{a(t)}{a_0}={(1-\frac{M_{\mathrm{acc}}}{M_{\mathrm{binary},0}})}^{2\chi {(1+q)}^2/q-3},\end{array}$$(12)

where a₀ is the initial orbital separation, M_acc(t) the accreted mass, M_binary, 0  =  M_core + M₂ the initial binary mass, q  =  M₂/M_core the binary mass ratio, and χ a dimensionless parameter quantifying the torque on the binary per unit accreted mass. We further express the accreted mass by the inflow rate from Sect. 3.3 (Fig. 5) such that M_acc ≈ Ṁ_acc t_disk. Here, we assume that the accretion rates onto the core binary are 1 % of the inflow rates (Sect. 3.3, given that most of the inflowing material is ejected through the magnetically driven outflow) and that the expected lifetime of the disk t_disk  =  10⁴ yr is set by photo-evaporation (e.g., Tuna & Metzger 2023). In this way, the total accreted mass amounts to 0.73 M_⊙ and the only unknown variable left is χ. Since χ depends on the mass accretion onto the central binary (preventing us from computing it directly) and the literature is currently lacking any precise values, we follow the results of Siwek et al. (2023a,b), who found that systems with q  ≳  0.3 tend to shrinkage, and assume χ  =   − 2χ_crit  =   − 3q/(1 + q)² as an approximation for all our simulations (χ_crit is the critical value dividing widening for χ >  χ_crit and broadening for χ <  χ_crit of the binary orbit, Tuna & Metzger 2023). It must be noted that these estimates of the further evolution are rather a qualitative approximation, given the large uncertainties. The impact on the final distance and, thus, the merging times of the systems are summarized in the third row in Table 2, where we assumed an equilibrium eccentricity of e  =  0.5 (Siwek et al. 2023a).

For the SN kick, we adhere to the analysis in Brandt & Podsiadlowski (1995, Sect. 2.1) assuming a Maxwellian distribution of kick velocities with σ  =  265 km s⁻¹ (Hobbs et al. 2005) and the formation of a 1.4 M_⊙ NS.⁸ Under these assumptions, however, the natal NS kick with the kick parameters as described above alters the final fate of the core binary system, as illustrated in Fig. 6. In this figure, the probability of the system to merge via GW emission (fraction f_merge) is shown in purple, the fraction of bound (but not merging) systems f_bound is shown in orange, and the fraction of systems disrupted by the SN kick f_disrupt is indicated in blue. As possible pathways for the systems, we consider the following two scenarios:

Fig. 6. Fractions in percent of systems that are considered to merge within the Hubble time fmerge (light purple), expected to stay in a bound binary configuration and are still too wide to merge within a Hubble time fbound (orange) and to be disrupted fdisrupt (light blue) by the natal NS kick. Considered are the following evolutionary pathways after the CE event: Either the post-CE binary first experiences an MT episode followed by the SN (MT+SN) or the CBD interaction and an MT episode transitions into the SN (CBD+MT+SN). All events are modeled separately, and the assumptions for the calculations are provided in Sect. 3.5.

• (i)

  After the CE evolution, an MT phase occurs, followed by the SN explosion of the primary core, excluding interaction with the CBD. The MT induces orbital widening in binaries with a BH companion, while for systems involving a NS companion, the separation decreases (Table 2). The merger rates for scenarios involving a BH companion amount to 0.9% and 10.7 % for those with a NS companion. Compared to Moreno et al. (2022), the merger rates increase by 0.3 %pt and 3.0 %pt for the BH and NS scenarios, respectively.

• (ii)

  We consider the CBD interaction and an MT episode ending with the SN event. Physically, the CBD and MT event can take place simultaneously (Wei et al. 2024; Tuna & Metzger 2023), but we choose to model them as subsequent phases. The CBD interaction then first hardens all separations and the MT further tightens the systems involving NS companions while widening those with BH companions (Table 2). As a result, the systems with NS companions increase their merger rates once more to 54.7%, while the merger fraction in the systems containing a BH enhances the merger fraction to 4.9%.

In summary, the interaction between the core binary and a CBD can significantly tighten the binaries and increase the likelihood of a merger event, but the individual results should be interpreted cautiously due to the considerable uncertainties arising from our simplified assumptions. A comprehensive exploration of the parameter space and the development of more sophisticated models, such as those proposed by, for example, Wei et al. (2024), are essential for obtaining more detailed insights.

4. Discussion

4.1. Evolution of the orbital separations

The tension between the orbital separations resulting from CE interactions according to expectations of the α formalism and the outcomes of our simulations (Sect. 3.4) is not unique to our modeled scenarios. In fact, the findings in lower mass systems (e.g., Passy et al. 2012; Ricker & Taam 2012; Ohlmann et al. 2016a; Iaconi et al. 2017; Prust & Chang 2019; Sand et al. 2020; Reichardt et al. 2020; Ondratschek et al. 2022) or more massive systems (e.g., Moreno et al. 2022; Lau et al. 2022a) reach similar conclusions. In contrast to the energy prescription of the CE evolution and in agreement with, for instance, Sand et al. (2020), after the initial perturbation of the envelope by the rapid spiral-in of the companion, the envelope ejection induced by the release of recombination energy barely affects the orbital separations (Sect. 3.1). Consequently, our final orbital distances significantly exceed the separations expected from the energy formalism (Sect. 3.4). Even under full envelope ejection, the distances only marginally decrease compared to Moreno et al. (2022).

Physically, the strongly decreased orbital contraction rates after the spiral-in of the core binary are caused by the ceasing of the dynamical drag force (e.g., Ostriker 1999), which is induced by the rapidly decreasing density of the gas in the vicinity of the core binary system, together with the subsonic and partially corotating motion of this gas (e.g., Moreno et al. 2022; Röpke & De Marco 2023). The reduced drag is then operating on the binary until the system dynamically equilibrates, and the further evolution is set by processes on longer timescales (e.g., Ivanova et al. 2013; Hirai & Mandel 2022). Given the occurrence of a CBD directly after the rapid in-spiral, this picture must be complemented by gravitational torques (or other stress-related angular momentum extraction processes) and the accretion of mass and angular momentum onto the primary core and companion, which could affect the post-spiral-in contraction. These processes could even proceed on longer timescales and are shown to impact the further evolution of the core binary significantly (Wei et al. 2024).

Our simulations with κ  =  1 show good agreement with the simulations involving κ  =  0 (Sect. 3.1 and Appendix A). These, in turn, show convergence in the (averaged) orbital evolution and envelope ejection with regard to the initial cut radius according to the resolution study presented in Fig. 9 of Moreno et al. (2022).⁹ It therefore appears plausible to conclude that the found (averaged) final orbital separations are indeed not only an upper limit but reflect the expected outcome given the considered timescales and physical model. The discrepancy to the α-CE formalism can thus only be mitigated by reporting updated envelope ejection efficiencies (e.g., Sect. 3.4).

Most numerical and physical uncertainties concerning the orbital separations are related to approximating the primary core and companion as gravity-only point particles (Sect. 2.1). We focus here on arguments that were not accounted for in previous publications, such as the adaptive reduction of the softening length (Moreno et al. 2022) and the reaction of the primary core to envelope loss (e.g., Wei et al. 2024, see also Bronner et al. 2024), and refer to these publications for a more detailed discussion.

As stated in Sect. 2.1, with the gravity-only point particle approximation, we neglect any accretion of angular momentum and mass onto the primary core and the companion during our simulation. However, this is expected to leave only a weak impact on the orbits over the dynamical CE interaction. Quantitatively, this can be shown by applying Eq. (12) and assuming the extreme case that the inner binary can accrete mass an angular momentum on a similar rate than the inflow rates. Assuming further, that this process proceeded from the onset of the magnetized outflows to the end of our simulations (i.e., M_acc  ≈  2 × 10⁻⁵ M_⊙d⁻¹[t_end − t₀], where t_end is the total simulated time and t₀ is set by the onset of the magnetized outflows) and χ  =   − 3q/(1 + q)² as in Sect. 3.5, this results in deviations from the final orbital separation on the order of 0.1%. Hence, the accretion onto the central binary seems to be unimportant for the overall evolution during the dynamical CE interaction.

We find changes in the eccentricity during the late stages of our simulations (as can be qualitatively inferred from the fluctuations in the orbital separations in Fig. 1 and seen in Fig. B.2). As already discussed in Moreno et al. (2022), our simulations do not converge in eccentricity, and the reported values in Sec. 3.5 must be taken with caution. The arising changes seem to be inherently connected to the point-particle approximation of the companion and primary star core, but there is currently no coherent explanation for this behavior. One potential cause of these changes might be the adaptive reduction of the softening length, which may introduce perturbations to the HSE (if present). However, the abrupt changes in the eccentricities do not correlate with the reduction steps in the softening length (compare Fig. 1 and Figs. A.1a and A.1b) and the total energy is conserved down to a few per mil throughout our simulations (relative to the initial total energy in the system the deviations amount to 0.11% for the BH and 0.16% for the NS companion case, respectively). Moreover, the eccentricity changes are mostly observed in the later stages of the simulations, where the softening length either changes only smoothly or is constant. Therefore, a direct connection between the changes in eccentricity and the changes in the gravitational softening length seems unlikely. While some other numerical artifacts cannot be excluded as a reason, the changes in eccentricity might have a physical origin, for example, from the complex interaction of the bound material and the central binary, similar to the effects associated with warped disks (e.g., Nixon & King 2016).

4.2. The common-envelope ejection and the disk-like bound material

In Sect. 3.2 and Sect. 3.3, we found evidence for the presence of a thick CBD surrounding the core binary after the rapid spiral-in. The material showed characteristics of a thick, pressure-supported sub-Keplerian disk. We further identified two envelope ejection channels in our simulations, and in contrast to the spherical expansion, the channels are tightly coupled to the CBD. The first channel is a disk wind driven by recombination and stripping the outer layers of the torus-like structure. The second channel is the ejection of material through a magnetically driven jet-like outflow in the polar directions, which is refueled by the CBD material.

Since the first channel is driven by recombination, it depends sensitively on the coupling of the radiation field to the hydrodynamic evolution. Our simulations are fully adiabatic and assume local thermalization of recombination energy. Ultimately, this leads to the envelope ejection, as shown in Fig. 1, predominantly via this channel. Physically, however, the photosphere is expected to be near the hydrogen recombination front (Fig. 4) and some of the energy may be lost by radiation without performing mechanical work. Hence, the ejection efficiency by the radial winds could be overestimated. Therefore, a more advanced model incorporating radiative cooling processes is necessary to address this issue comprehensively.

While the question of how much hydrogen recombination energy is accessible to the system is highly speculative, one could also turn the question around and ask how much of the hydrogen recombination energy is required to accelerate material at the hydrogen recombination front beyond the local escape velocity. For the BH companion, this question can be approximated by equating the specific energies G(M_c + M₂)/2R and Δe_{HII → HI} at the hydrogen recombination front (here denoted by R; c.f. argument in footnote 3). For R ≈ 4  × 10³ R_⊙ this implies that almost 21% of the recombination energy is required to eject the material according to the kinetic-energy criterion and as the recombination front recedes in the evolution, this number becomes progressively higher (e.g., at around t ≈ 13 698 d, R ≈ 10³ R_⊙ and thus 83% are needed). This leads us to the conclusion that, under the influence of cooling, some of the disk material can remain bound and form a long-lived CBD. Assuming further that the disk winds can indeed be lowered or even ceased by cooling (depending on the cooling timescale compared to the ejection timescale), the mass transport inward can still be maintained by the persistence of turbulent transport mechanisms (like the saturated MRI phase) or angular momentum extraction by gravitational torques inflicted by the engulfed binaries. In this scenario, the second mass-loss channel through the jet-like polar outflow could become the only way for the system to eject material.

While we analyze the properties of the magnetically driven outflow in detail in Vetter el al. (in prep.), we can already conclude from our analysis of the inflow rates and the disk wind in Sect. 3.3 that the envelope mass lost through this channel has to be subdominant during our simulations. So far, it is not resolved by which mechanisms angular momentum is transported, and a more elaborate analysis (as attempted in e.g., Gagnier & Pejcha 2023, 2024) is required to address this question. One prominent candidate in our case is the saturated phase of the MRI, but further investigations are required to support this hypothesis.

The mass inflow rates might be affected by various other parameters. For instance, the initial mass ratio q  =  M₂/M₁ of the binary can affect the magnetically driven outflows. In our scenario involving the NS companion, less material remains bound near the core binary at the onset of the magnetically propelled outflow (Fig. 1) and, consequently, less mass can be ejected through the magnetized outflows. Hence, this second envelope-loss channel is expected to be less important compared to the CE event with the BH. Likewise, the initial mass of the primary star (i.e., the amount of mass that is participating in the CE event) and its evolutionary stage (i.e., the binding energy of the envelopes, which are expected to be lower for less evolved stars) should directly impact the evolution. This might explain why Ondratschek et al. (2022) find the magnetically driven outflows to be unimportant for the overall CE ejection process.

The lower resolution in low-density regions because of the mass-adapted unstructured mesh generated by AREPO implies rather poor resolution of the pressure scale height (we find 3 − 60 cells per H) in the CBD material (Sect. 3.3). This can enhance numerical viscosity and may lead to overestimated α-viscosities.

5. Conclusions

In this paper, we turned our attention once more to the CE interaction of a 10 M_⊙ RSG with a 5 M_⊙ BH as well as a 1.4 M_⊙ NS companion, that was first simulated by Moreno et al. (2022). Utilizing a new refinement approach for the softened primary core and companion particles, we were able to follow the dynamic evolution until full kinetic envelope ejection. Our main results can be summarized as follows:

• After 1220 (1831) orbits, a final orbital separation of a_f,  BH  ≈  37 R_⊙ (a_f,  NS  ≈  14 R_⊙) is established for the BH (NS) companion. In terms of the energy prescription, these final distances translate to α_CE,  BH  =  0.50 and α_CE,  NS  =  0.57 for the BH and NS companions, respectively. Despite full envelope ejection, the orbital distances we obtain in our simulations exceed the expected separations of a_f,  BH  ≈  19.2 R_⊙ and a_f NS  ≈  5.9 R_⊙ by the α-formalism (here α_th = 1, α_CE  =  1). In agreement with simulations of CE interactions in low-mass systems (e.g., Sand et al. 2020), the envelope ejection induced by recombination in the post-spiral-in systems has little effect on the final orbital separations.

• The morphology of our system changes significantly during our conducted CE simulations. The initial perturbation of the envelope during the spiral-in of the companions breaks the spherical symmetry of the primary star and leaves behind a system imprinted by the symmetry of the orbital motion of the core binary. Consequently, after the spiral-in, the bound material forms a toroidal structure visually resembling a geometrically thick CBD and fulfills the dynamic requirements (i.e., gravitationally bound, circularized motion, Raleigh stability and vertical HSE) associated with stable disk structures. The bound material is found to be eroded by two mass loss channels: radial inward mass fluxes and recombination-driven winds, which drive the ejection of envelope material in the opposite direction.

• The first channel is refueling magnetically driven polar outflows appearing in both of our simulations¹⁰ similar to those observed in simulations of lower mass systems (Ondratschek et al. 2022). They accelerate high-entropy material to peak velocities of 300 km s⁻¹, which flows along low-density channels perpendicular to the orbital plane. When the polar outflows interact with the surrounding material, bow shocks are observed. Together with the toroidal structure, they do resemble the morphological appearance of bipolar pre-PN and the magnetically driven outflow can act like a pressure valve in our systems.

• From the radial (averaged) distributions of the bound material, we can infer the time evolution of physical properties such as pressure scale height, aspect ratio, temperature, column density and angular velocity. They show the characteristics expected from pressure-supported, sub-Keplerian thick disks with an aspect ratio of ∼0.6 − 1.0 (see Sect. 3.3, especially Fig. 4). We additionally estimate the radial mass fluxes through both mass loss channels (see Fig. 5). The inward-directed mass transport rates can be reproduced by a turbulent transport model with an α-viscosity of 0.06 −− 0.13. It appears likely that the stress-related transport is induced by the nonlinear regime of the MRI, and most of the material is expected to be launched by the magnetically driven outflow (Vetter et al. in prep.). Further, we show that the second channel of radial winds eroding the disk-like structure in radial direction in the outer part of the CBD is driven by the local thermalization of hydrogen recombination energy. This channel is found to exceed the inward flows significantly and is thus be expected to be the dominant mass ejection channel. It can potentially be affected by cooling processes, and further investigations, including cooling models, are required to refine our understanding of the properties of a post-CE CBD.

• We analyzed the possible final fates of the systems, similar to the analysis in Moreno et al. (2022). In addition to the default scenario of a subsequent MT and a natal NS kick, we also accounted for the interaction with a CBD. The fractions f_merge of mergers due to GW emission amount to 0.9–10.7% for the BH and NS systems when accounting for a stable MT after the CE interaction and a natal NS kick (see Sect. 3.5). The additional interaction of the engulfed binary with a CBD (Sect. 3.5, Tuna & Metzger 2023; Wei et al. 2024), however, increases the merger fractions significantly, and we found merger fractions of approximately 10(40)% for the system involving a BH (NS) companion. This agrees with the results obtained in Wei et al. (2024), and we refer to this publication for a more detailed analysis of the possible fate of the post-CE systems resulting from our simulations.

Acknowledgments

We thank the referee for the constructive and helpful comments, which helped us to improve the paper. M.V., F.K.R., F.R.N.S., M.Y.M.L., and R.A. acknowledge support by the Klaus-Tschira Foundation. M.Y.M.L. has been supported by a Croucher Foundation Fellowship. M.V., F.K.R., and R.A. acknowledge funding by the European Union (ERC, ExCEED, project number 101096243). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

References

Appendix A: A new refinement approach for gravity-only particles in CE interaction

In Sect. 2.3, we introduced a new refinement approach within the softened gravitational potential of the two point particles representing the core of the primary star and the companion. As mentioned, the softening length h is reduced over the course of our CE simulations (Fig. A.1a and A.1b) as the two point particles approach each other in order to avoid a potential overlap of the softened spheres. As a consequence of the implementation of this procedure, the softening length can be expressed as a monotonic decreasing function depending on the reduction given by a(t) (Fig. A.1). Combined with Eq. (4), it is immediately apparent that the cell sizes within the softened potential (R_cell  ≈  h(t)/N_CPS) must be constant for κ  =  1 and decrease in the κ  =  0 case (Fig. A.1c and A.1d).

The system time step Δt_sys in Fig. A.1e and A.1f seems to reflect the behavior of the smallest cell sizes according to the CFL condition (Courant et al. 1928) for the NS case. But we also observe an unexpectedly sharp decline in the system involving the BH companions for the simulations with κ  =  0 (black line at around 3000 d). This steep drop can be explained by the onset of the magnetized outflow (Vetter et al. in prep.), in which strongly magnetized and high velocity but low density material (i.e., large cell sizes) propagates through a rather weakly magnetized gas. Given the implementation of magnetohydrodynamics in AREPO (see Sect. 2), the relatively harsh difference in magnetization between the magnetically driven outflow and ambient material provokes a steep drop in the system time step. Additionally, for the NS companion and κ  =  0 (black line), there is a small recovery of Δt_sys at around, 500 d which is caused by a restart of the simulation with a higher C_CFL-factor. This adjustment increased the system time step but compared to κ = 0, the new refinement approach allows us to use time steps four times longer.

Fig. A.1. Temporal evolution of the softening length h in (a) and (b), the smallest cell radius within the softened spheres min(Rcell)|h in (c) and (d) as well as the system time step Δtsys in (e) and (f). The simulations conducted with κ = 1 (see Sect. 2.3) are shown in light blue, while the simulations with κ = 0 (Moreno et al. 2022) are shown in black. The plots on the left feature the system involving a BH companion and the plots on the right the CE models involving the NS companion, respectively.

Over the course of the simulations, the number of cells per pressure scale height in the BH (NS) case decreases from ≈11 to ≈7 (≈3) for κ  =  1 and ≈5 (≈7) for κ  =  0, respectively. These numbers may raise concerns and, for example, Ohlmann et al. (2017) suggests a linear resolution of the pressure scale height of 3.2 −− 6.8 cells for a stabilized atmosphere down to mach numbers of ℳ  =  0.01 −− 0.001 (with adiabatic index γ  =  5/3 and CFL factor of C_CFL  =  0.4). But this neither refers to the stellar remnant in our model, nor does it take the modified structure of the core into account. In conclusion, the new refinement approach mitigates the issue of decreasing time steps in the late stages of the modeled CE interactions, but the constant cell sizes may also introduce additional challenges in achieving convergence in orbital separation and envelope ejection with respect to the choice of the initial cut radius.

Appendix B: Influence of the new refinement criterion on the orbital separation

The new refinement criterion described in Sect. 2.3 enabled us to follow the entire envelope ejection in both analyzed systems. When comparing our simulations to the findings of Moreno et al. (2022) (see Fig. 1 and bottom row of Fig. B.1), the simulations seem to be consistent with the previous results, particularly in the earlier stages of the CE interaction. This affirmation gains further support when we plot the two orbital separations for κ  =  1 and κ  =  0 against each other (see top panels of Fig. B.1). The best linear fit yields a slope of ∂a_{κ  =  1}/∂a_{κ  =  0}  ≈  0.998 (≈0.995) for the simulations involving a BH (NS) companion. However, after the rapid spiral-in phase in both scenarios, the orbital separations begin to exhibit artificial eccentricities, resulting in deviations from our reference simulations. In total, the (averaged) orbital separations at the time when Moreno et al. (2022) terminated their simulations are a_{κ = 1 BH}  = 45.6 R_⊙ and a_{κ = 0 BH}  = 46.2 R_⊙ for the BH companion and a_{κ = 1 BH}  = 16.9 R_⊙ as well as a_{κ = 0 BH}  = 15.2 R_⊙ for the NS companion.

The variations in eccentricity are most likely due to the sudden change of the gravitational potential, induced by the decrease of the softening length. At the point of reduction, the cells within the former softened spheres might be pushed out of HSE and (depending on the decline in h and the position of the cells in the softened sphere) change their gravitational potential instantaneously, resulting in a restructuring of the density distribution in the vicinity of the point particle. In contrary to this argument, the total energy of the system is preserved down to per mil level and this effect can thus only have a minor effect on the system. Moreover, the realized cells must exhibit a significant, spatially asymmetric mass distribution at the moment of reduction to exert a kick on the binary. And last but potentially the strongest argument, the arising eccentricities seen in our simulations appear to be more pronounced in the late stages of the simulation. There, the temporal evolution of the softening length and the arising changes in eccentricity (e.g., Fig. 1) seem to show no direct correlation, and while the softening length remains constant, clear changes in eccentricity can be observed (see qualitatively in Fig. 1b). Regardless, the numerical approximation of the core region presents a source of uncertainty in current state-of-the-art simulations of massive CE interactions. Addressing this challenge is crucial for improving such simulations.

Fig. B.1. Comparison between the orbital separations in our simulation with κ = 1 and the fiducial model with κ = 0 (Moreno et al. 2022), for the simulation involving the BH (left panel) and NS (right panel) companion, respectively. In (a) and (b) a scatter plot of the orbital separations (aκ  =  1 and aκ  =  0) is shown. The black dashed lines represent the perfect correlation between the separations (i.e., aκ  =  1  =  aκ  =  0), while the best linear fit is aκ  =  1  ≈  0.998 × aκ  =  0 (≈0.995 × aκ  =  0) for the BH (NS) companion. In the bottom row (c and d), we show the zoomed in temporal evolution of the orbital separations for κ  =  0 (red) and κ  =  1 (blue) similar to Fig. 1.

Fig. B.2. Temporal evolution of the eccentricity. The simulation involving the BH companion is colored blue, while the system with the NS companion is plotted red. The eccentricity is computed via the Runge–Lenz vector.

Appendix C: Movies of the CE interactions

In Table C.1 we summarize the movies of the density evolution for both simulations shown in this publication.

All Tables

All Figures

Fig. 1. Temporal evolution of the orbital separation a (solid dark) and the fraction of unbound material fej according to the kinetic (dashed green, Eq. (1)), thermal (dashed rust, Eq. (2)), and internal energy (orange, Eq. (3)) criterion, respectively. The simulation involving the BH is shown on the left, while the NS scenario is on the right. The thin lines in the plots represent the simulations that contained the new refinement approach (see Sect. 2.3, κ  =  1), while the thick, faint lines correspond to the results of Moreno et al. (2022). On the left-hand side, we additionally show the unbound mass fraction only for cells contained within a radius of R  =  1.6 × 104 R⊙ around the central binary due to the artificial “rebinding” (see Sect. 3.1) of material in the pseudo-vacuum for t  ≥  4400 d (see Sect. 3.1). Furthermore, the hourglass markers represent the individual starting points of the magnetically driven outflow in the simulations with κ  =  1. The minor ticks in both plots in between 1 − 10 orbits are adapted for readability, and only the 5 orbit mark is shown.

In the text

Fig. 2. Slices of spherical radial velocity vr, sph. in the x − z plane in (a)-(f) and in the x − y plane in (g)-(j) for the BH (NS) companion in the left (right) column. The iso-surfaces in (a) and (b) enclose the bound mass according to the kin. energy criterion (Eq. (1)) in dashed gray. The gray streamlines in (c)-(j) represent the velocity field in the corresponding plane. The primary core (companion) is marked with a black cross (plus sign).

In the text

Fig. 3. Specific entropy s and density ρ slices of the simulation with the BH companion at t  =  8244 d. The quantities are shown to the left (x ≤ 0 for the spec. entropy) and right (x ≥ 0 for the density) in an edge-on view (x − z) in (b), respectively. Furthermore, the distribution of all cells along the z-direction in a cylindrical shell (indicated by the dashed rectangle) in (a) and (c). The gray streamlines in (b) show the velocity field in the x − z plane.

In the text

Fig. 4. Time evolution of averaged density ⟨ρ⟩ in (a), cumulative mass Mcum. in (b), specific angular momentum J2 in (c), angular velocity ⟨Ω⟩ in (d), isothermal pressure scale height H in (e), the aspect ratio H/r in (f), the temperature T in (g), the column density Σ in (h), the ionization fraction of HII (χHII), HeII (χHeII), and HIII (χHeIII) in (i), (j), and (k). Last, the average optical depth of cells τcell is shown in (l). The averages are taken over the range of r  ∈  [2a(t), 4 × 103 R⊙] (Sect. 3.3).

In the text

	Fig. 5. Radial distribution of mass transport rates $\dot{M}=2\pi r{\overline{v}}_r\mathrm{\Sigma }$ similar to Fig. 4. The faint bands over the entire radial domain in the background correspond to the overall envelope ejection rate Ṁej of the system at the corresponding time with the matching color.
In the text

Fig. 6. Fractions in percent of systems that are considered to merge within the Hubble time fmerge (light purple), expected to stay in a bound binary configuration and are still too wide to merge within a Hubble time fbound (orange) and to be disrupted fdisrupt (light blue) by the natal NS kick. Considered are the following evolutionary pathways after the CE event: Either the post-CE binary first experiences an MT episode followed by the SN (MT+SN) or the CBD interaction and an MT episode transitions into the SN (CBD+MT+SN). All events are modeled separately, and the assumptions for the calculations are provided in Sect. 3.5.

In the text

Fig. A.1. Temporal evolution of the softening length h in (a) and (b), the smallest cell radius within the softened spheres min(Rcell)|h in (c) and (d) as well as the system time step Δtsys in (e) and (f). The simulations conducted with κ = 1 (see Sect. 2.3) are shown in light blue, while the simulations with κ = 0 (Moreno et al. 2022) are shown in black. The plots on the left feature the system involving a BH companion and the plots on the right the CE models involving the NS companion, respectively.

In the text

Fig. B.1. Comparison between the orbital separations in our simulation with κ = 1 and the fiducial model with κ = 0 (Moreno et al. 2022), for the simulation involving the BH (left panel) and NS (right panel) companion, respectively. In (a) and (b) a scatter plot of the orbital separations (aκ  =  1 and aκ  =  0) is shown. The black dashed lines represent the perfect correlation between the separations (i.e., aκ  =  1  =  aκ  =  0), while the best linear fit is aκ  =  1  ≈  0.998 × aκ  =  0 (≈0.995 × aκ  =  0) for the BH (NS) companion. In the bottom row (c and d), we show the zoomed in temporal evolution of the orbital separations for κ  =  0 (red) and κ  =  1 (blue) similar to Fig. 1.

In the text

	Fig. B.2. Temporal evolution of the eccentricity. The simulation involving the BH companion is colored blue, while the system with the NS companion is plotted red. The eccentricity is computed via the Runge–Lenz vector.
In the text