Common Envelope Evolution on a Moving Mesh

TL;DR

This paper presents a new simulation approach for common envelope evolution using a moving-mesh code, demonstrating detailed results on envelope ejection and orbital changes in binary star systems.

Contribution

It introduces an extended moving-mesh code with multiple time-steps and realistic physics for simulating CEE, providing new insights into envelope ejection and orbital dynamics.

Findings

Significant envelope mass ejection (~63-66%) including recombination energy.

Corotation increases mass ejection and final orbital separation.

Potential for complete envelope ejection during the plunge phase.

Abstract

We outline the methodology of simulating common envelope evolution (CEE) with the moving-mesh code MANGA. We extend MANGA to include multiple time-steps. This provides substantial speedups for problems with large dynamic range. We describe the implementation of realistic equations of state relevant in stellar structure and the generation of suitable initial conditions. We then carry out two example simulations of a 2 $\rm{M}_{\odot}$ red giant with a 0.36 $\rm{M}_{\odot}$ core and a 1 $\rm{M}_{\odot}$ companion undergoing CEE for 240 days. In one simulation the red giant is set into corotation with the orbital motion and in the other it is non-rotating. We find that the separation between the companion and red giant core shrinks from 52 $\rm{R}_{\odot}$ to 3.6 $\rm{R}_{\odot}$ and 3.2 $\rm{R}_{\odot}$ respectively, ending with an eccentricity of 0.1. We also find that 66 and 63 per cent of the envelope mass is ejected. This is higher than in many previous works. Several reasons for this are discussed. These include our inclusion of recombination energy. Our simulations show that putting giants in corotation increases the fraction of mass ejected from the system and results in a larger final orbital separation. We conclude that the entire envelope of the red giant might be ejected during the plunge phase of CEE in this region of parameter space.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathbf{v}$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \rho \mathbf{v}\mathbf{v}+\nabla P$$

$$\frac{\partial \rho e}{\partial t}+\nabla \cdot \left(\rho e+P\right)\mathbf{v}$$

$$\frac{\partial \mathcal{U}}{\partial t}+{\int }_V\nabla \cdot \mathcal{F}dV=\mathcal{S},$$

$${\mathbf{U}}_i={\int }_{V_i}\mathcal{U}dV={\mathcal{U}}_i{V}_i,$$

$${\int }_{\partial V_i}\nabla \cdot \mathcal{F}dV={\int }_{\partial V_i}\mathcal{F}\cdot \hat{\mathbf{n}}dA,$$

$$\sum _{j\in \mathrm{neighbors}}{\mathbf{F}}_{ij}{A}_{ij}={\int }_i\mathcal{F}\cdot \hat{\mathbf{n}}dA,$$

