© The Authors 2023

1. Introduction

Common envelope evolution (CEE) of a binary star system occurs when one of the stars engulfs its companion, which then rapidly spirals in through the envelope (Paczynski et al. 1976). The drag experienced by the companion moving in the noncorotating envelope leads to energy and angular momentum deposition in the surrounding gas. One possible outcome of CEE is that the companion star is dissolved inside the primary and the two stars merge into one. Alternatively, the dynamical inspiral slows down and a quasi-steady spiral-in phase ensues. The reasons for the stalling of the inspiral are not completely understood, but the reduction of the drag does occur when the density decreases either due to envelope expansion or heating, or when the gas starts to locally corotate with the companion (Roepke & De Marco 2023). Simulations show that this stalled inspiral phase lasts for at least hundreds of orbits and that the two cores remain surrounded by a shared envelope (e.g., Ricker & Taam 2012; Passy et al. 2012; Ohlmann et al. 2016; Ivanova & Nandez 2016). It is believed that a self-regulating feedback loop of the weak local drag and its associated energy dissipation slowly brings the central binary together and eventually ejects the envelope on its thermal timescale leaving behind a post-CEE binary (e.g., Ivanova et al. 2013a; Clayton et al. 2017; Glanz & Perets 2018).

Common envelope evolution is responsible for a wide variety of binary systems such as cataclysmic variables (Paczynski et al. 1976), X-ray binaries (e.g., Kalogera & Webbink 1998; Taam & Ricker 2010; Chen et al. 2020), progenitors of Type Ia supernovae (e.g., Iben & Tutukov 1984; Belczynski et al. 2005; Ablimit et al. 2016), or planetary nebulae nuclei (e.g., De Marco 2009; Jones & Boffin 2017). CEE might be responsible for a substantial fraction of gravitational wave progenitors (e.g., Dominik et al. 2012; Klencki et al. 2021; Marchant et al. 2021), but CEE is also expected to emit gravitational waves on its own that are likely to be detected by space-based gravitational-wave detectors such as LISA (Thorpe et al. 2019) or TianQin (Huang et al. 2020) during the post-dynamical in-spiral CEE stage (Renzo et al. 2021). Binaries that do not survive CEE and merge can be observed as luminous red novae (e.g., Soker & Tylenda 2006; Ivanova et al. 2013b; Kochanek et al. 2014; Pejcha et al. 2016a; Blagorodnova et al. 2021).

Despite its importance, CEE is far from being fully understood. Great efforts have been made to confront numerical simulations’ outcomes to observational constraints over the last few decades. Three-dimensional hydrodynamical simulations have provided insight into the physical processes important in the dynamical inspiral (e.g., Passy et al. 2012; Ohlmann et al. 2016; MacLeod et al. 2018; Chamandy et al. 2020; Sand et al. 2020), the CEE ejecta dynamics and thermodynamics (e.g., Glanz & Perets 2018; Iaconi et al. 2019, 2020), or radiation hydrodynamics of the ejecta and the associated transients (e.g., Pejcha et al. 2016a,b, 2017; Metzger & Pejcha 2017; MacLeod et al. 2017a; Matsumoto & Metzger 2022). However, many ab initio simulations fail to eject the common envelope during dynamical plunge-in when only orbital energy injection by the secondary star is considered, and the obtained post-dynamical inspiral orbital separations are often larger than that observed in post-CE systems (e.g., Nebot Gómez-Morán et al. 2011; Iaconi & De Marco 2019; Politano 2021; Passy et al. 2012; Kruckow et al. 2021). A more realistic equation of state that takes ionization states into account seems to facilitate mass ejection, but might not affect the final separation (e.g., Nandez et al. 2015; Reichardt et al. 2020; Lau et al. 2022a,b). Because of the wide range of temporal and spatial scales that need to be resolved and the associated high numerical cost, 3D hydrodynamical simulations are often stopped soon after the end of the dynamical inspiral phase. The consecutive slow contraction of the orbit on a thermal timescale necessitates reverting to 1D models that cannot capture the multidimensional features (Taam et al. 1978; Meyer & Meyer-Hofmeister 1979; Fragos et al. 2019).

To facilitate rapid prediction of outcomes, the binary configurations preceding and following CEE are often linked using energy conservation with one free parameter, α_CEE (e.g., Webbink 1984; Livio & Soker 1988). The value of α_CEE can be estimated from simulations or from various observed binary populations. Yet, it is not clear whether the α_CEE formalism can truly encompass the complicated physics of CEE. In particular, if thermal timescale processes such as predynamical nonconservative mass transfer or post-dynamical self-regulated inspiral are important, then adiabatic energy conservation is violated. Of course, it is often possible to select a value of α_CEE to explain an observed population even if some assumptions of the formalism are not satisfied. This displeasing situation has motivated the development of alternative formalisms based on the conservation of angular momentum (Nelemans et al. 2000; Di Stefano et al. 2023) or a two-step prescription combining energy and angular momentum (Hirai & Mandel 2022).

The post-dynamical self-regulating inspiral plays an important role in many of the unsolved aspects of CEE. Due to numerical difficulties in studying late stages of CEE when the envelope has expanded and the central binary has tightened, very little is known about this phase, especially the duration, mechanism of orbital contraction and angular momentum transfer, and whether the thermal-timescale self-regulation is actually established. Clearly, even if the central binary orbits in a locally corotating gas and the gravitational drag, which is the prevailing source of orbital tightening during dynamical in-spiral, becomes very weak (e.g., Ostriker 1999; Ricker & Taam 2012; MacLeod & Ramirez-Ruiz 2015; MacLeod et al. 2017b; Reichardt et al. 2019; Chamandy et al. 2019b; De et al. 2020), the corotation cannot be maintained over arbitrary distances. As a result, the complex interaction between the binary and the gravitationally perturbed shared envelope can take over and drive the orbital separation evolution on a shorter timescale than the thermal one. An additional issue is the possibility that mass and angular momentum can reaccrete onto the binary.

The configuration of the post-dynamical inspiral resembles a very thick circumbinary disk (CBD), where a binary is embedded in a low density cavity surrounded by a disk with which it interacts by gravitational, advective, and viscous torques, mass accretion onto the central binary, and by binary eccentricity evolution (Sandquist et al. 1998). In the case of thin CBDs, such intricate interactions can lead to either orbital expansion or contraction and to excitation of the binary eccentricity (e.g., Artymowicz & Lubow 1994; MacFadyen & Milosavljević 2008; Shi et al. 2012; Tang et al. 2017; Miranda et al. 2017; Muñoz et al. 2019; Muñoz & Lithwick 2020; Duffell et al. 2020; D’Orazio & Duffell 2021; Dittmann & Ryan 2021; Penzlin et al. 2022). These results suggest an exciting possibility that the CEE post-dynamical evolution does not have to proceed as a simple monotonic contraction of a circular orbit, but there can be phases of orbital expansion or an eventual formation of an eccentric post-CEE binary. Our connection between post-dynamical CEE and CBDs is different from previous explorations of fallback CBDs around post-CEE binaries (De Marco et al. 2011; Kashi & Soker 2011).

There are also important differences between CBDs and the post-dynamical phase of CEE. CEE might not result in the formation of a cavity around the central binary, instead, the central binary could virialize the gas in its vicinity, which would provide pressure support of the envelope and prevent accretion. Therefore, we can identify two extreme regimes of zero or maximum accretion onto the binary. Which of the two regimes of accretion occurs depends on the absence or presence of a “pressure valve” inside the orbit, which allows the material to accrete onto the binary (Chamandy et al. 2018). An example of such a pressure valve could be jets (e.g., Soker & Livio 1994; Moreno Méndez et al. 2017; Shiber et al. 2019; López-Cámara et al. 2019, 2022). A realistic situation probably lies between these two extreme regimes of accretion.

In this paper, we aim to clarify the nature and dynamics of the post-dynamical inspiral of CEE by performing the first dedicated series of 3D hydrodynamical numerical simulations. To establish a well-controlled numerical experiment, we mimic the outcome of the dynamical inspiral phase by artificially injecting angular momentum in the primary envelope following the procedure of Morris & Podsiadlowski (2006). To follow the evolution of the system over long timescales, we excise an inner sphere containing the binary, but study the gravitational influence of the orbiting binary on the surrounding envelope by prescribing time-changing gravitational potential. The inner boundary condition at the excised sphere allows us to control the accretion on the central binary. To analyze our results, we employ techniques and diagnostics inspired by those commonly used in the context of CBDs (e.g., Shi et al. 2012; Miranda et al. 2017; Muñoz et al. 2019; Penzlin et al. 2022).

This work follows the following structure: in Sect. 2, we introduce our physical model and describe the numerical setup used in our common envelope simulations. In Sect. 3, we present the results of our simulations. In particular, we measure the timescale of binary separation evolution resulting from the various torques acting on the system, when accretion is turned on or off. We measure the typical frequencies associated with the short-term variability of mass accretion onto the binary, and we compare them with that from CBDs. We then study the formation of overdensities, the excitation of eccentricity, and the convective stability of the envelope. Finally, we analyze the angular momentum transport within the envelope. In Sect. 4, we discuss implications of our findings for CBDs and CEE. In Sect. 5, we summarize our results.

2. Physical model and numerical setup

We construct our post-dynamical inspiral model in the inertial frame at rest with the center of mass of the binary. We do not follow the previous evolution of the inspiraling binary, instead, we mimic its outcome following a procedure similar to Morris & Podsiadlowski (2006, 2007, 2009) and Hirai et al. (2021), where the envelope is spun-up until a satisfactory amount of total angular momentum is injected. This mimics angular momentum transfer from the secondary’s orbit into the envelope during the dynamical plunge-in (see Sect. 2.4). We set the gravitational constant G, the total binary mass M = M₁ + M₂, the primary’s initial radius R, and thus the angular velocity to unity. The orbital velocity is fixed to , where a_b = r₁ + r₂ is the fixed binary separation, M₁ is the mass of the primary’s core located at {r₁, θ₁, φ₁}, and M₂ is the mass of the secondary object (either a main-sequence star or a compact object) located at {r₂, θ₂, φ₂}. Orbital period is P_orb = 2π/Ω_orb. The two objects are not resolved and are considered as constant point masses. To simplify our model, we consider an equal mass binary (q ≡ M₂/M₁ = 1) on a fixed circular orbit. The mass of the envelope is M_env = 2 in our units. Because we are most concerned with the angular momentum transport within the common envelope in the two extreme regimes of mass and angular momentum accretion onto the binary rather than the specific details of the individual cores, we excise a central region encompassing the binary, which has a radius r_in = 0.625 a_b = R/10. This excised region represents the gas bubble virialized by the orbiting binary and the enforced conditions at its boundary determine whether the binary is accreting or not.

We compare our setup to several ab initio simulations of CEE in Table 1. The key quantity is the ratio of final separation to the initial radius of the primary, which we set in our model to a_b/R = 0.16. This comparison suggests that our choice of initial parameters to the binary and envelope broadly represents results of ab initio simulations across for a range of progenitor types.

In the rest of this Section, we describe the equations used for solving the problem (Sect. 2.1), boundary conditions (Sect. 2.2), initial conditions (Sect. 2.3), and initial deposition of angular momentum (Sect. 2.4). We then present our numerical setup for the averaging of the polar zones (Sect. 2.5), mesh refinement (Sect. 2.6), and equatorial symmetry of the simulations (Sect. 2.7).

