# Common envelope binary interaction simulations between a thermally-pulsating
AGB star and a low mass companion

Miguel González-Bolívar ,1,2 Orsola De Marco,1,2 Mike Y. M. Lau,3,4 Ryosuke
Hirai3,4 and Daniel J. Price3
1 Department of Physics and Astronomy, Macquarie University, Sydney, NSW 2109,
Australia
2 Astronomy, Astrophysics and Astrophotonics Research Centre, Macquarie
University, Sydney, NSW 2109, Australia
3 School of Physics and Astronomy, Monash University, Clayton, Victoria 3800,
Australia
4 OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery,
Australia
E-mail<EMAIL_ADDRESS>0000-0002-5939-9269
0000-0002-1126-869X 0000-0002-6592-2036 0000-0002-8032-8174
0000-0002-4716-4235

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

At least one in five of all planetary nebulae are the product of a common
envelope (CE) interaction, where the companion in-spirals into the envelope of
an asymptotic giant branch (AGB) star ejecting the nebula and leaving behind a
compact binary. In this work we carry out 3D smoothed particle hydrodynamics
simulations of the CE interaction between a 1.7 M⊙ AGB star and a 0.6 M⊙
companion. We model the AGB structure using a 1D stellar model taken at the
seventh thermal pulse. The interaction takes place when the giant is on the
expanding phase of the seventh thermal pulse and has a radius of 250 R⊙. The
post-CE orbital separations varies between 20–31 R⊙, with the inclusion of
recombination energy resulting in wider separations. Based on the observed
short in-spiral timescales, we suggest that thermal pulses can trigger CEs,
extending the ability of AGB stars to capture companions into CEs, that would
lead to the prediction of a larger population of post-AGB post-CE binaries.
Simulations that include a tabulated equation of state unbind a great deal
more gas, likely unbinding the entire envelope on short timescales. The shape
of the CE after the in-spiral is more spherical for AGB than RGB stars, and
even more so if recombination energy is included. We expect that the planetary
nebula formed from this CE will have different features from those observed by
Zou et al. 2020.

###### keywords:

hydrodynamics — methods: numerical — Stars: AGB — binaries: close

††pubyear: 2022††pagerange: Common envelope binary interaction simulations
between a thermally-pulsating AGB star and a low mass companion–18

##  1 Introduction

Common envelope (CE) interactions occur in binary stellar systems when one of
the stars, whether because of radial expansion or changes in the binary orbit,
engulfs the companion. This interaction is a critical step in the formation of
compact binaries, such as cataclysmic variables, the progenitors of type Ia
supernovae, at least one in five of all planetary nebulae (PN) and even of
binaries that later become gravitational waves sources. Since such
interactions have timescales of the order of months to years, the number of
systems observed during a CE is low (Ivanova et al., 2013; De Marco & Izzard,
2017).

At least one in five PN have a post-CE central star binary (e.g., Miszalski et
al., 2009; Jacoby et al., 2021). These close binary stars (typical periods
$\lesssim$3 days) owe their short periods to a CE in-spiral following unstable
mass transfer between an asymptotic giant branch (AGB) star and a companion
(main sequence or white dwarf). The PN around these binaries tend to have
axial symmetry or no apparent symmetry. PN in general tend to have short
lifetimes in virtue of their expansion, which leads to visibility times of
$\sim 50\,000$ years (Schönberner, 1983). When the central star of a PN is a
post-CE binary, this ensures that the PN is the ejected CE111Some of the PN
material could technically have been ejected before the CE interaction. This
would be the case in those instances when the AGB star has entered the short
superwind phase just ahead of the CE in-spiral, something that would require
some degree of fine tuning, but is not impossible.. This in turn means that
the observed nebular morphology and dynamics, as well as the central stars’
and binary parameters may provide observational constraints on CE interaction
energies and timescales. Additionally, the post-CE central stars of PN provide
an observational data point for the end state of a CE interaction.

Most CE simulations in the literature have used red giant branch (RGB),
instead of AGB donors, because they exhibit greater 3D hydrostatic stability
in the numerical domain. All giant stars are in an evolutionary phase with
relatively short timescales. RGB and AGB stars expand on timescales of
millions of years. Thermally-pulsating AGB stars, the last phase of AGB
evolution, expand significantly on timescales of decades-to-centuries. These
thermodynamically unstable objects tend to expand or contract on dynamical
timescales once they are mapped on the 3D computational domain, because of
typical numerical shortcuts, such as a lack of adequate resolution, or the use
of a simplistic equation of state. The problem is worse at the end of the AGB
phase, where even 1D stellar structure and evolution codes find it
increasingly difficult to find a hydrostatic solution for the envelope
(Ohlmann et al., 2017).

There are only a handful of AGB primary CE simulations in the literature.
Sandquist et al. (1998) and Staff et al. (2016) carried out simulations with
relatively massive AGB primaries (3-5 M⊙). More recently, Chamandy et al.
(2020) and Sand et al. (2020) have implemented CE simulations using low mass
AGB stars as donors ($1.78$ and $0.97$ M⊙, respectively, evolved from 2.0 and
1.2 M⊙ main sequence stars, respectively). These simulations were carried out
using the grid code astrobear (Carroll-Nellenback et al., 2013) and moving
mesh code arepo (Springel, 2010). Donors were in the non-thermally-pulsating
AGB phase (the stars had radii of 122 and 173 M⊙, respectively), when the star
is more dynamically stable as opposed to the later stage. Particularly, there
is no previous work done using a thermally pulsating AGB star as the CE donor.

In order to pursue the link between PN and CE interactions, García-Segura et
al. (2018, 2020, 2021), Frank et al. (2018) and Zou et al. (2020) carried out
hydrodynamics simulations of the “fast wind” (the wind that the post-AGB star
emits and that ploughs into the previously emitted, slow AGB wind) ploughing
into gas ejected during a CE. The first group used a simulation by Ricker &
Taam (2012), which, by virtue of its small computational domain, was halted
after $56$ days. The PN simulation was in 2.5D ensuring that symmetry would be
found. The latter two papers used instead a simulation by Reichardt et al.
(2019), which, having been carried out using a smoothed particle hydrodynamics
(SPH) technique, could be followed for a time of 14.5 years after the start,
resulting in more expanded ejecta. Their PN simulation was in 3D but covered a
much shorter PN evolution time than the one by García-Segura et al. (2018).
These initial efforts, however, showed promise to connect PN morphologies with
CE events.

Another possible shortcoming of the PN formation work of García-Segura et al.
(2018) and subsequent papers as well as Zou et al. (2020) is that the CE
simulations used to provide initial conditions for the simulation had not
fully ejected their envelopes. Hence the velocity of the gas into which the
super-wind ploughs was slower than is expected in a real PN formation
situation, where the AGB has ejected its envelope. Today, we have simulations
that can fully eject the envelope, by adopting a tabulated equation of state
including recombination energy that can be released into the envelope. While
these simulations may overestimate the dynamical effect of the recombination
energy by assuming an adiabatic gas, they allow one to test its effect on the
formation of the PN.

Finally Iaconi & De Marco (2019) found that the observed post-RGB CE systems
appear to be systematically more efficient in ejecting the envelope than the
post-AGB objects, or, in other words, they have similar final separations
despite having considerably more bound envelopes. At the same time they found
that the lack of simulations with low mass AGB stars made it impossible to
determine whether simulations retrieve this characteristic.

In this paper we therefore carry out the first 3D hydrodynamical simulation of
a CE interaction with a thermally-pulsating AGB star, carried out with an
ideal gas equation of state as well as with a tabulated equation of state that
includes the effects of radiation pressure and recombination. The initial
mass, 2.0 M⊙, was chosen because stars with masses lower value than this tend
to grow as large on the RGB as on the AGB, with the implication that most CE
interactions would take place on the RGB.

The structure of the paper is as follows: we outline the evolution through the
thermally pulsating AGB phase (TP-AGB) of the donor star in Section 2. In
Section 3, we discuss the simulation setup: the star relaxation procedure is
detailed in Section 3.1, the binding energy of the AGB stellar structures is
discussed in Section 3.2 and the binary setup is discussed in Section 3.3. The
results are presented in Section 4, starting with the orbital evolution in
Section 4.1, an analysis of the unbound mass in Section 4.2, a discussion on
the envelope rotation in Section 4.3, an appraisal of the morphology of the
common envelope in Section 4.4 and an projection of the likely envelope dust
content in Section 4.5. We finally present our discussion in Section 5 and
summarise in Section 6.

##  2 The thermally-pulsating AGB primary

Figure 1: Evolution during the TP-AGB phase. Top left panel: evolution of the
radius and luminosity. Top right panel: an enlargement of a thermal pulse; the
orange vertical lines indicate the star contraction period right after the
helium shell flash. Bottom left panel: evolution of the thermal and dynamical
timescales compared to luminosity outputs. The red symbol represents the mesa
stellar profile used for the common envelope simulation at $R=260$ R⊙. Bottom
right panel: Kippenhahn diagram of the star during the thermal pulse. The
contour represents the helium abundance clearly showing the boundary of the
hydrogen and helium shells. The blue dotted line at $\sim 0.562$ M⊙ is the
helium core limit. The square magenta pattern indicates the convective
regions. The time $t=0$ is set, arbitrarily, two pulses before the pulse we
analyse.

The AGB is one of the last evolutionary phases experienced by low-and
intermediate-mass stars ($\lesssim$ 8 M⊙). It is divided into an early and
late/thermally-pulsating phase. The early AGB phase is characterised by thick
helium shell burning above a mostly carbon-oxygen core and a slow and steady
radial expansion. Above this layer, there is a helium-rich shell, and above
that, a hydrogen-burning shell (Lattanzio & Karakas, 2001). After the helium
shell thins out, the late-AGB phase starts, which is underpinned by a series
of thermal pulses as the star alternates between explosive, runaway shell
helium burning and thermally stable shell hydrogen burning.

Figure 1 illustrates the relevant aspects of stellar evolution through this
phase, where we have used the Modules for Experiment in Stellar Astrophysics
(mesa, version 12778; Paxton et al., 2011; Paxton et al., 2013, 2015) to
calculate the structure of an AGB star with a zero age main sequence mass of 2
M⊙. Here we see a cyclical ignition of the helium shell with concomitant
variability of the hydrogen shell burning rate (green and cyan solid lines in
Figure 1, top left-panel, respectively). Helium shell flashes repeat
approximately every 100 000 years.

Thermal flashes are caused by the thermonuclear ignition of a thin layer of
helium, at the base of the helium shell just above the CO degenerate core and
below a large hydrogen-rich envelope. The instability that leads to the
flashes is explained by a combination of the thin-shell instability and
partial degeneracy (Kippenhahn & Weigert, 1990). Flashes follow long,
quiescent periods of hydrogen burning just above the helium shell. During
these periods, hydrogen burning shell accretes onto the helium shell until the
base of the shell ignites under partial degeneracy conditions.

