11institutetext: Heidelberger Institut für Theoretische Studien, Schloss-
Wolfsbrunnenweg 35, 69118 Heidelberg, Germany 22institutetext: Zentrum für
Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut,
Mönchhofstr. 12-14, 69120 Heidelberg, Germany 33institutetext: University of
Oxford, St Edmund Hall, Oxford OX1 4AR, UK 44institutetext: Zentrum für
Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik,
Philosophenweg 12, 69120 Heidelberg, Germany

# Going from 3D to 1D: A one-dimensional approach to common-envelope evolution

V.A. Bronner<EMAIL_ADDRESS>from 3D to 1D: A one-dimensional
approach to common-envelope evolutionGoing from 3D to 1D: A one-dimensional
approach to common-envelope evolution F.R.N. Schneider Going from 3D to 1D: A
one-dimensional approach to common-envelope evolutionGoing from 3D to 1D: A
one-dimensional approach to common-envelope evolutionGoing from 3D to 1D: A
one-dimensional approach to common-envelope evolutionGoing from 3D to 1D: A
one-dimensional approach to common-envelope evolution Ph. Podsiadlowski Going
from 3D to 1D: A one-dimensional approach to common-envelope evolutionGoing
from 3D to 1D: A one-dimensional approach to common-envelope evolutionGoing
from 3D to 1D: A one-dimensional approach to common-envelope evolutionGoing
from 3D to 1D: A one-dimensional approach to common-envelope evolution F.K.
Röpke Going from 3D to 1D: A one-dimensional approach to common-envelope
evolutionGoing from 3D to 1D: A one-dimensional approach to common-envelope
evolutionGoing from 3D to 1D: A one-dimensional approach to common-envelope
evolutionGoing from 3D to 1D: A one-dimensional approach to common-envelope
evolution

(Received xxx / Accepted yyy)

The common-envelope (CE) phase is a crucial stage in binary star evolution
because the orbital separation can shrink drastically while ejecting the
envelope of a giant star. Three-dimensional (3D) hydrodynamic simulations of
CE evolution are indispensable to learning about the mechanisms that play a
role during the CE phase. While these simulations offer great insight, they
are computationally expensive. We propose a one-dimensional (1D) model to
simulate the CE phase within the stellar evolution code MESA by using a
parametric drag force prescription for dynamical drag and adding the released
orbital energy as heat into the envelope. We compute CE events of a 0.97M_⊙
asymptotic giant-branch star and a point mass companion with mass ratios of
0.25, 0.50, and 0.75, and compare them to 3D simulations of the same setup.
The 1D CE model contains two free parameters, which we demonstrate are both
needed to fit the spiral-in behavior and the fraction of ejected envelope mass
of the 1D method to the 3D simulations. For mass ratios of 0.25 and 0.50, we
find good-fitting 1D simulations, while for a mass ratio of 0.75, we do not
find a satisfactory fit to the 3D simulation as some of the assumptions in the
1D method are no longer valid. In all our simulations, we find that the
released recombination energy is important to accelerate the envelope and
drive the ejection.

###### Key Words.:

hydrodynamics – methods: numerical – stars: AGB and post-AGB – binaries: close
– stars: mass-loss

##  1 Introduction

Common-envelope (CE) evolution, originally proposed by Paczyński (1976), is a
phase in binary star evolution during which a giant star engulfs its more
compact companion, leading to the formation of a shared envelope around the
core of the giant star and the companion. Frictional forces cause the binary
orbit to shrink. The released orbital energy expands the envelope, which can
result in complete or partial envelope ejection while forming a close binary
or a merger of the two stars if the envelope cannot be fully ejected (for
reviews, see Ivanova et al. 2013, De Marco & Izzard 2017 and Röpke & De Marco
2023). Today, it is believed that CE evolution is the main mechanism for
converting wide binary star systems into close binaries (Ivanova et al. 2013).
Therefore, CE evolution plays an important, often essential role in
understanding the formation channels of close binary systems: such as the CE
channel for gravitational-wave sources (Tutukov & Yungelson 1993; Belczynski
et al. 2002; Voss & Tauris 2003; Eldridge & Stanway 2016; Stevenson et al.
2017; Kruckow et al. 2018;Vigna-Gómez et al. 2018;Spera et al. 2019; for a
more complete summary of references see, e.g., Mandel & Broekgaarden 2022),
possibly the progenitors of SNe Ia (Iben & Tutukov 1984; Webbink 1984; Whelan
& Iben 1973), X-ray binaries (Tauris & van den Heuvel 2006), cataclysmic
variables (Warner 1995), short-period hot subdwarfs (Heber 2009), and even
potentially gamma-ray burst sources (Fryer & Woosley 1998; Izzard et al. 2004;
Detmers et al. 2008).

Soon after CE evolution was proposed, one-dimensional (1D) simulations were
carried out by Taam et al. (1978) and Meyer & Meyer-Hofmeister (1979). Using
1D hydrodynamic simulations, Podsiadlowski (2001) found that the CE evolution
can be divided into distinct phases. During the _initial phase_ , the co-
rotation between the binary star and the envelope of the giant star is lost.
Frictional forces then lead to the rapid spiral-in of the binary star during
the _plunge-in phase_ , which is the shortest phase during CE evolution,
happening on the dynamical timescale of the binary. As the CE expands, the
spiral-in slows down, which might then lead to the _self-regulated spiral-in
phase_ , its _termination_ and the _post-CE evolution_ , during which the
changes of the orbital separation are small compared to the plunge in phase
(see Ivanova et al. (2013) for the definitions of the final stages). For the
following discussion in this work, we follow the approach of Röpke & De Marco
(2023) and combine all stages after the plunge-in phase to the _post-plunge-in
phase_.

Three dimensional (3D) (magneto-)hydrodynamic simulations of CE events (Ricker
& Taam 2008, 2012; Passy et al. 2012; Nandez et al. 2014; Ohlmann et al.
2016a, b; Iaconi et al. 2017; Chamandy et al. 2018; Prust & Chang 2019;
Reichardt et al. 2019; Sand et al. 2020; Ondratschek et al. 2022; Lau et al.
2022a, b; Moreno et al. 2022) are currently the best tool to study the
physical mechanisms that cause the spiral-in of the binary and the ejection of
the envelope. These simulations are computationally expensive and easily need
more than e5 core-hours. They allow us to study CE events in a case-by-case
study, but they are infeasible when studying larger populations of systems
undergoing CE evolution.

To predict the outcome of a CE event, an energy formalism was introduced
(Webbink 1984; Livio & Soker 1988). The released orbital energy is compared to
the binding energy of the envelope, where the CE efficiency
$\alpha_{\mathrm{CE}}$ determines the fraction of the orbital energy that is
used to unbind the envelope. Classically, an efficiency
$\alpha_{\mathrm{CE}}\leq 1$ is assumed to account for energy conservation.
However, energy sources other than the orbital energy might be available,
e.g., the ionization energy from the recombination of hydrogen and helium (Han
et al. 1995; Ivanova et al. 2015), accretion on the companion (Chamandy et al.
2018), the formation of jets (Shiber et al. 2019), and dust formation (Glanz &
Perets 2018; Iaconi et al. 2020; Reichardt et al. 2020). These extra energy
sources allow physical scenarios with $\alpha_{\mathrm{CE}}>1$. Therefore, the
energy formalism cannot be used directly by itself to predict the outcome of a
CE event, since the efficiency parameter $\alpha_{\mathrm{CE}}$ varies from
system to system (Iaconi & De Marco 2019). Recent efforts have been made by
Marchant et al. (2021) and Hirai & Mandel (2022) to improve on the energy
formalism.

Even though 3D simulations are constantly improving and more computational
power is available, 1D simulations of the CE phase remain a useful tool (cf.
Ivanova & Nandez 2016; Clayton et al. 2017; Fragos et al. 2019; Trani et al.
2022). The comparatively low computational costs of 1D simulations compared to
3D simulations open the opportunity to study a larger number of CE systems and
explore the possible parameter space (Ivanova et al. 2013). This makes 1D
simulations more versatile than 3D simulations. In addition to 1D simulations,
semi-analytic models calibrated on 3D simulation (e.g., Trani et al. 2022)
have similar advantages.

As outlined above, it is still a difficult task to predict the outcome of CE
events, because the input physics of 1D methods are incomplete, prohibiting
accurate simulations of CE evolution, while 3D simulations are computationally
too expensive. We introduce a new 1D method, to efficiently compute CE events.
The assumptions made in 1D models, such as spherical symmetry and energy
deposition in shells, do not capture the physical processes involved in CE
events. To compensate for this, we introduce free parameters, which we
calibrate on the characteristics (e.g., orbital separation, the mass of
ejected envelope) of 3D simulations. We test to see how many parameters are
needed, such that the 1D simulation can reproduce the characteristics of 3D
simulations. Our method needs to be calibrated on 3D simulations, to obtain
the most accurate results while still retaining the computational advantages
of a 1D simulation. If possible, the calibrated method can be used to simulate
CE events and predict their outcome.

This paper is structured as follows. In Sect. 2, we describe our new 1D CE
method. Then, we apply this method to simulate CE events of a 0.97M_⊙
asymptotic giant (AGB) star and point-mass companions. Detailed results of the
simulation with a mass ratio of 0.25 are shown in Sect. 3, before presenting
the results for simulations with mass ratios of 0.50 and 0.75 in Sect. 4.
Finally, we discuss our proposed model and the results in Sect. 5, and
conclude in Sect. 6.

##  2 Methods

We model the 3D CE simulation of Sand et al. (2020) in the 1D stellar-
evolution MESA revision 12778 (Paxton et al. 2011, 2013, 2015, 2018, 2019). We
describe how to obtain the same pre-CE giant star as Sand et al. (2020) in
Sect. 2.1 and our general hydro setup in Sect. 2.2. For numerical reasons, we
first perform a relaxation run of the pre-CE model, which we explain in Sect.
2.3. Frictional forces are hindering the motion of the companion inside the
envelope and cause the orbit to decay. In Sect. 2.4, we introduce a parametric
prescription of this drag force and show how the drag force modifies the
equations of motion that we integrate to obtain the binary orbit. The mass of
the companion modifies the gravitational potential of the giant star, which we
describe in Sect. 2.5. In Sect. 2.6, we show how we model the back reaction of
the companion on the envelope by artificially heating the envelope layers
around the companion. The unbound layers at the surface of the CE are removed
from the simulation, as explained in Sect. 2.7. Within this model, we use two
calibration parameters; one parameter determines the strength of the drag
force and the other sets the size of the artificially heated zone. In Sect.
2.8 we show how we calibrate these parameters by comparing our 1D simulations
to 3D hydrodynamic simulations of CE events of Sand et al. (2020).

###  2.1 Initial model

We evolve an initially 1.2M_⊙ star with metallicity $Z=0.02$ from the zero-age
main sequence to the ABG phase and stop the evolution once the model reaches a
mass of 0.97M_⊙. For this run, we use the default MESA settings, except for
changing the wind-loss parameters in the Reimers prescription to
$\eta_{\mathrm{R}}=0.5$ (Reimers 1975) and in the Blöcker prescription to
$\eta_{\mathrm{B}}=0.1$ (Blöcker 1995). Additionally, we disable the MESA
‘gold tolerances’ to evolve the model through the helium flash, and we use a
mixing-length parameter of $\alpha_{\mathrm{MLT}}=2$.

Figure 1: The initial density profile $\rho_{\mathrm{ini,1D}}$ of the AGB star
in comparison to the initial density profile $\rho_{\mathrm{ini,3D}}$ of the
AGB star in Sand et al. (2020), as well as the 1D and 3D density profiles
after the relaxation run $\rho_{\mathrm{relax,1D}}$ and
$\rho_{\mathrm{relax,3D}}$. Relative deviations between the density profiles
are shown along the secondary ordinate.

Our initial model of the ABG star has a radius of 166R_⊙. This is similar to
the radius of 172R_⊙ for the initial model of Sand et al. (2020). The density
in the envelope of our initial model is ${\sim}\,5-10\,\%$ larger compared to
the initial model for the 3D simulations (Fig. 1), with deviations of up to
$45\,\%$ close to the surface. We attribute this to the older MESA revision
7624 used to construct the initial pre-CE model in Sand et al. (2020), which
makes it difficult to reproduce the exact same model.

###  2.2 Hydrodynamic setup in MESA

We simulate the CE evolution in MESA using its single star module star
together with an artificial heating source to account for the presence of a
companion star. We switch to the implicit hydrodynamics in MESA by using an
HLLC solver (Paxton et al. 2018). With implicit hydrodynamics, we can study
the dynamical evolution of the CE and also capture potential shocks. To better
resolve shocks near the surface of the envelope, we impose a vanishing
compression gradient for the surface boundary condition that sets the surface
density (Paxton et al. 2015; Grott et al. 2005). For the other surface
boundary condition, here the temperature, we use a black body.

