How negative feedback and the ambient environment limit the influence of recombination in common envelope evolution

We perform 3D hydrodynamical simulations to study recombination and ionization during the common envelope (CE) phase of binary evolution, and develop techniques to track the ionic transitions in time and space. We simulate the interaction of a $2\,M_\odot$ red giant branch primary and a $1\,M_\odot$ companion modeled as a particle. We compare a run employing a tabulated equation of state (EOS) that accounts for ionization and recombination, with a run employing an ideal gas EOS. During the first half of the simulations, $\sim15$ per cent more mass is unbound in the tabulated EOS run due to the release of recombination energy, but by simulation end the difference has become negligible. We explain this as being a consequence of (i) the tabulated EOS run experiences a shallower inspiral and hence smaller orbital energy release at late times because recombination energy release expands the envelope and reduces drag, and (ii) collision and mixing between expanding envelope gas, ejecta and circumstellar ambient gas assists in unbinding the envelope, but does so less efficiently in the tabulated EOS run where some of the energy transferred to bound envelope gas is used for ionization. The rate of mass unbinding is approximately constant in the last half of the simulations and the orbital separation steadily decreases at late times. A simple linear extrapolation predicts a CE phase duration of $\sim2\,\mathrm{yr}$, after which the envelope would be unbound.


INTRODUCTION
In the common envelope (CE) scenario, a star engulfs a typically much smaller companion and the companion and core of the primary star spiral in together until either the envelope is ejected or the core and companion merge (Paczynski 1976).The CE phase is short-lived and hence difficult to observe.However, understanding CE evolution (CEE) is crucial for understanding phenomena as wide-ranging as neutron star-neutron star mergers, type Ia supernovae and planetary nebulae (for recent reviews see Ivanova et al. 2020 andRoepke &De Marco 2022).While theoretical work on CEE has come a long way, simulating CEE up to a realistic end state (e.g.complete envelope ejection and stabilization of the orbit) has still not been achieved.The overarching reason is that simulating such a large dynamic range of spatial, temporal and density scales remains challenging.
Nevertheless, such calculations have been useful in constraining ★ lchamandy@niser.ac.in † jcarrol5@ur.rochester.edu‡ blackman@pas.rochester.edu§ afrank@pas.rochester.edu¶ yt2cr@virginia.edu∥ baowei.liu@rochester.edu★★ yzou5@ur.rochester.edu† † nordhaus@astro.rit.edu the effects of various physical processes during CEE.One such process is the recombination of ions and electrons as the envelope expands and cools, which releases recombination energy.As has become somewhat standard in the CE literature, we define recombination energy as the latent energy contained by a plasma that would be released upon recombination.
While the recombination energy content of the initial envelope is generally substantial and thus could potentially play an important role in envelope unbinding, its efficacy remains unclear.This is partly because some released recombination energy radiates away (Soker & Harpaz 2003;Ivanova et al. 2013b;Sabach et al. 2017;Ivanova 2018;Grichener et al. 2018;Soker et al. 2018;Reichardt et al. 2020;Lau et al. 2022a).How much is radiated is not addressed in the present work.Rather, we ask: what difference does recombination energy make to CEE, assuming (optimistically) that released recombination energy is thermalized locally?In particular: how much extra envelope mass does recombination unbind at any given time, where in the envelope does this extra unbinding happen, and what are the relative contributions of the various ionic species?Even in the optimistic case that assumes local absorption and thermalization, some of the recombination energy would be released into already unbound gas, and play no further role in envelope unbinding.Therefore, there is an efficiency associated with the transfer of released recombination energy to binding energy of the remaining envelope.
Sometimes authors count gas as unbound if its total energy density, including that of the recombination energy, is positive.This is misleading because it greatly overestimates the role of recombination energy in envelope unbinding.First, the recombination energy is still latent.Second, upon release, it might contribute little to unbinding due to the inefficiencies mentioned above.A more justifiable approach is to include only kinetic (bulk and perhaps thermal) and potential energy terms in assessing whether material is unbound, and compare the unbound mass in a simulation with a realistic tabulated EOS that includes recombination energy with a simulation that does not.If cooling is neglected, it is appropriate to use a  = 5/3 ideal gas EOS for the latter simulation.When such a comparison has been made (Ohlmann 2016;Reichardt et al. 2020;Sand et al. 2020;Lau et al. 2022a;González-Bolívar et al. 2022), it has been found that the unbound envelope mass is significantly higher in the tabulated EOS run.Moreover, the fractional difference in the unbound mass curves between the runs usually increases with time, indicating that the effect becomes more important at late times.
The details of how and how much recombination energy can really assist envelope unbinding are still emerging.Reichardt et al. (2020) included a spatial analysis where they mapped out the ionization fractions of various species using the Saha equation.This allowed them to identify recombination fronts and analyze their evolution.They then compared this data with data of where unbinding is happening in both  = 5/3 and tabulated EOS runs, allowing them to infer how recombination affects unbinding.In this work we in addition employ gas tracers to track envelope gas of a given ionization species at  = 0. Comparing the present ionization state calculated using the Saha equation with the original ionization state recorded by the tracers enables us to map out the net change in the ionization state of gas in the simulation, as a function of both space and time.This method also makes it possible to directly compute the net energy absorbed by ionization or released by recombination for some subset of the gas -for example, gas which is presently bound.
We also find that we need to assess the role played by the ambient medium in our simulations, i.e. extra gas that fills the simulation box.We thus compare two  = 5/3 simulations, one with a dense ambient medium and one with a low-density ambient medium.This helps to place bounds on effects that would stem from the presence of a circumbinary torus, formed in nature by mass loss that precedes the CE interaction.
The structure of the paper is as follows.In Section 2, we explain the methods used in the setup, running and analysis of the simulations, and describe the physical model, including the initial conditions and the parameters of the runs performed.Results are then presented in Section 3. In Section 3.3 we present the evolution of the unbound mass with time, in Section 3.4 we analyze the transfer of energy between the various components, in Section 3.5 we explore the ionic evolution and the release of recombination energy, and in Section 3.6 we discuss the role of the ambient medium.We summarize our results and conclude in Section 4.

