© The Authors 2025

1. Introduction

Common envelope evolution (CEE) is a phase in the life of binary systems during which one of the components expands and initiates a dynamically unstable mass transfer toward its more compact companion (e.g., Paczynski 1976; Meyer & Meyer-Hofmeister 1979; Ivanova et al. 2013; Röpke & De Marco 2023). The companion later becomes engulfed in the now shared envelope and the two stellar objects then rapidly spiral toward each other, depositing energy and angular momentum into the surrounding gas. The dynamical approach of the two cores can eventually slow due to the redistribution of the gas, and a phase of slow, quasi-stationary approach ensues. Ultimately, the common envelope will disperse and the two stellar cores will emerge as a post-common-envelope binary with properties often significantly different from the original binary. Because of the wide range of temporal and spatial scales that would need to be resolved, ab initio 3D (magneto)hydrodynamical simulations of this post-dynamical inspiral phase are often terminated well before the evolution naturally concludes.

In Gagnier & Pejcha (2023, 2024), we studied the post-dynamical inspiral phase of CEE by constructing idealized numerical experiments. In these simulations, a spherical envelope is first spun up to emulate angular momentum deposition by the dynamically inspiralling companion and then subjected to the evolving gravitational field of the two orbiting cores at the center. The (magneto)hydrodynamical simulations were performed in 3D using a spherical mesh with both adaptive and static mesh refinement. To reduce the computational cost of the simulations, we excised a spherical region around the two orbiting cores, removing them from the simulation domain. This also allowed us to control accretion onto the central regions by applying either an inflow or a reflecting inner boundary condition. With these simulations, we addressed issues such as angular momentum transfer between the binary and the envelope, the amplification of magnetic fields, and the effects of these fields on the envelope structure and dynamics, and on the binary's orbital evolution. We also showed that the structure and dynamics of the common envelope is remarkably similar to those of circumbinary disks. This similarity was also found in other recent ab initio 3D simulations (Ondratschek et al. 2022; Wei et al. 2024; Landri et al. 2025; Vetter et al. 2024) and used in 1D long-term models of CEE (Tuna & Metzger 2023; Valli et al. 2024).

Recent simulations of circumbinary disks in a more general context (e.g., Moody et al. 2019; Muñoz et al. 2019; Muñoz & Lithwick 2020; Duffell et al. 2020, 2024; Paschalidis et al. 2021; Li et al. 2021; Dittmann & Ryan 2021; Combi et al. 2022) have highlighted the crucial importance of adequate resolution of the cavity containing the binary in numerical simulations. This is essential in order to precisely quantify the torques that dictate the orbital evolution (e.g., Dittmann & Ryan 2021; Duffell et al. 2024), to analyze the potential mass transfer between the two binary components (e.g., Bowen et al. 2017; Avara et al. 2024), to study the evolution of their rotational velocity, or to understand the process of self-consistent launchings of jets (e.g., Paschalidis et al. 2021; Combi et al. 2022), for example.

In the context of the late phases of CEE, sufficient numerical resolution may become even more crucial than in circumbinary disks because of the potential absence of a low-density cavity encompassing the binary due to the more spherical geometry of the system. Consequently, our earlier simulations with an extracted central sphere presented in Gagnier & Pejcha (2023, 2024) provide limited information about processes occurring in the immediate vicinity of the central binary. In other CEE simulations, both the core of the primary star and the companion are generally represented as point masses, requiring the use of artificial softening of the gravitational potentials to avoid singularities. However, this approach inevitably results in reduced injection of orbital energy and angular momentum, as well as in nonphysical flow structures near the cores. This can potentially lead to erroneous envelope ejection rates and gravitational torques and diminished magnetic energy amplification. This problem is further exacerbated during the late phases of the CEE, when the stellar cores are close and the two softened regions constitute a non-negligible fraction of the intraorbital volume. Furthermore, the intraorbital region is often significantly under-resolved during this late phase.

To provide several examples of gravitational softening in CEE, in the magnetohydrodynamic (MHD) simulations of Ondratschek et al. (2022), the softening radius is fixed at ∼40% of the minimum orbital separation. In the hydrodynamical simulations of Kramer et al. (2020), Moreno et al. (2022), and Morán-Fraile et al. (2023), the softening radius is adaptively reduced to prevent the softened regions from overlapping. In Chamandy et al. (2020), the softening radius reaches up to ∼34% of the orbital separation, and up to 62% of the final separation in Prust & Chang (2019). Ohlmann et al. (2016a, b), Chamandy et al. (2019), Sand et al. (2020), Chamandy et al. (2024), and Vetter et al. (2024) maintain the softening radius below 20% of the binary separation. All of these works used Hernquist & Katz (1989)'s spline function to soften the gravitational potentials. Sandquist et al. (1998, 2000), Passy et al. (2012), Staff et al. (2016), Iaconi et al. (2018), and Shiber et al. (2019) instead used Ruffert (1993)'s softened potentials with softening radii often on the order of, or larger than, the final orbital separation. Perhaps even more alarming is the fact that, in some of these simulations, the final orbital separation only amounts to the width of a couple of mesh cells. In the context of circumbinary and protoplanetary disks, a Plummer-softened potential is almost exclusively used. However, Dong et al. (2011) showed that torques can be highly sensitive to potential shapes, and that low-order gravitational potentials such as Plummer's are particularly inaccurate. Such a wide spectrum of approaches and parameters raises the fundamental question of the potential impact of gravitational softening and intraorbital resolution on the ejection of the common envelope, on the amplification of magnetic energy within it, on the morphology of emerging planetary nebulae and the final separation of orbits, and on the reliability of current predictions.

We investigated the gas dynamics in the non-accreting binary's immediate vicinity and its effects on the orbital evolution with a special emphasis on the role of gravitational softening and intraorbital numerical resolution. To this end, we constructed a post-dynamical inspiral model similar to that of Gagnier & Pejcha (2023, 2024), emulating the preceding dynamical stage. In this work, however, we included both the binary and the intraorbital region in the numerical domain, using a statically refined Cartesian mesh to focus on the intraorbital region at the expense of the outer envelope.

In Sect. 2 we present our physical model and describe our numerical setup. In Sect. 3 we present our simulations and our analysis of gravitational softening formulations, gravitational torques, envelope ejection and morphology, and shear rates. In Sect. 4 we summarize and discuss our results in the context of the entire post-dynamical inspiral phase of CEE, as well as potential similarities with ordinary contact binaries.

2. Physical model and numerical setup

We constructed our post-dynamical inspiral model within the inertial frame centered on the binary's center of mass. M₁ is the mass of the primary's core located at {x₁, y₁, z₁} and M₂ is the mass of the secondary object, whether a main-sequence star or a compact object, located at {x₂, y₂, z₂}. Neither object is resolved and both are treated as point masses. We set the gravitational constant G, the total mass M=M₁+M₂ of the non-accreting binary, the orbital separation $a_{\mathrm{b}}^i$ at the onset of our simulations, and the associated orbital angular velocity ${\left(\mathit{GM}/a_{\mathrm{b}}^{i3}\right)}^{1/2}$ to unity. As a choice, we set the radius of the primary star R_*=50. To simplify our model and establish continuity with our prior works, we solely considered circular orbits ($e_{\mathrm{b}}={\dot{e}}_{\text{b}}=0$) and neglected mass accretion onto the individual cores. The mass of the envelope beyond the binary's orbit is set to M_env=2 in our units.

To convert the nondimensional units used in the simulations to dimensional cgs units, the reader should multiply the nondimensional values by the appropriate reference quantities. As an example, take a nondimensional density ρ, a reference mass of the binary M=M_⊙, and a reference length $a_{\mathrm{b}}^i=R_{*}/50$ with R_*=50 R_⊙. The dimensional density is then calculated as

$${\rho }_{\text{dim}}=5.9{\mathrm{g}\mathrm{cm}}^{-3}\rho \left(\frac{M}{M_{\odot }}\right){\left(\frac{a_{\text{b}}^i}{R_{\odot }}\right)}^{-3}.$$(1)

Instead of tracking the prior evolution of the inspiralling binary, we replicated its outcome through a method described by Morris & Podsiadlowski (2006, 2007, 2009), Hirai et al. (2021), and Gagnier & Pejcha (2023, 2024), where the envelope is spun-up until a satisfactory amount of total angular momentum,

$$J_z=(1+\beta )\frac{M_2(M_1+M_{\mathrm{env}})}{M+M_{\mathrm{env}}}\sqrt{G(M+M_{\mathrm{env}})m{a}_{\mathrm{b}}^i}-\mu \sqrt{\mathit{GM}{a}_{\mathrm{b}}^i},$$(2)

is imparted. Here, μ=M₁M₂/M is the binary's reduced mass and, in this work, we set β=0.1 and m=100. Our selection of initial parameters for the binary and envelope broadly mirrors the outcomes of ab initio simulations across various progenitor types. Table 1 summarizes the properties of our models. In the subsequent sections, we describe the equations used for solving the problem (Sect. 2.1), the initial conditions (Sect. 2.2), and the mesh structure (Sect. 2.3).

2.1. Equations of hydrodynamics

We used the Eulerian code Athena++ (Stone et al. 2020) to solve the compressible equations of hydrodynamics

$$\frac{\partial \rho }{\partial t}+\nabla ·\rho u=0,$$(3)