We use an adaptive mesh refinement in the velocity, optimized for hydrodynamic
simulations in MESA (Paxton et al. 2018), with uniform spacing in $\log\,r$
for the cells. The number of cells is chosen to be similar to the number of
cells in the initial hydrostatic model. As this mesh refinement does not
support rotation of the star, we simulate the CE without accounting for
rotation in the stellar structure equations.

The hydrodynamic equations with MESA’s HLLC scheme are solved implicitly in
time. Thus, the time-steps of our simulations are not limited by the Courant-
Friedrichs-Lewy (CFL) criterion (Courant et al. 1928; Paxton et al. 2018).
However, to resolve shock waves that travel through the envelope, we limit the
time-steps to 45% of the maximum time-step given by the CFL criterion.
Furthermore, we find that limiting the time-steps to $<\qty{e-4}{yr}$ helps
with the numerical stability of the simulations. All the 1D CE simulations are
run for as long as possible and stop for numerical reasons, e.g., too short
time-steps.

###  2.3 Relaxation run

Before starting the CE simulation, we perform a relaxation, where we switch
from the standard hydrostatic solver to the implicit hydrodynamic solver and
change the surface boundary conditions. The relaxation run lasts for 2yr or
about 15 dynamical timescales111We start the relaxation run with an initial
time-step of e-7yr and set the maximum time-step during the relaxation run to
e-3yr.. During the final 0.2yr of the relaxation run, the maximum time-step is
linearly decreased to a value of e-7yr. We find that this helps to start the
CE simulation from the relaxed model without numerical artifacts.

During the relaxation run, the envelope expands to 189R_⊙. The density profile
of the envelope of the AGB star changes by less than 5% during the relaxation
run (Fig. 1). At the end of the relaxation run, the density close to the
surface of the star decreased by 60%. This shows that mostly the surface
layers of the star expand during the relaxation run, leaving the envelope
structure unchanged. Fig. 1 also compares the density profile after the
relaxation run to the relaxed density profile of Sand et al. (2020).
Throughout most parts of the envelope, the density profiles deviate by less
than 10%, except for larger deviations towards to surface. The deviations
towards smaller radii originate from the cut-out core in the 3D simulation.

###  2.4 Orbital evolution

We integrate the equations of motion of the classical two-body problem
consisting of an extended sphere at location $\boldsymbol{x}_{1}$,
representing the giant star, and a point-mass particle at location
$\boldsymbol{x}_{2}$, representing the companion, to compute the orbital
evolution of the binary star. In general, quantities with index 1 refer to the
giant star, while index 2 refers to the companion. As the companion is
revolving within the giant star during the CE event, only the mass of the
giant star enclosed by the orbit $M_{1,a_{\mathrm{orb}}}$ needs to be
considered in the equations of motion,

$\displaystyle\boldsymbol{\ddot{x}}_{1}$
$\displaystyle=\dfrac{GM_{2}}{a_{\mathrm{orb}}^{2}}\boldsymbol{\hat{a}}_{\mathrm{orb}}$
(1) $\displaystyle\boldsymbol{\ddot{x}}_{2}$
$\displaystyle=-\dfrac{GM_{1,a_{\mathrm{orb}}}}{a_{\mathrm{orb}}^{2}}\boldsymbol{\hat{a}}_{\mathrm{orb}}-\dfrac{F_{\mathrm{d}}}{M_{2}}\boldsymbol{\hat{v}}_{\mathrm{rel}}$
(2)

where the orbital separation, i.e., the distance from the center of the giant
star to the point-mass companion, is given by
$\boldsymbol{a}_{\mathrm{orb}}=\boldsymbol{x}_{2}-\boldsymbol{x}_{1}$. The
mass of the companion is given by $M_{2}$ and $G$ is the gravitational
constant. Vector quantities in non-bold notation denote the absolute value of
that quantity and $\boldsymbol{\hat{x}}$ is a unit vector in the direction of
$\boldsymbol{x}$. Because the companion moves within the envelope of the giant
star, it is subject to a drag force
$\boldsymbol{F}_{\mathrm{d}}=-F_{\mathrm{d}}\boldsymbol{\hat{v}}_{\mathrm{rel}}$
caused by dynamical friction, opposing the relative velocity
$\boldsymbol{v}_{\mathrm{rel}}$ of the companion to the envelope.

We use the semi-analytic formalism of Kim (2010) and Kim & Kim (2007) to
calculate the drag force. Kim & Kim (2007) extended the classical model of
dynamical friction of Ostriker (1999) to perturbers moving in circular orbits
with radius $r$ through a homogeneous background medium of density $\rho$. Kim
(2010) then extended the model further to account for possible non-linear
effects arising from the circular motion of the perturber. They define the two
dimensionless parameters

$\mathcal{B}=\frac{GM_{2}}{c_{\mathrm{s}}^{2}r},\qquad\eta_{\mathcal{B}}=\frac{\mathcal{B}}{\mathcal{M}^{2}-1},$
(3)

with $M_{2}$ the mass of the perturber, i.e., the companion in our case,
$c_{\mathrm{s}}$ the sound speed, and
$\mathcal{M}=v_{\mathrm{rel}}/c_{\mathrm{s}}$ the Mach number. They find a
correction term $f_{\rho}$ for the density to account for non-linear effects,

$\rho_{\mathrm{nl}}=f_{\rho}\rho=\left(1+\frac{0.46\,\mathcal{B}^{1.1}}{{(\mathcal{M}^{2}-1)}^{0.11}}\right)\rho.$
(4)

Together, the drag force is then given by

$F_{\mathrm{d}}=\begin{cases}C_{\mathrm{d}}\dfrac{4\pi\rho_{\mathrm{nl}}{(GM_{2})}^{2}}{v_{\mathrm{rel}}^{2}}\left(\dfrac{0.7}{\eta_{\mathcal{B}}^{0.5}}\right)\quad&\mathrm{for}\,\eta_{\mathcal{B}}>0.1\wedge\mathcal{M}>1.01,\\\\[10.00002pt]
C_{\mathrm{d}}\dfrac{4\pi\rho{(GM_{2})}^{2}}{v_{\mathrm{rel}}^{2}}I&\mathrm{else},\end{cases}$
(5)

where the condition $\eta_{\mathrm{B}}>0.1$ and $\mathcal{M}>1.01$ corresponds
to a regime in which non-linear effects are present that modify the drag
force, as opposed to the classical linear regime for $\eta_{\mathcal{B}}<0.1$
or $\mathcal{M}<1.01$. The drag-force parameter $C_{\mathrm{d}}$ is introduced
as a free parameter in our model and modifies the strength of the drag force.
We modified the original condition by Kim (2010) for the non-linear regime
from $\mathcal{M}>1$ to $\mathcal{M}>1.01$ in Eq. (5) to avoid divergence
issues in the definitions of $\eta_{\mathcal{B}}$ and $\rho_{\mathrm{nl}}$.
According to Kim & Kim (2007), the Coulomb logarithm $I$ in the linear regime
is given by

$I=\begin{cases}0.7706\ln\left(\dfrac{1+\mathcal{M}}{1.0004-0.9185\,\mathcal{M}}\right)-1.473\,\mathcal{M}\quad\mathrm{for}\,\mathcal{M}<1.0,\\\\[10.00002pt]
\ln\left(330\dfrac{r}{r_{\mathrm{min}}}\dfrac{(\mathcal{M}-0.71)^{5.72}}{\mathcal{M}^{9.58}}\right)\quad\mathrm{for}\,1.0\leq\mathcal{M}<4.4,\\\\[10.00002pt]
\ln\left(\dfrac{r/r_{\mathrm{min}}}{0.11\,\mathcal{M}+1.65}\right)\quad\mathrm{for}\,4.4\leq\mathcal{M}.\end{cases}$
(6)

We choose $r_{\mathrm{min}}=3.1\,{\rm R}_{\odot}$, a measure of the size of
the companion, which is the same value as the softening-length of the
gravitational potential of the companion in the simulations of Sand et al.
(2020).

The initial orbital separation $a_{\mathrm{orb,ini}}$ in Sand et al. (2020) is
chosen such that the giant star overfills its Roche lobe, initiating the CE
event. In particular, the companion is initially outside of the giant star,
i.e., $a_{\mathrm{orb,ini}}>R_{1}$. This setup of the initial orbit would not
lead to a CE event in our 1D model on a dynamical timescale, because the
companion does not exert a tidal force on the giant star. Additionally, we
assume a vanishing density around the giant star which leads to a vanishing
drag force and no spiral-in behavior. Therefore, to initiate the CE event, we
choose $a_{\mathrm{orb,ini}}=\qty{160}{{\rm R}_{\odot}}$, i.e., we place the
companion well inside the giant’s envelope of 166R_⊙ before the relaxation.
The initial orbit is set up to be circular, following the setup in Sand et al.
(2020).

The giant star in the simulations of Sand et al. (2020) is initially set up to
rotate rigidly at 95% of the initial orbital frequency with the angular
velocity vector $\boldsymbol{\Omega}$ pointing perpendicular to the orbital
plane. We use the same $\boldsymbol{\Omega}$ for our 1D simulations as in Sand
et al. (2020), which means in particular that $\boldsymbol{\Omega}$ varies
with mass ratio $q=M_{2}/M_{1}$.

The relative velocity of the companion to the envelope of the giant star is
given by

$\boldsymbol{v}_{\mathrm{rel}}=\boldsymbol{\dot{x}}_{2}-\boldsymbol{\dot{x}}_{1}-\left(\boldsymbol{\Omega}\times\boldsymbol{a}_{\mathrm{orb}}+v_{r}(a_{\mathrm{orb}})\boldsymbol{\hat{a}}_{\mathrm{orb}}\right),$
(7)

where $v_{r}(a_{\mathrm{orb}})$ is the (radial) expansion velocity of the CE
at the location of the companion. While the giant star in the MESA simulation
does not rotate (Sect. 2.2), a rigidly rotating envelope with a constant
angular velocity $\boldsymbol{\Omega}$ in time is assumed for calculating the
relative velocity. The magnitude of the relative velocity enters the
calculation of the drag force via Eq. (5).

The equations of motion (Eqs. 1 and 2) are integrated in parallel with the
MESA simulation using a fifth-order explicit Runge-Kutta method with an
embedded fourth-order error estimation (Cash & Karp 1990). Relative and
absolute error tolerances are set to $10^{-10}$, and the maximum integration
time-step is set so that at least ten orbit-integration steps are taken during
one MESA time-step.

###  2.5 Modified gravitational potential

For the simulations of the CE evolution with our 1D model, we use the single-
star module star of MESA to model the giant star. However, because of the mass
of the companion, which is revolving inside the envelope of the giant star,
the gravitational potential can no longer be approximated by that of a single
star. To account for the change in the gravitational potential, we modify the
gravitational constant $G$ to a radius-dependent gravitational constant
$\tilde{G}(r)$ with

$\tilde{G}(r)=\begin{cases}G\qquad&\mathrm{for\,}r<a_{\mathrm{orb}},\\\
G\left(1+\dfrac{M_{2}}{M_{1,r}}\right)\qquad&\mathrm{for\,}r\geq
a_{\mathrm{orb}},\end{cases}$ (8)

where $M_{1,r}$ is the mass of the giant star within radius $r$ (Podsiadlowski
et al. 1992). This ensures that the layers outside the orbit are bound more
strongly to the system because of the additional mass of the companion. Layers
within the orbit are not affected. Changing the gravitational constant
modifies not only the binding energy of the outer layers but also the entire
envelope profile, e.g., density and pressure.

###  2.6 CE heating

The drag force acting on the companion inside the CE dissipates the orbital
energy of the binary system, causing the orbital separation to shrink. The
rate at which orbital energy is dissipated is given by
$\dot{E}_{\mathrm{dis}}=\boldsymbol{F}_{\mathrm{d}}\cdot\boldsymbol{v}_{\mathrm{rel}}$.
We add the dissipated energy from the binary orbit by artificially increasing
the internal energy of the envelope, i.e., heating with the same rate as the
orbital energy is lost, i.e.,
$\dot{E}_{\mathrm{heat}}=-\dot{E}_{\mathrm{dis}}$.

The range over which energy is injected is approximated by the accretion
radius of the Bondi-Lyttleton-Hoyle accretion model (Hoyle & Lyttleton 1941;
Bondi & Hoyle 1944). The accretion radius profile $R_{\mathrm{a}}(r)$ inside
the envelope is given by

$R_{\mathrm{a}}(r)=\frac{2\,GM_{2}}{{\Delta
v(r)}^{2}+{c_{\mathrm{s}}(r)}^{2}},$ (9)

