Issue 
A&A
Volume 674, June 2023



Article Number  A121  
Number of page(s)  26  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202346057  
Published online  13 June 2023 
Postdynamical inspiral phase of common envelope evolution
Binary orbit evolution and angular momentum transport
Institute of Theoretical Physics, Faculty of Mathematics and Physics, Charles University, V Holešovičkách 2, Praha 8, 180 00
Czech Republic
email: damien.gagnier@matfyz.cuni.cz
Received:
1
February
2023
Accepted:
13
April
2023
After the companion dynamically plunges through the primary’s envelope, the two cores remain surrounded by a common envelope and the decrease of the orbital period P_{orb} stalls. The subsequent evolution has never been systematically explored with multidimensional simulations. For this study, we performed 3D hydrodynamical simulations of an envelope evolving under the influence of a central binary star using an adaptively refined spherical grid. We followed the evolution over hundreds of orbits of the central binary to characterize the transport of angular momentum by advection, gravitational torques, turbulence, and viscosity. We find that local advective torques from the mean flow and Reynolds stresses associated with the turbulent flow dominate the angular momentum transport, which occurs outward in a disklike structure about the orbital plane and inward along the polar axis. Turbulent transport is less efficient, but can locally significantly damp or enhance the net angular momentum radial transport and may even reverse its direction. Shortterm variability in the envelope is remarkably similar to circumbinary disks, including the formation and destruction of lumplike overdensities, which enhance mass accretion and contribute to the outward transport of eccentricity generated in the vicinity of the binary. If the accretion onto the binary is allowed, the orbital decay timescale settles to a nearly constant value τ_{b} ∼ 10^{3} to 10^{4} P_{orb}, while preventing accretion leads to a slowly increasing τ_{b} ∼ 10^{5} P_{orb} at the end of our simulations. Our results suggest that the postdynamical orbital contraction and envelope ejection will slowly continue while the binary is surrounded by gas and that τ_{b} is often much shorter than the thermal timescale of the envelope.
Key words: binaries: close / methods: numerical / stars: kinematics and dynamics
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 quasisteady spiralin 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 selfregulating 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 postCEE 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), Xray 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 spacebased gravitationalwave detectors such as LISA (Thorpe et al. 2019) or TianQin (Huang et al. 2020) during the postdynamical inspiral 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. Threedimensional 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 plungein when only orbital energy injection by the secondary star is considered, and the obtained postdynamical inspiral orbital separations are often larger than that observed in postCE systems (e.g., Nebot GómezMorá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 & MeyerHofmeister 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 postdynamical selfregulated 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 twostep prescription combining energy and angular momentum (Hirai & Mandel 2022).
The postdynamical selfregulating 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 thermaltimescale selfregulation 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 inspiral, becomes very weak (e.g., Ostriker 1999; Ricker & Taam 2012; MacLeod & RamirezRuiz 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 postdynamical 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 postdynamical 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 postCEE binary. Our connection between postdynamical CEE and CBDs is different from previous explorations of fallback CBDs around postCEE binaries (De Marco et al. 2011; Kashi & Soker 2011).
There are also important differences between CBDs and the postdynamical 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ópezCá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 postdynamical inspiral of CEE by performing the first dedicated series of 3D hydrodynamical numerical simulations. To establish a wellcontrolled 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 timechanging 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 shortterm 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 postdynamical 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 spunup 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 plungein (see Sect. 2.4). We set the gravitational constant G, the total binary mass M = M_{1} + M_{2}, the primary’s initial radius R, and thus the angular velocity to unity. The orbital velocity is fixed to , where a_{b} = r_{1} + r_{2} is the fixed binary separation, M_{1} is the mass of the primary’s core located at {r_{1}, θ_{1}, φ_{1}}, and M_{2} is the mass of the secondary object (either a mainsequence star or a compact object) located at {r_{2}, θ_{2}, φ_{2}}. 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_{2}/M_{1} = 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.
Setup comparison to several ab initio simulations of CEE.
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
where E = e + ρu^{2}/2, e is the internal energy density, P = (Γ − 1)e, Γ = 5/3 is the adiabatic index, Φ(r) is the gravitational potential of the binary,
and T_{ij} is the symmetric viscous stress tensor,
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 = α_{1}H_{P} (e.g., Vitense 1953; Zahn 1989), and that the characteristic eddy velocity is a fraction of the local sound speed, v = α_{2}c_{s}. Following Shakura & Sunyaev (1973), we obtain ν(r, θ, φ) = α_{ν}c_{s}H_{P}, where α_{ν} = α_{1}α_{2}/3 is a free parameter. Assuming a typical effective kinematic viscosity of 𝒪(10^{15} cm^{2} s^{−1}), we take α_{ν} = 10^{−3} 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 RungeKuttaLegendre supertimestepping 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 spinup 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 firstorder 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
where is the cellface 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,
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 spinup 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 diodetype 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 diodetype 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 spinup 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 lowdensity 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 selfgravity and considering purely radial initial profiles, the equations governing the envelope initial structure read
where K is a constant related to the thermal conditions at the inner boundary. The Green’s function for Eq. (4) satisfies (Jackson 1975)
for r ≥ r_{in} ≥ max(r_{1}, r_{2}), 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
where ⟨ ⋅ ⟩_{θ} indicates a time and latitude average. The ℓ ≤ 2 latitude and time averaged binary potential finally reads
We use this latitude and timeaveraged potential in the momentum equation during spinup. We replace the averaged potential by its timedependent expression (Eq. (4)) after the spinup. Our choice of the initial potential facilitates the transition from the initial spinup 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, κ^{n} = ρ(R)/ρ(r_{in}), we obtain the initial density and pressure profiles
where n = 1/(Γ − 1) is the polytropic index and
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 lowdensity medium to which we apply our outer boundary conditions and in which the envelope will expand later on. To model this lowdensity 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
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_{1} such that the stellar surface is in hydrostatic equilibrium with the lowdensity atmosphere, which gives
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 GaussLegendre 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 polarspherical coordinates, especially near the polar axis because of the converging grid geometry.
2.4. Initial spinup
We aim to construct a model with an initial total angular momentum that is consistent with what is available in the system,
where a_{i} is the initial binary separation^{1}, a_{b} is the enforced separation at the end of the dynamical plungein, 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 plungein. is the orbital angular momentum of the binary at the end of the dynamical plungein, which coincides with the beginning of our simulations, and μ = M_{1}M_{2}/(M_{1} + M_{2}) 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 spinup rate to all cells in which the angular velocity is subKeplerian, . Simultaneously, the structure of the envelope slowly restructures.
We stop the spinup 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 spiralin 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 timestep constraints resulting from the converging grid geometry and the CourantFriedrichsLewy (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 postprocessing 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
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 secondorder spatial reconstruction procedure to the averaged values and we correct the minimum timestep 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 spinup 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
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 wellsuited 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,
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. Volumeaveraged 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, facecentered 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 valueunsafe optimizations of floatingpoint operations or the nonassociativity of floatingpoint 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^{−6} ρ(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 zoomedin 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 plungein, the density structure and flow morphology early in our simulations have striking resemblance with latetime snapshots from ab initio simulations of dynamical plungein (e.g., Ohlmann et al. 2016; Chamandy et al. 2020).
Fig. 2. Zoomedin 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 spinup 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 largescale 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 spinup 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 GM^{2}/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. 
Run parameters and simulations outcome.
To asses the importance of the bump, we estimate the CEE efficiency parameter α_{CEE} corresponding to injection of internal energy ΔE,
where a_{i} is the initial binary orbital separation before plungein 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 plungein 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 postdynamical 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 selfconsistently 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
where is the advective torque associated with the loss of angular momentum through the boundaries,
is the gravitational torque exerted by the binary,
and is the viscous torque,
Here, n_{⊥} is the outwardpointing 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 spinup, 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 timesteady so that we can qualitatively assume that all initial transients have decayed and that the flow properties have reached a quasisteady 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
If the central binary does not accrete from the shared envelope, Eq. (29) simplifies to
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
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
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/P_{orb} ≤ 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
Fig. 6. Evolution of key quantities relevant for the binary orbit after the initial envelope spinup 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 j_{crit} = 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 M_{env}(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^{3} 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^{5} 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^{−3} 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 spacetime 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. Spacetime diagram of the mass (a) and angular momentum (b) fluxes onto the binary during 21 orbital periods for model A. 
Fig. 9. Detailed view on the variability of mass flux for model A. Panel a: spacetime 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 lowdensity cavity encompassing the central binary. For t ≳ 156 P_{orb}, the mass and angular momentum accretion becomes more isotropic, suggesting the absence of such wellstructured 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 shortterm 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 spacetime evolution of quantity A_{1}(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 φ,
Fig. 11. Spacetime 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 P_{orb}) propagating outward and expanding radially, and a new lump building up from t ≃ 152 P_{orb}. 
First, we see a highfrequency variation of A_{1} 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,
Fig. 12. Surface density about the orbital plane (35) for model A (top) and A′ (bottom) at t = 158 P_{orb}. 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 equalmass 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 spacetime diagram of shellaveraged envelope eccentricity
Fig. 13. Spacetime diagram of shellaveraged envelope eccentricity. The meaning of symbols is the same as in Fig. 11. 
Similarly to A_{1}, 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
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 latitudeaveraged 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 quasiequilibrium that may be maintained throughout the entire postdynamical spiralin 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 postdynamical 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
where
ℓ_{z} = s^{2}Ω 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 latetime. 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 plungein from Ohlmann et al. (2016), where the theoretically stable and unstable layers alternate in a geometrically thick disklike 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 P_{orb} (top row), and t = 250 P_{orb} (bottom row). Negative values indicate convective instability according to Eq. (38). 
3.6.2. Local torque balance
While the various volumeintegrated 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
where
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 meanflow angular structure.
Fig. 16. Advective and gravitational contributions to the local angular momentum transfer rate for models A (full lines), model A′ (dashdotted lines), model B (dashed lines), and model C (dotted lines). The quantities are averaged in time interval 250 ≤ t/P_{orb} ≤ 275. 
3.6.3. Turbulent transport of angular momentum in the envelope
We can use hydrodynamic meanfield 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 zdirection, using Reynolds decomposition to define velocity fluctuations about their averages as , and ignoring density fluctuations from its mean value, we obtain the conservation law
where is the meridional mean velocity, overlined quantities indicate azimuthal average, and and S′ are the rateofstrain tensors of the mean flow and of the fluctuating flow defined as
The first term on the righthand 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,
This stress tensor is often separated into nondiffusive (“Λeffect”) and diffusive contributions described by turbulent viscosity,
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 disklike structure with an opening angle that is smaller for higher initial angular momentum. Outside of the disklike 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 disklike 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/P_{orb} ≤ 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/P_{orb} ≤ 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/P_{orb} ≤ 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 meanflow and turbulent contributions to the angular momentum transfer in terms of Reynolds stress,
Fig. 20. Mean and turbulent components of the advective contribution to local angular momentum transfer rate for models A (full lines), A′ (dashdotted lines), B (dashed lines), and C (dotted lines). The quantities are averaged over a time interval 250 ≤ t/P_{orb} ≤ 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 largescale meanflow, 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 autocorrelation 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}],
We also introduce the integral radial scale (e.g., Townsend 1976; O’neill et al. 2004; Mora & Obligado 2020),
Because the low amplitude tail of Eq. (51) may contain information extraneous to turbulent motion, we integrate Eq. (52) up to the first zerocrossing 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 energycontaining 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 autocorrelation 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 timeaveraged values, Λ^{tmin − tmax}, 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 postdynamical 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 postdynamical 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 spinup setup where the innermost layers are spunup to critical rotation, while they remain subcritical during spinup 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 ≤ R_{domain}/10), II (R_{domain}/10 < r ≤ R_{domain}/2), and III (r > R_{domain}/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 P_{orb} to t = 150 P_{orb}. 
4. Discussions
4.1. Comparison with CBDs
In this work, we draw many analogies between simulations of CBDs and postdynamical 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 supermassive 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 postdynamical shared envelope in CEE contains large amount of mass, which prevents cooling and keeps the geometry strictly threedimensional. Ultimately, the shared envelope disperses and whatever gas remains should cool to a thin disks, as is observed in postAGB binaries (e.g., Dermine et al. 2013; Kluska et al. 2022). The transition between these two regimes of postdynamical 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 shortterm 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 postdynamical 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^{3} to 10^{4} orbital periods of the binary when accretion is allowed or slowly decreases to τ_{b} ∼ 10^{5} 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 plungein 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 “selfregulating” 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^{3} to 10^{5} 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_{1}/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 lowmass red giants without fully developed convective envelopes, R ≲ 100 R_{⊙}, the ratio is t_{KH}/P_{orb} ≳ 10^{5}, 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 longlasting processes are strong pulsations or dustdriven wind (Clayton et al. 2017; Glanz & Perets 2018). For AGB and highmass stars, the thermal coupling between the inspiralled binary and the envelope seems more likely. However, unless t_{KH}/P_{orb} ≲ 10^{3}, 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 plungein 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 postdynamical phase, as Fig. 6 suggests. Furthermore, CEE is often preceded by strong thermaltimescale 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 longterm 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 energyconserving α_{CEE} formalism that is commonly used to predict CEE outcomes. Our results suggest that much of the orbital decay during the postdynamical phase might be over before thermal effects become important, which lends support to the energyconserving formalism. At the same time, the orbital decay is clearly separated into two regimes: a fast dynamical plungein and much slower postdynamical 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 postdynamical 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 twostep 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 plungein. This is in agreement with 3D hydrodynamical simulations, which typically find quasicircular 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 postdynamical 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 postdynamical spiralin 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 postAGB and postCEE 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 postdynamical inspiral phase of CEE. We used the procedure of Morris & Podsiadlowski (2006, 2007, 2009) to mimic the outcome of the preceding dynamical plungein 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 quasisteady value of 𝒪(10^{3} − 10^{4} 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^{5} 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 plungein. This suggests a significant reduction of orbital separation and more efficient envelope ejection is possible even after the dynamical plungein 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 postdynamical inspiral timescale is much shorter than the thermal timescale for primary stars with radius ≲100 R_{⊙}. Even for larger stars, the postdynamical decay does not have to be significantly influenced by the thermal response of the envelope and thermal “selfregulation” 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) twostep process.
Our second aim was to find the typical frequencies associated with the shortterm 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 postdynamical 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 shortterm variability. Finally, we found that envelope eccentricity is excited in the vicinity of the binary and propagates outward within and inbetween 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 plungein simulations, gravitational perturbations from the orbiting binary during postdynamical inspiral 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 disklike 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^{−3} 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/P_{orb} ≤ 150. 
Fig. 24. Ratio of thermal timescale t_{KH} to the orbital period of the binary P_{orb} inside the shared envelope in the Hertzprung–Russel diagram constructed for four solarmetallicity nonrotating evolutionary tracks from the MIST database (Dotter 2016; Choi et al. 2016). Here, t_{KH} = G(M_{1} + M_{env})^{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 binaryenvelope 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 binaryenvelope 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.
We ignore potential mass and angular momentum loss from the outer Lagrange point (L_{2}) preceding common envelope (e.g., Shu et al. 1979; Pejcha 2014; Pejcha et al. 2016b; Hubová & Pejcha 2019). Such additional angular momentum loss can be mimicked by lowering the value of β.
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 ‘CatInhAT’ (grant agreement no. 803158). This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the eINFRA 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 PHY1748958.
References
 Ablimit, I., Maeda, K., & Li, X.D. 2016, ApJ, 826, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [Google Scholar]
 Artymowicz, P., Clarke, C. J., Lubow, S. H., & Pringle, J. E. 1991, ApJ, 370, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
 Bate, M. R., & Bonnell, I. A. 1997, MNRAS, 285, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
 Belczynski, K., Bulik, T., & Ruiter, A. J. 2005, ApJ, 629, 915 [Google Scholar]
 Blagorodnova, N., Klencki, J., Pejcha, O., et al. 2021, A&A, 653, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cabot, W. 1996, ApJ, 465, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Chamandy, L., Frank, A., Blackman, E. G., et al. 2018, MNRAS, 480, 1898 [NASA ADS] [CrossRef] [Google Scholar]
 Chamandy, L., Tu, Y., Blackman, E. G., et al. 2019a, MNRAS, 486, 1070 [NASA ADS] [CrossRef] [Google Scholar]
 Chamandy, L., Blackman, E. G., Frank, A., et al. 2019b, MNRAS, 490, 3727 [NASA ADS] [CrossRef] [Google Scholar]
 Chamandy, L., Blackman, E. G., Frank, A., CarrollNellenback, J., & Tu, Y. 2020, MNRAS, 495, 4028 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, W.C., Liu, D.D., & Wang, B. 2020, ApJ, 900, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
 Clayton, M., Podsiadlowski, P., Ivanova, N., & Justham, S. 2017, MNRAS, 470, 1788 [NASA ADS] [CrossRef] [Google Scholar]
 De, S., MacLeod, M., Everson, R. W., et al. 2020, ApJ, 897, 130 [Google Scholar]
 De Marco, O. 2009, PASP, 121, 316 [NASA ADS] [CrossRef] [Google Scholar]
 De Marco, O., Passy, J.C., Moe, M., et al. 2011, MNRAS, 411, 2277 [CrossRef] [Google Scholar]
 Dermine, T., Izzard, R. G., Jorissen, A., & Van Winckel, H. 2013, A&A, 551, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Di Stefano, R., Kruckow, M. U., Gao, Y., Neunteufel, P. G., & Kobayashi, C. 2023, ApJ, 944, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Dittmann, A. J., & Ryan, G. 2021, ApJ, 921, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [Google Scholar]
 D’Orazio, D. J., & Duffell, P. C. 2021, ApJ, 914, L21 [CrossRef] [Google Scholar]
 D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, MNRAS, 436, 2997 [Google Scholar]
 Dotter, A. 2016, ApJS, 222, 8 [Google Scholar]
 Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020, ApJ, 901, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Fragos, T., Andrews, J. J., RamirezRuiz, E., et al. 2019, ApJ, 883, L45 [Google Scholar]
 Gagnier, D., & Rieutord, M. 2020, J. Fluid Mech., 904, A35 [NASA ADS] [CrossRef] [Google Scholar]
 Garaud, P., Gagnier, D., & Verhoeven, J. 2017, ApJ, 837, 133 [Google Scholar]
 Glanz, H., & Perets, H. B. 2018, MNRAS, 478, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Glanz, H., & Perets, H. B. 2021, MNRAS, 507, 2659 [NASA ADS] [CrossRef] [Google Scholar]
 Grichener, A., Sabach, E., & Soker, N. 2018, MNRAS, 478, 1818 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F. 2000, ApJ, 528, 462 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., & Krolik, J. H. 2001, ApJ, 548, 348 [CrossRef] [Google Scholar]
 Heath, R. M., & Nixon, C. J. 2020, A&A, 641, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirai, R., & Mandel, I. 2022, ApJ, 937, L42 [NASA ADS] [CrossRef] [Google Scholar]
 Hirai, R., Podsiadlowski, P., Owocki, S. P., Schneider, F. R. N., & Smith, N. 2021, MNRAS, 503, 4276 [NASA ADS] [CrossRef] [Google Scholar]
 Høiland, E. 1941, Avhandliger Norske VidenskapsAkademi i Oslo, i. math naturv. Klasse, 1 [Google Scholar]
 Huang, S.J., Hu, Y.M., Korol, V., et al. 2020, Phys. Rev. D, 102, 063021 [NASA ADS] [CrossRef] [Google Scholar]
 Hubová, D., & Pejcha, O. 2019, MNRAS, 489, 891 [Google Scholar]
 Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
 Iaconi, R., & De Marco, O. 2019, MNRAS, 490, 2550 [Google Scholar]
 Iaconi, R., Maeda, K., De Marco, O., Nozawa, T., & Reichardt, T. 2019, MNRAS, 489, 3334 [NASA ADS] [CrossRef] [Google Scholar]
 Iaconi, R., Maeda, K., Nozawa, T., De Marco, O., & Reichardt, T. 2020, MNRAS, 497, 3166 [NASA ADS] [CrossRef] [Google Scholar]
 Iben, I., Jr., & Tutukov, A. V. 1984, ApJS, 54, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Ivanova, N., & Nandez, J. L. A. 2016, MNRAS, 462, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Ivanova, N., Justham, S., Chen, X., et al. 2013a, A&ARv, 21, 59 [Google Scholar]
 Ivanova, N., Justham, S., Avendano Nandez, J. L., & Lombardi, J. C. 2013b, Science, 339, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, J. D. 1975, Classical Electrodynamics (New York: Wiley) [Google Scholar]
 Jones, D., & Boffin, H. M. J. 2017, Nat. Astron., 1, 0117 [Google Scholar]
 Jones, C. A., Kuzanyan, K. M., & Mitchell, R. H. 2009, J. Fluid Mech., 634, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Kalogera, V., & Webbink, R. F. 1998, ApJ, 493, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2019, A&A, 622, A195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Guerrero, G., Brandenburg, A., & Chatterjee, P. 2011, A&A, 531, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kashi, A., & Soker, N. 2011, MNRAS, 417, 1466 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Ruediger, G. 1995, A&A, 299, 446 [NASA ADS] [Google Scholar]
 Kitchatinov, L. L., Pipin, V. V., & Ruediger, G. 1994, Astron. Nachr., 315, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54 [EDP Sciences] [Google Scholar]
 Kluska, J., Van Winckel, H., Coppée, Q., et al. 2022, A&A, 658, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kochanek, C. S., Adams, S. M., & Belczynski, K. 2014, MNRAS, 443, 1319 [NASA ADS] [CrossRef] [Google Scholar]
 Kruckow, M. U., Neunteufel, P. G., Di Stefano, R., Gao, Y., & Kobayashi, C. 2021, ApJ, 920, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, M. Y. M., Hirai, R., GonzálezBolívar, M., et al. 2022a, MNRAS, 512, 5462 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, M. Y. M., Hirai, R., Price, D. J., & Mandel, I. 2022b, MNRAS, 516, 4669 [NASA ADS] [CrossRef] [Google Scholar]
 Livio, M., & Soker, N. 1988, ApJ, 329, 764 [Google Scholar]
 Lohner, R., Morgan, K., Peraire, J., & Vahdati, M. 1987, Int. J. Numer. Meth. Fluids, 7, 1093 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez Armengol, F. G., Combi, L., Campanelli, M., et al. 2021, ApJ, 913, 16 [NASA ADS] [CrossRef] [Google Scholar]
 LópezCámara, D., De Colle, F., & Moreno Méndez, E. 2019, MNRAS, 482, 3646 [Google Scholar]
 LópezCámara, D., De Colle, F., Moreno Méndez, E., Shiber, S., & Iaconi, R. 2022, MNRAS, 513, 3634 [CrossRef] [Google Scholar]
 Lubow, S. H. 1991a, ApJ, 381, 259 [Google Scholar]
 Lubow, S. H. 1991b, ApJ, 381, 268 [NASA ADS] [CrossRef] [Google Scholar]
 MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83 [Google Scholar]
 MacLeod, M., & RamirezRuiz, E. 2015, ApJ, 803, 41 [Google Scholar]
 MacLeod, M., Macias, P., RamirezRuiz, E., et al. 2017a, ApJ, 835, 282 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, M., Antoni, A., MurguiaBerthier, A., Macias, P., & RamirezRuiz, E. 2017b, ApJ, 838, 56 [Google Scholar]
 MacLeod, M., Ostriker, E. C., & Stone, J. M. 2018, ApJ, 863, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating stars (Berlin, Heidelberg: Springer) [Google Scholar]
 Marchant, P., Pappas, K. M. W., GallegosGarcia, M., et al. 2021, A&A, 650, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumoto, T., & Metzger, B. D. 2022, ApJ, 938, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, T., Saigo, K., & Takakuwa, S. 2019, ApJ, 871, 36 [CrossRef] [Google Scholar]
 Metzger, B. D., & Pejcha, O. 2017, MNRAS, 471, 3200 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, F., & MeyerHofmeister, E. 1979, A&A, 78, 167 [NASA ADS] [Google Scholar]
 Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2014, J. Comput. Phys., 257, 594 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
 Milosavljević, M., & Phinney, E. S. 2005, ApJ, 622, L93 [CrossRef] [Google Scholar]
 Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170 [NASA ADS] [CrossRef] [Google Scholar]
 Moody, M. S. L., Shi, J.M., & Stone, J. M. 2019, ApJ, 875, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Mora, D., & Obligado, M. 2020, Exp. Fluids, 61, 199 [CrossRef] [Google Scholar]
 Moreno Méndez, E., LópezCámara, D., & De Colle, F. 2017, MNRAS, 470, 2929 [CrossRef] [Google Scholar]
 Morris, T., & Podsiadlowski, P. 2006, MNRAS, 365, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Morris, T., & Podsiadlowski, P. 2007, Science, 315, 1103 [Google Scholar]
 Morris, T., & Podsiadlowski, P. 2009, MNRAS, 399, 515 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106 [CrossRef] [Google Scholar]
 Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84 [CrossRef] [Google Scholar]
 Nandez, J. L. A., Ivanova, N., & Lombardi, J. C. J. 2015, MNRAS, 450, L39 [Google Scholar]
 Nebot GómezMorán, A., Gänsicke, B. T., Schreiber, M. R., et al. 2011, A&A, 536, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelemans, G., Verbunt, F., Yungelson, L. R., & Portegies Zwart, S. F. 2000, A&A, 360, 1011 [NASA ADS] [Google Scholar]
 Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Ohlmann, S. T., Röpke, F. K., Pakmor, R., & Springel, V. 2016, ApJ, 816, L9 [Google Scholar]
 O’neill, P. L., Nicolaides, D., Honnery, D., & Soria, J. 2004, Autocorrelation Functions and the Determination of Integral Length with Reference to Experimental and Numerical Data (Sydney: University of Sydney) [Google Scholar]
 Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]
 Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, eds. P. Eggleton, S. Mitton, & J. Whelan, 73, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Passy, J.C., De Marco, O., Fryer, C. L., et al. 2012, ApJ, 744, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Pejcha, O. 2014, ApJ, 788, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Pejcha, O., Metzger, B. D., & Tomida, K. 2016a, MNRAS, 455, 4351 [Google Scholar]
 Pejcha, O., Metzger, B. D., & Tomida, K. 2016b, MNRAS, 461, 2527 [NASA ADS] [CrossRef] [Google Scholar]
 Pejcha, O., Metzger, B. D., Tyles, J. G., & Tomida, K. 2017, ApJ, 850, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Penzlin, A. B. T., Kley, W., Audiffren, H., & Schäfer, C. M. 2022, A&A, 660, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Politano, M. 2021, A&A, 648, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rüdiger, G. 2022, Differential Rotation and Stellar Convection: Sun and Solartype Stars (De Gruyter) [Google Scholar]
 Reichardt, T. A., De Marco, O., Iaconi, R., Tout, C. A., & Price, D. J. 2019, MNRAS, 484, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Reichardt, T. A., De Marco, O., Iaconi, R., Chamandy, L., & Price, D. J. 2020, MNRAS, 494, 5333 [NASA ADS] [CrossRef] [Google Scholar]
 Renzo, M., Callister, T., Chatziioannou, K., et al. 2021, ApJ, 919, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, P. M., & Taam, R. E. 2012, ApJ, 746, 74 [Google Scholar]
 Roedig, C., Dotti, M., Sesana, A., Cuadra, J., & Colpi, M. 2011, MNRAS, 415, 3033 [NASA ADS] [CrossRef] [Google Scholar]
 Roedig, C., Sesana, A., Dotti, M., et al. 2012, A&A, 545, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roepke, F. K., & De Marco, O. 2023, Liv. Rev. Comput. Astrophys., 9, 2 [CrossRef] [Google Scholar]
 Sabach, E., Hillel, S., Schreier, R., & Soker, N. 2017, MNRAS, 472, 4361 [NASA ADS] [CrossRef] [Google Scholar]
 Sand, C., Ohlmann, S. T., Schneider, F. R. N., Pakmor, R., & Röpke, F. K. 2020, A&A, 644, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sandquist, E. L., Taam, R. E., Chen, X., Bodenheimer, P., & Burkert, A. 1998, ApJ, 500, 909 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shi, J.M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118 [Google Scholar]
 Shiber, S., Iaconi, R., De Marco, O., & Soker, N. 2019, MNRAS, 488, 5615 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H., Lubow, S. H., & Anderson, L. 1979, ApJ, 229, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N. 1993, ApJ, 417, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N., & Livio, M. 1994, ApJ, 421, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N., & Tylenda, R. 2006, MNRAS, 373, 733 [NASA ADS] [CrossRef] [Google Scholar]
 Solberg, M. 1936, Union Géodésique et Géophysique internationale VIeme assemblée, Edinburg, 66 [Google Scholar]
 Stone, J. M., & Balbus, S. A. 1996, ApJ, 464, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Taam, R. E., & Ricker, P. M. 2010, New Astron. Rev., 54, 65 [CrossRef] [Google Scholar]
 Taam, R. E., Bodenheimer, P., & Ostriker, J. P. 1978, ApJ, 222, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Tang, Y., MacFadyen, A., & Haiman, Z. 2017, MNRAS, 469, 4258 [Google Scholar]
 Teyssandier, J., & Ogilvie, G. I. 2016, MNRAS, 458, 3221 [NASA ADS] [CrossRef] [Google Scholar]
 Thorpe, J. I., Ziemer, J., Thorpe, I., et al. 2019, Bull. Am. Astron. Soc., 51, 77 [Google Scholar]
 Townsend, A. A. 1976, The Structure of Turbulent Shear Flow, 2nd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Vitense, E. 1953, Zeit. Astrophys., 32, 135 [NASA ADS] [Google Scholar]
 Webbink, R. F. 1984, ApJ, 277, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, E. C., & Nordhaus, J. 2019, MNRAS, 485, 4492 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1989, A&A, 220, 112 [NASA ADS] [Google Scholar]
 Zhang, B., Sorathia, K. A., Lyon, J. G., Merkin, V. G., & Wiltberger, M. 2019, J. Comput. Phys., 376, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Zrake, J., Tiede, C., MacFadyen, A., & Haiman, Z. 2021, ApJ, 909, L13 [NASA ADS] [CrossRef] [Google Scholar]
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,
Using the fact the viscous stress tensor T is symmetric, and writing r = x_{j}e_{j},
where ϵ_{ijk} is the LeviCivita tensor of rank three, and
Eq. (A.1) can be rewritten as
Finally, using Gauss divergence theorem, Eq. (A.4) can be rewritten as
where n_{⊥} is is the outwardpointing 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,
where
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
that is the plane orthogonal to the orbital plane and intersecting the binary semimajor axis. In our simulations, the gravitational torque arises from the quadrupolar moment of the binary potential generating a twoarmed 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
Fig. B.1. Specific gravitational torque for simulation run A at t ≃ 25 P_{orb} (left) and t ≃ 250 P_{orb} (middle), and for run A’ at t ≃ 250 P_{orb} (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,
Here, 0 ≤ ϵ ≤ 1 quantifies the asymmetry of the mass distribution with respect to 𝒫, and ρ_{eff} is the weightedaveraged 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}. Turbulenceinduced 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) Weightedaveraged 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 spacetime 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_{1}) and core 2 (at t = t_{2}), with a period of P_{orb}/2 = t_{2} − t_{1}, that is, when the two cores switch position. This suggests that the source of the highlatitude 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_{1} and t = t_{2}, respectively in Figs. C.1 (B) and (C). We see that highlatitude accretion through the inner boundary results from accretion streams originating from lowlatitude 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: Spacetime 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 P_{orb}/2 and indicate two successive times, t_{1} and t_{2}, 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 = t_{1}. C: same as B but at t = t_{2}. D: Surface density about the orbital plane (35) averaged over t_{1} ≤ t ≤ t_{2}. 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 t_{1} ≤ t ≤ t_{2}. 
We show the surface density about the orbital plane (35) averaged over t_{1} ≤ t ≤ t_{2} 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 highlatitude accretion. Finally, we show the density cross section in the θφ plane at r = 0.7 averaged over t_{1} ≤ t ≤ t_{2} 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. Zoomedin 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 spinup phase. 

In the text 
Fig. 3. Kinetic, internal, gravitational, and total binding energy evolution in units GM^{2}/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/P_{orb} ≤ 158. 

In the text 
Fig. 6. Evolution of key quantities relevant for the binary orbit after the initial envelope spinup 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 j_{crit} = 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 M_{env}(r ≤ 10). Panel e: moving average of the orbital separation evolution timescale . 

In the text 
Fig. 7. Spacetime 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: spacetime 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. Spacetime 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 P_{orb}) propagating outward and expanding radially, and a new lump building up from t ≃ 152 P_{orb}. 

