Revisiting the common envelope evolution in binary stars: A new semianalytic model for $N$-body and population synthesis codes

We present a novel way of modeling common envelope evolution in binary and few-body systems. We consider the common envelope inspiral as driven by a drag force with a power-law dependence in relative distance and velocity. The orbital motion is resolved either by direct $N$-body integration or by solving the set of differential equations for the orbital elements as derived using perturbation theory. Our formalism can model the eccentricity during the common envelope inspiral, and it gives results consistent with smoothed particles hydrodynamical simulations. We apply our formalism to common envelope events from binary population synthesis models and find that the final eccentricity distribution resembles the observed distribution of post-common-envelope binaries. Our model can be used for time-resolved common-envelope evolution in population synthesis calculations or as part of binary interactions in direct $N$-body simulations of star clusters.

Figure 1

Evolution of semimajor axis ($a$), eccentricity $e$, argument of pericenter $\omega$ and true anomaly $\nu$ for a CE inspiral using the model in Eq. (9) with $l=2$, $k=0$ (top panel) and $l=1$, $k=0$ (bottom panel). The blue curves were obtained integrating Eqs. (10)–(13), while the orange curves are the result of direct $N$-body integration with the drag force applied as a perturbation. The two curves overlap to a great extent.

Figure 2

Evolution of semimajor axis $a$, eccentricity $e$, argument of pericenter $\omega$. For a CE inspiral using the model in Eq. (9) with $l=2$, $k=1$ (top panel), $l=1$, $k=2$ (middle panel) and $l=1$, $k=2$ (bottom panel). Each curve indicates a different starting eccentricity.

Figure 3

Evolution of semimajor axis (top) and eccentricity (bottom) as a function of time, assuming self-similar expansion of the envelope. The colors indicate different initial binding energy of the envelope, parametrized by the $\lambda$ parameter [Eq. (19)]. Blue: no envelope expansion. Orange: $\lambda =0.2$. Green: $\lambda =0.5$. Red: $\lambda =1$. Purple: $\lambda =2$. We set ${\chi }_a=0.05$.

Figure 4

Comparison of the binary separation evolution for four hydrodynamical simulations and the power-law drag force model. From top to bottom, left to right: simulations 8r2g5, 8r2g-0, 1r06p5, 1r06p7 from Glanz and Perets [78]. For clarity, we compared models in 8r2g-0 from the second apocenter approach, where the CE rapid plunge-in began. The blue line is the separation obtained from the hydrodynamical simulations, while the thick purple line is the power-law model. All the semianalytic curves use the quadratic drag force ($l=2$). The left panels use a drag force linearly decreasing with radius ($k=1$), while the right panel adopts a power-law of $k=3/2$. The top-right panel shows the difference between $k=1$ and $k=3/2$.

Figure 5

Semimajor axis (top), eccentricity (middle), and separation (bottom) of the inner binary of a triple stellar system as a function of time. Solid lines: simulation including a drag force term with $l$, $k=2$, 0 and ${\chi }_a=0.005$. Dotted line: simulation with no drag force. We stop the CE simulation as soon as the energy released by the drag equals the binding energy of the envelope. Evolved by means of direct $N$-body integration with an Hermite scheme [66, 67].

Figure 6

Pericenter distance (left) and eccentricity (right) distributions, pre-CE (blue) and post-CE (red). In the left panel, the post-CE semimajor axes distribution derived from bse is outlined in green. Using the drag force model with ${\chi }_a=0.05$, $l=2$, $k=1$, and stopping the inspiral to be consistent with the $\alpha \lambda$ model, as described in Sec. 3a. The common envelope events are taken from the BPS calculations of Tanikawa et al. [110]. In BPS codes, the final eccentricities are in practice always zero, while we obtain small but finite eccentricities in our model.

Figure 7

Comparison between phase-dependent [Eqs. (10)–(13)] versus orbit-averaged evolution [Eqs. (b5)–(b6)] for the case $l=2$, $k=0$. Top panel: semimajor axis. Bottom panel: eccentricity. Both integrations are carried out with a Dormand-Prince 8th order integrator. In both cases, ${\chi }_a=0.05$.

Figure 8

Values of the integral $-I(e)$ that arises when orbit averaging Eq. (c1), that is the rate of change in eccentricity for a dynamical friction force ($l=-2$). Each curve corresponds to a different value of $k$. For $k=3$, $⟨\dot{e}⟩=0$ and the derivative of the eccentricity changes sign, from positive ($k>3$) to negative ($k<3$).