where $\Delta v(r)$ is the relative velocity of the spherical shell at
coordinate $r$ with respect to the companion, and $c_{\mathrm{s}}(r)$ is the
sound speed. The relative velocity is given by

$\Delta v(r)^{2}=\Delta v_{r}(r)^{2}+\Delta v_{\varphi}(r)^{2},$ (10)

where

$\displaystyle\Delta v_{r}(r)$
$\displaystyle=v_{\mathrm{orb},r}-v_{r}(r)\quad\text{and}$ (11)
$\displaystyle\Delta v_{\varphi}(r)$
$\displaystyle=v_{\mathrm{orb},\varphi}-\Omega r.$ (12)

In the above equations,
$\boldsymbol{v}_{\mathrm{orb}}=\boldsymbol{\dot{x}}_{2}-\boldsymbol{\dot{x}}_{\mathrm{1}}$
is the orbital velocity and $v_{r}(r)$ is the (radial) expansion velocity of
the envelope. The quantities $v_{\mathrm{orb},r}$ and
$v_{\mathrm{orb},\varphi}$ are the projections of the orbital velocity vector
along the $r$-direction and the $\varphi$-direction.

The accretion radius defines the area of influence of the companion in the
giant’s envelope. We choose to heat all the envelope layers with a radial
coordinate between $r^{\mathrm{heat}}_{\min}$ and $r^{\mathrm{heat}}_{\max}$,
which are determined by solving

$\displaystyle r^{\mathrm{heat}}_{\min}$
$\displaystyle=a_{\mathrm{orb}}-C_{\mathrm{h}}R_{\mathrm{a}}(r^{\mathrm{heat}}_{\min}),$
(13) $\displaystyle r^{\mathrm{heat}}_{\max}$
$\displaystyle=a_{\mathrm{orb}}+C_{\mathrm{h}}R_{\mathrm{a}}(r^{\mathrm{heat}}_{\max}).$
(14)

The heating parameter $C_{\mathrm{h}}$ is introduced as a free parameter in
our models and determines the extent of the heating zone by modifying the
upper and lower boundary for the heating zone via Eqs. (13) and (14). If there
are multiple solutions222The accretion radius $R_{\mathrm{a}}(r)$ is not a
monotonically increasing function, allowing multiple solutions for Eqs. (13)
and (14). for Eqs. (13) and (14), we use the largest one for
$r^{\mathrm{heat}}_{\min}$ and the smallest one for
$r^{\mathrm{heat}}_{\max}$, i.e., the solutions closest to $a_{\mathrm{orb}}$.
This definition ensures that for all layers with $r^{\mathrm{heat}}_{\min}\leq
r\leq r^{\mathrm{heat}}_{\max}$ the radial separation from the companion is
smaller than the local accretion radius, i.e., these layers are
gravitationally deflected/focused by the companion. The heating rate
$\dot{\epsilon}_{\mathrm{heat}}$ then follows

$\dot{\epsilon}_{\mathrm{heat}}{\sim}\,\begin{cases}\exp\left[-\left(\dfrac{a_{\mathrm{orb}}-r}{a_{\mathrm{orb}}-r^{\mathrm{heat}}_{\min}}\right)^{2}\right]\quad&\mathrm{for}\,r^{\mathrm{heat}}_{\min}\leq
r<a_{\mathrm{orb}},\\\\[10.00002pt]
\exp\left[-\left(\dfrac{a_{\mathrm{orb}}-r}{a_{\mathrm{orb}}-r^{\mathrm{heat}}_{\max}}\right)^{2}\right]&\mathrm{for}\,a_{\mathrm{orb}}\leq
r\leq r^{\mathrm{heat}}_{\max},\\\\[10.00002pt] 0&\mathrm{else},\end{cases}$
(15)

which is normalized so that the total heating rate matches the dissipation
rate of the orbital energy, i.e.,
$\int\dot{\epsilon}_{\mathrm{heat}}\,\mathrm{d}m=\dot{E}_{\mathrm{heat}}$. In
the case where $r^{\mathrm{heat}}_{\max}$ is undefined by Eq. (14) because
$r\neq a_{\mathrm{orb}}+R_{\mathrm{a}}(r)$ for $a_{\mathrm{orb}}\leq r\leq
R_{1}$, we choose $r^{\mathrm{heat}}_{\max}=R_{1}$. The lower heating limit
$r^{\mathrm{heat}}_{\min}$ is always well defined because the sound speed
increases by orders of magnitude toward the giant’s core, causing the
accretion radius to decrease.

###  2.7 Dynamical envelope ejection

The heating perturbs the hydrostatic equilibrium of the envelope, causing it
to expand on the dynamical timescale. Envelope layers exceeding the local
escape velocity $v_{\mathrm{esc}}$ are formally unbound. The local escape
velocity is given by

$v_{\mathrm{esc}}(r)=\sqrt{\dfrac{2\,\tilde{G}(r)M_{1,r}}{r}},$ (16)

where $\tilde{G}(r)$ is the modified gravitational constant defined by Eq.
(8). We remove all surface layers with $v_{r}>v_{\mathrm{esc}}$ in our CE
simulations to avoid numerical difficulties in the unbound layers. If a
continuous layer of unbound mass333We only use the radial velocity $v_{r}$ in
the criterion for envelope ejection, i.e., we neglect the kinetic energy in
the rotation of the envelope. We find that excluding the kinetic energy in
rotation changes the fraction of the ejected envelope mass in the 3D
simulations of Sand et al. (2020) by less than 2% after the end of the
dynamical plunge-in phase. $M_{\mathrm{unbound}}$, i.e.,
$v_{r}>v_{\mathrm{esc}}$ for all shells in this layer, reaches the outer
boundary of the simulation domain, we remove the layer exponentially by
adopting a mass-loss rate

$\dot{M}=-\frac{M_{\mathrm{unbound}}}{\tau_{\dot{M}}},$ (17)

where $\tau_{\dot{M}}$ is the mass-loss timescale. In all our simulations, we
choose $\tau_{\dot{M}}=\qty{0.01}{yr}$, which is shorter than both the
dynamical timescale and the thermal timescale of the envelope. The applied
mass-loss prescription is the same as that in Clayton et al. (2017).

###  2.8 Comparison to 3D CE simulations

We compare and fit our 1D results to those of Sand et al. (2020). In
particular, we compare the time evolution of the orbital separation
$a_{\mathrm{orb}}(t)$, from hereon called the spiral-in curve, and the mass
fraction of the ejected envelope $f_{\mathrm{ej}}$,

$f_{\mathrm{ej}}=\frac{M_{\mathrm{env,ini}}-M_{\mathrm{env}}}{M_{\mathrm{env,ini}}},$
(18)

where $M_{\mathrm{env}}=M_{1}-M_{\mathrm{1,core}}$ is the envelope mass,
$M_{1}$ is the mass of the giant star and
$M_{\mathrm{1,core}}=\qty{0.545}{{\rm M}_{\odot}}$444The value of
$M_{\mathrm{1,core}}$ is taken from Sand et al. (2020) to ensure a meaningful
comparison.. In contrast to this, Sand et al. (2020) defined the total unbound
mass, and hence the mass fraction of the ejected envelope, by summing all
cells with positive total energy, i.e., the sum of the potential and kinetic
energy. To better match this quantity, we define a second envelope-ejection
fraction,

$\displaystyle
f_{\mathrm{ej,all}}=\frac{M_{\mathrm{env,ini}}-M_{\mathrm{env,bound}}}{M_{\mathrm{env,ini}}}\quad\mathrm{with}$
(19) $\displaystyle
M_{\mathrm{env,bound}}=M_{\mathrm{env}}-M_{\mathrm{unbound}},$ (20)

where $M_{\mathrm{unbound}}$ is the total mass of all layers in the envelope
with $v_{r}>v_{\mathrm{esc}}$. The difference between $f_{\mathrm{ej}}$ and
$f_{\mathrm{ej,all}}$ is that $f_{\mathrm{ej,all}}$ also includes the layers
with $v_{r}>v_{\mathrm{esc}}$ within the envelope.

We incorporate two free parameters in our model. The drag-force parameters
$C_{\mathrm{d}}$ modifies the strength of the drag force via Eq. (5) and the
heating parameter $C_{\mathrm{h}}$ modifies the extent of the heated layers
via Eqs (13) and (14). Both parameters are used to fit the spiral-in curves
and the mass fraction of the ejected envelope material from the 1D CE
simulations to the results of the corresponding 3D simulations in Sand et al.
(2020).

The initial orbital separations of our 1D simulations deviate from the initial
separations used by Sand et al. (2020) (see Sect. 2.4). Therefore, the spiral-
in curve of the 3D simulation is shifted by $\Delta t$ with respect to the
spiral-in curve of the 1D simulation. The optimal time shift $\Delta
t_{\mathrm{min}}$ for fixed $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ is found by
minimizing the mean relative deviation (MRD) between the spiral-in curves,

$\mathrm{MRD}(\Delta
t)=\frac{1}{T}\sum_{i}\frac{\left|a_{\mathrm{orb}}^{\mathrm{1D}}(t_{i})-a_{\mathrm{orb}}^{\mathrm{3D}}(t_{i}-\Delta
t)\right|}{a_{\mathrm{orb}}^{\mathrm{1D}}(t_{i})}\delta t_{i},$ (21)

where $T$ is the total simulated time, $\delta t_{i}$ is the time-step and
$a_{\mathrm{orb}}^{\mathrm{1D/3D}}$ are the orbital separations in the 1D and
3D simulation. For the calculation of the MRD, only the plunge-in phase is
considered, i.e., the phase during which the orbital separation changes
dynamically555We do not simulate the initial phase where the binary star loses
co-rotation, but rather start the simulation at the beginning of the plunge-in
phase. This is achieved because the stars have already lost co-rotation
according to our setup.. The end of the plunge-in phase is determined by

$\left\langle-\frac{\dot{a}_{\mathrm{orb}}}{a_{\mathrm{orb}}}P_{\mathrm{orb}}\right\rangle_{\mathrm{orbit}}<0.01,$
(22)

where the time average is computed over one full orbit. This defines the end
of the plunge-in phase when the average change in orbital separation during
one orbit is less than 1%. At the beginning of the CE phase in the 3D
simulations, a transition of the initially circular orbit towards the
dynamical spiral-in is observed. These first 100d are omitted in the MRD
calculations in Eq. (21).

The best fit of the 1D simulation to the 3D simulation is determined manually.
We find the best-fitting simulation by computing 1D CE models with varying
values for $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$. For each simulation, we find
$\Delta t_{\mathrm{min}}$ and then compare by eye the spiral-in curve and the
mass fraction of ejected envelope to the 3D simulation.

We tested the above described method at different spacial and temporal
resolutions and present the results in Appendix A.

##  3 Results for mass ratio $q=0.25$

In this section, we show the results of our 1D CE simulations with mass ratio
$q=0.25$ and compare our results with the 3D hydrodynamic simulations of Sand
et al. (2020). First, we present the best-fitting 1D CE simulations in Sect.
3.1. We then show how the drag-force parameter $C_{\mathrm{d}}$ and the
heating parameter $C_{\mathrm{h}}$ affect the outcome of the simulations in
Sect. 3.2. In Sect. 3.3, the dynamical processes taking place in the envelope
are discussed, before we describe recombination processes in Sect. 3.4. The
energy budget during the CE phase and the evolution of the drag force are
shown in Sects. 3.5 and 3.6 respectively.

###  3.1 Best fitting 1D simulation for $q=0.25$

Figure 2: Best fitting 1D CE simulation for $q=0.25$ when comparing to the 3D
CE simulation of Sand et al. (2020). Full lines show the orbital separation,
and the dotted and dash-dotted lines show the mass fraction of ejected
envelope $f_{\mathrm{ej}}$ and $f_{\mathrm{ej,all}}$ respectively. The results
of the 3D simulations are shown in orange, and the results from this work are
shown in blue. The vertical gray-dashed line indicates the end of the
dynamical plunge-in phase as determined by Eq. (22). All the results of the 3D
simulation are shifted by $\Delta t_{\mathrm{min}}=\qty{-58}{d}$ to compensate
for the difference in initial separation. The gray shaded region shows
$t<\qty{0}{d}$ for which we cannot compare the 1D simulation to the 3D
simulation.

The spiral-in curve of the best-fitting 1D simulation with parameters
$C_{\mathrm{d}}=0.23$ and $C_{\mathrm{h}}=4.0$ for a mass ratio of $q=0.25$ is
shown in Fig. 2. The spiral-in curve of the 3D benchmark model of Sand et al.
(2020) for the same mass ratio is also shown. We find $\Delta
t_{\mathrm{min}}=-\qty{58}{d}$ to best account for the difference in the
initial separation (see Sect. 2.8). The plunge-in phase is well reproduced by
the 1D simulation. The orbital separation of the 1D simulation after the
plunge-in phase is smaller than in the 3D simulation. In both the 1D and the
3D simulations, the post-plunge-in orbit shows a non-zero eccentricity. The
eccentricity $e$ of the orbit can be approximated by