METHODS AND MODELS
We first clarify the terminology we use, which might otherwise cause confusion."Envelope" or "envelope gas" refers to gas that was originally part of the stellar envelope, and "unbound envelope gas" refers to such gas which has become unbound (and hence is no longer technically part of the envelope).Envelope gas is separate from ambient gas, which is extra gas in the simulation not included in the original spherically symmetric stellar model used in our initial condition.The gravitating point particles are referred to as core particles, separate from the gas (even though they can represent gaseous astrophysical bodies).
Our study involves a binary system consisting of a 1.96 M ⊙ , 48.1 R ⊙ red giant branch (RGB) primary star with a 0.366 M ⊙ core and a 0.978 M ⊙ secondary representing a main sequence star, white dwarf, or even a very low-mass neutron star.The primary core and secondary are modeled as gravitation-only point particles.A spline function is used to model the potential around the point particles, which becomes Newtonian outside the spline softening sphere, with radius  soft .The initial density and pressure profiles of the primary are mapped to our 3D grid from a 1D Modules for Experiments in Stellar Astrophysics (MESA) snapshot (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2019)).The snapshot is taken from a MESA release 12778 simulation.It matches almost exactly the snapshot used for our previous RGB simulations using release 8845 (Chamandy et al. 2018), except that here we increased the radial resolution by a factor ∼ 20 to smooth the profile.We utilize the adaptive mesh refinement (AMR) code AstroBEAR (Cunningham et al. 2009;Carroll-Nellenback et al. 2013), and employ an HLLC Riemann solver.Simulations are carried out in a centre of momentum frame with the grid centred on the initial location of the centre of the primary star.The simulation box of side length  box = 1150 R ⊙ is discretized into 512 3 AMR level zero cells, corresponding to a base resolution of  0 ≈ 2.25 R ⊙ .For further numerical method details, see Chamandy et al. (2018Chamandy et al. ( , 2019aChamandy et al. ( ,b, 2020)).
Initially, the envelope and some of the ambient medium surrounding the RGB star are resolved at AMR level 4, or  4 ≈ 0.140 R ⊙ , and this refinement zone gradually reduces as the particle separation  decreases.However, unlike in our previous RGB simulations, AMR level 5, with  5 ≈ 0.070 R ⊙ , was added around the point particles out to slightly farther than the softening sphere.This extra level of refinement helps to conserve energy and to avoid artificial reduction of the central density and pressure during the simulation.Buffer zones with 16 cells were included to smoothly transition between AMR levels.At  = 25.23 d, the softening radius around the particles was halved to ≈ 1.2 R ⊙ and a sixth level of refinement with  6 ≈ 0.035 R ⊙ was added around each particle.At  = 50.46d the softening radius was again halved to ≈ 0.6 R ⊙ , and AMR level 7 was added, with  7 ≈ 0.018 R ⊙ .The enhanced refinement zones are approximately spherical in shape and contain uniform resolution.

Initial orbit
The binary is initialized in a circular orbit at a separation  i = 49 R ⊙ , slightly outside the stellar surface ( 1 = 48.1 R ⊙ ).The primary has no initial rotation.The initial conditions are almost identical to those in our previous works, and similar to those of Ohlmann et al. (2016) (see also Ohlmann 2016) and Prust & Chang (2019), which allows for quantitative comparison (see Chamandy et al. 2019a and App.A).

Ambient medium
The pressure scale height near the core and at the stellar surface of the MESA profile are too small for our 3D simulation to resolve.Therefore, we cut out the RGB core and replaced it with a spline-softened gravitating point particle, with softening length  soft = 2.41 R ⊙ , equal to the cut radius, along with a core density and pressure profile obtained by solving a modified Lane-Emden equation (Ohlmann et al. 2017).This solution also incorporates an iteration over the particle mass to ensure that the mass below the cut radius is preserved (Chamandy et al. 2018). 1 The softening radius of the secondary is the same as that of the RGB core particle.
The ambient medium has uniform density  amb = 1.0 × 10 −9 g cm −3 and uniform pressure  amb = 1.0 × 10 5 dyn cm −2 ; the ambient pressure is added everywhere on the grid to ensure a smooth transition in pressure at the stellar surface, while the density undergoes a sharp transition at the stellar surface.The value of  amb is chosen such that the RGB pressure profile effectively truncates (but smoothly) just inside the outer radius of the star to avoid the very small pressure scale height at that location in the MESA profile.The value of  amb is about 7 times smaller than the density at the outer radius of the star ( 1 ); smaller values are possible but would result in a higher ambient temperature and prohibitively small timesteps.We find that we require at least eight resolution elements per scale height to adequately resolve the initial stellar profile.This number was determined by studying the smoothness and stability of both the core and surface, where the scale heights are smallest, during the first ∼ 1 d of the simulations. 2 To ascertain the effect of the ambient medium on the results, we perform one run for which we drastically reduce the density and pressure of the ambient gas.We reduce the ambient medium near the beginning of an ideal gas run, at a time chosen as a compromise between reducing the ambient as early as possible and retaining the stability of the stellar surface.Also, the later the restart from the fiducial, the less the computational cost. 3This ambient reduction procedure is made possible by the tracer (Sec.2.3) assigned to the ambient gas during the initial setup at  = 0. Thus, between  1 = 19.9d and  2 = 24.5 d, the ambient density and pressure are reduced by three orders of magnitude according to the expressions and where  ≡ ( −  1 )/,  = 10 −3 ,  = 4 and  = 2.3 d.The parameter  determines how gradual the transition is.The code was stable for test runs with  an order of magnitude shorter, but we used the above value to be conservative.

Tracers
In the initial setup, tracers are added to track gas that is initially in the envelope (0 <  i <  1 ) or ambient ( i >  1 ).Tracking the ambient gas allows us to exclude it in post-processing (Zou et al. 2022) or to remove it during the simulation (Sec.2.2).By employing tracers, the relative mass of each component (e.g.envelope or ambient gas) in any given AMR cell can be accurately determined.Thus, the method remains effective if there is mixing between various components.Separate tracers are also added to track the initial hydrogen and helium ionization states of the gas.The number density profiles of the key ionic species, as computed from the Saha equation, are shown in Fig. B1.Then, by comparing the ionization state from the Saha equation at time  with the original ionization state from MESA, we deduce how the ionization state evolves in the simulation.

Equation of State
Our tabulated EOS is essentially the "DT2" part of the MESA tabulated EOS (see Paxton et al. 2019), which is based on the OPAL EOS (Rogers & Nayfonov 2002) and the SCVH EOS (Saumon et al. 1995).For the parameter space of the simulations, OPAL is much more relevant than SCVH.These tables have been adapted for use in AstroBEAR; in some cases, this requires inverting tables to change the independent variables.The tables cover the full parameter space encountered by the simulations presented. 4 For our 3D simulation, we modified the tabulated EOS by subtracting the radiation specific energy  rad  4 /, where  is temperature,  is gas density and  rad is the radiation constant.In the original MESA RGB profile, the ratio of this component to the local gas thermal specific energy 3 B /2 H , where  B is the Boltzmann constant and  is the mean molecular mass units of the hydrogen mass  H , reaches a maximum of 20% at  ≈ 0.6 R ⊙ and is < 7% outside  soft = 2.41 R ⊙ .The ratio of the total radiation energy to that of the total thermal energy over the entire star is < 0.1%.If included, the radiation energy would lead to a high internal energy density in the high-temperature ambient medium, which could artificially help unbind the envelope through mixing.The choice to remove radiation energy is also consistent with the fact that radiation pressure is not included in the tabulated EOS, although it is included separately in the stellar models computed by MESA. 5  That said, we noticed during post-processing that the subtraction of the radiation energy was imperfect, due to a slightly different value of the radiation constant  rad used for the subtraction compared to that implied by the tables.This resulted in 1.5 × 10 −4 of the radiation energy remaining in the simulation.This leftover radiation energy is rather negligible in the envelope, but makes up about two thirds of the energy of the ambient medium at  = 0; the gas thermal energy makes up about 90 per cent of the remaining one third, and the rest is recombination energy.Thus, the energy density of the ambient medium is a few times higher in the tabulated EOS run as compared to the ideal gas run.We account for this leftover radiation energy in the analysis below.

Runs
Our runs are summarized in Table 1.Model A employs an ideal gas EOS, with  = 5/3.This is an updated version of the fiducial RGB CE model from our previous works (Chamandy et al. 2018(Chamandy et al. , 2019a(Chamandy et al. ,b, 2020)), with higher resolution and with initial ambient density reduced by a factor of about seven.Model B makes use of the MESA EOS with radiation energy removed, as explained in Sec.2.4, but is otherwise the same as Model A. It assumes that recombination energy is thermalized locally when it gets released.Model C is like Model B except that the internal energy is replaced by the gas thermal energy alone; since both recombination and radiation energy are removed, the results are expected to be closer to those of Model A than Model B. Model D is like Model A except that the ambient density and pressure are reduced by a factor of 1000 early in the simulation (Sec.2.2).Models A and B were run for 115.7 d, Model C for 50.5 d, and Model D for 142.4 d. 4 To make the tables rectangular in log 10 of the independent variables  and ,  and specific energy , or  and , we extrapolate them using an ideal gas law with mean atomic weight of metals set equal to the value 17.4 amu.This is the value obtained from the MESA model.However, the extrapolated regions correspond to regions of the parameter space that do not arise in the simulations. 5When preparing the simulation initial condition, we chose to make the gas pressure equal to the total pressure (gas plus radiation) in the MESA 1D RGB profile.

RESULTS
This presentation of results is organized as follows.We first present 2D snapshots that illustrate the evolution of certain key quantities in Sec.3.1.We then discuss the orbital evolution of the core particles in Sec.3.2 and analyze the unbinding of the envelope with time in Sec.3.3.In Sec.3.4 we attempt to shed light on the envelope unbinding process and role of recombination energy by focussing on the energy transfer in the bound envelope gas.We then study the ionic transitions and their energetics in detail in Sec.3.5.In Sec.3.6 we try to determine what role the ambient medium is playing in the simulations.

Overall evolution and spatial dependence
In this section we focus on the tabulated EOS run, Model B. Snapshots are shown in Figs. 1 and 2. Slices are through the orbital plane and axes are those of the lab frame.From top to bottom, the first row shows the gas density , the second row the temperature , the third row the density of each helium ionization species (He I in green, He II in orange and He III in blue), the fourth row the density of helium ion tracer gas, indicating initial ionization state (using the same colour scheme), the fifth row a measure of the unbound mass, discussed below, and the sixth row is the same as the fifth but for

Density
The top row of Figs. 1 and 2 shows the mass density in g cm −3 for Model B. The primary core particle (secondary) softening sphere is visible in the first four columns as a purple (red) circle.Since the softening radius decreases by a factor of four during the simulation, the softening spheres are too small to be plotted in the final two snapshots.We do not include colour-plots of density for the other run, Model A, because the results are very similar to Model B, but density contours for Model A are shown in the bottom row.

Temperature
The second row shows the gas temperature for Model B, with density contours overplotted.At first (Fig. 1) the local temperature does not undergo much change as the star deforms.Then the outer envelope shows a modest amount of cooling (Fig. 2 left and middle columns), and eventually it heats up again owing to mixing with ambient gas (Fig. 2 right column).

Recombination and ionization
The third row of Figs. 1 and 2 shows the helium ionization state of the gas, as computed from the Saha equation.This is done by plotting the mass density of each ionic species: He III (blue), He II (orange) and He I (green).We focus on helium rather than hydrogen because, as we will see, helium releases more recombination energy during the simulation.For clarity, we plot only the ionic species whose density is highest at a given position.In other words, for ease of presentation, we do not attempt to plot overlapping colours.However, the calculations below include contributions from every hydrogen and helium ionic species at every location. 6o get a sense of where and when ionization and recombination are occurring, we need to understand how the ionization state of each fluid element changes over time.The fourth row shows the density of tracer gas for each of the helium ion species (Sec.2.3).These plots tell us where envelope gas that was originally dominated by one or the other helium ionic species at  = 0 is located at the present time .Visually comparing the third and fourth rows allows one to broadly understand where ionization and recombination have taken place. 7s an example, let us first consider the snapshot at  = 11.6 d (Fig. 1, right column), which is just before the first periastron passage time of  = 13.0 d.Clearly, much of the He III tracer gas (4th row, blue) has recombined to become He II (3rd row, orange) or He I (3rd row, green).Additionally, much of the He II tracer (4th row, orange) has recombined to become He I (3rd row, green).However, the bulk of the gas remains He III (3rd row, blue).Interestingly, ionization has also occurred in some places.In particular, it appears that some of the He II tracer gas gets dragged down by the secondary to hotter regions, where it gets ionized.This is visible in the thin wake of dense (first row) and hot (second row) material attached to the secondary (red circle) that is orange in the fourth row (originally dominated by He II) but blue in the third row (now dominated by He III).To the right of the primary core particle, there is another stream of material within which ionization of He II has occurred.
Ionization also happens in envelope gas that is in direct contact with the ambient medium.At the early times shown in Fig. 1, this seems to be caused by shock heating at the envelope-ambient interface.This can be seen along the outer edge of the structure, where blue He III has replaced orange He II and green He I.However, this gas is very low-density, so this ionization does not have a large influence on the energy budget at this stage.At a later stage (Fig. 2), mixing with the ambient medium heats and reionizes gas in the outer envelope, and we revisit this in Sec.3.6.

Envelope unbinding
The fifth row of Figs. 1 and 2 shows the differentiation of the envelope gas into bound (blue) and unbound (red) components, along with velocity vectors in the orbital plane and density contours, for Model B. The sixth row shows the same snapshots, but now for the comparable ideal gas run, Model A. Gas is defined to be unbound if (3) Here, E blk,env is the bulk kinetic energy density of envelope gas, E thm,env is the thermal energy density of envelope gas, 2E pot,env−1 is twice the potential energy density of envelope gas associated with the RGB core particle, 2E pot,env−2 is twice the potential energy density of envelope gas associated with the secondary, and E env−gas is the potential energy density due to self-gravity. 8The quantity plotted in the fifth row of Figs. 1 and 2 is where To begin with, comparing the panels in the fifth (Model B, tabulated EOS) and sixth (Model A, ideal gas EOS) rows of Figs. 1 and  2 shows that the envelope unbinding process is extremely similar in the two runs, but that Model B has a slightly faster expansion.In addition, the unbound mass is slightly higher in Model B; this is most evident when comparing the panels in the left column of Fig. 2 8 The factors of 2 are included to account for the half of the potential energy nominally located inside the core particles.The self-gravity term is labeled "env-gas" because it uses the gravitational potential from all the gas (en-velope+ambient) to determine whether envelope gas is bound or unbound.
Omitting the ambient from the potential would reduce its magnitude by ≲ 1%, and this small difference would be comparable for Models A, B and C. No definition of unbound predicts with accuracy whether fluid elements will ultimately become unbound, because the outcome will depend on other factors (e.g.Ivanova et al. 2013a).The definition used here is fairly conservative.We assess the sensitivity to the definition of "unbound" in Sec.3.3; see also Chamandy et al. (2019a).
( = 23.1 d).The faster expansion and greater unbinding are consistent with the expectation that the energy released into the gas owing to recombination assists envelope unbinding.However, the similarity between the runs shows that the release of orbital energy, rather than recombination energy, is the main driver of mass unbinding at this stage.
How much of the released recombination energy helps to unbind the envelope?Recombination in gas which is already unbound will likely not assist unbinding.By comparing the third, fourth and fifth rows of Figs. 1 and 2, we can infer how much recombination happens in bound versus unbound envelope gas.Recombination is significant in both, which confirms that only part of the released recombination energy is utilized to unbind the envelope.We address this further in Sec.3.5.The separation curve of Model C is quite similar to those of Models A and B, but shows less correspondence with Model A than does Model B. Recall that in Model C we replaced the internal energy of the EOS tables with the thermal energy alone, akin to Model A, so the distinct evolutions between between Model C and Model A must be due to other model differences.

Orbital Evolution of the core particles
The separation curve for the reduced ambient run Model D agrees very closely with that of Model A. This verifies that the orbital evolution is insensitive to the density and total mass of ambient gas.
Finally, note that the mean separation continues to decrease monotonically at the end of the simulations, which shows that the orbital inspiral has not ended.This is consistent with the fact that envelope gas is still present around the core particles.We expand on this point in Sec.3.7.

Evolution of the unbound mass
Fig. 4, shows the unbound mass as a function of time (see also Sec. 3.1.4).To aid interpretation and facilitate comparison with the literature, we present the unbound mass using different definitions of "unbound."Our fiducial definition is given by inequality (3), and the corresponding unbound mass is represented by solid lines in Fig. 4. Model A is shown in red, Model B in blue, Model C in green (or dotted black), and Model D in yellow.Spatial information to go with the blue and red solid lines are provided in the fifth and sixth rows, respectively, of Figs. 1 and 2.
The dashed lines in Fig. 4 do not include the term E thm .This is thus a more conservative definition of "unbound."The dotted lines show a more liberal definition, where E thm,env is replaced by E int,env .As discussed in Sec. 1 however, the latter definition is unrealistic because it effectively assumes that all of the recombination energy will be released into gas that is still bound, or into marginally unbound gas that would otherwise rebind.In reality, the process will not be perfectly efficient because some of this energy will be released in gas that is already unbound, or radiate away.This is also why orbital energy cannot be utilized with 100% efficiency to unbind the envelope (Chamandy et al. 2019a).An efficiency  rec < 1 can be assigned analogous to  CE < 1 for orbital energy (e.g.Zorotovic et al. 2010).In addition, recombination energy is latent not kinetic, so it is misleading to include it unless released as kinetic energy.
The solid and dashed curves of Fig. 4 tell a similar story.Model B (tabulated EOS) unbinds ∼ (10-20)% more mass than Model A (ideal gas EOS) starting at the first periastron passage until  ≈ 40 d.This is expected based on previous work (Sec.1), and because the envelope acquires latent energy released by recombination, some marginally bound envelope material unbinds and some marginally unbound material is prevented from rebinding.However, at  ≈ 40 d the difference in unbound mass between the runs declines, nearly vanishing by  ≈ 70 d and later.We will later address why this difference in unbound mass vanishes at late times even though recombination occurs in Model B but not Model A.
Next consider the solid and dashed green curves, and dotted black curves, showing Model C, where the internal energy was replaced with thermal energy in the tabulated EOS.As expected, the curve that includes the internal energy (dotted) precisely traces the curve that includes thermal energy only, giving confidence that the EOS was modified correctly.Also, the solid and dashed curves closely trace their Model A counterparts.The general agreement between Models A and C is expected, and suggests that the differences in the unbound mass between Models A and B are due to release of recombination energy.That said, the agreement is not perfect, and Model C unbinds slightly more mass than Model A.

Envelope energization, unbinding and expansion
We now turn to Fig. 5, which shows the evolution of the volumeintegrated energy components of the bound envelope gas, where the fiducial definition of unbound [equation (3) and solid lines in Fig. 4] has been used. 9The subscript 'bnd' refers to that fraction 9 As the bound envelope mass is not the same between runs (Sec.3.3), we also considered the specific energy, but that plot is qualitatively similar to Fig. 5 and leads to the same conclusions.
of the energy contained in gas that is bound at time .The opaque curves represent Model B and the semi-transparent curves represent Model A in Fig. 5.The particle energy  1−2 , in green, is the sum of the kinetic energies of the particles and their mutual gravitational potential energy. 10 The blue curves show the gas thermal energy, the red curves show the bulk kinetic energy, and the purple curves show the potential energy, including the gas self-gravity as well as the potential energy associated with the mutual attraction between the core particles and the bound gas.The orange curves show the sum of blue, red and purple curves, and the black curves are the sum of green and orange curves. 11 Let us compare the two runs at  ≈ 50 d, when the unbound mass is still significantly higher in Model B than in Model A (Fig. 4).The purple lines in Fig. 5 show that | pot,env,bnd | is smaller in Model B at this time, so the envelope has expanded more than in Model A. This extra expansion was discussed in Sec.3.1.4and is likely caused by the extra energy supplied by recombination.The total energy  env,bnd is correspondingly higher (i.e. less negative) in Model B. And although Model B incurs a slightly deeper inspiral up to this time as implied by the gap between the green lines, this deeper inspiral can only partially account for the wider energy gap between the orange lines.
However, the fact remains that the difference in the potential energy between the runs (purple) is larger than that of the total gas energy (orange) -why?The bulk kinetic energy, shown in red, is almost identical between the runs, but the thermal energy content is higher in the ideal gas run (blue), so this explains why the gap between the orange curves is smaller than that between the purple curves.Although recombination converts latent recombination energy into thermal energy, the extra heat from recombination is evidently used to expand the envelope, causing it to cool, in turn resulting in less overall thermal energy content compared to Model A.
Next we compare the energy terms at the end of the simulations at  ≈ 115 d.The total energy of the bound system of envelope and core particles is clearly higher in Model B (see black lines in Fig. 5).However, the energy of the gas is almost equal for the two runs (orange lines).This is because the inspiral is deeper in Model A, resulting in smaller (i.e. more negative) orbital energy of the system of core particles (green lines).The extra release of orbital energy in Model A has evidently led to more expansion and cooling of the envelope than in Model B.
The overall picture is that energy released by recombination leads to extra expansion and cooling of the envelope.This reduces the drag force on the core particles, their inspiral, and their transfer of orbital energy from the cores to the envelope.This picture implies the existence of a self-regulating mechanism that limits the effectiveness of the recombination energy for energizing and unbinding the envelope.

Spatially Integrated Recombination and Ionization
In Sec.3.1.3we saw that much recombination and ionization occur in Model B as the envelope expands, particularly at late times.To better understand these processes and their effect on unbinding, we now explore the volume-averaged evolution of the recombination 10 In addition to the kinetic energy associated with the quasi-Keplerian motion,  1−2 also includes the kinetic energy associated with the centre of mass motion.This contribution lessens as the system evolves and by  = 40 d it is negligible (Chamandy et al. 2019a). 11The black curves do not include the energy of unbound envelope gas or ambient gas and so need not remain constant.The discontinuities at  ≈ 25 d and  ≈ 50 d are due to the instantaneous reduction in the softening radius (Sec.2) and the associated extra negative potential energy.Although unphysical, it does not affect our results.energy associated with the various ionic transitions, where transitions are defined to include both single and double recombinations or ionizations.

Helium
Fig. 6 shows the cumulative energy release by the various ionic transitions of helium up to time , for all envelope gas (black) and bound envelope gas (blue).Recombination is shown in the top row and ionization in the bottom row.Positive values indicate net release of recombination energy and negative values indicate net gain of recombination energy due to ionization.Fig. 6 shows that the He III → He II and He III → He I transitions are the most important, with the former releasing the most recombination energy by the end of the simulation, and mostly from bound gas.Net ionization of tracer gas during the simulation occurs over time, but mostly in unbound gas, and much less commonly than recombination.Reionization of recombined tracer gas is however important, as discussed below.For reference, Fig. C1 in App.C shows the mass of each ionic species for each of the three helium ion tracers, and App.C also provides details of how the energy release is calculated.
The top right panel of Fig. 6 shows a peak in the energy released by He III → He I at  ≈ 30 d, and then a decline until  ≈ 80 d, followed by a local peak at  ≈ 100 d.At the first peak, about 40% of the released energy comes from unbound envelope gas, so much of the released energy is "wasted" in energizing already unbound gas.The decline after the first peak is caused by reionization of He I into He II and, to a lesser extent, He III.This is evident from the right column of Fig. C1.In the middle panel of this column, we see that between   5), and is shown by dashed lines.The green line shows the difference in the energy injected by the particles in the two runs, with the ideal gas run hosting a deeper inspiral and thus greater energy release.By the end of the simulation, the absence of recombination energy release in the ideal gas run is fully or partially compensated by the deeper inspiral.also shows that much energy is used to ionize He II tracer gas into He III.Ionization absorbs bulk kinetic, thermal and potential from the surroundings, hindering unbinding.However, comparing the black and blue lines reveals that most of the ionization happens in already unbound gas.
The net recombination energy release for helium is shown in the middle panel of Fig. 7.The peak occurs at  ≈ 25 d, when ∼ 2/3 of released energy has been released by bound gas, which can be used to unbind the envelope.

Hydrogen
The net energy release from the recombination and ionization of hydrogen is shown in the left panel of Fig. 7. Transitions between three ionic states (H 2 , H I and H II) are included in the calculation, but H 2 does not play a significant role.The energy released by hydrogen is small compared to that from helium, despite the greater hydrogen mass fraction.This is mainly because only a small fraction of the primarily H II envelope recombines into H I. Moreover, less than half of the hydrogen recombination energy release, at its peak at  ≈ 35 d, is from bound gas.Thus, hydrogen less efficiently releases recombination energy than helium up to this time because the former releases most of this energy in already unbound material.These results are qualitatively consistent with those of Lau et al. (2022a) and Lau et al. (2022b), who studied different binary systems from ours.Fig. 7, shows that most all of the H I formed from recombination of H II reionizes into H I, likely due to the interaction between envelope and ambient gas.

Net effect from hydrogen and helium
Fig. 7 shows the net recombination energy released by hydrogen and helium.The spatially integrated release continues up to  ≈ 25 d in both total and bound gas.Subsequently, He II, He I and H I reonizations decrease the net energy released.The net released recombination energy in the bound envelope gas plateaus higher than that of the total envelope gas at the end of the simulation, as some of the unbound helium was ionized (Fig. 6, bottom row).

Net release of recombination energy from all species
Fig. 8 compares the recombination energy released by the envelope gas according to the above analyses (solid black and blue lines, as in Fig. 7) with the same quantities obtained in a different way (dashed black and blue).To obtain the dashed black line, the total recombination energy was computed by differencing the total internal and thermal energy of the envelope gas, taking into account leftover radiation energy (Sec.2.4).Thus, where  rad,env =  rad,env  4 and  rad,env is the small difference in  rad of about 1.5 × 10 −4  rad discussed in Sec.2.4.The thermal energy density is given by where the mean molecular mass  is taken from the simulation, and calculated using the EOS tables.The dashed black line in Fig. 8 shows the released recombination energy, equal to negative the net change in  rec,env at time , denoted as − rec,env .The difference in the released orbital energy  1−2 between Models A and B is also plotted for comparison (green).Details of the calculation are provided in App. C. The solid and dashed black curves in Fig. 8 nearly agree, but with more recombination energy release obtained by the more direct method described just above.Metals likely account for some of the difference, as they have ionization energy of about 2 × 10 45 erg.Most of this would remain locked up in ions, so this cannot account for the full excess."Leftover" radiation energy (Sec.2.4) rises as gas is heated, and this energy must be taken from other forms of energy such as recombination energy.This energy absorption, − ∫   rad,env  4 with ∫  denoting volume integration over the simulation domain, is plotted as a dashed-dotted black line in Fig. 8.This extra energy component may also contributee to the discrepancy between solid and dashed black lines.The shape of the black dashed line is sensitive to the value of  rad most appropriate for the tabulated EOS, which is not precisely known.In any case, the small discrepancy between solid and dashed lines does not impact our main conclusions.
The quantity − rec,env,bnd is plotted as a blue dashed line in Fig. 8. Likewise − ∫   rad,env,bnd  4 is plotted as a blue dasheddotted line, but its contribution is negligible, which means that all of the leftover radiation energy is in unbound gas.Since it does not affect bound gas, the leftover radiation energy is likely inconsequential for envelope unbinding.

Comparison with orbital energy release
Comparing the green line with the blue solid and dashed lines in Fig. 8 shows that the release of recombination energy into bound envelope gas is compensated by extra orbital energy released in the ideal gas run, owing to deeper particle inspiral.As argued in Secs.3.3 and 3.4, this likely explains why recombination energy release in Model B does not result in more unbound mass than in Model A.   4), is a measure of binding energy density: red gas is unbound and blue gas is bound (see Sec. 3.1.4).In Model A, there is mixing with the ambient medium (the ambient gas is uncoloured), whereas in Model D the ambient medium is not present the slice because it has been completely displaced by envelope gas.

