{'fraction_non_alphanumeric': 0.030185688813644796, 'fraction_numerical': 0.020171820686922536, 'mean_word_length': 4.43553597650514, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating the Common Envelope Phase of a Red Giant Using SPH and Uniform Grid Codes

We use three-dimensional hydrodynamical simulations to study the rapid infall phase of the common envelope interaction of a red giant branch star of mass equal to 0.88 \msun and a companion star of mass ranging from 0.9 down to 0.1 \msun. We first compare the results obtained using two different numerical techniques with different resolutions, and find overall very good agreement. We then compare the outcomes of those simulations with observed systems thought to have gone through a common envelope. The simulations fail to reproduce those systems in the sense that most of the envelope of the donor remains bound at the end of the simulations and the final orbital separations between the donor's remnant and the companion, ranging from 26.8 down to 5.9 \rsun, are larger than the ones observed. We suggest that this discrepancy vouches for recombination playing an essential role in the ejection of the envelope and/or significant shrinkage of the orbit happening in the subsequent phase.

During the giant phases of the primary, companions closer than ∼ 5 AU enter a strong interaction phase with the primary and, under certain circumstances, a common envelope (CE) may form around the two stars. The secondary star spirals inside the envelope of the primary and may also fill its own Roche lobe because it cannot accrete all the matter coming from the donor star. This process is called a common envelope interaction and was originally described by Paczynski (1976). For a general review of the topic see, e.g., Iben & Livio (1993). There are two different processes leading to the onset of a CE phase: the start of unstable mass transfer from the expanding primary to the secondary (Hjellming & Webbink 1987;Hurley et al. 2002) and the development of a tidal instability that occurs if there is not enough angular momentum in the orbit to maintain the primary's envelope in synchronization (Darwin 1879). The post-CE system will be either a compact binary system, if there is enough energy to eject the primary's envelope, or a merger, if not.
The CE interaction is an essential ingredient for any binary population synthesis study of intermediate (e.g., Politano et al. 2010) or massive stars (e.g., Belczynski et al. 2008). Compact binaries are believed to be formed through at least one CE phase. Among them are symbiotic binaries, supersoft X-ray sources, cataclysmic variables and double white dwarfs, which are all possible supernova Type Ia progenitors. As Meng et al. (2010) pointed out, results deduced from population synthesis studies such as the Type Ia supernova birth rate are highly dependent on the physics of the CE phase. Therefore, it is paramount to understand more accurately the CE interaction in order to identify the formation channels of such supernovae and to compare observations with predictive models. Moreover, many substellar companions to evolved stars have recently been discovered with small orbital separation. Maxted et al. (2006) found a brown dwarf orbiting a white dwarf with a 116 min period, while Setiawan et al. (2010) discovered a system composed of a Jupiter-like object orbiting an horizontal branch star with a 16.2 days period. We therefore know that substellar companions can survive a CE interaction, but what is the minimum mass of the companion that can eject the envelope? Is the ejected envelope entirely unbound or will some of it eventually fall back and form a circumbinary disk? Were the substellar companions present before and survived the CE or were they formed later on in such a disk (Perets 2010)? Those questions remain unanswered.
Although the CE process was outlined more than 30 years ago, it is still far from understood quantitatively. Numerical simulations suggest that the typical duration of the entire CE phase is short -less than 10 3 years -which makes CE ejections unlikely to be observed. However, one can use observations of post-CE binaries to better understand CE evolution. With the use of stellar models, the initial configuration of such systems can be approximately determined from the final configuration. Using either the α-formalism (Webbink 1984, but see De Marco et al. 2011) or the γ-formalism (Nelemans et al. 2000) the relevant parameters can be constrained and the CE ejection efficiency can be predicted. Using this approach, De Marco et al. (2011) suggested an anti-correlation between α, the CE efficiency parameter, and the secondary to primary mass ratio.
The entire CE evolution can be divided into three different phases (Podsiadlowski 2001) with different timescales, length scales and physics involved. These differences are the reasons why reproducing the entire CE evolution of a given system accurately is challenging. Therefore, one usually treats one phase after the other with different methods. In this paper, we focus only on the rapid infall phase, which has a short timescale (∼ 1 − 10 years), and in which the evolution is driven by drag forces. Several numerical hydrodynamic studies of the CE interaction have been carried out in the past (for an exhaustive list, see , including a series of ten papers starting with the two-dimensional calculation of the interaction of a 16 M supergiant and a 1 M neutron star (Bodenheimer & Taam 1984), and most recently treating three-dimensional simulations of the CE interaction between 3 or 5 M giant stars and 0.4 or 0.6 M main sequence (MS) companions (Sandquist et al. 1998). The latter study has been extended first by Sandquist et al. (2000) to 1 M and 2 M red giant branch (RGB) stars with companion masses ranging from 0.1 to 0.45 M , then by De Marco et al. (2003) to a 1 M asymptotic giant branch star with a 0.1 or 0.2 M companion. Ricker & Taam (2008) computed high resolution simulations of the CE phase between a 1.05 M RGB star and a 0.6 M compact companion, and concluded that the gravitational component of the drag dominates over the hydrodynamical component (also see Taam & Ricker 2010;Ricker & Taam 2011).
A direct comparison of the results obtained using different numerical methods has however never been carried out. Although analytical/empirical work has included discussion regarding observational data, there are only a couple of publications that connect simulations and observations in a meaningful way (see e.g., Sandquist et al. 2000). Those are, as we will explain in §2 and §4, key steps to better understand the implications of CE interactions and the physical processes driving them. In this paper we therefore present numerical simulations with two different algorithms of the CE interaction of a 0.88 M RGB star with a MS companion. Different companion masses from 0.1 M to 0.9 M are considered. The simulations are carried out with both an Eulerian code (Enzo in uniform-grid mode, O'Shea et al. (2004) and enzo.googlecode.com) and a Lagrangian code (SNSPH, Fryer et al. 2006), and for different resolutions. We describe the numerical methods and the initial conditions of our 15 simulations in §2 and §3. We describe and discuss the results in §4 and §5, and finally conclude and summarize in §6.

Description of codes
In this section we describe the numerical methods we use. We first compare the code algorithms and explain why a code-to-code comparison is necessary. Then, we describe both codes in detail and finally discuss different ways to compare resolution.

Eulerian vs Lagrangian 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 under unity. On the other hand, SPH codes don't 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.

Input physics
Both codes solve the fully compressible hydrodynamics equations with self-gravity included. These equations can be written using an Eulerian formulation: where ρ, v, p, Φ, u, u int , γ 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: 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 (Eq. 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 §5.2.2. Finally, the gravitational potential is calculated using the Poisson equation (Eq. 5).

The Enzo code
Enzo is a three-dimensional, adaptive mesh refinement hybrid (hydrodynamics + Nbody) 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 (Eqs. 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 better resolves shocks and turbulence, our tests show that it slows down the computation by a factor 2. As we will point out in §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 a RGB star and a 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 Fig. 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 that each particle represents many actual particles, such as in cosmological simulations with dark matter. In that case their mass is deposited in the 8 nearest cells and added to the gas density of those cells to find the total density for use in solving the Poisson equation (Eq. 5). In a simple two-body interaction between 1 M and 0.1 M objects in a one 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 6 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): where M P M 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: where C 1 = 0.4 is the Courant factor, C 2 = 0.4 is the particle Courant factor, c s is the sound is the velocity of a particle and 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.

The SNSPH code
SNSPH (Fryer, Rockefeller, & Warren 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 a SPH particle, i, is smoothed using the following formula: 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 Fig. 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.
SNSPH uses the fast multipole method to calculate gravitational accelerations (Warren & Salmon 1993). The SPH particles are also advanced using an leapfrog algorithm. Finally, in order to keep the same overall spatial coverage, the smoothing length varies according to the formula from Benz (1989): Plotted are the theoretical potential (solid line), the Ruffert (1993) potential used in Enzo (dashed line) and the Monaghan (1992) one used in SNSPH (dash-cross line).

Resolution comparison
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: where N tot , V 1 , V G , R 1 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 Fig. 2, the smoothing length at small radii does not vary, so an Enzo run with a 128 3 grid will be under-resolved compared to our canonical 500 000 (roughly 80 3 ) particle SPH run no matter how deep the companion penetrates while a run with a 256 3 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. -Resolution comparison between the SNSPH smoothing length field (dots) for a run with 500 000 particles, and the Enzo size of a grid cell for the 128 3 (dash line) and the 256 3 (solid line) runs for the initial (left) and final (right) particle distributions.
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 under-resolve 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.

The simulations
We perform 5 SNSPH and 12 Enzo simulations of CE interactions with a 0.88 M RGB primary that are summarized in Table 1. The SNSPH simulations are computed using 500 000 particles whose initial smoothing length follows the radial profile shown in Fig. 2. The Enzo simulations are performed using either a 128 3 or a 256 3 grid. In both cases, the linear size of the computational domain is L = 3 × 10 13 cm. We consider companion masses of 0.9, 0.6, 0.3, 0.15 and 0.1 M . Giant stars are slow rotators with rotational velocities of the order of a few km s −1 (de Medeiros & Mayor 1999). Although it is expected that a close companion will, through the action of tides and the transport of angular momentum in the primary envelope, spin up the envelope during the pre-CE phase, the actual rotation of the primary at the onset of the CE interaction is hard to quantify. Moreover, even if the primary was uniformly rotating at 50 km s −1 , its rotational energy would be where ω, r g , M 1 and R 1 are the angular velocity, the radius of gyration, the mass and the radius of the primary, respectively. For RGB stars r g is typically about 0.1 . This rotational energy does not affect the energetics of the system since it is more than two orders of magnitude smaller than the binding energy of the primary (see below). Consequently, we assume that the primary is initially non-rotating. Finally, the companion is at the start placed at the surface of the primary in a circular orbit. We thus have three different simulations for each initial companion mass -one with SNSPH, and two with Enzo on 128 3 and 256 3 grids. Additionally, we also run two Enzo simulations in order to study the dependency of the final parameters on the initial conditions. We consider the 128 3 Enzo simulation with a 0.3 M companion (Enzo3) as the reference and run identical simulations increasing, by 5 %, either the initial velocity of the companion (Enzo11) or the initial separation (Enzo12). All the runs follow the evolution of the system for about 1 000 days.  ), the companion mass, the initial orbital separation (A 0 ), the initial orbital period (P 0 ), the ratio of the initial orbital velocity of the companion (v 0 ) to the velocity required for a circular orbit (v circ ) and the final orbital separation (A f ) taken at the end of the rapid infall phase ( § 4.1).
As a primary, we use a one-dimensional model of a star with a MS mass of 1 M . Using the stellar evolution code EVOL (Herwig 2000), this progenitor was evolved to the RGB phase until the core reached M c = 0.392 M . At that time, the radius of the star was 83 R and its total mass was M 1 = 0.88 M due to mass loss, which was treated using the Reimers formalism with η = 0.5. We adapt this model by using the density and pressure profiles, but computing the internal energy using Eq. 4. A sample of relevant profiles are plotted in Fig. 3.
We now explain how this stellar model is modified in order to be compatible with an input suitable for each of our codes. For the SNSPH simulations, the initial particle configuration is a weighted Voronoi tesselation (WVT) similar to that described by Diehl & Statler (2006). As we have explained in §2.1 one limitation of SPH codes is the large number of particles required by dense regions such as the core of the primary. Since the time step  (blue), the SPH initial model computed with 500 000 particles (black) and the Enzo initial models for a 128 3 (green) and a 256 3 (purple) unigrid. The vertical line represents the core-envelope boundary according to the criterion of De Marco et al. (2011).
induced by a particle i can be roughly estimated by h i /c s,i where c s,i is the local sound speed, a small smoothing length will require a small time step resulting in a high computational cost. Since the equation of state changes significantly around the helium core, we represent the core by a particle with mass M c . The associated smoothing length is h c = 0.1 R . We add SPH particles in the region around the core such that the density values and gradient profiles connect smoothly at the core/envelope boundary (r = 2h c = 0.2 R ). In this way, we obtain the profile shown in Fig. 3. Since the density profile has been changed, one must modify the gravitational acceleration accordingly. Assuming hydrostatic equilibrium in spherical symmetry, we integrate the pressure gradient choosing the integration constant to match the true profile outside the core (at r = 0.2 R ). The specific energy profile is computed using Eq. 4. Finally, the acceleration of a SPH particle is due either to gravity or to gas pressure. These two components are computed using the same particle mass for all particles except the core and the companion, for which we distinguish between the gravitational and SPH masses. The gravitational mass of the core is M c and its SPH mass is set to balance the gravitational acceleration of the envelope and prevent the star from collapsing. As for the companion, we treat it as an N-body particle so its SPH mass is 0 M .
For the Enzo simulations, the grid is initialized using the stellar model of the primary with the addition of a PM particle that represents its core. We fill the computational domain with a constant background density to prevent the star from expanding and set the ratio between the background density and the minimum density of a cell that belongs to the primary to 10 −4 . This setup is not initially numerically stable. The star tends to expand, so we let the initial configuration evolve for a few dynamical times in the absence of the companion, while damping the velocity field by a factor of 2 after each cycle. Finally, we evolve this relaxed model normally for another few dynamical times to obtain a numerically stable model. As a side effect of the relaxation to hydrostatic equilibrium, the Enzo models are a little bit bigger -the lower the resolution, the larger the radius of the primary isthus the initial orbital separations between the models are slightly different (Table 1).

Results
In this section, we describe the results obtained from our 15 simulations. Since the qualitative behavior is the same in all of them, we detail the 0.6 M case (SPH2, Enzo2 and Enzo7).

Description of the rapid infall phase
As explained in §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 3 Enzo simulations is shown in Fig. 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 (Fig. 5). As expected, the orbital decay is initially quite high (∼ 0.01 day −1 ), 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 (Fig. 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. 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 (Fig. 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 (Fig. 7) with the local sound speed of the gas (Fig. 3, bottom left panel). The former does not exceed 50 km s −1 while the highest sound speed encountered is about 60 km s −1 . 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. 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 Fig. 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, no more additional mass gets unbound (see § 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 more exchange of angular momentum between the unbound mass and the rest of the system. After ∼ 150 days, there is no more angular momentum exchange in the system. The primary's core and the companion -which are the main contributors to the calculation of the centre 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.
We plot the various energy components in Fig. 9. We start by explaining the different components of potential, thermal and kinetic energies represent 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, 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 §3): 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 Finally, we subtract the contribution of the secondary from the total potential energy in order to calculate the binding energy of the envelope: where the subscript "2" stands for the secondary.
During the first 200 days when most of the in-spiral 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 Fig. 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. Fig. 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.

Code comparison
The fact that a code solves the equations in an accurate and precise way in a particular situation does not necessary 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 Fig. 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 (Fig. 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 −2 M between the Enzo and the SNSPH runs. This difference is large compared with the mass of a SPH particle (∼ 10 −6 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 Fig. 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 −3 g cm −3 . 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 3 runs is about 1.6 R so the accretion zone cannot be resolved although it is still better than for the Enzo 128 3 simulations as can be seen from comparing the different density profiles at 50 days (Fig. 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 3 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. 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 Fig. 7 that for the 0.6 M companion, the velocity of the companion stays below 50 km s −1 and therefore, the gas flows are subsonic except in the outer layers. This conclusion was also reached by Ricker & Taam (2008).

The impact of initial conditions
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 Fig. 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 0 /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 outwards 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 (Fig. 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.

Gravitational vs Hydrodynamic drag
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: where R 2 is the radius of the secondary, 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): where the accretion radius R A is defined as: where c s is the sound speed of the medium. Choosing |v| = 2c s = 80 km s −1 with an 0.6 M companion yields R A ∼ 30 R . Assuming R 2 ∼ 1 R , we conclude that the hydrodynamical drag is of the order of almost 1 000 times smaller than the gravitational drag, 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 (Fig. 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).

Comparison of simulations and observations
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).

Final separations
For a given companion mass (or alternatively mass ratio q) we obtain 3 values for the final separation A f , one for each simulation carried out with that companion mass (Table 1 and Fig. 15). One can distinguish between these values at high q (q ≥ 0.34), which correspond to "heavy" companions (M 2 ≥ 0.3 M ), and the ones at low q (q < 0.34) corresponding to "light" companions (M 2 < 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 3 Enzo simulations is not sufficient. However, as one increases the resolution to 256 3 cells, the final separations converge to the solutions given by the SNSPH simulations. 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 Fig. 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 2 = 0.9 and 0.6 M leave the companion far out. Systems with lower mass companions (M 2 ≤ 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 Fig. 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 §5.2.

The state of the envelope at the end of the simulations
As shown on 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 Fig. 18. It first confirms that some envelope mass is unbound only during the first 50 days, after which neither angular momentum (Fig. 8) nor kinetic energy (Fig. 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 Fig. 19), well outside the orbit of the primary core and the companion but remaining bound. 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. Fig. 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 3 R is bound, which corresponds to more than 85% of the envelope mass. The remaining mass is found at radii between 10 3 and 6 × 10 3 R , where the radial velocities are typically between 25 and 75 km s −1 . 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 in-spiral, the shock was slightly supersonic (V 2 ∼ 35 km s −1 and c s ∼ 20 km s −1 ). 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 (Fig. 7) but the sound speed of the medium does (Fig. 3).
We can measure how much extra energy would be required to unbind the envelope at each radius, using the definition where v e,i and v r,i are the escape velocity at the location of the i-th shell and its average radial velocity, respectively. One finds E extra ∼ 8.4 × 10 45 ergs 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 the envelope ejection and a significant reduction of the orbit actually happens during the later subsequent slow in-spiral phase. We discuss both possibilities in the following section.

Reproducing the observations
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 in-spiral phase and investigate whether the envelope can be ejected and the separation significantly reduced during this subsequent phase.

Rotation of the primary
The envelope of the progenitor is initially non-rotating and although the calculation done in §3 shows that, regardless of the initial rotation velocity of the envelope, its rotational energy is negligible in comparison with its binding energy, we suspected at first that the absence of rotation might be the reason for most of the envelope to remain bound. However, Sandquist et al. (1998) carried out two identical simulations where they modified the initial rotation state of the primary from a giant star in synchronization with the orbit to a nonrotating one (their simulations 1 and 2). In both cases, the evolution of the bound mass and the final orbital parameters are similar. It thus does not seem that changing the initial rotation of the primary leads to a different CE outcome.

Physics not included in the simulations
The hydrodynamics codes use an ideal gas equation of state ( §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 0 = 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 2 = 1. The partition function for the hydrogen atom at temperature T is where k = 8.6173 × 10 −5 eV.K −1 is the Boltzmann constant. We truncate the sum in Eq. 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 −6 cm, i.e. a 0 n 2 max > l max where a 0 = 5.2918×10 9 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): 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 where N A is the Avogadro number and M H the atomic mass of hydrogen. One finds E recomb ∼ 1.18 × 10 46 ergs, which is slightly higher than the extra energy E extra required to eject the envelope in our canonical model ( §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 4 /3, where a is the radiation constant) is negligible compared to gas pressure (Eq. 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 in-spiral 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: P rad P gas after ∝ P rad P gas before (M 2 ) 3 = 6.4 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.

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 in-spiral 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. Fig. 18 shows that less than 10 −2 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 3 years during which the star will have a luminosity between 300 and 1 000 L (Iben & Tutukov 1993, their Fig. 1), which is consistent with the more recent work of Driebe et al. (1998) (their Fig. 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 where r is the distance between the gas particle and the core, κ = 0.4 cm 2 g −1 is the opacity for Thompson 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 -mainlyas 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 for.

Summary and Future Work
In this work we have carried out three-dimensional hydrodynamic simulations of the CE interaction between a 0.88 M RGB star and companions with mass ranging from 0.1 to 0.9 M . We have used both an Eulerian grid code (Enzo) and a Lagrangian SPH code (SNSPH) with various resolutions. They both have advantages and disadvantages and can be used for different purposes: while one might rather use SPH to study the accretion around the secondary, even a uniform grid code is more suitable in resolving the low-density extended envelope. Of course, adaptive mesh refinement combines the advantages of both of these methods at the cost of increased code complexity.
We first compared the outcomes of those simulations with each other. We found that the results are very similar for companion masses M 2 0.3 M . We thus conclude that in this regime, the resolutions used are sufficient to study the global evolution of the system during the rapid infall phase of the interaction, which is driven mainly by gravitational drag. For lower companion masses (M 2 0.3 M ) that penetrate deeper in the giant's envelope, the 128 3 Enzo runs are under-resolved but the Enzo results converge to the solutions from the SNSPH simulations.
We then compared the outcomes of our simulations with observed post-CE systems. The final separations are found to be systematically higher than those deduced from observations, as is the case for the past simulations by Sandquist et al. (1998), De Marco et al. (2003 and Ricker & Taam (2008). Moreover, mass is only unbound during the early stages of the interaction (∼ 50 days for the 0.6 M companion) and most of the envelope remains bound at the end of the simulations, as was the case for the earlier simulations of Sandquist et al. (1998). We investigated whether there might be additional processes that were not accounted for in the simulations. We found that recombination can contribute significantly, but stellar rotation and radiation pressure play only marginal roles. Finally, we wondered whether the bound envelope is a result of imprecise simulations or a real physical feature. If the latter, then one would have to follow the subsequent evolution of the system to determine the actual outcome of the CE. Fall back disks may form and even have an impact on the inner binary (Artymowicz et al. 1991, Kashi & Soker 2011. After the submission of this paper, Ricker and Taam made their paper Ricker & Taam (2011) available. This paper continues the work introduced in Ricker & Taam (2008). In their simulation, only about 25 % of the primary's envelope is unbound. Although this value is slightly higher than ours, it is in agreement with our work in the sense that most of the envelope remains bound. They also claim that the ejection occurs mostly in the orbital plane, as it is the case in our simulations. However, the extended envelope at the end of their simulation is rotating much faster than it is expanding which is in contradiction with our results ( § 5.1.2) but might be due to the fact that their primary is initially rotating.

Acknowledgments
J-CP, ODM and M-MML acknowledge funding from NSF grant 0607111. FH acknowledges funding from an NSERC Discovery grant. J-CP thanks Colin McNally for useful comments and discussions. J-CP is grateful to both the Enzo and yt (Turk et al. 2011) communities for their precious help, especially Matthew Turk's dedication. The authors acknowledge computer time provided by Westgrid and Compute Canada as well as the Los Alamos National Laboratory. The authors thank an anonymous referee for his useful comments that helped to improve this manuscript.  Fig. 9.-Energy components for the SPH2 simulations. Plotted are the total energy (E tot ), the total gravitational potential energy Φ tot , the internal energy of the system (U tot ), the gravitational potential energy of the envelope (Φ env ), the gravitational potential energy from the core-companion interaction (Φ c2 ), the kinetic energy of the core and the companion (K c2 ), the kinetic energy of the bound mass (K b ), the kinetic energy of the unbound mass (K u ), the orbital energy of the core-companion system (E c2 ) and the total energy of the envelope (E env ≡ Φ env + U tot + K b ). The beat frequency seen on K c2 and Φ c2 are due to the non-synchronization between the orbital period and the data dumping frequency.   Fig. 11.-Separation between the core of the primary and the 0.6 M companion as a function of time for the SPH2 (left), Enzo2 (middle) and Enzo7 (right) simulations. Again, the beat frequency seen in the SPH simulation is due to the non-synchronization between the orbital period and the dumping frequency.