Kinematics of Mass Loss from the Outer Lagrange Point L2

TL;DR

This paper analyzes the kinematics of mass loss from the L2 point in binary systems, revealing how initial conditions influence whether material escapes, forms a disk, or impacts the binary, aiding interpretation of simulation results.

Contribution

It provides a detailed characterization of particle trajectories near L2, highlighting the effects of initial velocities and offsets on mass loss outcomes in binary systems.

Findings

Unbound outflows can occur even with initial velocities below corotation.

Trajectory outcomes depend on initial position offsets and velocity vectors.

Results help interpret hydrodynamic simulation morphologies.

Abstract

We investigate kinematics of mass loss from the vicinity of the second Lagrange point L2 with applications to merging binary stars, common envelope evolution and the associated transient brightenings. For ballistic particle trajectories, we characterize initial velocities and positional offsets from L2 which lead to unbound outflow, fall back followed by a formation of a decretion disk, collision with the binary surface, or a hydrodynamic shock close to the binary, where some particle trajectories loop and self-intersect. The latter two final states occur only when the trajectories are initiated with offset from L2 or with velocity vector different from corotation with the binary. We find that competition between the time-dependent and steeply radially decreasing tidal torques from the binary, Coriolis force and initial conditions lead to a non-trivial distribution of outcomes in the vicinity of L2. Specifically, even for initial velocities slower than corotation, we find that a set of initial position offsets lead to unbound outflows. Our results will aid in the interpretation of the morphology of mass loss streams in hydrodynamic simulations.

Equations

$$\varphi =-\frac{\mu }{\sqrt[]{(x-1+\mu {)}^2+y^2}}-\frac{1-\mu }{\sqrt[]{(x+\mu {)}^2+y^2}}-\frac{1}{2}\left(x^2+y^2\right),$$

$$f{x}_{L2}-\frac{\mu }{{\left(x_{L2}-1+\mu \right)}^2}-\frac{(1-\mu )}{{\left(x_{L2}+\mu \right)}^2}=0.$$

$$\ddot{x}$$

$$\ddot{y}$$

$$E$$

$$J$$

$$C=\varphi +\frac{1}{2}\left({\dot{x}}^2+{\dot{y}}^2\right)$$

Full text

††pubyear: 2019††pagerange: Kinematics of Mass Loss from the Outer Lagrange Point L2–LABEL:lastpage

1 Introduction

The evolution of binary stars, including exoplanets and their host stars, can be significantly influenced by the mass loss from the system. When the mass loss from a binary with semi-major axis $a$ and total mass $M$ occurs with velocities comparable to or lower than the binary escape velocity $v_{\rm esc}=\sqrt{2GM/a}$ the morphology of the outflow is significantly affected by the gravitational forces of the system. In particular, the outer critical Lagrange point L2 plays a prominent role. The L2 point is the saddle point of the Roche potential located on the binary axis just outside the less massive component (Fig. 1). At L2, the gravitational and centrifugal accelerations in the corotating reference frame cancel out. This means that a only a small perturbation is required to eject an initially corotating test particle. Typical situations where this might arise are gradual expansion of the contact binary envelope, which was studied by Webbink (1976, 1977) for low-mass stars and Flannery & Ulrich (1977) for high-mass stars. Mass loss from L2 can also result from rapid mass transfer between the binary components (e.g., Livio et al., 1979; MacLeod et al., 2018a). In many cases, mass and angular momentum loss from L2 causes rapid shrinking of the orbit so that the binary soon experiences a common envelope evolution episode and might even merge into a single object (e.g., Paczynski, 1976).

Recently, Pejcha (2014) and Pejcha et al. (2017) argued that L2 mass loss lasting many hundreds and thousands of orbital periods was responsible for the pre-merger behaviour of V1309 Sco – a contact binary that exhibited accelerating orbital period decay accompanied by changes in the light curve shape and overall brightening preceding its luminous red nova outburst (Tylenda et al., 2011). L2 mass loss or a similar process could lead to an equatorially-concentrated mass distribution that is later overtaken by a faster and more spherical explosion as is implicated in many supernovae, classical novae, luminous blue outbursts, and other transients (e.g., Li et al., 2017; Andrews & Smith, 2018; Smith et al., 2018). Furthermore, mass loss from the vicinity of L2 can feed the circumbinary disk of microquasar SS433 (Fabrika, 1993) and may significantly affect the gravitational wave frequency and amplitude of extreme mass-ratio inspirals (Linial & Sari, 2017).