2.1. Equations of hydrodynamics

We use Athena++ (Stone et al. 2020) to solve the equations of hydrodynamics

(1)

(2)

(3)

where E = e + ρu²/2, e is the internal energy density, P = (Γ − 1)e, Γ = 5/3 is the adiabatic index, Φ(r) is the gravitational potential of the binary,

(4)

and T_ij is the symmetric viscous stress tensor,

(5)

which is nonzero when we prescribe a kinematic viscosity ν.

For runs with nonzero viscosity, we prescribe an isotropic effective kinematic viscosity of turbulent nature, , where v is the velocity of the turbulent eddies, and l is their vertical mean free path or the mixing length. We further assume that the mixing length is proportional to the local pressure scale height, l = α₁H_P (e.g., Vitense 1953; Zahn 1989), and that the characteristic eddy velocity is a fraction of the local sound speed, v = α₂c_s. Following Shakura & Sunyaev (1973), we obtain ν(r, θ, φ) = α_νc_sH_P, where α_ν = α₁α₂/3 is a free parameter. Assuming a typical effective kinematic viscosity of 𝒪(10¹⁵ cm² s⁻¹), we take α_ν = 10⁻³ in this work. The nature of such effective viscosity is unknown and has been debated lively in the context of astrophysical accretion flows. We further impose zero kinematic viscosity radial gradient in ghost cells at both boundaries. In our viscous simulation, we use the Runge-Kutta-Legendre super-time-stepping algorithm in Athena++ (Meyer et al. 2014; Stone et al. 2020), which integrates diffusive terms forward with hyperbolic timesteps. Although this algorithm dramatically reduces the timestep constraints, we still could not evolve our viscous runs for as long as we could when α_ν = 0.

2.2. Boundary conditions

2.2.1. At the inner edge of the domain

During the initial spin-up of the star (see Sect. 2.4), we assume that the inner boundary supports the weight of the primary’s envelope and we forbid the fluid to flow through it. To achieve that, we assume that ρ is constant in ghost cells, which we initialize with the value of the density in the first interior cell i and which we assume to be in equilibrium with ghost cells. Considering a simple first-order finite volume integration algorithm with a number of ghost cells N_g = 2, the pressure and density boundary conditions in the inner ghost cell of index j read

(6)

(7)

where is the cell-face value of the variable f evaluated at , Δr = r_j − r_j + 1, and ⟨Φ⟩_θ is the time and latitude average of the binary gravitational potential (see Sect. 2.3). We deal with the horizontal velocity by applying zero radial gradient in the adjacent ghost cells and we impose reflecting radial velocity to enforce that u_r = 0 at the inner boundary,

(8)

(9)

(10)

We note that in some cases Eq. (7) gives negative P in ghost zones. When that is the case, equilibrium cannot be enforced at the boundary, and we impose a zero pressure gradient instead.

Once the spin-up phase is terminated, we either allow or forbid accretion onto the central binary. When accretion is allowed, we open the inner boundary to angular momentum and mass flow by imposing zero radial gradient of ρ, P, u_θ, and angular momentum in ghost cells, and a diode-type radial infall only condition, u_{r, i − j} = min(u_r, i, −u_{r, i + j − 1}). When accretion is forbidden, we impose purely reflecting boundary conditions.

2.2.2. At the outer edge of the domain

We use diode-type boundary conditions at the outer edge of the domain and we impose zero density and pressure gradient in the outer ghost zones. However, this condition implies that the ambient medium of our initial model is out of equilibrium and there is an inflow near the outer boundary during the initial spin-up phase. In this region, ρ and P are initially very low (see Sect. 2.3) and thus there is negligible influx of mass.

2.3. Initial conditions and outer low-density medium

We assume that the gas in the envelope is initially in hydrostatic equilibrium and that it can be described by a polytropic equation of state, as is often done in stellar physics (e.g., Maeder 2009; Jones et al. 2009; Gagnier & Rieutord 2020). Ignoring the gas self-gravity and considering purely radial initial profiles, the equations governing the envelope initial structure read

(11)

where K is a constant related to the thermal conditions at the inner boundary. The Green’s function for Eq. (4) satisfies (Jackson 1975)

(12)

for r ≥ r_in ≥ max(r₁, r₂), where are the usual normalized scalar spherical harmonic functions of degree ℓ and order m. The parity properties of the spherical harmonic function and of the time and latitude average of the binary potential, that is ⟨Φ(r)⟩_θ = ⟨Φ(−r)⟩_θ, imply that

(13)

where ⟨ ⋅ ⟩_θ indicates a time and latitude average. The ℓ ≤ 2 latitude and time averaged binary potential finally reads

(14)

We use this latitude- and time-averaged potential in the momentum equation during spin-up. We replace the averaged potential by its time-dependent expression (Eq. (4)) after the spin-up. Our choice of the initial potential facilitates the transition from the initial spin-up by preventing a large injection or removal of gravitational energy, which could lead to a sudden envelope ejection or collapse. By combining Eqs. (11) and (14) and by specifying the ratio between ρ at the stellar surface and at the inner boundary surface, κⁿ = ρ(R)/ρ(r_in), we obtain the initial density and pressure profiles

(15)

(16)

where n = 1/(Γ − 1) is the polytropic index and

(17)

Density at the inner boundary ρ(r_in) can be calculated from the prescribed total mass of the envelope.

The primary star is initially embedded in a low-density medium to which we apply our outer boundary conditions and in which the envelope will expand later on. To model this low-density medium, we consider an atmosphere in hydrostatic equilibrium with constant ambient sound speed c_s, amb (see also MacLeod et al. 2018) and we obtain analytical ρ and P profiles assuming , which yields

(18)

We note that it is not possible to transition from the envelope to the ambient region without a discontinuity, either in ρ, P, or both. To accommodate this discontinuity, we derive the constant C₁ such that the stellar surface is in hydrostatic equilibrium with the low-density atmosphere, which gives

(19)

Here, P⁻ and ρ⁻ are the pressure and density in the last radial cell of the envelope from Eq. (15), r⁺ is the radius of the first cell of the ambient medium, Δr(R) is the difference between r⁺ and the radius of the last radial cell of the envelope, and . We illustrate our initial conditions in Fig. 1.

Fig. 1. Initial density, pressure, and sound speed profiles used in all our simulation runs.

In order to minimize the effects of the nonexact numerical hydrostatic equilibrium resulting from the finite grid resolution, we use Gauss-Legendre quadrature to map the initial profiles onto the mesh as volume averaged variables at the volume averaged center of each cell, which is different from geometric center in polar-spherical coordinates, especially near the polar axis because of the converging grid geometry.

2.4. Initial spin-up

We aim to construct a model with an initial total angular momentum that is consistent with what is available in the system,

(20)

where a_i is the initial binary separation¹, a_b is the enforced separation at the end of the dynamical plunge-in, M_env = 2 is the total mass of the envelope, and β is the ratio between the primary’s envelope angular momentum and the orbital angular momentum before the plunge-in. is the orbital angular momentum of the binary at the end of the dynamical plunge-in, which coincides with the beginning of our simulations, and μ = M₁M₂/(M₁ + M₂) is the reduced mass. We require β ≤ 1/3 to ensure Darwin stability (e.g., Hut 1980). To impart angular momentum to the envelope, we use the procedure of Morris & Podsiadlowski (2006, 2007, 2009) and we apply a fixed spin-up rate to all cells in which the angular velocity is sub-Keplerian, . Simultaneously, the structure of the envelope slowly restructures.

We stop the spin-up once a satisfactory amount of total angular momentum is injected in the envelope. After a short adjustment phase, we replace the latitude and time averaged potential ⟨Φ⟩_θ with its real expression (Eq. (4)). As a result, there is a small bump of internal energy that is exclusively due to the increase of the gravitational energy density in the inner envelope. Though it has no physical origin beyond the sudden anisotropy of the gravitational potential and despite the fact that its amplitude cannot be easily constrained, it has the benefit of mimicking a small gravitational energy deposition by the spiral-in of the secondary star. We discuss this more in Sect. 3.1.

2.5. Polar averaging

It is well known that the use of spherical coordinates leads to strong time-step constraints resulting from the converging grid geometry and the Courant-Friedrichs-Lewy (CFL) condition. To mitigate this issue, we use a polar averaging technique based on the Ring Average technique of Zhang et al. (2019), which is conservative and computationally inexpensive. This technique consists of a post-processing treatment of the variables in cell “chunks” adjacent to the polar axis, which is applied after the cells have been updated by the Riemann solver. Hence, this technique does not involve the modification of the grid nor of the solver. For nonuniformly spaced spherical coordinates, we compute the appropriate number of chunks N_c = 2^k per latitudinal ring of index m in each mesh block, where

(21)

Here, square brackets indicate rounding to the nearest integer and Δφ_block is the azimuthal extent of the mesh block. Then, we average conserved variables in the azimuthal direction within each chunk of each ring and in each mesh block. We subsequently apply second-order spatial reconstruction procedure to the averaged values and we correct the minimum time-step within a mesh block to account for the coarsened effective mesh.

2.6. Mesh refinement

Our initial models are statically refined to properly resolve regions with strong initial gradients. These are regions close to the central binary and to the initial surface of the star. Specifically, our initial mesh is refined two levels above the base in the regions r_in ≤ r ≤ 0.25 and 0.95 ≤ r ≤ 1.05. After the first timestep of the initial spin-up of the envelope, we switch from static to adaptive mesh refinement. We adopt a criterion based on the second derivative error norm of a function σ of a variable (Lohner et al. 1987). This criterion measures the smoothness of the solution for a given refinement variable. Similarly to the PLUTO code (Mignone et al. 2012), a mesh block is tagged for refinement whenever

(22)

Here, Δ_r, ±1/2 = ±(σ_i ± 1 − σ_i) and σ_r, ref = |σ_i + 1|+2|σ_i|+|σ_i − 1|. The value of the threshold χ_r is problem dependent and also depends on the chosen refinement variable σ. Finally, ϵ acts as a filter preventing refinement in regions of small ripples. We find that for our simulations, σ = ρ|u| tends to capture the flow contrasts the best with and ϵ = 0.01. We note that the criterion in Eq. (22) does not include cross derivatives, unlike the original work of Lohner et al. (1987). We find little difference when those terms are included and we opt not to include them.

2.7. Equatorial symmetry

Since the orbiting cores are aligned in the equatorial plane at all times, our setup should be exactly symmetric about the equator. In practice, such symmetry can be difficult to enforce despite Athena++’s integration method being well-suited to preserving it Stone et al. (2020). For instance a mesh symmetric about the double precision rounding accuracy of π/2 in the θ-direction yields asymmetric volume averaged cell colatitudes (Mignone 2014). For example,

(23)

is asymmetric because of the asymmetry of the trigonometric functions about the rounded value of π/2. This introduces asymmetry in the theoretically symmetric source terms in our problem. Volume-averaged colatitudes are thus computed about the double precision rounding of π/2, hereafter noted . We obtain

where is the latitude measured from . Although this change considerably improves symmetry, face-centered cell colatitudes are, in practice, not symmetric in the last place precision. Such tiny asymmetry of dθ leads to the asymmetry of the physical and geometric sources terms and to a residual asymmetric flow that may amplify when it is linearly unstable. Furthermore, additional sources of asymmetries may amplify the problem, such as compiler value-unsafe optimizations of floating-point operations or the nonassociativity of floating-point arithmetic. In order to control such inevitable perturbations, we impose ad hoc initial random weak seed perturbation to the initial density profile with maximum amplitude 10⁻⁶ ρ(r, θ, φ), which is orders of magnitude larger than the amplitude of perturbations resulting from grid asymmetries.