$$\frac{\partial \rho u}{\partial t}+\nabla ·(\rho u\otimes u+PI)=-\rho \nabla \mathrm{\Phi },$$(4)

$$\frac{\partial E}{\partial t}+\nabla ·\left((E+PI)u\right)=-\rho \nabla \mathrm{\Phi }·u,$$(5)

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

$$\mathrm{\Phi }(r)={\displaystyle \sum _{i=1}^2}{\mathrm{\Phi }}_i(r)=-{\displaystyle \sum _{i=1}^2}{\mathit{GM}}_{i}f(\Vert r-r_i\Vert ,ϵ),$$(6)

where ϵ is the gravitational softening length, $r_i=\sqrt{x_i^2+y_i^2+z_i^2}$ is the distance to the center of mass of either the primary's core or the secondary object, and $f(\Vert r-r_i\Vert ,ϵ)$ is the shape of the softened potential. In the Plummer model we have

$$f(\Vert r-r_i\Vert ,ϵ)=\frac{1}{\sqrt{{\Vert r-r_i\Vert }^2+{ϵ}^2}}·$$(7)

In the Ruffert (1993) model we have

$$f(\Vert r-r_i\Vert ,ϵ)=\frac{1}{\sqrt{{\Vert r-r_i\Vert }^2+{ϵ}^{2}exp\left(-{\Vert r-r_i\Vert }^2/{ϵ}^2\right)}}·$$(8)

In the Hernquist & Katz (1989) spline model we have