Role of the ambient medium
To assess whether the ambient medium influences envelope unbinding, we compare Models A and D, which are identical except for ambient gas removal in the latter.Fig. 9 shows the evolution of various bound and unbound envelope gas energy components.The subscript "box" means that only gas remaining in the simulation domain is counted.The total energy of envelope gas including that lost through the simulation boundary is shown in black, while that which remains is shown dashed.The difference becomes significant after  ≈ 50 d as only unbound gas (equation 3) is lost. 12etween  ≈ 25 d, soon after the ambient reduction, and  ≈ 40 d, the bulk kinetic energy in Model A reduces relative to that of Model D. This is caused by a reduction in the unbound component.There is a corresponding deficit in the gas energy, shown in orange, and in the total energy, shown in black.In Model A, this is explained by a transfer of energy from the outgoing ejecta to the ambient medium (Chamandy et al. 2019a).
Between  ≈ 40 d and  ≈ 60 d, the deficit in  blk,env,box between the runs reduces to zero, and the thermal energy becomes larger in Model A as compared to Model D. From the middle and right columns of Fig. 2, it can be seen that mixing occurs between ambient and envelope gas during this time, likely from the Rayleigh-Taylor instability (Chamandy et al. 2019a).This mixing accelerates and compresses envelope gas, raising thermal and bulk kinetic energy components relative to Model D.
After  ≈ 60 d,  blk,env,box is about the same in the two runs, but the envelope continues to gain thermal energy at the expense of potential energy in Model A relative to Model D. In Model D, ejected envelope gas can freely expand, but in Model A, the inertia and weight of the ambient stalls the expansion.The expanding envelope collides with previously ejected overlying envelope, resulting in more shock heating, and limiting outward flow.Density and E env snapshots at  = 92.6 d are presented in Fig. 10, with Model A and Model D in the left and right columns respectively.The structure of the bound, inner envelope is similar between the two runs, but that of the outer envelope is very different.