These ignitions of the helium layer produce a thermal energy increase of the
inner layers, which leads to its expansion and to a decreased burning rate. At
this point there is relatively low energy input in any part of the star and
the surface of the star begins to contract. Eventually, the energy radiated
from the helium shell produces a radial expansion of the envelope, which
reaches a larger radius compared to the hydrogen shell burning phase before
the flash.

For stars with main sequence masses of 2 M⊙, both hydrogen and helium burning
layers are at $R<0.1$ R⊙ by the time the TP-AGB starts. We define those radii
as the (radial) cell with the highest burning rate of hydrogen or helium.

At the time of the flash, the stellar radius is 234 R⊙. After a period of
contraction, the radius returns to its initial value (900 years after the
flash peak), expands to 304 R⊙ (1050 years after the flash peak) and
eventually returns to its initial value 37 000 years after the peak of the
flash. The interval between the helium flash and the maximum stellar radial
expansion is $\sim 400$ years. This lag is of the order of the thermal
timescale of the layers above the helium shell (bottom plot, black dashed line
in Figure 1, calculated from the 1D mesa model as the time for the layers
above the helium burning shell to radiate all their thermal energy), and this
time is about two orders of magnitude slower than the dynamical timescale of
these regions.

The upper layers ($R\geq 0.2$ R⊙) are convective during the TP-AGB phase.
Figure 1, bottom-right panel, shows the Kippenhahn diagram during the same
time span as the previous plot. The orange vertical lines are the same as for
the top-right panel and are plotted here for reference. The non-burning
hydrogen envelope is convective, while convection extends all the way to the
core just preceding the maximum helium shell luminosity output, during the
early part of the thermonuclear runaway of the helium shell. This is the
convective event responsible for the third dredge-up. This combination of
events creates a unique situation for the star. The thermally-pulsating star
is a mostly convective, rapidly expanding star. Thermal pulses expand the star
into volumes never reached before creating new opportunity for interaction.
Accounting for tidal interactions between the star and its companion shows
that the thermal pulses are able to progressively draw in the companion such
that even if one thermal pulse expansion does not result in a CE, the
following pulse can generate a CE even before the star expands past the
maximum radius of the previous pulse (Madappatt et al., 2016, figure 5). The
next question is whether the short duration of the pulse radial expansion is
such that a CE would be stymied by the contraction of the star post-pulse.

##  3 Simulation setup

Figure 2: The mesa and softened density profiles of the four models used. The
black dashed vertical line is at $R=2$ R⊙, which is the value set for $r_{\rm
soft}$ in all the models except for R-hi. The density and pressure profiles
outside the softening length are identical for all three softened models.
I-hi, R-hi and M-hi are the stellar profile mapped in the 3D computational
domain using the ideal gas, ideal gas plus radiation and tabulated equations
of state, respectively.

In this section we describe the initial conditions for our 3D CE simulation,
starting with the single star 3D mapping method, the relaxation procedure to
obtain a stable model of the TP-AGB star and the subsequent binary star setup.
We use the mesa model of the TP-AGB star discussed in Section 2. The zero age
main sequence mass is 2 M⊙ with solar metallicity ($Y=0.28$; $Z=0.02$). The
stellar wind parameter implemented had the code default values: cool wind RGB
scheme was set to “Reimers”, with a Reimers scaling factor of $0.1$; the cool
wind AGB scheme was “Blocker” with a scaling factor of 0.5 and the RGB to AGB
wind switch was set at $10^{-4}$. We let the model evolve until it reaches the
seventh thermal pulse expansion. At this point, the stellar radius is 260 R⊙,
the carbon-oxygen core mass is 0.45 M⊙ and the total stellar mass is 1.7 M⊙.
The peak radius of 304 R⊙ would be reached $150$ years later (Figure 1). We
took the 1D stellar profile and mapped it into the 3D SPH code phantom (Price
et al., 2018) using the procedure explained in Section 3.1. We then proceeded
to create four different, stable 3D setups of the same mesa model: two with an
ideal gas equation of state at low and high resolutions of $1.37\times 10^{5}$
and $1.37\times 10^{6}$ SPH particles, respectively. The other two models used
an equation of state including ideal gas plus radiation pressure and a
tabulated equation of state implemented by Reichardt et al. (2020) and taken
from mesa, both with the same number of particles as the high resolution
simulation for the ideal gas case. Both equations of state assume local
thermodynamic equilibrium, which is a good approximation in the optically
thick regions of the star.

###  3.1 Core excision and star relaxation procedure

Model | Equation | $N_{\rm part}$ | Relaxation | $r_{1,\rm soft}$ | $r_{2,\rm soft}$ | $M_{\rm core}$ | $t_{\rm model}$ | $E_{\rm bin}$ | Artificial
---|---|---|---|---|---|---|---|---|---
| of state |  | technique | (R⊙) | (R⊙) | (M⊙) | ($t_{\rm dyn}$) | ($10^{46}$ erg) | conductivity
I-lo | ideal gas | $1.37\times 10^{5}$ | damp | 2 | 2 | 0.5633 | 5.2 | $3.12$ | $1$
I-hi | ideal gas | $1.37\times 10^{6}$ | damp | 2 | 2 | 0.5634 | 33.9 | $3.11$ | $1$
R-hi | ideal + rad | $1.37\times 10^{6}$ | relax-o-matic | 3 | 3 | 0.5642 | 17.3 | $2.92$ | $0.8$
M-hi | tabulated | $1.37\times 10^{6}$ | relax-o-matic | 2.5 | 2.5 | 0.5639 | 5.2 | $-0.11$ | $0.1$
† For our TP-AGB star, $t_{\rm dyn}\approx 42$ days |  |  |  |
Table 1: Parameters of the stellar models used as starting conditions for our
common envelope simulations. The eighth column, $t_{\rm model}$ is the number
of dynamical times that the stellar structure was evolved in isolation, before
implementing it in the binary setup. These times correspond to the colored
squares in Figure 3.

A procedure to generate a reasonable stellar structure that is in equilibrium
over sufficiently long timescales is an integral and important component of
all hydrodynamics simulations of stars. Typically stars generated by 1D
implicit codes have high resolution and sophisticated equations of state and
are not stable once they are mapped in the 3D computational domain at lower
resolution, a different equation of state, grid mismatch and less physics
included, such as, for example, the lack of radiation transport.

Most approaches followed today are similar to those in the early simulations
of Rasio & Livio (1996, who used an SPH approach) and Sandquist et al. (1998,
who used a grid-based approach), who replaced the cores of their giant stars
with a point mass particle, because the extremely high central densities are
orders of magnitude larger than near the surface, resulting in the
Courant–Friedrichs–Lewy condition near the centre imposing impractically short
time steps. In order to include a central point mass particles, these authors
implemented a smoothed gravitational potential over a certain distance from
the center. This technique solves the time stepping problem in the core region
and similar solutions have been adopted by stellar numerical simulations ever
since. They also added a velocity-dependent friction term (e.g. Gingold &
Monaghan, 1977) for the gas applied over a certain time, that aims to dampen
abnormal speeds due to the mapping process or lack of hydrostatic equilibrium.
Running the code with this damping term for a few dynamical timescales relaxes
stellar profiles mapped from hydrostatic equilibrium solutions.

Variations on this technique were used before in CE simulations with the
phantom code (Iaconi et al., 2018; Reichardt et al., 2019, 2020), when the
primary was an RGB star. In this work we have used two different techniques to
excise the core and relax the stellar profile. The first is a slightly
modified version of the technique presented by Reichardt et al. (2019), which
was very similar to what was devised by Ohlmann et al. (2017). The second
procedure was introduced by Lau et al. (2022b) in order to accurately
reproduce the density profile in their work with a more massive star and which
resulted to be superior in the stability of our TP-AGB star with an equation
of state other than ideal gas. Below we give an account of these procedures.
All relevant parameters are listed in Table 1. The first prescription was used
to create an input stellar model that adopted an ideal gas equation of state
(in Table 1 we nickname this first procedure damp). We start with the
selection of a softening radius, $r_{\rm 1,soft}$, for the potential of the
point mass particle that replaces the core. The potential is softened within a
volume that has twice this radius. The mesa density profile is modified inside
$2\times r_{\rm 1,soft}$, such that the density is lower and almost constant
but it joins smoothly with the original mesa profile at $2\times r_{\rm
1,soft}$.

The value of $r_{\rm 1,soft}$ needs to be smaller than the final core-
companion orbit, but large enough that the time step in the simulation is not
impractically short. Once a value of $r_{\rm 1,soft}$ is selected, phantom
then calculates the mass for the core point mass particle ($M_{\rm core}$ in
Table 1). After the softened profile is constructed, a geometric distribution
of particles, with an approximate density similar to that of the softened
profile, is generated (see Figure 2, where the three profiles mapped using
three equations of state are displayed alongside the original 1D mesa
profile). After finalising the mapping process into the phantom computational
domain, we evolve the 3D star with the velocity-dependent damping factor
implemented by Reichardt et al. (2019) for 5 dynamical times, $5t_{\rm dyn}$.
We calculate $t_{\rm dyn}$ as:

$t_{\rm dyn}=\frac{R}{v_{esc}}=\sqrt{\frac{R^{3}}{2GM}},$ (1)

where G is the gravitational constant and R and M are the photosphere radius
and mass of the original mesa stellar profile, respectively. We then let the
simulation evolve for an additional 5 dynamical times with no velocity damping
to assess the numerical stability of the model. In Appendix A we give further
details of this process and show the particle distributions and density cross
sections for the stellar models taken at all stages in the relaxation process
just described.

Figure 3: Radial evolution of the implemented models during and after
relaxation. The “particle radius” is shown in the left panel and the “density
radius” evolution is showing in the right panel. The square symbol represents
the profile we have used in the binary simulations. The black dashed line is
the radial evolution of the 1D, mesa profiles during the thermal pulse.

The second relaxation technique used in this paper (called relax-o-matic in
Table 1) similarly follows from mapping the 1D profile and creating a density
distribution inside the core region based on the selected softening length.
This time, however, we set the density and pressure within
$r<r_{\mathrm{1,soft}}$, such that the entropy is constant. In the case of an
ideal gas the specific entropy in the central region is given by:

$s=\frac{k_{B}}{m_{h}}\frac{1}{\mu}\ln\left(\frac{T_{\rm
soft}^{3/2}}{\rho_{\rm soft}}\right),$ (2)

whereas for ideal gas with radiation pressure and for the tabulated equation
of state we use:

$s=\frac{k_{B}}{m_{h}}\frac{1}{\mu}\ln\left(\frac{T_{\rm
soft}^{3/2}}{\rho_{\rm soft}}\right)+\frac{4a_{\rm rad}}{3}\frac{T_{\rm
soft}^{3}}{\rho_{\rm soft}},$ (3)