$$\begin{array}{cc}\hfill & f(\Vert r-r_i\Vert ,ϵ)\hfill \\ \hfill & =\frac{1}{h}\{\begin{array}{cc}\left(-\frac{16}{3}{u}^2+\frac{48}{5}{u}^4-\frac{32}{5}{u}^5+\frac{14}{5}\right),\hfill & \text{if}\hspace{0.5em}0\le u<\frac{1}{2},\hfill \\ \left(\frac{-1}{15u}-\frac{32}{3}{u}^2+16u^3-\frac{48}{5}{u}^4+\frac{32}{15}{u}^5+\frac{16}{5}\right),\hfill & \text{if}\hspace{0.5em}\frac{1}{2}\le u<1,\hfill \\ \frac{1}{u},\hfill & \text{if}\hspace{0.5em}u\ge 1.\hfill \end{array}\hfill \end{array}$$(9)

Here, $u=\Vert r-r_i\Vert /h$ and h≡14ϵ/5 so that the minimum of the potential has the same depth with the three formulations. We note that this definition of h is arbitrary, it has no physical implication, and it does not establish equivalence between the different formulations (see Appendix A). In our simulations, we gradually decreased the value of ϵ from 0.5 to its target value over a span of 20 orbital periods after the initial spin-up of the envelope.

2.2. Initial conditions

As in Gagnier & Pejcha (2023, 2024), we assumed that the gas in the envelope is initially in hydrostatic equilibrium and that it can be described by a polytropic equation of state ignoring the gas self-gravity and considering purely radial initial profiles. Because the binary's orbit is now included within the numerical domain, we also initialized the gas within the orbit as a polytrope in hydrostatic equilibrium. This initial stratification is marginally convectively stable. Since the equilibrium within and near the orbit is immediately perturbed by the binary motion, the specific details of the initial conditions in the inner regions are not important and we thus provide the corresponding analytical expressions for initial density and pressure in Appendix B.

2.3. Mesh structure and boundary conditions

At the root level, we set 256 × 256 × 256 active cells in {x, y, z}. The domain extends from −100 to 100 in the three directions. On top of the root level, we added multiple levels of static mesh refinement up to level 8 (see Table 1), which corresponds to a grid spacing $\delta \simeq 0.003a_{\mathrm{b}}^i$. The highest refinement level is concentrated within a box of dimensions $1.6a_{\mathrm{b}}^i\times 1.6a_{\mathrm{b}}^i\times 0.2a_{\mathrm{b}}^i$, centered on the binary's center of mass. We used no inflow boundary conditions in the three directions (“diode” boundary conditions).

3. Results

Here, we describe results of our simulations. In Fig. 1, we show a density cross section of one of our simulations, which illustrates a typical output of our runs. We then address different gravitational softening formulations (Sect. 3.1), gravitational torques (Sect. 3.2), envelope ejection and its morphology (Sect. 3.3), and shear rate (Sect. 3.4).

Fig. 1. Snapshot of the gas density cross section in the xy plane for simulation E005.S.l.q1 after ∼2550 orbits.

3.1. Overview of gravitational softening formulations

In Appendix A we provide a thorough comparison of different gravitational softening formulations. Here, we simply summarize the most important findings. We show that the Ruffert and spline formulations (Eqs. (8) and (9)) yield a more accurate representation of the true gravitational potential of the binary than the Plummer formulation (Eq. (7)), at a distance $\Vert r-r_i\Vert \gtrsim ϵ/5$ from the core. All three formulations underestimate the gravitational potential by more than a factor of 10 at a distance $\Vert r-r_i\Vert \lesssim ϵ/10$.

The accuracy of the softened potential itself, however, is not relevant. Instead, it is the gradient of the potential, which appears in the equations of hydrodynamics, that is crucial for meaningful comparison. In Fig. A.2, we show that, as expected, the three considered softening formulations effectively dampen the gravitational potential gradient within the softening spheres. While both the Plummer and spline methods achieve effective gradient reduction, the Ruffert method results in a more dramatic dampening without necessarily further relaxing the time step constraint. This variation in the degree of gradient suppression within softening spheres among the different considered formulations may have important implications for the outcome of our simulations.

In Appendix A we further show that while the spline-softened potential gradient becomes exactly Newtonian at $\Vert r-r_i\Vert \ge h$, the Plummer softening formulation still significantly underestimates the potential gradient at distance a $\Vert r-r_i\Vert =O(10ϵ)$ from the point mass. Hence, the adopted definition h≡14ϵ/5 is far from implying any form of relevant equivalence between the Plummer and spline formulations, and we suggest an alternative definition of h that implies conditional equivalence in Appendix A. For the sake of consistency with previous studies however, we still adopted h≡14ϵ/5 in this work, while acknowledging its limitations. Finally, we show that an asymptotic regime of small ϵ can be achieved for the total torque and gravitational body force provided any non-axisymmetric structures (e.g., mini-disks in circumbinary disks) significantly contributing to the total torque and body force are properly resolved outside of the softened regions, that is, provided the necessary but insufficient condition (Eq. (A.6)) is satisfied.

3.2. Gravitational torque

Here we address the impact of the gravitational softening formulation and radius and of the intraorbital resolution on the gravitational torque exerted by the binary on the surrounding envelope. The z-directed gravitational torque exerted by the binary on the envelope reads

$$\begin{array}{cc}\hfill {\dot{J}}_{z,\mathrm{grav}}& =-\int s\rho \frac{\partial \mathrm{\Phi }}{\partial \phi }\mathrm{d}V\hfill \\ \hfill & =\int s\rho {\displaystyle \sum _{i=1}^2}{\mathit{GM}}_i\frac{\partial f(\Vert r-r_i\Vert ,ϵ)}{\partial \phi }\mathrm{d}V,\hfill \end{array}$$(10)

where s is the radial cylindrical coordinate, and

$$\frac{\partial f(\Vert r-r_i\Vert ,ϵ)}{\partial \phi }=-\frac{(x_{i}y-y_{i}x)}{{\left({\Vert r-r_i\Vert }^2+{ϵ}^2\right)}^{3/2}},$$(11)

using a Plummer model,

$$\frac{\partial f(\Vert r-r_i\Vert ,ϵ)}{\partial \phi }=-\frac{(x_{i}y-y_{i}x)\left(1-exp\left(-{\Vert r-r_i\Vert }^2/{ϵ}^2\right)\right)}{{\left({\Vert r-r_i\Vert }^2+{ϵ}^{2}exp\left(-{\Vert r-r_i\Vert }^2/{ϵ}^2\right)\right)}^{3/2}},$$(12)

using Ruffert (1993)'s softening method, or, using the spline defined in Eq. (9)

$$\begin{array}{cc}\hfill & \frac{\partial f(\Vert r-r_i\Vert ,ϵ)}{\partial \phi }\hfill \\ \hfill & =\frac{(x_{i}y-y_{i}x)}{h^3}\{\begin{array}{cc}-\frac{32}{3}+\frac{192}{5}{u}^2-32u^3,\hfill & \text{if}\hspace{0.5em}0\le u<\frac{1}{2},\hfill \\ \frac{1}{15u^3}-\frac{64}{3}+48u-\frac{192}{5}{u}^2+\frac{32}{3}{u}^3,\hfill & \text{if}\hspace{0.5em}\frac{1}{2}\le u<1,\hfill \\ -\frac{1}{u^3},\hfill & \text{if}\hspace{0.5em}u\ge 1.\hfill \end{array}\hfill \end{array}$$(13)

In the absence of core accretion, the time derivative of the binary's z-directed angular momentum is ${\dot{J}}_{z,\mathrm{b}}=-{\dot{J}}_{z,\mathrm{grav}}$. Accurately quantifying the gravitational torque is thus crucial for understanding the long-term evolution of the binary system. Achieving this requires careful selection of the gravitational softening formulation and radius, along with sufficient numerical resolution in the binary's vicinity where the torque density is maximum.

3.2.1. Gas distribution near the two cores

In Figs. 2, 3, and 4 we show snapshots of the gas density on the xy and xz planes, as well as gravitational equipotentials, for our q=1 and q=1/3 runs using the spline softening method. We see that isopycnic surfaces tend to coincide with gravitational equipotentials in the binary's vicinity, in other words, ∇ρ×∇Φ≃0. Consequently, the structure around the two cores visually resembles contact binary stars, which share a hydrostatic envelope (e.g., Lucy 1968; Lombardi et al. 2011). To further illustrate the quasi-hydrostatic structure corotating with the two cores, in Fig. 5 we show a snapshot of the normalized deviation from hydrostatic equilibrium, $\Vert \nabla \rho \times \nabla \mathrm{\Phi }\Vert /(\Vert \nabla \rho \Vert \Vert \nabla \mathrm{\Phi }\Vert )$. We see that the equilibrium condition is satisfied in the immediate vicinity of the two cores.

Fig. 2. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼300 orbits for ϵ=0.05, 0.1, 0.2, 0.3, and 0.5 (from left to right) for q=1 with fixed orbital separation and using the spline softening method. The fourth snapshot corresponds to the low-resolution simulation run E02.S.f.q1.vlr. Solid black lines indicate equipotentials of the smoothed binary potential. The dashed black line marks the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

Fig. 3. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼2100 orbits for ϵ=0.05, 0.1, 0.2, 0.3, and 0.5 (from left to right) for q=1/3 with live orbital separation and using the spline softening method. The dashed black line illustrates the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

Fig. 4. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼2100 orbits for ϵ=0.05, 0.2, and 0.5 (from left to right) for q=1 with live orbital separation and using the spline softening method. The dashed black line illustrates the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

Fig. 5. Close-up view of the normalized deviation from hydrostatic equilibrium in the xy plane after 1500 orbits for simulation run E005.S.l.q03. The cross and plus signs respectively indicate the position of the primary's core and of the companion. The white circles mark the extent of the softened spheres of radius ϵ=0.05.

Hydrostatic equilibrium implies ρ=ρ(Φ), that is, ρ has the same symmetry properties as Φ. In particular,

$$\rho (x_i+x,y_i-y,z)=\rho (x_i-x,y_i+y,z).$$(14)

Conversely, ∂Φ/∂φ is antisymmetric with respect to the plane orthogonal to the orbital plane and parallel to the r_i−r_j vector:

$${\frac{\partial \mathrm{\Phi }}{\partial \phi }|}_{(x_i+x,y_i-y,z)}=-{\frac{\partial \mathrm{\Phi }}{\partial \phi }|}_{(x_i-x,y_i+y,z)}.$$(15)

Consequently, the gravitational torque density, sρ∂Φ/∂φ, is also antisymmetric with respect to such a plane, resulting in a near-zero gravitational torque contribution from the quasi-hydrostatic regions. We emphasize that the reduction in gravitational torque arises from the formation of a corotating quasi-hydrostatic equilibrium structure, rather than from simply gas corotation, as is often assumed.

Maintaining the hydrostatic equilibrium near the cores largely depends on the gravitational softening formulation and on the grid resolution. For instance, while the spline and Ruffert softening methods behave very similarly for $\Vert r-r_i\Vert >ϵ$, the Ruffert softened potential gradient is much weaker than the spline-softened potential gradient within the softening spheres (see Fig. A.2). Such a shallower potential can be resolved with fewer cells. However, the potential gradients for both the Ruffert and spline methods reach their maximum values outside the softening radius, at a distance of about ∼1.2ϵ from the core (see Fig. A.2). Hence, the supposed advantage of the less steep Ruffert potential inside the softening sphere may be rendered ineffective, as the region just outside the softening radius, where the potential gradient is at its maximum, likely dictates the grid resolution requirements. In fact, the use of Ruffert's softening method may be problematic because the associated shallower potential within the softening spheres implies weaker hydrostatic support with less concentrated mass, potentially compromising the stability and mass distribution of the gas in the surrounding regions.

The hydrostatic equilibrium of the gas surrounding the binary is a result of mass accumulation about the point masses. Such a mass accumulation and equilibrium prevent the formation of non-axisymmetric structures, such as circumstellar disks, and is permitted due to the absence of core mass accretion in our simulations: it can be circumvented by the use of numerical mass sinks, which we discuss in Sect. 4.3. However, the ability of the primary's core and/or of the secondary star to accrete matter is likely strongly dependent on their nature and evolutionary stage as well as on the thermodynamical and radiative properties of the gas at the core-envelope interfaces (e.g., Hjellming & Taam 1991; Webbink 2008).

3.2.2. Torques with fixed orbital separation

In Fig. 6 we show the gravitational torque exerted by the binary on the envelope, averaged over the last 50 orbits, as a function of ϵ for our simulation runs with q=1 and fixed orbital separations (see Table 1). This allows us to isolate the effect of the softening radius using a spline softening method. We observe that smaller softening radii generally lead to reduced gravitational torque at quasi-steady state.

Fig. 6. Gravitational torque averaged over the last 50 orbital periods as a function of the softening radius, ϵ, for all of our q=1 simulations with fixed orbital separation.

This trend occurs because increasing ϵ increases the critical distance to the cores, beyond which the gravitational potential gradient is significantly underestimated (Eq. (A.2)). As a result, the hydrostatic support of the gas surrounding the binary weakens, making it more vulnerable to perturbations from the surrounding dynamics. These perturbations disrupt the symmetry of the gas distribution, resulting in a net gravitational torque. Therefore, although a larger ϵ reduces the contribution from existing density asymmetries to the gravitational torque (Fig. A.1), it also promotes gas distribution asymmetries by artificially weakening hydrostatic equilibrium. This effect can be seen in the slight misalignment of isopycnic lines and equipotentials near the cores in Fig. 2 for ϵ=0.5.

However, our simulation run E005.S.f.q1 deviates from this trend, indicating that an asymptotic regime of small ϵ is not reached. This is because the maximum refinement level of 7 is insufficient to resolve the pressure gradient within the softening sphere of radius 0.05 at all times. As a result, hydrostatic equilibrium is not consistently maintained, and the gas distribution fails to match the symmetry of the gravitational potential, leading to artificial gravitational torques and impacting the large-scale density structure (see Sect. 3.3). By increasing the maximum refinement level to 8 (run E005.S.f.q1.hr), the simulation produces a quasi-steady gravitational torque and a mass distribution similar to that observed in run E01.S.f.q1, where the softening radius ϵ is set to 0.1.

3.2.3. Torques with “live” orbital separation

We now aim to understand the role of the gravitational softening formulation and radius, as well as of the intraorbital resolution on the secular evolution of binaries’ orbits. In this work, we fixed the orbital eccentricity, e_b, to zero and neglected core accretion. In this simplified setup, the time derivative of the binary's angular momentum can be written as the orbital separation evolution equation

$$\frac{\dot{a_{\mathrm{b}}}}{a_{\mathrm{b}}}=2\frac{{\dot{J}}_{z,\mathrm{b}}}{J_{z,\mathrm{b}}}·$$(16)

We used this equation to evolve the binary semimajor axis directly in the simulation.

In Fig. 7 we show the orbital separation evolution and the orbital contraction timescale averaged over the last 500 orbital periods as a function of the softening radius, for all of our “live” binary simulations (see Table 1). We find that simulations E03.S.l.q03, E02.S.l.q03, E01.S.l.q03, and E005.S.l.q03 yield relatively similar outcome, and that simulations with ϵ≤0.1 result in a final orbital evolution timescale of $\sim {10}^5{P}_{\mathrm{orb}}^i$. For q=1/3, we find that the orbital evolution of our simulation with the largest softening radius (simulation E05.S.l.q03) stands out from the other well-resolved simulations. That is because this simulation does not satisfy criterion (Eq. A.6). In particular, we find that the contribution from the (extended) softening regions to the total gravitational torque exceeds the total gravitational torque measured in the simulation (see Fig. A.4). This simulation thus cannot be used to obtain any physical insight nor to make predictions.

Fig. 7. Panel (a): Orbital separation evolution for all “live” binary simulations (see Table 1). Panel (b): Orbital separation evolution timescale, ${\tau }_{a_{\mathrm{b}}}=\langle \dot{a_{\mathrm{b}}}/a_{\mathrm{b}}{\rangle }_{500N_{\mathrm{orb}}}^{-1}$, averaged between orbits 1625 and 2125 as a function of the softening radius, ϵ.

Simulations E02.S.l.q03 and E02.S.l.q03.lr are initialized identically, except that the binary's vicinity and the intraorbital region are less well resolved in the latter (see Table 1). We find that, for the same softening radius and method, the less well-resolved run results in a final orbital separation approximately 20% larger than its more resolved counterpart, and an orbital contraction timescale about four times larger. We note that despite having lower resolution in the binary's vicinity than simulation E02.S.l.q03, simulation E02.S.l.q03.lr is still likely to be better resolved than most CEE simulations in the literature. Simulation R.l.q03.lr uses a Ruffert softening method with ϵ=3δ with a minimum grid cell width $\delta \simeq a_{\mathrm{b}}^f/5$. These parameters broadly represent the late phases of CEE simulations of Sandquist et al. (1998, 2000), Passy et al. (2012), Staff et al. (2016), and Iaconi et al. (2018), for example. Here, we find a final orbital separation approximately 30% larger and, crucially, an orbital contraction timescale more than ten times larger than in our well-resolved simulations. These results suggest that high spatial resolution is crucial for accurate binary separation evolution, as it is required to capture the relevant spatial scales associated with the physical processes involved in mass redistribution. Achieving a sufficiently small softening radius also requires high resolution to resolve the pressure gradients within and around the softening spheres. Conversely, simulations with q=1 using three different softening radii, where pressure gradients within the softening spheres are well resolved, yield similar final orbital separations and contraction timescales. The final separation is 50% larger compared to binaries with q=1/3, while the orbital contraction timescale is approximately ten times longer, around ${10}^6\hspace{0.5em}{P}_{\mathrm{orb}}^i$.

We further note that the final separation increases for increasing mass ratio. This has been previously noted by Passy et al. (2012) and Nandez & Ivanova (2016). However, unlike our work where the total binary mass M₁+M₂ is conserved as the mass ratio varies, their results were based on a constant primary mass M₁.

Finally, we note that the orbital separation evolution is likely significantly impacted by the absence of core accretion in our simulations. Indeed, mass and angular momentum accretion can have a significant impact on the orbital evolution and require the addition of an extra term to Eq. (16), which would then read

$$\frac{\dot{a_{\mathrm{b}}}}{a_{\mathrm{b}}}=\frac{\dot{M}}{M}\left(2\frac{M{\dot{J}}_{z,\mathrm{b}}}{\dot{M}{J}_{z,\mathrm{b}}}-3\right),$$(17)

where in this case ${\dot{J}}_{z,\mathrm{b}}=-{\dot{J}}_{z,\mathrm{grav}}-{\dot{J}}_{z,\mathrm{accr}}$ and $-{\dot{J}}_{z,\mathrm{accr}}$ is the total z-directed angular momentum accretion rate onto the binary. Accretion processes are typically implemented as numerical mass sinks (Dittmann & Ryan 2021; Dempsey et al. 2022), which can also significantly affect the gravitational torque by hampering the hydrostatic equilibrium of the gas surrounding the binary resulting from mass accumulation (see the discussion in Sect. 4.3).

3.3. Envelope ejection and morphology

We show the evolution of gas mass within the numerical domain¹ for all our “live” simulations in Fig. 8. We see that although the orbital contraction timescale plateaus around 10⁵–10⁶ orbits, most of the gas leaves our numerical domain on a much shorter timescale of a few hundred orbits. For reasonable properties of the system, this timescale is a few decades and is much shorter than the asymptotic giant branch lifetime and inferred ages of (proto)planetary nebulae (Bujarrabal et al. 2001; Chamandy et al. 2020; Sand et al. 2020; Ondratschek et al. 2022; Moreno et al. 2022; Vetter et al. 2024). The softening radius introduces significant nonlinearity in envelope ejection outcomes. A larger softening radius reduces the available gravitational energy that can be transferred from the binary to the shared envelope. Specifically, a larger ϵ decreases the gravitational energy available for envelope ejection by reducing the depth of the binary potential. Conversely, a larger ϵ decreases the binding energy density of the gas near the binary, potentially making the remaining envelope easier to expel. Additionally, changes in orbital separation influence the amount of gravitational energy injected into the envelope. However, the orbital evolution is also directly affected by the mass of the remaining envelope and its spatial distribution. This interplay between the softening radius ϵ and the orbital separation makes predicting the precise impact on envelope ejection challenging.

Fig. 8. Enclosed gas mass for various values of ϵ using the spline-softening formulation.

However, extreme cases are more straightforward to interpret. A large softening radius, such as ϵ=0.5 (e.g., run E05.S.l.q03), significantly limits the total gravitational energy available for envelope ejection, resulting in reduced mass ejection. Similarly, low grid resolution runs (e.g., runs R.l.q03.lr and E02.S.l.q03.lr) fail to capture the relevant spatial scales involved in mass redistribution and do not accurately resolve pressure gradients near the cores. This leads to an inaccurate representation of gravitational forces and energy transfer and, in runs R.l.q03.lr and E02.S.l.q03.lr, in the ejection of (almost) the entire envelope.

In addition to the significant differences in the amount of mass ejected from the computational domain, the large-scale morphology of the remaining envelope strongly depends on the binary's mass ratio, softening radius, and intraorbital numerical resolution. Notably, in all simulations with a mass ratio q=1/3, the gas distribution remains nearly isotropic, regardless of resolution or softening radius. In contrast, simulations with q=1 result in the formation of a circumbinary disk, driven by the higher amount of angular momentum injected into the envelope. This is illustrated in Fig. 9, where we show snapshots of the gas density cross-section in the xz-plane at the end of simulations E005.S.l.q033 and E005.S.l.q1, as well as in Figs. 10, C.1, and C.3.

Fig. 9. Snapshot of the gas density cross section in the xz plane after 2150 orbits for simulations E005.S.l.q033 (panel a) and E005.S.l.q1 (panel b).

Fig. 10. Snapshots of the gas density, radial velocity, normalized Bernoulli parameter, temperature ratio T/Tvir, specific entropy, and Mach number in the xz plane for simulation E05.S.l.q1 after ∼2250 orbits.

Finally, despite the absence of magnetic fields in our simulations, we observe jet-like polar outflows in simulations with q=1, where sufficient angular momentum has been injected into the envelope to clear, at least partially, a low-density polar funnel and form a circumbinary disk. These polar outflows are present for all three simulations with live orbits and q=1. The outflows can alternate between unipolar and bipolar configurations and sometimes take the form of rising hot bubbles when the polar funnel is not yet cleared. This aligns with the findings of Zou et al. (2020), who showed that an injected isotropic wind can be collimated into jet-like polar outflows by the common envelope ejecta, even in the absence of magnetic fields. In our simulations, such outflows are self-consistently launched. In Figs. 10, C.1, and C.3, we show snapshots of the gas density, the radial velocity, and the normalized Bernoulli parameter

$$\begin{array}{cc}\hfill \frac{\mathit{Be}}{|\mathrm{\Phi }|}& =\frac{1}{2}\frac{|u{|}^2}{|\mathrm{\Phi }|}+\frac{\gamma P}{(\gamma -1)\rho |\mathrm{\Phi }|}-1\hfill \\ \hfill & =\frac{1}{2}\frac{|u{|}^2}{|\mathrm{\Phi }|}+\frac{H}{|\mathrm{\Phi }|}-1,\hfill \end{array}$$(18)

where H is the enthalpy. The Bernoulli parameter Be represents the total specific energy of a gas parcel. A negative value of Be indicates that the gas parcel is gravitationally bound to the binary, while a positive value of Be indicates that it is unbound. In these figures, we also show the gas temperature ratio T/T_vir=3P/(ρ|Φ|), the Mach number, and the specific entropy of the gas lnP/ρ^γ.

We find that these polar outflows are characterized by positive Bernoulli parameters, suggesting that they have sufficient kinetic and thermal energy to overcome the binary's gravitational potential, and therefore be effectively unbound. We further find that these polar outflows are very hot compared to the circumbinary disk and have high specific entropy. This suggests that the accreting disk material is strongly heated in the binary's vicinity, and that part of this heating is irreversible. The heated material then may expand along the polar axis, either freely as a pressure-driven jet-like outflow (e.g., Nobuta & Hanawa 1999) or as buoyant hot bubbles. In our in simulations, the main source of gas heating in the binary's vicinity is shock-heating, which can be decomposed into two main contributions, adiabatic compression, −P∇·u, and irreversible heating from viscous dissipation at the shock front, the latter being the source of specific entropy in our simulations.

In our three live binary simulations with a mass ratio of q=1, we observe a predominantly unipolar outflow that persists for most of the simulations’ duration. This outflow arises as part of the kinetic energy from the cold inflowing material, typically coming from the opposite pole, is converted to thermal energy through shock heating as it collides with the hydrostatic structure surrounding the cores. This results in a fast, hot, high-entropy jet-like structure propagating outward along the polar axis (Fig. C.1). The strong vertical pressure gradient resulting from the shock heating in the binary's vicinity can eventually quench polar accretion, and, doing so, limit the polar outflow in the opposite hemisphere (Fig. C.3).

In the absence of strong polar accretion, the existence of a polar outflow depends on the binary's ability to clear a centrifugally supported low-density cavity around itself. This process inhibits mass accretion onto the binary and leads to the generation of additional shocks. These shocks occur where accretion streams, drawn from the inner edge of the disk, collide with the hydrostatic structure surrounding the two cores (e.g., Fig. C.2), as well as at the cavity edge where the streams collide (e.g., Farris et al. 2015). The increase in temperature and entropy associated with these shocks is dependent on the pre-shock Mach number of the accretion streams. This Mach number is directly related to the size of the cavity, the boundary of which corresponds to the location of the centrifugal barrier. In our simulations, the cavity size is linked to the prescribed gravitational softening length around the cores, with larger softening radii resulting in a more underestimated gravitational torque in the vicinity of the binary (see Appendix A). This can be seen by comparing Figs. 11, C.2, and C.4, for instance. These figures also show that, regardless of the softening radius, the disk is off-centered and the cavity is elliptical. This has already been observed in a large number of hydrodynamical simulations of circumbinary disks around equal mass ratio binaries in a circular orbit (e.g., Shi et al. 2012; Penzlin et al. 2022, 2024) and is associated with the disk's eccentricity (see also Gagnier & Pejcha 2023). Such eccentricity is thought to be excited by stream impacts at the cavity edge (e.g., Shi et al. 2012) and/or by resonant Lindblad excitation (e.g., Lubow 1991; Ogilvie 2007), and to be a long-lived mode trapped by the steep density gradient at the cavity edge (Muñoz & Lithwick 2020). The disk eccentricity can in turn excite the orbital eccentricity of the binary (e.g., Papaloizou 2001). Understanding the excitation mechanisms and long-term evolution of orbital eccentricity is crucial, as eccentricity in post-CEE binaries strongly affects both the merger rates and the properties of their gravitational wave signals. Investigating this in more detail should be the focus of future studies.

Fig. 11. Snapshots of the gas density, Mach number, and specific entropy in the xy plane for simulation E02.S.l.q1 after ∼2250 orbits.

Simulation run E05.S.l.q1 represents a special case. The large softening radius of 0.5 leads to the formation of a nearly spherical, quasi-hydrostatic, and gravitationally bound gas bubble surrounding the binary (Sect. 3.2.1). Within this bubble, the gravitational potential gradient is significantly smoothed, which inhibits strong shocks or compression. Consequently, the bubble maintains low entropy and low Mach number. As the binary begins to clear a low-density cavity, the bubble becomes isolated from the disk. Direct collisions between supersonic accretion streams and such an isolated bubble results in the formation of strong shocks at the bubble's boundary, dramatically increasing local heat and specific entropy production. This process generates a hot bipolar outflow that efficiently clears the polar funnel in both hemispheres. The bubble eventually erodes due to mixing with hotter, higher-entropy material in the cavity and destabilization caused by compression from the accretion streams, leading to its dissolution. The collapse of the bubble results in reduced entropy production, weakening the polar outflow and leading the funnel and cavity to partially refill. In simulation run E05.S.l.q1, this bipolar outflow persists from approximately the 1900th to the 2400th orbit. The bubble can eventually reform, potentially leading to the recurrence of jet-like outflows.

While we observe jet-like outflows in our hydrodynamical simulations, high-momentum outflows from protoplanetary nebulae likely cannot be explained without additional energy transfer mechanisms, such as accretion (Blackman & Lucchini 2014) or the transfer of magnetic energy through the work of Lorentz force (Ondratschek et al. 2022; Gagnier & Pejcha 2024; Vetter et al. 2024). In our simulations, we do not measure significant mean kinetic helicity in either envelope hemisphere. Kinetic helicity is considered a source of the α effect in mean-field dynamo theory (Steenbeck et al. 1966; Brandenburg et al. 2012; Moffatt 2014). Its absence suggests that an envelope dynamo – whether of the αΩ, α², or another type – would not be driven by kinetic helicity.

To conclude, both magnetic and nonmagnetic jets are interesting phenomena, but a more detailed study of jet-launching mechanisms is beyond the scope of this work and will be addressed in future studies.

3.4. Shear rate

In the context of ideal MHD, the magnetic energy density evolution equation reads

$$\frac{\partial E_B}{\partial t}=\left(B\otimes B\right):\nabla u-E_B\nabla ·u-\nabla ·\left(E_{B}u\right),$$(19)

where E_B=B·B/2 is the magnetic energy density, and the terms on the right hand side correspond to the change in magnetic energy from the stretching of the magnetic field lines by the fluid elements ($\left(B\otimes B\right):\nabla u$), the changes in magnetic energy due to expansion or compression effects (E_B∇·u), and changes in local magnetic energy from field advection ($\nabla ·\left(E_{B}u\right)$). The first two terms, associated with magnetic energy amplification, are proportional to the local shear rate,

$$\left(B\otimes B\right):\nabla u=\frac{1}{2}{B}_i{B}_j{\dot{\gamma }}_{\mathit{ij}}\mathrm{and}{E}_B\nabla ·u=\frac{1}{4}{B}_i{B}_i{\dot{\gamma }}_{\mathit{jj}},$$(20)

where

$${\dot{\gamma }}_{\mathit{ij}}=\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}$$(21)

is the local shear rate.

In Figs. 12 and 13 we show snapshots of the absolute horizontal shear rate, $|{\dot{\gamma }}_{\mathit{xy}}|$, on the orbital plane after 300 orbital periods for our fixed q=1 runs and ∼650 orbital periods for our “live” q=1/3 runs, respectively. In Fig. 14 we show the absolute shear rate components, averaged within a cubic region of edge length 2 encompassing the binary system, over last 100 orbital periods for our fixed q=1 runs, as a function of the softening radius ϵ. We find an increase in local shear rate with higher grid resolution and smaller softening radius, each effect observed independently while the other parameter is held constant. This shows the importance of sufficiently high resolution and small softening radii to understand angular momentum transport and magnetic field amplification in common envelopes. Increased resolution indeed allows the capture of more detailed turbulence features and velocity gradients, while smaller softening radii render the binary's gravitational body force more accurate. The often lower resolution and larger softening radii commonly used in ab initio CEE simulations may not accurately capture these effects and may thus lead to an incomplete and potentially erroneous depiction of the dynamics, evolution, and outcome of common envelopes. Unfortunately, achieving such high resolution and small softening radii is often computationally prohibitive for ab initio simulations. Of course, the presence of magnetic fields and the associated Lorentz force in our simulation would act back on the fluid, damping the local shear rate, and promoting a more ordered flow.

Fig. 12. Absolute horizontal shear rate in the orbital plane after ∼300 orbits with fixed separation, for q=1 with a spline softening formulation, and for ϵ=0.5, 0.2, 0.1, and 0.05, going from left to right and top to bottom in a clockwise direction. The bottom-left and bottom-center snapshots respectively correspond to simulation runs E02.S.f.q1.lr and E02.S.f.q1.vlr (see Table 1). The cross and plus signs respectively indicate the position of the primary's core and of the companion. The black circles mark the extent of the softened spheres of radius ϵ.

Fig. 13. Absolute horizontal shear rate on the orbital plane after ∼650 orbits, for q=1/3 with a spline softening formulation, and for ϵ=0.5, 0.3, 0.2, 0.1, and 0.05, going from left to right and top to bottom in a clockwise direction. The cross and plus signs respectively indicate the position of the primary's core and of the companion. The black circles mark the extent of the softened spheres of radius ϵ.

Fig. 14. Volume- and time-averaged absolute shear rate within a cube of side length 2ab centered on the binary's center of mass, averaged over the last 50 orbital periods for simulations with q=1 and fixed orbital separation, as a function of the softening radius, ϵ.

4. Summary and discussions

We conducted a detailed investigation of the dynamics of two point-mass cores orbiting inside a shared envelope with application to the post-dynamical inspiral phase of CEE. Special attention was given to the effects of the gravitational potential softening and the numerical resolution of the central region. We find that the spline-softened potential offers the best overall performance, and we used it for most of our simulations (Table 1). In all our simulations the two cores become surrounded by a shared envelope, which is to a high degree hydrostatic and corotating with the orbit (Figs. 2, 3, and 5). We investigated torques on the central binary both when the semimajor axis was held fixed and when it was allowed to respond to simulation events. We find that for sufficiently resolved simulations, the orbital evolution timescales converges to $\sim {10}^5{P}_{\text{orb}}^i$ for q=1/3 and $\sim {10}^6{P}_{\text{orb}}^i$ for q=1 (Figs. 6 and 7). The softening length and resolution affect the efficiency of mass ejection in complex ways (Fig. 8), but higher binary mass ratios typically result in more flattened ejecta (Fig. 9). In certain situations, our simulations develop clear polar outflows escaping through the low-density polar funnel even in the absence of the collimating effect of a magnetic field. These outflows are pressure-driven and result from the intense shock-heating of the gas in the binary's vicinity. Using kinetic helicity as an indicator of α dynamo, we found that large-scale magnetic fields were unlikely to develop through the usual α effect in our simulations. Finally, we used the local shear rate as another indicator of turbulence-driven magnetic field amplification and observed that it increases with higher grid resolution and smaller softening radius. Each of these effects was observed independently while the other parameter was held constant. We find that quantities such as the volume-averaged shear rate – measured within a small central region surrounding the binary – and the orbital evolution timescale appear to converge for $ϵ\le 0.1a_{\text{b}}^i$ and a grid spacing inside the binary's orbit $\delta \le 6\times {10}^{-3}{a}_{\text{b}}^i$. Smaller values of ϵ require correspondingly smaller sizes of δ to ensure the resolution of pressure gradients within the softening spheres. We now discuss the broader implications of our findings for CEE.

4.1. Implications for current simulations

One of the open questions in CEE theory is what drives the transition from the fast dynamical inspiral to the much slower post-dynamical inspiral phase. The key physical effects influencing this transition are typically considered to be a reduction in gas density or a decrease in the Mach number. Each of these factors would decrease the drag felt by the two cores (e.g., Ostriker 1999; Kim & Kim 2007; Röpke & De Marco 2023).

For instance, the 3D hydrodynamic simulations from Reichardt et al. (2019) suggest that the end of the dynamical inspiral coincides with the corotation of the gas (i.e., to zero relative to the Mach number), consistent with the idea of Bondi-Hoyle-Lyttleton accretion (e.g., Edgar 2004), where accreted mass transfers momentum to the cores. This theory assumes a steady, stable, and homogeneous medium and a low-mass perturber. These assumptions are generally not met in the context of CEE. Our simulations reveal a different mechanism. We find that it is not the corotation of the gas that stalls the inspiral, but rather the fact that the gas surrounding the cores reaches a quasi-hydrostatic equilibrium. The associated symmetry properties of the gas distribution imply (almost) no net gravitational torque and lead to the end of the dynamical inspiral phase. This is in line with the results of Kim (2010).

Building upon the findings presented here and in Gagnier & Pejcha (2023, 2024), we propose a coherent picture for the post-dynamical inspiral phase, applicable to systems with high to moderate mass ratios and cores experiencing weak accretion feedback. The dynamical inspiral ends when the cores become enveloped by a corotating, nearly hydrostatic structure resembling a contact binary, which is itself eventually embedded in a low-density cavity. Unless the outer envelope is fully ejected during the dynamical phase, the subsequent evolution is driven by interactions between the central binary and the outer envelope layers. These outer layers cannot maintain corotation with the central binary, giving rise to spiral waves, shocks, and turbulence. Angular momentum transport is primarily driven by local advective torques from the mean flow and Reynolds stresses from turbulence. This transport occurs outward in a disk-like structure around the orbital plane and inward along the polar axis. The variability in mass accretion within the envelope is remarkably similar to that of circumbinary disks. Apart from brief phases of orbital expansion, the binary's orbit shrinks over a timescale of ∼10⁵–10⁶ P_orb, though significant core accretion can shorten this timescale to about 10³ P_orb. Depending on the properties of the remaining envelope, the contraction timescale of the central binary can be either longer or shorter than the envelope's thermal timescale. Consequently, different processes, such as energy input from the contracting binary, dust- and recombination-driven winds, and feedback from core accretion, can more or less contribute to the ultimate ejection of the envelope across different stars. With current 3D simulations, often limited by both their short duration and range of included physics, it remains challenging to assess the relative importance of these processes.

Magnetic field amplification and the potential launching of jets or similar collimated outflows is another topic of importance for theory and observations. Our results on kinetic helicity, shear rate, and our earlier ideal MHD simulations with excised central region tend to suggest that the envelope with a central binary does not efficiently produce large-scale organized magnetic fields. However, this statement has a limited validity due to us not simulating earlier CEE stages, where magnetic field amplification can occur around the plunging companion (Ohlmann et al. 2016b; Ondratschek et al. 2022; Vetter et al. 2024), and because we have not yet performed MHD simulations that explicitly include dynamics around the two cores. At the same time, we have shown that thermally launched polar outflows collimated by the centrifugally evacuated polar funnel can occur for at least some amount of time without the presence of magnetic fields. Our earlier simulations indicate that the existence and properties of the polar funnel can vary considerably due to the turbulence in the envelope (Fig. 3 in Gagnier & Pejcha 2024). This implies that magnetic fields are not a necessary condition for launching at least intermittent time-variable jets, but it is not yet clear whether this is also the case for jets that are stronger and sustained over longer times.

4.2. Relation to contact binaries

The corotating hydrostatic structure surrounding the two cores bears a remarkable similarity to the shared envelope of ordinary contact binaries. Despite decades of observations and theoretical efforts, the internal structure and evolution of contact binaries remain unknown. The fundamental problem is the Kuiper (1941) paradox arising from the observational fact that the two cores, with very different masses and luminosities, are able to maintain a nearly constant effective temperature across the shared surface. Proposals to solve this problem typically involve strong subsurface heat-redistributing flows and the absence of global thermal equilibrium (e.g., Lucy 1968; Shu & Lubow 1981; Stepien 2006).

Based on the analogy with contact binaries, we can now speculate about two potential effects that this shared, corotating hydrostatic envelope can have on the two cores inside CEE. First, the shared envelope can develop internal flows, which can redistribute heat from one core to the other. If this transport is sufficiently efficient and long-lasting, it might be incorrect to treat the cores as single stars, for example by emulating the effect of CEE as a simple rapid envelope removal. Second, after the dynamical inspiral, it is unlikely that the spin of the individual cores would corotate with the new binary orbit. As a result, shear will necessarily develop between the corotating hydrostatic envelope and the two cores, potentially leading to the generation and amplification of magnetic fields. Indeed, contact binaries typically show higher levels of magnetic activity than single stars of the same type (e.g., Vilhu & Walter 1987). Neither of these two effects would be captured by existing 3D simulations due to insufficient resolution.

It is also important to discuss how long the corotating hydrostatic structure can survive around the two cores. If one of the cores is capable of accreting gas, the mass of the envelope will decrease. However, the accretion rate might be very low, for example when one of the cores is a stellar-mass black hole accreting at the Eddington limit, resulting in long accretion timescales. More disruptive is the potential feedback from the accreting core in the form of disk winds or jets. Even if neither of the cores accrete, the shared envelope should eventually disperse due to ordinary stellar winds once the outer envelope is cleared.

Finally, it is interesting to note that the corotating hydrostatic envelope formed very rapidly within the first few tens of orbits. This suggests that similar structures might be formed in events such as rapid mass transfer preceding CEE and thus make binaries appear to be in contact, similar to what was observed in the pre-explosion orbital variability in V1309 Sco and other similar transients (Tylenda et al. 2011; St¸epień 2011; Pejcha 2014; Pejcha et al. 2017; Blagorodnova et al. 2020, 2021).

4.3. Future work

There are several areas where our work can be improved in the future. First, we did not account for accretion onto the cores. This process, often implemented as numerical mass sinks, prevents potentially unrealistic mass accumulation around the point masses by removing excess mass at a set rate. The absence of such sinks hinders the formation of asymmetric structures around the point masses. Their implementation could lead to an increase in the total gravitational torque, potentially resulting in smaller final orbital separations as well as differential accretion toward q=1. Additionally, accretion would complexify the gas dynamics and thus affect magnetic energy amplification. However, whether accretion occurs on the primary, secondary, or both cores, and at what rate, remains highly uncertain. These factors heavily depend on the nature of the cores and the thermodynamic and radiative properties at their interface with the envelope. Furthermore, feedback from the accretion could similarly fundamentally change the gas dynamics; however, resolving these effects would require much higher resolution. Second, our calculations could be repeated in the limit of ideal MHD, shedding further light on the amplification of magnetic fields. Finally, implementing self-gravity of the gas would improve the accuracy of the shared envelope's internal dynamics and its interaction with the outer envelope. This may, in turn, influence the eccentricity of the binary orbit, which could be excited or damped by the effects of the surrounding gas.

Acknowledgments

The research of D.G. and O.P. has been supported by Horizon 2020 ERC Starting Grant ‘Cat-In-hAT’ (grant agreement no. 803158). D.G. acknowledges support by the Klaus-Tschira Foundation. D.G. acknowledges funding by the European Union (ERC, ExCEED, project number 101096243). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254).