Effect of ambient on envelope unbinding and ionic state
The unbound mass evolution of Models A and D can be compared using the red and yellow curves of Fig. 4. Beginning at  ≈ 50 d, the unbound mass curve of Model A (solid red) exceeds that of Model D (solid yellow).There is no corresponding surge in red relative to yellow in the dashed lines, which omit thermal energy in the definition of unbound.Therefore, Model A unbinds more mass than Model D due to extra heating of envelope gas.This extra heating, discussed earlier, is visible in Fig. 9.
In Fig. 4, the yellow solid and dashed curves of the reduced ambient run, Model D, tell an interesting story.Until  ≈ 70 d, these curves nearly overlap and are flat after the first periastron passage.Thus Model D unbinding up to the first periastron passage at  ≈ 13 d is driven by the bulk kinetic energy and thermal energy is unimportant.Only when a dense ambient medium is present, as in Model A, does some of this bulk kinetic energy convert to thermal energy, helping to unbind the envelope (compare red solid, yellow dashed and red dashed lines).The solid yellow line increases again at  ≈ 70 d, while the dashed yellow line stays flat and then increases at  ≈ 85 d, though not as fast as the solid line.This indicates a second phase of unbinding begins at  ≈ 70 d kickstarted by thermalisation of core particle orbital energy.That the solid red and yellow lines have almost the same slope after  ≈ 70 d shows that the rate of mass unbinding is roughly independent of the ambient density, as long as the thermal energy is accounted for in the definition of unbound.
The additional thermal energy helps to unbind envelope gas in Model A between  ≈ 50 d and  ≈ 70 d (red solid line in Fig. 4), but not in Model B (blue solid line), even though the ambient is equally dense and high-pressured with even more energy owing to the leftover radiation energy (see Sec. 2.4).Moreover, the envelope expands slightly faster at early times in Model B owing to the release of recombination energy (Sec.3.1), suggesting a stronger interaction between envelope and ambient gas.Why does  env,unb remain flat during this time for Model B (tabulated EOS), but increase for Model A? In Model B, the extra energy transferred from ambient gas may be redirected into recombination energy (i.e.latent energy of reionized gas).This is supported by the analysis of Sec.3.5.The left panel of Fig. 7, shows that all of the recombined hydrogen is reionized between  ≈ 35 d and  ≈ 85 d.Similarly, from the top-right panels of Fig. C1 and Fig. 6, between  ≈ 30 d and  ≈ 80 d, most of the He I formed by recombination of He III tracer gas reionizes into He II or He III.But this reionization cannot account for all of the He II formed, which means there is still recombination of He III into He II taking place (compare Fig. C1 middle-right panel).But the release of this recombination energy happens almost exclusively in bound gas, likely deep inside the envelope.In short, the extra unbinding that happens between  ≈ 50 d and  ≈ 70 d in the ideal gas run (Model A) but not the tabulated EOS run (Model B) likely occurs because the latter transfers ambient energy to reionize the gas rather than heat it.In Model B, ambient energy is also likely being used to increase leftover radiation energy of unbound envelope gas, Ambient energy is also likely being used to increase leftover radiation energy of unbound envelope gas, as discussed in Sec.3.5.

