SIMULATING THE COMMON ENVELOPE PHASE OF A RED GIANT USING SMOOTHED-PARTICLE HYDRODYNAMICS AND UNIFORM-GRID CODES

Although they are meant to simulate similar astrophysical situations, high-order Eulerian grid codes and Lagrangian smoothed-particle hydrodynamics (SPH) codes differ fundamentally, with each having advantages and disadvantages. Among other studies, Davies et al. (1993), Frenk et al. (1999), Agertz et al. (2007), Tasker et al. (2008), and Heitsch et al. (2011) aim at identifying these differences. On the one hand, high-order Eulerian grid codes have a better wavenumber resolution than SPH codes for an equal number of cells and particles and are more accurate at resolving the rarefied regions since, unlike SPH, the resolution does not depend on the density of the gas; Eulerian codes also better resolve shocks (Tasker et al. 2008) compared to SPH codes; and finally, SPH noise dominates subsonic flows and therefore makes it difficult for SPH codes to follow perturbations in flows with Mach numbers below unity. On the other hand, SPH codes do not diffuse material properties, and inherently conserve mass, momentum, and energy (Rosswog 2009). While the treatment of boundary conditions can be challenging in grid-based codes when the flow expands beyond the computational domain, SPH easily handles vacuum conditions. It is still unclear which method is the most appropriate to simulate CE interactions. Therefore, we use both methods and confront the results from both codes in order to draw conclusions about their physical relevance.

Both codes solve the fully compressible hydrodynamics equations with self-gravity included. These equations can be written using an Eulerian formulation:

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \mathbf {v}) = 0 \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial \mathbf {v}}{\partial t} + (\mathbf {v} \cdot \nabla) \mathbf {v} = - \frac{1}{\rho } \nabla p - \nabla \Phi \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial u}{\partial t} + \mathbf {v} \cdot \nabla u = - \frac{1}{\rho } \nabla \cdot (p\mathbf {v}) - \mathbf {v} \cdot \nabla \Phi \end{equation} \tag{ 3 }$$

$$\begin{equation} u_{{\rm int}} = \frac{1}{\gamma - 1} \frac{p}{\rho } \end{equation} \tag{ 4 }$$

$$\begin{equation} \Delta \Phi = 4 \pi G \rho, \end{equation} \tag{ 5 }$$

where $\rho, \mathbf {v}, p, \Phi, u, u_{{\rm int}}, {\rm and}\ \gamma$ are the density, velocity, pressure, gravitational potential, specific total energy, specific internal energy, and adiabatic index of the gas, respectively. The total energy is the sum of the internal energy and the macroscopic kinetic energy:

$$\begin{equation} u = u_{{\rm int}} + \mathbf {v}^2/2. \end{equation} \tag{ 6 }$$

Equations (1), (2), and (3) express mass continuity, conservation of momentum, and conservation of energy, respectively. Both codes evolve the internal energy rather than the total energy. An ideal gas equation of state (Equation (4)) for a monoatomic gas (γ = 5/3) closes the system composed by Equations (1)–(3). Such an equation of state represents an adequate approximation of the deep convective envelope of RGB stars (Hjellming & Webbink 1987) although it ignores some physical processes such as radiation pressure and ionization. We discuss this point in detail in Section 5.2.2. Finally, the gravitational potential is calculated using the Poisson equation (Equation (5)).