$${\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n+\Delta t\sum _j{\hat{\mathbf{F}}}_{ij}{A}_{ij}+\Delta t{\mathbf{S}}_i,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho {\mathbf{v}}^{\prime }$$

$$\frac{\partial \rho {\mathbf{v}}^{\prime }}{\partial t}+\nabla \cdot \rho {\mathbf{v}}^{\prime }{\mathbf{v}}^{\prime }+\nabla P$$

$$\frac{\partial \rho e^{\prime }}{\partial t}+\nabla \cdot \left(\rho e^{\prime }+P\right){\mathbf{v}}^{\prime }$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathbf{v}-\mathbf{w}\cdot \nabla \rho$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \rho \mathbf{v}\mathbf{v}-\mathbf{w}\cdot \nabla \rho \mathbf{v}+\nabla P$$

$$\frac{\partial \rho e}{\partial t}+\nabla \cdot \rho e\mathbf{v}-\mathbf{w}\cdot \nabla \rho e+\nabla \cdot P\mathbf{v}$$

$$v_{\varphi }=f_{\mathrm{cr}}{\Omega }_{\mathrm{orb}}{R}_{\mathrm{cyl}},$$

$$v_{\varphi }=\frac{f_{\mathrm{cr}}{\Omega }_{\mathrm{orb}}{R}_{\ast }^2{\sin }^2(\theta )}{R_{\mathrm{cyl}}}.$$

$$E_{\mathrm{tot},\mathrm{i}}=m_i[\frac{1}{2}({\mathbf{v}}_i-{\mathbf{v}}_{\mathrm{CM}}{)}^2+{\varphi }_i+I_{e,i}]$$

$${\mathbf{v}}_{\mathrm{CM}}=\frac{m_{\mathrm{c},1}{\mathbf{v}}_1+m_{\mathrm{c},2}{\mathbf{v}}_2+m_{\mathrm{g},\mathrm{b}}{\mathbf{v}}_{\mathrm{g},\mathrm{b}}}{m_{\mathrm{c},1}+m_{\mathrm{c},2}+m_{\mathrm{g},\mathrm{b}}},$$

$$f_{\mathrm{unb}}=\frac{m_{\mathrm{g},\mathrm{u}}}{m_{\mathrm{env}}},$$

$$e=\frac{r_{\mathrm{apo}}-r_{\mathrm{peri}}}{r_{\mathrm{apo}}+r_{\mathrm{peri}}},$$

Full text

††pubyear: 2019††pagerange: Common Envelope Evolution on a Moving Mesh–Common Envelope Evolution on a Moving Mesh

1 Introduction

Common envelope evolution (CEE) is the name given to a brief but important phase in the evolution of a binary system where two stars – one compact, one giant – share a common envelope (for a review see Ivanova et al., 2013). As the core of the giant and the compact star orbit each other, loss of orbital energy and angular momentum drive the two stellar centres close to one another and eject the envelope from the system. CEE is a critical process in the lives of these binary systems and is responsible for progenitors of Type Ia supernova (potentially), X-ray binaries, double white dwarfs, double neutron stars and possibly the merging double black hole and double neutron star systems recently discovered by Advanced LIGO (Ivanova et al., 2013; Belczynski et al., 2016).

In spite of its importance, CEE is not well understood. Our knowledge of CEE mainly comes from population necessity rather than theory or direct observation. Part of this problem is the diverse physics, including hydrodynamics, convection, turbulence, accretion, nuclear burning, radiation, self-gravity, ionization/recombination, jets and magnetism involved. In addition, the physical time-scales span the dynamical time of compact remnants or cores (milliseconds to seconds) to the many thermal times of an envelope (years) (Webbink, 2008; Ivanova et al., 2013).

In the presence of the different physics and time-scales, astronomers have mainly employed conserved quantities – energy and angular momentum – to roughly model CEE. These are known as the energy formalism (Webbink, 1984; Livio & Soker, 1988) and the angular momentum (or $\gamma$) formalism (Nelemans et al., 2000; Nelemans & Tout, 2005). To make progress beyond these formalisms, fully 3-D numerical tools must be brought to bear. Here the defining numerical studies are those carried out by Passy et al. (2012), Ricker & Taam (2012), Nandez et al. (2015) and Ohlmann et al. (2016). These simulations generally find a wide variation in efficiencies from about 5 to 100 per cent. Nandez et al. (2015) use a smooth particle hydrodynamics (SPH) code, Passy et al. (2012) used both a SPH and an Eulerian code, Ricker & Taam (2012) used an adaptive mesh refinement Eulerian code and Ohlmann et al. (2016) used perhaps the most interesting code of all – the moving-mesh (MM) code, AREPO.

AREPO is a class of arbitrary Lagrangian-Eulerian (ALE) or MM schemes that have been devised as an effort to capture the best characteristics of both Lagrangian and Eulerian approaches, combining superior conservation properties of Lagrangian schemes with the superior shock capturing of Eulerian schemes. Springel (2010, hereinafter S10) described a usable ALE scheme that has proven successful and is implemented in the code AREPO. The scheme relies on a Voronoi tessellation (Okabe et al., 1992) to generate a well-defined, unique, and continuously deformable mesh for an arbitrary distribution of points.

S10 argued that the use of ALE schemes is important to maintain the Galilean invariance of Eulerian schemes. It has also been argued that these schemes are superior at capturing boundary layer instabilities such as Kelvin-Helmholtz instabilities to SPH and Eulerian grid schemes (S10, but also see Lecoanet et al. 2016). In any case, MM methods do seem ideal to model colliding galaxies or stars. In particular, AREPO has been used in a number of different problems including cosmological galaxy formation (see for instance Vogelsberger et al., 2014), and stellar mergers (Zhu et al., 2015; Ohlmann et al., 2016).

Recently, we have developed a moving-mesh hydrodynamic solver for the N-body simulation code CHANGA (Charm N-body GrAvity solver) (Jetley et al., 2008, 2010; Menon et al., 2015). We call this moving-mesh solver MANGA (Chang et al., 2017, hereinafter C17). The solver is based on the scheme described by S10, but also utilizes advances in gradient estimation (Steinberg et al., 2016) and limiters (Duffell & MacFadyen, 2011). We also use an alternative method to construct the Voronoi tessellation (Rycroft 2009; C17). More recently, we have added magnetohydrodynamics (Chang, in preparation) and radiation hydrodynamics (Chang, Davis, & Jiang in preparation) to MANGA.

n addition, CHANGA uses the CHARM++ language and run-time system (Kale & Krishnan, 1996) for parallelization rather than a custom message-passing interface design.111Other astrophysical codes that use CHARM++ include ENZO-P and FVMHD3D (Gaburov et al., 2012). The use of CHARM++ promises that CHANGA will be much more scalable than previous astrophysical codes. In particular, CHANGA has demonstrated strong scaling on single time-stepping problems with $12$ billion particles to 512K cores (with 93% efficiency) and on multi-time-stepping problems with 52 million particles to 128K cores (Menon et al., 2015).

Given the importance of CEE and the recent advances in numerical computation, we apply MANGA to study CEE using MM techniques. In this initial work, we describe our methodology, highlighting the algorithmic improvements to MANGA to enable this study. In particular, we describe the implementation of individual time-steps, the incorporation of a realistic equation of state and the generation of realistic stars and initial conditions. We also show one example CEE evolution computed both with and without the giant set into corotation with the binary orbit.

We have organized this initial paper as follows. In section 2.1, we summarize the algorithm of MANGA, highlighting improvements made since C17. We also clarify the half time-step prediction that was not well described by C17. We describe the implementation of multistepping in section 2.2 and realistic equations of state in section 2.3. We describe how we use the stellar evolution code MESA to produce 1-D models of red giants that we map into 3-D unstructured grids in section 2.4 and use these models as initial conditions for CEE. In section 3, we describe the results of our simulations of a 2 $\textrm{M}_{\odot}$ red giant and a 1 $\textrm{M}_{\odot}$ companion. We discuss our results and close in section 4.

2 Methodology

2.1 Numerical Approach

We begin with a brief summary of MANGA. This is mainly covered by C17 but the algorithm has been modified since C17 which we highlight below. MANGA solves the the Euler equations, e.g., the conservation of mass equation, conservation of momentum and conservation of energy, which written in conservative form are

[TABLE]

and

[TABLE]

where $\rho$ is the density, $\boldsymbol{v}$ is the velocity, $\Phi$ is the gravitational potential, $e=\epsilon+v^{2}/2$ is the specific energy, $\epsilon$ is the internal energy, and $P(\rho,\epsilon)$ is the pressure. Equations (1) to (3) can be written in a compact form by introducing a state vector $\boldsymbol{\mathcal{U}}=(\rho,\rho\boldsymbol{v},\rho e)$:

[TABLE]

where $\boldsymbol{\mathcal{F}}=(\rho\boldsymbol{v},\rho\boldsymbol{v}\boldsymbol{v},(\rho e+P)\boldsymbol{v})$ is the flux function, $\boldsymbol{\mathcal{S}}=(0,-\rho\boldsymbol{\nabla}\Phi,-\rho\boldsymbol{v}\cdot\boldsymbol{\nabla}\Phi$) is the source function and $V$ is an elemental volume.

To solve equation (4), we adopt the same finite volume strategy as S10. We refer the interested reader to S10 for a more detailed discussion of the scheme. Here, we only briefly describe the scheme to document the algorithm we have implemented and to highlight the differences between our scheme and that of S10.

For each cell, the integral over the volume of the $i$th cell defines the charge of the $i$th cell, $\boldsymbol{U}_{i}$, to be

[TABLE]

where $V_{i}$ is the volume of the cell. As do S10, we then use Gauss’ theorem to convert the volume integral over the divergence of the flux in equation (4) to a surface integral

[TABLE]

where $\partial V_{i}$ is the boundary of the cell. We now take advantage of the fact that the volumes are Voronoi cells with a finite number of neighbours to define an integrated flux

[TABLE]

where $\boldsymbol{F}_{ij}$ and $A_{ij}$ are the average flux and area of the common face between cells $i$ and $j$. The discrete time evolution of the charges in the system is given by

[TABLE]

where $\hat{\boldsymbol{F}}_{ij}$ is an estimate of the half time-step flux between the initial $\boldsymbol{U}_{i}^{n}$ and final states $\boldsymbol{U}^{n+1}_{i}$ and $\boldsymbol{S}_{i}^{(n+1/2)}=\int_{i}\boldsymbol{\mathcal{S}}dV$ is the time-averaged integrated source function.

We estimate the flux $\hat{\boldsymbol{F}}_{ij}$ across each face as follows.

1. (i)

Use the gradient estimates at the initial time-step to predict the half time-step cell centered values. 2. (ii)

Drift the cells a half time-step and rebuild the Voronoi tessellation at the half time-step. This step is new compared to C17. 3. (iii)

Estimate the half time-step state vector (in the rest frame of the moving face) at the face centre ($\boldsymbol{\tilde{r}}_{ij}$) between the neighboring $i$ and $j$ cells by linear reconstruction. 4. (iv)

Estimate the (half time-step) velocity $\boldsymbol{\tilde{w}}_{ij}$ of the face, following S10, and boost the state vector from the “lab” frame (the rest frame of the simulation box) to the rest frame of the face to find the flux along the normal of the face. I.e., in the direction from $i$ to $j$. 5. (v)

Estimate the flux $\hat{\boldsymbol{F}}_{ij}$ across the face using an HLL or HLLC (or HLLD for MHD, Chang, in preparation) Riemann solver implemented following Toro (2009). 6. (vi)

Boost the solved flux back into the “lab” frame.

We can then use the estimated fluxes to time-evolve the charges $\boldsymbol{U}_{i}$ following equation (8) and using the full time-step $\delta t$ and apply changes owing to the source terms.

2.1.1 Half Time-Step Predictions on a Moving Voronoi Mesh

We note from the above that we compute the half time-step Voronoi cells to achieve true second-order time integration. This is a major change from C17, who used the initial time-step cells. Here we derive the estimate for the time derivative for the half time-step prediction step as this estimate was insufficiently described by C17.

The Euler equations written in conservative form are equations (1) to (3). Following C17, let’s assume that the mesh generating point moves with velocity $\boldsymbol{w}$. Boosting (in a Galilean sense) into the reference frame of the moving point, equations (1) to (3) become

[TABLE]

and

[TABLE]

where $\boldsymbol{v}^{\prime}=\boldsymbol{v}-\boldsymbol{w}$ and $e^{\prime}=\epsilon+v^{\prime 2}/2$. Expanding this out, we can write

[TABLE]

and

[TABLE]

Equations (12) to (14) provide estimates for the time derivatives of the conserved variables, $\rho$, $\rho\boldsymbol{v}$ and $\rho e$ with the estimated gradients for the conserved variables calculated as described by C17. This allows us to formally achieve the same second order level of accuracy in time integration as Pakmor et al. (2016); (see also Duffell & MacFadyen, 2011), but there are some differences which we discuss below.

2.2 Multistepping

MANGA as described by C17 imposes a universal time-step. Large speedups for problems with large dynamic range are possible with individual time-steps. Hence, to greatly enhance the applicability of MANGA to the problem of CEE, we have implemented an individual time-step scheme for MANGA.

The basic (universal) time-stepping algorithm for a second order accurate (in time) integrator can be broken down into 3 stages, the initial time (a), the half time-step (b) and the full time-step (c).

1. (a)

Initial time $t=0$: determine Voronoi cells using current positions of mesh-generating points. Estimate gradients and construct half time-step predictions. Zero out the changes to the charges, e.g., $\delta\boldsymbol{U}=0$. 2. (b)

Half time-step $t=0.5\delta t$: drift Voronoi cells to half time-step positions. Construct Voronoi cells at half time-steps. Perform linear reconstruction to the half time-step cell faces and compute fluxes. Perform a Riemann solution and incorporate source terms with the full time-step. Place these changes into $\delta\boldsymbol{U}$. 3. (c)

*Full time-step $t=\delta t$: * drift Voronoi cells to full time-step positions. Advance charges to be $\boldsymbol{U}^{n+1}=\boldsymbol{U}^{n}+\delta\boldsymbol{U}$. Reset the state to be the new initial time.

In comparison, Pakmor et al. (2016) achieve the same formal second order accuracy in two steps, an initial and full time-step. Here, they make a full time-step prediction for the primitives, perform a Riemann solution using the initial Voronoi cells and faces for a half time-step and perform a closing Riemann solution (again a half time-step) using the final Voronoi cells and faces at the full time-step and the full time-step prediction for the primitives. The two methods differ in the comparative numbers of Voronoi tessellations and Riemann solutions. In particular, we only need to solve the Riemann problem once, whereas Pakmor et al. (2016) require two solutions of the Riemann problem (with associated linear reconstructions). However, we require an extra Voronoi calculation at the half time-step. Depending on the relative cost of Riemann solution and the Voronoi calculation – expensive for tabulated EOSs – the two schemes have their relative merits.

To adapt our scheme for multistepping, we note that the full time-step for any cell involves actions at a full time-step (a,c) and the half time-step (b). Let us then consider the time advancement of three separate levels or rungs at 0,1,2, where the time-step decreases by a factor of 2 at every level or rung in Fig. 1. In the interest of keeping terminology between the code and this paper synchronized, we drop the term level and refer to the different time-step levels as rungs.

In Fig. 1, the full time-steps are represented by filled circles, while the half time-steps are represented by empty circles. At step 0, all three rungs are synchronized at their respective full time-steps: all the circles are black and step (a) is executed on all three rungs. The particles are then drifted $0.25\delta t$ to step 1. The half time-step at rung 2, and the half time-step reconstruction and Riemann solution is performed for rung 2. Note that the Riemann solution is for the full time-step $0.5\delta t$ at rung 2 so that when the $\delta\boldsymbol{U}$ is summed at the next step, the $\boldsymbol{U}$s for rung 2 are updated to the full time-step. Now at step 2, the full time-steps (c) and initial time-steps (a) are performed at rung 2 and the half time-step (b) is performed for rung 1. In addition, cells at rung 0 that share a face with a cell at rung 1 (denoted by red-dashed circles) also have a Riemann solution at the common face at the rung 1 half time-step. We skip step 3 as it is the same as step 1 and now proceed to step 4. Here, we do the full time-steps for rungs 1 and 2 and the half time-step for rung 0.

From this discussion, we can conclude that for the half time-step for rung $i$, full time-steps are executed from all rungs $j>i$. Conversely, there can only be one half time-step for only one rung per (drift) time-step. This is the half time-step of the highest rung. In addition, cells only switch rungs during their full time-steps to ensure consistency of the second-order time integration and they can only reduce their rung when their new reduced rung is at its full time-step.

One additional feature that we have included in our implementations is the smoothing of time-steps by large propagating velocity fluctuations, or shocks. This is a well known issue in SPH, for which Saitoh & Makino (2009) show that large localized changes in the time-step can lead to incorrect results. They also showed that smoothing out time-stepping changes over a SPH kernel such that the time-step variation is limited to a factor of 4 mitigates these issues. We have adopted a similar approach here in that we limit the cell’s time-step to be less than $\sqrt{2}$ of the minimum time-step of its neighbours. This results in time-step changes of a factor of two occurring over two cells away from a minimum cell. We find that this smoothing allows for stable integrations in MANGA.

This multistepping algorithm was applied to a simple CEE simulation to test the speedup factor. At each step, we find speedups ranging from 3 to 9 times that of the universal time-stepping algorithm, resulting in an average speedup of 4 to 5 times.

2.3 Realistic Equation of State

In our original implementation of MANGA described by C17, we only used an ideal, adiabatic equation of state. However, it is straightforward to extend MANGA to arbitrary equations of state. For CEE, we have elected to use the equations of state implemented in the open source stellar evolution code, MESA (Paxton et al., 2011; Paxton et al., 2013, 2015). The MESA equation of state relies on several equations of state relevant over the different regimes of stellar structure, allowing MESA to span the full range from giant planets and white dwarfs to massive stars (Paxton et al., 2011).

For a hydrodynamic solver, the equation of state can be called many times – during the update of cell centred variables, during the reconstruction step, during the flux calculation and during the Riemann solution (to determine sound speeds). We have found that tabulated equations of state are much too slow to be called so many times during a single solve. In particular, an average Voronoi cell has about 20 faces over which a Riemann solver must be executed. This means at least 20 calls to the equation of state to find the sound speed (for linear or higher order reconstruction). Instead, we only make one call per step in MANGA – during the update of cell centred variables. To find values of these variables on the faces, we estimate the relevant thermodynamic quantities using linear reconstruction of the density, energy and adiabatic index. We have found that this results in a factor of 25 speedup in the Riemann solver.

2.4 Initial Conditions

We now discuss the development of appropriate initial conditions for MANGA. As stated in the Introduction, Ohlmann et al. (2016) previously performed a similar moving-mesh simulation of CEE. Here, we adopt similar initial conditions. We first use MESA to evolve a 2 $\textrm{M}_{\odot}$ star with metallicity $Z=0.02$ from the pre-main sequence to the red giant phase. We stop when the star reaches 52 $\textrm{R}_{\odot}$ with a He core mass of 0.36 $\textrm{M}_{\odot}$. Compared with that of Ohlmann et al. (2016) our red giant is slightly larger and its He core mass is slightly smaller. This could be due to the updates in the MESA code since their original result or discrepancies in our MESA inlists. From the MESA output, we take the entropy and hydrogen fraction. For the core, we take the total mass at a density that is $50$ times greater than the mean density of the red giant, giving a core mass $M_{\rm{c}}=0.379\textrm{M}_{\odot}$. This corresponds to a core radius $R_{\rm{c}}=1.99\textrm{R}_{\odot}$, which we use as the core gravitational softening length. Because of the great difference in density between the He core and the H envelope, we model the core as a dark matter particle with mass $M_{\rm{c}}$ and softening radius of $R_{\rm{c}}$ similar to the initialization of Ohlmann et al. (2016)’s simulations. We then take the entropy profile and construct a star of mass $M-M_{\rm{c}}$, with an entropy profile which matches that of the original star and contains a dark matter particle core. This yields a radial profile of density, temperature and H-fraction that can be mapped to a particle (mesh generating point) profile.

We construct an appropriate particle mesh for the star from a precomputed glass distribution of points embedded in a 3-D cube. We periodically replicate this glass distribution to produce sufficient numbers of particles. We assume that each particle is of equal mass and rescale them to the appropriate radial position based on the computed $M(r)$ from MESA. These particles are also endowed with the radially interpolated temperature and H-fraction. The total number of particles representing the star is $3\times 10^{5}$, which is smaller than the $2\times 10^{6}$ particles used by Ohlmann et al. (2016). Outside of the star, we include a low density atmosphere of $10^{-13}$ g cm*-3* with temperature $10^{5}$ K that extends out to the total box size of $3.5\times 10^{14}$ cm ($5000~{}\textrm{R}_{\odot}$), with periodic boundary conditions at its edges. The total number of particles in the simulation is $8\times 10^{5}$.

To lower the computational cost, we use a mesh refinement algorithm to decrease the number of gas particles in the atmosphere far from the star. We define a scale factor $\mathcal{S}(r)=(r/R_{*})^{n}$ where $R_{*}$ is the radius of the star, $r$ is the spherical radius and $n$ is an adjustable parameter which we have set to $n=2/3$ in this case. Starting with the same uniform glass distribution as for the star, the linear spacing between particles is increased by $\mathcal{S}$ and their mass is increased by $\mathcal{S}^{3}$, preserving the external density.

The profiles produced from the mapping from 1-D stellar evolution models to fully 3-D hydrodynamics simulations are not in perfect hydrostatic balance owing to discretization errors. Previously C17 used this fact to compare (roughly) hydrostatic stars experiencing oscillations using the MM and SPH solver. Here we are interested in the equilibrium solution in the MM solvers. Thus, following the mapping from 1-D to 3-D, we damp spurious velocities and energies and allow the star to relax to a final configuration. We note that this is the same procedure that we have used previously for merging white dwarfs (Zhu et al., 2013) and was also used by Ohlmann et al. (2016) for giant star initial conditions.

We show the result of this relaxation in Fig. 2 as blue circles and compare this to the radial profile from MESA (orange circles). Here we downsample the number of points from MANGA by a factor of 100. The horizontal line denotes the point where the density in the MESA profile is 50 times greater than the mean density. The mass inside of this radius is represented by a dark matter particle with a gravitational softening length equal to the radius. The star is then recomputed to follow the same entropy–mass relation (corrected for the inclusion of the dark matter particle) and then relaxed in MANGA. The envelope follows the original MESA profile quite faithfully in spite of the differences in the physical structure. The inner profile approaches a fixed density asymptotically owing to the softened gravitational forces from the dark matter core.

We represent the 1 $\textrm{M}_{\odot}$ companion as a dark matter particle that is in an initially circular Keplerian orbit at the red giant radius $a=52\textrm{R}_{\odot}$. This follows the initial conditions of Ohlmann et al. (2016) to facilitate a direct comparison. Although this neglects the evolution of the binary prior to CEE, we compensate by altering the dynamics of the giant. During the phase leading up to CEE, tidal forces from the companion on the envelope are expected to spin up the rotation of the envelope (Soker, 1996) to a significant fraction of its breakup velocity. We have implemented this corotation between the envelope rotation and orbital motion into our simulations by following the scheme of MacLeod et al. (2018). Within the envelope, we assume rigid body rotation and initialize the velocity as

[TABLE]

but give the atmosphere a velocity

[TABLE]

Here, $f_{\rm{cr}}$ is an adjustable parameter, $\phi$ is the azimuthal angle, $\theta$ is the polar angle, $\Omega_{\rm{orb}}$ is the orbital frequency of the red giant and companion, $R_{*}$ is the radius of the giant and $R_{\rm{cyl}}$ is the cylindrical radius from the rotation axis of the giant. We note that (15) and (16) ensure that the velocity is continuous at the surface of the giant.

Here, we run two simulations in which the giant is in different states of corotation, 95 per cent corotation ($f_{cr}=0.95$) and 0 per cent corotation ($f_{cr}=0$). We note that the 95 per cent corotation case follows the simulations of Ohlmann et al. (2016) and Ricker & Taam (2012). The choice to set $f_{cr}$ to less than unity is also motivated by MacLeod et al. (2018) who found that orbital angular momentum lost by the companion prior to the onset of a common envelope phase does not necessarily go into spinning up the envelope, leading to a desynchronization between the orbital frequencies of the companion and envelope.

3 Results

Starting with the above initial conditions, we simulate the binary for 240 d and show several density projections at $t=$ 1, 10, 30, 75, 120 and 240 d in Figs. 3 and 4. We note that the simulation period is about 10 times the initial orbital period of 25 d.

In both simulations, the companion experiences an initial plunge into the envelope of the giant that lasts about 15 days, throwing off a large tidal tail (Figs. 3 and 4 upper centre) and greatly decreasing the separation between the companion and core. The companions then continue to spiral in as their orbital energy is transferred to the gas; the spiral shocks facilitating this transfer can be seen in the projections (Figs. 3 and 4 upper right). We see smaller tidal tails thrown off at $t=75$ d (lower left) and $t=120$ d (lower centre), as well as several others. The simulation ends well before the outflow reaches the edge of the simulation box.

3.1 Envelope Ejection

The total energy of each gas particle is given by

[TABLE]

(Nandez et al., 2015), where $m_{i}$, $\boldsymbol{v}_{i}$, $\boldsymbol{v}_{\rm CM}$, $\phi_{i}$ and $I_{e,i}$ are the mass, velocity, centre of mass velocity of the bound material, gravitational potential and specific internal energy of each particle. Gas particles with a negative total energy are bound to the binary, while those with a positive total energy are unbound. Fig. 5 shows a density slice in the $z=0$ plane at $t=15$ d in the 0 per cent case. The black contour encloses the matter bound to the binary.

The kinetic energy is computed relative to the velocity of the centre of mass (CM) of the bound matter; that is, the bound gas as well as both dark matter particles. However, because the total energy is needed in order to determine which gas particles are bound, we used an iterative scheme to find the velocity of the CM.

1. (i)

The CM velocity is assumed to be zero. 2. (ii)

The velocities of the gas particles are calculated relative to the CM velocity. 3. (iii)

These velocities are used to find the kinetic energies of the gas particles, which are used to determine which particles are bound and which are not. 4. (iv)

A new CM velocity $\boldsymbol{v}_{\rm CM}$ is found for everything except the unbound particles,

[TABLE]

where $m_{\rm c,1}$ and $m_{\rm c,2}$ denote the dark matter core and companion, respectively, $m_{\rm g,b}$ denotes the bound gas, m_{\rm c,1}$$\boldsymbol{v}_{1}, m_{\rm c,2}$$\boldsymbol{v}_{2}, and m_{\rm g,b}$$\boldsymbol{v}_{\rm g,b} and $\boldsymbol{v}_{1}$, $\boldsymbol{v}_{2}$, and $\boldsymbol{v}_{\rm g,b}$ are the momenta and velocities of the dark matter core and companion and bound gas, respectively. 5. (v)

Steps (ii) to (iv) are repeated until $\boldsymbol{v}_{\rm CM}$ reaches convergence.

This is similar to the method of Guillochon & Ramirez-Ruiz (2013). With the total energies of all particles known, we can find the mass of the unbound gas as a fraction of the total mass of the envelope $m_{\rm env}$,

[TABLE]

which we refer to as the ejection efficiency. This is shown in Fig. 6. Qualitatively, the simulations are very similar in this regard. They show a large ejection of gas during the initial plunge as well as during the ejection of additional tidal tails. Unsurprisingly, the system set into corotation consistently exhibits a higher ejection efficiency owing to the inclusion of additional energy and angular momentum. We find that matter continues to be ejected up until the end of the simulation period, ending with an efficiency of 66 per cent with corotation and 63 per cent for no corotation.

We also show in Fig. 6 the ejection efficiency that would be obtained if the internal energies were neglected. We find that only 8 per cent of the envelope is unbound in both cases, which is similar to the reported ejection efficiency in Ohlmann et al. (2016). Such a situation can be reached if the envelope expands sufficiently such that it loses energy rapidly by radiation. This indicates that the internal energy plays a large role in unbinding the envelope. It also suggests that further work on CEE may not be able to ignore radiative cooling because of its possible large effects.

3.2 Orbital Parameters

The separations between the stellar core and companion are shown in Fig. 7. For ease of reading, we also show the separations on a logarithmic scale and smooth them over a period of 15 d. That is, each separation $r_{i}$ is taken to be the average of all separations within an interval of 15 d centred on $r_{i}$.

After one orbit, both binaries have reduced their orbital separation to less than half the initial separation, and they continue to spiral in at a slower rate. Although the simulation with corotation initially falls to a smaller separation, it is eventually surpassed by its counterpart. The final (smoothed) separations are 3.6 $\textrm{R}_{\odot}$ with corotation and 3.2 $\textrm{R}_{\odot}$ with no corotation.

The eccentricity $e$ can be found from the periapsis and apoapsis distances,

[TABLE]

where $r_{\rm apo}$ and $r_{\rm peri}$ denote the distances of closest and furthest approach. Here $r_{\rm apo}$ and $r_{\rm peri}$ are determined by interpolating the maxima and minima of the separation. The alert reader will notice that this definition is meaningful in this situation because the potential is not Keplerian owing to the distribution and nonaxisymmetry of the unbound gas. However, this method allows the determination of the eccentricity from only the orbital separations. Both binaries circularize their orbits over time, seeing drops in eccentricity from $e\approx 0.2$ down to $e=0.1$ by the end of the simulation period.

In addition, each binary receives a kick velocity owing to the ejection of gas from the system, as shown in Fig. 8. Interestingly, although corotation results in a larger mass ejection and a large initial spike in CM velocity during the plunge, it quickly drops below the non-rotating case by about 1 km s*-1*.

4 Discussion and Conclusions

In this initial work, we describe the methodology for simulations of CEE using the moving-mesh code MANGA. First, we detail a new multistep scheme to speed up problems with large dynamic range. In the case of CEE, we find speedups of up to 4 to 5 times over our previous single-stepping scheme. We then implement a realistic equation of state using the libraries in the stellar evolution code MESA and discuss the construction of realistic initial conditions.

Using these initial conditions, we run two simulations of a 2 $\textrm{M}_{\odot}$ red giant with a 1 $\textrm{M}_{\odot}$ companion (thought to be a main sequence star), where we take the rotation of the envelope to be either 95 or 0 per cent of the (initial) orbital corotation. After 240 d, we find that the orbital separation of the red giant core and the companion has been reduced to 3.6 $\textrm{R}_{\odot}$ and 3.2 $\textrm{R}_{\odot}$, respectively. This is a reduction in the separation by a factor of 15. Most of the envelope is ejected (66 and 63 per cent respectively).

Our envelope ejection efficiencies are much larger than many of the results of previous work. However, three-dimensional numerical computations of CEE have found a wide variance in the ejection efficiency of the envelope. Early work (for example, de Kool, 1987; Livio & Soker, 1988; Terman et al., 1994; Sandquist et al., 1998; Sandquist et al., 2000; Terman et al., 1995; Rasio & Livio, 1996) on 3-D numerical simulations of CEE found envelope ejection efficiencies from about 10 to 80 per cent. However, the resolutions of these early simulations were poor. More recent work Passy et al. (2012); Ricker & Taam (2012); Nandez et al. (2015); Ohlmann et al. (2016) also find a large variance in the efficiencies (about 5 to 100 per cent). For instance, Passy et al. (2012) found that only 15 per cent of the mass was unbound in their simulations. Ricker & Taam (2012) found that about 25 per cent of the mass is unbound, but the unbound mass is still rising at the end of their simulations (see their Fig. 9). Nandez et al. (2015) found that between 50 and 100 per cent of the envelope mass is ejected in their simulations, depending on the inclusion of recombination energy. Ohlmann et al. (2016) found that about 8 per cent of the envelope mass is ejected. Our results are on the higher side and are similar to what Nandez et al. (2015) find. One difficulty in comparing these results is that the initial conditions and simulation parameters are not similar in these different groups. In this paper, we have adopted initial conditions and parameters similar to those used by Ohlmann et al. (2016) but our results are dramatically different. The difference may be attributed to different simulation parameters such as the gravitational softening, the parameterization of the dark matter particle that represents the core of the red giant and the companion, or our inclusion of the recombination energy though the use of the MESA EOS. Indeed, we have found that the internal energy has a large impact on the ejection efficiency. This suggests that the incorporation of radiation physics may be an important next step.

It has also been suggested that inclusion of recombination energy would result in more efficient expulsion of the envelope (see the review by Ivanova, 2017). We also note that this point is controversial, as Soker and collaborators have claimed in a series of papers (Sabach et al., 2017; Soker, 2017; Soker et al., 2018; Grichener et al., 2018) that the energy from recombination is easily transported away either by radiation or convection. However, Ivanova (2018) reached the opposite conclusion, that the fraction of recombination energy transported away from the regions in which recombination takes place is negligible. The importance of radiative cooling in this case highlights the need for a more careful accounting of radiation effects in future calculations.

Our work here reinforces the results of Nandez et al. (2015) which suggests that the entire envelope of the red giant might be ejected during the plunge phase of CEE or at least can be numerically simulated. Thus, the effects or relevance of the other phases of CEE such as the self-regulated spiral-in phase is unclear. It may be the case that these phases are important or relevant for different system parameters, though this requires a more detailed exploration of parameter space.

Finally, we note here that other processes may also aid the expulsion of the envelope. Much recent work have been focused on energy and momentum injection from jets produced by accretion onto the companion (Papish et al., 2015; Moreno Méndez et al., 2017; Shiber & Soker, 2018; Chamandy et al., 2018). In addition, the transition from laminar flow to accretion flows around the companion star can be complicated and also provide another mechanism by which orbital energy is dissipated (MacLeod & Ramirez-Ruiz, 2015a, b; MacLeod et al., 2017). Much of the work have focused on Bondi-Littleton-Hoyle type accretion around a point mass but a typical accretor has a physical scale of a solar radius, which gives it a nontrivial cross section in a red giant envelope. Typically, we can show that a 1 $R_{\odot}$ accretor encounters 10 to 50 per cent of the envelope in one orbit. This hydrodynamic interaction may play a significant role in the evolution of the common envelope and the binary orbit. Finally, we mention that it is also suggested that dust formation during CEE can drive winds to expel the envelope (Soker, 1994; Glanz & Perets, 2018).

Acknowledgements

We acknowledge useful conversations with Thomas Quinn and Volker Springel during the development of the multistepping algorithm and refinement on the half time-step prediction. We also are grateful to Bill Paxton for tracking down a MESA bug in release 10398, supplying a workaround, and correcting it in future versions of MESA. We thank S. Ohlmann and R. Pakmor for useful discussions. We thank the reviewer, Christopher Tout, for a useful and constructive referee report. We are supported by the NASA ATP program through NASA grants NNH17ZDA001N-ATP, NSF CAREER grant AST-1255469, and the Simons Foundation. We also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (NSF) grant No. ACI-1053575. We also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper (URL: http://www.tacc.utexas.edu). We also use the YT software platform for the analysis of the data and generation of plots in this work (Turk et al., 2011).