3. Results

We used a total of 4.6 million CPU hours on the Karolina cluster at IT4Innovations to perform our simulation runs. In Table 2, we summarize the parameters of the runs. In Fig. 2, we present zoomed-in snapshots of density cross section in the xy and xz planes at different times and for three inviscid simulation runs. The inviscid runs A, B, and C only differ by the initial size of the envelope’s angular momentum reservoir. Run A is computed with β = 0.3, that is close to the limit of Darwin instability, run B is computed with β = 0.1, and run C with β = −0.3. Negative value of β implies that the total z component of angular momentum is smaller than the initial orbital angular momentum of the binary orbit. Although our setup only approximates the process of angular momentum transfer from the orbit to the primary’s envelope during the dynamical plunge-in, the density structure and flow morphology early in our simulations have striking resemblance with late-time snapshots from ab initio simulations of dynamical plunge-in (e.g., Ohlmann et al. 2016; Chamandy et al. 2020).

Fig. 2. Zoomed-in snapshots of density cross section in the xy and xz planes at different times and for our three inviscid simulations runs A (left), B (middle) and C (right). The snapshots on the first line are taken shortly after the end of the initial spin-up phase.

Our simulations show that overall the envelope is destabilized by the central binary gravitationally torquing the inner envelope, exciting spiral density waves, and shearing the fluid flow. Energy from such flow is transferred to large-scale turbulence, and angular momentum is then transported by mean and turbulent flows. In the rest of this section, we investigate these processes in detail. We address the initial jump in energy (Sect. 3.1), binary evolution and mass accretion (Sect. 3.2), short timescale dependence of accretion (Sect. 3.3), presence and origin of the lump (Sect. 3.4), eccentricity of the envelope (Sect. 3.5), and convective stability and angular momentum transport (Sect. 3.6).

3.1. Energy injection

In Fig. 3, we show the kinetic, internal, gravitational, and total binding energy evolution for model A. We first discuss the bump in energy, which occurs at the end of the initial spin-up when we replace the latitude and time averaged binary potential with its real expression.

Fig. 3. Kinetic, internal, gravitational, and total binding energy evolution in units GM2/R for model A. The vertical black dashed line indicates the replacement of the averaged binary potential with its full expression (Eq. (4)), and the bump in energy is indicated with the double arrow.

To asses the importance of the bump, we estimate the CEE efficiency parameter α_CEE corresponding to injection of internal energy ΔE,

(24)

where a_i is the initial binary orbital separation before plunge-in that we assume to be equal to 10 a_b (e.g., Passy et al. 2012; Ohlmann et al. 2016; Chamandy et al. 2020). That is, the amplitude of the internal energy bump corresponds to a gravitational energy deposition during plunge-in of 46% of the difference between initial and final total orbital energy. Furthermore, because the difference between averaged and real gravitational potentials is only significant in the vicinity of the central binary, the energy is almost exclusively injected in the inner part of the envelope, which agrees with numerical simulations of Chamandy et al. (2019a). Because both amplitude and location of the energy injection are consistent with orbital energy deposition during CEE, we do not add more.

3.2. Binary evolution and mass accretion

Here, we address the evolution of the orbit of the central binary. So far, CEE theory has assumed that the binary separation decreases almost monotonically in time. However, recent studies of CBDs (e.g., Muñoz et al. 2019; Penzlin et al. 2022) suggest that for equal mass binaries there is a wide range of viscosity and disk thickness that leads to the expansion of the orbit. Therefore, finding out what actually happens to the binary separation in post-dynamical CEE inspiral is of fundamental importance. In our setup, we keep the orbital parameters fixed, but we can measure how much angular momentum was exchanged between the binary and the envelope and therefore assess what would happen to the binary if it was self-consistently coupled to the envelope. The advantage of our setup is that by excising the inner region, we can run the simulations for more orbits of the central binary.

3.2.1. Torques and angular momentum conservation

In order to predict the secular evolution of the binary separation, it is necessary to evaluate the torques in the common envelope. Such torques originate from the quadrupolar component of the gravitational potential as well as the advective (and perhaps viscous) angular momentum fluxes through the domain boundaries. The angular momentum conservation equation reads

(25)

where is the advective torque associated with the loss of angular momentum through the boundaries,

(26)

is the gravitational torque exerted by the binary,

(27)

and is the viscous torque,

(28)

Here, n_⊥ is the outward-pointing unit vector at the boundaries’ surface and s = r sin θ is the radial cylindrical coordinate. We give more details in Appendix A.

In Fig. 4a, we show the evolution of these torques for runs A and A′. We also perform consistency check for angular momentum conservation by comparing the time evolution of the individual torques with the evolution of the total angular momentum budget for all of our models after the initial spin-up, and we show the result for model A in Fig. 4b. We find that the angular momentum is conserved to within about 0.1 − 1% margin for all of our models. We also see that for all our models the total angular momentum evolution is dominated by the outflow at the outer boundary, which results from the expansion of the envelope and the finite radial extent of our numerical domain. When the inner boundary is open to angular momentum and mass flow toward the binary, angular momentum accretion dominates over the gravitational torque, which only weakly contributes to the injection of the angular momentum in the envelope. Because we choose zero radial gradient of angular momentum and viscosity at the inner boundary, the contribution of viscous torque remains weak even for eccentric flows in the binary close vicinity or for larger values of α_ν. After 140 orbital periods, we consider the various torques to be sufficiently time-steady so that we can qualitatively assume that all initial transients have decayed and that the flow properties have reached a quasi-steady state.

Fig. 4. Angular momentum evolution and conservation in our simulations. Panel a: time evolution of the advective, viscous, and gravitational torques for runs A and, A′ (dotted lines). Panel b: relative difference between the sum of the torques and the measured time derivative of the total angular momentum showing the angular momentum conservation for run A.

3.2.2. Binary orbital evolution

In this work, we set the orbital eccentricity e_b to zero and impose the binary mass ratio q = 1. We thus assume that mass and angular momentum accretion through the inner boundary distribute equally between the two cores, . Furthermore, we assume that accretion does not excite orbital eccentricity, as suggested by CBD simulations, and we therefore fix ė_b = 0 (Muñoz et al. 2019; Heath & Nixon 2020; Penzlin et al. 2022). The validity of this assumption will be discussed in Sect. 3.5. The time derivative of the binary’s angular momentum can be written as the orbital separation evolution equation

(29)

If the central binary does not accrete from the shared envelope, Eq. (29) simplifies to

(30)

and the binary orbit contracts (ȧ/a < 0) if , that is if the gravitational torque transfers angular momentum from the binary orbit to the envelope. If the inner boundary is open to mass and angular momentum flow onto the binary, it is useful to consider the specific angular momentum transfer rate

(31)

which yields a critical value j_crit = 3/8 for q = 1 (e.g., Miranda et al. 2017; Moody et al. 2019; Dittmann & Ryan 2021; Penzlin et al. 2022). Above this value, the binary orbit expands (ȧ/a > 0) and below it contracts (ȧ/a < 0). In Eq. (31), Ṁ is the measured mass accretion through the inner boundary

(32)

and .

Contrary to 3D simulations of CBDs where mass and angular momentum accretion only occur within a limited angle about the orbital plane dictated by the geometrical thickness of the disk, mass accretion could span the whole solid angle in our simulations. Hence, we may observe accretion along the polar axis with very small j, which could favor the contraction of the orbit according to Eq. (31). To diagnose this issue, we show in Fig. 5 the time average over 10 P_orb of the angular distribution of mass and angular momentum accretion fluxes through the inner boundary for model B. We see that mass and angular momentum accretion mostly occur within an annular ring centered on the orbital plane. Above and below such annular ring, mass accretion is accompanied with weak angular momentum accretion. Hence, the geometry of the CEE problem favors contraction of the binary orbit unlike what is the case for CBDs.

Fig. 5. Time average of the angular distribution of mass (panel a) and angular momentum (panel b) accretion fluxes through the inner boundary for ten orbital periods of model A at 148 ≤ t/Porb ≤ 158.

In the top panel of Fig. 6, we show the evolution of the moving average of the specific angular momentum transfer rate over one orbital period

(33)

Fig. 6. Evolution of key quantities relevant for the binary orbit after the initial envelope spin-up and adjustment. Panel a: moving average of the specific angular momentum transfer ⟨j⟩P (Eq. (32)) for models A, B, C, and D. The black dashed line indicates the critical specific angular momentum transfer jcrit = 3/8. Panel b: moving average of the gravitational torque for simulation runs A and A′. The black dashed line indicates . Panel c: moving average of the mass accretion rate through the inner boundary ⟨Ṁ⟩P. Panel d: envelope mass in our numerical domain Menv(r ≤ 10). Panel e: moving average of the orbital separation evolution timescale .

We see that the combined effects of gravitational torque and mass and angular momentum accretion lead to the contraction of the orbit for all the considered values of β and for both viscous and inviscid fluids. In the second panel of Fig. 6, we show the evolution of the moving average of the gravitational torque for all our models. We see that rapidly settles to a value that closely oscillates around zero and thus does not contribute to the orbital evolution when the inner boundary is open to mass and angular momentum flow toward the binary. Conversely, when accretion is prevented by reflecting boundary conditions decreases much slower, remains positive, and thus drives orbital contraction. This crucial difference comes from the stabilizing effect of higher density in the vicinity of the binary when reflecting boundary conditions are enforced. This is discussed in more depth in Appendix B. In the bottom panel of Fig. 6, we show the orbital separation evolution timescale . We find that τ_ab reaches a statistically steady value of O(10³ P_orb) for all models allowing mass and angular momentum accretion onto the binary. Conversely, when accretion is forbidden (simulation run A′), the gravitational torque is exclusively responsible for the orbital contraction and the slow decrease of implies a slow increase of τ_ab. The orbital separation evolution timescale eventually reaches a value of O(10⁵ P_orb) at the end of simulations run A′ at t ≈ 450 P_orb. It is possible that τ_ab would continue increasing if we were able to run our model for more orbits.

Finally, we address the influence of envelope viscosity on binary evolution. Unfortunately, we could not run simulation run D for as long as the inviscid ones. Still, we can see that α_ν = 10⁻³ does not significantly affect τ_ab. Because of the very variable nature of mass and angular momentum accretion rates, it is not clear whether the limited impact of viscosity would eventually lead to a slower or faster contraction of the orbit. Similarly, higher values of α_ν should be investigated as well.

3.3. Time variability of mass and angular momentum accretion

Now that we have investigated the secular binary separation evolution, we more thoroughly analyze the gas dynamics in the vicinity of the central binary, in particular, the time variability of mass and angular momentum accretion in the simulations that permit accretion. In Figs. 7 and 8, we show the latitudinal space-time diagram of the mass and angular momentum fluxes onto the binary, normalized by their maximum value in the considered time interval, for runs A and B. In the top panel of Figs. 9 and 10, we show a more detailed view of a shorter time interval. To construct these plots, we increased the simulation output rate to 80/P_orb. We see that mass and angular momentum fluxes exhibit periodic variability at all colatitudes. This variability is manifold: we observe a high frequency variability that is modulated by a lower frequency, at least near the orbital plane.