Lack of Stalling and Extrapolation to Termination
Despite Models A, B and D having completed over 100 orbits, the orbit-averaged separation steadily decreases, as seen in the inset of Fig. 3 and the green lines in Figs. 5 and 9.There is no evidence for "stalling" of the inspiral, as is often reported in the CE simulation literature.From  ≈ 70 d onward, the envelope unbinds at almost the same rate of ≈ 0.91 M ⊙ yr −1 in all three simulations (slope of solid lines in Fig. 4).Linearly extrapolating this constant rate, we estimate that the envelope would be ejected (i.e.fully unbound) at  ≈ 1.7 yr (where  = 0 corresponds to our initial condition).This is about four times shorter than the estimate obtained using the same method for a similar simulation involving a more evolved (asymptotic giant branch; AGB) primary (Chamandy et al. 2020).Interestingly, this ratio of 1/4 happens to be equal to the ratio of the initial orbital periods of the RGB and AGB simulations.
Likewise, we can extrapolate the orbital inspiral down to the final separation  f , with the latter estimated using the CE energy formalism where the change in binding energy  bind equals   1 ( 1 −  1,c )/ 1 , with  ≈ 1.3 (Chamandy et al. 2019a).Let us focus on the longest run, Model D. From  = 115 d onward, the nearly constant rate of core particle orbital energy loss is  1−2 ≈ −3.9×10 44 erg d −1 , where  1−2 ≈ −  1,c  2 /2.Assuming that  1−2 remains constant, we can calculate the duration of the CE phase for different values of  CE .The energy formalism yields for the final separation  f = 0.3, 0.8, 1.5 or 2.6 R ⊙ for  CE = 0.1, 0.25, 0.5 or 1, respectively (Chamandy et al. 2019a).We then estimate 12.5, 4.6, 2.0 or 0.7 yr for the CE duration.Conversely, for a duration  ≈ 1.7 yr as estimated above, we obtain  CE ≈ 0.55.We can take the same approach using the orbital angular momentum of the particles in the particle centre-of-mass frame.Approximating the orbit as circular,  ,1−2 ≈ √︁   1,c  2  red , where  red =  1,c  2 /( 1,c +  2 ) is the reduced mass.From  = 115 d onward, the torque  ,1−2 is approximately −3.5 × 10 43 g cm 2 s −2 .Assuming that this rate of angular momentum transfer remains constant and that the CE terminates at  ≈ 1.7 yr, we obtain  CE ≈ 0.27.Thus, extrapolations of energy and angular momentum both yield comparable values of  CE < 1 which are broadly consistent with observation-based estimates (e.g.Scherbak & Fuller 2023).The same procedure can be applied to Models A and B.  1−2 and  ,1−2 are approximately constant from 100 d onward in Model A and 92 d onward in Model B. Extrapolating  1−2 yields  CE = 0.4 for both runs, but extrapolating  ,1−2 yields much smaller values.This mismatch suggests that the estimates are less meaningful for Models A and B than they are for Model D. Finally, we note that  CE ultimately depends on the details of physical processes not yet included in numerical simulations, such as convective and radiative transport (Wilson & Nordhaus 2019, 2020, 2022).
During the revision of this work, Valsan et al. (2023) reported an ideal gas simulation run using the moving mesh code MANGA with parameters very similar to our Model D. However, their simulation ran for much longer, about 13 yr.They found that the envelope was fully ejected after about 4 yr, and that 80% of the envelope mass was ejected after 1.1 yr.The estimate of ∼ 2 yr obtained from our naive extrapolation of the unbound mass thus aligns quite well with their results.However, Valsan et al. (2023) also reported that the separation of the central binary plateaus at about 5 R ⊙ after about 200 d, whereas in our Model D the mean separation is 3.1 R ⊙ at the end of the simulation ( = 142 d) and is still decreasing.This difference is unlikely to be caused by the slightly different initial conditions, and suggests possible limitations in the various methods, which needs to be explored.In Appendix A, we provide comparisons between various models in the literature that use parameters similar to the ones used in the present work.In any case, extending Model D up to envelope ejection would be interesting future work.