References

Appendix A: Plummer versus Ruffert versus Spline

We show the ratio between the non-softened gravitational potential and the Plummer (Eq. 7), Ruffert (Eq. 8), and spline-softened (Eq. 9) gravitational potentials as a function of the softening radius ϵ and the distance to the core in Fig. A.1. This figure also shows the ratio between the non-softened and softened gravitational potential gradients as a function of the softening radius and the distance to the core for the three softening methods. The same quantities for a fixed softening radius (ϵ=0.2) are displayed in Fig. A.3. Using the definition h≡14ϵ/5, the three softening methods yield very similar gravitational potentials with a difference of less than or equal to 1% up to a distance of ϵ/5. Beyond this distance, the Plummer-softened potential converges to Newton's law more slowly than the Ruffert- and spline-softened potentials, with the spline-softened potential becoming exactly Newtonian at a distance greater than or equal to h.

Fig. A.1. Top row: Ratio between the non-softened gravitational potential and the Plummer-, Ruffert-, and spline-softened gravitational potentials, as a function of the softening radius and of the distance to the core. Bottom row: Same as the top row but for the gravitational torque density associated with one of the two binary components. The white full lines indicates $\Vert r-r_i\Vert =ϵ$ and the dashed white line $\Vert r-r_i\Vert =h=14ϵ/5$.