Material initially positioned in corotation at L2 is always bound. However, after leaving L2 the time-changing gravitational field (i.e., binary tidal torques) transfer angular momentum and energy from the binary orbit to the L2 outflow (Kuiper, 1941). When considering purely ballistic motions of test particles, Shu et al. (1979) found that the final fate of L2 mass loss depends only on the binary mass ratio $q\leq 1$. For $0.064\leq q\leq 0.78$, the tidal torqueing is efficient enough that the L2 outflow achieves positive energy and leaves the binary to infinity with asymptotic velocities $v\lesssim 0.25v_{\rm esc}$. For $q<0.064$, tidal torqueing is inefficient owing to the small mass of the secondary component, and for $q>0.78$, the initial energy of test particles corotating at L2 is too negative and the particle remains bound. For such a $q$, the L2 ballistic outflow does not achieve positive energy, falls back to the binary and establishes a viscous decretion disk. The picture becomes more complex when the ballistic particles are replaced by a realistic gaseous medium with shocks and radiative processes such as diffusion and cooling. Pejcha et al. (2016a, b) found that, depending on the radiative cooling efficiency and the initial temperature, L2 mass loss leads to an equatorial outflow with a narrow or wide opening angle, a nearly isotropic outflow, a decretion disk or an inflated envelope.

Despite these advances, most of the work focusing specifically on L2 mass loss has so far been done assuming corotation and point-like injection at L2. This assumption must eventually become invalid when the timescale of orbital decay becomes shorter than the timescale for maintaining stellar rotation synchronized with the orbit111Since corotation is easier to maintain only in a relatively thin surface layer, the synchronization time-scale relevant for L2 mass loss can be much shorter than the tidal synchronization time-scale of the entire star., as is likely the case for common envelope evolution and stellar mergers (e.g., Paczynski, 1976; Meyer & Meyer-Hofmeister, 1979). Recently, MacLeod et al. (2018b) studied the hydrodynamics of Roche-lobe overflow of a giant star and the associated mass loss from the vicinity of L2. For their binary mass ratio $q=0.3$, the L2 outflow should become unbound in ballistic calculations. However, MacLeod et al. (2018b) found that the outflow actually remains bound and forms a decretion disk because the initial tangential velocities of the stream are on average slower than corotation. This implies that the kinematics of the L2 mass loss stream is potentially sensitive to the actual initial conditions near L2.

In this paper, we extend the previous ballistic calculations and investigate the trajectories and asymptotic properties of mass loss launched from a range of spatial positions around L2 and with a varying direction and magnitude of the initial velocity. In Section 2, we describe our setup for integrating ballistic particle trajectories, which we supplement for selected cases with smoothed particle hydrodynamics. In Section 3, we identify regions of positive and negative stream energy, stream collisions with the binary and loops leading to self-collisions within the stream. In Section 4, we summarize our findings and discuss implications for catastrophic binary interactions.

2 Methods

Here we describe the geometry of the Roche potential and our numerical setup for integrating ballistic particle trajectories (Section 2.1) as well as our implementation of smoothed particle hydrodynamics, which we use for investigation of hydrodynamical shocks (Section 2.2).

2.1 Ballistic trajectories

We study the restricted three-body problem in the non-inertial reference frame corotating with the binary circular orbit with angular frequency $\omega$. The positions of the two point-mass stars are fixed and the test particles launched from around L2 do not affect the gravitational field around the binary. The binary components have masses $M_{1}$ and $M_{2}$ and are positioned on the $x$-axis with barycentre at $x=0$ and the less massive star $M_{2}$ with $x>0$. We restrict our analysis to trajectories in the binary orbital plane $(x,y)$, because the problem has mirror symmetry with respect to the orbital plane and the motions in the vertical direction are trivial. Using $a$ as the unit of length, $a\omega$ as the unit of velocity, the net potential $\phi$ in units of $a^{2}\omega^{2}$ for the particle motion in the frame corotating with the binary is