Fig. 7. Space-time diagram of the mass (a) and angular momentum (b) fluxes onto the binary during 21 orbital periods for model A.

Fig. 8. Same as Fig. 7 but for model B during 26 orbital periods.

Fig. 9. Detailed view on the variability of mass flux for model A. Panel a: space-time diagram of the local mass flux through the inner boundary. Panel b: time evolution of the mass accretion rate onto the binary. Panel c: power spectral density of the total mass accretion rate onto the binary (black line) and of the mass flux at each colatitude (colored lines). Green lines correspond to the range 0 ≤ θ ≤ π/3, blue lines to the range 2π/3 ≤ θ ≤ π, and orange lines to the range π/3 ≤ θ ≤ 2π/3.

Fig. 10. Detailed view on the variability of mass flux for model B and its two different accretion regimes (left and right panels). Meaning of symbols and lines in each panel is the same as in Fig. 9.

To identify the modes associated with mass and angular momentum accretion variability, we use Fourier transform to compute the power spectral density of the mass accretion rate Ṁ and mass flux ṁ(θ). We distinguish four colatitude ranges: 0 ≤ θ ≤ π, 0 ≤ θ ≤ π/3, 2π/3 ≤ θ ≤ π, and π/3 ≤ θ ≤ 2π/3 for the mass flux. We show our results in Figs. 9c and 10c for simulation runs A and B. In both simulations we identify two main peaks and their harmonics: one located at ω_b ≃ 2 Ω_orb and the other one at ω_ρ ≃ Ω_orb/5. Here, ω_b corresponds to the forcing angular frequency of the quadrupolar moment contribution to the binary potential for a binary mass ratio q = 1, which is the frequency at which material is pulled toward the central binary. The frequency ω_ρ is also often seen in CBDs, where it corresponds to an overdensity in the envelope orbiting at the angular frequency ω_ρ. This overdensity is often called “lump” and typically forms when some of the accreting material is strongly torqued in the vicinity of the inner boundary, which flings it back into the envelope where it locally accumulates. The interaction of the binary forcing frequency and the orbital angular frequency of the overdensity materializes as a modulation with beat angular frequency ω_beat = ω_b − ω_ρ.

We note a dramatic change in the latitudinal distribution from t ≃ 156 P_orb in simulation run B (Fig. 8). For t ≲ 156 P_orb, we find that mass accretion shows inclined and periodic stripes spanning all latitudes. We interpret this as an indication of the presence of a tilted lump, successively feeding the individual binary components through accretion streams. We give more details in Appendix C. In the context of CBDs, accretion onto the binary results exclusively from analogous accretion streams propagating in a low-density cavity encompassing the central binary. For t ≳ 156 P_orb, the mass and angular momentum accretion becomes more isotropic, suggesting the absence of such well-structured latitudinally extended and tilted lump. Simultaneously, the complexity of the variability increases.

Interestingly, we see that while the ω_b mode appears in all three latitudinal regions in both runs and in both regimes of run B, such peak is not present in the power spectral density of the total mass accretion rate in the first regime of simulation run B. We can explain this by the asynchronocity of mass accretion between colatitudes, which results from the migrating accretion stream and which suggests that local latitudinal analysis is necessary when studying short-term evolution of accretion in CEE. The presence of a power spectral density peak at ω_ρ in the three latitudinal regions in the first phase of simulation run B suggests that there is a large latitudinal extent of an overdense region amplifying the accretion. However, this peak frequency is not present for the total mass accretion rate in simulation run A in the same time interval. This difference is likely due to the eccentric structure of overdensities above and below the orbital plane (see Sect. 3.5), which splits ω_ρ and its harmonics about their original value. Similar phenomenon was identified in CBD simulations (e.g., Shi et al. 2012; Noble et al. 2012; D’Orazio et al. 2013). In simulation run A, we can identify two peaks at around 0.18 Ω_orb and 0.28 Ω_orb, which correspond to the splitting of Δω ≃ 0.05 about an unsplit lump angular frequency ω_ρ ≃ 0.23. For the second phase of simulation run B, additional peaks appear for 2π/3 ≤ θ ≤ π and π/3 ≤ θ ≤ 2π/3 at around 0.13 Ω_orb and 0.255 Ω_orb, which correspond to an angular frequency splitting of Δω ≃ 0.0625 about an unsplit lump angular frequency ω_ρ ≃ 0.1925.

3.4. The lump

Although we saw signatures of the lump in the power spectra, the density snapshots in Fig. 2 do not make the existence of a lump glaring. To better visualize the lump and to assess its potential effects on the inner envelope dynamics and accretion onto the central binary, in Fig. 11 we examine the space-time evolution of quantity A₁(r, t) (e.g., Roedig et al. 2011, 2012; Shi et al. 2012; Noble et al. 2012; Lopez Armengol et al. 2021), which is the θ-integrated m = 1 mode of the Fourier transform of the density with respect to the azimuth φ,

(34)

Fig. 11. Space-time diagram of the m = 1 mode of the Fourier transform of the θ-integrated density with respect to the azimuth φ for model A. Panel b is obtained with a much larger time resolution than panel a, such that high frequency fluctuations are well resolved. On panel b, we see the overdensity generated at early time (t ≃ 147 Porb) propagating outward and expanding radially, and a new lump building up from t ≃ 152 Porb.

First, we see a high-frequency variation of A₁ in the inner envelope (r ≲ 0.8), which is caused by the forcing angular frequency, ω_b. A small fraction of such overdensities contribute to the increase of mass and angular momentum accretion shown in Figs. 7 and 8 while the rest of the material is strongly gravitationaly torqued by the binary and is flung back into the envelope. These outflowing streams collide and accumulate in a large range of colatitudes starting from r ≃ 0.8. The resulting overdense region dilutes and propagates radially far into the envelope, as we can see from the outward propagating overdensity in Fig. 11. The inner part of this overdense lump feeds the inner envelope, but eventually the lump propagates far enough into the envelope that it no longer interacts with the inner region and a new lump begins to form again. We illustrate this process in Fig. 12, where we show the surface density averaged in the z direction for a thin region of opening angles ±π/8 about the orbital plane,

(35)

Fig. 12. Surface density about the orbital plane (35) for model A (top) and A′ (bottom) at t = 158 Porb. We see the outward propagation of a lump and the formation of a new one for model A, and the absence of a structured lump in model A′ resulting from the absence of accretion streams. Black dots indicate the position of the two cores.

While overdensities we see in our simulations are in many aspects remarkably similar to the lump present in CBDs simulations, they exhibit fundamental differences. In CBD simulations, ω_ρ is the orbiting frequency of a single lump that is fed by accreting material flung back into the envelope and that typically remains near the cavity edge. In contrast, in our CEE simulations, ω_ρ characterizes the formation frequency of nonaxisymmetric overdensities that propagate far into the envelope.

3.5. Eccentricity growth and evolution

Throughout this work, we assume that the binary orbit remains circular and fixed. As a result, we can only speculate about implications of our findings for binary eccentricity (see Sect. 4.2), but we can directly study the related eccentricity of the envelope. We would expect that the initially noneccentric envelope encompassing an equal-mass binary on a circular orbit will not become eccentric. However, the frequency splitting of ω_ρ observed in Figs. 9 and 10 suggests that envelope eccentricity develops in our simulations. To illustrate this more thoroughly, in Fig. 13 we show the space-time diagram of shell-averaged envelope eccentricity

(36)

Fig. 13. Space-time diagram of shell-averaged envelope eccentricity. The meaning of symbols is the same as in Fig. 11.

Similarly to A₁, e is subject to a high frequency variation in the inner envelope (r ≲ 1) according to the forcing frequency ω_b. The dynamics of accretion and of the lump is tightly linked to the generation and propagation of eccentricity in the envelope. Eccentricity is excited by the amplification of small asymmetries in the interaction between accretion flows and the central binary either by stream impact on the inner boundary (e.g., Shi et al. 2012) or by resonant Lindblad excitation (e.g., Lubow 1991a,b; Papaloizou et al. 2001; Muñoz & Lithwick 2020). We see that while a fraction of the newly generated eccentricity is contained in the colliding outflowing streams forming lumps, the rest is trapped in between successively created lumps, where e grows over time. Consequently, as the lumps propagate outward, the eccentricity follows.

In Fig. 14, we show the evolution of the mean envelope eccentricity within our simulation domain, e_env, which is defined as

(37)

Fig. 14. Evolution of the mean envelope eccentricity within the numerical domain. The black dashed line shows a linear fit yielding a growth rate λenv ≃ 0.022 Ωorb.

We find that similarly to CBDs (e.g., Shi et al. 2012), the e_env initially increases very rapidly in response to the quadrupole perturbation associated with the replacement of the time- and latitude-averaged binary potential with the true expression. For 25 ≲ t/P_orb ≲ 50, eccentricity grows exponentially with a growth rate λ_env ≃ 0.022 Ω_orb. The growth rate does not depend on the initial angular momentum of the envelope nor on the presence or absence of viscosity. Interestingly, the value of λ_env is of the same order as the eccentricity saturation growth rate of ∼0.018 Ω_orb obtained by Shi et al. (2012) in the context of CBDs, which could suggest a common physical origin. After t ≃ 50 P_orb, the exponential growth saturates and e_env reaches a statistically stationary state with a mean value e_env, f ≃ 0.12, which is independent of β and of the presence or absence of viscosity. Such eccentricity saturation likely results from nonlinear effects, which suggests that eccentricity excitation and damping reach a quasi-equilibrium that may be maintained throughout the entire post-dynamical spiral-in phase (e.g., Shi et al. 2012; Teyssandier & Ogilvie 2016; Miranda et al. 2017; Muñoz & Lithwick 2020).

3.6. Convective stability and angular momentum transport in the envelope

During the post-dynamical CEE, the central binary interacts with the surrounding gas and a complex interplay between the torques and internal stresses continuously injects, removes, and redistributes angular momentum within the envelope. In this section, we investigate the stability of the envelope and analyse its dynamics by characterizing the various angular momentum transport processes.

3.6.1. Solberg–Høiland criterion for convective stability

In the absence of viscosity, thermal diffusion, and radiation pressure, Solberg (1936) and Høiland (1941) proposed the following necessary but not sufficient condition for convective stability, which for a stratified and rotating fluid with Ω = Ωe_z reads

(38)

where

(39)

ℓ_z = s²Ω is the specific angular momentum, S is the specific entropy, and c_p is the heat capacity at constant pressure. In Fig. 15, we show the Solberg–Høiland criterion for convective stability in the xy and xz planes soon after replacing the averaged binary potential with its full expression and at late-time. We see that as soon as they are present, gravitational perturbations from the central binary destabilize the flow according to the Solberg–Høiland criterion. In practice, this translates into small scale turbulent mixing between spiral arms (see Fig. 2, first row), which is initially not strong enough to destroy the spiral structure. As the envelope expands and the stabilizing effect of density stratification is reduced, the vertical size of the turbulent eddies increases and the spiral structure is partially destroyed. We observe behavior resembling the ab initio simulation of dynamical plunge-in from Ohlmann et al. (2016), where the theoretically stable and unstable layers alternate in a geometrically thick disk-like structure about the orbital plane. The radial spatial frequency decreases outward as the stabilizing effect of stratification becomes weaker.

Fig. 15. First Solberg–Høiland criterion for convective stability for run B at t = 20 Porb (top row), and t = 250 Porb (bottom row). Negative values indicate convective instability according to Eq. (38).