Comparing the softened potentials, however, is meaningless as the equations of hydrodynamics only involve their gradient. The softened gradients are underestimated by orders of magnitude for $\Vert r-r_i\Vert \lesssim ϵ$, with the three softening formulations (see Figs. A.1, A.2, and A.3). Hence, as intended, the binary's body force density is dramatically underestimated within the softening regions. We further note that

$$\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Soft}}\Vert }=\frac{{\dot{\ell }}_{z,i}}{{\dot{\ell }}_{z,i}^{\mathrm{Soft}}}\hspace{0.5em},$$(A.1)

that is the gravitational torque density is equally underestimated. While dramatically underestimated, the Plummer and spline-softened gravitational potential gradients are better at reproducing the true gravitational torque and gravitational force densities than the Ruffert alternative, provided $\Vert r-r_i\Vert \lesssim ϵ$. As for the softened potential, the Plummer-softened potential gradient converges more slowly to Newton's law than the Ruffert and spline-softened potential gradients, with the spline-softened potential gradient becoming exactly Newtonian at $\Vert r-r_i\Vert \ge h$. In particular, one can derive that, using the spline softening formulation, the relative difference

$$\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert -\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Spline}}\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Spline}}\Vert }\{\begin{array}{cc}\le 0.8\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.295ϵ.\hfill \\ \le 0.5\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.448ϵ.\hfill \\ \le 0.1\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.911ϵ.\hfill \\ \le 0.01,\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 2.325ϵ.\hfill \\ =0,\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \ge h=14ϵ/5.\hfill \end{array}$$(A.2)

Equivalently, using the Plummer-softened potential

Fig. A.2. Comparison between the Euclidian norm of the gradient of the different potentials as a function of the distance to the core for a softening radius ϵ=0.2, as indicated by the vertical line.

Fig. A.3. Left: Ratio between the non-softened gravitational potential and the Plummer- and spline-softened gravitational potentials, as a function of the distance to the core, for ϵ=0.2. Right: Same but for the gravitational torque density associated with one of the two binary components. Horizontal black lines indicate a relative difference of 1% (full line) and 10% (dashed line) with the non-softened value.

$$\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert -\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Plummer}}\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Plummer}}\Vert }\{\begin{array}{cc}\le 0.8\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.444ϵ,\hfill \\ \le 0.5\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.795ϵ,\hfill \\ \le 0.1\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 3.904ϵ,\hfill \\ \le 0.01,\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 12.258ϵ,\hfill \\ \to 0,\hfill & \text{as}\hspace{0.5em}\Vert r-r_i\Vert \to \infty ,\hfill \end{array}$$(A.3)

And the Ruffert-softened potential

$$\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert -\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Ruffert}}\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Ruffert}}\Vert }\{\begin{array}{cc}\le 0.8\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.149ϵ,\hfill \\ \le 0.5\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.264ϵ,\hfill \\ \le 0.1\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 1.672ϵ,\hfill \\ \le 0.01,\hfill & \text{if}\hspace{0.5em}\Vert r-r_i\Vert \gtrsim 2.209ϵ,\hfill \\ \to 0,\hfill & \text{as}\hspace{0.5em}\Vert r-r_i\Vert \to \infty .\hfill \end{array}$$(A.4)

Hence, for a given softening radius ϵ, the Ruffert and spline-softened potential gradients converge similarly to the real potential with distance to the core. In particular, Ruffert softened potential gradient is more accurate than the spline-softened one up to $\Vert r-r_i\Vert \simeq 2.459ϵ$. At $\Vert r-r_i\Vert =h$, Ruffert's softened potential gradient's deviation from the real potential is ∼0.047%.

Comparing Eqs. A.2 and A.3, it is clear that the adopted definition h≡14ϵ/5 is far from implying any form of relevant equivalence between the different softening methods. For example, taking ϵ=0.2 implies that a spline-softened potential gradient will be underestimated by less than 1% at a distance $\Vert r-r_i\Vert \gtrsim 0.465$ from the core. On the other hand, for the same value of ϵ, a Plummer-softened potential gradient will be underestimated by less than 1% at a distance $\Vert r-r_i\Vert \gtrsim 2.452$ from the core, that is, far beyond the companion core. Similarly, while the body force density and gravitational torque density is underestimated by more than 10% within $\Vert r-r_i\Vert \lesssim 0.3822$ with a spline softening, is goes up to a distance ∼0.781 with a Plummer softening. One can instead define h by imposing, for instance,

$$\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert -\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Spline}}\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Spline}}\Vert }=\frac{\Vert \nabla {\mathrm{\Phi }}_i\Vert -\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Plummer}}\Vert }{\Vert \nabla {\mathrm{\Phi }}_i^{\mathrm{Plummer}}\Vert }=0.01,$$(A.5)