Enzo is a three-dimensional, adaptive mesh refinement hybrid (hydrodynamics + N-body) grid-based code (Bryan et al. 1995; O'Shea et al. 2004) that we use in uniform-grid mode only. It is primarily designed to simulate cosmological structure formation (Norman et al. 2007). However, its numerous features make it useful for reproducing many different astrophysical situations, including CE interactions.

The Euler equations (Equations (1)–(3)) are solved using the van Leer (1977) second-order advection method also implemented in Zeus (Stone & Norman 1992). Although those equations can also be solved in Enzo by a third-order piecewise parabolic method that resolves shocks and turbulence better, our tests show that it slows down the computation by a factor of two. As we will point out in Section 4, there are neither strong shocks nor important turbulence in our simulations so we favor efficiency and use the van Leer solver. The Poisson equation is solved using fast-Fourier transforms.

In the case of a CE interaction between an RGB star and an MS companion, the radius of the secondary—typically 0.5 R_☉—is small compared to the primary's radius (∼100 R_☉), so we can legitimately model the companion as a point mass particle. Furthermore, as shown in Figure 3, the primary's core is also small (∼0.01 R_☉) and dense, so it can also be modeled as a point mass.

Enzo usually models collisionless particles as a continuous mass field appropriate for computing the gravitational potential in the case where each particle represents many actual particles, such as in cosmological simulations with dark matter. In that case their mass is deposited in the eight nearest cells and added to the gas density of those cells to find the total density for use in solving the Poisson equation (Equation (5)). In a simple two-body interaction between 1 M_☉ and 0.1 M_☉ objects in a 1 year circular orbit without gas, this method does not provide the accuracy required by our problem because of the spreading out of the mass of the point source, leading to an inaccurate gravitational potential. Indeed, a 1% error in the orbit is reached after only six orbits. Consequently, we implemented, as a new type of particle, point mass particles. These particles create a potential that is added analytically to the gas potential calculated using the Poisson equation. Using an analytic potential yields an accuracy of the orbit more than two orders of magnitude better than with the default particles. The gravitational potential created by a point mass particle is smoothed according to the prescription of Ruffert (1993), used in Sandquist et al. (1998):

$$\begin{equation} \Phi _{{\rm PM}}(r) = \frac{-GM_{{\rm PM}}}{ \sqrt{r^2 + \epsilon ^2 \delta ^2 \exp {[-r^2/(\epsilon \delta)^2]}}}, \end{equation} \tag{ 7 }$$

where M_PM is the mass of the particle, r is the distance from the particle, δ is the size of a cell, and = 1.5. The point mass particles are advanced using a leapfrog algorithm. Time stepping is determined by taking the minimum time step between the Courant conditions for the gas, the particles, and the acceleration field:

$$\begin{equation} \delta t_{{\rm gas}} = \min _{\rm cells}\left(\frac{C_1 \delta }{c_s + \max (|v_x|,|v_y|,|v_z|)}\right), \end{equation} \tag{ 8 }$$

$$\begin{equation} \delta t_{{\rm part}} = \min _{\rm particles}\left(\frac{C_2 \delta }{\max (|V_x|,|V_y|,|V_z|)}\right), \end{equation} \tag{ 9 }$$

$$\begin{equation} \delta t_{{\rm accel}} = \min _{\rm cells}\left(\sqrt{\frac{\delta }{\max (|g_x|,|g_y|,|g_z|)}}\right), \end{equation} \tag{ 10 }$$

where C₁ = 0.4 is the Courant factor, C₂ = 0.4 is the particle Courant factor, c_s is the sound speed, $\mathbf {v} = (v_x,v_y,v_z)$ is the velocity of the gas, $\mathbf {V} = (V_x,V_y,V_z)$ is the velocity of a particle, and $\mathbf {g} = (g_x,g_y,g_z)$ is the acceleration field.

Finally, we remark that the current Enzo Poisson solver prevented us from using nested or adaptive grids that would have allowed us to increase resolution locally. The inaccurate treatment of boundary conditions within the refined grids prevented us from stabilizing the RGB progenitor in a multi-grid initial setup. We are currently developing a new Poisson solver that will allow us to use nested grids as well as adaptive mesh refinement and carry out better-resolved simulations.

SNSPH (Fryer et al. 2006) is a three-dimensional, parallel SPH code using tree gravity. It uses a regular Monaghan cubic spline kernel (Monaghan 1992). For the artificial viscosity we use the sum of a bulk viscosity and a von Neumann and Richtmyer viscosity (Rosswog 2009). The particles are organized into a parallel hashed oct-tree as described in Warren & Salmon (1993). The gravitational potential of an SPH particle, i, is smoothed using the following formula:

$$\begin{eqnarray} &&\Phi _i (x_i = r_i/h_i)\nonumber\\ && = \left\lbrace \begin{array}{@{}ll@{}}-G m_i/h_i \times \big(\frac{2}{3} x_i^3 - \frac{3}{10} x_i^4 + \frac{1}{10} x_i^5 - 1.4\big) & {\rm if} \ 0 \le x \le 1 \\ \ms -G m_i/r_i \times \big[\big(\frac{4}{3} x_i^2 - x_i^3 + \frac{3}{10}x_i^4 &\\ \ms \quad -\, \frac{1}{30} x_i^5 - 1.6\big)\big/h_i + 1/15r_i \big] & {\rm if} \ 1 \le x \le 2 \\ \ms -G m_i/r_i & {\rm otherwise} \end{array}, \right. \nonumber\\ \end{eqnarray} \tag{ 11 }$$

where h_i, m_i, and r_i are the smoothing length, the particle mass, and the distance from the particle, respectively. We compare both numerical potentials to the theoretical potential in Figure 1. For a given smoothing length, h_i, the Monaghan (1992) potential used in our SNSPH simulations is deeper than the Ruffert (1993) one used in our Enzo simulations. Also, the Monaghan (1992) potential is exact at distances larger than 2h_i whereas the Ruffert (1993) potential only asymptotically tends to the exact potential.

Zoom In Zoom Out Reset image size

SNSPH uses the fast multipole method to calculate gravitational accelerations (Warren & Salmon 1993). The SPH particles are also advanced using a leapfrog algorithm. Finally, in order to keep the same overall spatial coverage, the smoothing length varies according to the formula from Benz (1989):

$$\begin{equation} \frac{h_i(t)}{h_i(0)} = \left(\frac{\rho _i(0)}{\rho _i(t)}\right)^{1/3}. \end{equation} \tag{ 12 }$$

There is no ideal way to compare the resolution between SPH and uniform-grid codes. However, a few criteria can give us a general idea of how to relate them.

As mentioned by Davies et al. (1993), a first global criterion would be to compare the total number of SPH particles N_part with the total number of cells originally inside the progenitor:

$$\begin{equation} N_{{\rm cells}} = \frac{V_1}{V_G} \times N_{{\rm tot}} \sim 4.19 \times \left(\frac{R_1}{L}\right)^3 N_{{\rm tot}}, \end{equation} \tag{ 13 }$$

where N_tot, V₁, V_G, R₁, and L are the total number of cells, the volume of the primary, the volume of the grid, the radius of the primary, and the linear dimension of the grid, respectively. As time goes by, the gas will however fill a larger fraction of the numerical grid and thus increase the number of relevant cells, but not the real resolution of the simulation.

A more local criterion is to compare the size of an Enzo grid cell, δ, with the SPH smoothing length, which varies in space and time. Indeed, if the companion does not sink much into the primary's envelope and does not modify the inner part of the smoothing length distribution too much, then the resolution deep inside the progenitor does not matter. Therefore, we compare the smoothing length distribution of the SPH model to the cell size of the Eulerian grid. As shown in Figure 2, the smoothing length at small radii does not vary, so an Enzo run with a 128³ grid will be underresolved compared to our canonical 500,000 (roughly 80³) particle SPH run no matter how deep the companion penetrates while a run with a 256³ grid would be equivalent to our SPH runs if the separation between the primary core and the companion always exceeds 20 R_☉. This local criterion for the resolution is not perfect either since it does not take into account the variation of the smoothing length throughout the SPH simulation.

Zoom In Zoom Out Reset image size

Again, comparing the resolution between uniform-grid and SPH codes is quite challenging and both methods have, in the situation we are interested in, strengths and weaknesses: SPH will underresolve the low-density outer parts of the envelope, where the smoothing length dramatically increases, while it will be more accurate in the later phase of the evolution when the separation between the primary's core and the secondary will typically sink below a few cells. Therefore, comparing SPH and grid-based simulations is paramount in order to state which one is more adapted to our problem, and the combination of both global and local criteria is the best way to compare the resolutions of both methods.

As explained in Section 3, the companion is placed at the surface of the primary. Thus, the primary extends beyond its Roche lobe and unstable mass transfer starts immediately. The companion, surrounded by stellar matter, exchanges momentum and energy with this gas through drag. The orbital separation shrinks on a dynamical timescale and its evolution for the 256³ Enzo simulations is shown in Figure 4. Although the orbit is initially circular, it quickly develops eccentricity due to the geometry of the gas ejection. In order to define quantitatively the end of the rapid infall phase and the final orbital separation ad hoc, we consider the evolution of the orbital decay (Figure 5). As expected, the orbital decay is initially quite high (∼0.01 day⁻¹), decreases as less gas is available for the companion to exchange energy with, and eventually reaches a plateau. We decide to define the end of the rapid infall phase to be at the start of this plateau, which occurs at about 280 days for the 0.6 M_☉ companion (Figure 5). All the simulations show the same trend and the lighter the companion, the deeper it falls, and the longer it needs to reach its final orbital separation. The duration of the rapid infall phase is 260, 280, 280, 300, and 340 days, for the 0.9, 0.6, 0.3, 0.15, and 0.1 M_☉ companion, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As orbital energy is transferred to the envelope, the latter is ejected, initially in the orbital plane; at later phases there is an almost equal distribution of matter into the polar direction as well (Figure 6). Overall, almost 90% of the envelope is ejected within an angle of 30° on each side of the equatorial plane. We compare the orbital velocity of the companion (Figure 7) with the local sound speed of the gas (Figure 3, bottom left panel). The former does not exceed 50 km s⁻¹ while the highest sound speed encountered is about 60 km s⁻¹. The companion moves only slightly above or below the local sound speed. We therefore conclude that the SPH noise could not significantly influence the solution. Also, since the motion of the companion is not highly supersonic, the shocks are not strong and we can use Enzo with the faster Zeus solver.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unlike the SPH computational domain, the Enzo grid is spatially limited. Thus, the evolution of the gas that leaves the grid cannot be followed. Therefore, we use the SPH2 simulation to study the global evolution of the angular momentum and the energy of the system.

We compute the angular momentum using the center of mass of the SPH particles as the center of reference. As shown in Figure 8 for the 0.6 M_☉ companion case, the total angular momentum of the system is conserved to less than 1%. Since the ejection of the gas is asymmetric, the center of reference is eventually located outside the orbit. Consequently, studying the orbital components individually is irrelevant as the sign of each component changes during a single orbit. Therefore, we study their sum J_orb instead. During the first 50 days, angular momentum from the orbit almost equally spins up the envelope and unbinds mass from the outer layers. Later on, additional mass no longer gets unbound (see Section 5.1.2) and the angular momentum lost from the orbit spins up the bound envelope only. Since the unbound mass is located at large distances from the primary's core, there is no longer exchange of angular momentum between the unbound mass and the rest of the system. After ∼150 days, there is no longer angular momentum exchange in the system. The primary's core and the companion—which are the main contributors to the calculation of the center of mass—switch positions twice per orbit, which leads to small periodic motions of the center of mass. These periodic displacements are the causes for the small angular momentum fluctuations of the orbital components and the bound mass occurring after 100 days.

Zoom In Zoom Out Reset image size

We plot the various energy components in Figure 9. We start by explaining the different components of potential, thermal, and kinetic energies represented and how they are computed. Among numerous other attributes, each particle i possesses a specific gravitational potential energy ϕ_i, a specific thermal energy u_i, and a specific macroscopic kinetic energy k_i. By definition,

$$\begin{equation} \phi _i = \sum _{j \ {\rm particles}, \ j \ne i} -G \frac{M_j^{\rm grav}}{r_{ij}}, \end{equation} \tag{ 15 }$$

where G is the gravitational constant, M^grav_j is the gravitational mass of particle j, and r_ij is the distance between particles i and j. We compute these different components using the gravitational mass of the particle for the gravitational potential energy and the macroscopic kinetic energy, and the SPH mass for the thermal energy (see Section 3):

$$\begin{eqnarray} &&\Phi _i = M_i^{\rm grav} \phi _i \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&K_i = M_i^{\rm grav} k_i \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&U_i = M_i^{\rm sph} u_i, \end{eqnarray} \tag{ 18 }$$

where M^sph_i is the SPH mass of particle i used to compute its acceleration due to pressure. We recall that both masses are identical for all particles except the primary's core and the secondary. Thus, the total gravitational potential energy of the system is

$$\begin{equation} \Phi _{\rm tot} = \frac{1}{2} \sum _{i \ {\rm particles}} \Phi _i. \end{equation} \tag{ 19 }$$

Zoom In Zoom Out Reset image size

Finally, we subtract the contribution of the secondary from the total potential energy in order to calculate the binding energy of the envelope:

$$\begin{equation} \Phi _{\rm env} = \Phi _{\rm tot} - M_2^{\rm grav} \phi _2, \end{equation} \tag{ 20 }$$

where the subscript "2" stands for the secondary.

During the first 200 days when most of the inspiral happens, the total internal energy of the system decreases by more than a factor of two: the envelope expands and therefore cools. The energy released is transferred mostly into macroscopic kinetic energy of the gas: the envelope is lifted up, accelerated, and the outermost part of the envelope becomes unbound in the first 50 days. At later times, more energy is transferred from the orbit to the envelope but no more material becomes unbound. One can easily note in Figure 9 how the variations of the orbital energy of the core–secondary system and of the total energy of the envelope balance each other. The total energy of the envelope remains negative throughout the simulation. We follow the evolution of the unbound particles and determine their initial position in the envelope. Figure 10 shows the cumulative mass of the particles that will eventually get unbound as a function of their initial distance from the core. It confirms that the unbound mass was initially located in the outer part of the envelope and that almost all gas located initially closer than 40 R_☉ from the primary's core remains bound at the end of the simulations.

Zoom In Zoom Out Reset image size

The fact that a code solves the equations in an accurate and precise way in a particular situation does not necessarily mean it will do so in another regime. Thus, a direct comparison of simulations of the CE interaction using two different numerical methods is a good solution for testing the ability of the two methods to model this problem. One can see in Figure 11 and Table 1 that for each binary system, the final separations in the Enzo simulations are very close to those obtained with the equivalent SNSPH simulations. We may then compare the mass evolution of the material in the volume defined by the Enzo grid, the matter within the initial volume of the primary, and within the current separation. For the 0.6 M_☉ companion (Figure 12), both the mass within the Enzo grid and the mass within the initial volume of the progenitor agree well between the Enzo and the SNSPH runs. For the mass within the orbit we notice a difference of ∼10⁻² M_☉ between the Enzo and the SNSPH runs. This difference is large compared with the mass of an SPH particle (∼10⁻⁶ M_☉) and is due to how accurately accretion of the gas by the core and the companion is resolved by the two codes. We have plotted, in Figure 13, density profiles at different times along the line joining the primary core and the secondary, for the three simulations with the 0.9 M_☉ companion. Accretion onto the secondary is better resolved in the SPH simulations in which the maximum density of the matter accreted by the companion is about 10⁻³ g cm⁻³. This maximum value depends on the resolution of the runs. In the single-grid Enzo runs, accretion is poorly resolved due to the low number of cells resolving the local region around each particle. Although mass is still accreted around the particles, it eventually becomes dispersed. For the SPH runs, around 60 particles interact within a smoothing length so the accretion zone is well resolved. On the other hand, the cell width of the Enzo 256³ runs is about 1.6 R_☉ so the accretion zone cannot be resolved although it is still better than for the Enzo 128³ simulations as can be seen from comparing the different density profiles at 50 days (Figure 13). However, the accurate simulation of accretion onto the secondary is not crucial for the global evolution of the system: as we mentioned earlier, the evolution is not driven by accretion but by drag forces. Although the density of the matter accreted by the companion differs by up to three orders of magnitudes between the two methods, the accreted mass is negligible compared with the companion mass and the final orbital separations are very similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Ricker & Taam (2008) used the FLASH code (Fryxell et al. 2000) to study the CE evolution of a binary system consisting of a 1.05 M_☉ RGB star having a 0.36 M_☉ core and a 0.6 M_☉ companion. Their implementation is somewhat different from ours since they treat the red giant core and the companion as spherical clouds of particles. In spite of those differences, their progenitor is almost identical to ours and they find a final separation of 20 R_☉ which falls within the range of the results given by our simulations SPH2, Enzo2, and Enzo7. Moreover, one can see in Figure 7 that for the 0.6 M_☉ companion, the velocity of the companion stays below 50 km s⁻¹ and therefore, the gas flows are subsonic except in the outer layers. This conclusion was also reached by Ricker & Taam (2008).

In order to determine the sensitivity of the final state of the system to the initial parameters, we start with the Enzo3 simulation and increase by 5% either the initial velocity of the secondary (Enzo11) or the initial separation between the two particles (Enzo12), which correspond to initial eccentricities of 0.10 (Enzo11) and 0.05 (Enzo12). The evolution of the separation for those three simulations is compared in Figure 14. For Enzo11 and Enzo12, the ratio of the initial velocity of the companion to the velocity required for a circular orbit is higher than one (v₀/v_circ > 1), so the separation must first increase. The larger the orbital separation, the more delayed the rapid infall phase is, and the later the system reaches its final separation. The final separations for Enzo3, Enzo11, and Enzo12 are 11.7, 12.0, and 12.2 R_☉, respectively, and the final eccentricities are 0.09, 0.17, and 0.18, respectively. As expected, the companion that moves outward the farthest initially, sinks into the envelope with a higher orbital decay velocity. Therefore, it attains a more eccentric orbit and completes fewer revolutions around the primary core (Figure 14). However, the standard deviation of the final separation between the three simulations (σ ∼ 0.2 R_☉) is more than 10 times smaller than the width of a cell. Consequently, we conclude that the final results are quite insensitive to the initial conditions at the level tested.

Zoom In Zoom Out Reset image size

The drag exerted on the companion has two components: gravitational and hydrodynamical. The former is due to gravitational forces from matter flowing past the companion and colliding with its wake (Bondi & Hoyle 1944; Iben & Livio 1993), while the latter is due to ram pressure forces on the companion. The hydrodynamical contribution can be estimated as

$$\begin{equation} F_{\rm hydro} \sim \rho \mathbf {v}_2^{2} \times \pi R_2^2, \end{equation} \tag{ 21 }$$

where R₂ is the radius of the secondary, $\mathbf {v}_2$ is the relative velocity between the secondary and the envelope, and we have taken the coefficient of drag to be unity for simplicity. In a similar manner, the gravitational drag is approximated by (Iben & Livio 1993)

$$\begin{equation} F_{\rm grav} \sim \rho \mathbf {v}_2^{2} \times \pi R_A^2, \end{equation} \tag{ 22 }$$

where the accretion radius R_A is defined as

$$\begin{equation} R_A = \frac{2GM_2}{\mathbf {v}_2^2 + c_s^2}, \end{equation} \tag{ 23 }$$

where c_s is the sound speed of the medium. Choosing $|\mathbf {v}| = 2 c_s = 80$ km s⁻¹ with an 0.6 M_☉ companion yields R_A ∼ 30 R_☉. Assuming R₂ ∼ 1 R_☉, we conclude that the hydrodynamical drag is of the order of almost 1000 times smaller than the gravitational drag, and thus negligible.

This conclusion is also confirmed by the outcomes of our simulations. Indeed, the primary's core and the companion are treated as point masses and are not pressure sources, except for the primary's core in the SNSPH simulations. Instead of being caused by the finite size of the particles, hydrodynamical drag in the models is thus due to the matter accreted around them. We pointed out earlier that the accuracy with which accretion was treated was different between the two different models because of the different finest resolutions and softenings used: accretion is poorly modeled in the Enzo simulations whereas in the SNSPH simulations, the companion builds up a sphere of accreted matter about a few R_☉ wide around itself (Figure 13). This should lead to differences in the magnitude of hydrodynamic drag forces. Nevertheless, the consistency of the results suggests that the hydrodynamic drag is unimportant in the evolution of the system, confirming the results of Ricker & Taam (2008).

We now compare the numerical results with a sample of 61 observed post-CE systems listed in Zorotovic et al. (2010) and De Marco et al. (2011).

5.1.1. Final Separations

For a given companion mass (or alternatively mass ratio q) we obtain three values for the final separation A_f, one for each simulation carried out with that companion mass (Table 1 and Figure 15). One can distinguish between these values at high q (q ⩾ 0.34), which correspond to "heavy" companions (M₂ ⩾ 0.3 M_☉), and the ones at low q (q < 0.34) corresponding to "light" companions (M₂ < 0.3 M_☉). At high q, the values of A_f are very similar and the standard deviation is more than 20 times smaller than the average value of A_f. At low q, the companion sinks deeper and as a consequence, the resolution used in the 128³ Enzo simulations is not sufficient. However, as one increases the resolution to 256³ cells, the final separations converge to the solutions given by the SNSPH simulations.

Zoom In Zoom Out Reset image size

Figure 16 shows the distribution of orbital separations reached by the 61 post-CE systems. For all these systems, there has been no substantial orbital shrinkage due to phenomena such as magnetic braking or radiation of gravitational waves (see discussion in Schreiber & Gänsicke 2003). Although they cover a significant range in secondary masses, going from a 1.1 M_☉ MS star down to a 0.05 M_☉ brown dwarf, all of them have separations smaller than 11 R_☉. Furthermore, 87% of those systems have separations smaller than 4 R_☉, which is smaller than any value obtained in our simulations. This is even more obviously shown in Figure 17, where the final separations for simulations presented here and in the literature are compared to the orbital separations of the observed post-CE systems. Although a couple of observed systems have q ⩾ 0.5, one clearly sees that the simulations with M₂ = 0.9 and 0.6 M_☉ leave the companion far out. Systems with lower mass companions (M₂ ⩽ 0.3 M_☉) have by and large lower orbital separations than in our simulations. Simulations of Sandquist et al. (1998) and Ricker & Taam (2008) shown in Figure 17 give results consistent with ours. All these numerical simulations suggest that the separations between the secondary and the primary's remnant at the end of the simulated rapid infall phase are too large to explain the orbital separation of the currently observed post-CE systems. This suggests that further evolution of the orbital separation must occur during the phase immediately following the rapid infall phase. We discuss this point further in Section 5.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1.2. The State of the Envelope at the End of the Simulations

As shown in Table 2, most of the primary's envelope remains bound in all of our simulations. We study the situation in detail for our canonical model with the 0.6 M_☉ companion here. The evolution of the mass for different components is plotted in Figure 18. It first confirms that some envelope mass is unbound only during the first 50 days, after which neither angular momentum (Figure 8) nor kinetic energy (Figure 9) are exchanged between the unbound mass and the rest of the system. It also shows that more than 85% of the mass remains bound at the end of the simulation. This outcome, already pointed out by Sandquist et al. (1998), is quite intriguing, since the post-CE binaries observed must have succeeded in ejecting their envelope. After about 400 days, most of the envelope mass in our models has been moved to a larger radius (∼100 R_☉; see bottom panel in Figure 19), well outside the orbit of the primary core and the companion but remaining bound.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Amount of the Envelope Mass Still Bound at the End of the SNSPH Simulations

Name	M2	Mbounda
	(M☉)	(M☉)
SPH1	0.9	0.44
SPH2	0.6	0.44
SPH3	0.3	0.45
SPH4	0.15	0.46
SPH5	0.1	0.48

Note. ^aAt the start of the simulations, M_bound equals the total envelope mass M_e ≡ M₁ − M_c = 0.49 M_☉.

Download table as:  ASCIITypeset image

We now investigate how bound the final system is. We consider the center of mass of the system composed by the secondary and the mass within the current orbit as the center of our frame of reference. Then, we partition the domain into concentric shells with identical thickness, calculate the average radial velocity of each shell, and compare it to the escape velocity at that location. Figure 19 shows the escape velocity and the average radial velocity of the shells. The radial velocity is always positive and is similar to the space velocity at radii larger than 600 R_☉, as expected for envelope ejection. At radii smaller than 300 R_☉, the radial velocity is much smaller than the space velocity, suggesting that orbital motions dominate at those radii. All the mass within 10³ R_☉ is bound, which corresponds to more than 85% of the envelope mass. The remaining mass is found at radii between 10³ and 6 × 10³ R_☉, where the radial velocities are typically between 25 and 75 km s⁻¹. Those particles were initially in the outer parts of the giant star, and were the first to encounter the secondary. At that time of the inspiral, the shock was slightly supersonic (V₂ ∼ 35 km s⁻¹ and c_s ∼ 20 km s⁻¹). This regime of evolution is thus different from later phases when the secondary sinks deeper into the primary's envelope, where its velocity does not really increase (Figure 7) but the sound speed of the medium does (Figure 3).

We can measure how much extra energy would be required to unbind the envelope at each radius, using the definition

$$\begin{equation} E_{\rm extra} = \sum _i \frac{1}{2} M_i (v_{e,i} - v_{r,i})^2, \end{equation} \tag{ 24 }$$

where v_{e, i} and v_{r, i} are the escape velocity at the location of the ith shell and its average radial velocity, respectively. One finds E_extra ∼ 8.4 × 10⁴⁵ erg which represents just over 10% of the initial binding energy of the primary envelope. Thus, a relatively small additional input of energy could be sufficient to completely unbind the remaining envelope material.

We have compared here the final separations deduced from observations and those determined from the simulations. We have purposefully stayed away from calculating the ejection efficiency α (Webbink 1984; De Marco et al. 2011). Indeed, we question what the relevance of calculating α is when the envelope has not yet been fully ejected, true both in the Sandquist et al. (1998) and our simulations. We therefore defer for the moment the task of calculating α from simulation—a long-term goal of this project—until the simulations are more advanced.

In conclusion, the hydrodynamic simulations do not reproduce the post-CE systems in the sense that the system is left at too large separations and the envelope is not unbound at the end of the rapid infall phase. This means that either physical processes that are not accounted for in the simulations are responsible for the envelope ejection, or that the envelope ejection and a significant reduction of the orbit actually happens during the later subsequent slow inspiral phase. We discuss both possibilities in the following section.

In this section we first study and quantify physical processes that are not taken into account in our hydrodynamic simulations and that might be responsible for ejecting the envelope. Then, we focus on the subsequent slow inspiral phase and investigate whether the envelope can be ejected and the separation significantly reduced during this subsequent phase.

5.2.1. Rotation of the Primary

5.2.2. Physics not Included in the Simulations

The hydrodynamics codes use an ideal gas equation of state (Section 2.1) which, by definition, does not include variable abundances and the different ionization layers of the envelope. Han et al. (1995) suggested that recombination might play a role in CE interactions. As the outer parts of the envelope expand and cool, ions recombine with electrons, releasing energy that could aid in unbinding the envelope. Although it is unclear how efficient this process is and how much of the initial recombination budget can be used, one can calculate an upper limit on how much energy can be injected into the envelope by recombination.

According to our stellar evolution model, the hydrogen fraction within the convective envelope of our RGB star is X ∼ 0.68. The mass of the envelope is M_e = 0.49 M_☉ and each proton recombining with an electron produces an energy E₀ = 13.6 eV. We also have to calculate how much of the envelope is ionized. Therefore, we calculate the partition functions Z for hydrogen. The hydrogen ion has no degeneracy so Z₂ = 1. The partition function for the hydrogen atom at temperature T is

$$\begin{equation} Z_1 = \sum _{n=1}^{\infty } 2n^2 \exp {\frac{E_0(1/n^2-1)}{kT}}, \end{equation} \tag{ 25 }$$

where k = 8.6173 × 10⁻⁵ eV K⁻¹ is the Boltzmann constant. We truncate the sum in Equation (25) at the first integer n_max such that the distance at which the electron orbits the proton for this quantum number is larger than l_max = 10⁻⁶ cm, i.e., a₀n²_max > l_max, where a₀ = 5.2918 × 10⁹ cm is the Bohr radius (Miranda 2001). We then use the Saha formula to calculate the ratio of ionized to neutral hydrogen (Carroll & Ostlie 1996):

$$\begin{equation} N_2/N_1 = \frac{2Z_2}{n_e Z_1}\left(\frac{2\pi m_e k T}{h^2}\right)^{3/2} \exp {(-E_0/kT)}, \end{equation} \tag{ 26 }$$

where n_e is the number density of free electrons and m_e is the electron mass. We find that 91% of the envelope is ionized. Consequently, the recombination of the whole ionized envelope would produce an extra energy

$$\begin{equation} E_{\rm recomb} = 0.91 \times XM_e \frac{N_A}{M_H} \times 13.6\ {\rm eV}, \end{equation} \tag{ 27 }$$

where N_A is the Avogadro number and M_H is the atomic mass of hydrogen. One finds E_recomb ∼ 1.18 × 10⁴⁶ erg, which is slightly higher than the extra energy E_extra required to eject the envelope in our canonical model (Section 5.1.2). Thus, we conclude that recombination in the envelope could substantially aid in unbinding it.

Another source of energy could be radiation pressure. For low- and intermediate-mass giants in hydrostatic equilibrium, radiation pressure (P_rad ≡ aT⁴/3, where a is the radiation constant) is negligible compared to gas pressure (Equation (4)): for our primary, P_rad/P_gas ≲ 0.01 except in a small zone (0.1 R_☉ ⩽r ⩽ 10 R_☉), where P_rad/P_gas ≲ 0.1. However, the deep inspiral of the companion within the primary's envelope will induce local shock heating. The increase of temperature is proportional to the square of the Mach number (Tarbell et al. 1999), so even if the companion is orbiting at twice the local sound speed, the radiation pressure to gas pressure after the shock becomes

$$\begin{equation} \left(\frac{P_{{\rm rad}}}{P_{{\rm gas}}} \right)^{\rm after} \propto \left(\frac{P_{{\rm rad}}}{P_{{\rm gas}}} \right)^{\rm before} (M^2)^3 = 6.4. \end{equation} \tag{ 28 }$$

Therefore, including radiation pressure in the equation of state will increase the total pressure locally and might reduce the energy required to eject the envelope. However, it is possible that this effect is globally small, since this extra heating source is probably very localized around the companion.

5.2.3. The Post-rapid-infall Phase

At the end of the rapid infall phase, the orbit is stable until the end of the simulations (a few more years). Consequently, there is no further hydrodynamical coupling between the extended envelope and the surviving binary. We now investigate whether the envelope is likely to be ejected during this slower inspiral phase.

Although the resolution of the simulations prevents us from quantifying how much envelope will be left around the core of the primary, one can still describe qualitatively what the evolution of the primary's remnant will be. Figure 18 shows that less than 10⁻² M_☉ is left around the primary's core, so the primary will depart the giant branch (Bloecker 1995; but see also the discussion in De Marco et al. 2011). Then, two scenarios might occur depending on how long the partially ejected envelope will take to fall back.

If the star is given enough time to transit to the blue due to hydrogen burning at the base of the envelope before the lifted envelope falls back, the star will readjust on its thermal timescale of the remaining envelope, and eventually end its life as a helium white dwarf. This transition will last ∼10³ years during which the star will have a luminosity between 300 and 1000 L_☉ (Iben & Tutukov 1993, see their Figure 1), which is consistent with the more recent work of Driebe et al. (1998, see their Figure 1). If we assume the remnant to have a luminosity L_c ∼ 500 L_☉, we can compare the gravitational acceleration of a gas particle with the radiation acceleration defined by

$$\begin{equation} a_{{\rm rad}} = \frac{L_c}{4\pi r^2} \frac{\kappa }{c}, \end{equation} \tag{ 29 }$$

where r is the distance between the gas particle and the core, κ = 0.4 cm² g⁻¹ is the opacity for Thomson scattering for hydrogen, and c is the speed of light. We still find the radiation acceleration to be overall almost two orders of magnitude smaller than the gravitational acceleration.

If on the contrary, the envelope falls back before the primary's remnant had crossed the Hertzsprung–Russell diagram, a circumbinary disk will form (Kashi & Soker 2011). They refer to the numerical work done by Artymowicz et al. (1991), which suggests that in such a configuration, the binary separation will decrease due to Lindblad resonances—mainly—as well as viscous tides. Although this mechanism has the advantage of explaining how the orbital separation will diminish during the subsequent phase, the ability of radiation to eject the gas will even be reduced in comparison with the previous situation, so it is not clear how the latter will eventually be unbound.

In conclusion, radiation acceleration alone does not seem to be responsible for unbinding the remaining gas, regardless of the time the partially ejected envelope will remain suspended.