3.6.2. Local torque balance

While the various volume-integrated torques presented in Sect. 3.2.1 trace the evolution of the total angular momentum reservoir of the common envelope in our numerical domain, it is also important to examine the spatial variation of such torques. As we show in detail in Appendix A, the local angular momentum transfer rate across the common envelope reads

(40)

where

(41)

(42)

(43)

In Fig. 16, we show the contributions to as a function of r for simulation runs A, A′, B, and C. We see that for accreting models A, B, and C, gravitational torque only plays a minor role in the redistribution of angular momentum in the inner envelope and essentially no role far from the central binary. This occurs because the density is globally a decreasing function of r and because lim_{r → ∞}∂Φ/∂φ = 0. Instead, it is the advective torque that transports angular momentum. Up to r ≃ 1.5, the angular momentum is transported inwards , but for larger r the angular momentum flows outward. However, when accretion is prevented by reflecting boundary conditions (simulation run A′), gravitational torque dominates advective torque up to r ≃ 0.2, which results in inward angular momentum transport. At larger r, the advective torque transports angular momentum outward. The main differences in the advective torque profiles between the four models result from different contributions of turbulence and different mean-flow angular structure.

Fig. 16. Advective and gravitational contributions to the local angular momentum transfer rate for models A (full lines), model A′ (dash-dotted lines), model B (dashed lines), and model C (dotted lines). The quantities are averaged in time interval 250 ≤ t/Porb ≤ 275.

3.6.3. Turbulent transport of angular momentum in the envelope

We can use hydrodynamic mean-field theory to asses the turbulent fluxes of angular momentum (e.g., Käpylä 2019; Rüdiger 2022). Taking the azimuthal average of the angular momentum equation in the z-direction, using Reynolds decomposition to define velocity fluctuations about their averages as , and ignoring density fluctuations from its mean value, we obtain the conservation law

(44)

where is the meridional mean velocity, overlined quantities indicate azimuthal average, and and S′ are the rate-of-strain tensors of the mean flow and of the fluctuating flow defined as

(45)

(46)

The first term on the right-hand side of Eq. (44) corresponds to the advective transport by large scale meridional flow, the second represents the turbulent meridional flux, and the third term is the viscous transport. Turbulent angular momentum fluxes are often described using Reynolds stress,

(47)

This stress tensor is often separated into nondiffusive (“Λ-effect”) and diffusive contributions described by turbulent viscosity,

(48)

where is the nondiffusive part and 𝒩_ijkl is the (turbulent) viscosity tensor (e.g., Kitchatinov et al. 1994; Kitchatinov & Ruediger 1995; Käpylä et al. 2011; Rüdiger 2022). Even in the case where we do not prescribe subgrid viscosity, an effective (convective) turbulent viscosity still exists and can be derived from the expression of the turbulent viscosity tensor 𝒩_ijkl. Conversely, when we prescribe ν > 0 (simulation run D), an additional effective viscosity associated with the simulation’s intrinsic turbulence still exists. In this case, the total effective viscosity is given by the sum of the two contributions. The disentangling of the two contributions to the Reynolds stress and the measurement of the associated simulation’s intrinsic effective turbulent viscosity is however beyond the scope of this work. Instead, we focus on the total stress.

We assume that the turbulent velocity u′, which is the deviation from the mean velocity, is the nonaxisymmetric component of the fluid flow velocity. However, because the large scale flow resulting from the gravitational torque exerted by the binary orbit is itself nonaxisymmetric, the contribution of turbulence to angular momentum transport is likely overestimated, especially in the close vicinity of the binary. Unfortunately, there is no straightforward way to establish what the mean flow is in our simulations. This is an issue also in the context of accretion and CBDs, where Hawley (2000), Hawley & Krolik (2001), Shi et al. (2012), and Lopez Armengol et al. (2021) use departure from density weighted shell average to compute velocity perturbations. Still, one could extract the actual turbulent flow with reasonable accuracy by filtering out the large scale flow using Fourier and inverse transforms (e.g., Käpylä et al. 2011). This is however beyond the scope of this work and we refer to the nonaxisymmetric perturbation u′ as the turbulent fluid flow velocity, though one has to keep in mind that this may be inaccurate in the binary close vicinity.

In Fig. 17, we show the meridional components of the azimuthally averaged advective fluxes of the total angular momentum . In Figs. 18 and 19, we show the mean and turbulent flow contributions, and . We see that for accreting models, the mean axisymmetric flow results in outward angular momentum advective flux in a equatorial disk-like structure with an opening angle that is smaller for higher initial angular momentum. Outside of the disk-like structure, the angular momentum advective flux points inward. The morphology of the radial turbulent transport of angular momentum is more complicated, because it changes sign in both cylindrical directions s and z. Such disk-like structure is not present in our nonaccreting model A′. Indeed, the inward flow is deflected by the inner boundary and any polar mass flux asymmetry between northern and southern hemispheres, however small, is amplified and can even lead to a polar outflow in one of the hemispheres.

Fig. 17. Azimuthally averaged advective radial and latitudinal fluxes of angular momentum, and averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).

Fig. 18. Mean flow contribution to the azimuthally averaged advective angular momentum radial and latitudinal fluxes averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).

Fig. 19. Turbulent flow contribution to the azimuthally averaged advective angular momentum radial and latitudinal fluxes averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).

In Fig. 20, we show the radial profile of the mean-flow and turbulent contributions to the angular momentum transfer in terms of Reynolds stress,

(49)

(50)

Fig. 20. Mean and turbulent components of the advective contribution to local angular momentum transfer rate for models A (full lines), A′ (dash-dotted lines), B (dashed lines), and C (dotted lines). The quantities are averaged over a time interval 250 ≤ t/Porb ≤ 275. Neglecting density perturbations, .

We find that for our three inviscid and accreting runs, the net radial angular momentum transport is essentially dominated by the contribution from the mean flow, which is directed inward for r ≲ 1.5 and outward for r ≳ 1.5.

The mean axisymmetric flow also leads to angular momentum advective transport in the θ direction. Specifically, such mean flow advects angular momentum toward the orbital plane in the inner part of the envelope and away from it further out. Conversely, turbulent flow advects angular momentum away from the midplane in the close vicinity of the binary and toward it in the rest of the envelope. Overall, the structure of the total angular momentum flux follows the mean flow contribution, where the angular momentum is advected toward the orbital plane in the inner envelope and away from the orbital plane far from the binary.

3.6.4. Vertical eddy scales

Since we are interested in the ability of turbulent structures to transport angular momentum radially in the envelope, we aim to estimate the typical vertical scale of turbulent convective eddies exchanging angular momentum with one another. To make sure that we properly isolate turbulent flow, we focus on its latitudinal component in the orbital plane. This is because the nonaxisymmetric contribution of the large-scale mean-flow, which pollutes the inferred turbulent velocity, results from the envelope’s response to the gravitational perturbations exerted by the binary and is zero in the θ-direction in the orbital plane. Let us first introduce the normalized auto-correlation of the turbulent latitudinal velocity on the orbital plane, which we azimuthally average and we integrate over an arbitrary radial domain [r_min, r_max],

(51)

We also introduce the integral radial scale (e.g., Townsend 1976; O’neill et al. 2004; Mora & Obligado 2020),

(52)

Because the low amplitude tail of Eq. (51) may contain information extraneous to turbulent motion, we integrate Eq. (52) up to the first zero-crossing of R_θθ(r′,t), as is commonly done in experimental and numerical fluid dynamics (O’neill et al. 2004). In homogeneous turbulence, Λ can be interpreted as the typical radial scale of the energy-containing turbulent eddies. Because the mean density stratification becomes weaker as the distance from the central binary increases, turbulence is not homogeneous in our simulations, and the interpretation of Λ is more ambiguous. Here, Λ represents the weighted average of all turbulent radial scales in the flow with a dominant contribution from eddies containing higher energy. In an effort to mitigate this ambiguity, we integrate the auto-correlation of the azimuthally averaged turbulent latitudinal velocity on the orbital plane in three arbitrary regions of strong, moderate, and weak stratification labeled I (r ≤ R_domain/10), II (R_domain/10 < r ≤ R_domain/2), and III (r > R_domain/2). We interpret the resulting integral scale as the typical eddy scale in each region.

In Fig. 21, we show the time evolution of Λ in the three regions for our three inviscid simulation runs. In Table 2, we provide the time-averaged values, Λ^{t_min − t_max}, where t_min and t_max denote the time interval of averaging in the units of P_orb. We see that in all three models the vertical extent of turbulent eddies decreases outward as the stabilizing effect of density stratification decreases. However, depending on the amount of angular momentum present in the envelope at the onset of the post-dynamical inspiral phase (parameter β), the radial dependence of turbulent eddy vertical scales varies substantially. The more angular momentum the secondary star injects into the shared envelope during dynamical inspiral, the more the envelope gets deformed by centrifugal forces. Hence, the injected angular momentum modifies density stratification by expanding the envelope anisotropically and affects envelope stability through the sign of angular momentum gradient (Eq. (38)). As a result, we find that envelopes with lower angular momentum content at the onset of the post-dynamical phase end up being less effectively stratified, which leads to a reduction of their ability to limit the vertical extent of turbulent eddies. However, it is important to note that in region I, for both simulation runs A and B while for run C. The fact that we observe similar eddy scales for runs A and B with different initial total angular momentum in this inner region is due to our initial spin-up setup where the innermost layers are spun-up to critical rotation, while they remain subcritical during spin-up for run C. Additionally, we further note that far away from the binary in region III, simulation runs B and C yield similar vertical eddy scale. This roughly constant value likely constitutes a limit at very low stratification that may depend on the numerical size of the domain (see also Garaud et al. 2017). This limit formally implies that the assumption of a vertical eddy scale that is proportional to the local pressure scale height (e.g., Vitense 1953; Zahn 1989) may fail in the limit of low stratification as the effective viscosity would locally tend to infinity.

Fig. 21. Time evolution of the integral radial scales Λ (Eq. (52)) in three regions I (r ≤ Rdomain/10), II (Rdomain/10 < r ≤ Rdomain/2), and III (r > Rdomain/2) for simulation runs A (full lines), B (dashed lines), and C (dotted lines). The discontinuity of ΛIII for run A is a result of not changing sign on the same radial scales for all φ, that is Rθθ(r′,t) does not cross zero during a few orbital periods.

3.6.5. The role of viscosity

In simulation run D, we prescribed an isotropic viscosity whose effects on the envelope dynamics add to those from the effective viscosity associated with Reynolds stresses. In Fig. 22, we show the mean, turbulent, and viscous contributions to local angular momentum transfer rate across the common envelope for model D and in Fig. 23 we show the mean flow, Reynolds stress, and viscous contributions to the azimuthally averaged advective radial and latitudinal angular momentum fluxes. We see that for our value of α_ν, the azimuthally averaged total angular momentum viscous flux is essentially directed outward, and only plays a minor role in the radial and latitudinal transport of angular momentum in the shared envelope. Hence, simulation run D retains the same angular momentum transport features as our inviscid runs. This result reflects the efficiency of the transport by the mean flow and by the Reynolds stress transport and suggests that the effective viscosity associated with the latter is likely several orders of magnitude larger than ν. Still, despite being small, the prescribed viscosity has a stabilizing effect on the shear flow in the inner envelope, and the ratio between turbulent and mean components of the advective angular momentum transfer rate is, at a given time, smaller in simulation run D. Such stabilizing effect of viscosity also locally delay the onset of turbulence.