[TABLE]

where $\mu\equiv M_{2}/(M_{1}+M_{2})$. The L2 point is located at $x_{\rm L2}>0$ defined by solving

[TABLE]

The parameter $f={\omega^{\prime}}/\omega$ is introduced to define an equilibrium point analogous to L2 (zero net of gravitational and centrifugal accelerations) in situations when the initial particle angular frequency $\omega^{\prime}$ differs from $\omega$ (see Sepinsky et al., 2007, for more details on Lagrange points in eccentric, non-corotationg orbits). Unless explicitly stated otherwise, $f=1$ throughout the paper.

The particle trajectory is obtained by solving the equations of motion

[TABLE]

using the scipy.integrate.solve_ivp function from the SciPy library for Python. We applied the Dormand-Prince method of order 4/5 (RK45), which is an explicit method of the Runge-Kutta family with adaptive stepsize. During our computations, relative and absolute tolerances controlling the accuracy of $x$ and $y$ in units of $a$ were set to $10^{-12}$.

We explored the trajectories for a range of initial positions $\mathbf{r}_{0}=(x_{0},y_{0})$ and velocities $\mathbf{v}_{0}=(\dot{x}_{0},\dot{y}_{0})$. For the sake of brevity, we express $\mathbf{r}_{0}$ relative to the L2 point so that particles starting in corotation at L2 have $\mathbf{r}_{0}=(0,0)$ and $\mathbf{v}_{0}=(0,0)$. The integration was terminated when one of the following events occurred. (i) The particle crossed the inner equipotential line passing through the L2 point. We approximated the equipotential with a polygon and tested whether the trajectory entered inside every $0.25/\omega$ time units. (ii) The trajectory intersected with itself. We determined this by approximating the trajectory with line segments and checking whether a newly-calculated segment intersected any of the previous ones. (iii) The distance of the particle from system’s centre of mass has exceeded $200a$. At this distance, the tidal torque from the central binary is sufficiently weak to not significantly change the asymptotic outcome.

The energy $E$ and angular momentum $J$ per unit mass measured with respect to a stationary inertial observer can be computed from the position $(x,y)$ and the velocity $(\dot{x},\dot{y})$ of the particle in the corotating frame using the relations

[TABLE]

The specific energy computed from the above-mentioned expression has units $a^{2}\omega^{2}$, the specific angular momentum is given in units $a^{2}\omega$. We denote initial and final energies of particles by $E_{\rm i}$ and $E_{\rm f}$, respectively. Along the particle’s trajectory the value of Jacobi constant

[TABLE]

is conserved. We see from equations (4) and (5) that along the particle’s trajectory, the difference between its energy and angular momentum remains constant because $C=E-J$. Consequently, we focus on the behaviour of $E$, because $J$ can always be recovered from $C$.

2.2 Smoothed particle hydrodynamics

Self-intersections of ballistic particle trajectories lead to hydrodynamic shocks. To investigate the gas dynamics in these cases, we utilize our smoothed particle hydrodynamics (SPH) code (Pejcha et al., 2016a, b; Pejcha et al., 2017), which injects particles in a small volume around L2 and follows their dynamics in the external gravitational field of the binary star, approximated as two point masses. There is no gravitational interaction between the gas particles. Shocks are captured with standard artificial viscosity (Monaghan & Gingold, 1983; Balsara, 1995). Radiative diffusion, cooling and irradiation by the central binary are not included. We employ the realistic equation of state of Tomida et al. (2013).

To approximately take into account re-accretion of particles by the binary star, we introduce an absorbing inner boundary. The natural choice for the inner boundary would be the outer critical surface of the Roche potential going through the L2 point. However, it would be too time-consuming to evaluate whether each particle crossed inside this three-dimensional implicitly-defined surface at every timestep. Instead, we approximate this surface with a prolate ellipsoid that exactly touches L2 and the corresponding opposite point on the $x$ axis (Fig. 2). The other two axes are manually set to approximate the shape of the Roche equipotential near $M_{1}$. Although a proper way of implementing the inner boundary takes into account the properties of the stellar surface, we simply remove any particle that crosses inside the ellipsoid. As a result, we cannot reliably model interaction of the mass loss with the star and our results are thus only indicative of a more complete calculation.