$e=\frac{A-P}{A+P}$ (23)

where $A$ is the apastron distance and $P$ is the periastron distance. This
approximation is valid as long as the orbital separation is not changing
significantly during one orbit, e.g., during the post-plunge-in phase of the
CE evolution. In the 1D simulation, we find an eccentricity of 0.005, which is
approximately equal to the eccentricities of 0.006 in the 3D simulation (for a
summary, see Table 1). In both cases, the eccentricity is evaluated at the end
of the plunge-in phase as determined by Eq. (22), and averaged over 5 orbits.

The fraction of unbound envelope material $f_{\mathrm{ej}}$ is also shown in
Fig. 2. For the 1D simulation, two descriptions ($f_{\mathrm{ej}}$ and
$f_{\mathrm{ej,all}}$) are used to determine the fraction of the unbound
envelope (Sect. 2.8). In the 3D simulation, the envelope ejection starts at
$t<\qty{0}{d}$ and reaches 0.1 at $t=\qty{0}{d}$ (Fig. 2). As the companion
starts to enter the envelope, it causes one large spiral arm where mass is
immediately ejected (Sand et al. 2020). In the 1D simulation, envelope
ejection sets in later during the simulations. When considering
$f_{\mathrm{ej,all}}$, the envelope ejection starts at $t=\qty{700}{d}$
compared to $t=\qty{850}{d}$ based on $f_{\mathrm{ej}}$. This is expected
since $f_{\mathrm{ej,all}}$ captures all the unbound mass that is considered
for $f_{\mathrm{ej}}$ and, additionally, all layers with
$v_{r}>v_{\mathrm{esc}}$ inside the envelope. At the end of the 1D simulation,
about $70\,\%$ of the envelope is unbound, as well as $f_{\mathrm{ej}}\approx
f_{\mathrm{ej,all}}$, i.e., there are only unbound layers at the outer
boundary of the CE but not within the CE. It is important to note that the
slopes of the fraction of the ejected envelope, i.e., the envelope-ejection
rate, at the end of the simulations are similar in both the 1D and 3D case,
possibly suggesting a similar ejection mechanism. The envelope ejection in the
1D simulations relies on the energy released from the recombination of
hydrogen and helium and will be analyzed in more detail in Sect. 3.3.

The 1D simulation ended because of too short time-steps. However, it seems
appropriate to assume that the envelope ejection will be sustained if the
simulation is run for longer, possibly leading to a full envelope ejection. At
the end of the 3D simulation at 4000d, 91% of the envelope is ejected (cf.
Table 4 in Sand et al. 2020).

Although there are small quantitative differences between the 1D CE simulation
in this work and the 3D CE simulations of Sand et al. (2020), it is remarkable
that it is possible to achieve such similar results with our 1D model.

###  3.2 Role of $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ for the orbital spiral-
in and envelope ejection

Figure 3: Effects of $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ on the spiral-in
curves $a_{\mathrm{orb}}$ and the envelope-ejection fraction
$f_{\mathrm{ej}}$. In panel (a), the heating parameter is held constant at
$C_{\mathrm{h}}=4.00$ while varying the drag-force parameter $C_{\mathrm{d}}$.
In panel (b), the drag-force parameter is held constant at
$C_{\mathrm{d}}=0.23$ while varying the heating parameter $C_{\mathrm{h}}$.
The thick black lines in both panels show the results from the 3D hydrodynamic
CE simulation of Sand et al. (2020), which are shifted by $\Delta
t_{\mathrm{min}}=-\qty{58}{d}$, i.e., the time shift determined using the
best-fitting 1D simulation with $C_{\mathrm{d}}=0.23$ and
$C_{\mathrm{h}}=4.00$.

In Fig. 3, we show how $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ affect the
orbital separation and the mass fraction of the ejected envelope by varying
$C_{\mathrm{d}}$ between $0.1$ and $1.0$ at the best-fit, constant
$C_{\mathrm{h}}=4.00$ (Fig. 3a), and varying $C_{\mathrm{h}}$ between $1.0$
and $10.0$ at the best-fit, constant $C_{\mathrm{d}}=0.23$ (Fig. 3b). For
simplicity, we only show $f_{\mathrm{ej}}$ in Fig. 3. However, we find for all
simulations that $f_{\mathrm{ej}}\approx f_{\mathrm{ej,all}}$ during the late
phases of the simulations, where the envelope-ejection rate is constant (Figs.
2 and 9).

The slope of the spiral-in is mostly determined by the strength of the drag
force via $C_{\mathrm{d}}$ (Fig. 3a). When varying $C_{\mathrm{d}}$ between
$0.1$ and $1.0$ at fixed $C_{\mathrm{h}}$, the spiral-in timescale decreases
from more than 1500d to less than 200d. Additionally, the envelope ejection
starts earlier for larger $C_{\mathrm{d}}$ while also ejecting a larger
fraction of the envelope. The start of the rapid envelope ejection occurs at
similar orbital separations of about $\qty{30}{{\rm R}_{\odot}}$ throughout
the different simulations. The post-plunge-in separation changes only slightly
when varying the drag-force parameter $C_{\mathrm{d}}$.

We find that the heating parameter $C_{\mathrm{h}}$ plays an important role in
determining the mass fraction of the ejected envelope, but not so much in
setting the post-plunge-in separation (Fig. 3b). Increasing the heating
parameter causes a slightly deeper spiral-in of the companion and a higher
fraction of ejected envelope. The behavior of $f_{\mathrm{ej}}$ for
$C_{\mathrm{d}}=0.23$ and varying $C_{\mathrm{h}}$ between 3.2 and 4.2 is
noteworthy. In this range, a changing $C_{\mathrm{h}}$ results in a large
change in $f_{\mathrm{ej}}$ while the orbital separation at the end of the
plunge-in changes only between 22.9R_⊙ for $C_{\mathrm{h}}=3.2$ and 21.3R_⊙
for $C_{\mathrm{h}}=4.2$. For the same models, the envelope-ejection curves
show a rapid ejection event at the beginning, after which a more steady
envelope ejection settles in. In this second phase, the envelope-ejection
rates seem to be the same for all the simulations, which are also in agreement
with the 3D simulation. This is another indicator, that the envelope ejection,
especially in the later phase of the simulation, is mainly caused by
recombination processes because the recombination rate is comparable between
the models (Sect. 3.4). Only the timing and strength of the initial rapid
ejection event is changing when varying $C_{\mathrm{h}}$. Therefore, the
heating parameter $C_{\mathrm{h}}$ is important to set the envelope ejection
in the 1D CE simulation.

From these experiments, we conclude that the drag-force parameter
$C_{\mathrm{d}}$ mostly affects the orbital separation during the plunge-in
phase by determining the spiral-in timescale. The final mass-fraction of the
ejected envelope is determined by both the heating parameter $C_{\mathrm{h}}$
and the drag-force parameter $C_{\mathrm{d}}$. The orbital separation after
the plunge-in phase changes only slightly when varying $C_{\mathrm{d}}$ and
$C_{\mathrm{h}}$. This suggests, that there is a more fundamental principle
that sets the post-plunge-in separation and that does not depend on the
details of the plunge-in-phase. When comparing the 1D model to the 3D
simulation, we cannot match the results with only one of the two parameters.
Therefore, 2 parameters are needed in 1D CE models similar to ours to
reproduce the spiral-in curve as well as the envelope ejection.

###  3.3 Dynamical evolution of the envelope

Figure 4: Kippenhahn diagrams of the simulation with $q=0.25$,
$C_{\mathrm{d}}=0.23$ and $C_{\mathrm{h}}=4.0$ showing the radial velocity in
terms of the local escape velocity (panel a), where layers colored in red have
a radial velocity larger than the local escape velocity, the ionization
fractions of H ii, He ii and He iii (panels b, c and d respectively), the
heating kernel as described in Eq. (15) (panel e), and the specific heating
rate (panel f). The red line indicates the orbital separation between the
companion and the core of the giant star. The gray-dotted region is heated
during the CE simulation. The cyan dashed line shows $R_{\tau=10}$ which
traces the hydrogen-recombination radius. The gray dashed lines represent
envelope-mass-fractions of 0.95, 0.9, 0.8, 0.6, 0.4, 0.2, 0.1, 0.06, 0.04,
0.03 from the surface toward the center, and visualize the expansion of the
envelope. The white patches in panel (a) show layers that have an inward
(negative) radial velocity. Figure 5: Kippenhahn diagrams similar to Fig. 4 of
the simulation with $q=0.25$, $C_{\mathrm{d}}=0.23$ and $C_{\mathrm{h}}=4.0$
showing the Mach number $\mathcal{M}=v_{r}/c_{\mathrm{s}}$, velocity
divergence $\mathrm{d}v/\mathrm{d}r$, specific entropy $s$, and opacity
$\kappa$. The white hatched zones in panel (c) indicate convective mixing. For
better visualizing the velocity divergence, the log modulus transformation,
$\mathrm{logmod}(x)=\mathrm{sgn}(x)\log_{10}(|x|+1)$, is used (John & Draper
1980). The arrows in panel (b) roughly indicate the recombination zones of
${H\textsc{i}}$, ${He\textsc{i}}$ and ${He\textsc{ii}}$ (see also panels (b),
(c) and (d) in Fig. 4).

The dynamical evolution of the envelope is shown by several Kippenhahn
diagrams in Figs. 4 and 5. During the first ${\sim}\,\qty{200}{d}$, the
envelope oscillates radially. These oscillations are artifacts of the change
to the hydrodynamic mode in MESA as well as the change in the boundary
condition, which are not yet fully damped during the relaxation run (see Sect.
2.3). As the velocities in these oscillations are small compared to the
expected expansion velocity which is comparable to the escape velocity, they
do not affect the further evolution of the CE.

After this initial phase, the envelope begins to expand, as more and more heat
is injected into the outer layers (see Fig. 9). The partial ionization zones
of hydrogen and helium are expanding together with the envelope, as can be
seen in Figs. 4b, c, and d, and no recombination is happening yet.

For $t<\qty{700}{d}$, we find $r_{\mathrm{max}}^{\mathrm{heat}}=R_{1}$, i.e.,
the upper boundary of the heating zone coincides with the outer boundary of
the simulation domain. Once $r_{\mathrm{max}}^{\mathrm{heat}}<R_{1}$ for
$t>\qty{700}{d}$, the specific heating rate $\dot{\epsilon}_{\mathrm{heat}}$
increases by about one order of magnitude because the mass of the material
that is heated is decreasing significantly (Figs. 4e and f). As the upper
heating radius moves deeper in the envelope, the partial ionization zones of
singly ionized hydrogen and helium are no longer heated and recombination sets
in. The beginning of recombination and the end of heating in these zones are
causally connected, as the heating provides an energy source that keeps the
atoms ionized. This is also the reason why $C_{\mathrm{h}}$ determines the
ejection of the envelope so strongly.

The layer in which hydrogen recombination takes place moves outwards in radius
at a typical optical depth of $\tau=10$ ($R_{\tau=10}$), hereafter referred to
as the hydrogen-recombination radius, below the outer boundary of the
simulation. We find that $R_{\tau=10}$ stays always close to the hydrogen
recombination front because, above the hydrogen recombination layer, the
opacity drops dramatically (Fig. 5d). In fact, this layer expands faster than
the local escape velocity (Fig. 4a). This demonstrates that the energy
released from hydrogen recombination provides an important acceleration
mechanism in expanding and ejecting the envelope. Even though the
recombination spatially takes place close to the photosphere (with $\tau\sim
1$), it is still in the optically thick region where the recombination photons
cannot freely escape666At an optical depth of $\tau=10$, only
$e^{-10}\approx\qty{0.005}{\%}$ of the photon are expected to escape without
any scattering event, compared to $\,\approx\qty{37}{\%}$ at the photosphere
around $\tau\sim 1$. Initially, $R_{\tau=1}$ and $R_{\tau=10}$ are spatially
separated by $\lesssim\qty{10}{{\rm R}_{\odot}}$ which increases to
$\gtrsim\qty{1000}{{\rm R}_{\odot}}$ around $t=\qty{1500}{d}$. This
illustrates the significant difference between $\tau=1$ and $\tau=10$. It is
worth noting that in 3D simulations, this transition region is difficult to
resolve. (also see Ivanova & Nandez 2016; Clayton et al. 2017). On top of this
fast layer at $\qty{700}{d}\leq t\leq\qty{1000}{d}$, there is a slower-moving
layer that prevents the fast-expanding layer from becoming unbound and being
removed. As the fast material crashes in the slower layer on top, a shock wave
is produced, seen by the large negative velocity divergence
$\mathrm{d}v/\mathrm{d}r$ in Fig. 5b.