where $k_{B}$ is the Boltzmann constant, $a_{\rm rad}$ is the radiation
constant, $m_{h}$ the hydrogen atomic mass, and $T_{\rm soft}$ and $\rho_{\rm
soft}$ are the temperature and density values at $r=2\times r_{\rm soft}$,
respectively (we only used the constant entropy prescription for R-hi and
M-hi; see Hirai et al. 2020, appendix B for details of this procedure). At
this point, a random distribution of SPH particles is generated. Particles are
generated in pairs such that their positions are symmetric with respect to the
origin of the domain. This initially removes any net gravitational force on
the core. Instead of using velocity damping, we apply an artificial particle
shifting method where the thermal energy of the particles is updated after
each iteration so that it follows the given $P(m)/\rho(m)$ profile (for a more
detailed explanation see appendix C of Lau et al. 2022b). The relaxation
procedure iterates until the particle distribution resembles the original
softened density profile within a certain tolerance, the total kinetic energy
falls below $10^{-7}$ of the potential energy, or after 1000 iterations - for
our case, after 1000 iterations $E_{\rm kin}/E_{\rm pot}\sim 2\times 10^{-7}$.
Once again, after the relaxation procedure was complete, we evolved the star
in isolation for a few dynamical times to test the model’s hydrostatic
stability. We used the relax-o-matic technique to relax the stellar structures
with radiation pressure and a tabulated equation of state. We adopted default
artificial conductivity parameter values (see section 2.2.8 of Price et al.,
2018) of unity for the ideal gas models (I-hi and I-lo), 0.8 for R-hi and 0.1
for M-hi (see Table 1). These values are the same for the relaxation phase and
for the CE simulations. The default value is set to ensure accurate treatment
of contact discontinuities (Price, 2008). The conductivity is second order in
phantom, with the effective thermal conduction $\kappa_{\rm
cond}\propto\alpha_{u}h^{2}(\nabla\cdot{\bf v})$, where $h$ is the resolution
length (c.f. Price et al., 2018). It thus vanishes when velocity gradients are
zero and the resolution is high, but we found that the remaining heat transfer
could nevertheless lead to a slight undesirable expansion of the star over the
timescales of interest. The simple solution was to lower the pre-factor for
the R-hi and M-hi models to $\alpha_{u}=0.8$ and $0.1$, respectively.

The fidelity of the stellar structures was further assessed by inspecting the
density distributions as well as by assessing the individual velocities of the
SPH particles. This analysis leads to the conclusion that the stellar
structures are sufficiently in hydrostatic equilibrium to proceed with the
binary interaction simulations. This numerical stability analysis is fully
presented in Appendix A.

In Figure 3, we plot the stellar radius over time. Time zero is the time when
we start running the simulation of the single star to check its degree of
stability. The square symbols indicate the times used to set up the binary
simulations. As can be seen in Figure 3, we run the single star simulations
for longer to determine the stability of the star over a time that is almost
as long as the entire CE simulation. In Figure 3, “particle radius” is the
average distance of the 0.5 per cent of all particles that are farthest from
the center of mass, whereas the “density radius” is defined as the volume-
equivalent radius (Nandez et al., 2014) of all the SPH particles whose density
is higher than the least dense SPH particle at $t=0$ ($2.29\times 10^{-9}$,
$1.81\times 10^{-9}$, $3.16\times 10^{-9}$, and $3.16\times 10^{-9}$ g cm-3
for the I-lo, I-hi, R-hi and M-hi, respectively).

For the I-lo model, the stellar “particle radius” is still growing by 14
years, while the “density radius” seems to be asymptotically reaching a value
of $\sim$280 R⊙. On the other hand, I-hi remains reasonably constant in both
radius definitions for the entire time. The R-hi model’s radius is gently
growing over the 5 years for which the single model was run, while the M-hi
model rapidly expands but then just as rapidly decreases and stabilises
thereafter at a value of $\sim$285 R⊙, depending on the radius definition.

In Figure 3 we also plot the natural radial growth of the 1D mesa model which,
being at the top of a thermal pulse is still growing on relatively short
timescales. The radial growth of our single star models is due to residual
motion that we have not been able to eliminate. Later on we discuss the slow
expansion of our numerical stellar models in relation to the expected
expansion of the actual TP-AGB star during the same time as the CE
interaction.

###  3.2 The binding energy of the AGB star

To calculate the binding energy of our models we perform a numerical
integration using two terms: the gravitational potential energy, which is the
same for all the models, and the internal energy, which depends on the
equation of state implemented. This second term represents one of the sources
of energy that will be used to accelerate the gas during the common envelope
simulation.

In Figure 4, we plot the the cumulative binding energy as a function of radius
$r$, defined as:

$E_{\rm
bind}(r)=\int_{M}^{m(r)}\left[u(r^{\prime})-\frac{Gm^{\prime}}{r^{\prime}}\right]dm^{\prime}$
(4)

where M is the total stellar mass, $G$ is the gravitational constant and
$m^{\prime}$ is the mass coordinate. The integration is carried out from the
surface inward. The internal energy, $u(r)$, depends on the equation of state.
As expected, when we only account for the gravitational energy, the
integration has the highest binding energy. Adding the internal energy
corresponding to an ideal gas with the density of the stellar profile (green
and orange curves in Figure 4), the full integration gives a total binding
energy equivalent to half of the gravitational binding energy, satisfying the
Virial theorem. This result also applies when adding the radiation pressure,
since this does not contribute significantly for low mass stars. Lastly,
calculating $u(r)$ with the tabulated equation of state, the contribution of
the latent energy of ionised hydrogen and helium decreases the binding energy
in the inner layer to such level that the central region ($<75$ R⊙) has a
negative integrated binding energy value.

Figure 4: The cumulative binding energy for the TP-AGB star. The blue curve includes only the gravitational potential energy (second term in the integrand in Equation 4). Orange, green and red lines include the $u(r)$ term in Equation 4 which is the internal energy for an ideal gas (equivalent to the thermal energy), ideal gas with radiation pressure (equivalent to thermal energy plus radiation energy) and the mesa equations of state (equivalent to thermal, radiation and recombination energies), respectively. All of the quantities correspond to the softened profile generated during the core excision of the original mesaprofile. Model |  $E_{\rm bin}$ (phantom) | Equivalent radius | Equivalent mass
---|---|---|---
| ($10^{46}$ erg) | (R⊙) | (M⊙)
I-lo | 3.118 | 0.084 | 0.5603
I-hi | 3.105 | 0.084 | 0.5603
R-hi | 2.918 | 0.086 | 0.5604
M-hi | -0.1062 | 2.004 | 0.5645
Table 2: Envelope binding energy of the implemented models. The second column
is the binding energy calculated using Eq. 4 and the phantom profiles. The
third column is the radius coordinate in the mesa profile that has the same
binding energy value. The fourth column is the mass inside that radial
location. The binding energy, including only gravitational potential energy of
the model, is $6.036\times 10^{46}$ erg.

This analysis works as a sanity check: in Table 2 we see that the mass inside
the location with the binding energy calculated using the phantom file
(“Equivalent mass”, in Table 2), is less than 1 per cent different from the
value of the degenerate core mass in the mesa profile ($M_{\rm core}$ in Table
1) or the value of the point mass we used.

###  3.3 Binary setup

Each of our four stellar models, the two models with ideal gas equation of
state at low and high resolutions (I-lo and I-hi), the high resolution model
with radiation pressure added (R-hi) and the model with a tabulated equation
of state (M-hi) are combined with a companion with a mass of 0.6 M⊙ at a
separation of 550 R⊙. For simplicity, we set the softening length of the
companion, $r_{\rm 2,soft}$, to be the same as that of the primary’s core. The
initial separation is such that the donor with a radius of 260 R⊙ fills its
Roche lobe, which was calculated using the formula of Eggleton (1983). The
primary and companion are set to orbit one another with zero eccentricity and
Keplerian velocities of $35.3$ km s-1. The initial period is 2.7 years. None
of our simulations were relaxed in the potential of the binary. For a
discussion of the differences of such a setup see Reichardt et al. (2019).

The primary is not rotating, and we do not relax the primary in the potential
of the companion. When the companion is introduced in the computational domain
the primary suffers a “plucking” effect, where a sudden tide is induced that
sets off oscillations. This results in a slight additional expansion of the
star.

Table 3 shows the binary setup parameters for the four models and for the
simulation conducted by Reichardt et al. (2019, hereafter Rei19) that we use
as reference.

Model | $q$ | $a_{\rm i}$ | $a_{\rm max,pl}$ | $a_{\rm f,pl}$ | $a_{\rm f}$ | $a_{\rm f,fit}$ | $t_{\rm max,pl}$ | $t_{\rm f,pl}$ | $t_{\rm f}$ | $e_{\rm f,pl}$ | $M_{\rm b,PN}$ | $M_{\rm b,PN}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  | (R⊙) | (R⊙) | (R⊙) | (R⊙) | (R⊙) | (yr) | (yr) | (yr) |  | (M⊙) | (%)
I-lo | $0.35$ | $550$ | $166$ | $54$ | $25.9$ |  $25.0$ – $27.4$ | $13.3$ | $14.2$ | $17.7$ | $0.06$ | $0.97$ | 85
I-hi | $0.35$ | $550$ | $161$ | $32$ | $20.6$ |  $19.1$ – $19.4$ | $15.3$ | $16.8$ | $27.4$ | $0.03$ | $0.96$ | 84
R-hi | $0.35$ | $550$ | $162$ | $44$ | $25.8$ |  $23.3$ – $24.8$ | $7.5$ | $8.7$ | $15.2$ | $0.09$ | $0.97$ | 85
M-hi | $0.35$ | $550$ | $250$ | $35$ | $32.7$ |  $31.8$ – $32.2$ | $4.6$ | $8.7$ | $9.8$ | $0.005$ | $0.69$ | 60
Rei19 | $0.68$ | $218$ | $61$ | $23$ | $18.5$ |  $17.9$ – $18.6$ | $13.3$ | $14.2$ | $20.1$ | $0.01$ | $0.46$ | 94
Table 3: Simulation parameters for our simulations and for the simulation of
Reichardt et al. (2019). Columns are as follows: $q$ is the mass ratio
$M_{2}/M_{1}$; $a_{1}$, $a_{\rm max,pl}$, $a_{f,pl}$ and $a_{f}$ are the
separations at the start of the simulation, at the time of fastest in-spiral,
at the end of the in-spiral, as defined in the text, and at the end of the
simulation. Column $a_{\rm f,fit}$ contains the range of the extrapolated
final separations (periastron and apastron separations). Columns $t_{\rm
max,pl}$, $t_{f,pl}$ and $t_{f}$, are the corresponding times. The last two
columns list the bound mass at the time of the morphology comparison described
in Subsection 4.4.