Fig. 22. Mean, turbulent, gravitational, and viscous contributions to local angular momentum transfer rate across the common envelope for model D, averaged in time from t = 140 Porb to t = 150 Porb.

4. Discussions

4.1. Comparison with CBDs

In this work, we draw many analogies between simulations of CBDs and post-dynamical stages of CEE and it is thus instructive to briefly discuss differences and similarities between these situations. First, there are significant differences in the origin of the gas surrounding the binary. CBDs can often occur as remnants of star formation out of molecular clouds or they are thought to accompany orbiting super-massive black holes (Begelman et al. 1980; Bate & Bonnell 1997; Milosavljević & Phinney 2005; Matsumoto et al. 2019). In this situation, the density distribution and angular momentum content of the disk depends not only on the properties of the binary, but perhaps more significantly on the accretion for larger distances. Conversely, CEE can often be regarded as an isolated object, where the formation process inextricably links together the distribution of density, energy, and angular momentum in the envelope with the properties of the central binary. Second, CBDs are often observed and simulated as relatively optically and geometrically thin, 2D objects. Instead, the post-dynamical shared envelope in CEE contains large amount of mass, which prevents cooling and keeps the geometry strictly three-dimensional. Ultimately, the shared envelope disperses and whatever gas remains should cool to a thin disks, as is observed in post-AGB binaries (e.g., Dermine et al. 2013; Kluska et al. 2022). The transition between these two regimes of post-dynamical CEE should be a subject of future study. Finally, thermal convection is weak or even absent in CBDs, where it is the magnetorotational instability (MRI, Balbus & Hawley 1991) that is instead often recognized as the main source of turbulence (e.g., Cabot 1996; Stone & Balbus 1996; Balbus & Hawley 1998) modeled with a turbulent effective viscosity using the α ansatz (Shakura & Sunyaev 1973). Conversely, common envelopes are expected to be vigorously convective, making thermal convection inevitable (e.g., Soker 1993; Ohlmann et al. 2016; Sabach et al. 2017; Grichener et al. 2018; Wilson & Nordhaus 2019).

Despite these fundamental differences, we have shown that there are similarities and even commonalities between these two systems. We found that mass and angular momentum accretion onto the central binary (when allowed) has the same temporal variability with two characteristic frequencies. The first frequency is associated with the quadrupolar moment contribution to the binary potential, while the second one with the formation and propagation of overdensities, which share many common characteristics with the lump located near the cavity edge in CBD simulations. Because of the complicated geometry of common envelopes, we have shown that a local analysis of the accretion flux is necessary to understand its short-term variability. The behavior of the orbital separation evolution is dictated by the same condition for the two problems, specifically, j < 3/8 gives orbital contraction when q = 1 and . This condition suggests predominant orbital contraction in CEE simulations, while orbital expansion is possible for a wide orbital parameter range in CBD simulations (e.g., Miranda et al. 2017; Muñoz et al. 2019). Finally, we found that the shared envelope develops eccentricity, which grows with an exponential growth rate that is of the same order as that obtained by Shi et al. (2012) in the context of CBDs and which saturates to reach a statistically stationary value as is also seen in CBD simulations (e.g., Miranda et al. 2017; Muñoz & Lithwick 2020). To summarize, the abundant literature and ongoing work in the field of accretion and CBDs can be of precious help to better understand the post-dynamical inspiral phase of CEE.

4.2. Implications for CEE

Our findings have a number of implications for CEE. We find that the orbital evolution of the central binary does not stall even when the gas in the immediate vicinity of the binary corotates. Instead of the commonly assumed drag, the binary transfers angular momentum to the envelope by generating spiral waves and turbulence. We find that the associated timescale of orbital contraction is τ_b = |a_b/ȧ_b| ∼ 10³ to 10⁴ orbital periods of the binary when accretion is allowed or slowly decreases to τ_b ∼ 10⁵ P_orb at t ≃ 450 P_orb when accretion is prevented. These timescales are similar to what is typically found in CBDs (e.g. Artymowicz et al. 1991). We emphasize that this timescale refers to the inner binary orbit, which is much smaller than the outer extent of the envelope. Our results suggest that while there is gas in the shared envelope, the binary should continue to spiral in due to nonlocal interactions with the nearby gas, albeit much slower than in the preceding dynamical plunge-in phase. Many ab initio works on CEE find that the central binary orbits continue to slowly shrink at the end of the simulations. Based on our results, we suggest that the orbital contraction rate should not gradually approach zero but always remains at a rather small but finite value. We also suggest that achieving complete envelope ejection and final orbital separations compatible with expectations could be possible simply by following the evolution for much longer time than what is currently done.

No thermal energy transfer through the envelope is required to reduce the orbital separation. However, at some point in time, the energy diffusion timescale through the envelope should become comparable to τ_b. When that happens, it is possible that thermal coupling between the binary and the envelope is established, which might affect the orbital decay. For example, efficient removal of energy deposited in the vicinity of the binary could increase the envelope density around the binary, which would lead to higher torques and smaller τ_b. In order for a thermal “self-regulating” process to have a chance to accelerate the orbital decay, the thermal timescale of the envelope of the primary star, t_KH, has to be shorter, or of the order of τ_b.

In order to evaluate t_KH/τ_b for different stars, we need to make three approximations. First, we make use of the fact that τ_b/P_orb ∼ 10³ to 10⁵ and we study instead quantity t_KH/P_orb. Second, P_orb is evaluated assuming that parameters of our simulation described in Sect. 2, a_b/R = 0.16 and M₁/M = 0.2, are applied uniformly to all progenitor primary stars. Third, stellar quantities R, M, and L represent the values of the stellar model before the binary interaction. In Fig. 24, we show the ratio t_KH/P_orb evaluated along evolutionary tracks of four single stars with masses 1, 5, 12, and 20 M_⊙. We see that for low-mass red giants without fully developed convective envelopes, R ≲ 100 R_⊙, the ratio is t_KH/P_orb ≳ 10⁵, which implies that the envelope would not have enough time to thermally couple to the inspiralled binary before its orbit significantly decays. For these stars, the time window for any additional processes acting to remove the envelope on long timescales might be severely restricted, because the lifetime of the envelope is set by the fast orbital decay timescale rather than the thermal timescale. An example of such possible long-lasting processes are strong pulsations or dust-driven wind (Clayton et al. 2017; Glanz & Perets 2018). For AGB and high-mass stars, the thermal coupling between the inspiralled binary and the envelope seems more likely. However, unless t_KH/P_orb ≲ 10³, which occurs only near the maximum expansion of the stars, the orbital decay could still remain unaffected by the thermal processes if the binary is able to efficiently accrete from the envelope due to a “pressure valve” such as launching of jets (Soker & Livio 1994; Chamandy et al. 2018; Shiber et al. 2019).

Naturally, our simple estimates in Fig. 24 have a number of caveats. For example, the binary might dynamically plunge-in to much lower values of a_b/R, as was seen in several recent simulations (Ohlmann et al. 2016; Lau et al. 2022a), or the binary might relatively quickly shrink its orbit early in the post-dynamical phase, as Fig. 6 suggests. Furthermore, CEE is often preceded by strong thermal-timescale mass transfer, which leads to a significant decrease of primary’s luminosity due to thermal restructuring of the envelope. The luminosity could decrease by a up to a factor of 10, which would enlarge t_KH (e.g., Blagorodnova et al. 2021). All of these effects would tend to increase t_KH/P_orb and make thermal timescale influences less likely. Conversely, as the envelope expands, t_KH and the diffusion timescale decrease. Since there are number of effects working in opposite directions, the overall importance of thermal effects on the inspiral is not immediately obvious. This long-term evolution cannot be easily studied by direct multidimensional simulations, but some insight can be obtained with 1D models with prescriptions calibrated to include physical processes studied in this work. In the future, we aim to enlarge our grid of parameters to make such parameterization possible.

In the light of our results, it is also interesting to discuss the energy-conserving α_CEE formalism that is commonly used to predict CEE outcomes. Our results suggest that much of the orbital decay during the post-dynamical phase might be over before thermal effects become important, which lends support to the energy-conserving formalism. At the same time, the orbital decay is clearly separated into two regimes: a fast dynamical plunge-in and much slower post-dynamical inspiral. Even if energy is conserved in both of them, the value of the α_CEE parameter might be different, because both types of orbital decay depend differently on binary properties and envelope structure. Furthermore, the post-dynamical inspiral depends on the efficiency of accretion, which could be influenced, among other effects, by jets. All of this could lead to different effective values of α_CEE for different populations of binary stars, which is not surprising, but perhaps also to a spread of α_CEE among a single population. In any case, our results generally motivate the development of two-step or multistep CEE formalisms (e.g., Hirai & Mandel 2022).

Although one of our original motivations for this work was to see whether the binary could reaccrete some of its angular momentum and expand its orbit, our results suggest that this is unlikely. The shared envelope is very thick and accretion near the polar regions brings in gas with very low specific angular momentum. The situation could change at later phases when the remaining envelope is able to cool to a thinner disk. Investigating this transition should be a subject of future study.

Finally, in our work we have made the assumption that the binary orbital motion has completely circularized after the dynamical plunge-in. This is in agreement with 3D hydrodynamical simulations, which typically find quasi-circular orbits at the end of this phase provided that the initial eccentricity is low, (e.g., Ricker & Taam 2012; Passy et al. 2012; Ohlmann et al. 2016; Glanz & Perets 2021). However, this does not necessarily imply that the orbit remains circular throughout the post-dynamical inspiral phase. A variety of physical processes can lead to the growth or decrease of binary eccentricity, such as accretion streams impact on binary components or the gravitational interaction between the central binary and its nonaxisymmetric eccentric envelope such as the one we find in Sect. 3.5. While such phenomena can potentially lead to the binary eccentricity growth, binary eccentricity may generate new resonances that can in turn damp eccentricity (e.g., Lubow 1991a,b). Orbital eccentricity may therefore be nonzero during the post-dynamical spiral-in phase, and stabilize or oscillate about a fixed value, similarly to what is seen in CBDs (e.g., Roedig et al. 2011; Zrake et al. 2021), and could help to explain nonzero eccentricities seen in some post-AGB and post-CEE binaries (e.g., Dermine et al. 2013; Kruckow et al. 2021).

5. Conclusions

In this work, we performed a series of 3D hydrodynamic numerical simulations of the post-dynamical inspiral phase of CEE. We used the procedure of Morris & Podsiadlowski (2006, 2007, 2009) to mimic the outcome of the preceding dynamical plunge-in and to establish controlled initial conditions for our simulations (Fig. 2). Our first aim was to determine the timescale of binary separation evolution in response to the various torques acting on the system when accretion is turned on or off. We have computed the various torques acting on the binary and we found that they always result in the contraction of the orbit, regardless of whether accretion is allowed or not (Fig. 6, Sect. 3.2). When accretion is allowed, mass and angular momentum accretion drive the orbital contraction and the orbital contraction timescale rapidly reaches a quasi-steady value of 𝒪(10³ − 10⁴ P_orb). Without accretion, orbital contraction is solely driven by the gravitational torque. Because of the envelope expansion, the amplitude of the gravitational torque slowly decreases, leading to a slow increase of the orbital contraction timescale. After 450 P_orb, this timescale reaches a value of 𝒪(10⁵ P_orb).