that is, h′≃14.7612ϵ. This alternative definition becomes a physically relevant equivalence between the two formulations provided any non-axisymmetric structure (e.g., mini-disks in circumbinary disks) significantly contributing to the total torque and body force, are properly resolved outside of the softened regions. In other words, an asymptotic regime of small ϵ can be considered reached for softening radii ∼ 2–5 times larger using a spline-softened potential compared to using a Plummer-softened potential. A necessary but insufficient condition for convergence to an asymptotic regime of small ϵ is, for example,

$$-{\displaystyle {\int }_{\Vert r-r_i\Vert \le \alpha ϵ}}s\rho {\displaystyle \sum _{i=1}^2}\frac{\partial {\mathrm{\Phi }}_i}{\partial \phi }\mathrm{d}V\ll {\dot{J}}_{z,\mathrm{grav}}\hspace{0.5em}.$$(A.6)

where α can be chosen from Eqs. A.2, A.3, and A.4, depending on desired inaccuracy tolerance. We note that if such a condition is not satisfied, the Plummer and spline softening methods may provide a less inaccurate representation of the binary's body force than the Ruffert's potential. However, such simulations cannot be used to obtain any physical insight nor predictions. We check criterion (A.6) for different runs in Fig. A.4. Finally, we note that azimuthal grid asymmetry relative to the plane orthogonal to the orbital plane and parallel to the r_i−r_j vector induces torque artificially, independent of gas density distribution and coordinate system. This anomalous gravitational torque is primarily driven by the gravitational torque density near the cores, hence specifically within softening spheres. This again suggests the importance of satisfying condition (A.6).