##  4 Simulation results

Figure 5: Density cross sections in the orbital plane (first, third and fifth
columns) and the perpendicular plane (second, fourth and last columns) of all
models at the start (first two columns) in the middle (third and fourth
columns) and at the end (last two columns - note these slices are zoomed out)
of the plunge-in phase. Columns 1 to 4 have panels of $10^{4}$ R⊙on each side.
Columns 5 and 6 have panels of $1.5\times 10^{4}$ R⊙on each side. The rows,
from top to bottom, represent the four models: I-lo, I-hi, R-hi and M-hi,
respectively. A complete movie of the I-hi, R-hi and M-hi models is available
by following the link: https://miguelglezb.github.io/mgb/simulations/2msun-tp-
agb.html.

Here we discuss the simulations’ results using the high resolution, ideal gas
equation of state simulation (I-hi) as our baseline. In Section 4.1 we analyse
the orbital evolution of all our simulations, including that of Reichardt et
al. (2019). In Section 4.2 we discuss the mass unbound during the simulations.
We analyse the rotational properties of the stars in Section 4.3. Finally, the
morphology of the post-CE system is explained in Section 4.4. All our
simulations conserve total energy and total angular momentum to 0.1 per cent
and 0.05 per cent, at worst, respectively (more information is presented in
Appendix B).

###  4.1 Orbital evolution

Figure 6: Left panel: orbital evolution of the common envelope models. The
black circles indicate the end of the plunge-in. The magenta square symbols on
the I-hi and M-hi orbital evolution curves correspond to the times of the
slices in Figure 11. Right panel: orbital evolution of the models, including
that of Reichardt et al. (2019). The x-axis have been shifted so that $t=0$
coincides with the steepest plunge-in moment, $t_{\rm max,pl}$ (Table 3), in
every model. The orbital separation is normalised to the initial separation
between point mass particles. The square symbols show the times that are
compared in Section 4.4.

Before we can describe the orbital evolution and compare simulations we must
define a quantitative criterion to determine the time of the end of the
dynamical part of the in-spiral. We define the dynamical in-spiral phase as
the time at which $|\dot{a}P|/a>0.01$, where $a$ is the orbital separation and
$P$ is the orbital period. This criterion was used by Ivanova & Nandez (2016)
as well as by Sand et al. (2020); the former used a value of $0.1$ while the
latter applied $0.01$. Because of the eccentricity that develops towards the
end of the in-spiral, we apply a Savitzky–Golay filter to the point mass
particles’ separation. The filter applies a piecewise polynomial with a user-
defined order, to a number of the data points within a window which is also
determined by the user. We set a first order polynomial and a window length
such that the residuals between the smoothed and the original data was at most
the sum of the point mass particles’ softening radii. The end of the plunge-in
separation for all simulations, determined in this way, along with the
separation at the steepest plunge-in point are listed in Table 3, alongside
the times when they happen. In Figure 5 we show density slices in the orbital
and perpendicular planes for all our models taken at three times: (i) an
arbitrary time shortly after the beginning of the simulation, when the
dynamical in-spiral has not yet started, (ii) at a time close to when,
visually, the orbital separation starts to decline more rapidly, and (iii) at
the end of the in-spiral as quantified by the criterion above (called $t_{\rm
pl,f}$ in Table 3). For simulation M-hi there is no clear delineation between
pre-plunge-in and plunge in, so for time (ii) we selected a somewhat arbitrary
moment in the early in-spiral.

What is immediately evident in Figure 5 is that the envelope expansion and
mass ejection are preferentially along the orbital plane. This ejection is
primarily through the $L_{\rm 2}$ point. By the end of the plunge-in, a spiral
structure is observed around the close binary, including a ring of material
ejected at time $16.8$ yrs in the case of I-hi. Striking is that with the
tabulated equation of state the distribution is more spherical and smooth. It
is possible that as the envelope expands and recombines the recombination
energy is delivered in a more spherically symmetric way and the mass is
accelerated accordingly.

At the beginning of each simulation, the primary star fills its Roche lobe
such that mass transfer starts immediately. The orbital period of the binary
at the start of the simulations is $2.7$ years. Although there is some
variability between models, the orbital separation decreases relatively slowly
for the first few years, after which the plunge-in phase takes over. This is
followed by a relatively sudden slowing down of the in-spiral for all models.
A plot of the separation between core and companion as a function of time is
shown in Figure 6.

Our baseline model, I-hi, spends $\sim$5 years in a pre-in-spiral phase.
Eventually gas starts to outflow from the $L_{2}$ and $L_{3}$ Lagrangian
points, at which point the in-spiral commences. Envelope gas expands at this
point, some of which becomes unbound.For I-hi, the pre-plunge-in phase lasts
visibly the longest of all models. For this model, the time of steepest
plunge-in is at $t=14.2$ years from the start of the simulation. After the
plunge-in, at about 17 yrs, the point mass particles move in an orbit with a
mild eccentricity, $e\approx 0.02$, and slowly decreasing average separation.
The separation at the end of the plunge-in as defined above (16.8 yrs) is 32
R⊙, while at the end of the simulation (27.4 yrs), the separation is 20 R⊙.
This is evident in Figure 6, right panel, where the sink separation curves are
normalised by their initial values, and are shifted so that $t=0$ is at the
time of steepest plunge-in, $t_{\rm max,pl}$ (Table 3), for every model. Here,
the slopes of the fast in-spiral phases are more easily compared, as is the
end of the in-spiral (Fig. 6, right panel).

At low resolution the time before the plunge-in is shorter. This was also
observed by Reichardt et al. (2019) and it is due to a higher mass transfer
rate through $L_{1}$ that hastens the plunge-in. It is possible that at even
higher resolution the pre-plunge-in time would be even longer as the early
mass transfer rate would be lowered further. Under steady conditions (a star
that is not naturally expanding or contracting over short timescales) it is
possible that the pre in-spiral time would lengthen considerably. The TP-AGB
star that we chose to simulate is, however, expanding over short timescales,
because of its thermal pulse, as described in Section 2 and as can clearly be
seen in Figure 3. In this situation it is likely that the early mass transfer
at the time of first Roche lobe overflow, may be hastened by the natural
expansion of the star. The plunge-in slope and the post-plunge-in separation
are very consistent between I-lo and I-hi (Figure 6, right panel), showing
that resolution does not play a large role in these quantities.

The values of the separations for I-lo and I-hi at the end of plunge-in,
$a_{\rm f,pl}$ in Table 3, are 54 and 32 R⊙, respectively, a difference of 51
per cent. Aligning I-hi and I-lo at the steepest in-spiral point (Figure 6,
right panel) and reading the separation at 4.4 years after that point, close
to the end of the simulations, the two values are very similar: 24.2 R⊙ and
25.8 R⊙, for I-hi and I-low, respectively. The eccentricity of the orbit after
the plunge-in is comparatively larger for the low resolution simulation, I-lo
(0.04) than for the high resolution one, I-hi (0.02).

As we have discussed in Section 3.1, the R-hi initial stellar model is more
extended and has higher overall particle speeds at its surface than the model
I-hi (Figures 3 and 14). This may be responsible for a shorter Roche lobe
overflow phase. The separation after the plunge-in, $a_{\rm pl,f}=44$ R⊙
(Table 3), can be compared to the equivalent value for I-hi of 32 R⊙. The
orbital separation 4.4 yrs after the time of maximum plunge-in is 25.3 R⊙ for
R-hi and 25.8 R⊙ for I-hi. We therefore conclude that including radiation
pressure in the equation of state does not lead to large differences.

The M-hi model has an in-spiral timescale shorter than for all the other
simulations, including the one at lower resolution. The reason for this
behaviour are (i) the larger size of the star at the time of the start of the
simulation, (ii) the slow radial expansion of the initial model and (iii) the
delivery of recombination energy to the gas in the outer layers which promotes
its expansion (and some early envelope unbinding - see Section 4.2 for
additional details). The exact role played by these three factors on the
faster in-spiral on-set, in-spiral shape and final separation is not well
quantified. However, we rerun simulation M-hi but this time starting from a
stellar structure that was larger (285 R⊙), but more stable and the in-spiral
on-set, shape and final separation were remarkably similar to those associated
with M-hi, showing that a star that is smaller, but slowly expanding results
in a similar CE as one simulated using a stellar structure that is not
expanding, but is already somewhat larger (see Appendix C for more details).
This demonstrates that the CE interaction characteristics are not strongly
dependent on the fine details of the initial stellar structure.

We note that M-hi does not just in-spiral sooner, but also has a completely
different in-spiral shape compared to the other three simulations, with an
overall shallower slope throughout (Figure 6, right panel). We will discuss
this result in Section 5 where we will contextualise it with other work where
similar comparisons are carried out.

The post in-spiral separation of M-hi is clearly larger as can be seen in
Figure 6. The values at the end of the plunge-in is 35 R⊙, 42 per cent larger,
than the corresponding value for I-hi. Once I-hi and M-hi are aligned at the
time of steepest in-spiral, the separation of M-hi 4.4 years after the
steepest in-spiral time is 34.6 R⊙ compared to the equivalent value for I-hi
of 24.2 R⊙. As is also clear from Figure 6, right panel, M-hi stands out in
having larger final separation compared to the other three simulations.

Since the orbits are still decreasing by the end of the simulation, we
extrapolate the final separations using a exponential decay fit (Table 3,
column 7). The details of the extrapolation method are described on Appendix
D. The two values correspond to the upper and lower limits for the final
separations. While these extrapolated values indicate that the orbits will
further shrink, the differences are not significant to influence the envelope
or considering a merger scenario. There is one caveat about this extrapolation
technique. The number of oscillations in the orbital evolution after the
plunge-in have to be large enough to find a decay factor that fits the decay
of the orbit at later times, otherwise the upper limit of the extrapolation
may be larger than the final separation in the simulation. This is the case
for the orbital evolution in I-lo and Rei18. However, the rest of the AGB
models have extrapolated (upper and lower) values that are below the values in
column 6.

Lastly, the post-plunge-in orbit for the M-hi model does not present an
eccentricity as large as the other models (0.005, Figure 6, right panel and
Table 3).

For all the simulations, the time between Roche lobe overflow and orbital
stabilisation is at least one order of magnitude smaller than the thermal
pulse of the TP-AGB star, from the moment the star reaches 260 R⊙ until it
contracts to that same size, $400$ years later. This suggests that the common
envelope evolution can easily take place during a thermal pulse. The only
question is whether, in nature, the time between the on-set of Roche lobe
overflow and the on-set of the in-spiral could be long enough that the CE
would not take place, because the star contracts before in-spiral could start
in earnest. After all, we do not have a good grasp on this pre-plunge-in time.
We discuss this further in Section 5.