Layers exceeding $v_{\mathrm{esc}}$ reach the outer simulation boundary at
$t=\qty{1000}{d}$. Then, a rapid mass-loss event removes
${\sim}\,\qty{30}{\%}$ of the envelope mass in a short period of time,
followed by continuous steady envelope ejection. The sudden mass-loss event
can also be seen in Fig. 2, as the mass fraction of ejected envelope
$f_{\mathrm{ej}}$ increases rapidly at $t=\qty{1000}{d}$.

We find that $R_{\tau=10}$ stays constant for $700-\qty{1200}{d}$, after which
$R_{\tau=10}$ decreases as the hydrogen recombination front moves to smaller
radii. At about $t=\qty{1400}{d}$, the singly ionized helium-recombination
front is located at the hydrogen-recombination radius, which means that
hydrogen and singly ionized helium recombine at approximately the same
physical location. The released recombination energy accelerates the envelope
material, which can be seen by the increase in the slope of the lines of
constant mass in Fig. 4.

Along the hydrogen-recombination radius, there is a layer of positive velocity
divergence $\mathrm{d}v/\mathrm{d}r$ (Fig. 5b). This is another indicator that
the energy from hydrogen recombination contributes to accelerate the envelope.
The recombination of singly and doubly-ionized helium also causes layers of
positive velocity divergence but with a smaller magnitude compared to the
recombination of hydrogen. Because the partial ionization zones of singly and
doubly-ionized helium are more radially extended than the partial ionization
zone of hydrogen, the increase in the velocity divergence is less pronounced.
Once the upper heating radius decreases at $t=\qty{700}{d}$, the velocity
divergence inside the heating zone, i.e., all layers with
$r_{\mathrm{min}}^{\mathrm{heat}}\leq r\leq r_{\mathrm{max}}^{\mathrm{heat}}$,
increases significantly (Fig. 5b). This is an immediate consequence of the
increase in the specific heating rate as a result of the decrease in the total
envelope mass that is heated (Fig. 4f). Therefore, the localized heating
source causes the envelope to expand rapidly. At the same time, there is an
outward-traveling feature defined by a lower velocity divergence compared to
the surroundings (Fig. 5b). This feature is launched at $t=\qty{750}{days}$ at
the top of the heating zone and reaches the outer boundary at
$t=\qty{1000}{d}$. The cause of this feature, possibly a pressure/acoustic
wave, is likely connected to the decrease in the size of the heating layers.

The neutral layers outside the hydrogen-recombination radius expand
supersonically (Fig. 5a). The sound speed of the outer layers is significantly
lower due to adiabatic and photon cooling in the neutral layers (Fig. 5c).
These outward-moving layers show alternating positive and negative velocity
divergence (Fig. 5b), i.e., there are alternating layers with faster and
slower expansion velocities compared to the average expansion velocity of the
envelope. The details of the origin of this pattern are unclear. The large
negative velocity divergence at the end of the simulation between
$200-\qty{300}{{\rm R}_{\odot}}$ indicates a shock wave, but we cannot observe
any direct consequences arising from the shock.

The white hatching in Fig. 5c shows layers in the CE that are unstable against
convection. Initially, the entire envelope is convective, i.e., the entropy
gradient is less than or equal to zero. Heating stops the convective energy
transport after about 200d, at the same time as the envelope begins to expand.
For $t<\qty{200}{d}$, the injected energy is transported almost
instantaneously to the outer boundary, without affecting the entropy structure
of the envelope. As the heating rate increases with time (Fig. 4f) because of
the increase in the drag force, the injected energy stops convection.
Consequently, the energy cannot be transported away by convection but rather
leads to an expansion of the envelope. As energy is injected into the heating
zone, the specific entropy in the heating zone increases. The gradient of the
heating rate $\dot{\epsilon}_{\mathrm{heat}}$ is positive between
$r_{\mathrm{min}}^{\mathrm{heat}}$ and $a_{\mathrm{orb}}$ (Eq. 15). Hence, the
entropy gradient is also expected to be positive for
$r_{\mathrm{min}}^{\mathrm{heat}}\leq r\leq a_{\mathrm{orb}}$, implying that
these layers become stable against convection. While for $a_{\mathrm{orb}}\leq
r\leq r_{\mathrm{max}}^{\mathrm{heat}}$ the heating rate decreases (Eq. 15)
and most of these layers are also stable against convection for
$t>\qty{750}{d}$. They have a positive entropy gradient because the
heat/entropy of the inner layers is transported outwards. Outside the
hydrogen-recombination radius, the transport of energy by radiation is more
efficient than convection due to the lower opacity. Some smaller convective
zones are visible close to the companion and the outer boundary of the CE at
the end stages of the simulations. They are not expected to affect the
evolution of the CE because of their limited radial extent.

At the end of the simulation, the opacity close to the outer boundary of the
CE increases (Fig. 5f). This increase is caused by the appearance of hydrogen
molecules that can form at such low temperatures and densities. At these
temperatures, dust formation close to the outer boundary might be possible as
well, but it is not included in our simulations (dust formation is expected at
a temperature ${\sim}\,1200-\qty{2100}{K}$; Iaconi et al. 2020).

###  3.4 Recombination energy

Figure 6: Fraction of released recombination energy $f_{\mathrm{rec}}$ and
contributions of hydrogen (H ii), singly ionized helium (He ii) and doubly
ionized helium (He iii) for the simulation with $q=0.25$,
$C_{\mathrm{d}}=0.23$ and $C_{\mathrm{h}}=4.0$. The total fraction of released
recombination energy is shown in black. The horizontal colored-dashed lines
show the initially available recombination energy in each species. The
vertical gray-dashed line indicates the end of the dynamical plunge-in phase
as determined by Eq. (22).

To study the effects of recombination on envelope ejection, we calculate the
amount of released recombination energy. For this, the ionization potentials
of hydrogen (H ii), singly-ionized helium (He ii) and doubly-ionized helium
(He iii) are taken from Kramida et al. (2021),

$\displaystyle E_{H\textsc{ii}}$ $\displaystyle\equiv
E_{{H\textsc{i}}\rightarrow{H\textsc{ii}}}=13.5984\,\mathrm{eV}$ (24)
$\displaystyle E_{{He\textsc{ii}}}$ $\displaystyle\equiv
E_{{He\textsc{i}}\rightarrow{He\textsc{ii}}}=24.5874\,\mathrm{eV}$ (25)
$\displaystyle E_{{He\textsc{iii}}}$ $\displaystyle\equiv
E_{{He\textsc{ii}}\rightarrow{He\textsc{iii}}}=54.4178\,\mathrm{eV}$ (26)

and the atomic masses from Prohaska et al. (2022),

$\displaystyle m_{\mathrm{H}}$ $\displaystyle=1.0080\,\mathrm{amu}$ (27)
$\displaystyle m_{\mathrm{He}}$ $\displaystyle=4.0026\,\mathrm{amu}.$ (28)

The potentially available energy from the recombination of hydrogen and
helium, usually referred to as the recombination energy, is then given by

$E_{\mathrm{rec}}=\int\limits_{M_{\mathrm{1,core}}}^{M_{1}}\left[\frac{X}{m_{\mathrm{H}}}E_{{H\textsc{ii}}}f_{{H\textsc{ii}}}\right.\\\
\left.+\frac{Y}{m_{\mathrm{He}}}\Big{(}E_{{He\textsc{ii}}}f_{{He\textsc{ii}}}+(E_{{He\textsc{iii}}}+E_{{He\textsc{ii}}})f_{{He\textsc{iii}}}\Big{)}\right]\mathrm{d}m,$
(29)

where $X$ and $Y$ are the mass fractions of hydrogen and helium, and
$f_{{H\textsc{ii}}}$, $f_{{He\textsc{ii}}}$ and $f_{{He\textsc{iii}}}$ are the
ionization fractions of H ii, He ii and He iii respectively777The envelope
material that is removed during the simulation has fully recombined before
reaching velocities larger than the escape velocity (Fig. 4). Thus, the
envelope ejection does not remove any potential recombination energy from the
system.. The recombination energies of elements more massive than helium are
not taken into account, because their contribution to the total recombination
energy is negligible (Ivanova & Nandez 2016).

In Fig. 6, we show the fractions $f_{\mathrm{rec,total}}$ and
$f_{\mathrm{rec},i}$ of released recombination energy

$\displaystyle f_{\mathrm{rec,total}}$
$\displaystyle=\frac{E_{\mathrm{rec,ini}}-E_{\mathrm{rec}}}{E_{\mathrm{rec,ini}}},$
(30) $\displaystyle f_{\mathrm{rec},i}$
$\displaystyle=\frac{E_{\mathrm{rec,ini},i}-E_{\mathrm{rec},i}}{E_{\mathrm{rec,ini}}}\;\mathrm{with}\;i={{H\textsc{ii}},\,{He\textsc{ii}},\,{He\textsc{iii}}},$
(31)

where $E_{\mathrm{rec,i}}$ is computed by integrating each term in Eq. (29)
individually.

The released recombination energy is negligible for $t<\qty{700}{d}$. This can
be understood with the help of Fig. 5, which shows that the partial ionization
zones of hydrogen and helium remain at constant mass coordinates during this
time. At the end of the simulation, more than $98\,\%$ of the initially
available recombination energy is released, and about $70\,\%$ of the released
energy is due to hydrogen recombination. The reason for this is the high
abundance of hydrogen in the envelope ($X_{\mathrm{ini}}=0.7$). For
$\qty{600}{d}\leq t\leq\qty{750}{d}$, we find
$f_{\mathrm{rec},{He\textsc{ii}}}<0$, because doubly-ionized helium recombines
to form singly-ionized helium such that the recombination energy available
from singly-ionized helium is larger than the initial recombination energy
stored in singly-ionized helium.

Our calculations show that recombination energy is likely to provide an
important contribution to the envelope ejection process. Hydrogen
recombination occurs well below the photosphere (at a typical optical depth
$\tau=10$), where photons cannot escape directly. This results in a positive
velocity divergence at the location where hydrogen recombines, causing further
acceleration of the envelope (Fig. 5b). The recombination front of singly-
ionized helium is deeper in the envelope compared to the recombination front
of hydrogen. Around this recombination zone, the velocity divergence is
positive. It is not as large as in the case of hydrogen recombination but much
more extended radially. A similar behavior can be observed for the
recombination zone of doubly-ionized helium. There is a layer of negative
velocity divergence around the He iii recombination zone for
$t>\qty{1300}{d}$. As most of the doubly-ionized helium has already recombined
at this point (see Fig. 6), the energy released might not be enough to further
accelerate the envelope.

At the end of the plunge-in phase, about $50\,\%$ of the available
recombination energy is released, while only ${\sim}\,40\,\%$ of the envelope
is ejected. Hereafter, the energy that drives the envelope ejection is mostly
from recombination and no longer from the heating, as by definition the
orbital separation decreases slowly in the post-plunge-in phase which means
that the drag force and hence the heating are much lower than during the
plunge-in phase (see also Sects. 3.5 and 3.6). At the end of the simulation,
more than $70\,\%$ of the envelope is ejected, which suggests that the
recombination energy contributes significantly to the ejection, especially in
the post-plunge-in phase. For a comparison, at the end of the 3D simulation at
$t=\qty{2500}{d}$, $\qty{91}{\%}$ of the envelope mass is ejected (Sand et al.
2020).

###  3.5 Energy budget

Figure 7: Major energy sources (dashed lines) and sinks (full lines) during
the CE simulation with $q=0.25$, $C_{\mathrm{d}}=0.23$ and
$C_{\mathrm{h}}=4.0$. Shown are the total injected energy $\Delta
E_{\mathrm{heat}}$, the change in internal energy $\Delta E_{\mathrm{int}}$,
the total released recombination energy $\Delta E_{\mathrm{rec}}$ (the
recombination energy is part of the internal energy and is only shown
separately to highlight the contribution from recombination), the change in
orbital energy $\Delta E_{\mathrm{orb}}$ according to (34), the extra term
$\Delta E_{m}$ defined by Yarza et al. (2022), the change in the gravitational
energy $\Delta E_{\mathrm{grav}}$, and the energy that is carried away by
radiation $\Delta E_{\mathrm{rad}}$. The vertical gray-dashed line indicates
the end of the dynamical plunge-in phase as determined by Eq. (22). The
energies are shown in units of the initial binding energy of the envelope with
$E_{\mathrm{bind,ini}}=\qty{-4.03e45}{erg}$999The definition of the core mass
$M_{\mathrm{1,core}}$ is arbitrary and can affect the value of the binding
energy hugely (Tauris & Dewi 2001). For consistency, we use the definition of
$M_{\mathrm{1,core}}=\qty{0.545}{{\rm M}_{\odot}}$ of Sand et al. (2020)
throughout this paper. When using $M_{\mathrm{1,core}}=\qty{0.540}{{\rm
M}_{\odot}}$, the helium core mass defined by a hydrogen mass-fraction less
than 0.1, the binding energy of the envelope increases to
$\qty{-5.67e46}{erg}$. Hence, the fractions $\Delta E/|E_{\mathrm{bind,ini}}|$
in Fig. 9 should not be used for absolute comparisons.. The lines showing
contributions from $\Delta E_{\mathrm{orb}}$ and $\Delta E_{\mathrm{m}}$ are
dash-dotted, as these show source terms for the heating $\Delta
E_{\mathrm{heat}}$, but are not direct energy sources in the MESA simulation.