3 Results

In this section, we describe the types of trajectories we identified in our work (Section 3.1) followed by investigation of how these trajectory types depend on the initial position (Section 3.2) and velocity (Section 3.3).

3.1 Types of particle trajectories

In Fig. 1, we show four types of particle trajectories identified in our work. The upper row shows previously known trajectory types, which were investigated in both ballistic and hydrodynamic approximations (e.g. Shu et al., 1979; Pejcha et al., 2016a, b). In panel (a), the particle stream becomes unbound and forms a thin equatorial outflow heated by relatively weak internal spiral shocks. In panel (b), the particles travel to a certain maximum distance, fall back to the binary, and eventually form a decretion disk or a more complicated object222For trajectories initially in corotation at L2, we find that the upper dividing line between bound and unbound trajectories is $q=0.792$, while Shu et al. (1979) found $q=0.78$. We were not able to trace the origin of this discrepancy.. We detect this outcome as self-intersecting trajectories which occur with $E<0$. The bottom row shows trajectories newly identified in this work, which occur for more general initial conditions. Panel (c) shows a particle colliding with the binary star. This results in hydrodynamical shocks and at least partial re-accretion of the gas stream. The nature of this interaction depends on the detailed structure of the binary star and it is thus difficult to assess the fractions of re-accreted and ejected material. Finally, panel (d) shows a trajectory which loops and intersects with itself near the binary. Most of the loops occur with $E<0$ and we thus do not distinguish them from decretion disk formation. The trajectory in panel (d) is an exception, where the loop occurs with $E>0$. This required an initial velocity vector in a particular direction and of relatively high magnitude. We do not find similar cases in the results shown below. We expect that trajectory loops lead to formation of hydrodynamical shocks in the stream and outside of the binary.

We used our SPH code to model the self-intersecting trajectory from panel (d) of Fig. 1. We show a representative snapshot of the final steady-state configuration in Fig. 2. After the simulation begins, the centre of the SPH particle stream follows the ballistic stream trajectory. As the stream turns back to the binary, around the third Lagrange point, some particles cross the inner boundary and are removed from the simulation without causing a hydrodynamic shock333When we made the inner absorbing boundary smaller, these particles formed a narrow accretion funnel similar to what is seen in accreting circumbinary disks (e.g. Artymowicz & Lubow, 1996; Muñoz & Lai, 2016).. As the stream centre gets to the turning point of the ballistic particle loop, SPH particles collide and the inner rim of the stream becomes shocked. This is visible as a dense and hot arc immediately below and to the right of the ballistic particle loop in Fig. 2. However, part of the shocked rim is now on the leading side of the binary, loses energy and angular momentum from the tidal torques and soon disappears below the inner boundary. However, some of the shock-heated gas also expands away from the orbital plane and fills the binary surroundings with a hot medium. Some fraction of the gas becomes unbound and flows away. To summarize, self-intersecting trajectories deposit gas with low angular momentum and/or high temperature in the vicinity of the binary star. This is somewhat similar to decretion disk formation although, for a decretion disk and moderate $q$, the shock is farther away from the binary.

3.2 Dependence on initial position

In Fig. 3, we show the trajectory type and final energy as a function of the starting position $\mathbf{r}_{0}=(\Delta x,\Delta y)$, relative to L2, and the binary mass ratio. All trajectories were initiated with $\mathbf{v}_{0}=(0,0)$. Trajectories launched very close to the binary surface end up eventually colliding with the binary (grey hatched areas). The hatched area on the $\Delta y>0$ side is larger than for $\Delta y<0$, because trajectories need to pass around the physical barrier near L2, where the star extends to largest $x$. Trajectories launched further away from the binary can escape to infinity ($E_{\rm f}>0$) or return to the binary and form a decretion disk ($E_{\rm f}<0$). We see in Fig. 3 that negative energies always occupy a lobe-like region, which is either confined to $\Delta y>0$, for intermediate $q$, or extends to L2 and $\Delta y<0$, for small and high $q$. Interestingly, there is always a narrow region of $E_{\rm f}>0$ located at $\Delta y>0$ and adjacent to initial conditions that lead to collisions with the binary.