Our results imply that while the binary is embedded in gas, the orbit contracts even if the gas immediately surrounding the binary is corotating. The orbital decay timescale is much slower than what is seen during the dynamical plunge-in. This suggests a significant reduction of orbital separation and more efficient envelope ejection is possible even after the dynamical plunge-in and that current simulations have not been carried out over a sufficiently long period to observe this effect. Since the orbital separation is very small compared to the outer extent of the envelope, the post-dynamical inspiral timescale is much shorter than the thermal timescale for primary stars with radius ≲100 R_⊙. Even for larger stars, the post-dynamical decay does not have to be significantly influenced by the thermal response of the envelope and thermal “self-regulation” does not seem to be unavoidable, but this is contingent on the dynamics very close to the binary such as the presence or absence of accretion or jets. The short inspiral timescales lend support to adiabatic treatment of CEE, but motivate viewing CEE as an (at least) two-step process.

Our second aim was to find the typical frequencies associated with the short-term variability of mass accretion onto the binary and to compare the results to CBDs. We found that the main features of mass accretion variability in the context of post-dynamical CEE are similar to that of CBDs. Specifically, the variability is connected to the forcing angular frequency of the quadrupolar moment contribution to the binary potential ω_b = 2 Ω_orb and to the frequency associated with the formation of nonaxisymmetric overdensities in the inner part of the envelope ω_ρ = Ω_orb/5. Such overdensities result from accreting material being flung back into the envelope, and accumulating at a distance of roughly six binary separations from the binary center of mass (Figs. 7–11, Sect. 3.3). We found that the resulting overdense and eccentric “lump” then propagates far into the envelope, contrary to the case of CBDs, and feeds and enhances mass accretion (Fig. 12 and Sect. 3.4). Because of the spherical shape of the CEE problem contrasting with flat CBDs, such frequencies do not necessarily characterize the global mass accretion rate. Instead, the presence or absence of latitudinally migrating accretion streams leads to the synchronicity or asynchronocity of the mass flux at all colatitudes. When asynchronous, time variability of the latitudinal integrated mass flux may be smoothed out. A local analysis of the accretion flux is therefore necessary to understand its short-term variability. Finally, we found that envelope eccentricity is excited in the vicinity of the binary and propagates outward within and in-between successive lumps. During this process, the eccentricity amplifies and builds up in the envelope, which leads to the splitting of the accretion frequencies (Fig. 11 and Sect. 3.5). The envelope mean eccentricity grows exponentially with a growth rate λ_e ≃ 0.022 Ω_orb until it reaches a statistically stationary value of e_env, f ≃ 0.12 independent of the initial angular momentum of the envelope and of the presence or absence of viscosity (Fig. 14). Growth of eccentricity in the envelope could in turn excite eccentricity in the inner binary.

Our third aim was to understand how angular momentum is transported within the envelope (Sect. 3.6). We showed that, similarly to previous dynamical plunge-in simulations, gravitational perturbations from the orbiting binary during post-dynamical in-spiral phase trigger the destabilization of the envelope. This destabilization results in turbulent convection contributing to the transport of energy and angular momentum throughout the shared envelope (Fig. 15 and Sect. 3.6.1). However, we showed that this contribution is rather small. Instead, the angular momentum flux is dominated by large scale axisymmetric fluid flows, which consist of an outward transport in a geometrically thick disk-like structure about the orbital plane, except in the close vicinity of the central binary where the flux points inward. However, because the net (latitudinally integrated) radial angular momentum transport by the mean flow at a given radius is rather weak, turbulent transport measured by Reynolds stresses in fact is also important and can have an effect of the same order as that of the mean flow on the global transport of angular momentum. In particular, we showed that Reynolds stresses can locally strongly damp or enhance the outward transport of angular momentum by the mean flow (Figs. 16–20).

Our final aim was to characterize the role of viscosity originating both from turbulent motions and from unspecified processes acting on subgrid scales. We showed that because the stabilizing effect of density stratification decreases outward, vertical convective eddy scales increase with radial distance from the central binary. We showed that envelopes with higher initial angular momentum content are more strongly stratified and thus contain smaller convective cells, however, their typical size is limited to a fraction of the domain radius in the limit of low stratification, suggesting that α-type viscosity models would fail in outer layers (Figs. 21–23). We further found that prescribing a background kinematic viscosity with α_ν = 10⁻³ does not significantly affect the binary separation evolution, nor the transport of angular momentum within the shared envelope (Sect. 3.6.5).

Fig. 23. Mean flow (a), turbulent flow (b), and viscous (c) contributions to the azimuthally averaged advective radial and latitudinal angular momentum fluxes for the viscous model D. We averaged the quantities over ten orbital periods for 140 ≤ t/Porb ≤ 150.

Fig. 24. Ratio of thermal timescale tKH to the orbital period of the binary Porb inside the shared envelope in the Hertzprung–Russel diagram constructed for four solar-metallicity nonrotating evolutionary tracks from the MIST database (Dotter 2016; Choi et al. 2016). Here, tKH = G(M1 + Menv)2/(2RL), where L is the luminosity of the star.

Our new way of studying late stages of CEE has its limitations. For example, it is important to keep in mind that we have excised a central region encompassing the binary and that we only considered two extreme regimes of accretion (maximum or none). We have further considered fixed binary orbit, assuming that the contraction or expansion timescale is much longer than the duration of our simulation. Our estimates of the orbital contraction timescale suggest that such assumption might not be completely valid when accretion is allowed. This limitation cannot be easily lifted: if we allowed the orbit to shrink or expand we would have to change the position of the inner boundary with time in order to not discard the flow dynamics in the vicinity of the binary or to prevent the binary from entering the numerical domain. This would affect the conservations of mass and angular momentum and the numerical cost of our simulations. Another solution would be to use Cartesian grid, but in such a case we would lose the advantages of spherical geometry. We have further assumed that binary eccentricity remains zero throughout our simulations. However, a variety of physical processes may lead to the binary eccentricity growth or decrease, although they may balance each other out. Finally, we do not include gas self gravity which could affect the binary-envelope interaction.

More sophisticated initial parameters could complicate the our results and need to be explored in future works. We plan to investigate eccentric binary orbits, binaries with mass ratios different from unity, more sophisticated inner boundary conditions, or the effect of magnetic fields. Changing these parameter could dramatically impact the binary-envelope interaction, resulting in very different variability, amplitude, and angular distribution of mass and angular momentum accretion onto the binary, binary separation evolution, and angular momentum transport within the envelope.

Acknowledgments

We thank the anonymous referee for comments that improved this paper. We thank Kengo Tomida for discussions about Athena++. The research of D.G. and O.P. has been supported by Horizon 2020 ERC Starting Grant ‘Cat-In-hAT’ (grant agreement no. 803158). This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90140). O.P. thanks the KITP program “Bridging the Gap: Accretion and Orbital Evolution in Stellar and Black Hole Binaries” for hospitality and inspiration: this research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

References

Appendix A: Angular momentum conservation

We take the cross product of r with the momentum Eq. (1), multiply by e_z and integrate it over the domain’s volume,

(A.1)

Using the fact the viscous stress tensor T is symmetric, and writing r = x_je_j,

(A.2)

where ϵ_ijk is the Levi-Civita tensor of rank three, and

(A.3)

Eq. (A.1) can be rewritten as

(A.4)

Finally, using Gauss divergence theorem, Eq. (A.4) can be rewritten as

(A.5)

where n_⊥ is is the outward-pointing unit vector at the boundaries’ surface, is the time derivative of the z−component of the gas angular momentum, and , T_z, grav, and are the advective, gravitational, and viscous torques on the system, respectively. Replacing the viscous stress tensor by its expression finally yields,

(A.6)

where

(A.7)

Similarly, local torques are obtained by taking the cross product of r with the momentum Eq. (1), multiplying it by e_z, taking the radial derivative, integrating the result over the volume comprised between r and r_out, and using the fact that the gravitational torque vanishes at the outer boundary, that is ∂Φ/∂φ tends to zero far from the central binary (see e.g., Miranda et al. 2017).

Appendix B: Effect of boundary conditions on the gravitational torque

Fig. 6(b) shows that the amplitude of the gravitational torque largely depends on the enforced inner boundary conditions. In particular, while rapidly settles to statistically (almost) zero when the inner boundary is open to mass and angular momentum flow toward the binary, decreases much slower and thus drives the orbital contraction of the binary when accretion is turned off. To understand this fundamental difference, let us first investigate the conditions for . For the gravitational torque to be zero, mass distribution must satisfy the symmetry property of |∂Φ/∂φ| with regard to the plane 𝒫 defined as

(B.1)

that is the plane orthogonal to the orbital plane and intersecting the binary semi-major axis. In our simulations, the gravitational torque arises from the quadrupolar moment of the binary potential generating a two-armed spiral density wave and breaking the mass distribution symmetry with respect to the plane 𝒫 (see Fig. B.1, left). We decompose the amplitude of the gravitational torque as

(B.2)

Fig. B.1. Specific gravitational torque for simulation run A at t ≃ 25 Porb (left) and t ≃ 250 Porb (middle), and for run A’ at t ≃ 250 Porb (right). Black crosses indicate the location of the two cores, and the dashed black line indicates the intersection of 𝒫 with the orbital plane.

where ∫|∂Φ/∂φ|dV is constant in our simulations,

(B.3)

Here, 0 ≤ ϵ ≤ 1 quantifies the asymmetry of the mass distribution with respect to 𝒫, and ρ_eff is the weighted-averaged density with a dominating contribution from the gas subject to larger |∂Φ/∂φ|.

Fig. B.2 shows the time evolution of ⟨ϵ⟩_P and ⟨ρ_eff⟩_P for models A and A′, where ⟨ ⋅ ⟩_P indicates a temporal smoothing over one orbital period. We see that closing the inner boundary to mass and angular momentum flow toward the binary typically leads to larger ϵ and ρ_eff. Turbulence-induced alteration of the asymmetric spiral density wave structure likely explains the difference in ϵ between the two models, simulation run A being more turbulent in the vicinity of the binary (see Fig. B.1). Reflecting boundary conditions indeed lead to the accumulation of material at the inner boundary leading to both larger ρ_eff and increased stabilizing density stratification. However, because the energy transferred from the binary orbit to the envelope as kinetic and internal energy is higher for higher ρ_eff, the envelope expands more rapidly in simulations run A’. This leads to a more rapid decrease of ρ_eff which in turn leads to the reduction of the stabilizing effect of density stratification and thus to a slow decrease of ϵ. Inevitably, the gravitational torque in simulation run A’ eventually also reaches statistically (almost) zero (see Fig. 6(b)). Because the gravitational torque is exclusively responsible for the orbital contraction when the binary does not accrete, |a/ȧ| does not reach a steady value in simulation run A’, and keeps increasing.

Fig. B.2. (a) quantifying the asymmetry of the mass distribution with respect to the plane 𝒫 as a function of time. (b) Weighted-averaged density with a dominating contribution from the gas subject to larger |∂Φ/∂φ|, , as a function of time.

Appendix C: Migrating accretion streams