Major energy sources and sinks during the 1D CE simulation are shown in Fig.
9. The energies are shown in units of the initial binding energy of the
envelope

$E_{\mathrm{bind}}=\int_{M_{1,\mathrm{core}}}^{M_{1}}\left(-\frac{Gm}{r}+u\right)\,\mathrm{d}m,$
(32)

where $u$ is the internal energy. The total energy injected via heat is
determined from the drag force and the relative velocity by

$\Delta
E_{\mathrm{heat}}=-\int\boldsymbol{F}_{\mathrm{d}}\cdot\boldsymbol{v}_{\mathrm{rel}}\,\mathrm{d}t.$
(33)

The source of the heating energy is the orbital energy of the two stars in the
CE phase. The released orbital energy is classically given by

$\Delta
E_{\mathrm{orb}}=-\frac{GM_{1,a_{\mathrm{orb}}}M_{2}}{2\,a_{\mathrm{orb}}}+\dfrac{GM_{1}M_{2}}{2\,a_{\mathrm{orb,ini}}},$
(34)

where $M_{1,a_{\mathrm{orb}}}$ is the mass of $M_{1}$ enclosed by the orbit
with separation $a_{\mathrm{orb}}$. From Fig. 9, it is clear that $\Delta
E_{\mathrm{heat}}$ is approximately a factor 2 larger than $-\Delta
E_{\mathrm{orb}}$. This discrepancy is expected in our CE formalism because we
mix both 1D and 3D treatments of physical mechanisms. For example, we evolve
the binary orbit in 3D but are bound to use 1D approximations to calculate the
change in potential energy. Therefore, we cannot expect to conserve energy in
our CE formalism, i.e., we do expect $\Delta E_{\mathrm{heat}}\neq-\Delta
E_{\mathrm{orb}}$.

In addition, Eq. (34) only applies to systems, where the envelope is
completely ejected, and the ejecta has zero velocity at infinity. In our
simulations, the above-mentioned criteria are not fulfilled, and, therefore,
Eq. (34) is insufficient to capture all the relevant effects. This means, that
there is a need for extra terms in Eq. (34) to correctly calculate $\Delta
E_{\mathrm{orb}}$ in 1D CE treatments similar to ours. One such term arises
from the change in potential energy because of the change in the enclosed
mass. In the limit $M_{2}\ll M_{1,a_{\mathrm{orb}}}$, Yarza et al. (2022)
suggest this term to be

$\Delta
E_{m}=GM_{2}\int_{a_{\mathrm{orb,ini}}}^{a_{\mathrm{orb}}}\frac{1}{r}\frac{\mathrm{d}M_{1,r}}{\mathrm{d}r}\,\mathrm{d}r,$
(35)

where $M_{1,r}$ is the mass of the giant star enclosed by the radius $r$. We
test the effect of $\Delta E_{m}$ to Eq. (34) on our simulation, to see if it
can resolve the observed mismatch between $\Delta E_{\mathrm{heat}}$ and
$-\Delta E_{\mathrm{orb}}$. Even when considering the contribution of $\Delta
E_{m}$ to the change in orbital energy, there is more energy injected in the
envelope than is released (Fig. 9). We note that in our simulation the
condition $M_{2}\ll M_{1,a_{\mathrm{orb}}}$ is not fulfilled. Further reasons
for the discrepancy are that the change in the orbital energy is not only
caused by the change in the orbital separation but also by the envelope
expansion. As the envelope expands, the mass enclosed by the orbit
$M_{1,a_{\mathrm{orb}}}$ decreases and thus the potential energy of the
companion. This means that the orbital energy of a companion at a constant
$a_{\mathrm{orb}}$ inside an expanding envelope increases because of the
decrease in the enclosed mass by the orbit. This effect is not included in the
term $\Delta E_{m}$ proposed by Yarza et al. (2022).

Additionally, the energy loss via radiation at the outer boundary of the
envelope is shown in Fig. 9 as well as the energy release from the
recombination of hydrogen and helium. The losses from radiation are about half
of the energy injected via heat. Both the losses via radiation and the energy
gain from recombination start to become significant for $t>\qty{700}{d}$,
after which they increase almost synchronously. This does not mean that the
energy released from recombination is immediately radiated away, but rather
contributes to accelerate the envelope material, e.g., seen by the kink in the
lines of constant mass at the hydrogen recombination front (Fig. 4b) as well
as by the positive velocity divergence (Fig. 5b). This acceleration can even
cause the layers to expand faster than the local escape velocity (Fig. 4a).

###  3.6 Drag force evolution

Figure 8: Time evolution of the drag force and its components for $q=0.25$,
$C_{\mathrm{d}}=0.23$ and $C_{\mathrm{h}}=4.0$. The visualization of these
quantities represents how each of them enters the drag force according to Eq.
(5). During the whole simulation, the drag force is in the non-linear regime,
i.e., $\mathcal{M}>1.01$ and $\eta_{\mathcal{B}}>0.1$. Additionally, the ratio
of the drag force and gravitational force $F_{\mathrm{g}}$ is shown. The
vertical gray-dashed line indicates the end of the dynamical plunge-in phase
as determined by Eq. (22).

The drag force acting on the companion is calculated by Eq. (5) following the
results of Kim (2010) and Kim & Kim (2007). In Fig. 8, we show the individual
components that enter the drag force as well as the ratio of the drag force
and the gravitational force $F_{\mathrm{g}}$. The drag force increases for
$t<\qty{750}{d}$. From Fig. 8 it is apparent, that the drag force is mostly
influenced by the density $\rho$ and the relative velocity $v_{\mathrm{rel}}$
as they vary the most compared to the other components of the drag force. For
$t<\qty{750}{d}$, the density increases which causes the drag force to
increase. The relative velocity increases throughout the entire simulation
because the orbital velocity increases with ${\sim}\,a_{\mathrm{orb}}^{-1/2}$
and the rotational velocity decreases with ${\sim}\,a_{\mathrm{orb}}$ as the
orbit separation $a_{\mathrm{orb}}$ decreases (compare Eq. (7)). Hence, the
drag force, which is proportional to $v_{\mathrm{rel}}^{-2}$, decreases. For
$t<\qty{750}{d}$, the increase of $\rho$ dominates over the decrease of
$v_{\mathrm{rel}}^{-2}$, causing the drag force to increase. For
$t>\qty{750}{d}$, the drag force decreases, because both the $\rho$ and
$v_{\mathrm{rel}}^{-2}$ decrease.

The ratio $F_{\mathrm{d}}/F_{\mathrm{g}}$ is an indicator for the relative
strength of the drag force compared to the gravitational force, e.g., a larger
value of $F_{\mathrm{d}}/F_{\mathrm{g}}$ shows a high drag force and hence we
expect a larger orbital decay compared to a lower value of
$F_{\mathrm{d}}/F_{\mathrm{g}}$. At $t=\qty{250}{d}$,
$F_{\mathrm{d}}/F_{\mathrm{g}}$ starts to decrease, indicating that the
relative strength of the drag force decreases (Fig. 8). This is at
approximately at same time as the envelope starts to expand (Fig. 4).
Initially, we find values for $F_{\mathrm{d}}/F_{\mathrm{g}}$ between
$10^{-1.8}$ and $10^{-1.5}$. At $t=\qty{750}{d}$, the drag force peaks and
then decreases. This change in behavior of $F_{\mathrm{d}}$ translates into
the change in slope of $F_{\mathrm{d}}/F_{\mathrm{g}}$. At the end of the
plunge-in phase at $t=\qty{1000}{d}$, we find
$F_{\mathrm{d}}/F_{\mathrm{g}}=\qty{e-3}{}$, which decreases to e-4 at the end
of the simulation. This shows, that indeed the drag force becomes small
compared to gravity during the post-plunge-in phase, which is the main reason,
why the dynamical plunge-in ends.

##  4 Results for mass ratios $q=0.50$ and $q=0.75$

Figure 9: Similar to Fig. 2, but showing the best fitting 1D simulation
compared to the 3D simulation of Sand et al. (2020) for mass ratios $q=0.5$ in
panel (a) and $q=0.75$ in panel (b).

The spiral-in and envelope-ejection curves for our 1D simulations with mass
ratios $q=0.5$ and $q=0.75$ are shown in Fig. 9. For a mass ratio of $q=0.5$,
the spiral-in curve as well as the mass fraction of the ejected envelope
capture well the results of the 3D simulation of Sand et al. (2020). Similar
to the case with mass ratio $q=0.25$, the envelope-ejection rate in the later
phases of the simulation is very similar to the 3D simulation. Additionally,
both the 1D and the 3D simulation show a comparable post-plunge-in
eccentricity of 0.020 and 0.017, respectively (see Table 1).

For a mass ratio of $q=0.75$, however, it is more difficult to find values of
$C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ such that the 1D simulation reproduces
well the spiral-in and the envelope ejection of the 3D simulation. Both the
post-plunge-in orbital separation and the final envelope-ejection rate are
lower than in the 3D simulation. This suggests that we have reached the
limitations of our model as the companion for mass ratio $q=0.75$ is no longer
a small perturber, but rather comparable in mass to the initial AGB star and
even more massive than its helium core of 0.54M_⊙. The eccentricity of the
post-plunge-in orbit decreases after about $1250\,\mathrm{d}$ (Fig. 9) because
parts of the envelope material fall back below the orbit of the companion,
causing a spike in the drag force that, in turn, alters the orbit.

Table 1: Summary of the fitting parameters of the three 1D CE simulations. The orbital separations at the end of the plunge-in phase $a_{\mathrm{orb,pl}}$ are given for both the 1D and the 3D simulations. The eccentricities $e_{\mathrm{pl}}$ are determined for 1D and 3D simulations via Eq. (23) at the end of the plunge-in phase. Both $a_{\mathrm{orb,pl}}$ and $e_{\mathrm{pl}}$ are averaged over 5 orbits. $q$ | $C_{\mathrm{d}}$ | $C_{\mathrm{h}}$ | $a_{\mathrm{orb,pl}}^{\mathrm{1D}}$ | $a_{\mathrm{orb,pl}}^{\mathrm{3D}}$ | $e_{\mathrm{pl}}^{\mathrm{1D}}$ | $e_{\mathrm{pl}}^{\mathrm{3D}}$
---|---|---|---|---|---|---
|  |  | $\left[\mathrm{R}_{\odot}\right]$ | $\left[\mathrm{R}_{\odot}\right]$ |  |
0.25 | 0.23 | 4.00 | 21.2 | 25.2 | 0.005 | 0.006
0.50 | 0.30 | 1.30 | 40.1 | 47.6 | 0.020 | 0.017
0.75 | 0.30 | 0.75 | 57.2 | 76.7 | 0.034 | 0.018

The best-fit parameters $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$ as well as the
post-plunge-in orbital separations and eccentricities of the three simulations
with mass ratios $q=0.25$, 0.5 and 0.75 are summarized in Table 1. In all
cases, the post-plunge-in orbital separation of the 1D simulations is smaller
than that of the 3D simulation. Possible reasons are discussed in Sect. 5. For
$q=0.25$ and 0.5, the eccentricity of both the 1D and the 3D simulations are
comparable while for $q=0.75$ the eccentricity in the 1D simulations is twice
as large as in the 3D simulation. The drag-force parameter $C_{\mathrm{d}}$
stays almost the same for all simulations with $C_{\mathrm{d}}{\sim}\,0.25$.
The heating parameter $C_{\mathrm{h}}$ decreases with increasing mass ratio
$q$. However, the results for a mass ratio $q=0.75$ may not be accurate, as
the simulation does not fit the 3D simulation as well as for the lower mass
ratios.

The results for $q=0.25$ discussed in Sect. 3 regarding the role of the
recombination energy in ejecting the envelope, the behavior of the drag force
and the energy budget qualitatively also apply for the simulations with higher
mass ratios. We find that the recombination energy is important in driving the
envelope ejection. The behavior of the drag force is mostly determined by the
change in relative velocity and density. Both of which cause the drag force to
drop at the end of the plunge-in phase.

##  5 Discussion