To better understand these features, we show in Fig. 4 the initial energy $E_{\rm i}$ (left panel) and the energy gain $\Delta E=E_{\rm f}-E_{\rm i}$ (right panel) as a function of positional offset. We also show excerpts of three trajectories with different starting points. Furthermore, Fig. 1d shows time derivative of $E$ as a function of position for particles initially in corotation with the binary. In Fig. 4, $E_{\rm i}$ is higher away from the binary, because both potential and kinetic energy increase. The potential energy increases, because the particle is farther away from the two point masses. The kinetic energy increases, because higher velocities are needed to ensure corotation with the binary. The initial energy is symmetric with respect to the $x$ axis. The energy gain is set by the efficiency of tidal torquing, which is illustrated in Fig. 1d. In the vicinity of L2 tidal torques increase $E$ for particles lagging behind ($\Delta y<0$) and decrease $E$ for particles leading ahead ($\Delta y>0$). We also see that there is always a positive energy gain adjacent to the region of collisions with the binary surface, even for $\Delta y>0$. As is illustrated by the middle particle trajectory (Fig. 4, green line), trajectories launched from this region pass around the obstacle of the binary surface as they are accelerated by the Coriolis force and Roche potential gradient. After crossing to $\Delta y<0$, the trajectory is still very close to the binary and experiences strong tidal torqueing (Fig. 1d). The energy gain in this region more than compensates for $E_{\rm i}<0$ and for energy losses while the particle was leading ahead of the binary.

3.3 Dependence on initial velocity

We now investigate changes to the picture outlined in Section 3.2 when the trajectory starts with non-zero velocity in the corotating frame. It is instructive to first investigate purely radial initial velocities, for which we show the results in Fig. 5. As the velocity in the radial direction increases, the trajectories reach farther away from the binary. As a result, the area of initial conditions leading to collisions with the binary surface shrinks. At the same time, for small initial velocities, $|\mathbf{v}_{0}|\equiv v_{0}\lesssim 0.3$, the area with $E_{\rm f}<0$ increases. This occurs because the efficiency of tidal torqueing sharply drops with the distance from the binary (Fig. 1d) and the trajectories do not gain enough energy to become unbound despite higher $E_{\rm i}$. For even higher $v_{0}$, the initial kinetic energy is so large that the area of $E_{\rm f}>0$ increases again (bottom row of Fig. 5).

An alternative way to view our results is to investigate the range of $q$ leading to $E_{\rm f}>0$. We studied trajectories starting from $\mathbf{r}_{0}=(0,0)$ on a fine grid of $q$ and radial velocity magnitude $v_{0}$. As $v_{0}$ increases the range of $q$ leading to $E_{\rm f}>0$ slightly increases and peaks around $v_{0}\approx 0.04$. For higher $v_{0}$, the range shrinks and achieves a minimum around $v_{0}\approx 0.3$, where $E_{\rm f}>0$ is achieved only for $0.15\lesssim q\lesssim 0.27$. For $v_{0}\gtrsim 0.7$ all $q$ lead to $E_{\rm f}>0$.

In Fig. 6, we show $E_{\rm f}$ for trajectories initially on a circular orbit but with $\omega^{\prime}=f\omega$, $f>1$. The results are qualitatively similar to the radial velocity: for $f\lesssim 1.1$, it is easier to evade collision with the binary surface, especially for the initial conditions of the leading side of the binary motion, but often at the price of less efficient tidal torqueing, because the trajectories move further away from the binary. As a result, the contour of $E_{\rm f}=0$ gets closer to L2. For $f\gtrsim 1.1$, the high initial velocity results in an increase of the area with $E_{\rm f}>0$. We also note a thin line of initial positions with $E_{\rm f}<0$ passing through the hatched area at $\Delta y>0$. These trajectories looped and self-intersected before they collided with the binary (Fig. 1d).