CONCLUSIONS
We performed 3D hydrodynamic simulations involving a 2 M ⊙ RGB primary and 1 M ⊙ point particle companion to study the role of ionization and recombination in envelope unbinding and orbital evolution.We also studied the role of the ambient gas.The runs were performed using the code AstroBEAR, and employed adaptive mesh refinement and spline softening for the point particles.We compared four simulations, as summarized in Tab. 1, to isolate the effects of ionization/recombination and ambient-envelope interaction.Our main conclusions are as follows: 1) Released recombination energy somewhat helps to unbind the envelope, with a 10-20 per cent increase in unbound mass compared to the ideal gas run during the first half of the simulation.
2) However, release of recombination energy also expands the envelope, reducing the drag force, and reducing the orbital energy release relative to the ideal gas case.Thus, we find a stabilizing feedback which limits the gain in unbound mass caused by the release of recombination energy.
3) Interaction between bound envelope gas, unbound ejecta and the ambient medium leads to transfer of ambient energy to the envelope.Some of this energy is used to reionize gas, which prevents it from being used to directly unbind the envelope.This is a second stabilizing feedback that limits the effectiveness of recombination energy for envelope unbinding when a dense circumstellar environment is present.
4) By employing tracers in the code, we separately tracked gas initially dominated by each ionic species (i.e.H I, H II or H 2 for hydrogen and He I, He II or He III for helium).Computing the ionization state using the Saha equation then enabled us to follow the net ionic transitions that had occured up until a given time.Given that tracers are computationally inexpensive, future work could, in addition, include tracers for metal species and even introduce new ionic tracers at every snapshot to better resolve the time evolution.5) We also tracked the various ionization transitions for bound and unbound gas separately, and found that at the time of peak cumulative release of recombination energy ( ≈ 25 d), ∼ 2/3 was released by envelope gas that is still bound, suggesting that this recombination energy is used efficiently for unbinding.
6) Recombination energy release was dominated by helium.Recombination energy release by hydrogen contributed at first, but hydrogen was later almost fully reionized.Moreover, at the time of maximum recombination energy release of hydrogen, more than half of this release happened in already unbound gas, so was not used as efficiently as helium recombination energy for envelope unbinding.Radiative losses, not included in our simulation, likely reduce this efficiency further.7) At late times (> 100 orbits) we find a steady decrease in the period-averaged orbital separation, not a "stalled" inspiral as commonly reported in the literature.The approximately constant rates of mass unbinding, core particle orbital energy release and core particle orbital angular momentum release at late times motivate a simple linear extrapolation to the termination of the CE phase.This predicts a CE phase duration of about 2 yr, and a best-guess 0.2 ⩽  CE ⩽ 0.6.
The second stabilizing effect summarized in item (3) depends on the ambient gas.The uniform ambient medium of density 10 −9 g cm −3 in our Models A, B and C was imposed for numerical expediency but also helps to place an upper bound on its influence.
In nature, the environment likely has material concentrated in the orbital plane, formed from pre-dynamical mass loss during a Roche lobe overflow phase and by the early part of the CE plunge phase which was not included in our simulations.These earlier phases of mass loss have been studied in previous simulations, where they result in torus-like mass distributions with densities of order 10 −9 g cm −3 for systems similar to the one simulated in this work (MacLeod et al. 2018a,b;Reichardt et al. 2019). MacLeod et al. (2018b), simulated a case with secondary-to-primary mass ratio  = 0.3, finding that the primary loses ∼ 10% of its mass up to when the orbital separation equals the original primary radius, and about half of this remains bound.MacLeod & Loeb (2020) find that about 25% of the secondary mass is unbound during this time, roughly independent of .Taken together, this suggests that for the system we simulated ( = 0.5), ∼ 5-10% of the original primary mass, or ∼ 0.1-0.2M ⊙ , could have been present in the circumstellar environment in a bound torus-like structure.
Indeed, observations suggest that interaction between CE ejecta and pre-CE ejecta may be common.Matsumoto & Metzger (2022), citing Metzger & Pejcha (2017), model the light curves of luminous red novae (LRN) by considering both thermal emission and emission owing to recombination.However, they find that several LRNe are too bright to be explained by their model unless shocks between CE ejecta and predynamical ejecta are invoked.They point out that such predynamical mass loss likely explains certain observations, including the slow early rise in the light curves of some LRNe.Furthermore, similar shock-heating models have been proposed for other transient events involving high-mass progenitors.In many cases, narrow emission lines are observed, which may stem from the ejecta-circumstellar material interaction leading to ionization and subsequent recombination (e.g Smith 2013; Smith et al. 2015;Pearson et al. 2023).Our simulations (Model B) suggest that something analogous may happen during CE events involving lowto intermediate-mass progenitors.
Overall, our models suggest that the potency with which recombination assists envelope unbinding is blunted by two effects.First, released recombination energy appears to expand the envelope, creating a negative feedback loop which opposes the release of orbital energy.Secondly, reionization of recombined gas can act as an energy sink, preventing heating by the ambient environment that would otherwise assist envelope unbinding.Future work should explore these self-regulation effects in more detail.
Pathways allocation AST20034.This work also used the computational resources of the Rosen Center for Advanced Computing (RCAC) at Purdue University under XSEDE and ACCESS allocation TG-AST120060.Financial support for this project was provided by the Department of Energy grants DE-SC0020432 and DE-SC0020434, the National Science Foundation grants AST-1813298, AST-2009713 and AST-2319326, and the National Aeronautics and Space Administration grant 80NSSC20K0622.