In this section, we discuss our 1D CE method as well as the results that we
described above. First, we show the limitations as well as the advantages of
our 1D CE method in Sects. 5.1 and 5.2 respectively. Then, we compare our 1D
CE method to other proposed 1D methods for simulations of the CE phase within
the context of 3D simulations (Sect. 5.3). Finally, we physically motivate the
necessity of the two free parameters used in our model (Sect. 5.4).

###  5.1 Limitations of the 1D CE model

In our 1D CE model, the drag force acts only on the companion and not on the
core of the giant star (Eqs. 1 and 2). For low mass ratios, this assumption is
valid, as the center of mass (CM) of the binary is close to the core of the
giant star. For higher mass ratios, the CM is located in between the companion
and the core of the giant star. Hence, both the core of the giant star and the
companion are orbiting inside the CE around the CM and experience a drag
because of dynamical friction. To some extent, this effect is compensated by
using a free parameter for the drag-force calculations. For our 1D CE
simulation with a mass ratio of 0.75, it is no longer valid to use this
assumption, because the mass of the companion ($0.73\,{\rm M}_{\odot}$) is
larger than the mass of the helium core of the AGB star ($0.54\,{\rm
M}_{\odot}$), i.e., the CM is located closer to the companion than to the core
of the AGB star.

Additionally, the CE is simulated as the perturbed envelope of a giant single
star. Similarly to the argument above, the CE is not centered on the core of
the giant star, but rather on the CM of the binary. Therefore, the pressure,
density, and sound speed at the location of the companion are different from
the predictions in our model, which consequently causes the drag force to be
different as well. Again, for low mass ratios, this assumption is valid,
because the CM is located closer to the core of the giant star. However, this
simplification breaks down for larger mass ratios.

In 3D simulations of CE events similar to the simulations in Sand et al.
(2020), spiral arms emerge as the companion plunges into the CE and the
companion and the core of the giant star orbit each other inside the CE. There
are usually two spiral arms observed, one originating at the companion at the
other originating at the core of the giant star. These spiral arms might
transport orbital energy from the companion and the giant’s core to the CE as
they induce shocks in the envelope. This hypothesis needs further testing
using 3D CE simulations. If true, an energy transport mechanism where the
orbital energy is transferred to the envelope via shocks in the spiral arms
would imply that the entire envelope could be heated in a 1D CE model.
Additionally, there might be a time delay as the energy needs to be
transported via the spiral arms first, before thermalizing in the envelope.
Within our 1D model, we assume that all the released energy from the orbit is
converted into heat. However, some of the orbital energy might also be
converted into rotational energy, i.e., spinning up the envelope. We ignore
this effect and only allow energy conversion into heat. Ivanova & Nandez
(2016) argue that in 1D CE simulations, energy should not be added as heat but
rather as kinetic energy. The spiral arms which are present in the simulations
in Sand et al. (2020) convert some of the kinetic energy into thermal energy.
Therefore, it seems reasonable to inject heat into the envelope to effectively
simulate the response of the envelope to the spiral-in instead of kinetic
energy.

During the early plunge-in phase, there is a difference in the envelope
ejection mechanism between our 1D CE method and the 3D CE simulations in Sand
et al. (2020). In the 1D model, a large part of the envelope is slowly heated
up. The injected internal energy is converted into kinetic energy and almost
the entire envelope starts to expand. It takes a certain time until the
surface of the envelope exceeds the escape velocity and becomes unbound. In
the 3D simulations, the outer layers are dynamically flung out and become
immediately unbound, i.e., orbital energy is directly converted to kinetic
energy rather than heat. As our 1D model does not include energy injection as
kinetic energy, the initial differences in the mass fractions of the ejected
envelope are expected (see Figs. 2 and 9).

To account for the mass of the companion inside the envelope, we modified the
gravitational constant outside the orbit of the companion (Sect. 2.5). This is
equivalent to modeling the companion as a thin shell with radius
$a_{\mathrm{orb}}$, as such a thin shell has a constant potential for
$r<a_{\mathrm{orb}}$ and a ${\sim}\,1/r$ potential for $r>a_{\mathrm{orb}}$.
When comparing this approximation to 3D CE simulations, Ivanova & Nandez
(2016) find a deeper potential in 3D simulations close to the orbit and a more
shallow potential far outside, with a relative deviation of up to $50\,\%$
between the 1D potential with the thin shell approximation and the 3D
simulations for a similar CE configuration. This approximation could be
improved with a different potential form for the companion, but it works well
to first order. Inconsistencies might also arise as we sometimes model the
companion as a point particle, e.g., for evaluating the drag force and
integrating the orbits (Sect. 2.4), and sometimes as a shell, e.g., for
modeling the influence of the companion on the structure of the CE (Sect.
2.5).

###  5.2 Advantages of the 1D CE model

The main reason for creating a 1D CE model is to reduce the computational cost
of simulations, which enables larger parameter studies. Our 1D CE simulations
take ${\sim}\,10$ core-hours to simulate a physical time of 2000d, while the
computational cost for the 3D simulation is of the order of $10^{5}$ core-
hours. This means that the 1D CE simulation can be run easily on a desktop
computer and does not rely on high-performance computing facilities.
Therefore, a larger parameter space of systems that evolve through a CE phase
can be explored with such 1D models.

In the simulation of Sand et al. (2020), the core of the giant star is cut out
and replaced by a point mass to allow for sufficiently large time-steps
(Ohlmann et al. 2016a; Sand et al. 2020). The cut-out core of the giant star
cannot respond to the dynamic changes in the envelope. In Sand et al. (2020),
the core is cut at $5\,\%$ of the initial radius of the giant star, i.e., the
cut-out core is larger than the helium core and contains hydrogen-rich layers.
As the CE expands and is subsequently ejected, the core in the 3D simulation
is expected to expand on the dynamical timescale to react to the loss of
pressure by the CE. Our 1D simulations can resolve the core throughout the
whole CE simulation and can capture the expansion of the layers that are
considered to be inside the excised core in the 3D simulations. We find that
the mass of the core, as defined in Sand et al. (2020), decreases by
${\sim}\,3\,\%$ during our 1D CE simulations because the core expands. If the
core is allowed to expand, we expect a deeper spiral-in as more envelope
material needs to be ejected. This is a possible explanation for why we
consistently find lower post-plunge-in separations compared to the 3D models
Sand et al. (2020).

In addition, the 1D simulations in MESA include energy transport via photons
and photosphere cooling which is more difficult to implement 3D simulations.
This also means that the release and transport of the recombination energy is
better handled in the 1D model, which is another possible reason why we find
lower post-plunge-in separations compared to the 3D simulations.

###  5.3 Comparison to other 1D and 3D CE simulations

Several 1D CE models were proposed in the past to simulate CEs. The model of
Fragos et al. (2019) is similar to our model in the sense that they assume a
drag force acting on the companion, which removes orbital energy that is then
injected as heat into the envelope. In contrast to our model, they assume a
circular orbit for the companion and the released energy is injected into a
region defined by one local accretion radius, i.e., heating within
$a_{\mathrm{orb}}\pm R_{\mathrm{a}}(a_{\mathrm{orb}})$ using the terminology
of our model. We find that the extent of the heating zone directly affects the
envelope ejection and we need to tune the extent of the heated zone via
$C_{\mathrm{h}}$ to obtain the same amount of envelope ejection as in the 3D
simulations. A similar model compared to Fragos et al. (2019) is used by
O’Connor et al. (2023) to study the envelope’s response to a planetary
engulfment. In our model, we have the advantage of comparing the predictions
from the 1D simulation to 3D simulations of the same initial CE setup and
tuning the free parameters of our model accordingly.

Sand et al. (2020) find in their 3D CE simulations that the energy from
recombination is necessary to eject significant amounts of the envelope during
the simulation by comparing different equations of state, which either include
or do not include recombination energy. In our 1D simulations, we also find
that the envelope ejection is only triggered once hydrogen recombines. A
significant fraction of the recombination energy accelerates the envelope,
which is subsequently ejected. Once recombination starts, the hydrogen
recombination front stays at a constant radius ($\qty{700}{d}\leq
t\leq\qty{1300}{d}$ in Fig. 4b). This process is similar to the _steady
recombination outflow_ described by Ivanova & Nandez (2016). The recombination
energy expands and ejects the envelope. Therefore, the inner layers adjust to
the decrease in pressure by expanding. As these inner layers expand, they cool
down until recombination expands them further, causing the even deeper layers
to expand, and so on. Ivanova & Nandez (2016) argue that this leads to the
recombination at a constant radial coordinate. In the late stages of the
simulations, we find, however, that the radius of the recombination front of
hydrogen and also helium decreases. This might be caused by the envelope
running out of material, since more than $97\,\%$ of the initial envelope
recombined at the end of the simulations.

###  5.4 Physical motivation for two free parameters

We use the drag force models of Kim (2010) and Kim & Kim (2007) for all of our
1D CE simulations. Both models assume a perturber moving on a circular orbit
through a gaseous background medium. While Kim & Kim (2007) assume a low-mass
perturber, the model of Kim (2010) applies also to high-mass perturbers. In
both cases, the background medium is assumed to be homogeneous in density. Kim
(2010) argue that density gradients can be ignored if the Bondi radius is much
smaller than the pressure scale height, i.e.,
$\mathcal{M}\ll(M_{1,a_{\mathrm{orb}}}/M_{2})^{1/4}$. In our 1D CE
simulations, we find that this condition is not always fulfilled, especially
during the plunge-in phase. Additionally, Kim (2010) assumes that the
centrifugal force and the Coriolis force can be ignored. This assumption is
valid if $\mathcal{B}/\mathcal{M}^{2}\ll 1$. We also find that this condition
is not satisfied during our simulations. Because our simulations disagree with
these assumptions, we do not expect the drag force given by the models of Kim
(2010) and Kim & Kim (2007) to provide the correct values for our simulation
without considering a calibration factor. Hence, it seems appropriate to
introduce a calibration factor $C_{\mathrm{d}}$ for the drag force, which we
adjust by comparing our 1D to the 3D CE simulations. From the three
comparisons between 1D and 3D simulations, we find that $C_{d}\approx 0.25$
(Table 1).

The Bondi-Lyttleton-Hoyle model for accretion (Hoyle & Lyttleton 1939; Bondi &
Hoyle 1944) also assumes a homogeneous background density. We use the
accretion radius of the Bondi-Lyttleton-Hoyle model to estimate the heating
zone, i.e., the layers where the material is gravitationally deflected and
focused by the companion. Because we assume spherical symmetry in the 1D CE
model, we always heat a spherical shell defined by the accretion radius. In a
real CE scenario, only a region around the companion might be gravitationally
affected, and not a spherical shell centered on the core of the giant star.
Introducing a second calibration factor $C_{\mathrm{h}}$ for the accretion
radius and therefore also for the size of the heating zone seems reasonable.
From the results of our simulations, we find the trend of decreasing
$C_{\mathrm{h}}$ for increasing $q$ (Table 1). We also analyzed the quantities
$qR_{\mathrm{a}}(a_{\mathrm{orb}})/a_{\mathrm{orb}}$ and
$q(r^{\mathrm{heat}}_{\mathrm{max}}-r^{\mathrm{heat}}_{\mathrm{min}})/a_{\mathrm{orb}}$
to see if we can find a calibration for the extent of the heating zone. For
the three 1D models we computed, we could not see a better calibration with
these quantities. Therefore, we kept using the heating parameter
$C_{\mathrm{h}}$ as a free parameter in our model.

From the simulations discussed in Sect. 3.2, it is necessary to have two free
parameters in our 1D CE model. If we restrict our model to only one free
parameter, we cannot reproduce the results of the 3D CE simulations of Sand et
al. (2020). Many approximations of our 1D model are captured by the free
parameters $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$. Therefore, the best-fitting
values do not directly carry any physical meaning and cannot be immediately
used to constrain the physics during CE events.

##  6 Conclusions

We presented a 1D approach to simulate the CE evolution within the stellar
evolution code MESA using its hydrodynamic capabilities. The CE is modeled by
the envelope of a giant star, in which a point-mass companion is placed on a
circular orbit. The viscous forces hinder the motion of the companion because
of the dynamical drag by the CE. These extra forces are included in the
equations of motion of the companion via a parametric drag-force prescription,
allowing us to integrate the orbital evolution of the binary star. The energy
lost because of the drag force is added as heat to the envelope. We trace the
layers expanding faster than the escape velocity and remove them from the
simulated domain. To fit our results to the 3D CE simulations, we include two
free parameters, which we use as optimization parameters for the fits.

We simulate the CE phase of a 0.97M_⊙ AGB star and a point-mass companion with
mass ratios $q=0.25$, 0.50, and 0.75. When comparing and fitting these
simulations to the 3D CE simulations of Sand et al. (2020), we find the
following conclusions:

  * •