In the text 
Fig. 12. Surface density about the orbital plane (35) for model A (top) and A′ (bottom) at t = 158 P_{orb}. 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. Spacetime diagram of shellaveraged 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 P_{orb} (top row), and t = 250 P_{orb} (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′ (dashdotted lines), model B (dashed lines), and model C (dotted lines). The quantities are averaged in time interval 250 ≤ t/P_{orb} ≤ 275. 

In the text 
Fig. 17. Azimuthally averaged advective radial and latitudinal fluxes of angular momentum, and averaged over time interval 250 ≤ t/P_{orb} ≤ 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/P_{orb} ≤ 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/P_{orb} ≤ 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′ (dashdotted lines), B (dashed lines), and C (dotted lines). The quantities are averaged over a time interval 250 ≤ t/P_{orb} ≤ 275. Neglecting density perturbations, . 

In the text 
Fig. 21. Time evolution of the integral radial scales Λ (Eq. (52)) in three regions I (r ≤ R_{domain}/10), II (R_{domain}/10 < r ≤ R_{domain}/2), and III (r > R_{domain}/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 P_{orb} to t = 150 P_{orb}. 

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/P_{orb} ≤ 150. 

In the text 
Fig. 24. Ratio of thermal timescale t_{KH} to the orbital period of the binary P_{orb} inside the shared envelope in the Hertzprung–Russel diagram constructed for four solarmetallicity nonrotating evolutionary tracks from the MIST database (Dotter 2016; Choi et al. 2016). Here, t_{KH} = G(M_{1} + M_{env})^{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 P_{orb} (left) and t ≃ 250 P_{orb} (middle), and for run A’ at t ≃ 250 P_{orb} (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) Weightedaveraged density with a dominating contribution from the gas subject to larger ∂Φ/∂φ, , as a function of time. 

In the text 
Fig. C.1. A: Spacetime 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 P_{orb}/2 and indicate two successive times, t_{1} and t_{2}, 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 = t_{1}. C: same as B but at t = t_{2}. D: Surface density about the orbital plane (35) averaged over t_{1} ≤ t ≤ t_{2}. 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 t_{1} ≤ t ≤ t_{2}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.