C2 Masses
Fig. C1 shows the amount of mass involved in the various ionic transitions of helium during the simulation.The left, middle and right columns show the ionic transitions of He I, He II and He III tracer gas, respectively.Recall that, for example, the He I tracer follows envelope gas for which He I was the most common ionization state in the 1D MESA snapshot of the RGB star.From top to bottom, the rows indicate the mass of helium gas in a given tracer that is respectively He I, He II, or He III, at time .First, note that the He III tracer (right column) contains more than 10 times as much mass as the He II tracer, which, in turn, contains more than 10 times as much mass as the He I tracer.Note also that the sum of the masses of all ionic states of a given tracer is equal to the tracer mass in the simulation domain at time  (dashed green).Finally, note that the same quantities are plotted for both the total mass (black) and only that part of it which is bound (blue), where by "bound" we mean that it does not satisfy the local condition (3) to qualify as unbound gas at time .That part which is unbound is the difference between black and blue, which is not plotted to avoid clutter.Comparing the blue and black curves in a given panel gives an indication of what fraction of the recombination or ionization has occurred in material that was bound at the time of said recombination or ionization.We have no way of verifying whether gas that was bound (unbound) then is still bound (unbound) now, so the correspondence is imperfect.Nevertheless, given that only a small fraction of the envelope mass is unbound during the simulation (Fig. 4, solid lines) the difference is likely to be small.

C3 Energetics
Given that we are interested in explaining envelope unbinding, we are mainly interested in the energetics of the part of the envelope that is still bound at the time of recombination/ionization.Because the release (acquisition) of recombination energy during recombination (ionization) is local in space, we are interested in recombination/ionization that takes place in bound envelope gas.Therefore, we often focus on the blue curves in the discussion, but keep in mind that we technically mean gas that is bound at time , not necessarily at the time of recombination/ionization.From the right column of Fig. C1, we see that most of the recombination has happened in gas that is bound.By contrast, while most of the He I and He II tracer gas ionizes into He III (left and middle panels of bottom row), there is almost no net ionization of the bound component.
The energy release shown in Fig. 6 is calculated in the following way.For each cell in the simulation domain, we calculate the ionization state at time , including partial ionization, by making use of the Saha equation.This allows us to compute the total recombination energy for all the cells of each ionization state.Importantly, we do this for each ion tracer separately, adding up the recombination energy from all the cells of that tracer only.This gives us the total recombination energy of a given ionic species  that started out in ionic state .This is similar to what was done by eye by comparing the third and fourth rows of Figs. 1 and 2 in Sec.3.1.3.However, now we are integrating over the full volume of envelope gas (rather than just a slice) and including partial ionization.Recall that we were also able to estimate, by visual inspection of the fifth row of those figures, how much recombination takes place in envelope gas that remains bound.Here we can choose to include only envelope gas that is bound at time  in the integration, allowing us to plot the bound envelope gas separately.Values are calculated with respect to the MESA model.Initial values are nonzero for at least two reasons.First, the high temperature near the surface caused by imposing a higher pressure there (Sec.2.2) causes extra ionization in the initial condition, as mentioned in Sec.3.1.3.Second, partial ionization regions are not accounted for by the tracers, which are assumed to represent gas of a single ionization state of H and He (Sec.2.3).This effect apparently causes a slight overestimate of the ionization of the MESA star, and hence gives slightly positive values of the released recombination energy at  = 0.

C4 Direct calculation of total recombination energy release
Here we explain how the dashed curves in Fig. 8 are obtained.To make a fair comparison with the solid curves which considered only tracer gas that is still present in the simulation domain at time , we define  rec,env to be the difference between the recombination energy in the simulation domain and that of the same gas at  = 0.This initial value (which depends on ) is estimated as follows.We first calculate the mass of each ion tracer at time .Then we multiply by the hydrogen or helium mass fraction, divide by the hydrogen or helium atomic mass to obtain the number of atoms, multiply by the appropriate ionization energy, and sum the contributions from all hydrogen and helium ionic species.For example, the initial recombination energy of He III tracer gas at time  is estimated as  rec,He III (0) =  (He III)  m (He)  He (  He I→He II +  He II→He III ), where  (He III) is the mass of He III tracer gas at time ,  m (He) is the mass fraction of He (= 0.29),  He is the mass of the helium atom, and  He I→He II and  He II→He III are the ionization energies associated with the relevant ionic transitions.This gives us an estimate which neglects metals.This estimate is then subtracted from  rec,env ( ) to obtain  rec,env ( ).The initial value of − rec,env then turns out to be negative rather than zero, implying the presence of extra recombination energy at  = 0.This initial value matches quite well with the recombination energy of metals in the star at  = 0 that is calculated from the MESA model. rec,env is then redefined to equal zero at  = 0.The part of  rec,env associated with unbound gas is calculated using condition (3).The same procedure is done for the unbound gas only to obtain  rec,env,unb , and then we calculate the bound component using  rec,env,bnd =  rec,env −  rec,env,unb .

APPENDIX D: ROLE OF RADIATION
The leftover radiation energy that is present on account of the slightly different value of  rad mistakenly used when subtracting  rad  4 from the MESA EOS tables (see Sec. 2) is accounted for in the analysis.Considering that this leftover energy is just 0.015 per cent of the actual radiation energy that would be present assuming thermodynamic equilibrium, Fig. 8 implies that  rad  4 can become exceedingly large in Models A and B, where a dense ambient is present.In fact, this term can greatly exceed the gas thermal energy density 3 gas /2 in regions where envelope gas and ambient gas are mixing.In Fig. D1 we plot the ratio . 14From the figure, we see that even in the reduced ambient run, Model D,  rad  4 exceeds 3 gas /2 at the very centre where the core particles are located and in outward moving spiral density wakes.Thus, even though radiation was negligible in the initial envelope, it may become important during CEE.However, note that if radiation energy and pressure had been included in the simulation, then regions of such large  rad  4 would be less likely because (i) they would require extra energy to form, and (ii) radiation pressure would resist the concentration of radiation.Moreover, for low enough optical depth radiation transport would also impede the formation of such regions.

APPENDIX E: DEVIATION FROM PERFECT ENERGY AND ANGULAR MOMENTUM CONSERVATION
The evolution of the total energy and angular momentum in the simulation, including fluxes through the boundaries, are shown in Fig. E1 and Fig. E2, respectively.The total energy and angular momentum variation are small, and the differences between the runs smaller.In some cases these differences (over the entire domain) can be comparable to differences discussed in this work (restricted to bound gas or the core particle orbit).However, we note that the differences seen in Figs.E1 and E2 would, in any case, have the wrong sign to explain the deeper inspiral in Model A as compared to Model B. Model D is not plotted because there is no way to account for the energy and angular momentum change caused by the removal of ambient gas.Some non-conservation of energy and angular momentum is caused by  the exclusion by the multipole Poisson solver of the potential associated with gas with  >  box /2.Note that the sum of envelope gas and particle orbital energy (black in Fig. 9) eventually increases with time in both simulations.
In Model D the ambient energy is insignificant so transfer from the ambient cannot be causing the rise.The rise may be caused by the potential energy gain (toward less negative values) due to gas moving out of the sphere  =  box /2.
Since more gas leaves the box and hence the sphere  =  box /2 in Model D than Model A, the extra energy gain from this effect in Model D would then be comparable to the extra energy gained from the ambient in Model A. Imperfect energy and angular momentum conservation during the simulation could significantly affect the estimate for  CE of Section 3.7 if the particle energy and angular momentum are significantly impacted, which may or may not be the case.To obtain a conservative estimate of how much this could change our estimate for  CE , we considered the unlikely scenario that all of the deviation from non-conservation is from the particle orbital energy.We then recalculated the  CE estimates for Model D assuming energy and angular momentum non-conservation comparable to Models A and B. The gain in energy seen in Fig. E1 would then imply that  1−2 was in reality more negative, which leads to a smaller estimate for  CE (0.43 instead of 0.55).The loss of angular momentum seen in Fig. E2 would then imply that  ,1−2 was in reality less negative, which leads to a larger estimate for  CE (0.51 instead of 0.22).In both cases, the estimates remain within the range 0.27 <  CE < 0.55 mentioned in Section 3.7.