Figure 7: Bound mass as a fraction of the envelope mass for all simulations,
shifted so that $t=0$ corresponds to the time of plunge-in, $t_{\rm max,pl}$.
The square symbols show the times compared in Section 4.4. For these plots,
unbound mass is defined as $e_{\rm pot}+e_{\rm K}>0$.

Figure 8: Bound mass (coloured curves) and binary separation (black curves)
plots of I-hi (left panel) and M-hi (right panel) models. We plot the bound
mass using the mechanical (solid lines) and thermal (dashed lines) criteria in
each panel.

###  4.2 Analysis of the unbound mass

In Figure 7 we present the bound mass as a function of time for our four
simulations as well as for that of Reichardt et al. (2019), using the most
stringent criterion for calculating unbound mass: $e_{K}+e_{\rm pot}\geq 0$,
where $e_{K}$ and $e_{\rm pot}$ are the kinetic and potential energies of all
the gas particles at each moment in time. The potential energy for each SPH
particle was calculated with respect to the rest of the SPH particles and each
point mass particle. As extensively discussed in several publications (e.g.,
Iaconi et al., 2017; Staff et al., 2016), including the thermal energy (and
even more importantly the recombination energy) in this criterion assumes that
the entire thermal energy payload of the stellar envelope (including the
thermalised recombination energy) will, at some point, be transformed into
bulk kinetic energy. Since this is not necessarily the case, it is good
practice to retain a more stringent “mechanical” criterion as well as noting
the amount of additional unbound gas when adopting a thermal energy criterion.

In order to use a criterion that includes the thermal energy, we follow the
procedure described below. The thermal energy calculation depends on the
equation of state. For I-hi the thermal energy is provided by phantom in the
form of the internal energy. For the mesa and ideal plus radiation pressure
equations of state, we calculate the thermal energy in the following way:

$E_{\rm th}=\frac{3}{2}\frac{k_{\rm B}T_{i}}{m_{\rm h}\mu}+\frac{a_{\rm
rad}T_{i}^{4}}{\rho_{\rm i}}.$ (5)

where $k_{\rm B}$ is the Boltzmann constant, $m_{\rm h}$ is the mass of the
hydrogen atom, $\mu$ is the mean molecular weight, $T_{i}$ is the temperature
of the i-th particle, $\rho_{i}$ is its density and $a_{\rm rad}$ is the
radiation constant. The temperature and mean molecular weights are output by
the computation.

Figure 8 compares the unbound mass calculated with the criteria that omit, or
include the thermal energy component, for the I-hi and M-hi simulations. As
expected, the mechanical criterion has less unbound gas than the thermal
energy one. The bound mass at the end of I-hi is $0.83$ and $0.82$ M⊙, using
the mechanical and thermal criteria, respectively. For M-hi, these values are
$0.2$ and $0.08$ M⊙, respectively. Therefore, I-hi unbinds $27-28$ percent of
the envelope, while M-hi unbinds $82-93$ percent. The curves in Figure 8 are
all trending down, implying additional unbinding. However, as we are about to
see, one has to be cautious about the later part of the unbinding.

(a) I-hi model, mechanical criterion.

(b) I-hi model, thermal criterion.

(c) M-hi model, mechanical criterion.

(d) M-hi model, thermal criterion.

Figure 9: Newly unbound particles for the I-hi (top row) and M-hi (bottom row)
models using the mechanical (left column) or thermal (right column) unbound
mass criteria. The color grading indicates the amount of gas that is newly
unbound, approximately every 5 days. The black curve represents the binary
separation. The y-axis measures distance from the primary’s point mass core
particle.

Figure 9 shows four 2D histograms representing the radial position and time of
newly unbound SPH particles for I-hi (top panels) and M-hi (bottom panels),
using a heat map to represent the amount of mass that becomes unbound (in
arbitrary units). We can divide the unbinding into four stages. The first
stage is only clearly seen in M-hi: some particles become unbound starting at
$\sim$1 yr, likely promoted by some early expansion and the delivery of
recombination energy in those surface layers. This would promote additional
expansion and precipitate the in-spiral, as explained in Section 4.1.

The second phase, more prominent in I-hi, is due to the transfer of orbital
energy to the gas. It appears as a band of particles right above the black
curve that denotes the orbit. This unbinding starts before the plunge-in and
it continues and intensifies in the early phase of the plunge-in, but it stops
at approximately the time of maximum in-spiral at $t\sim 15$ years for I-hi
and $t\sim 5$ years for M-hi. Gas unbound near the companion and moving
outward transfers energy to gas further out leading to more unbinding – this
is seen as the curving streaks of unbound particles in Figure 9.

The third phase is denoted by a very localised region of unbound mass between
$\sim$15 years ($\sim$5 years) and $\sim$22 years ($\sim$10 years) for the
I-hi (M-hi) simulation in Figure 9. Particles are no longer becoming unbound
near the companion because the companion is closer to the core, where higher
density, more bound gas resides. Instead, the gas unbound in the previous
phase collides and transfers some of its kinetic, and to a lesser degree
thermal, energy to the bound gas in the upper layers of the envelope,
effectively unbinding them. The M-hi simulation differs from I-hi in this
phase: the “streak" of unbound gas is very dense (more unbound mass) with most
of the mass being unbound here rather than in the previous phase as is the
case for I-hi. This is likely due to the release of recombination energy.

The fourth and last phase occurs, for the I-hi simulation at $\sim 22$ years
and for M-hi at $\sim$10 years. Here the region of newly unbound particles
extends again from deeper layers to near the surface. This fourth phase
corresponds to the second drop in the bound mass for I-hi in Figure 7 (left
panel). For M-hi the decline in bound mass in the same Figure is more gradual
and this fourth phase is not as distinct. The nature of this late unbinding is
related to a decrease in resolution near the cores. This event is resolution
dependent as can be seen by comparing I-hi and I-lo in Figure 7 and was
already described by Reichardt et al. (2020). When the SPH particles’
smoothing lengths become comparable to the softening length of the point mass
particles, the gas particles fail to resolve the core — due to the modified
point mass gravitational potential — and the pressure can become negative.
This causes an acceleration of particles near the point masses, which at this
point have a separation of $\lesssim 30$ R⊙. These particles may or may not
become unbound, but even if not, they transfer their energy to upper layers of
the envelope that, being only loosely bound, can become unbound. This
resolution-dependent unbinding is far more obvious in ideal gas simulations
and the extent to which the late unbinding in M-hi is resolution dependent has
not yet been established.

The panels including thermal energy (right panels in Figure 9) are similar,
but in the fourth phase show additional regions where gas is becoming newly
unbound. For I-hi this difference is marginal and shock heating may be the
contributing agent. For M-hi the extra unbinding seen using the thermal
criterion is more widespread and likely indicates the release of recombination
energy (which is assumed to be immediately thermalised).

We explore further the resolution-dependent unbinding in Figure 10, where we
show perpendicular slices of total energy (mechanical and thermal energy) at
times when I-hi and M-hi display resolution-dependent unbinding. A bi-conic
region of unbound gas becomes evident at the approximate time when the phase
of additional, resolution-dependent unbinding starts. The geometry of the CE
ejection is such that these regions have lower density than the density closer
to the equatorial plane. This supports our interpretation that the particles
near the core acquire higher kinetic energies but only those in the lower
density regions can more readily escape as their path is not as impeded by
high density overlaying particles. As reference, the specific potential energy
at the surface of the models at time zero is $4.1\times 10^{12}$ erg g-1. The
lack of one of the two cones in I-hi (Figure 10 \- left panel) is due to the
lack of perfect symmetry in the density distribution.

Figure 10: Slices in the xz-plane for the I-hi (left) and M-hi (right) models
at times 27.4 and 12.6 years, when the resolution-dependant unbinding is under
way. The rendering correspond to the total specific energy of the gas, defined
as the sum of the specific kinetic, potential, and thermal energies, with red
and blue colours representing unbound and bound gas, respectively (using the
thermal energy criterion).

###  4.3 Rotation and angular momentum

The $z$ component of the initial angular momentum of the system with respect
to the centre of mass is $9.55\times 10^{52}$g cm2 s-1 and it mostly belongs
to the orbital angular momentum of the stellar cores. By the end of the
interaction, most of the orbital angular momentum is transferred to the
envelope, some of which becomes unbound (Figure 16): $41\%$, $41\%$ and $88\%$
of the $z$-angular momentum unbound, for the I-hi, R-hi and M-hi, respectively
(using the mechanical energy criterion).

If the envelope is only partly ejected by the dynamical plunge-in, the
remaining angular momentum will play a key role in the adjustment of the
remaining envelope into a new equilibrium configuration. Indicatively, if the
remaining bound angular momentum at the end of the I-hi simulation (which
unbinds only 27-28 per cent of the envelope) were given to the initial I-hi
model, and the star were to rotate as a solid body, its surface velocity would
be $65$ km s-1, well above the local escape speed. This demonstrates that as
CE mergers (that have not ejected a great amount of gas) attempt to readjust
thermally, they have to contend with a great deal of angular momentum
redistribution which will greatly impact the future structure of the star.

In Figure 11 we analyse the gas density and velocity field (top panels), as
well as the angular velocity (bottom panels) for M-hi. In the top panels, the
Roche lobes are shown - the primary core is surrounded by higher density. At
$t=5.8$ and $6$ years the gas in the primary’s Roche lobe moves approximately
in the same direction as the orbit as expected, while the gas in the
secondary’s lobe orbits the point mass companion, as also expected after
entering the lobe via L1. The corresponding bottom panels show an angular
velocity contrast pattern that is in line with what we see in the top panels.
At $6.7$ years, right at the end of the plunge-in, the region within and
around the lobes is rotating with the orbit and is close to corotation (green
colour). Any remaining patterns of differential rotation are quickly
disappearing. By $7.41$ years the gas in and around the orbit is mostly in
corotation. We conclude that the end of the in-spiral in these simulations is
primarily the result of the inner gaseous region being brought into corotation
by the binary cores.

Figure 11: Slices of density including velocity arrows (top row) and
‘corotation’ (bottom row) for the M-hi simulation. We show the Roche lobes
corresponding only to the point mass particles system in the top panels.
Corotation is defined as $\Delta\omega/\omega_{\rm bin}$, where
$\Delta\omega\equiv\omega_{\rm bin}-\omega_{\rm gas}$ and $\omega_{\rm bin}$,
$\omega_{\rm gas}$ are the angular velocity of the point mass particles and
the of the gas, respectively, with respect to the centre of mass of the two
point masses, only. The angular velocity, $\omega$, is defined as positive for
anti-clockwise motion. The white dots represent the point masses. The times on
the top panels correspond to the times labelled on the bottom panels. These
times are also marked by magenta symbols in Figure 6.

###  4.4 The morphology of the ejected envelope

