1. Introduction
Mass transfer in binary systems involving a giant donor and a more compact companion remains one of the most significant ‘black boxes’ in modern astrophysics, with implications on our understanding of any evolved binary population (e.g. De Marco & Izzard Reference De Marco and Izzard2017). When a giant star fills its Roche lobe, the stability of the resulting mass transfer is dictated by the response of the donor’s radius to mass loss compared to the response of the Roche lobe itself (see e.g. Soberman et al. Reference Soberman, Phinney and van den Heuvel1997 and Hjellming & Webbink Reference Hjellming and Webbink1987). If the star expands faster than the lobe, the process becomes unstable, leading to a common envelope (CE) phase – a brief, catastrophic episode where the companion is engulfed, leading to dramatic orbital shrinkage or stellar mergers (Ivanova et al. Reference Ivanova2013).
Unlike stars with radiative exteriors, which tend to contract or stay nearly the same size upon mass loss, convective donors have an adiabatic response that causes them to expand as they lose their outer layers.
This expansion creates a feedback loop: as the star grows, it overfills its Roche lobe even further, leading to even higher mass-loss rates. If the Roche lobe cannot expand fast enough to accommodate the growing star – which typically happens if the donor is significantly more massive than its companion – runaway mass transfer ensues, leading to a CE phase.
In this paper we consider a class of post-red giant branch and post-asymptotic giant branch (post-RGB and post-AGB) stars with main sequence companions with semi-major axes of
$\sim$
$0.5$
–4 au (periods of
$\sim$
100–
$2\,000$
d), eccentricities up to 0.6 (Van Winckel Reference Van Winckel2025) and surrounded by a circumbinary disc (Oomen et al. Reference Oomen, Van Winckel, Pols and Nelemans2019).
The orbits of these post-RGB and post-AGB binaries are typically small enough to imply a recent phase of Roche lobe overflow (RLOF), because their progenitor RGB and AGB stars were larger than the orbital periastron distance. However, a RLOF with a giant should have resulted in a CE interaction with a resulting, much reduced orbital separation, which is not observed for these stars.
A possible solution was envisioned by Soker (Reference Soker2015), who suggested that if during the interaction the companion accretes gas via RLOF and blows jets, it can avoid the CE inspiral and remain in a perpetual phase of non-conservative mass transfer that results in the complete loss of the giant star’s envelope. Kashi & Soker (Reference Kashi and Soker2018) later suggested that this ‘grazing envelope’ mechanism can even produce eccentric binaries because of the enhancement of the mass accretion rate near periastron. These scenarios, however, have never been investigated further.
Another way in which the RLOF, mass transfer phase could have been stabilised and have resulted in the relatively wide separations of post-RGB/AGB binaries was suggested by the discovery that the mean companion mass for these binaries is
$1.09\pm0.62$
M
$_{\odot}$
(Oomen et al. Reference Oomen, Van Winckel, Pols and Nelemans2019), larger than the typical
$\sim$
0.5 M
$_{\odot}$
seen in equivalent post-CE, post-AGB binaries such as central stars of PN (Iaconi et al. Reference Iaconi, Maeda, De Marco, Nozawa and Reichardt2019). These larger companion masses imply larger
$q\equiv M_\textrm{2} /M_\textrm{1}$
values at the time of interaction. Larger q values are known to lead to wider post-CE separation in simulations (e.g. Passy et al. Reference Passy2012, although always below
$\sim$
30 R
$_{\odot}$
). It is also known that larger q values may stabilise mass transfer and lengthen the pre-inspiral mass transfer phase. One can then speculate that this may avoid a CE inspiral or lead to a weakened CE.
This said, the exact mass ratio,
$M_2/M_1$
, above which stability can be assumed is not well known. Tout (Reference Tout1991) found analytically that for a non rotating star, stable mass transfer would take place for a value of
$q\equiv M_2/M_1 \gtrsim 1.43$
Footnote
a
(here
$M_2$
is the accretor, and
$M_1$
is the donor).
Ge et al. (Reference Ge, Webbink, Chen and Han2015) mapped the response of donor stars to rapid mass loss and provided stability criteria for CE precursors. Using adiabatic mass-loss sequences from the zero-age main sequence to the base of the giant branch, they determined that the critical mass ratio for dynamical instability varies strongly with envelope structure. For fully convective donors, this corresponds to
$q \approx 1.5$
–
$1.6$
, whereas for radiative envelope donors
$q \approx 0.8$
, down to
$q \approx 0.3$
–
$0.5$
for more massive radiative donors. Approaching the asymptotic giant branch, the critical mass ratio rises back towards
$q \approx 1.3$
–
$1.7$
, making evolved stars far more prone to a CE inspiral.
Following on from this, Ge et al. (Reference Ge, Webbink, Chen and Han2020) showed that superadiabatic convective envelopes can enhance donor expansion and lower the dynamical stability, implying that purely adiabatic stability criteria likely underestimate the tendency of giants to undergo a CE inspiral. However, Woods & Ivanova (Reference Woods and Ivanova2011) found that this instability depends on assuming the superadiabatic surface layer is entirely removed; if mass transfer proceeds on a timescale comparable to the local thermal timescale, the layer can thermally readjust, yielding a more stable response than the adiabatic assumption predicts (see also Temmink et al. Reference Temmink, Pols, Justham, Istrate and Toonen2023). The question of where the stability-instability boundary sits is clearly not yet answered.
In this paper we set out to explore CE interactions between RGB stars and companions with a range of q values, monitoring the phase of RLOF before the CE inspiral, to determine its impact on the final separation and morphology of the interaction. We use the same star as adopted by Passy et al. (Reference Passy2012) and Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019), a 0.88 M
$_{\odot}$
, 90 R
$_{\odot}$
, RGB star. The latter work used a compact companion of mass 0.6 M
$_{\odot}$
(
$q=0.68$
), starting the simulation at the separation that just allows Roche lobe overflow. They found that material leaving the outer Lagrangian points before the inspiral appears to remain mostly bound in a circumbinary disc (see also MacLeod et al. Reference MacLeod, Ostriker and Stone2018a and MacLeod & Loeb Reference MacLeod and Loeb2020). This hinted at the possibility that with a larger companion mass and therefore mass ratio, one may obtain a wider separation as well as the formation of the observed circumbinary discs. In this paper we therefore extend their investigation to higher values of q to increase mass transfer stability, determine whether a circumbinary disc can form, whether a delayed inspiral occurs and whether the final separation will be wider.
This paper is structured as follows. In Section 2 we outline the simulations used and justify the parameters chosen. In Section 3, we analyse the binary orbital decay and unbound mass. We then look at the
$L_1$
mass transfer rates as a function of both resolution and mass ratio, q. We also investigate the mass lost during the pre-inspiral phase through the
$L_2$
and
$L_3$
points and examine the kinematic identity of early ejecta versus that of the ejecta produced during the inspiral, to determine if a disc formed during the early phase can survive the CE. Finally, we look at the possibility of the formation of a post-CE circumbinary disc around these surviving binaries by fall-back of bound gas. In Section 4 we discuss our findings and explore the implications for the formation of post-RGB binary systems, particularly those documented by Oomen et al. (Reference Oomen2018) and Kluska et al. (Reference Kluska2022). We summarise, conclude, and give details of future investigations in Section 5.
2. Method
In order to explore the early pre-CE phase of mass transfer and the characteristics of the post-CE ejecta, we performed 12 simulations using the smoothed particle hydrodynamics (SPH; e.g. Monaghan Reference Monaghan1992; Price Reference Price2012) code phantom (Price et al. Reference Price2018). Across these simulations, we have varied
$q \equiv M_2/M_1 = 0.68$
–
$1.5$
, the resolution (76 000 and 531 000 SPH particles), and the simulation’s equation of state (EoS; Table 1). Our donor star is the same as that used by Passy et al. (Reference Passy2012), Iaconi et al. (Reference Iaconi2017) and Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019): a 1 M
$_{\odot}$
main sequence star that was evolved to the RGB using the 1D code mesa (Paxton et al. Reference Paxton2010), till it had a total mass of
$M_1 = 0.88$
M
$_{\odot}$
with a core mass
$M_{c} = 0.392$
M
$_{\odot}$
and a radius of approximately 90 R
$_{\odot}$
. We model both our companion star and the primary’s core as point mass particles with a 3 R
$_{\odot}$
softened potential. The 1D stellar structure was mapped into the 3D computational domain. This stellar structure was shown to be stable in the 3D computational domain by Passy et al. (Reference Passy2012), Iaconi et al. (Reference Iaconi2017), and Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) by evolving the star in isolation and showing that it remains in reasonable hydrostatic equilibrium for several dynamical times. Later, González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022) and Lau et al. (Reference Lau2022) improved the stabilisation method, and this is the one we used for both EoSs, obtaining an even more stable star.
Simulations’ inputs. The second column is the number of SPH particles each simulation uses. The value
$q = M_2/M_1$
, where
$M_2$
is the companion mass and
$M_1=0.88$
M
$_{\odot}$
. The initial separation at Roche lobe contact is given by
$a_0$
. The value of
$h_{s}$
is the smallest SPH smoothing length at time
$t=0$
(noting that the gravitational softening radius for both the primary and secondary core is 3 R
$_{\odot}$
).
$m_{p}$
is the mass of all SPH particles in the simulation. The time
$t_\textrm{end}$
is the total simulation physical time. The last two columns are the artificial conductivity and the minimum artificial shock viscosity, respectively (the maximum artificial shock viscosity is
$\alpha_\textrm{max}=1$
). The differences in these last two columns between the ideal and tabulated EoS simulations are due to the different requirements for stellar stability.
L and H denote low and high resolution, respectively (see
$n_\textrm{part}$
).
I and M denote the use of an ideal and tabulated (mesa) EoS, respectively.
Our lowest mass ratio is
$q = 0.68$
, the same as used by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019). Additional simulations were carried out up to
$q=1.5$
. A binary with
$q=1.5$
implies a 1.32 M
$_{\odot}$
companion, more massive than the main sequence progenitor of the primary. Purely from an evolutionary point of view, we would have to assume the companion in this simulation to be a massive white dwarf star. Observed post-RGB and post-AGB binaries have main sequence, not white dwarf companions, so this simulation is not quite modelling one of those systems. We therefore consider our
$q=1.5$
simulation as a technical experiment to push the ratio to a higher value.
In each simulation, we start our binary with a circular
$e=0$
orbit with a semi-major axis separation such that the radius of the star is approximately equal to its Roche lobe radius. In doing so, the process of mass transfer from the donor onto the accretor begins immediately in our simulations. This also ensures each simulation begins at Roche lobe contact to model as much of the pre-inspiral phase of evolution as is computationally feasible.
Table 1 lists a summary of parameters for our 12 simulations, where the names are self-explanatory: the number is the mass ratio
$\times 100$
, ‘L’ and ‘H’ are low and high resolution, respectively, and ‘M’ stands for mesa EoS as opposed to ideal gas. For each of the four values of q we carried out a low- and a high-resolution simulations with an ideal gas EoS, as well as an equivalent high-resolution simulation with a mesa, tabulated EoS (Paxton et al. Reference Paxton2011; Paxton et al. Reference Paxton2013; Paxton et al. Reference Paxton2015; Reichardt et al. Reference Reichardt, De Marco, Iaconi, Chamandy and Price2020). This additional EoS simulation is key to providing an upper limit to the amount of gas unbound, due to its inclusion of recombination energy, which delivers additional thermal energy to the envelope. The effects of this EoS have previously been studied by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Chamandy and Price2020), González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022), and Lau et al. (Reference Lau2022).
We adopted default artificial conductivity parameter values (see Section 2.2.8 of Price et al. Reference Price2018) of unity for the ideal gas models and 0.1 for tabulated EoS simulations (see Table 1). The default value is set to ensure accurate treatment of contact discontinuities (Price Reference Price2008). The conductivity is second order in phantom, with the effective thermal conduction
$\kappa_\textrm{cond} \propto \alpha_u h^2 (\nabla \cdot \textbf{v})$
, where h is the resolution length (c.f., Price et al. Reference Price2018). It thus vanishes when the velocity divergence is zero and the resolution is high, but we found that the remaining heat transfer could nevertheless lead to a slight undesirable expansion of the star over the timescales of interest. The simple solution was to lower the pre-factor for the mesa EoS simulations models to
$\alpha_u = 0.1$
(see also González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022).
3. Results
Here, we present the orbital evolution, the characteristics of the ejected mass, and the state of the binary ejecta environment at the conclusion of the CE interaction. In Figure 1, we show an overview of the evolution for three of our simulations, 68IH, 85IH, and 100IH, from left to right respectively. We take four snapshots in time, shown vertically, where the first two show progressing RLOF, leading into the start and end of the dynamical inspiral in the latter two panels, respectively.
Cross sections of density in the orbital plane of the 68IH (left), 85IH (centre), and 100IH (right) simulations. Each column is a time sequence starting with two moments before the inspiral (top two rows), and ending with the start (
$t_i$
) and end (
$t_f$
) of the inspiral (bottom two rows). Each box is approximately 7 au in size.
3.1. Orbital evolution and unbound mass
Figure 2 shows the orbital evolution and bound mass of our systems as a function of time. We define the beginning and end of the inspiral in the same manner as was done by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019), that is, the upper and lower bounds of the criterion:
$|\frac{\dot{a}}{a}| \geq \frac{1}{15} \textrm{max}|\frac{\dot{a}}{a}|$
(circles and triangles in Figure 2), where the limiting value, here
$1/15$
, for the timescale is chosen by reasonable inspection. Due to the shallow inspiral of both the 150MH and 150IL simulations this criterion does not capture, even with modifications to the limiting value, satisfactory points that could be considered the start and end of the inspiral. Instead, we have chosen an alternative method for these simulations.
Top panel: binary core separation as a function of time for the twelve simulations (see Table 2). The circles and triangles are the start and end of the inspiral, respectively, as determined by the criterion
$|\frac{\dot{a}}{a}| \geq \frac{1}{15} \textrm{max}|\frac{\dot{a}}{a}|$
(Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019). Note that this criterion is not adopted for the
$q=1.5$
simulations, either due to a very shallow inspiral (150IL and 150MH), or the lack of inspiral (150IH; see text). Extrapolating from the time taken for the 100IL and 100IH simulations to inspiral, the computational cost for continuing the 150IH simulation is currently unfeasible. Bottom panel: the evolution of the bound mass for each simulation. Circles and triangles have the same meaning as in the upper panel, while the stars denote the time at which the resolution-dependent mass unbinding is estimated to start.
For the other simulations the average rate of orbital decline,
$\langle|\frac{\dot a}{a}|\rangle$
is less than 5% of the simulations
$\textrm{max}~|\frac{\dot{a}}{a}|$
. For the 150MH and 150IL simulations; however, we estimate it to be approximately 25%. Given its rate of decline is more linear than the other simulations shown in Figure 2, we find the start and end of the inspiral, respectively, as follows:
$a_{i} = a_0 - 0.25(a_0-\bar{a})$
and
$a_{f} = a_\textrm{end} + 0.25(a_0 - \bar{a})$
, where
$\bar{a}$
is the mean orbital separation over the whole simulation,
$a_0$
is the initial separation, and
$a_\textrm{end}$
is the orbital separation at the end of the simulation. An attempt to determine a universal algorithm that would select meaningful points for all our simulations proved fruitless, so we retain these somewhat arbitrary but quantitative criteria even if for some simulations they do not appear to select reasonable times by visual inspection.
A large value of q leads to a shallower inspiral and larger final separations (see
$\tau_\textrm{steep}$
and
$a_f$
in Table 2), with the 100IH and 100MH having
$a_f \sim 45$
–50 R
$_{\odot}$
. The 150IH simulation shows a decreasing orbital separation, but by 40 yr no sign of the dynamical inspiral is present, and we stopped this simulation due to excessive computational times. The 150MH simulation, aided by a slightly expanded stellar structure has a gentle inspiral (
$a_f\sim 150$
R
$_{\odot}$
) and plateau, even if the final separation is still declining at up to approximately 0.1 R
$_{\odot}$
/yr by the time we stop the simulation. This also true of the 150IL simulation.
Summary of parameters relating to the CE inspiral. Here
$a_0$
is the initial separation at Roche lobe contact. The beginning and end of the CE inspiral are found using the criterion from Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) and are denoted
$t_{i}$
and
$t_{f}$
, respectively, with their associated separations,
$a_{i}$
, and
$a_{f}$
. The parameters with the subscript ’steep’ refer to the time, separation, and timescale of the point of fastest inspiral in the interaction. The column
$t_{*}$
denotes the star in Figure 2, approximately the last moment before the resolution-dependent unbinding takes place.
$\dot{M}_{{{L}}_1, i}$
is the rate of mass transfer onto the companion one year after the start of the simulation.
$^{\#}$
This separation is at the end of the simulation.
In the bottom panel of Figure 2 we show the evolution of the bound mass for each simulation (defined as mechanical energy
$\lt 0$
). The circles and triangles represent the start and end of the CE inspiral as explained above, while the star symbol in the IL simulations represents the beginning of the resolution dependent unbinding, discussed below.
Variable amounts of pre-CE mass transfer occur before the rapid inspiral. This phase of early mass transfer typically increases in duration with increasing resolution, and with an increasing value of
$q=M_2/M_1$
(
$t_{i}$
in Table 2). Note this is not the case for the 100IL, and 150IL simulations due to their significantly smoother inspirals, for which
$t_{i}$
by our criteria begins to have difficulty. The duration of this early phase is not converged with resolution and is generally agreed to be substantially longer in reality (Iaconi et al. Reference Iaconi2017; Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019; González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). Simulations carried out with the tabulated mesa EoS also exhibit a shorter duration for this early phase of mass transfer compared to that of their ideal EoS counterparts, due to a greater degree of expansion of the star surface layers resulting from a reduced level of surface stability (González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022).
Higher resolution simulations unbind less mass overall than their lower resolution counterparts (as systematically observed across codes; e.g. González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). As expected due to the inclusion of recombination energy, the tabulated EoS simulations unbind most of the envelope. This behaviour has been reported before (e.g. Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019; González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). The additional unbinding seen immediately after inspiral’s conclusion completely disappears for high q simulations (present in 68IH but not in 85IH and 100IH), presumably because the envelope is expanded by the deposition of orbital energy but not unbound, which we will continue to investigate shortly.
The star symbols in Figure 2 denote the beginning of unbinding due to the density in proximity to the core dropping sufficiently that the SPH particle smoothing lengths become larger than their distance to the core point mass. This leads to additional unbinding already described by González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). This phenomenon is particularly pronounced in simulations with low resolution and an ideal EoS, clear in all the IL simulations (except 150IL) as a pronounced decrease in the bound mass after the star symbol in Figure 2. At high resolution, this does not occur (see Appendix A for further details).
In Figures 3 and 4, we plot the bound mass distribution as a function of time for the 68IH, 85IH, 100IH, and 150IH, and 68MH, 85MH, 100MH, and 150MH simulations, respectively. The vertical black lines correspond to the circle and triangle in Figure 2, the beginning and end of the inspiral. The total energy of each SPH particle was computed as the sum of its kinetic, core–gas potential, and gas–gas potential energies (making up the mechanical energy), adding the values for all particles, k, in each bin, i, as follows:
\begin{equation} \langle E_\textrm{tot}(r_i)\rangle = \frac{1}{N_i}\sum_{k\in i}\left[\frac{1}{2}m_{ p}v^2_{k}+U_{k}^\textrm{core-gas}+U_{k}^\textrm{gas-gas}\right],\end{equation}
where
$N_i$
is the number of particles per bin, and
$m_{p}$
is the mass of each SPH particle. The particles were binned radially into 5 R
$_{\odot}$
intervals, and the mean energy of the particles in each bin was rendered as colour. Note that we do not include thermal energy and only use the more stringent mechanical criterion. This is because the inclusion of thermal energy assumes that the entire thermal energy payload of the stellar envelope will be transformed into bulk kinetic energy. However, as has been shown in other work, this is not necessarily the case (e.g. Staff et al. Reference Staff2016; Iaconi et al. Reference Iaconi2017).
Only bound material is plotted (
$K+U\lt0$
) where we note there is very little bound material with energies larger than the minimum plotted energy of
$-10^{38}$
erg. Because each SPH particle has a constant mass, this plotted average energy per particle is equivalent to an average specific energy, since the division by mass introduces only a scaling factor.
Distribution of bound mass (
$K+U\lt0$
) throughout simulations 68IH (top left), 85IH (top right), 100IH (bottom left), and 150IH (bottom right). The pixels are binned at approximately 10 d in width, and 5 R
$_{\odot}$
in height, where we calculate the average energy of the gas within that radial bin, at that time step. Top panel: normalised orbital separation (blue) and the bound envelope (red). The vertical lines spanning the two plots denote, from left to right, the start (solid) and end (dashed) of the inspiral. These lines correspond respectively to the circle and triangle in Figure 2.
As in Figure 3, but for the 68MH, the 85MH, the 100MH, and the 150MH simulations.
An initial expansion of the stellar structure upon the start of RLOF also sees material at the surface becoming loosely bound. In the ideal gas EoS simulations the inspiral delivers orbital energy at the base of the envelope, which causes a wave of unbound material with a time delay that depends on the simulation (the trough between two peaks). Very little material is unbound in these high-resolution simulations because the unbound mass at the base of the envelope collides with bound gas above it.
The features in Figures 3 and 4 show just how complex the dynamics of the unbinding is. The deep troughs in the upper panels of Figure 3 are unbinding events, not too dissimilar to those displayed by the low-resolution simulations (Appendix A), showing that mass unbinding is not converged property but that, for an ideal gas EoS, less mass is unbound at higher resolution. For the tabulated EoS simulations, that are only carried out at higher resolution (Figure 4), the unbinding events are less pronounced as unbinding here happens more homogeneously due to the inclusion of recombination energy. For clarity, we note that the high potential energy seen as the horizontal red lines in Figures 3 and 4 represent the gas in close proximity of the core particles.
In Figure 5, we plot mechanical energies as slices in both the
$x-y$
plane (top) and
$x-z$
plane (bottom) for the 68IH simulation across four different time steps. These time steps are chosen to show the early, pre-inspiral phase of mass transfer, the start of the inspiral, and then two more time steps shortly after the conclusion of the inspiral. These latter two time steps show the unbinding ‘wave’ of material that is ejected from the region of the cores as the inspiral ends, which is similarly seen in the top right of Figure 3 as the spike of unbound material from the base of the envelope between the two peaks of bound material.
Slices of energy (
$E_\textrm{tot}$
, calculated as in Figures 3 and 4 but for both positive and negative energies) in the
$x-y$
plane (top) and the
$x-z$
plane (bottom) for the 68IH simulation. The selected times reflect the early mass transfer period (top left), the start of the inspiral (top right), whereas the bottom two panels depict the unbinding that occurs shortly after the inspiral concludes (as seen after the dashed line in Figure 3).
This wave of unbinding travels outwards from the cores and crosses the previously bound material at higher radii, but at the same time leaving bound material in its wake. Due to the primarily equatorial ejection throughout the pre-inspiral mass transfer phase, this later ejecta is also partially obstructed in the
$x-y$
plane, favouring polar ejection, as seen in the
$x-z$
slices of Figure 5. This unbound outflow through the poles was similarly seen in González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022).
3.2. The mass transfer rate through
$L_1$
Figure 6 shows the mass transfer rate through the inner Lagrange point,
$L_1$
. It begins at the start of the simulation, much earlier than the start of the inspiral and is calculated by counting the number of SPH particles that have passed from the Roche lobe of the donor into that of the accretor between each timestep. For this calculation the mass of the donor and accretor are the sum of their respective sink particle mass, plus any SPH gas particles within their respective Roche lobes. Note that mass transfer may not be conservative as mass may flow through
$L_2$
and
$L_3$
. We end this calculation at the moment of steepest inspiral.
Numerically-derived
$L_1$
mass transfer rates as a function of time for each simulation, where the low, high, and tabulated EoS simulations, are the dotted, solid, and dashed lines, respectively, in each panel. The calculation is then stopped at the point of steepest inspiral (
$t_\textrm{steep}$
in Table 2). Inserts zoom in on the first year of mass transfer.
The initial mass transfer rate (
$t \lesssim 1$
yr) appears erratic as it begins with a few particles crossing
$L_1$
, giving a minimum measurable mass transfer rate (Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019). We estimate that we can only calculate a reliable mass transfer rate for
$t\gt 1$
yr (Table 2). The mass transfer rate at 1 yr for all simulations is in the range
$\sim$
$0.5$
–
$1 \times 10^{-3}$
M
$_{\odot}$
yr
$^{-1}$
. The value of the mass transfer rate at the moment of steepest inspiral decreases by approximately an order of magnitude as the value of q increases, from 0.7 M
$_{\odot}$
yr
$^{-1}$
for
$q=0.68$
down to 0.08 M
$_{\odot}$
yr
$^{-1}$
for the
$q=1$
simulations.
The behaviour of the
$q=1.5$
simulations is distinctive in that the mass transfer rate remains relatively low with values of a few
$\times 10^{-3}$
M
$_{\odot}$
yr
$^{-1}$
. In simulations where
$q\gt1$
, mass is transferred from the least to the more massive binary component. In a conservative evolution, this leads to a widening of the orbit. The simulations 150IL/IH/MH show a modest inspiral likely caused by a combination of physical mechanisms, including mass loss through
$L_2$
and
$L_3$
, and tides.
Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) found (their Figure 7) that the analytical mass transfer rate prescription by Paczyński & Sienkiewicz (Reference Paczyński and Sienkiewicz1972) agreed with their derived values, but this conclusion critically depended on their estimated stellar radius. From 3D SPH simulations, it is difficult to calculate stellar radii accurately and even more so when the CE interaction brakes the donor’ spherical symmetry. It may be meaningless to carry out a comparison with analytical theory in this case but for more details of this problem see Appendix B.
Density slices in the
$x-y$
plane of simulations 68IL (left), 68IH (middle), and 68MH (right) at the beginning of the dynamical inspiral (circles in Figure 2). The red circle is centred on the centre of mass and passes through
$L_2$
. Significant mass ejection from behind the accretor (
$L_2$
, right) is accompanied by a slightly less pronounced ejection from behind the donor (
$L_3$
, left). At high mass ratios such as those used in this work, the difference between
$L_2$
and
$L_3$
is small. Plot was generated using Splash (Price Reference Price2007).
3.3. Pre-CE mass loss through
$L_2$
and
$L_3$
Mass loss through
$L_2$
and
$L_3$
can be responsible for the formation of a circumbinary disc and for the reduction of the orbital separation seen during the dynamical inspiral. To estimate the total mass loss from the binary before the onset of the inspiral, we first determine the distance
$D_{L_2}$
between the primary star and
$L_2$
using the formulation of Misra et al. (Reference Misra, Fragos, Tauris, Zapartas and Aguilera-Dena2020):
\begin{equation} \frac{D_{L_2}}{R_{L}} = \begin{cases} 3.334 q^{0.514}e^{-0.052q} + 1.308 & \text{for $q \lt 1$},\\[2pt] -0.040q^{0.866}e^{-0.040q} + 1.883 & \text{for $q \geq 1$},\\ \end{cases} \end{equation}
where
$R_{L}$
is the Roche lobe radius of the donor given by Eggleton (Reference Eggleton1983) formula
\begin{equation} {\frac{R_{{L}}}{a}} = \frac{0.49q^{2/3}}{0.6q^{2/3} + \textrm{ln}(1+q^{1/3})}.\end{equation}
We then select all SPH particles farther than
$L_2$
from the centre of mass (
$r_{{{L}}_2}$
) before the time of inspiral,
$t_i$
, as denoted by the circles in Figure 2. This pre-inspiral mass loss occurs primarily through
$L_2$
(as also seen for varying q values by Scherbak et al. Reference Scherbak, Lu and Fuller2025), but will also increasingly include material ejected through
$L_3$
as we approach the onset of the dynamical inspiral. As such, using
$r_{L_2}$
captures ejected material through both
$L_2$
and
$L_3$
, providing an upper limit to the total mass lost from the binary. Circles with radius
$r_{{{L}}_2}$
are shown on density slices in Figure 7.
It is possible that the onset of
$L_3$
mass outflow may also exert a torque on the binary that results in an additional drain of angular momentum, displayed as a small increase in eccentricity traced by the wiggles on the separation curve present in all simulations (particularly obvious at high resolution in Figure 2) prior to inspiral.
Table 3 presents a summary of the bound and unbound gas mass for the material ejected through
$L_2$
and
$L_3$
. Higher resolution, ideal gas EoS (IH vs. IL simulations) increases slightly the amount of mass outside of
$L_2$
by the start of the inspiral (
$t_i$
), as seen in Column 4, increases dramatically the fraction of that mass that is unbound (three-quarters of the mass is unbound vs. half for lower resolution; Column 6), and results in approximately the same mass remaining unbound by the end of the simulations (
$t_f$
) – meaning no more of the material outside
$L_2$
at the start of inspiral is unbound by the end of it (Column 9).
Data describing mass lost through
$L_2$
. Columns are as follows: (2):
$r_{{{L}}_2}$
– distance of
$L_2$
from the centre of mass at the onset of the dynamical inspiral, (3):
$t_{i}$
; (4):
$M_{\gt L_2, {i}}$
– percentage of the envelope mass outside
$r_{{{L}}_2}$
at time
$t_{ i}$
; (5)
$M_{\gt L_2, \textrm{UB}, i}$
– unbound mass outside of
$r_{{{L}}_2}$
at time
$t_{i}$
as a percentage of envelope mass, and of the mass exterior to
$L_2$
(6), respectively; (7):
$t_{f}$
– end of the CE inspiral (triangles in Figure 2);
$M_{\gt L_2,\textrm{UB}, f}$
– mass outside
$r_{{{L}}_2}$
at
$t_i$
that is unbound at
$t_f$
, as a percentage of envelope mass (8), and as a percentage of mass exterior to
$L_2$
(9); (10):
$M_{\textrm{Tot, UB}, i}$
, and (11)
$M_{\textrm{Tot,UB},f}$
– total unbound mass in the simulation at
$t_i$
and
$t_f$
, respectively (triangles in Figures 2); (12)
$M_\textrm{Tot,UB,*}$
– total envelope mass unbound at the end of the simulation.
$M_\textrm{env} = 0.49$
M
$_{\odot}$
.
$^*$
Did not go through a CE inspiral, percentages are instead taken at the end of the simulation.
$^\#$
Percentage is instead taken at the time denoted by the stars in Figure 2.
The tabulated EoS (MH vs. IH simulations) has distinctly less mass outside of
$L_2$
by
$t_i$
than the ideal gas EoS (but this could be a timing issue), the fraction of that mass that is unbound is higher (three quarters vs. half), but by the end of the inspiral the entirety of the mass outside
$L_2$
has been unbound – an effect of the recombination energy delivery, as expected.
Finally, looking at the dependency on q for the IH and MH simulations only: the mass outside of
$L_2$
by
$t_i$
slightly increases along the
$q=0.68, 0.85, 1.00$
sequence, between 10–20% and 20–30%, and the fraction of that mass that is unbound at
$t_i$
increases somewhat between 50–75% and 60–80%. We notice, however that these increases are between
$q=0.68$
and 0.85, while there is no real change for the
$q=1.00$
simulation compared to the
$q=0.85$
one. These effects balance out: the more massive
$L_2$
outflows for high q values, are overall less bound, making the potential circumbinary disc mass similar.
The extreme
$q = 1.5$
simulations, with
$L_2$
appearing now on the outside of the donor (primary) not accretor (secondary) star, seem to exhibit their own unique behaviour. These simulations regardless of EoS or resolution are extremely efficient at unbinding almost all of the 5–10% of the envelope they do eject. Even simulation 150IH shows these trends, despite the fact that it has progressed only 23 yr in total and has not gone through the inspiral at all. Beyond unity, simulations ultimately eject less mass from the binary, and since virtually all of this material is unbound, we find circumbinary disc formation unlikely.
In columns 10, 11, and 12, we present the total amount of envelope mass unbound at the start of the inspiral, the end of the inspiral, and the end of the simulation, respectively. Again we find that higher resolution unbinds less mass (Column 11) and the tabulated EoS unbinds generally more (Columns 11 and 12). Low-resolution simulations display a strong resolution-dependent unbinding pattern, as discussed previously (González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). Interestingly for the 150MH simulations the fraction of unbound mass at the end of the simulation is only
$\sim$
80%, although it is decreasing. Fundamentally this information indicates that retaining a fraction of disc mass from
$L_2$
ejecta is unlikely.
3.4 Kinematic properties of mass lost through
$L_2$
and
$L_3$
Here, we present a short analysis of the angular momentum that is ejected through
$L_2$
and
$L_3$
before the inspiral. Caution is needed in carrying out comparisons between simulations, due to the arbitrary definition of the time of inspiral start. In particular the definition for the 150IL/IH/MH simulations is substantially different, as is the nature of the inspiral.
Similar to what was found in Section 3.3, each simulation ejects around 4–8% of the total binary’s mass through
$L_2$
by the start of the inspiral (Table 4, column 2; where we express percentages of total mass instead of envelope mass as was the case in Table 3). Of the mass lost through
$L_2$
and
$L_3$
, about half is unbound in the IH simulations, and almost all is unbound in the MH simulations, as already discussed in Section 3.3. In column 3 of Table 4 we calculate the total angular momentum of the mass that is lost from
$L_2$
. The IH simulations eject
$\sim$
30% of the binary’s total angular momentum leading up to the inspiral, with the IL simulations ejecting somewhat less and likely demonstrating a lack of convergence. IL simulations, but not IH, exhibit a trend of decreasing fraction of ejected angular momentum for larger values of q. The MH simulations eject a smaller,
$\sim$
20%, with the 150MH showing a much larger value of 55%, the inverse trend to the L simulations (the 150IH cannot be considered due to the short time of the simulation).
Properties of the
$L_2/L_3$
ejecta material outside
$r_{L_2}$
at the start of inspiral. Columns 2 and 6 are taken from columns 4 and 5 of Table 3, but now show the percentage of the total binary mass (including the core of the primary and companion) outside
$L_2$
at the start of the inspiral; columns 3 and 7 are the angular momentum outside
$L_2$
at this time, as well as the angular momentum that is outside
$L_2$
and unbound, as a percentage of total angular momentum of the system. The parameter
$\gamma_\textrm{loss}$
is defined in Nelemans et al. (Reference Nelemans, Verbunt, Yungelson and Portegies Zwart2000) and shows the ratio between the specific angular momentum of the material lost from the binary, and the initial specific angular momentum of the binary. We also calculate
$\gamma_{L_2} = h_{L_2}/h_\textrm{bin}$
, where
$h_{L_2}$
is the initial specific angular momentum of
$L_2$
. In the final two columns we provide the binary’s specific and total angular momenta, respectively.
$^*$
Did not go through a CE inspiral, value is instead taken at the end of the simulation.
A convenient way to express the amount of angular momentum that is lost (
$\Delta J$
) from a binary due to mass loss (
$\Delta M$
) is with the parameter
$\gamma_\textrm{loss}$
(Nelemans et al. Reference Nelemans, Verbunt, Yungelson and Portegies Zwart2000), which can be defined as:
\begin{equation} \frac{\Delta J}{J_\textrm{bin}} = \gamma_\textrm{loss}\frac{\Delta M}{M_\textrm{tot}}.\end{equation}
The quantity
$\gamma_\textrm{loss}$
quantifies the impact of mass loss on the evolution of the orbit. Larger values of
$\gamma_\textrm{loss}$
imply each percentage of mass lost carries with it a higher percentage of the binary’s angular momentum, thus eliciting a more significant change in the binary’s orbit.
In Table 4 (Column 4) we list
$\gamma_\textrm{loss}$
(calculated with values with subscript ‘i’, hence at the start of in-spiral) to be in the range of 3–4 for most of our simulations. This indicates that the escaping material carries away 3–4 times the mean specific orbital angular momentum of the binary, implying relatively efficient removal of angular momentum for a given amount of mass loss. We also see a slight decrease of
$\gamma_\textrm{loss}$
as mass ratio increases for MH simulations, though the IH simulations have a pretty constant value. This may imply that, for a given amount of mass loss, higher mass ratio binaries remove angular momentum slightly less efficiently, experiencing slower orbital evolution. Our extreme
$q=1.50$
simulations tell a complex story. Both ideal gas simulations have considerably smaller
$\gamma_\textrm{loss}=2$
–3.
Using binary evolution reconstruction techniques based on three observed double white dwarf systems, Nelemans et al. (Reference Nelemans, Verbunt, Yungelson and Portegies Zwart2000) determined that for the first of two interactions that generated those systems today, a stable mass transfer phase between a
$\sim$
2 M
$_{\odot}$
giant and a similar mass main sequence companion (with no CE inspiral) was characterised by
$\gamma_\textrm{loss} \sim 1.7$
. Larger values of around
$\gamma_\textrm{loss} \approx 3$
–7 (modulated by orbital evolution – a decline from
$\sim$
7 to
$\sim$
3 for the initial slow orbital decline, followed by an increase back to
$\sim$
7 during the fast inspiral) were instead determined by MacLeod et al. (Reference MacLeod, Ostriker and Stone2018b) for three binary simulations with
$q=0.3$
. MacLeod & Loeb (Reference MacLeod and Loeb2020) calculated overall larger values, with simulations with
$q=0.03$
–
$0.3$
yielding values between
$\sim$
21 and
$\sim$
6 (lower
$\gamma_\textrm{loss}$
for larger q, similar to what we observe).
We also calculate
$\gamma_{L_2} = h_{L_2}/h_\textrm{bin}$
, where
$h_{L_2}$
is the specific angular momentum of
$L_2$
at the start of the simulation when it is at its maximum and
$h_\textrm{bin}$
is the total specific angular momentum of the binary at time zero (Table 4, Column 8). We find that in all cases,
$\gamma_\textrm{loss} \lt\gamma_{L_2}$
, consistent with both MacLeod et al. (Reference MacLeod, Ostriker and Stone2018b) and MacLeod & Loeb (Reference MacLeod and Loeb2020). We note that
$\gamma_\textrm{loss}$
need not be constant, as shown by MacLeod & Loeb (Reference MacLeod and Loeb2020), while we have calculated what is effectively an average value over the slow in-spiral before the CE fast inspiral takes place.
Figure 8 shows the velocity of each SPH particle as a function of distance from the centre of mass for the 68IH (left panels) and 68MH (right panels) simulations at the start (top row) and end (middle row) of the inspiral, and at the final timestep of the simulation (bottom row). The particles ejected from
$L_2$
and
$L_3$
before the onset of the dynamical inspiral are contoured in red and correspond to those outside the red circle in Figure 7. The black curve in these plots represents the escape velocity, separating bound from unbound gas.
Velocity of each SPH particle as a function of distance from the binary’s centre of mass for the 68IH (left) and 68MH (right) simulations. The times shown in chronological order from top to bottom are the start of the inspiral, the end of the inspiral, and the last timestep of the simulation. For simplicity the black line is an approximate escape velocity that assumes the central mass is the primary star and the companion – providing an upper limit for bound material. The particles within the red contour are located outside the radius of
$L_2$
at the onset of the inspiral as shown is Figure 7. To give an indication of how the material is distributed we construct a 2D histogram of mass with 300 bins in each axis, where the mass shown is the mass per bin. We have also marked the core and companion particles in orange.
The material with (
$v\gt v_\mathrm{esc}$
) that is unbound at small radii at the time of the inspiral’s onset (top row), collides with, and pushes out the material above it. This causes some of the inner ejecta to slow down and possibly become bound once again. In doing so it is also accelerating the outer ejecta, where a peak of unbound gas develops at the boundary between the
$L_2$
ejecta and material that follows after, ejected during the inspiral.
Unbound material expands nearly homologously (the
$L_2$
outflow is arranged in a straight line in a linear-linear plot). The MH simulations unbind almost all the
$L_2/L_3$
ejecta, while the ideal gas EoS simulations do not but the material is swept by the ejecta that comes later, following the inspiral. As such, the material that was ejected through the outer Lagrange points does not have the opportunity to form a disc. We conclude that disc formation via
$L_2$
and
$L_3$
ejection is unlikely if a CE ejection takes place. Next we investigate whether a fall-back disc is possible.
3.5. A fall-back disc
When the orbit stabilises after the interaction, the fraction of the envelope that remains bound may fall back. Depending on how much envelope falls back onto the binary, and where it falls back to, it may cause a second interaction (as simulated by Kuruwita et al. Reference Kuruwita, Staff and De Marco2016) that could result in a further change of the binary orbit. If, however, the material has sufficient angular momentum, it could form a circumbinary disc instead.
We start by assuming that the ejected, but bound material is in a ballistic orbit, such that both angular momentum and orbital energy are conserved. This material will then reach the pericentre of its orbit within an orbital period. We take the fall back time for this initial phase to be half of an orbital period, on average:
\begin{equation} t_\textrm{fb} = \pi\sqrt{\frac{a^3}{GM}}, \end{equation}
where M is the total central mass (primary core and companion), G is the gravitational constant, and a is the semi-major axis of the orbit (derived from each SPH particle velocity vector, distance from the centre of mass and central mass). The fall back radius coincides with the circularisation radius and is given by
\begin{equation} R_{c} = \frac{h^2}{GM}, \end{equation}
h is the specific angular momentum.
As gas streams back towards the centre of mass and approaches the fall-back radius, collisions between particles generate shocks and dissipation of kinetic energy. The orbit circularises, but angular momentum is conserved. Circularisation timescales are of the order of the orbital period at this new semi-major axis given by
$R_{c}$
\begin{equation}\tau_\textrm{circ} \sim 2\pi \sqrt{\frac{R_c^{3}}{GM} }.\end{equation}
Eventually the system will further evolve through viscous interactions, which will lead to the disc spreading – some material falling into the potential well and accreting to the binary while some material moves outwards. The associated timescale is
\begin{equation} \tau_\textrm{visc}=\frac{R_{c}^2}{\nu}, \end{equation}
where
$\nu$
is the gas viscosity which can be written as
\begin{equation} \nu = \alpha c_{s}^2\sqrt{\frac{R_{c}^3}{GM}}, \end{equation}
where
$c_s$
is the sound speed and
$\alpha$
the disc viscosity parameter (Shakura & Sunyaev Reference Shakura and Sunyaev1973). We now use these estimates of fall-back radius,
$R_{c}$
(henceforth
$R_\textrm{fb}$
), and timescale,
$t_\textrm{fb}$
, to approximate the size and distribution of any post-interaction disc that might form. For the mass of this disc, we take the bound mass present at the end of the simulation.
Ideal gas simulations likely underestimate the amount of unbound mass (and overestimate the amount of bound mass), while tabulated EoS simulations would do the opposite since all liberated recombination energy is allowed to do work and none escapes. In this way IH and MH simulations likely bracket the available mass to make a fall-back disc.
To calculate the fall-back timescale given by Equation (5), we take a to be the semi-major axis of the orbit for each SPH particle, calculated from the specific orbital energy,
$a = -GM/2E$
. Here we calculate E in a similar manner to what was done previously, by summing the specific mechanical energies for the bound material
$E = 0.5v_{p}^2 + U_{p, \textrm{core}} + U_{p,p} \lt 0$
, where
$v_{p}$
is the velocity of the SPH particle, and
$U_{p, {\rm core}}$
and
$U_{p, p}$
are the specific potential energies of the particle due to both the core particles and the other SPH gas particles. In doing so the current motion of each bound SPH particle is considered, as opposed to just its distance from the centre of mass. For our 68IH/MH simulations, we plot
$t_\textrm{fb}$
as a function of current location from the centre of mass in Figure 9. In this plot we also logarithmically bin the mass of each particle according to these axes, thus giving an indication of how the mass is distributed at this point in time when we might expect material to begin falling back.
Fall-back time as a function of distance from the centre of mass for our 68IH (top), and 68MH (bottom) simulations. The colour bar indicates the amount of mass within each logarithmically sized bin (pixel), to show the distribution of bound ejecta at
$t=t_\textrm{end}$
. To calculate the fall back time we take half the orbital period, where the semi-major axis, a, is derived from the orbital energy of the gas particle.
As shown for simulations 68IH and 68MH, we find that the mass available to fall back is initially distributed up to a distance of a few
$\times 1\,000$
R
$_{\odot}$
at the end of the simulation. Material with longer fall-back times can be seen at larger radii and corresponds to gas that is moving away from the binary at speeds that are close to escape velocity. Although not shown here, a similar distribution of fall back material is found for all simulations.
Figure 10 shows the distribution of the bound gas with respect to its fall-back radius, as given by Equation (6). The blue lines are for the IH simulations, that unbind only part of the gas, leading to more massive discs, while the orange lines are for the MH simulations that unbind most of the gas leading to a significantly less massive disc. We find that the distribution of fall-back material is narrower for the tabulated EoS simulations and depends somewhat on the mass ratio of the binary, with increasing q leading to slightly larger inner disc’s edge. This is due to simulations with larger values of q having more angular momentum because they start mass transfer when the companion is farther out. Our simulation with the lowest mass ratio (68IH/MH) results in a fall-back disc with a radial span of 10–3 000 R
$_{\odot}$
and 20–750 R
$_{\odot}$
for the ideal gas and tabulated EoS simulations, respectively. The secondary, smaller peak near 3–10 R
$_{\odot}$
, is due to the gas that effectively remains bound around the stellar cores (seen in Figure 9 as the two upside-down triangle shapes at small values of r). For the 85 and 100 simulations similar values are observed (15–4 500 R
$_{\odot}$
and 35–750 R
$_{\odot}$
for the 85IH and 85MH simulations, respectively and 20–4 500 R
$_{\odot}$
and 40–750 R
$_{\odot}$
for the 100IH and 100MH simulations respectively.)
Histograms of bound mass as a function of fall-back radii for
$q= 0.68$
(top),
$q= 0.85$
(middle), and
$q= 1$
(bottom), where the blue and orange lines are the high and tabulated EoS simulations for each value of q, respectively. Plots are generated at
$t=t_\textrm{end}$
. The circles and diamonds in each plot represent the orbital distance from the centre of mass of the primary and companion cores, respectively.
Using this range of fall-back radii we now calculate an approximate viscous timescale from Equation (8). The temperature of the gas surrounding our post-interaction binary is in the range
$T \approx 4\times10^3 - 1\times10^4$
K for all simulations. This then corresponds to a sound speed in the range of
$c_{s} = 5$
–12 km s
$^{-1}$
. If we take
$\alpha = 0.01$
–
$0.1$
, then the viscous timescale, and thus the expected circularisation timescale of our circumbinary disc is
$\tau_\textrm{visc}\approx10$
–100 yr.
Table 5 summarises the amount of mass expected to fall back onto the binary (the bound mass at the end of the simulations) for our simulations. In the case of the tabulated EoS simulations, fall-back masses are in line with the observed masses of post-giant circumbinary discs (0.02–0.05 M
$_{\odot}$
; Pourmand et al. Reference Pourmand, Kamath, De Marco and Wardle2025). How much more massive these discs may have been at the time of formation is unclear, but with the observed mass range in mind, we can state that only up to
$\sim$
5% of our envelope would need to remain bound. In Table 5 we also list the average location,
$R_\textrm{fb} \sim R_{c}$
, of this mass from the centre of mass of the system (binary plus bound gas), and the distances from the centre of mass of each of the two stars in the binary,
$d_\textrm{orb,1/2}$
.
Mass-weighted average fall-back times (Equation 5) and total mass of the bound envelope at end of inspiral (
$t=t_{\textrm{end}}$
) for each simulation. The total mass of the fall-back material is calculated from the area underneath the curve for the regions exterior to the orbit, that is, the outer set of peaks for each simulation in Figure 10.
$R_{\textrm{fb}}$
is the disc radius from the centre of mass, (the peak of each histogram in Figure 10, note for the peak plateau in the 68IH simulation we have marked the inner and outer radii), and the distance from the centre of mass of the core and companion, respectively (
$d_{\textrm{orb,1/2}}$
), is given to show the location of the disc with respect to the orbit of the binary (the circles and diamonds in Figure 10, respectively).
4. Discussion
In this section we discuss the impact of the mass ratio on CE interactions between RGB stars and compact companions. The CE interaction likely plays an important role in the evolution of binaries with periods less than approximately 2 000 d (Krynski et al. Reference Krynski, Siess, Jorissen and Davis2025). As such we discuss the weakening of the CE interaction in the context of post-RGB (and post-AGB) binaries. In particular we consider the final orbital separation, and whether a stable circimbinary disc can form.
4.1. Mass transfer stability and final orbital separation
Our simulations demonstrate that the duration of the stable mass transfer phase preceding the CE inspiral is sensitive to the mass ratio
$q = M_2/M_1$
, the resolution, and the equation of state (EoS). For example, ideal gas EoS simulations 68IH, 85IH, and 100IH, with increasing q see the pre-inspiral time go from 39 to 64 to 72 yr (shorter, but with the same trend for the tabulated EoS simulations; see
$t_i$
in Table 2). Noting that the pre-inspiral time is not converged, such that this value must be a lower limit, our simulations support the analytic prediction (Tout Reference Tout1991; Frank et al. Reference Frank, King and Raine2002; Dan et al. Reference Dan, Rosswog, Guillochon and Ramirez-Ruiz2011) that mass transfer begins to stabilise for
$q\gtrsim1$
.
Mass transfer rate through
$L_1$
early in the simulation has a slight correlation with q (
$\dot{M}_{L_1,i} = 4\ \textrm{to}\ 9 \times 10^{-4}$
M
$_{\odot}$
yr
$^{-1}$
for the ideal gas and tabulated EoS simulations, similarly), but this is not a converged quantity, with the low-resolution simulations showing values that are about twice those of the respective high-resolution simulations. The mass transfer rate remains approximately constant during the pre-inspiral time for the tabulated EoS simulations. A longer the pre-inspiral mass loss phase achieved by increasing resolution may not increase the total amount of mass transferred as the mass transfer rate also decreases with increased resolution. It is therefore possible that the lack of convergence of the length of time before inspiral may not change the amount of total mass transferred before inspiral.
The final separation depends on the mass ratio with values of 27, 38, and 49 R
$_{\odot}$
for the high resolution, ideal gas EoS 68IH, 85IH, and 100IH simulations at the end of the inspiral, with slightly decreasing values afterwards. The final separation after the inspiral is converged, and it is also very similar for the different EoS. Contributing factors to this trend are both the higher initial angular momentum of the higher q simulations (due to the larger initial orbital separations and mass) and the larger orbital energy inherent to the more massive companions (as observed by Passy et al. Reference Passy2012, who started different q simulations at the same initial orbital separation). This said, the largest value, at
$\sim$
50 R
$_{\odot}$
(or a period of
$\sim$
40 d) is still smaller than the smallest separation of the post-RGB/AGB binaries (
$\sim$
100 R
$_{\odot}$
; or a period of
$\sim$
100 d; Kluska et al. Reference Kluska2022), showing that entering the CE will tend to cause too much orbital shrinkage to explain the observed separations, no matter what the mass ratio.
The
$q = 1.5$
is not representative of post-RGB binaries because the companion could not be so much more massive than the primary and still be on the main sequence. (These massive companions could, however, represent massive WDs). Our 150IL/IH/MH simulations allowed us to investigate an even more extreme value of q. Our ideal gas simulation, 150IH, shows limited inspiral within the limited simulation runtime, and the orbital separation begins to plateau along with the mass transfer rate itself. The tabulated EoS simulation, 150MH has a modest, shallow inspiral, leading to a final separation at the end of the simulation of
$\sim$
140 R
$_{\odot}$
(a period of
$\sim$
130 d) at the low end of the observed separations (plateauing but still decreasing gently). This confirms the trend and suggests that, in real systems where mass transfer proceeds more gradually than in simulations, higher q values likely allow the binary to avoid inspiral and achieve wider post-interaction separations that are more inline with observed post-giant systems.
With increasingly stable mass transfer enabled by larger values of q, heat loss from the surface of the star may be sufficient to stabilise and slow down the mass transfer rate and possibly avoid the inspiral altogether. The stellar envelope could be removed in this way. Alternatively, the binary could widen and mass transfer could stop altogether. In such case, the RGB star envelope removal would eventually complete by regular winds, leaving a post-RGB binary systems with uncertain orbital parameters.
4.2. Unbinding the envelope
Higher resolution simulations typically unbind less mass than their lower-resolution counterparts unless aided by a tabulated EoS that includes recombination energy. This agrees with findings by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Chamandy and Price2020) and González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022) and emphasises the need for caution when interpreting the total unbound mass in simulations. While the use of a tabulated EoS has become common practice in CE simulations, the retention of all energy by our adiabatic simulations means that the unbinding observed with tabulated EoS must be considered an upper limit. Hence, when discussing the circumbinary disc mass, the ideal gas simulation, which unbind less mass, would generate a more massive disc, while the tabulated EoS simulations, which unbind almost the entire envelope, would generate a much lower mass disc.
Our 150MH simulation illustrates that nearly the entire envelope can be unbound without a deep dynamical inspiral. This suggests that, in physical systems, a slow process of unbinding enabled by recombination energy – combined with more stable mass transfer – could be sufficient to eject the envelope.
4.3. Post-interaction circumbinary disc formation
There are at least two ways to form a circumbinary disc: via
$L_2$
and
$L_3$
mass outflow and by fallback of bound mass after a CE ejection. In the current context of CE interactions, the former method works if the CE ejection that follows the formation of a circumbinary disc via
$L_2/L_3$
does not sweep the disc away. The latter method works if there is sufficient bound mass, with sufficient angular momentum to form a disc around the binary.
We observe almost no increase in the mass outside
$L_2$
by the start of the inspiral as a function of q for all simulations. Of the mass outside
$L_2$
, only about half is still bound at that time for all simulations. That is the mass that would be destined to form a disc. The 150IH/MH simulations shows that pushing q to higher values increases the amount of mass lost through
$L_2$
only for the MH simulation, but for all 150 simulations it increases the fraction of that mass that is unbound, leaving almost no gas to make a disc. This, once again, argues against the hypothesis by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) that higher q simulations could form a more massive (bound) disc.
This said, some gas (whether ejected via
$L_2$
before the inspiral, or ejected later on via the CE ejection) remains bound to the system (more so for the ideal gas EoS simulations and less for the tabulated EoS ones, as explained in Section 4.2). By modelling bound gas as ballistic particles, we find fall-back timescales of centuries to about one millennium. The fall-back/circularisation radii range from 0.5 to about 1.5 au, large enough to be stable around the binary, but much smaller than the discs observed around post-RGB/AGB binaries. It is possible that viscous spreading will take place on relatively short timescales thereafter, resulting in larger discs.
4.4. Circumbinary material and orbital elements
The fall back of material that remains bound to the system post-inspiral presents a potential opportunity for secondary interactions that may alter the orbital elements. In Section 4.3 we discussed how the fall back of material in these simulations should form a disc outside the orbit of our binaries. This said, some ejected material may still return close enough to the orbit of the binary such that it may re-accrete before the disc formation.Footnote b
Even if some gas had so little angular momentum that it would fall back onto the orbit, it is unlikely that an interaction between relatively little gas and the binary would lead to a significant decrease of the binary’s orbital separation: at orbital separations of
$\sim$
50 R
$_{\odot}$
, the energy required to shrink the orbit by half, down to 25 R
$_{\odot}$
, is approximately the same amount of energy it takes to shrink the orbit from the original
$\sim$
200 R
$_{\odot}$
to
$\sim$
50 R
$_{\odot}$
. This said, any gas that falls back to within the binary orbit should continue to interact with the binary through both accretion and resonances, until it is redistributed leaving a cavity of 2–3 times the orbital separation (Artymowicz & Lubow Reference Artymowicz and Lubow1994). We will consider this phase of disc adjustment in a future work.
Gagnier & Pejcha (Reference Gagnier and Pejcha2023, Reference Gagnier and Pejcha2024, Reference Gagnier and Pejcha2025) have investigated the post interaction environment using 3D hydrodynamic simulations. Using a setup similar to Morris & Podsiadlowski (Reference Morris and Podsiadlowski2006, Reference Morris and Podsiadlowski2007, Reference Morris and Podsiadlowski2009) and Hirai et al. (Reference Hirai, Podsiadlowski, Owocki, Schneider and Smith2021) they initialised their systems with a distribution of circumbinary material endowed with angular momentum consistent with that expected from material spun up during the interaction. They found that the material forms a thick, bound, out-flowing disc in the equatorial plane, with inward falling material around the poles. This orbiting material then grows eccentric with local over-densities present despite the central binary’s orbit being held circular (Gagnier & Pejcha Reference Gagnier and Pejcha2023). It is speculated this eccentric material could influence the eccentricity of the binary, which is interesting since many observed post-RGB/AGB binaries feature non-zero orbital eccentricities. We suspect, however, that tides, particularly after the CE, are much too efficient at circularising the orbit, such that binary-gas interactions are unlikely to meaningfully change this. This does still leave room for the decay of the orbital separation even after the inspiral, which Gagnier & Pejcha (Reference Gagnier and Pejcha2023) found may be possible on timescales of
$10^3$
–
$10^4 P_\textrm{orb}$
in the presence of accretion. The presence of magnetic fields also seems to have little effect on the future binary evolution, though can act to modestly shape the resulting disc-like structure (Gagnier & Pejcha Reference Gagnier and Pejcha2024).
Previous simulations (Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019; González-Bolívar et al. Reference González-Bolvar, De Marco, Lau, Hirai and Price2022) similar to ours have also shown that material surrounding the binary is typically dragged into co-rotation which is in part responsible for halting the inspiral due to a quenching of net torques. Recently Gagnier & Pejcha (Reference Gagnier and Pejcha2025) similarly found that any further evolution of the binary in the post-interaction environment is likely due to interactions between the inner, co-rotating layers around the binary and those outside this region that are not co-rotating. However, as previously mentioned, the resulting influence on the binary’s orbit is expected to be minimal.
5. Conclusions and summary
We performed a suite of 3D hydrodynamic simulations of CE evolution between a 0.88 M
$_{\odot}$
RGB star and compact companions with varying mass ratios,
$q \equiv M_2/M_1 = 0.68$
–
$1.5$
and two equations of state, to investigate regimes of increasing mass-transfer stability and their impact on the CE interaction outcome. We find that:
-
• Increasing q extends the duration of the pre-inspiral mass transfer phase, weakens the inspiral, and leads to wider final orbital separations. Some of these quantities, with the exception of final separation, are not converged with respect to simulation resolution. However, increasing the resolution meaningfully is excessively expensive computationally.
-
• We expect that binary interactions with high q values would be even more stable in nature, with even longer pre-inspiral timescales, although the final separation would be similar to the one obtained here.
-
• For sufficiently high q, the classical CE inspiral might be entirely avoided, suggesting that the envelope of the donor could be unbound through sustained, stable mass loss rather than dynamical orbital decay.
-
• Final orbital separations, even for large q values are of the order of 50–100 R
$_{\odot}$
, or periods of 30–120 d, still generally smaller than the range of semi-major axes and periods of post-RGB and post-AGB binaries. -
• We do not witness the formation of a post-CE, circumbinary disc through mass ejection via
$L_2$
. Tabulated equation of state simulations unbind all the
$L_2$
ejecta, while ideal gas simulations result in only about half of the
$L_2$
ejecta remaining bound to the system. In both cases, however, the mass ejected during and after the inspiral, tends to sweep away any low mass disc that may have formed just before. -
• In all simulations a fraction of the envelope remains bound to the system and we calculate that it will fall back towards the binary forming a disc that spans 15–4 000 R
$_{\odot}$
(0.1–20 au) for ideal gas simulations, and 25–750 R
$_{\odot}$
(0.1–3.5 au), for tabulated EoS simulations, in a time frame between a tenth of a year and 1 000 yr, with the bulk of the bound mass returning within a year. A further 10–100 yr may be needed to viscously spread the disc further. Such discs may develop characteristics that are in line to those surrounding post-RGB/AGB binaries. -
• Our findings suggest that high mass ratios at the time of the interaction may help to provide ideal conditions for a weakened CE that results in wider binaries with circumbinary discs; however, this alone is unlikely to be the solution to the formation of such observed systems. In a second paper, we will investigate the effect of mass ratio on AGB, instead of RGB stars which, with their lighter envelope may exhibit even weaker interactions.
-
• Finally, we start to see a behaviour that is critically different from a CE inspiral, only for our
$q = 1.5$
simulation. This ratio may be reduced if we could simulate a more realistic stellar surface, where rapid cooling may stabilise the mass transfer and avoid the CE for lower values of q.
Acknowledgements
This work was supported by resources awarded under Astronomy Australia Ltd’s ASTAC merit allocation scheme on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR programme receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government. This research was supported by the Commonwealth through an Australian Government Research Training Program Scholarship [DOI: https://doi.org/10.82133/C42F-K220]. LS is research director at the F.R.S-FNRS. We acknowledge discussions with Mike Lau, Taissa Danilovich, Kayla Martin, Ana Juarez, Stephen Neilson, Chunliang Mu, and Ali Pourmand.
Appendix A. Resolution-dependent unbinding for low-resolution simulations
In Figure A1 we show the distribution of bound gas in ideal gas EoS simulations demonstrating that the low-resolution simulations unbind gas at the base of the envelope (the white zone that develops after the dotted vertical lines in three of the panels in Figure A1), a behaviour typical of the resolution-dependent unbinding discussed by González-Bolívar et al. (Reference González-Bolvar, De Marco, Lau, Hirai and Price2022). This is not observed in high-resolution simulations (see Figure 3). Interestingly, simulation 150IL does not display the wave of unbinding at the base of the envelope. By the time we stop that simulation, there has been only little, gentle inspiral, which however has unbound 73% of the envelope mass (Table 3). Yet the modality of this unbound mass is clearly distinct from that of the other three simulations in Figure A1.
Distribution of bound gas for low-resolution simulations 68L (top left), 85IL (top right), 100IL (bottom left), and 150IL (bottom right). A demonstration of resolution-dependent gas unbinding at the base of the envelope is present in all simulations that undergo a strong inspiral.
Appendix B. On the comparison of numerical and analytical mass transfer rates
Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) carried out a comparison between the mass transfer rate before the CE inspiral in their simulations, versus values derived using the analytical approximation of Paczyński & Sienkiewicz (Reference Paczyński and Sienkiewicz1972), using variables measured from the simulations. In that work the comparison was satisfactory, but, we argue here, it was also fortuitous and due to the specific value of the stellar radius they measured from simulations, which should have been considered very uncertain.
In Figure A2 we show a recreation of their Figure 7, where their numerical mass transfer rate was compared with the analytical prescription using a stellar radius of 107 R
$_{\odot}$
. This value was measured from the core of the primary in a direction perpendicular to the orbital plane (the direction with the least amount of distortion due to the binary interaction). The radius of the star in that direction was defined as the distance between the core and the least dense particles, then adding twice their SPH smoothing lengths as has been standard practice for SPH simulations (e.g. Nandez, Ivanova, & Lombardi Reference Nandez and Ivanova2013). The radius value is used in the analytical expression to determine by how much the stellar radius exceeds the
$L_1$
point. This ‘excess’ is then taken to the third power to determine, along with other variables, the mass transfer rate. As it happens the radius value determined leads to an analytical mass transfer rate that closely matches the numerical mass transfer rate! However, a radius value as little as 20% smaller (88 R
$_{\odot}$
) changes the fit considerably (Figure A2). This goes to show that when a satisfactory answer is found, one tends not to investigate matters further.
Given the difficulty of defining radii in SPH even for spherical stars, and given the complexity of measuring a distorted star, we conclude that it is uninstructive to carry out this comparison.
A recreation of figure 7 of Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019). The orange and blue lines are those used in that work, representing the mass transfer rate calculated using the analytical equation for mass transfer, along with quantities measured using the simulation (orange line) and the mass transfer measured directly from the simulation (blue line). The green line shows the mass transfer rate calculated analytically, but this time using a somewhat smaller value for the stellar radius, which is however still a reasonable estimate, resulting in a significantly lower rate of mass transfer and a mismatch between analytical theory and the simulation.
Appendix C. Energy and angular momentum conservation
Energy and angular momentum are conserved to excellent precision in all of our simulations. Figure A3 shows the 68IH simulation as an example. In all of our simulations momentum is conserved to within 0.2%, and total energy is conserved to within 0.1%.
Plots of angular momentum (left) and energy (right) evolution for the 68H simulation, as an example. In each plot the subscripts tot, b, u, c, env, refer to the total, bound, unbound, cores, and gas, respectively. In the plot of energy, K, U, and
$\phi$
, are the kinetic energy, thermal energy, and the potential energy, respectively, while E is used for the combination of potential, kinetic, and thermal as appropriate to the cores, envelope or total.
Open access