Fig. A.4. Simulation average of the ratio between the gravitational torque contribution from the softened region and the total gravitational torque, for various values of ϵ and α, for our live binary simulations, using the spline formulation.

Appendix B: Initial density and pressure profiles

The radial density and pressure profiles outside of the binary orbit (i.e., at r≥max(r₁, r₂)=r_>) are

$$\frac{\rho (r\ge r_{>})}{\rho (r_{>})}={\left(1+C\left(\frac{B}{3}\left(\frac{1}{r^3}-\frac{1}{r_{>}^3}\right)-A\left(\frac{1}{r}-\frac{1}{r_{>}}\right)\right)\right)}^n,$$(B.1)

$$P(r\ge r_{>})=K{\rho }^{\mathrm{\Gamma }}\hspace{0.5em},$$(B.2)

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

$$\begin{array}{cccc}\hfill A& =\mathit{GM}\hspace{0.5em},\hfill & \hfill B& =\frac{3A{a}_{\mathrm{b}}^{2}q}{8(1+q{)}^2}\hspace{0.5em},\hfill \\ \hfill A^{\prime }& =A\left(\frac{1}{R_{*}}-\frac{1}{r_{>}}\right)\hspace{0.5em},\hfill & \hfill B^{\prime }& =\frac{B}{3}\left(\frac{1}{R_{*}^3}-\frac{1}{r_{>}^3}\right)\hspace{0.5em},\hfill \\ \hfill C& =\frac{\kappa -1}{B^{\prime }-A^{\prime }}\hspace{0.5em},\hfill & \hfill K& =\frac{(1-\mathrm{\Gamma })(B^{\prime }-A^{\prime })}{\mathrm{\Gamma }(\kappa -1)\rho (r_{>}{)}^{1/n}}\hspace{0.5em},\hfill \\ \hfill {\kappa }^n& =\rho (R_{*})/\rho (r_{>})\hspace{0.5em}.\hfill & \hfill & \hfill \end{array}$$(B.3)

The ρ(r_>) is calculated from the prescribed mass of the envelope outside of the orbit, M_env.

Unlike Gagnier & Pejcha (2023, 2024), we did not excise an inner sphere containing the binary. We therefore needed to derive the radial density and pressure profiles inside of the binary orbit. The non-softened gravitational potential of the two point-mass objects at r≤min(r₁, r₂)=r_< reads

$$\mathrm{\Phi }(r\le r_{<})=-4\pi {\displaystyle \sum _{i=1}^2}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{m=-2k}^{2k}}\frac{G{M}_i}{4k+1}\frac{r^{2k}}{r_i^{2k+1}}(Y_{2k}^m({\theta }_i,{\varphi }_i){)}^{*}{Y}_{2k}^m(\theta ,\varphi )\hspace{0.5em},$$(B.4)

where $Y_{\ell }^m$ are the usual normalized scalar spherical harmonic functions. The latitude- and time-averaged binary potential for spherical harmonics degrees ℓ≤2 (or equivalently k≤1) within the inner binary orbit (r≤r_<) reads

$$\langle \mathrm{\Phi }(r\le r_{<}){\rangle }_{\theta }=-A\left(\frac{1+q^2}{a_{\mathrm{b}}q}-B^{\prime \prime }{r}^2\right)\hspace{0.5em},$$(B.5)

where 〈·〉_θ indicates a time and latitude average, and

$$B^{\prime \prime }=\frac{{\left(1+q\right)}^2}{8a_{\mathrm{b}}^3}\left(\frac{1+q^4}{q^3}\right)\hspace{0.5em}.$$(B.6)

Injecting Eq. B.5 into the equation of hydrostatic equilibrium,

$$\nabla P=-\rho \nabla \mathrm{\Phi }\hspace{0.5em},$$(B.7)

yields

$$\frac{\rho (r\le r_{<})}{\rho (r_{<})}={\left(1+B^{\prime \prime \prime }\left(r^2-r_{<}^2\right)\right)}^n,$$(B.8)