In Fig. 7, we show the results for $f<1$. As expected, the area of initial conditions leading to collisions with the binary is significantly larger owing to lower kinetic energy. However, if we relax the assumption of corotation, the Roche potential is modified and so are the positions of Lagrange points (Sepinsky et al., 2007). Specifically, for $f<1$, $x_{\rm L2}$ increases as we show in Fig. 7. Binary stars with surface layers rotating slower than the orbital frequency would lose mass from the vicinity of the modified L2 point. In such cases, there still exist initial positions near the modified L2 which eventually achieve $E_{\rm f}>0$, at least for moderate $q$. However, trajectory self-intersections (Fig. 1d) also become prominent especially for smaller $f$. Furthermore, in generating Fig. 7, we still checked for collisions with the binary surface assuming the original corotating equipotential. It is likely that taking into account the modified equipotential would increase the number of colliding trajectories.

Finally, we also investigated initial velocity vectors pointing to an arbitrary direction. In some cases, the final outcome of the trajectory sensitively depends on the initial conditions. This suggests a possibility of deterministic chaos. Typically, these features occur for trajectories closely approaching some of the other Lagrange points. Around these stationary points, the direction of the Roche potential gradient varies significantly and the trajectories effectively scatter to a wide range of directions. However, this behaviour would likely result in hydrodynamic shocks which remain to be explored in more detail in the future.

4 Discussion and conclusions

We have performed detailed investigations of asymptotic behaviour of ballistic trajectories launched from the vicinity of the L2 point. In addition to previously known outcomes of unbound outflow and fallback forming a decretion disk, we identified a relatively narrow space of initial conditions where the ballistic trajectories self-intersect close to the binary star (Fig. 1d). This leads to formation of hydrodynamic shocks, as we illustrate with smoothed particle hydrodynamics (Fig. 2).

We studied how the final outcome depends on the initial position (Fig. 3) and velocity vector (Figs. 5, 6 and 7). We find that all types of outcomes are present in a relatively narrow extent of physical offsets from L2. For example, there generally exists a thin arc of positive final energies on the leading side of the binary, where tidal torqueing is surprisingly efficient (Fig. 4, right panel, and Fig. 1d). Trajectories initiated with a small offset from this arc either collide with the binary or do not achieve positive final energy. Furthermore, even for trajectories initiated about $20$% slower than corotation, there still exist regions where the trajectories end up unbound.

Because a realistic stream leaving the binary from L2 would start from a range of initial positions and velocities, our results imply that interpreting hydrodynamical simulations of binary mass loss with ballistic trajectories needs to take into account the sensitivity to initial conditions. This agrees with findings of MacLeod et al. (2018b), who also point out that hydrodynamical interactions within the stream, which we cannot study with ballistic trajectories, contribute to the final outcome. However, our Fig. 3 suggests that very small shifts in the initial $\Delta x$ of ballistic trajectories with $\Delta y>0$, such as those used in Fig. 5 of MacLeod et al. (2018b), could potentially result either in collisions with the binary or even positive asymptotic energies. Furthermore, the mass loss from the vicinity of L2 is affected by binary surface perturbations owing to convective motions, meridional and zonal circulations, and magnetic fields, all of which intensify as the mass-loss rate runs away. These physical processes are not resolved in most of the binary hydrodynamical simulations.

Despite these uncertainties, our results already offer clues how to interpret observations of suspected L2 mass loss. In their study of luminous red nova V1309 Sco, Pejcha et al. (2017) found that the initial temperature of the L2 stream needs to dramatically increase in the last about 75 days to explain the observed upward curvature of luminosity. Their results imply that $f\approx 0.95$ at 75 days before the merger. Our results in Figure 7 show that even for $f=0.95$ there is an increased fraction of initial conditions that lead to collision with the binary surface and that $E_{\rm f}$ of unbound particles is lower than for $f=1$. Although we do not know the distribution of initial positions of the stream, it is likely that the fraction of L2 stream remaining near the binary increases with decreasing $f$. This material will shock-heat to temperatures similar to the kinetic temperature of the binary, which could naturally explain the increase of initial temperature of the L2 stream.

Acknowledgements

We thank our referee, Prof. Christopher Tout, for comments and suggestions, which improved the paper. Our research is supported by Horizon 2020 ERC Starting Grant “Cat-In-hAT” (grant agreement #803158), award PRIMUS/SCI/17 from Charles University and INTER-EXCELLENCE grant LTAUSA18093 from the Czech Ministry of Education, Youth, and Sports.