Figure 12: Density slices of Rei19 (left column), I-hi (central column) and
M-hi (right column) at three zoom levels in the x-z plane: top and middle rows
show two zoom levels using the same density colour scale, while the bottom row
zooms in further using a colour scale that shows only the densest parts of the
structure. The first row has a box size of $1.6\times 10^{4}$ R⊙ per side. The
second row has sides of 4000R⊙on each box. The bottom row is has dimensions of
$625$ R⊙.

Approximately 20 per cent of all planetary nebulae (PNe) have a surviving post
CE binary in the middle (e.g., Jacoby et al., 2021), indicating that the
nebulae are the ejected envelopes. The great range of post-CE PN shapes is in
stark contrast to the relative homogeneity of the post-CE binaries at the
centres of the nebulae. This suggests that post-CE PN morphology is very
sensitive to small changes in stellar and binary parameters.

García-Segura et al. (2018), Frank et al. (2018) and Zou et al. (2020) were
first to simulate PN formation in the aftermath of a CE interactions. However,
they used, as initial conditions, CE simulations with RGB, instead of AGB
stars; in addition, those CE simulations had only partly ejected their
envelopes. One of the aims of this study is to generate an AGB CE simulation
that fully ejects the envelope to be used as input for a new PN simulation.
Here, we therefore compare the I-hi and M-hi simulations’ gas distributions
with that of the simulation of Reichardt et al. (2019) that was already used
to simulate the PN formation by Frank et al. (2018) and Zou et al. (2020). We
use this comparison to predict whether our new simulation is expected to
result in substantially different nebular shapes.

We start by selecting a post-CE time for I-hi and M-hi at which we can compare
the CE structure with the simulation of Reichardt et al. (2019) adopted by Zou
et al. (2020). We do not know the time after the CE ejection when the fast,
post-AGB wind that ploughs into the CE gas should start blowing. We can
presume that the fast wind will start when the phtospheric temperature reaches
a certain value ($\sim$20 000 K). This time, in turn, depends on the
adjustment of the core once the envelope has departed. None of these processes
are accurately modelled by the CE simulation. Zou et al. (2020) therefore used
a CE simulation that was stopped about one year after the time of maximum in-
spiral.

Somewhat arbitrarily, we select the dynamical time of the primary star prior
to the start of the CE interaction as a “yardstick”. We then assume that the
CE ejection takes place at the time of maximum in-spiral and call that time-
zero: $t_{0}=t_{\rm max,pl}$ (Table 3). The Reichardt et al. (2019) simulation
was used to generate the PN at 14.1 years of simulation time and 0.8 years
after $t_{0}$. With a dynamical time for that star of $\tau_{\rm dyn,RGB}$ =
11 days, the simulation was stopped $26\times\tau_{\rm dyn,RGB}$ after
$t_{0}$. Using the same factor, together with the AGB star dynamical time of
$\tau_{\rm dyn,AGB}=41$ days, we stop the simulation for the purpouse of
comparison, 2.9 years after $t_{0}$, or $18.2$ and $7.58$ yr from the start of
the I-hi and M-hi simulations, respectively. These times are marked with
square symbols in Figures 6 (right panel) and 7.

We finally compare the three simulations in Figure 12. We observe that our AGB
CE simulations lead to a somewhat more extended and more spherical post-CE gas
distribution (top row).

At higher magnification (Figure 12, middle row) all three simulations show a
peanut shaped overdensity, with the simulation of Reichardt et al. (2019)
being the most pronounced. This shape leaves a reasonably evacuated corridor
in the vertical direction, which is responsible for hydrodynamically
collimating the fast wind (Zou et al., 2020). In Figure 12, bottom row, we
zoom further into the centre of the gas distribution and change the density
colour bar to highlight the structure there. Clearly, at the highest densities
($\log\rho/{\rm g~{}cm}^{-3}=7-7.3$) the M-hi is far more spherical than the
two ideal gas simulations.

The velocity field shows some striking differences between the three
simulations. Overall, the M-hi simulation seems to show faster outflow
velocities with an inflow at the poles at smaller scales (Figure 12, lower
row). The RGB simulation of Reichardt et al. (2019) has fast outflow
velocities (including at the poles) only at intermediate scales (Figure 12,
middle row), while the I-hi AGB simulation has outflow velocities comparable
with M-hi at large scale, but at smaller scale the velocity field is slow and
confused. Overall, the simulation that unbinds the envelope (M-hi) appears to
have a faster outflow.

In conclusion we propose that M-hi and other simulations like it, that include
recombination energy, are the most realistic and hence the most appropriate
simulations from which to derive PN shapes. It will be then interesting to see
whether the strong collimation observed by Zou et al. (2020) is still a
feature, something that we leave for future work.

###  4.5 Dust and the expanding common envelope

The dust formation that is likely to happen in the dense, yet expanding and
cooling common envelope has never been calculated consistently. Iaconi et al.
(2019) and Iaconi et al. (2020) have, however, post-processed the dust
formation in the RGB, CE simulations of Reichardt et al. (2019) and found that
a CE event like that would produce approximately 10-3 M⊙, of dust
independently of whether the dust is carbon- or oxygen-rich.

For AGB donor stars, the formation of dust in the expanding CE is likely to
play an even greater role than for RGB giants, impacting the objects’ optical
properties as well as the dynamics of the expansion. In our simulations we
have not considered dust formation. However, we can infer where the dust
particles would likely form based on the gas temperature and density from
which our tabulated equation of state tables provide knowledge of the expected
opacity. In Figure 13 we show horizontal and vertical slices of the opacity
for M-hi at $t=7.58$ years. This figure is similar to one by Reichardt et al.
(2020, figure 10) for their RGB star. The two black contours represent the
recombination temperature for hydrogen ($6000$ K, inner contour) and an
approximate temperature for dust condensation ($1500$ K, outer contour).
Inside the $6000$ K contour, hydrogen is ionised and the opacity is high and
dominated by electron scattering. Outside of the $T=6000$ K region, hydrogen
is recombined and the medium becomes transparent. However, as we move
outwards, the temperature decreases and when it reaches $\sim$1500 K,
molecular and dust opacity contribute to increasing the overall opacity once
again.

Finally we observe that in the 0.88 M⊙ RGB star of Reichardt et al. (2019),
the high opacity shell is concentrated around the equator, with lower opacity
along the polar direction. On the other hand, their 2 M⊙ RGB giant model shows
higher opacity in the direction of the polar caps, with lower opacity at the
equator. In our case, M-hi shows a thick shell of higher opacity all around
the star, with a somewhat oblate shape similar to the 1-M⊙ simulation of
Reichardt et al. (2019). In future work we will calculate the dust mass
associated with this CE and infer the dust driving properties of the dusty
envelope.

Figure 13: Top row: opacity slices in the orbital (left panel) and
perpendicular (right panel) directions for the simulation M-hi. Bottom row:
equivalent slices of density with, overlaid, projected velocity arrows.

##  5 Discussion

We have presented four simulations of a common envelope interaction using a
TP-AGB star donor. We have used the 1D, stellar structure and evolution code
mesa to calculate a 2 M⊙, ZAMS star that we have evolved until the seventh AGB
thermal pulse, when it had a radius of 260 R⊙ and a mass of 1.7 M⊙. We have
implemented two different relaxation procedures, one used by Ohlmann et al.
(2017) and the other designed by Lau et al. (2022b), to map and relax the
stellar model in the 3D smoothed particle hydrodynamics (SPH) code phantom.
The relaxed star was set up with a $0.6$-M⊙ companion at an orbital distance
of 550 R⊙, such that the primary was Roche lobe-filling. Moreover, we have
simulated the system using three different equations of state: ideal gas,
ideal gas plus radiation and a tabulated form that takes recombination energy
into account. All simulations were run with a resolution of $1.3\times 10^{6}$
SPH particles. For the ideal gas model, we also run a simulation at low
resolution, with $1.37\times 10^{5}$ SPH particles.

The final orbital separation ranges between $20$ (I-hi) and $31$ R⊙ (M-hi),
depending on the equation of state. The adoption of a tabulated equation of
state results in a final separation that is $55$ per cent larger than for the
simulation with an ideal gas equation of state. Reichardt et al. (2020) found
almost no difference in the final separations between the two equations of
state for their $0.88$ M⊙, RGB simulation and for their $2$ M⊙, RGB
simulation. Sand et al. (2020) found instead an increase between $11$ and 25
per cent (depending on companion mass) in final separation when using a
tabulated equation of state in their 1.2 M⊙ AGB star.

In addition to the final separation, the shape of the inspiral curve in the
separation vs. time plot (Figure 6), changes dramatically between the
simulations I-hi and M-hi: not only does the M-hi simulation in-spiral sooner,
but the slope of the in-spiral is overall less steep. Reichardt et al. (2020)
observed virtually no difference in the in-spiral shape between ideal gas and
tabulated equation of state simulations. Sand et al. (2020), on the other
hand, found that the early in-spiral is similar for the two equations of
state, but that the late in-spiral is different, with the tabulated equation
of state resulting in a less steep curve, as we also find.

We interpret these differences as due to early expansion of the stellar
structure because of a slight lack of hydrostatic equilibrium, as explained in
Section 4.2 and to the recombination energy delivered upon early expansion
aiding that expansion. Our starting the simulations from the time of Roche
lobe overflow, affords more time before the in-spiral than for the simulations
of Reichardt et al. (2020), which were started on the surface. During this
time the M-hi star “expresses" its tabulated equation of state by enacting
additional expansion driven by the outer layers recombining (Section 4.2 and
Figure 9). This leads to a relatively faster in-spiral on-set, compared to
I-hi. As the companion in-spirals and further recombination energy is
delivered, a gentler inspiral slope is observed, as the densities encountered
are lower than for I-hi. Finally, we see a wider separation, as the star has
been effectively “helping" with its own unbinding.

Sand et al. (2020) started their simulations with the companion farther away
from the primary surface, but closer than Roche lobe overflow, such that they
observed more similarity in the in-spiral between ideal and tabulated
equations of state, but the shape of the in-spiral diverged later on with the
tabulated equation of state giving rise to a gentler slope, in a manner not
dissimilar to us222This is generally true for their 6 simulations, but there
is some variability across the three companion masses that complicated the
analysis further.. The simulations of Reichardt et al. (2019) were instead
started with the companions on the surface of the stellar structures, such
that the in-spiral was initiated immediately, giving no time to the star to
expand and for the recombination energy to play a role. Their envelope
ejection took place afterwards. This likely explained the similarity of their
in-spirals and final separations. We conclude preliminarily that the equation
of state does play a role on the envelope dynamics and even on the final
separation, contrary to what was concluded by Reichardt et al. (2020).

Our final orbital separations (20-31 R⊙) are wider than encountered in the
bulk of observed post-CE central binaries in PN (Iaconi & De Marco, 2019).
This is so also for the 1 M⊙, AGB star simulations of Sand et al. (2020, 22-69
R⊙, depending on the companion’s mass) and even for the 1.78 M⊙, AGB star
simulation with a 0.98 M⊙ companion by Chamandy et al. (2020, 15.5 R⊙, and we
note they used an ideal gas equation of state, suggesting that the separations
would be even larger had they adopted a tabulated equation of state).