$$P(r\le r_{<})=K{\rho }^{\mathrm{\Gamma }}\hspace{0.5em},$$(B.9)

where

$$B^{\prime \prime \prime }=\frac{B^{\prime \prime }(\kappa -1)}{B^{\prime }-A^{\prime }}{\left(\frac{\rho (r_{>})}{\rho (r_{<})}\right)}^{1/n}\hspace{0.5em}.$$(B.10)

In simulations initialized with q=1 (and e=0), $r_{<}=r_{>}=a_{\mathrm{b}}^i/2$. However, when we set q<1, we further need to compute the density and pressure profiles in the spherical shell located between the orbits of the two point masses. The ℓ≤2 latitude and time averaged binary potential at r_<≤r≤r_> reads

$$\langle \mathrm{\Phi }(r_{<}\le r\le r_{>}){\rangle }_{\theta }=$$(B.11)

$$-A\left(\frac{q}{a_{\mathrm{b}}}+\frac{1}{(1+q)r}-\frac{\mathit{qB}}{3(1+q){\mathit{Ar}}^3}-\frac{B^{\prime \prime }{q}^4{r}^2}{2(1+q^4)}\right)\hspace{0.5em}.$$

Injecting Eq. B.11 into the equation of hydrostatic equilibrium Eq. B.7 finally yields

$$\frac{\rho (r_{<}\le r\le r_{>})}{\rho (r_{>})}=(1+\frac{\kappa -1}{B^{\prime }-A^{\prime }}(\frac{\mathit{qB}}{3(1+q)}\left(\frac{1}{r^3}-\frac{1}{r_{>}^3}\right)-$$(B.12)

$$\begin{array}{cc}\hfill & \frac{A}{1+q}\left(\frac{1}{r}-\frac{1}{r_{>}}\right)+\frac{A{B}^{\prime \prime }{q}^4}{2(1+q^4)}\left(r^2-r_{>}^2\right))),\hfill \\ \hfill P(r_{<}\le r\le r_{>})& =K{\rho }^{\mathrm{\Gamma }}\hspace{0.5em}.\hfill \end{array}$$(B.13)

In Fig. B.1 we show the initial density and pressure profiles for q=1, q=2/3, and q=1/3. We note that while we have ensured the continuity of ρ and P throughout the envelope, the first derivative of Green's function, and hence the gravitational potential gradient, is discontinuous at r=r_> and r=r_<. This naturally leads to discontinuities in the initial density and pressure gradients at these radii, which slightly disrupt hydrostatic equilibrium. We find such perturbations to be sufficiently weak to not require special attention.

Fig. B.1. Initial density and pressure profiles of the envelope for q=1 (full lines) and q=1/3 (dashed lines), obtained from Eqs. B.1, B.8, and B.12. The vertical lines represent the values of r<=min(r1, r2) and r>=max(r1, r2) for q=1 (full lines) and for q=1/3 (dashed lines). When q=1, r1=r2, and thus r>=r<. See Eq. 1 for an example how to scale plotted quantities to physical systems.

Appendix C: Extra figures for Sect. 3.3

Fig. C.1. Same as Fig. 10 but for simulation run E02.S.l.q1.

Fig. C.2. Same as Fig. 11 but for simulation run E02.S.l.q1.

Fig. C.3. Same as Fig. 10 but for simulation run E02.S.l.q1.

Fig. C.4. Same as Fig. 11 but for simulation run E005.S.l.q1.

All Tables

All Figures

	Fig. 1. Snapshot of the gas density cross section in the xy plane for simulation E005.S.l.q1 after ∼2550 orbits.
In the text

Fig. 2. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼300 orbits for ϵ=0.05, 0.1, 0.2, 0.3, and 0.5 (from left to right) for q=1 with fixed orbital separation and using the spline softening method. The fourth snapshot corresponds to the low-resolution simulation run E02.S.f.q1.vlr. Solid black lines indicate equipotentials of the smoothed binary potential. The dashed black line marks the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

In the text

Fig. 3. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼2100 orbits for ϵ=0.05, 0.1, 0.2, 0.3, and 0.5 (from left to right) for q=1/3 with live orbital separation and using the spline softening method. The dashed black line illustrates the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

In the text

Fig. 4. Zoomed-in snapshot of the gas density cross section in the xy (top row) and xz planes (at y=0, bottom row) after ∼2100 orbits for ϵ=0.05, 0.2, and 0.5 (from left to right) for q=1 with live orbital separation and using the spline softening method. The dashed black line illustrates the intersection between the orbital plane and the plane orthogonal to the orbital plane and parallel to the ri−rj vector.

In the text

	Fig. 5. Close-up view of the normalized deviation from hydrostatic equilibrium in the xy plane after 1500 orbits for simulation run E005.S.l.q03. The cross and plus signs respectively indicate the position of the primary's core and of the companion. The white circles mark the extent of the softened spheres of radius ϵ=0.05.
In the text

	Fig. 6. Gravitational torque averaged over the last 50 orbital periods as a function of the softening radius, ϵ, for all of our q=1 simulations with fixed orbital separation.
In the text

	Fig. 7. Panel (a): Orbital separation evolution for all “live” binary simulations (see Table 1). Panel (b): Orbital separation evolution timescale, ${\tau }_{a_{\mathrm{b}}}=\langle \dot{a_{\mathrm{b}}}/a_{\mathrm{b}}{\rangle }_{500N_{\mathrm{orb}}}^{-1}$, averaged between orbits 1625 and 2125 as a function of the softening radius, ϵ.
In the text

	Fig. 8. Enclosed gas mass for various values of ϵ using the spline-softening formulation.
In the text

	Fig. 9. Snapshot of the gas density cross section in the xz plane after 2150 orbits for simulations E005.S.l.q033 (panel a) and E005.S.l.q1 (panel b).
In the text

	Fig. 10. Snapshots of the gas density, radial velocity, normalized Bernoulli parameter, temperature ratio T/Tvir, specific entropy, and Mach number in the xz plane for simulation E05.S.l.q1 after ∼2250 orbits.
In the text

	Fig. 11. Snapshots of the gas density, Mach number, and specific entropy in the xy plane for simulation E02.S.l.q1 after ∼2250 orbits.
In the text

Fig. 12. Absolute horizontal shear rate in the orbital plane after ∼300 orbits with fixed separation, for q=1 with a spline softening formulation, and for ϵ=0.5, 0.2, 0.1, and 0.05, going from left to right and top to bottom in a clockwise direction. The bottom-left and bottom-center snapshots respectively correspond to simulation runs E02.S.f.q1.lr and E02.S.f.q1.vlr (see Table 1). The cross and plus signs respectively indicate the position of the primary's core and of the companion. The black circles mark the extent of the softened spheres of radius ϵ.

In the text

Fig. 13. Absolute horizontal shear rate on the orbital plane after ∼650 orbits, for q=1/3 with a spline softening formulation, and for ϵ=0.5, 0.3, 0.2, 0.1, and 0.05, going from left to right and top to bottom in a clockwise direction. The cross and plus signs respectively indicate the position of the primary's core and of the companion. The black circles mark the extent of the softened spheres of radius ϵ.

In the text

	Fig. 14. Volume- and time-averaged absolute shear rate within a cube of side length 2ab centered on the binary's center of mass, averaged over the last 50 orbital periods for simulations with q=1 and fixed orbital separation, as a function of the softening radius, ϵ.
In the text

Fig. A.1. Top row: Ratio between the non-softened gravitational potential and the Plummer-, Ruffert-, and spline-softened gravitational potentials, as a function of the softening radius and of the distance to the core. Bottom row: Same as the top row but for the gravitational torque density associated with one of the two binary components. The white full lines indicates $\Vert r-r_i\Vert =ϵ$ and the dashed white line $\Vert r-r_i\Vert =h=14ϵ/5$.

In the text

	Fig. A.2. Comparison between the Euclidian norm of the gradient of the different potentials as a function of the distance to the core for a softening radius ϵ=0.2, as indicated by the vertical line.
In the text

Fig. A.3. Left: Ratio between the non-softened gravitational potential and the Plummer- and spline-softened gravitational potentials, as a function of the distance to the core, for ϵ=0.2. Right: Same but for the gravitational torque density associated with one of the two binary components. Horizontal black lines indicate a relative difference of 1% (full line) and 10% (dashed line) with the non-softened value.

In the text

	Fig. A.4. Simulation average of the ratio between the gravitational torque contribution from the softened region and the total gravitational torque, for various values of ϵ and α, for our live binary simulations, using the spline formulation.
In the text

	Fig. B.1. Initial density and pressure profiles of the envelope for q=1 (full lines) and q=1/3 (dashed lines), obtained from Eqs. B.1, B.8, and B.12. The vertical lines represent the values of r<=min(r1, r2) and r>=max(r1, r2) for q=1 (full lines) and for q=1/3 (dashed lines). When q=1, r1=r2, and thus r>=r<. See Eq. 1 for an example how to scale plotted quantities to physical systems.
In the text

	Fig. C.1. Same as Fig. 10 but for simulation run E02.S.l.q1.
In the text

	Fig. C.2. Same as Fig. 11 but for simulation run E02.S.l.q1.
In the text

	Fig. C.3. Same as Fig. 10 but for simulation run E02.S.l.q1.
In the text

	Fig. C.4. Same as Fig. 11 but for simulation run E005.S.l.q1.
In the text