It is possible to reproduce the spiral-in curve and the mass fraction of the
ejected envelope when using both free parameters for $q=0.25$ and 0.50. With
one parameter alone, this is not possible. Therefore, 1D CE models similar to
ours probably require at least two free parameters to reproduce 3D CE
computations. We find that the spiral-in timescale is mostly determined by
$C_{\mathrm{d}}$, while the mass-fraction of the ejected envelope is
determined by both $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$.

  * •

We are unable to find satisfactory fits with the 1D CE model for $q=0.75$.
This might be an extreme case where the assumptions and approximations of our
1D CE model are no longer valid and do not result in a physical solution.

  * •

For all mass ratios, we find the post-plunge-in separation in the best-fitting
1D simulation to be smaller than those of the 3D simulation.

  * •

Regardless of $C_{\mathrm{d}}$ and $C_{\mathrm{h}}$, the post-plunge-in
separation is almost always similar and only depends on the mass ratio. This
suggests a deeper physical mechanism that determines the post-plunge-in
separation.

  * •

The recombination energy released from hydrogen and helium likely plays an
important role in accelerating and ejecting the envelope.

In a future study, we plan to simulate CE events with different initial
configurations, i.e., different masses, and evolutionary stages of the giant
star, and compare them to 3D CE simulations to investigate whether a global
calibration of the two free parameters based on the mass and evolutionary
stage of the giant star as well as the mass ratio is possible. If we can find
such a calibration, this CE model can be used to predict the outcome of CE
events, i.e., the post-plunge-in orbital separations and the ejecta mass if
the envelope is ejected successfully.

###### Acknowledgements.

We thank the anonymous referee for their feedback, which helped to improve the
quality of the paper. We thank C. Sand for providing the data for the
simulations described in Sand et al. (2020). VAB, FRNS, PhP and FKR
acknowledge support from the Klaus Tschira Foundation. This work has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (Grant agreement No. 945806).
This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948
(the Heidelberg STRUCTURES Excellence Cluster). VAB acknowledges support from
the International Max Planck Research School for Astronomy and Cosmic Physics
at the University of Heidelberg (IMPRS-HD).

## References

  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407
  * Blöcker (1995) Blöcker, T. 1995, A&A, 299, 755
  * Bondi & Hoyle (1944) Bondi, H. & Hoyle, F. 1944, MNRAS, 104, 273
  * Cash & Karp (1990) Cash, J. R. & Karp, A. H. 1990, ACM Trans. Math. Softw., 16, 201–222
  * Chamandy et al. (2018) Chamandy, L., Frank, A., Blackman, E. G., et al. 2018, MNRAS, 480, 1898
  * Clayton et al. (2017) Clayton, M., Podsiadlowski, P., Ivanova, N., & Justham, S. 2017, MNRAS, 470, 1788
  * Courant et al. (1928) Courant, R., Friedrichs, K., & Lewy, H. 1928, Mathematische Annalen, 100, 32
  * De Marco & Izzard (2017) De Marco, O. & Izzard, R. G. 2017, PASA, 34, e001
  * Detmers et al. (2008) Detmers, R. G., Langer, N., Podsiadlowski, P., & Izzard, R. G. 2008, A&A, 484, 831
  * Eldridge & Stanway (2016) Eldridge, J. J. & Stanway, E. R. 2016, MNRAS, 462, 3302
  * Fragos et al. (2019) Fragos, T., Andrews, J. J., Ramirez-Ruiz, E., et al. 2019, ApJ, 883, L45
  * Fryer & Woosley (1998) Fryer, C. L. & Woosley, S. E. 1998, ApJ, 502, L9
  * Glanz & Perets (2018) Glanz, H. & Perets, H. B. 2018, MNRAS, 478, L12
  * Grott et al. (2005) Grott, M., Chernigovski, S., & Glatzel, W. 2005, MNRAS, 360, 1532
  * Han et al. (1995) Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1995, MNRAS, 272, 800
  * Heber (2009) Heber, U. 2009, ARA&A, 47, 211
  * Hirai & Mandel (2022) Hirai, R. & Mandel, I. 2022, ApJ, 937, L42
  * Hoyle & Lyttleton (1939) Hoyle, F. & Lyttleton, R. A. 1939, Proceedings of the Cambridge Philosophical Society, 35, 405
  * Hoyle & Lyttleton (1941) Hoyle, F. & Lyttleton, R. A. 1941, MNRAS, 101, 227
  * Iaconi & De Marco (2019) Iaconi, R. & De Marco, O. 2019, MNRAS, 490, 2550
  * Iaconi et al. (2020) Iaconi, R., Maeda, K., Nozawa, T., De Marco, O., & Reichardt, T. 2020, MNRAS, 497, 3166
  * Iaconi et al. (2017) Iaconi, R., Reichardt, T., Staff, J., et al. 2017, MNRAS, 464, 4028
  * Iben & Tutukov (1984) Iben, I., J. & Tutukov, A. V. 1984, ApJS, 54, 335
  * Ivanova et al. (2013) Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59
  * Ivanova et al. (2015) Ivanova, N., Justham, S., & Podsiadlowski, P. 2015, MNRAS, 447, 2181
  * Ivanova & Nandez (2016) Ivanova, N. & Nandez, J. L. A. 2016, MNRAS, 462, 362
  * Izzard et al. (2004) Izzard, R. G., Ramirez-Ruiz, E., & Tout, C. A. 2004, MNRAS, 348, 1215
  * John & Draper (1980) John, J. A. & Draper, N. R. 1980, Journal of the Royal Statistical Society: Series C (Applied Statistics), 29, 190
  * Kim & Kim (2007) Kim, H. & Kim, W.-T. 2007, ApJ, 665, 432
  * Kim (2010) Kim, W.-T. 2010, ApJ, 725, 1069
  * Kramida et al. (2021) Kramida, A., Ralchenko, Y., Reader, J., & and NIST ASD Team. 2021, NIST Atomic Spectra Database (ver. 5.9), [Online, 2022, August 10]. National Institute of Standards and Technology, Gaithersburg, MD
  * Kruckow et al. (2018) Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908
  * Lau et al. (2022a) Lau, M. Y. M., Hirai, R., González-Bolívar, M., et al. 2022a, MNRAS, 512, 5462
  * Lau et al. (2022b) Lau, M. Y. M., Hirai, R., Price, D. J., & Mandel, I. 2022b, MNRAS, 516, 4669
  * Livio & Soker (1988) Livio, M. & Soker, N. 1988, ApJ, 329, 764
  * Mandel & Broekgaarden (2022) Mandel, I. & Broekgaarden, F. S. 2022, Living Reviews in Relativity, 25, 1
  * Marchant et al. (2021) Marchant, P., Pappas, K. M. W., Gallegos-Garcia, M., et al. 2021, A&A, 650, A107
  * Meyer & Meyer-Hofmeister (1979) Meyer, F. & Meyer-Hofmeister, E. 1979, A&A, 78, 167
  * Moreno et al. (2022) Moreno, M. M., Schneider, F. R. N., Röpke, F. K., et al. 2022, A&A, 667, A72
  * Nandez et al. (2014) Nandez, J. L. A., Ivanova, N., & Lombardi, J. C., J. 2014, ApJ, 786, 39
  * O’Connor et al. (2023) O’Connor, C. E., Bildsten, L., Cantiello, M., & Lai, D. 2023, ApJ, 950, 128
  * Ohlmann et al. (2016a) Ohlmann, S. T., Röpke, F. K., Pakmor, R., & Springel, V. 2016a, ApJ, 816, L9
  * Ohlmann et al. (2016b) Ohlmann, S. T., Röpke, F. K., Pakmor, R., Springel, V., & Müller, E. 2016b, MNRAS, 462, L121
  * Ondratschek et al. (2022) Ondratschek, P. A., Röpke, F. K., Schneider, F. R. N., et al. 2022, A&A, 660, L8
  * Ostriker (1999) Ostriker, E. C. 1999, ApJ, 513, 252
  * Paczyński (1976) Paczyński, B. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 75
  * Passy et al. (2012) Passy, J.-C., De Marco, O., Fryer, C. L., et al. 2012, ApJ, 744, 52
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Podsiadlowski (2001) Podsiadlowski, P. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 229, Evolution of Binary and Multiple Star Systems, ed. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, 239
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246
  * Prohaska et al. (2022) Prohaska, T., Irrgeher, J., Benefield, J., et al. 2022, Pure and Applied Chemistry, 94, 573
  * Prust & Chang (2019) Prust, L. J. & Chang, P. 2019, MNRAS, 486, 5809
  * Reichardt et al. (2020) Reichardt, T. A., De Marco, O., Iaconi, R., Chamandy, L., & Price, D. J. 2020, MNRAS, 494, 5333
  * Reichardt et al. (2019) Reichardt, T. A., De Marco, O., Iaconi, R., Tout, C. A., & Price, D. J. 2019, MNRAS, 484, 631
  * Reimers (1975) Reimers, D. 1975, Memoires of the Societe Royale des Sciences de Liege, 8, 369
  * Ricker & Taam (2008) Ricker, P. M. & Taam, R. E. 2008, ApJ, 672, L41
  * Ricker & Taam (2012) Ricker, P. M. & Taam, R. E. 2012, ApJ, 746, 74
  * Röpke & De Marco (2023) Röpke, F. K. & De Marco, O. 2023, Living Reviews in Computational Astrophysics, 9, 2
  * Sand et al. (2020) Sand, C., Ohlmann, S. T., Schneider, F. R. N., Pakmor, R., & Röpke, F. K. 2020, A&A, 644, A60
  * Shiber et al. (2019) Shiber, S., Iaconi, R., De Marco, O., & Soker, N. 2019, MNRAS, 488, 5615
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889
  * Stevenson et al. (2017) Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nature Communications, 8, 14906
  * Taam et al. (1978) Taam, R. E., Bodenheimer, P., & Ostriker, J. P. 1978, ApJ, 222, 269
  * Tauris & Dewi (2001) Tauris, T. M. & Dewi, J. D. M. 2001, A&A, 369, 170
  * Tauris & van den Heuvel (2006) Tauris, T. M. & van den Heuvel, E. P. J. 2006, in Compact stellar X-ray sources, Vol. 39, 623–665
  * Trani et al. (2022) Trani, A. A., Rieder, S., Tanikawa, A., et al. 2022, Phys. Rev. D, 106, 043014
  * Tutukov & Yungelson (1993) Tutukov, A. V. & Yungelson, L. R. 1993, MNRAS, 260, 675
  * Vigna-Gómez et al. (2018) Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MNRAS, 481, 4009
  * Voss & Tauris (2003) Voss, R. & Tauris, T. M. 2003, MNRAS, 342, 1169
  * Warner (1995) Warner, B. 1995, Cataclysmic variable stars, Vol. 28
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Whelan & Iben (1973) Whelan, J. & Iben, Icko, J. 1973, ApJ, 186, 1007
  * Yarza et al. (2022) Yarza, R., Everson, R. W., & Ramirez-Ruiz, E. 2022, arXiv e-prints, arXiv:2210.00010

##  Appendix A Resolution study

Figure 10: Similar to Fig. 2, but also showing the simulation at a lower
spatial resolution and at a higher time resolution.

To prove the robustness of the results we obtained with our 1D CE method
described in Sect. 2, we perform a small resolution study. We rerun the
simulations for $q=0.25$, presented in Sect. 3. For one simulation, we
increase the time resolution by a factor $2$, for a second simulation we
decrease the spatial resolution by a factor of $1.6$. The results of these
simulations are shown in Fig. 10.

The spiral in the behavior of the companion is mostly unchanged when compared
to the original simulation. There are however some differences in the envelope
ejection between the different simulations. The exact timing of the main
envelope-ejection event varies by about 100d. In addition, the shape of the
ejection curve also depends on the resolution. Once a more continuous ejection
is established after the main ejection events at $t=900-1000\,\mathrm{d}$, the
values of $f_{\mathrm{ej}}$ only vary by 5%, i.e., within our expected
tolerances.

The reason for the differences in the envelope ejection originates mostly from
differences in the surface properties between the simulations. As described in
Sect. 2.3, we perform a relaxation run before the CE simulations where we
switch the outer boundary condition and turn on the hydro mode. These changes
cause small amplitude oscillations close to the surface of the CE (see also
Fig. 4a). The oscillations are stochastic and therefore also some of the
surface properties, e.g., surface velocity. This means that the exact timing
of when the surface velocity exceeds the escape velocity has a minor
dependence on the initial oscillations. This causes the difference in the
starting time of the ejection. The total ejected mass towards the end of the
simulations, the quantity in which we are mostly interested varies within our
tolerances.

Based on this resolution study, we conclude that the final values for the
orbital separation do not depend on the resolution. The total ejecta mass has
a minor dependence on the resolution but varies within our tolerances of a few
percent.