The in-spiral itself is very brief, taking place on the stellar dynamical
timescale and about an order of magnitude shorter than the thermal pulse.
However, we could ask whether the period of Roche lobe overflow preceding the
in-spiral, which could be longer than we have modelled, may be long enough
that an in-spiral may not happen at all before the star contracts at the end
of the thermal pulse. We propose here that the information gathered so far
indicates that once a thermal pulse leads to Roche lobe overflow the expanding
envelope of the plusating star will force a CE in-spiral. The natural
expansion of the thermally pulsating AGB star radius is not dissimilar to our
gently expanding structures (that expand due to a slight violation of
hydrostatic equilibrium; Figure 3). Additional simulations with increasing
accuracy of the stellar structure will be needed to finalise this conclusion.

In Figure 8, it is clear that even adopting the most stringent, mechanical
criterion, M-hi unbinds a great deal more envelope than I-hi or R-hi.
Including thermal energy in the criterion, results in an even greater amount
of unbound gas. On the other hand, the inclusion of thermal energy makes
almost no difference to the unbound mass calculated for I-hi. This is likely
because there is much less marginally bound material in the I-hi model that
would be marked as unbound with the addition of a small amount of thermal
energy. The overall trend in the bound gas curve over time for M-hi indicates
that the entire envelope could easily be unbound in a short time.
Interestingly, this trend is very different to what was observed by Reichardt
et al. (2020) in their RGB star simulations: their 0.88-M⊙ simulation with
ideal gas and tabulated equation of state behaves similarly to ours. However,
their 1.8 M⊙ simulation unbound very little gas, even when using a tabulated
equation of state showing that their recombination energy had not been fully
released by the end of that simulation.

There is relatively little difference between I-hi and R-hi, where the latter
includes radiation energy. While this is not unexpected for low mass giants,
we note here some subtle differences. R-hi starts to in-spiral sooner than
I-hi, likely due to the same reasons argued above, namely the initial model
was ever so sightly more out of hydrostatic equilibrium than I-hi. The shape
of the in-spirals, however, are very similar. The separation at the end of the
plunge-in, $a_{\rm f,pl}$, and the separation at the end of the simulations,
$a_{\rm f}$, are both larger for R-hi (44 and 27 R⊙) than for I-hi (32 and 20
R⊙) and take place 1.2 years after the time of maximum plunge in, $t_{\rm
max,pl}$, compared to 1.5 years for I-hi. The unbound mass at the end of the
in-spiral (which excludes the resolution-dependent unbinding at the end) is
almost identical for I-hi and R-hi.

Towards the end of all simulations, but particularly in I-lo/I-hi and R-hi,
there appears to be a second envelope unbinding phase (Figure 7), no matter
what criterion is used (Figure 9). Reichardt et al. (2019) described this
resolution-dependent unbinding as caused by the decrease in density (and
concomitant increase in the local SPH particles smoothing lengths) in the
proximity of the core. When the local SPH particles no longer “resolve" the
core, the pressure gradient inverts and they evacuate the core region. When
they confer their kinetic energy to surrounding particles, these particles may
become unbound. Figure 10 indicates that this late, resolution-dependent
ejection occurs in the direction of z-axis. This zone has a lower density and
would facilitate the escape of particles. For M-hi, this phase of resolution-
dependent unbinding is less obvious from Figure 7, but we determined that it
does happen at about 10 years from the start of the simulation. While this
reduces the fraction of unbound envelope that we can measure with certainty,
the gas unbinding ratio observed in Figure 7 for M-hi still indicates that the
envelope will soon become fully unbound: $15.6$ and $16.7$ years with the
final mass-loss rate for mechanical and thermal criteria, respectively.

Morphologically and kinematically, the ejected CE of this AGB star is somewhat
different from that of the RGB star of Reichardt et al. (2019, Figure 12). The
innermost structure of our M-hi simulation is less elongated than the
equivalent region in the RGB simulation of Reichardt et al. (2019) \- this
inner mass distribution influences greatly the nascent PN. In future work we
will use the ejected CE of M-hi to simulate the PN and will further analyse
differences at that time. The opacity within the envelope also indicates that
dust formation would take place. This would likely result in radiation driving
and unbinding further envelope gas.

Finally, some considerations on the common envelope unbinding efficiency.
Following Lau et al. (2022b) we can determine the common envelope efficiency
parameter based on our three high resolution simulations, which have
individual values for the orbital energy delivered, $\Delta E_{\rm orb}$
(where we use the average of the apastron and periastron fitted values from
Table 3 to calculate it). We also use a definition of the binding energy that
includes thermal energy. If thermal energy were omitted, the values of
$\alpha$ would increase by a factor of two. In Table 4 we list $\alpha$ for
each of the three high resolution simulations. The values of $\alpha$ larger
than unity tell us that the relatively meager orbital shrinkage in the
simulations did not deliver enough energy to unbind the envelope, according to
the two binding energy definitions. The reason why the in-spiral stalls before
delivering sufficient energy to unbind the envelope is that our simulations
also contain recombination energy, which aids the envelope ejection and which
is not accounted for in the formalism.

Lau et al. (2022b) used three definitions of the binding energy. One without
thermal energy, only includes gravitational potential energy. The second, is
the same we have used, while the third includes recombination energy ($E_{\rm
bin}=E_{\rm grav}+E_{\rm th}+E_{\rm rec}$). Adopting a binding energy
definition that excludes thermal energy would increase $\alpha$ by a factor of
two. We cannot use the third definition here: the inclusion of recombination
energy in the definition of binding energy leads to a formally unbound
envelope and a negative value of $\alpha$. The negative value of $\alpha$
simply states that given that the binding energy including recombination is
negative (the envelope is formally unbound) the orbit needs to expand instead
of contracting. Clearly this means that the alpha formalism does not work
here. The reason is that while the envelope might be formally unbound at time
zero, when the in-spiral commences, the thermal and recombination energies
have not yet turned into bulk kinetic energy so the companion can still feel a
drag force, spiral in, and deliver orbital energy. This formalism does not
include this time sequence and fails.

Sim (EoS) | $\alpha$
---|---
I-hi (Ideal) | 1.03
R-hi (Gas+rad) | 1.24
M-hi (Full) | 1.76
Table 4: The common envelope efficiency parameters derived from the three high
resolution simulations. The binding energy includes thermal energy. The
delivered orbital energy uses the mean of periastron and apastron final
orbital separations as listed in Table 3.

##  6 Conclusions

We simulated the common envelope interaction between a 1.7 M⊙, thermally
pulsating AGB star with a radius of 250 R⊙, and a 0.6-M⊙ companion started
approximately at the time of Roche lobe overflow using an ideal gas, ideal gas
plus radiation and tabulated equations of state. Our conclusions are:

  1. 1.

The thermalisation of recombination energy released during common envelope
expansion is sufficient to unbind most of the envelope. This is consistent
with previous simulations that include recombination energy (Nandez & Ivanova,
2016; Sand et al., 2020; Ivanova & Nandez, 2016; Lau et al., 2022b). We assume
the gas opacity to be sufficiently high to absorb the photons emitted from the
recombination process before they can escape the envelope. While this may not
be true for all of the hydrogen recombination photons (albeit some fraction of
them may be thermalised), there is reason to believe that the entirety of
helium recombination energy can be transferred to the envelope through this
process (Ivanova, 2018; Lau et al., 2022a).

  2. 2.

The inclusion of recombination energy leads to a dynamical, but shallower in-
spiral and a wider final separation, by about 58 per cent (or 66 per cent
using the lower limit extrapolated values). This is a similar trend to that
obtained by Sand et al. (2020), who modelled a low mass AGB star and three
companions, where the final separations were larger by between 20 and 30 per
cent with a tabulated EoS, and by Lau et al. (2022b) who modelled a 12 M⊙ red
supergiant and obtained separations that were wider by 16% when recombination
energy was included. Reichardt et al. (2020) modelled two RGB donor stars and
obtained separations that were wider for their lower binding energy star (by
16 per cent), but ever so slightly smaller for their more bound RGB star (by 4
per cent). Clearly this conclusion does depend on the specific parameters of
the interaction, but the reason for a wider separation when recombination
energy is included is likely to be that recombination energy expands the
stellar envelope more and earlier than is the case when using an ideal gas
EoS, which contributes to less in-spiral.

  3. 3.

Thermal pulses are likely to trigger a common envelope in-spiral, so long that
the primary giant overfills its Roche lobe when the star is in the expanding
phase of the thermal pulse. Doubling the resolution just increases the RLOF
timescale by 18 per cent or by 2 years. While a third simulation with another
increment in resolution is needed to do a proper convergence test, the RLOF
timescale is still under one order of magnitude smaller that the duration of
the pulse. In addition, in Nature, the star is expanding between 250 and 300
R⊙ in less than 100 years, forcing a CE even if, under a constant radius, the
stable mass transfer would last longer. We conclude that thermal pulses
increase the ability of AGB stars to capture companions into common envelopes,
thereby increasing the relative number of post-AGB, post-CE systems relative
to a situation where the thermal pulse does not cause a CE to take place.

  4. 4.

At the large scale, the CE shape shortly after the in-spiral is somewhat more
spherically distributed for AGB CEs than for RGB CEs and even more so for AGB
CEs calculated with the inclusion of recombination energy.

  5. 5.

While formally calculating the CE efficiency fails, we note that the wide
final separations obtained in our AGB CE simulation are in disagreement with
the smaller separations seen in observations (Iaconi & De Marco, 2019). This
leaves wide open the question of how AGB CE result in final separations of the
order of a few solar radii.

## Acknowledgements

M. G. B. acknowledges funding support from Macquarie University through the
International Macquarie University Research Excellence Scholarship (‘iMQRES’)
and to Ilya Mandel for his insightful commentaries during the course of this
research. OD acknowledges funding from Australian Research Council Discovery
Project, DP210101094. The authors thank to the Referee for the constructive
comments and recommendations which helped to improve the quality of the paper.
Parts of this research work were performed on the Gadi supercomputer of the
National Computational Infrastructure (NCI), which is supported by the
Australian Government. This research has made use of the NASA’s Astrophysics
Data System (ADS) and NASA Exoplanet Archive, which is operated by the
California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.