In Sect. 3.3 and Fig. 8, we discuss that for t ≲ 156 P_orb, mass accretion onto the binary exhibits inclined and periodic stripes spanning a wide range of latitudes. To understand the cause of this latitudinally migrating accretion, we show a space-time diagram of the mass flux onto the individual components of the binary for model A during this time period in Fig. C.1 (A). We see that the mass flux peaks at a maximum latitude successively on core 1 (at t = t₁) and core 2 (at t = t₂), with a period of P_orb/2 = t₂ − t₁, that is, when the two cores switch position. This suggests that the source of the high-latitude accretion orbits much slower than the binary, that is Ω ≪ Ω_orb. We show the mass flux toward and away from the individual cores, azimuthally averaged within a π/6 opening angle about the position of the cores at t = t₁ and t = t₂, respectively in Figs. C.1 (B) and (C). We see that high-latitude accretion through the inner boundary results from accretion streams originating from low-latitude regions. Such streams initially propagate radially, then toward the orbital plane when approaching the binary. Accretion streams approaching the orbital plane carry a certain amount of momentum. Their inertia allows them to travel a certain distance above the orbital plane before the material is pulled back toward it by the gravitational pull of the core. In addition, in the vicinity of the binary, the enhanced accretion streams deflect the outflow associated with spiral density waves to the north.

Fig. C.1. A: Space-time diagram of the mass flux onto the individual components of the binary for model A, azimuthally averaged in the range φ1 − π/12 ≤ φ1 ≤ φ1 + π/12 (a) and in the range φ2 − π/12 ≤ φ2 ≤ φ2 + π/12 (b). The horizontal dashed lines indicate the orbital plane θ = π/2. The two vertical dashed lines are separated by Porb/2 and indicate two successive times, t1 and t2, when the mass flux is maximum at the highest latitude. B: Mass flux toward and away from the individual cores, azimuthally averaged in the range φ1 − π/12 ≤ φ1 ≤ φ1 + π/12 (a) and in the range φ2 − π/12 ≤ φ2 ≤ φ2 + π/12 (b) as a function of colatitude and radius at t = t1. C: same as B but at t = t2. D: Surface density about the orbital plane (35) averaged over t1 ≤ t ≤ t2. The dashed circle indicates the approximate radial location of the overdensity responsible for the migrating accretion streams (r = 0.7). E: Density cross section in the θφ plane at r = 0.7 averaged over t1 ≤ t ≤ t2.

We show the surface density about the orbital plane (35) averaged over t₁ ≤ t ≤ t₂ in Fig. C.1 (D). We see that the overdensity located around r ≃ 0.7 is characterized by an angular frequency Ω_lump ≪ Ω_orb and is therefore a candidate for the source of high-latitude accretion. Finally, we show the density cross section in the θφ plane at r = 0.7 averaged over t₁ ≤ t ≤ t₂ in Fig. C.1 (E). We find that this overdensity is spatially extended over a wide range of θ. Its presence and its shape explains the migrating accretion streams shown in Fig. 8 as follows. As individual orbiting cores pass in front of the overdensity, they first pull its material from high colatitudes, enhancing accretion streams (Figs. B(a) and C(b)). Due to their high inertia, such accretion streams reach the inner boundary in the northern hemisphere. The material impacting the inner boundary is partially accreted by the binary, the rest is pulled back toward the orbital plane. As the core now faces the part of the overdensity that is located close to the orbital plane, mass accretion becomes symmetric about the orbital plane. Finally, the core eventually approaches the northern hemisphere end of the overdensity, analogously the accretion streams pulled from this part of the lump reach the inner boundary in the southern hemisphere. As one core orbits away from the lump and becomes subject to symmetric mass accretion (Figs. B(b) and C(a)), the companion approaches it, and the phenomenon is repeated.

All Tables

All Figures

	Fig. 1. Initial density, pressure, and sound speed profiles used in all our simulation runs.
In the text

	Fig. 2. Zoomed-in snapshots of density cross section in the xy and xz planes at different times and for our three inviscid simulations runs A (left), B (middle) and C (right). The snapshots on the first line are taken shortly after the end of the initial spin-up phase.
In the text

	Fig. 3. Kinetic, internal, gravitational, and total binding energy evolution in units GM2/R for model A. The vertical black dashed line indicates the replacement of the averaged binary potential with its full expression (Eq. (4)), and the bump in energy is indicated with the double arrow.
In the text

	Fig. 4. Angular momentum evolution and conservation in our simulations. Panel a: time evolution of the advective, viscous, and gravitational torques for runs A and, A′ (dotted lines). Panel b: relative difference between the sum of the torques and the measured time derivative of the total angular momentum showing the angular momentum conservation for run A.
In the text

	Fig. 5. Time average of the angular distribution of mass (panel a) and angular momentum (panel b) accretion fluxes through the inner boundary for ten orbital periods of model A at 148 ≤ t/Porb ≤ 158.
In the text

Fig. 6. Evolution of key quantities relevant for the binary orbit after the initial envelope spin-up and adjustment. Panel a: moving average of the specific angular momentum transfer ⟨j⟩P (Eq. (32)) for models A, B, C, and D. The black dashed line indicates the critical specific angular momentum transfer jcrit = 3/8. Panel b: moving average of the gravitational torque for simulation runs A and A′. The black dashed line indicates . Panel c: moving average of the mass accretion rate through the inner boundary ⟨Ṁ⟩P. Panel d: envelope mass in our numerical domain Menv(r ≤ 10). Panel e: moving average of the orbital separation evolution timescale .

In the text

	Fig. 7. Space-time diagram of the mass (a) and angular momentum (b) fluxes onto the binary during 21 orbital periods for model A.
In the text

	Fig. 8. Same as Fig. 7 but for model B during 26 orbital periods.
In the text

Fig. 9. Detailed view on the variability of mass flux for model A. Panel a: space-time diagram of the local mass flux through the inner boundary. Panel b: time evolution of the mass accretion rate onto the binary. Panel c: power spectral density of the total mass accretion rate onto the binary (black line) and of the mass flux at each colatitude (colored lines). Green lines correspond to the range 0 ≤ θ ≤ π/3, blue lines to the range 2π/3 ≤ θ ≤ π, and orange lines to the range π/3 ≤ θ ≤ 2π/3.

In the text

	Fig. 10. Detailed view on the variability of mass flux for model B and its two different accretion regimes (left and right panels). Meaning of symbols and lines in each panel is the same as in Fig. 9.
In the text

Fig. 11. Space-time diagram of the m = 1 mode of the Fourier transform of the θ-integrated density with respect to the azimuth φ for model A. Panel b is obtained with a much larger time resolution than panel a, such that high frequency fluctuations are well resolved. On panel b, we see the overdensity generated at early time (t ≃ 147 Porb) propagating outward and expanding radially, and a new lump building up from t ≃ 152 Porb.

In the text

	Fig. 12. Surface density about the orbital plane (35) for model A (top) and A′ (bottom) at t = 158 Porb. We see the outward propagation of a lump and the formation of a new one for model A, and the absence of a structured lump in model A′ resulting from the absence of accretion streams. Black dots indicate the position of the two cores.
In the text

	Fig. 13. Space-time diagram of shell-averaged envelope eccentricity. The meaning of symbols is the same as in Fig. 11.
In the text

	Fig. 14. Evolution of the mean envelope eccentricity within the numerical domain. The black dashed line shows a linear fit yielding a growth rate λenv ≃ 0.022 Ωorb.
In the text

	Fig. 15. First Solberg–Høiland criterion for convective stability for run B at t = 20 Porb (top row), and t = 250 Porb (bottom row). Negative values indicate convective instability according to Eq. (38).
In the text

	Fig. 16. Advective and gravitational contributions to the local angular momentum transfer rate for models A (full lines), model A′ (dash-dotted lines), model B (dashed lines), and model C (dotted lines). The quantities are averaged in time interval 250 ≤ t/Porb ≤ 275.
In the text

	Fig. 17. Azimuthally averaged advective radial and latitudinal fluxes of angular momentum, and averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).
In the text

	Fig. 18. Mean flow contribution to the azimuthally averaged advective angular momentum radial and latitudinal fluxes averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).
In the text

	Fig. 19. Turbulent flow contribution to the azimuthally averaged advective angular momentum radial and latitudinal fluxes averaged over time interval 250 ≤ t/Porb ≤ 275 for models A (panel a), A′ (panel b), B (panel c), and C (panel d).
In the text

	Fig. 20. Mean and turbulent components of the advective contribution to local angular momentum transfer rate for models A (full lines), A′ (dash-dotted lines), B (dashed lines), and C (dotted lines). The quantities are averaged over a time interval 250 ≤ t/Porb ≤ 275. Neglecting density perturbations, .
In the text

Fig. 21. Time evolution of the integral radial scales Λ (Eq. (52)) in three regions I (r ≤ Rdomain/10), II (Rdomain/10 < r ≤ Rdomain/2), and III (r > Rdomain/2) for simulation runs A (full lines), B (dashed lines), and C (dotted lines). The discontinuity of ΛIII for run A is a result of not changing sign on the same radial scales for all φ, that is Rθθ(r′,t) does not cross zero during a few orbital periods.

In the text

	Fig. 22. Mean, turbulent, gravitational, and viscous contributions to local angular momentum transfer rate across the common envelope for model D, averaged in time from t = 140 Porb to t = 150 Porb.
In the text

	Fig. 23. Mean flow (a), turbulent flow (b), and viscous (c) contributions to the azimuthally averaged advective radial and latitudinal angular momentum fluxes for the viscous model D. We averaged the quantities over ten orbital periods for 140 ≤ t/Porb ≤ 150.
In the text

	Fig. 24. Ratio of thermal timescale tKH to the orbital period of the binary Porb inside the shared envelope in the Hertzprung–Russel diagram constructed for four solar-metallicity nonrotating evolutionary tracks from the MIST database (Dotter 2016; Choi et al. 2016). Here, tKH = G(M1 + Menv)2/(2RL), where L is the luminosity of the star.
In the text

	Fig. B.1. Specific gravitational torque for simulation run A at t ≃ 25 Porb (left) and t ≃ 250 Porb (middle), and for run A’ at t ≃ 250 Porb (right). Black crosses indicate the location of the two cores, and the dashed black line indicates the intersection of 𝒫 with the orbital plane.
In the text

	Fig. B.2. (a) quantifying the asymmetry of the mass distribution with respect to the plane 𝒫 as a function of time. (b) Weighted-averaged density with a dominating contribution from the gas subject to larger |∂Φ/∂φ|, , as a function of time.
In the text

Fig. C.1. A: Space-time diagram of the mass flux onto the individual components of the binary for model A, azimuthally averaged in the range φ1 − π/12 ≤ φ1 ≤ φ1 + π/12 (a) and in the range φ2 − π/12 ≤ φ2 ≤ φ2 + π/12 (b). The horizontal dashed lines indicate the orbital plane θ = π/2. The two vertical dashed lines are separated by Porb/2 and indicate two successive times, t1 and t2, when the mass flux is maximum at the highest latitude. B: Mass flux toward and away from the individual cores, azimuthally averaged in the range φ1 − π/12 ≤ φ1 ≤ φ1 + π/12 (a) and in the range φ2 − π/12 ≤ φ2 ≤ φ2 + π/12 (b) as a function of colatitude and radius at t = t1. C: same as B but at t = t2. D: Surface density about the orbital plane (35) averaged over t1 ≤ t ≤ t2. The dashed circle indicates the approximate radial location of the overdensity responsible for the migrating accretion streams (r = 0.7). E: Density cross section in the θφ plane at r = 0.7 averaged over t1 ≤ t ≤ t2.

In the text