APPENDIX F: ADIABATIC INDEX
In stars, low values of the adiabatic index  1 = ( log / log ) ad can lead to dynamical instability.For example, in the idealized case of homologous motion in a constant density, constant  1 sphere, instability sets in for  1 < 4/3 (e.g.Hansen et al. 2004, chapter 8).In Figure F1, we present  1 in a

Figure 2 .
Figure 2. Continuation of Fig. 1 (but zoomed out to show the slice for the entire simulation domain) for times  = 23.1, 46.3 and 92.6 d.
Model A instead of Model B, for comparison.The left column of the first figure shows snapshots at  = 0 and the middle column shows  = 5.79 d.Each subsequent column, starting with the right column of Fig. 1 and moving left to right in Fig. 2, has a time equal to double that of the previous column: 11.57, 23.15, 46.30 and 92.59 d.Note that the panels of Fig. 1 are more zoomed in than those of Fig. 2.

Figure 3 .
Figure 3. Evolution of the orbital separation of the core particles for the three runs.The softening radius is shown for reference.The inset shows the evolution of the separation at late times.

Fig. 3
Fig.3shows the evolution of the core particles' separation .The inset shows a zoom-in of the late stage of the simulation.For the first ∼ 10 orbits (∼ 40 d), the particle separations are very similar in Models A ( = 5/3 EOS; red) and B (tabulated EOS; blue).Then Model B acquires a slightly smaller one period-averaged mean separation than Model A. This can be seen at  = 50 d.After that, the difference in the mean separations gradually reduces and flips sign near  = 65 d.By the end of the simulation, Model B has a slightly larger mean separation than Model A, as seen in the inset.The larger separation in the tabulated EOS run as compared to the ideal gas run is qualitatively consistent withSand et al. (2020),Lau et al. (2022a), González-Bolívar et al. (2022), which report on systems with different binary parameter values.The results ofReichardt et al. 2020 were less conclusive on this point.The separation curve of Model C is quite similar to those of Models A and B, but shows less correspondence with Model A than does Model B. Recall that in Model C we replaced the internal energy of the EOS tables with the thermal energy alone, akin to Model A, so the distinct evolutions between between Model C and Model A must be due to other model differences.The separation curve for the reduced ambient run Model D agrees very closely with that of Model A. This verifies that the orbital evolution is insensitive to the density and total mass of ambient gas.

Figure 4 .
Figure 4. Evolution of the unbound envelope mass for different definitions of "unbound."Model A (ideal gas EOS) is plotted in red, Model B (MESA EOS) is plotted in blue, and Model C (MESA EOS with recombination energy removed) is plotted in green or dotted black.Unbound mass according to fiducial definition of equation (3) (solid), as equation (3) but excluding the thermal energy density (dashed), or as equation (3) but also including the recombination energy density (dotted).Note the excess unbound mass in Model B as compared to Model A between  ≈ 13 d and  ≈ 70 d, caused by thermalization of released recombination energy.

Figure 5 .
Figure 5. Contributions to the energy of bound envelope gas where equation (3) defines unbound; Model B (MESA EOS, opaque lines) and Model A (ideal gas EOS, semi-transparent lines).

Figure 6 .
Figure 6.Energy released by each ionic transition of helium.Recombination results in positive values (top row) whereas ionization results in negative values (bottom row).The recombination of He III tracer gas into He II releases the most energy into bound envelope gas, and the recombination of He III tracer gas into He I also contributes significantly.From the top row, we also see that reionization of He I into He II absorbs substantial energy during the simulation.The non-zero values at  = 0 in some panels are explained in the main text.

Figure 7 .
Figure 7. Release of energy from recombination (left column) and ionization (middle column) of H (top row), He (middle row), and their combination (bottom row).The net energy release is plotted in the right column.Note that the ambient gas is excluded (as it should be).
d the mass of bound He II increases.Comparing the panels of this column, this increase is caused by both ionization of He I and recombination of He III.Unbound He I gas at  ≈ 30 d is mostly reionized to He III.As discussed below, this helps to explain why Model B does not unbind more mass than Model A at late times.Comparison with Fig.2suggests that the reionized envelope gas is mainly in the outer envelope and energized by the ambient medium, whereas the recombining bound gas is located deeper in the expanding envelope.The middle panel of the bottom row of Fig.6

Figure 8 .
Figure 8.Comparison of recombination energy release calculated using two different methods.The first uses the Saha equation and is shown as a solid black line, or a solid blue line for just the component which is bound at time .These lines are identical to those in the right panel of Fig. 7.The second method uses equation (5), and is shown by dashed lines.The green line shows the difference in the energy injected by the particles in the two runs, with the ideal gas run hosting a deeper inspiral and thus greater energy release.By the end of the simulation, the absence of recombination energy release in the ideal gas run is fully or partially compensated by the deeper inspiral.

Figure 9 .
Figure 9. Similar to Fig. 5 except now comparing Model A (ideal gas EOS, semi-transparent lines) with Model D (ideal gas EOS with reduced ambient, opaque lines), and including both bound and unbound envelope gas.

Figure 10 .
Figure 10.Comparison of Model A (left) and Model D (right), which restarts from Model A at  ≈ 20 d but drastically reduces the ambient gas after a few days.Panels show slices through the orbital plane.Top row: Mass density in g cm −3 .In Model A, the gas is confined by the ambient, whereas in Model D it expands freely.Bottom row: E env , defined in equation (4), is a measure of binding energy density: red gas is unbound and blue gas is bound (see Sec. 3.1.4).In Model A, there is mixing with the ambient medium (the ambient gas is uncoloured), whereas in Model D the ambient medium is not present the slice because it has been completely displaced by envelope gas.

Figure A1 .
Figure A1.Same as Fig. 4 but omitting the factors of 2 in the terms E pot,env−1 and E pot,env−2 .The horizontal black shows the total mass of envelope gas, accounting for the flux through the boundary.
= 92.6 d for Model A (left) and Model D (right)

Figure C1 .
Figure C1.Mass evolution for the various ionic transitions of helium in Model B, for all envelope gas (black) and envelope gas that is bound at time  (blue).Most of the helium gas is in the state He III at  = 0, and almost half of this He III recombines during the simulation to become He II or He I.The green line showing the total mass of the tracer reduces with time because envelope gas leaves the simulation domain through the boundary; this outgoing envelope gas is almost completely comprised of unbound gas, by our fiducial definition of equation (3).

Figure D1 .
Figure D1.Ratio of the quantity  rad  4 to the gas thermal energy density 3 gas /2 for the ideal gas run Model A (left) and ideal gas reduced ambient run Model D (right).The radiation energy is not explicitly included in any of the simulations.The temperature is estimated by assuming  = 0.62, though the results are insensitive to this choice.

Figure E1 .
Figure E1.Evolution of the total energy, accounting for energy contained in gas that passes through the domain boundary, for Model A (red), Model B (blue) and Model C (yellow).Negative jumps in the energy at  ≈ 25 d and  ≈ 50 d are caused by the sudden halving of the softening radius of both core particles.

Figure E2 .
Figure E2.Evolution of the total angular momentum, accounting for angular momentum contained in gas that passes through the domain boundary, for Model A (red), Model B (blue) and Model C (yellow).

Figure F1 .
Figure F1.The quantity  1 = ( log / log ) ad plotted in the orbital plane at  = 23.1 d (top) and  = 92.6 d (bottom).Purple contours show  1 = 4/3 and vectors show the gas velocity in the lab frame, projected onto the orbital plane.

Table 1 .
List of runs