This is a pre-copyedited, author-produced PDF of an article accepted for
publication in Monthly Notices of the Royal Astronomical Society following
peer review. The version of record is available online at:
https://doi.org/10.1093/mnras/stac2301.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Carroll-Nellenback et al. (2013) Carroll-Nellenback J. J., Shroyer B., Frank A., Ding C., 2013, Journal of Computational Physics, 236, 461
  * Chamandy et al. (2020) Chamandy L., Blackman E. G., Frank A., Carroll-Nellenback J., Tu Y., 2020, MNRAS, 495, 4028
  * De Marco & Izzard (2017) De Marco O., Izzard R. G., 2017, Publ. Astron. Soc. Australia, 34, e001
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * Frank et al. (2018) Frank A., Chen Z., Reichardt T., De Marco O., Blackman E., Nordhaus J., 2018, Galaxies, 6, 113
  * García-Segura et al. (2018) García-Segura G., Ricker P. M., Taam R. E., 2018, ApJ, 860, 19
  * García-Segura et al. (2020) García-Segura G., Taam R. E., Ricker P. M., 2020, ApJ, 893, 150
  * García-Segura et al. (2021) García-Segura G., Taam R. E., Ricker P. M., 2021, ApJ, 914, 111
  * Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
  * Hirai et al. (2020) Hirai R., Sato T., Podsiadlowski P., Vigna-Gómez A., Mandel I., 2020, MNRAS, 499, 1154
  * Iaconi & De Marco (2019) Iaconi R., De Marco O., 2019, MNRAS, 490, 2550
  * Iaconi et al. (2017) Iaconi R., Reichardt T., Staff J., De Marco O., Passy J.-C., Price D., Wurster J., Herwig F., 2017, MNRAS, 464, 4028
  * Iaconi et al. (2018) Iaconi R., De Marco O., Passy J.-C., Staff J., 2018, MNRAS, 477, 2349
  * Iaconi et al. (2019) Iaconi R., Maeda K., De Marco O., Nozawa T., Reichardt T., 2019, MNRAS, 489, 3334
  * Iaconi et al. (2020) Iaconi R., Maeda K., Nozawa T., De Marco O., Reichardt T., 2020, MNRAS, 497, 3166
  * Ivanova (2018) Ivanova N., 2018, The Astrophysical Journal Letters, 858, L24
  * Ivanova & Nandez (2016) Ivanova N., Nandez J. L., 2016, MNRAS, 462, 362
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Jacoby et al. (2021) Jacoby G. H., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 5223
  * Kippenhahn & Weigert (1990) Kippenhahn R., Weigert A., 1990, Stellar Structure and Evolution. Springer
  * Lattanzio & Karakas (2001) Lattanzio J. C., Karakas A. I., 2001, Memorie della Societa Astronomica Italiana, 72, 255
  * Lau et al. (2022a) Lau M. Y., Hirai R., Price D. J., Mandel I., 2022a, arXiv preprint arXiv:2206.06411
  * Lau et al. (2022b) Lau M. Y., Hirai R., González-Bolívar M., Price D. J., De Marco O., Mandel I., 2022b, Monthly Notices of the Royal Astronomical Society, 512, 5462
  * Madappatt et al. (2016) Madappatt N., De Marco O., Villaver E., 2016, MNRAS, 463, 1040
  * Miszalski et al. (2009) Miszalski B., Acker A., Moffat A. F. J., Parker Q. A., Udalski A., 2009, A&A, 496, 813
  * Nandez & Ivanova (2016) Nandez J. L., Ivanova N., 2016, Monthly Notices of the Royal Astronomical Society, 460, 3992
  * Nandez et al. (2014) Nandez J. L., Ivanova N., Lombardi Jr J. C., 2014, ApJ, 786, 39
  * Ohlmann et al. (2017) Ohlmann S. T., Röpke F. K., Pakmor R., Springel V., 2017, A&A, 599, A5
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Price (2008) Price D. J., 2008, J. Comp. Phys., 227, 10040
  * Price et al. (2018) Price D. J., et al., 2018, Publ. Astron. Soc. Australia, 35, e031
  * Rasio & Livio (1996) Rasio F. A., Livio M., 1996, ApJ, 471, 366
  * Reichardt et al. (2019) Reichardt T. A., De Marco O., Iaconi R., Tout C. A., Price D. J., 2019, MNRAS, 484, 631
  * Reichardt et al. (2020) Reichardt T. A., De Marco O., Iaconi R., Chamandy L., Price D. J., 2020, MNRAS, 494, 5333
  * Ricker & Taam (2012) Ricker P. M., Taam R. E., 2012, ApJ, 746, 74
  * Sand et al. (2020) Sand C., Ohlmann S. T., Schneider F. R., Pakmor R., Röpke F. K., 2020, A&A, 644, A60
  * Sandquist et al. (1998) Sandquist E. L., Taam R. E., Chen X., Bodenheimer P., Burkert A., 1998, ApJ, 500, 909
  * Schönberner (1983) Schönberner D., 1983, ApJ, 272, 708
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Staff et al. (2016) Staff J. E., De Marco O., Macdonald D., Galaviz P., Passy J.-C., Iaconi R., Low M.-M. M., 2016, MNRAS, 455, 3511
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Zou et al. (2020) Zou Y., et al., 2020, Monthly Notices of the Royal Astronomical Society, 497, 2855

##  Appendix A Stellar relaxation

In this Appendix we provide additional details of the stellar relaxation
procedures adopted in this paper and of the quantification of the hydrostatic
stability achieved.

To assess the fidelity of our stellar structures we plot, in Figure 14 the
density profiles of the four models, taken right after mapping, after the
damping (for the I-lo and I-hi models) or after the relax-o-matic for the R-hi
and M-hi models), time and after evolving with no damping.

The last part of the relaxation procedure, however, sees some particles moving
outward in both the R-hi and M-hi models. Also, both I-hi and R-hi have a
slightly higher central density than the cored mesa models.

To gauge the stability of the 3D hydrostatic equilibrium of the stellar
structures we checked the distribution of the SPH particle velocities at the
time after relaxation, when each profile was used as input to the binary
simulation. From there it is clear that only 1 per cent of all particles are
at velocities above 7.5 km s-1 for I-lo, R-hi and M-hi, while for I-hi it is
0.50 km s-1. While this is a reasonable number, and the best we can obtain
with the current method, we should be cognisant that at 7.5 km s-1 a particle
would cross the stellar diameter in 1.5 years, a time that is comparable with
the simulation time.

Figure 14: Density profiles of the four stellar models at various phases of
setup. The first column (blue) is for I-lo, the second (red) for I-hi, the
third (green) for R-hi and the last column (magenta) for M-hi. The coloured
profiles are the cored mesa profiles used in the initial mapping. The vertical
dashed lines show the softening length, $r_{\rm soft}$, of the core point mass
particles. The top row is right after mapping, the middle row is after damping
(for I-lo and I-hi) or relax-o-matic (for R-hi and M-hi) and the last row is
at the end of the 5 dynamical times ($\approx 0.56$ years).

##  Appendix B Conservation properties of our simulations

Figures 15 and 16 show the energy and angular momentum conservation properties
of the simulations. The conservation of energy and angular momentum are
excellent, with energy conserved to 0.05 per cent, and angular momentum
conserved to 0.1 per cent in I-hi and conservation in the other simulations
being similarly good. The energy and angular momentum exchanges between the
different components of our systems are typical for CE simulations (e.g.,
Iaconi et al., 2017; Sandquist et al., 1998). The simulation M-hi unbinds
almost the entire envelope, justifying the observed differences between the
left and right panels of Figures 15 and 16.

Figure 15: Energy evolution for the I-hi (left) and M-hi (right) simulations.
The plots show total energy (Tot energy), kinetic energy (Kin energy),
potential energy (Pot energy), thermal energy (Therm energy), orbital energy
of the point mass particles (Cores pot), potential energy of the gas without
point mass particles (Gas pot) and potential energy of the gas with the point
mass particles (Env pot).

Figure 16: Evolution of the $z$ component of the angular momentum for the I-hi
(left) and M-hi (right). In both cases, the total angular momentum, $J_{\rm
tot}$, is conserved. $J_{\rm orb}$, $J_{\rm b}$ and $J_{\rm unb}$ are the
orbital, bound and unbound (using mechanic criterion for all curves) angular
momenta, respectively.

##  Appendix C Sensitivity of the in-spiral to initial stellar structure

To test whether the speed of in-spiral on-set, in-spiral shape and final
separation were due to the fact that the initial stellar model in M-hi was
somewhat expanding, we calculated an additional model but where the initial
stellar model was taken at a later time, 2.5 years after application of the
relax-o-matic relaxation technique. At this time, the stellar radius had
expanded and then contracted, to reach a stable size of 285 R⊙.

In Figure 17 we see that the orbital separation over time for M-hi and for the
additional simulation are extremely similar. M-hi takes ever so sightly longer
to inspiral, but overall the two simulations have a similar in-spiral shape
and final separation. This demonstrates that the expansion of the single star
used for the M-hi simulation ($R=260$ R⊙, but expanding to 300 R⊙in $\sim$1
year; Figure 3) is not inducing a systematically different behaviour compared
to a more stable, albeit slightly larger (285 R⊙) stellar structure.

Figure 17: A comparison of M-hi (red) and an identical simulation started with
a stellar model that was evolved in isolation in the phantom computational
domain for an additional period of time, such that it reached a more stable,
albeit larger radius of 285 R⊙(black).

##  Appendix D Extrapolation of the final separation

In order to obtain a value of the final orbital separation that would more
closely resemble that which we would obtain by running the simulation further,
we extrapolated the apastron and periastron orbital distances observed at the
end of the simulations. To do so, we used the curve_fit function of the scipy
module in Python 3 (Virtanen et al., 2020) to carry out a non-linear fit. For
the fitting curve, we used an exponential function, $s(t)$:

$s(t)=\gamma e^{-\kappa t}+\omega,$ (6)

where $\gamma$, $\kappa$ and $\omega$ are the fitting parameters.
Particularly, $\kappa$ and $\omega$ are the decay factor and the extrapolated
final separation. We decided to fit perastron and apastron separations
separately. For each fit, we select the last $n$ data points, such that $n$ is
the number of points that minimizes the element $\sigma_{\kappa,\kappa}$ in
the covariance matrix $C(\gamma,\kappa,\omega)$:

$C(\gamma,\kappa,\omega)=\begin{bmatrix}\sigma_{\gamma,\gamma}&\sigma_{\gamma,\kappa}&\sigma_{\gamma,\omega}\\\
\sigma_{\kappa,\gamma}&\sigma_{\kappa,\kappa}&\sigma_{\kappa,\omega}\\\
\sigma_{\omega,\gamma}&\sigma_{\omega,\kappa}&\sigma_{\omega,\omega}\\\
\end{bmatrix},$ (7)

where the diagonal elements are the variance values for each of the three free
parameters, compared to the data points. The result of the fitting for
simulation I-hi is shown in Figure 18.

Figure 18: Extrapolation of the orbital separation for simulation I-hi (orange
curve). The green and blue solid curves are fits to the apastron and
periastron separations, respectively, using the method described in appendix
D. These curves have two different asymptotic values ($\omega$ in Equation 6),
indicated by the dashed lines.