Runaway Coalescence at the Onset of Common Envelope Episodes

[ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA School of Natural Sciences, Institute for Advanced Study, Princeton, NJ, 08540, USA [ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA [ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

Luminous red novae transients, presumably from stellar coalescence, exhibit long-term precursor emission over hundreds of binary orbits leading to impulsive outbursts, with durations similar to a single orbital period. In an effort to understand these signatures, we present and analyze a hydrodynamic model of unstable mass transfer from a giant-star donor onto a more-compact accretor in a binary system. Our simulation begins with mass transfer at the Roche limit separation and traces a phase of runaway decay leading up to the plunge of the accretor within the envelope of the donor. We characterize the fluxes of mass and angular momentum through the system and show that the orbital evolution can be reconstructed from measurements of these quantities. The morphology of outflow from the binary changes significantly as the binary orbit tightens. At wide separations, a thin stream of relatively high-entropy gas trails from the outer Lagrange points. As the orbit tightens, the orbital motion desynchronizes from the donor’s rotation, and low-entropy ejecta trace a broad fan of largely-ballistic trajectories. An order-of-magnitude increase in mass ejection rate accompanies the plunge of the accretor with the envelope of the donor. We argue that this transition marks the precursor-to-outburst transition observed in stellar coalescence transients.

binaries: close, methods: numerical, hydrodynamics

0000-0002-1417-8024]Morgan MacLeod \altaffiliationNASA Einstein Fellow

0000-0002-0509-9113]Eve C. Ostriker

0000-0001-5603-1832]James M. Stone

1 Introduction

Common envelope episodes occur in binary systems when a giant star engulfs its companion within its gaseous envelope (Paczynski, 1976). Drag forces tighten the orbit of the stellar cores during this phase, depositing orbital energy and momentum into the envelope gas (e.g. Iben & Livio, 1993; Taam & Sandquist, 2000; Ivanova et al., 2013b). The outcome of orbital tightening during this phase can be a more-compact binary and an ejected envelope, or a merged, single object with some retained envelope material (Taam et al., 1978).

Such a phase of orbital tightening is thought to be crucial in the transformation of binaries into compact systems that can merge within the age of the universe under the influence of gravitational radiation (e.g. Belczynski et al., 2002; Kalogera et al., 2007; Belczynski et al., 2008). As merging compact binaries are discovered in increasing numbers by gravitational wave interferometers like the LIGO-VIRGO network (Abbott et al., 2016a, b, 2017a, 2017c, 2017b) there is an increased significance to tracing the evolutionary history that leads to the formation of these objects (see Postnov & Yungelson, 2014, for a review). Despite the centrality of the common envelope phase to the formation of compact binaries (Belczynski et al., 2002; Toonen & Nelemans, 2013; Tauris et al., 2017), the details of these encounters, and which binaries map to which outcomes, have remained major sources of theoretical uncertainty for decades (Ivanova et al., 2013b).

New, empirical, evidence is now offering an avenue to better understand mass ejection in common-envelope-like encounters. It has become clear that the class of transients known as luminous red novae represent emission from cooling, expanding material ejected in either common envelope or stellar-merger encounters (Ivanova et al., 2013a). A key object in identifying this connection was the galactic transient V1309 Sco (Mason et al., 2010), which was identified by Tylenda et al. (2011) in pre-outburst data as an eclipsing binary with a decreasing orbital period. Several extragalactic systems, including M31 LRN 2015 (Kurtenkov et al., 2015; Williams et al., 2015), M101OT 2015-1(Blagorodnova et al., 2017), and NGC 4490 OT (Smith et al., 2016), have added extra evidence and posed new questions.

These transients exhibit outbursts defined by rapid rise to peak in the optical followed by fading and reddening emission. On timescales of months, dust formation leads to optical extinction and extensive reprocessing into the infrared (e.g. Nicholls et al., 2013; Tylenda & Kamiński, 2016). Another shared feature among recent red novae transients is precursor emission in the months to years leading to outbursts themselves (e.g. Tylenda et al., 2011; Blagorodnova et al., 2017), strongly suggesting binaries interacting in circular orbits rather than single, impulsive passages (e.g. Smith, 2011). The phase of rise and peak is particularly well sampled in the M31 LRN 2015 transient, in which the rise of the light curve from transient discovery to peak is similar to one orbital period of a test mass on the surface of the observed giant star progenitor (MacLeod et al., 2017).

The occurrence of impulsive transients might seem in tension with hints at circular binary orbits: how can a binary, which is stable over long-term, stellar-evolution timescales, produce a transient with timescale similar to the orbital period? Perhaps not surprisingly, the optical photosphere traces rapidly-expanding ejecta in these transients, rather than the enveloped binary (e.g. Pejcha, 2014; MacLeod et al., 2017). Rather than expanding freely, internal shocks may play a significant role in the ejecta’s thermodynamics (Metzger & Pejcha, 2017). To decipher how the underlying binary dynamics emerge into various components of the light curve, we need to trace these ejecta with the maximum fidelity possible in our simulation counterparts to these encounters.

In an effort to better understand the evidence about mass ejection in common envelope episodes offered by luminous red novae transients, we develop and present a new simulation method to study the lead-in to common envelope events. We study the phase from the onset of mass transfer from a more-massive, evolved primary star onto a lower mass accretor. Because of the mass ratio, mass transfer shrinks the binary separation, leading to an unstable increase in mass-loss rate from the donor. The secondary star is modeled only gravitationally, in reality, it may be either an unevolved main-sequence star or a compact remnant. We validate our approach through a series of careful convergence studies. We trace this runaway mass transfer and the resulting dynamics and forces that lead to the engulfment of the secondary star within the envelope of the giant-star primary and the onset of a common envelope phase.

We begin with some analytic background for the evolution of binary orbits under the influence of mass exchange in Section 2. We describe our numerical simulation method in Section 3; a number of tests are presented in Appendix A. Using this method, we present and analyze a simulation of the coalescence of a binary with mass ratio 1:0.3 that begins with mass transfer from a more-massive, giant-star primary to a less-massive and more-compact accretor in Section 4. In Section 5, we argue that the features of runaway orbital inspiral at the onset of common envelope episodes, as seen in our simulations, explain the presence of both precursor emission and impulsive outbursts in luminous red novae transients. In section 6, we summarize and conclude.

2 Analytic Context for Binary Mass Transfer and Orbital Evolution

In a binary system that evolves into contact, mass transfer begins from a donor star onto the accretor. The accreting object may accept this mass, in which case the mass transfer is described as conservative in that the system conserves total mass and angular momentum. In this case, there is a simple relationship that describes the rate of change of orbital separation, , in terms of the mass transfer rate and binary mass ratio (Huang, 1956).

If mass loss and the associated changes in orbit and donor-star structure cause the donor to become less contained within its Roche lobe, the mass-loss rate intensifies. Paczyński & Sienkiewicz (1972) estimate the mass loss rate of a polytropic donor star of index to be


where and are the mass and radius of the mass-losing donor star, is the radius of the Roche lobe (Eggleton, 1983), and is the binary orbital period (see Shore et al., 1994, for a derivation). Therefore for , the exponent is 3. Later, Edwards & Pringle (1987) validated this proportionality across a small range of mass-loss rate in two-dimensional hydrodynamic simulations. These scalings indicate that when orbital evolution leads to increased Roche lobe overflow, the mass transfer rate is subject to a runaway to ever-increasing rates.

In these cases of unstable mass transfer, the mass transfer rate is very likely to increase past the rate at which material may be accepted by the accreting object. This implies that some mass, and with it, some angular momentum, is carried away from the binary. In the extreme limit the majority of mass transferred is lost from the system, and mass transfer is described as non-conservative. The orbital evolution equation (Huang, 1963) for the case of fully non-conservative mass loss from a binary is


where is the mass of the donor-star primary, is the mass of the accretor, and parameterizes the specific angular momentum of the lost mass as


which is the ratio of the specific angular momentum of lost material, , to the specific angular momentum of the binary, , which is its total angular momentum divided by its mass.111A useful summary is presented in O. Pols’ binary-evolution notes, Chapter 7.

Studying equation (2) reveals that the evolution of the orbit under non-conservative mass transfer depends critically on the specific angular momentum of the material lost, through the parameter . A critical value separates solutions that imply orbital expansion from those that imply orbital tightening, . Because , situations with lead to widening orbits, while those with lead to shrinking orbits.

To illustrate the possibilities, consider a scenario where a more-massive donor transfers material onto a less-massive accretor (). We might imagine that material lost from the donor can leave the binary system with (Huang, 1963):

  1. the specific angular momentum of the donor star, as if being blown off in a spherical wind. In this case, .

  2. the specific angular momentum of the accretor star, as if being pulled from the donor onto the accretor and re-emitted. In this case, .

  3. the specific angular momentum of the outer Lagrange point near the secondary, , which as a saddle point in the potential might allow outflow in the plane of the orbit. In this case, (Pribulla, 1998).

These scenarios lead to major differences in orbital evolution. In case (1), and the binary widens its separation. In case (2), , and the binary orbit tightens. In case (3), again , but with larger magnitude because the point represents a larger lever arm from the system center of mass.

Our knowledge of and intuition for the geometry of mass transferring binaries (e.g. Lubow & Shu, 1975) makes case (1) appear unlikely – mass is typically taken from the donor to the vicinity of the accretor, which implies an exchange of angular momentum (Huang, 1963). Qualitatively, then, the prediction of these analytic models (cases 2 and 3) for unstable mass transfer is orbital tightening driven by increasingly-rapid mass transfer (e.g. Webbink, 1977). By the time the orbital separation shrinks to be similar to the size of the donor, , we enter a common envelope phase – during which it is difficult to discuss the stars as separate entities.

While these parameterized models provide a framework to understand binary evolution toward coalescence, the do not offer key quantitative details. In particular, while the scaling in (1) comes from dimensional analysis, the coefficient is unknown. Similarly, the precise magnitude of the order-unity quantity determines the rate and even direction of evolution, but is also unknown. Motivated by this uncertainty, we use hydrodynamical simulations of mass exchange in a binary system to study the runaway coalescence of merging binaries in the remainder of this paper.

3 Simulation Models and Method

Simulating the common envelope phase of binary systems has been a decades-long effort in astrophysics, from 1D approximations (Meyer & Meyer-Hofmeister, 1979; Taam et al., 1978), to the “Double-core” series of multi-dimensional efforts which revealed the lack of spherical symmetry in the problem (Bodenheimer & Taam, 1984; Taam & Bodenheimer, 1989, 1991; Terman et al., 1994; Taam et al., 1994; Terman et al., 1995; Yorke et al., 1995; Terman & Taam, 1996; Sandquist et al., 1998), to recent three dimensional hydrodynamic (e.g. Ricker & Taam, 2008, 2012; Nandez et al., 2014, 2015; Iaconi et al., 2017; Ohlmann et al., 2016a) and magnetohydrodynamic (Ohlmann et al., 2016b) calculations with both Eulerian and Lagrangian methods (see a review of methods in De Marco & Izzard, 2017, section 4). The majority of this work has focused on the outcome of the phase and the final separation of the tightened orbit of the stellar cores (see Ivanova et al., 2013b, for a review).

Here we focus on the lead-in to the common envelope phase itself. Early interaction determines the observable properties of common envelope ejecta as well as the initial conditions for the bulk of the common envelope interaction. Previous authors of have either neglected the early interaction entirely in hydrodynamic calculations, starting the calculation with the stars already in contact (e.g. Ohlmann et al., 2016a), or focused on other characteristics of the system in their analysis. There are several notable exceptions. Iaconi et al. (2017) examined the effect of initial separation on simulation progression and final binary separation. Nandez et al. (2014) tracked ejecta in simulated systems thought to be similar to V1309 Sco, but using a Lagrangian method that, in tracing mass rather than volume, lacks the resolution to follow the low-mass, low-density early interaction. Relatedly, Galaviz et al. (2017) attempted to synthesize observables from hydrodynamic simulations but found that poor resolution and thermodynamics near the photosphere considerably limited current capabilities. Shiber & Soker (2017) and Shiber et al. (2017) study the interaction of jets with tenuous envelope material during an early, grazing, phase of binary interaction, but these calculations again begin with the binary components in contact and, further, neglect their gravitational interaction.

Other work of methodological relevance includes simulations of mass transferring and merging binaries. Rasio & Shapiro (1995) discuss initial data and the hydrodynamics of merger of polytropic binaries. Motl et al. (2002) and later D’Souza et al. (2006), and Motl et al. (2007), developed techniques for the simulation of mass-transferring binaries on Eulerian meshes, and emphasized the importance of frame-choice in conserving key quantities (see also related work by Katz et al., 2016). Motl et al. (2017) have recently reviewed this work, and compare to particle-based methods. Dan et al. (2011), Lombardi et al. (2011), and Pakmor et al. (2012) emphasize the importance of relaxed initial conditions (in the binary potential) for the subsequent evolution using smoothed-particle-hydrodynamic approaches.

The approach presented here does not represent a dramatic departure from previous efforts, but we do make a series of numerical choices such that our method is particularly suited to studying this phase of lead-in to engulfment that has become so interesting in light of recently-observed transient outbursts. Of particular importance, we adopt an Eulerian method to solve the equations of gas dynamics, which allows us to follow low-density flow in (and outside) stellar envelopes. Because the distribution of this material contributes to torques in the system, following this low mass of material is critical to resolving the early interaction. Additionally, we simulate in spherical polar coordinates surrounding the giant-star donor, which allows us to preserve the hydrostatic structure of the donor envelope to particularly high fidelity.

We describe our method in the sections that follow and validate it through a series of convergence studies in Appendix A.

3.1 Simulated System and Reference Frame

We consider a binary star system composed of a giant star donor (primary, object 1) and its accreting companion (secondary, object 2). A cartoon of the simulated system is shown in Figure 1. Our simulations are conducted in spherical polar coordinates in three dimensions.

Here we consider cases where the companion is a main sequence star or compact object, and is not resolved on the scale of our simulations. We therefore represent this secondary object only gravitationally with a point mass. Similarly, the core of the giant star is excised from the computational domain and replaced with a point mass. The simulated system therefore contains two point masses as well as the gaseous envelope of the giant star.

The total mass of the primary star is , where is the mass of point mass 1, representing the core of the giant star and is the mass of the gaseous envelope. The primary star has total radius . The secondary has mass , it is made up entirely of a point mass, so . The total system mass , which generates an orbital frequency of , where is the binary semi-major axis. Equivalently, the orbital period is given by .

We adopt an orbiting, but non-rotating frame of reference. The coordinate origin centered on (and orbiting with) the giant star donor. The coordinate frame of the calculation is therefore non-inertial, in that it accelerates with the orbital motion. After testing both rotating and non-rotating coordinate frames, we chose the non-rotating frame for the calculations presented here because we found that it more-consistently preserved hydrostatic equilibrium of stars of arbitrary rotation. In a rotating frame, when the rotation of the star does not match that of the mesh, hydrostatic equilibrium implies a numerical balance of the centrifugal and Coriolis accelerations in addition to those from gravity and pressure gradients.

Figure 1: Coordinate system and domain of numerical simulations. Two point masses, and interact with the simulated gas. The orbiting, non-rotating simulation frame orbits with and the mesh is allocated in spherical polar coordinates. The secondary object, labeled orbits in the () plane. The simulation domain extends from to in azimuthal angle , from a minimum radius to in radial coordinate and does not necessarily extend to the poles, with the polar angle, , ranging from 0 to .

3.2 Equations Solved

Here we outline the equations solved in our hydrodynamic evolution, which are the basic expressions of inviscid gas dynamics with source terms related to the binary system and choice of reference frame. Our simulation is developed within the Athena++222Stone et. al. (in preparation), publicly available at: code, which is an Eulerian (magneto)hydrodynamic code descended from Athena (Stone et al., 2008).

3.2.1 Inviscid Gas Dynamics

The coupled equations of inviscid gas dynamics that we solve are


expressing mass continuity, the evolution of gas momenta, and the evolution of gas energies. In the above expressions, is the mass density, is the momentum density, and is the total energy density and is the internal energy density. Additionally, is the pressure, is the identity tensor, and is an external acceleration (source term). The equations above are closed by an ideal gas equation of state, , in which is the gas adiabatic index or ratio of specific heats.

3.2.2 Source Terms

The source term on the right-hand side of equations (4), , represents the forces associated with the binary, as well as fictitious accelerations associated with our choice to perform the calculation in a non-inertial frame of reference.

The orbital plane in our simulations is the (or ) plane, and the angular momentum vector of the orbit points in the -direction (). The terms that make up are


representing the acceleration from the gravity of the primary core particle, , the secondary, , the acceleration from the gradient of the gas self-gravitational potential, and the orbiting-frame correcting acceleration of the primary measured in the inertial frame, respectively.

Within the code, we derive these accelerations in cartesian coordinates, then convert them to spherical polar. With the accelerations defined, we convert back to the local spherical polar coordinates using


Finally, these source accelerations enter the gas momentum and energy equations in (4).

The acceleration from is


where is the vector displacement from the coordinate origin – where the donor-core particle is located. The acceleration from the secondary particle is


where is the vector position of the secondary particle. The approximate version is applied numerically and includes a softening spline at small distances from the particle,


where (Hernquist & Katz, 1989, equation A2). As compared, for example, to a Plummer softening, this spline softening has the advantages of converging much more quickly to the Newtonian force when and (more importantly) being exact outside two softening lengths, .

We treat the self-gravitational interaction of the gas only approximately in our simulation. We adopt


where is the enclosed mass of the original primary-star’s (initial) spherical profile at distance . That is, we apply an approximate, “static self-gravity”, to the mass of envelope gas. This choice most-affects our solution when significant mass is lost from the donor star, or under dramatic departures from spherical symmetry. We discuss the impacts of this choice for total angular momentum conservation further at the end of Appendix A.

The acceleration of the donor star relative to the inertial frame is


where the first and second terms represent the acceleration due to the accretor (particle 2) and the acceleration due to the gas is


We compute the integral above once per hydrodynamic time step. We apply the inverse of this acceleration to the gas in equation (5) to correct for the non-inertial frame of the simulation.

3.3 Secondary Point-Mass Particle

As described earlier, to simulate the effect of a secondary in the binary system we add a perturbing particle with mass to the system. We integrate the motion of this particle along an arbitrary orbit, which responds to the primary particle, , and gas-distribution accelerations. We assume that the mass of the secondary is constant in time – which implies that it neither gains nor loses mass due to interaction with the surrounding gas. This is likely an oversimplification, but without much higher spatial resolution it is difficult to even determine the qualitative outcome of accretion versus mass loss.

3.3.1 Initial Condition

The secondary is initialized based on a choice of semi-major axis and eccentricity. It is initialized in the -direction, at apocenter, such that the initial x-coordinate (and separation) is , the initial coordinates are set to zero.

The initial velocity is set along the direction (such that the plane becomes the orbital plane) and with a positive velocity such that the angular momentum is oriented in the -direction. Together, we have,


Finally, on the first call to the particle integrator, we kick the particle velocities backward a half step to achieve staggered position and velocity for the leapfrog integrator described below.

3.3.2 Integration

We integrate the motion of the secondary’s particle in response to the gravitational force from the primary-star particle and from the gas distribution. Because of our choice of a non-inertial frame, the primary star particle is fixed at the origin and the equation of motion for the secondary particle includes extra terms.

The particle step is currently synchronized with the hydrodynamic time step, and is controlled by an integer number of substeps per hydrodynamic step. We use a leapfrog algorithm, with the velocities offset a half step behind the positions in time, in order to be centered across the time step. An evolution step therefore consists a kick and a drift. The kick is

and the drift is,

which advances the particle positions to the next step.

3.3.3 Accelerations on Particle 2

The acceleration acting on the secondary particle in the orbiting and rotating frame is


where is given by (11). The acceleration of point mass 2 from point mass 1 is


The contribution from the gas distribution within the simulation domain to the secondary particle’s acceleration is


which is summed once per time step. If substeps are used in integrating the particle motions the gas contributions are not re-integrated each substep, only each hydrodynamic time step.

3.4 Domain and Mesh

As described above, we solve the equations of gas dynamics on a spherical polar mesh centered on the giant-star core particle . We employ zones which are uniformly spaced in and coordinates, but logarithmically spaced in so that zone shapes remain roughly cubic across a wide range of radial scale.

For parallelization, the mesh is broken into meshblocks. In our fiducial cases, we employ 12 () x 8 () x 16 () mesh blocks of zones each at our base level. We additionally employ two nested levels of static mesh refinement opening about the binary orbital plane. Mesh refinement subdivides each zone by a factor of two in each direction, creating 8 smaller zones. We achieve these nested levels with three refinement patches. One level above the base, we refine within for one meshblock above and below the orbital plane, at all azimuthal angles, . We also refine the region and at all to one level above the base level. Finally, we refine the region and to two levels above the base.

The spherical polar coordinate system has singularities at and . Zones in these regions can become very elongated, and narrow in the direction, unless the resolution is very coarse. This carries accompanying, and often severe, restrictions on the simulation time step. One way to deal with this problem is to have fewer zones in the direction near the poles than at the equator. Mesh blocks in Athena++ are fixed in dimension, but we create effectively larger zones by averaging conserved quantities across neighboring zones near the poles, and modifying the time step criteria accordingly. In the zones nearest the pole, we average all of the zones within a mesh block, creating one effective zone per meshblock in the direction. As we move further from the pole the number of zones averaged over drops linearly with , rounding to the next power of two, such that the zone shape is approximately constant.

3.5 Units, Initial, and Boundary Conditions

3.5.1 Units

Quantity Code Example System A Example System B Example System C
grav. constant
time 0.583 d 1.91 d 14.2 d
velocity 138 km s 148 km s 125 km s
density  g cm  g cm  g cm
pressure  dyn  dyn  dyn
mass transfer rate 627  yr 766  yr 464  yr
angular momentum  g cm s  g cm s  g cm s
torque  g cm s  g cm s  g cm s
secondary mass
initial orbital period 16.3 9.49 d 31.1 d 231 d
initial orbital frequency 0.386 0.662 d 0.202 d  d

Note. – Together the unit mass, length, and gravitational constant define the dimensionless unit system we adopt in our simulations, labeled “Code” above. Derived unit values for various quantities are shown in the middle section. Finally, the lower section lists some properties of the binary system simulated – the mass of the secondary and the initial orbital period, , and frequency, , at the Roche-limit separation, . We convert our code units above to three example systems. Example System A is a , donor with a companion, which has an orbital period at the Roche limit of 9.49 days. Example System B is a , donor with a companion, which has an orbital period at the Roche limit of 31.1 days (these are roughly the inferred donor-star properties for M31 LRN 2015, though the secondary-star mass is unkown (MacLeod et al., 2017)). Example System C is a , donor with a companion, which has an orbital period at the Roche limit of 231 days (these are the approximate donor-star properties for M101 OT2015-1, where again the secondary-star properties are unconstrained (Blagorodnova et al., 2017)).

Table 1: Description of Units

The simulation is performed in units in which . This implies that the unit velocity is and the unit time is the dynamical timescale of the primary star . Quantities can be rescaled based on these units to physical systems. In Table 1 we list the units of primary and derived quantities in the dimensionless simulation, and illustrate their rescaling to several physical systems.

3.5.2 Initial Condition

The gas initial condition is a polytropic envelope in hydrostatic equilibrium. The equations governing its structure are


where and is the structural polytropic index of the envelope. In general, one must iterate to find a solution to these equations with a given total mass and radius. Outside the stellar radius, we continue in hydrostatic equilibrium, but create a higher sound speed, low density “corona” with constant sound speed of using


We compute and tabulate a profile separately and interpolate onto the mesh at runtime. We allocate the giant-star core point mass, , based on this profile and the mesh inner radius.

The initial and components of the gas velocity are set to zero. In some cases we want the envelope to rotate about the axis of the orbital angular momentum. We specify a rotational frequency, . Then,

within the envelope . Outside the envelope ,

implying constant specific angular momentum (of the surface of the star at that polar angle to the material).

3.5.3 Relaxation Scheme

Once the initial conditions are allocated on the mesh, we relax toward the fully dynamical portion of the simulation in two ways, by damping gas velocities at early times to allow the structure to settle into a numerically-hydrostatic equilibrium, and by turning on the force of the secondary object slowly following this relaxation period.

In the case of a rotating envelope, the hydrostatic structure is different (as a function of polar angle ) from the non-rotating initial profile. Injecting non-equilibrium conditions quickly leads to shock formation and the destruction of the envelope without some sort of damping and relaxation scheme.

We adopt a damping source term, which is applied (included in of (5)) for a time at the beginning of the simulation,


where is a damping timescale, and the damping is applied in the and directions only. The damping timescale is set as for and


so that the damping time increases to 100 by .

During the secondary mass is not active. We then turn on the source term from the secondary mass slowly after the relaxation of the initial hydrostatic profile. For times , we apply a fractional factor to the source terms associated with the gravity of the secondary. This factor increases linearly from 0 at to 1 at . By turning on this source term over significantly greater than the envelope dynamical time, we reduce the effects of any initial transient due to impulsively applying the gravitational force of the secondary to the envelope. We typically adopt such that the simulation is fully active after 30 donor-star dynamical times.

3.5.4 Boundaries

In the -direction, we impose a “diode” boundary condition on the outer boundary that is an outflow-only boundary that does not allow inflow onto the grid . The inner boundary condition is reflecting. In the -direction, we apply “polar” boundary conditions which wrap the mesh appropriately to allow free passage of flow through the pole. In the -direction, we apply periodic boundaries so that flow wraps fully around the to domain.

3.6 Diagnostics

We record several diagnostic quantities at runtime to facilitate analysis of the binary system. We record the position velocity of the secondary particle. We also record the gas accelerations on each of the particles, equations (12) and (17). We calculate the (cartesian) position and velocity of the particle and gas center of mass in the simulation frame


where all quantities are in the simulation frame (which implies ). These quantities allow correction of simulation quantities to the inertial frame in analysis.

With the center of mass position and velocity computed, we also compute quantities related to the exchange of angular momentum between the simulated components of gas and particles. The particles’ net angular momentum relative to the center of mass is


where all of the initial system orbital angular is in the component. The net gas angular momentum is


Finally, gas can leave the simulation domain through the outer boundary, carrying with it angular momentum. We record the instantaneous rate of angular momentum loss through the cell faces at the outer boundary,


where is the mass flux through the surface element of area dA. The total system angular momentum at time is, therefore,


The total angular momentum is a conserved quantity in a closed physical system. There are, however, reasons this may not be precisely the case numerically. First, we do not treat the backreaction of the inner- boundary condition on particle-1’s motion (loss of angular momentum from the outer- boundary of the domain is recorded through (26)). Second, we treat gas-particle gravitational interaction directly, but gas self-gravity only approximately with a static potential. Finally, our equations of gas dynamics are written in a form which explicitly conserves linear momenta along the coordinate directions, but not angular momenta, so some departure from perfect conservation is to be expected. The measurement of is therefore a useful diagnostic of the degree to which our model conserves system angular momentum as it is exchanged between the particle and gas components. We typically find that angular momentum is conserved to the percent level as particle and gas angular momenta exchange by roughly a factor of two (see Appendix A for further tests and discussion).

4 Results: Model Merging System

In this section, we follow the properties of gas flow and orbital evolution in a system that is initiated with mass transfer from a giant-star donor to a less-massive accretor. The orbit first evolves slowly, but enters into a phase of runaway coalescence. We introduce the simulated system in Section 4.1 and the overall evolution in 4.2. We follow the exchange of angular momentum through gravitational torques in Section 4.3. Finally, we discuss the evolution of outflow morphology as the orbital separation shrinks in Section 4.4.

4.1 Simulation Parameters and Initial Conditions

Figure 2: 1D initial condition for a polytrope with a point mass core that is mapped onto the simulation domain. The envelope has polytropic structural coefficient . Within , pressure, density, and sound speed are those of the stellar envelope. Outside, we impose a low-density, constant sound-speed atmosphere in hydrostatic equilibrium. The total mass of this background (the mass outside ) is roughly .

We simulate a system with a mass ratio , which implies and in our code units (similarly ). The primary-star mass, is made up of the excised core, and the simulated envelope gas. We use a 1D model (shown in Figure 2) in which 25% of the mass is within the inner 10% of radius – similar to a red giant branch star of in which the core mass is roughly one quarter the total mass. Our model adopts identical polytropic coefficient and gas adiabatic exponent, . We excise an inner region of , and include this mass in . The outer radius of the domain is . This implies and an envelope gas mass of approximately . In a grid-based hydrodynamic calculation we cannot embed the stars in vacuum, therefore, in our 1D profile, we apply a transition to a low-density, high-sound-speed medium at a pressure of , implying a total background mass (outside the star) of approximately .

Our fiducial model adopts an initial separation of , the analytic Roche limit separation for our mass ratio under the Eggleton (1983) approximation. The initial orbit is circular, . We set the donor envelope in synchronous rotation with the orbital motion such that and . The 1D profile is relaxed onto the computational mesh for 15 code time units (roughly 15 dynamical times of the donor star). Next, the gravity of is linearly turned on for 15 time units such that after 30 dynamical times, the simulation is fully active.

4.2 A Model of Runaway Binary Coalescence

Figure 3: Evolution of the orbit of the two stellar core particles, and , for the duration of our simulation. The upper panels show the orbital evolution in the simulation and inertial frames, respectively. The lower panel tracks the evolution of the orbital separation as a function of simulation time. In this panel the size of the star is marked with a dashed grey line, and the size of the inner boundary is marked with a solid grey line. The Roche Limit separation is marked with a blue line. Time is when the separation is equal to .

Our model system’s evolution begins at the Roche limit separation , and terminates at . Here we trace the evolution of the orbit and gas flow.

4.2.1 Orbit

As we see in Figure 3, the orbital evolution begins with many tightly-wound orbits with minimal change in separation, but eventually leads to a phase of dramatic inward spiral toward coalescence. The upper panels of Figure 3 show the path of particles and throughout the simulation. In the simulation frame, the donor-core particle, , is always at rest at the origin, and orbits, as shown in the black trajectory. In the inertial frame, both objects orbit around their mutual center of mass.

The evolution of orbital separation is most clearly seen in the lower panel of Figure 3. While the separation initially evolves slowly in time, it exhibits a runaway decrease to smaller separations and binary coalescense as the orbit shrinks (see similar behavior in simulations of Nandez et al. (2014), their Figure 2, and Iaconi et al. (2017), Figure 4). At separations , the rate of evolution of the orbit is rapid. We denote a reference time, , defined by the time when the separation equals the initial radius of the donor,


where we approximate the binary semi-major axis with the separation of the two core particles in what follows . This is a reasonable approximation because the orbit remains relatively circular.

Figure 4: Orbital decay rate, expressed as , the number of orbital periods over which the semi-major axis changes by of order itself. As the separation decreases from near the Roche separation to coalescence, the orbit goes from evolving very slowly to evolving rapidly, on the same order as the orbital period, .

Figure 4 further explores the accelerating rate of orbital decay noted in Figure 3. Here we parameterize this as , the orbital decay timescale divided by the instantaneous orbital period. The rate of change of the orbit becomes increasingly rapid as the binary comes closer to coalescence. The orbit first evolves slowly, with , for separations (90% of the Roche limit separation). When the separation is , we find , indicating that the orbit is rapidly plunging to tighter separations even prior to the engulfment of the secondary within the envelope of the donor. Finally, for separations , the orbital decay timescale is similar to the orbital period of the binary.

4.2.2 Gas flow

The consequences of this runaway decay of the binary orbit are dramatic. The binary interaction begins quite gently (and, in fact, the realistic system would include a much more gradual lead-in than is possible to simulate here). Even so, the system enters a phase of dynamical coalescence where the orbit is shrinking on timescales similar to the orbital period. This rapid coalescence indicates transition from a very long phase of quasi-hydrostatic evolution for the binary components to one where the evolution timescale is so short that equilibrium cannot be maintained and subsequent interaction is truly hydrodynamic in nature.

Figure 5: Slices of gas density through the orbital plane (top) and perpendicular plane (bottom) for the duration of our simulation. For clarity, we plot along and axes, which are rotated about the -axis relative to the simulation frame, such that the line of separation between and lies along . The origin in and is at the particle-gas system’s center of mass. Note that snapshots are much more closely spaced in the time series near the time of binary coalescence, .

The evolution to an increasingly dynamic encounter is visualized in Figure 5 with a time series of slices of gas density in the orbital midplane (upper panel) and perpendicular plane (lower panel). Slices are rotated relative to the frame in which they are computed (see Figure 1 and section 3.1) such that the -axis is parallel to the orbital separation vector, and shifted to have origin located at the particle-and-gas system center of mass. We note that the snapshots here are not separated equally in time, and we see that evolution in the gas properties evolve largely with binary separation, rather than time, during the system’s runaway coalescence.

Common features in the flow include exchange of mass from the donor envelope toward the accretor (Lubow & Shu, 1975) and then into tails which extend from near the outer (near the accretor) and (near the donor) Lagrange points (Shu et al., 1979; Pejcha et al., 2016b, a). This material forms a spiral pattern in the frame that rotates with the binary because, in moving radially outward, it no longer has sufficient angular momentum to maintain corotation with the binary (Shu et al., 1979). The changing pitch angle of this spiral reflects varying radial to orbital velocity ratios. At early times, the morphology of Roche lobe overflow and disk-formation near the secondary is reminiscent of simulations of mass transfer by low-velocity winds (e.g. Chen et al., 2017). As the system trends toward merger, the outflow morphology changes from a low-density stream originating largely near , to a broad, mass-rich fan of ejecta at . In slices perpendicular to the orbital plane, we observe that the outflow is concentrated toward the plane of the orbit. Plumes of trailed-off material from the and points can be seen to self-intersect and interact in this plane (Pejcha et al., 2016b, a), with those interactions shaping and broadening the outflow tails.

Based on the snapshots of Figure 5, we take a moment to re-consider the appropriateness of our choice of an adiabatic () equation of state for the simulated gas. A first useful property is the optical depth of material, , for a column of length . The magnitude of the optical depth depends on the rescaling of our dimensionless system to physical coordinates. If we adopt , then, in terms of donor properties,


where the first term represents density in simulation-scaled units. We note that varies by orders of magnitude over the density and temperature conditions that may be encountered. Nonetheless, for , densities and solar-mass donors, in code units, of represent optically-thick material when . This nominal optical depth scales with stellar properties, if , only densities of are optically thick. These scalings indicate, that, for more compact giant-star donors, much of the simulated material is optically thick, suggesting inefficient cooling, and that the adiabatic treatment of the gas thermodynamics is reasonable as a first approximation.

4.3 Mass Transfer and Angular Momentum Evolution

In order to understand the runaway orbital evolution described in the previous section, we examine the angular momentum budget of our simulated system and its flow between the particles and gas, mediated by gravitational torques.

4.3.1 System Angular Momentum

Figure 6: Decomposition of angular momentum in the particle-gas system throughout the simulation duration. Particle angular momenta decrease through exchange with the gas. The total angular momentum is conserved to within approximately 2% as the particles and gas exchange roughly a factor of two in angular momentum.

We begin by examining the total angular momentum of the particle and gas components of our simulated binary, the -component of of which is plotted in Figure 6. The lines plotted in this diagram show angular momenta relative to the system center of mass. The particle angular momentum represents the net angular momenta of the primary and secondary particles relative to the center of mass, as described by equation (24). The gas angular momentum is the total angular momentum of the gas on the computational mesh, equation (25), while the “gas off grid ” is the cumulative angular momentum carried by gas leaving the computational domain, which results from a time integral of equation (26). Finally, the total angular momentum of the system is the sum of these components, equation (27).

The particle angular momenta are observed to decrease, at first gradually, then rapidly, representing the orbital decay discussed previously. The net gas angular momentum mirrors this behavior – as angular momentum is lost from the particle pair it is given to the gas. Some of this gas leaves the domain during the simulation, carrying angular momentum with it, but the majority of the gas remains on the mesh. Importantly, we note that the total system angular momentum is nearly conserved during our simulation (conservation of angular momentum is not guaranteed in our method of solving the evolution equations for the particle-gas system, as discussed in Sections 3.6 and A.4). We find an approximately 2% change in total angular momentum as the angular momenta of the particles and gas exchange by a factor of approximately two.

4.3.2 Mass Flow and Gas Specific Angular Momenta

Figure 7: Mass loss rate from the donor star increases by orders of magnitude as the system tightens. Mass loss rate is measured from the simulation as the instantaneous flux of mass through a spherical surface with radius . Lines marked “simulation”, are smoothed; the instantaneous values are highly variable, especially at large separation, and are shown in grey in the background. The dashed line compares to the analytic formula for the expected mass-loss rate, equation (1) with constant of proportionality equal to one, and an assumed donor radius of , which is the equatorial radius of the rotating stellar profile following relaxation (see Figure 21). With these choices, the analytic expression produces a remarkable fit for the simulated result.

As the particle pair’s orbital angular momentum decays, angular momentum is clearly exchanged with the gas. It is interesting to begin to examine the details of this flow of angular momentum. Following the analytic description of mass transfer in binary systems outlined in Section 2, we break this down into a flux of mass and its distribution of specific angular momenta.

In Figure 7, we show the rate at which material is lost from the donor star. As the encounter proceeds toward a common envelope phase, where the stellar cores are engulfed in a shared envelope, the mass of the individual stars becomes less well defined. Here, we estimate the mass loss rate by computing the instantaneous flux of mass through a spherical surface at . In Figure 7, grey lines in the background show the instantaneous mass transfer rate, which is highly variable at large separations. The solid blue line shows a moving average of the instantaneous values (of varying width in time – widest at early times such that the bin width encompasses roughly equal mass at each separation).

The most remarkable feature of Figure 7 is how much the mass transfer rate increases as the binary proceeds toward coalescence. Over the simulation duration, approximately 34 initial orbital periods, the mass transfer rate increases by four orders of magnitude. The late-time rise is especially rapid; over the several dynamical times it takes for the binary separation to shrink from 1.5 to , the mass transfer rate increases by an order of magnitude.

Figure 7 also shows a comparison to the analytic formula to estimate mass transfer rate, equation (1). Under several assumptions, we achieve a remarkable fit between this expression and our simulation results. First, we adopt a constant of proportionality of 1 in equation (1). Next, we find that the fit is best represented if we use the volume-averaged radius of the oblate, rotating donor-star profile, which is as shown in Figure 21. With these calibrations, we find that the simple expression of equation (1) provides a very good description of the simulation result, especially at separations larger than .

Figure 8: Specific angular momentum of material in the orbital midplane, at the same times and rotated - frame as Figure 5. Contours show decades of density from to  . In each frame, material is drawn from low specific angular momentum near the donor to higher specific angular momentum near the accretor. In the earliest snapshot, material is trailed off with relatively high specific angular momentum, roughly from the outer Lagrange point (). As the separation tightens slightly in the subsequent snapshots the specific angular momentum of ejecta decreases, and we see an associated broadening of the outflow stream. In the lower panels at times near binary coalescence, even as the orbit shrinks rapidly, the specific angular momentum of ejecta are approximately constant as the evolving outflow morphology approximately compensates for the tightening orbit.

As outlined in Section 2, the determining factor in orbital evolution with mass loss is not just the rate at which mass is lost, but the angular momentum content it carries with it. Here we examine the specific angular momentum of material lost from the donor star’s envelope. We begin in Figure 8 by mapping the specific angular momentum of material in the orbital plane, in the same series of snapshots as Figure 5. A few useful features emerge from this diagram, particularly when compared to the over-plotted density contours. Across all the snapshots, we see that material making up the bulk of the donor envelope has relatively low specific angular momentum, typically less than . Material streaming from the donor toward the accretor acquires specific angular momentum through that gravitational interaction. This transition can be observed in the gradient of specific angular momenta leading up to the accretor position. Finally, some material is trailed off from the the binary, especially near the outer Lagrange point, . This outflow contains the highest specific angular momenta in the slices.

Figure 9: Specific angular momentum, , and dimensionless specific angular momentum, , as a function of binary orbital separation. The solid blue line shows smoothed simulation results, as in Figure 7, and the dashed lines show reference quantities for context. Material lost from the binary system carries specific angular momentum similar to the outer Lagrange point at the widest separations. This reduces to near the specific angular momentum of the accretor at separation of approximately , but then increases even as the orbit tightens, exceeding the angular momentum of even the point at separation of . Details on the calculation of these quantities are given in the text.

We quantify the changing specific angular momentum distributions that can be seen in Figure 8 further in Figure 9, where we estimate the average specific angular momentum of material lost as a function of binary orbital separation. We also compute the dimensionless quantity , which relates the specific angular momentum to the specific angular momentum of the binary, see equation (2). The average specific angular momentum of material carried away from the binary system is computed as , where is the component of the gas angular momentum about the center of mass, shown in Figure 6. We use smoothed versions of the mass transfer rate and angular momentum transfer rate, as described in Figure 7. The specific angular momenta of the accretor and points are also shown for context.

The lower panel of Figure 9 shows the conversion of the same diagram into the unitless quantity used to describe orbital evolution in equation (2), (the specific angular momentum of the material lost in units of the specific angular momentum of the binary). To make this conversion, we must approximate and , even when the system is coalescing and the separation of mass into two stars is increasingly unclear. In Figure 9, we adopt the mass within as the instantaneous value , and treat the angular momentum as the angular momentum of a circular binary at the current separation.

Despite the approximate nature described above of the quantities displayed in Figure 9, useful trends emerge. Near the Roche limit, material is trailed off from the binary with specific angular momentum similar to, or slightly larger than, the angular momentum of the outer Lagrange point (e.g. Shu et al., 1979; Pribulla, 1998). As the binary tightens, we see that and drop, reaching minima that are nearly a factor of two lower near a separation of or roughly 87% of the Roche limit separation. As the binary trends toward tighter separations we observe that the specific angular momenta of the escaping material actually increases, despite the tightening binary. Relative to the binary’s angular momentum, as shown with , this implies an even steeper increase. As might be inferred from Figure 8, the change in specific angular momenta of mass loss from the binary is reflected in the evolving morphology of outflows from the binary pair, a point to which we will return in Section 4.4.

Figure 10: Reconstructed orbital evolution from the analytic expression of equation (2) compared to the simulation result. In the analytic expression, we adopt the parameters described in Figure 7 for the mass-loss rate and the specific angular momentum of for material lost. With these choices adopted, the analytic description very successfully reproduces the runaway orbital evolution observed in the simulation.

In the figures and discussion above, we have measured each of the parameters one needs to assign to integrate the parameterized differential equation of orbital evolution under the influence of mass loss from a binary, equation (2). We use these ingredients from the simulation to reconstruct the orbital evolution of the binary, and compare to our full simulation result in Figure 10. In particular, we adopt the analytic mass transfer rate and associated constants as plotted in Figure 7, equation (1). Motivated by the results of Figure 9, we apply . The resulting orbital separation evolution is strikingly similar to the simulation result.

The similarity of these results offers insight into the physical mechanism of runaway coalescence in the binary system. The fact that the mass loss rate from the donor is an extremely steep function of separation combined with the fact that the typical angular momentum of ejecta is large – similar to that of the outer Lagrange point – indicates that runaway mass loss leads to runaway orbital decay.

4.3.3 Driving Orbital Decay: Torques

Figure 11: Torque volume density in slices through the rotated (-) orbital plane (see Figure 5). Positive torques indicate material pulling the binary components forward in their orbit, negative torques indicate material pulling backward. Imbalance in positive and negative torques generates a negative net torque on the system which drives it, increasingly, toward coalescence. In later snapshots, , as the donor’s envelope is distorted, the particle-gas system center of mass (black cross) moves away from the line of centers between the two particles, generating curved contours of positive and negative torque. At late times, the gas distribution around the accretor is highly asymmetric, leading to the net negative torque that drives orbital decay.

In this section, we turn to the gravitational forces within the binary system in order to understand how low specific angular momentum gas near the donor is pulled away and then expelled with higher specific angular momentum in the vicinity of the accretor. The backreaction of these forces drives the orbital evolution of the pair of stellar cores within their gaseous surroundings.

A key quantity is the torque on the pair of particles, which describes their rate of change of angular momentum. The gas distribution applies gravitational force to the two particles, the force per unit volume on particle from a parcel of fluid with density is


The torque on particle relative to the system center of mass as a result of this force also depends on its lever arm relative to the center of mass,


Integrating equation (31) over volume yeilds the net torque on particle . The net torque on the binary is the sum of the torques on each component, .

In Figure 11, we map out the -component of torque density on particles 1 and 2 relative to the particle-gas center of mass (torque per unit volume)333See Tang et al. (2017) for a similar analysis in the context of disk accretion onto binary supermassive black holes. We slice in the orbital midplane for the same snapshots as shown in Figures 5 and 8. The color bar in Figure 11 includes positive and negative values. Gaseous material colored in green is pulling the stellar cores forward in their orbital motion. By contrast, material colored in pink drags the system backward. Because the net torque (and thus rate of change of particle angular momenta – as show in Figure 6) is given by an integral over the full volume, much of the positive and negative torques within the domain cancel out. If the gas distribution (and quantity of pink versus green material) is perfectly symmetric, there is no net torque.

Figure 11 shows increasing asymmetry of positive and negative torquing material as the binary advances to smaller separations. This asymmetry is particularly pronounced about the accretor. Looking at the magnitude of torque density, we see strongly increasing magnitudes of torque densities near the accretor as the mass flux through the binary system increases. Both increased asymmetry and larger magnitudes, therefore, contribute to the increasing net torque, and accelerating binary coalescence.

Inference of the net torque from Figure 11 requires an “integral by eye” and does not include material above and below the orbital midplane. In Figures 12 and 13 we integrate the net torques and separate them into their components. We begin our analysis of net torques with Figure 12, which shows the -component of gravitational torque on particles 1 (donor core) and 2 (accretor) by the gas distribution with respect to the full system center of mass. The panels of this figure show the torque with respect to time (left panel) and orbital separation (right panel).

There are several immediate conclusions to draw from Figure 12. First, the magnitude of net torque is much smaller than the maximum torque densities represented in Figure 11. This indicates that the cancelation of positively versus negatively torquing gas is relatively efficient. Second, the majority of the net torque is on the accretor particle (2). Referencing the midplane maps of torque density in Figure 11, we see that this is not because of lack of torquing material near the donor, but because the donor envelope remains relatively symmetric for much of the coalescence and thus generates little net torque. Finally, since the net torque represents the rate of change of angular momentum, we see the imprint of the runaway orbital decay here as well. As the binary separation tightens, the torque increases, indicating an ever-more-rapid runaway of orbital decay.

Figure 12: Integrated torques from the gas distribution on the two particles as a function of time and separation. The majority of the net torque is applied to the accretor particle.

In Figure 13, we examine the spatial distribution of torques on the accretor particle versus binary separation. In addition to the net torque on particle 2, , we show the torque within spherical regions centered on the accretor particle of radii between 0.1 and 0.6. This decomposition indicates that much of the net torque on the accretor particle (and therefore also on the binary as a whole) originates in its immediate vicinity. By comparison to Figure 11, we see that this is the region that encompasses much of the asymmetry in positively and negatively-torquing gas. A similar conclusion is reached by Tang et al. (2017) in studying disk accretion onto binary black holes – significant torques are generated in the immediate vicinity of the black holes – implying that resolving these scales is important. As the separation decreases to be similar to , a larger fraction of the torque is generated outside , by the broad fan of outflowing ejecta. Nonetheless, even in this, most extreme, case, the majority of the torque does not come from the largest-scale flow distribution.

Figure 13: Contribution to torque on accretor particle, 2, from within radii from the accretor between 0.1 and 0.6. The bulk of the torque on the accretor, and therefore also on the binary system, is generated in the vicinity of the accretor.

The trend of the binary system toward runaway merger, studied in the orbital properties in Section 4.2 and in the angular momentum budget in Section 4.3.1, is a consequence of the forces acting on the stellar cores due to the distribution of gas originating from the donor envelope. As shown in Section 4.3.2, a dramatically-increasing flux of mass transfers from the donor to the vicinity of the accretor. Through gravitational torques originating primarily near the accretor particle, this material acquires specific angular momentum similar to (or somewhat larger than) that of the accretor particle (4.3.2), and is expelled from the binary. This accelerating loss of angular momentum leads to the runaway plunge of the accretor within the envelope of the donor and the onset of the common envelope phase of the binary system.

4.4 Changing Flow Morphology During Coalescence

We have noted above a changing outflow distribution as the binary orbit tightens, seen particularly clearly in Figure 5. Here we examine this changing morphology in more detail, and trace its origin to desynchronization of the orbital motion with donor-star rotation.

4.4.1 Evolving Flow Morphology

Figure 14: A comparison of flow properties at early and late stages in the binary coalescence in the - rotated orbital plane. The left panels show density, the center pressure, and the right specific entropy. The left and center panels overplot the instantaneous Roche lobes, while the right adds velocity vectors in the instantaneously corotating frame. A dramatic transformation in the morphology of flow through and outward from the binary can be seen in this comparison. At early times, material in the donor’s envelope is at rest in the corotating frame, and transfers to the accretor in a thin stream. Later, material in the donor’s envelope is in motion in the corotating frame, and fans outward following its interaction with the accretor.

We examine the transformation of flow properties in Figure 14. The upper panels show an early snapshot, approximately 300 dynamical times prior to coalescence, the lower panels show a snapshot near , when the accretor is plunging within the envelope of the donor. The left panels show density, the center show pressure, and the right show gas specific entropy. In the left and center panels we overlay the instantaneous Roche lobe. The Roche lobes are approximated by assuming the full donor mass and are shown centered on the system center of mass, with contours corresponding to the potential of the , , and Lagrange points, respectively – clearly, these are a more meaningful guide in the early snapshot in which the binary is still separated. In the right panels, we plot gas velocities in the instantaneously corotating frame.

The biggest distinction in flow morphology between early and late snapshots comes in the nature of gas flow between and outward from the two binary components. In the early frame (upper panels of Figure 14), the highest density gas is largely confined to the Roche lobes of the two stars, and gas density and pressure fall off roughly tracing contours of equipotential. A relatively thin stream of stellar material transfers from the donor to the accretor. This material assembles into a rotating flow with a high pressure and high specific entropy core within the accretor-particle softening radius. Surface layers of higher-specific entropy gas surround both stars, similar to the description of Shu et al. (1976). Material outflowing from the binary has relatively high specific entropy, and exits near the point.

In the later snapshot (show in the lower panels of Figure 14), material is flung outward in a broad fan of ejecta. The width of this fan is similar to its size. Comparison with Figure 5 shows that that ejecta also fan out in opening angle in the perpendicular plane, and that the leading edge of this material quickly catches the trailing. Material is unconfined by the instantaneous Roche lobes of the pair, and its specific entropy shows that it is largely unshocked or mildly shocked in its ejection, by contrast to the much higher entropies observed near in the early snapshot.

4.4.2 Desynchronization

A crucial clue to understanding the transformation of flow morphology lies in the velocity vectors plotted in Figure 14. Early in the simulation, material within the donor envelope is nearly at rest in the instantaneously-corotating frame, and it acquires velocity in falling toward the accretor. The later snapshot shows that the donor envelope is in motion in the corotating frame. Although the donor star started out in synchronous, solid-body rotation with the orbit, orbital decay has desynchronized the angular frequencies of the orbit and the donor’s envelope. In the corotating frame, donor envelope material appears to counterrotate. It is ejected on a range of trajectories that seem to depend primarily on impact parameter relative to the accretor.

Figure 15: Angular frequency of material about the donor-star origin, , relative to the initial, , or instantaneous, , orbital velocity in the rotated orbital midplane. Contours show density in half-dex intervals. The gas is initialized in synchronous rotation with the orbit, but desynchronizes with the orbit over time: the gas retains close to its initial rotation even as the orbit shrinks and the orbital frequency increases. Thus, relative to the instantaneous orbital frequency (lower panel), material in the donor envelope appears increasingly counter-rotating.

Figure 15 explores desynchronization of orbit and donor envelope further through mapping the angular velocity of material in the orbital plane. The color scale shows material’s angular velocity about the donor, , relative to the initial (, upper panel), or instantaneous, (, lower panel) orbital angular velocity. In the earliest snapshots, we observe that the donor envelope starts in synchronous rotation with the binary orbit. The upper panels show that the angular velocity of material in the orbital plane is largely unchanged, even as the orbit decays. With the exception of strong flow near the accretor, the primary feature in angular velocity maps at the time of coalescence is a quadrupolar perturbation of the velocity field due to the gravitational pull of the donor. When we shift our reference frequency to that of the instantaneous orbital frequency (lower panels), Figure 15 shows donor envelope material counter-rotating in the instantaneous orbital frame.

Figure 16: Angular frequency of orbital motion and donor star gaseous envelope rotation. The envelope gas initial rotates at the orbital frequency at the Roche limit. As the separation decreases, the orbital frequency goes up, roughly as , but the envelope still rotates at approximately its original frequency, increasing slightly only when . Orbital angular momentum lost during orbital decay does not go into spinning up the envelope prior to coalescence. This implies increasing desynchronization of the accretor and the neighboring envelope gas, such that by the time , the secondary skims across the surface of the donor star’s envelope.

By maintaining close to its original rotation rate, the donor envelope loses its initial synchronization with the shrinking orbit. Such behavior is qualitatively described as a “loss of corotation” by Podsiadlowski (2001) and Ivanova et al. (2013b). Figure 16 shows this desynchronization quantitatively. We compare the orbital angular frequency to the donor envelope angular frequency as a function of orbital separation. To do so, we compute the mass-averaged angular velocity of envelope material, defined crudely as material with . We observe that the orbital angular frequency increases as the separation shrinks over the course of the simulation. This behavior is close to the scaling of two point masses, but is slightly shallower because of the decreasing enclosed mass within the orbital separation (both due to mass ejection and as ). This increasing orbital frequency is not matched by the angular frequency of the envelope, which is relatively constant. Therefore, while the system begins mass transfer in synchronous rotation, this state is rapidly lost as the orbital separation shrinks toward .

4.4.3 Outflow Morphology in Desynchronized System: Comparison to Ballistic Trajectories

Figure 17: Comparison of ballistic trajectories in binary potential to the gas flow in the rotated, -, orbital midplane. We plot randomly sampled free-particle trajectories in the corotating frame on top of gas specific entropy (color) and velocity in the corotating frame (vectors). Free trajectories are initialized from the flow near the point. Free trajectories mirror the evolving flow morphology in the binary. Early on, ballistic trajectories are confined to the accretor’s Roche lobe. In later snapshots, desynchronization of the orbital motion relative to the donor’s rotation implies that material crossing from donor to accretor is flung out in a broad fan that closely mirrors free trajectories.

Loss of synchronization between donor rotation and orbit carries implications for gas crossing from the donor toward the accretor, and by consequence, for the morphology of outflow from the binary. The manner in which material overflows from the binary is particularly affected. Figure 17 compares ballistic trajectories integrated in the potential of point masses of mass and in circular orbit at the orbital separation with the fluid flow (showing specific entropy in color scale and velocity with vectors). Ballistic trajectories are initialized at randomly-sampled positions near the point, with velocities interpolated from the simulation mesh.

In early snapshots, Roche lobe overflow leads to the formation of rotating structures around the accretor. Free trajectories self-intersect within the Roche lobe. In the gas dynamic calculation, gas streams collide and circularize, increasing the specific entropy of the gas through shocks. In either case, as discussed in the context of Figure 14, material is largely confined by the Roche contours of equipotential. Material that is lost is trailed off from the point after being shock-heated to higher entropy and pushed over the potential barrier by the pressure gradient away from the accretor.

At times approaching , the stream of gaseous ejecta becomes broader and also progressively lower-entropy. Relatedly, an increasing number of free-trajectories exit the Roche lobes. By all of the free trajectories are lost and carried away in the tail, and the free particle streamlines nearly mirror the gas velocities.

These features indicate that the mode of mass shedding from the binary transforms as the binary desynchronizes: from being pushed out by pressure gradients to being flung out along ballistic trajectories. This changing mode of mass loss leads to the change in outflow morphology from a narrow stream to a broad fan.

4.5 Non-Synchronous Initial Rotation

Figure 18: Comparison of orbital evolution under different cases of initial donor rotation. We initialize simulations with and , in addition to our fiducial case of . This comparison shows that non-synchronous initial rotation at the onset of Roche lobe overflow leads to somewhat more rapid orbital decay. The difference is most pronounced at large separations where the initially-synchronous cases remain synchronized and vanishes as the orbits tighten and all cases naturally desynchronize.

We have discussed the role of desynchronization of the orbital motion with donor rotation in giving rise to highly asymmetric flow around the accretor, which enhances the net torque on the binary and drives it toward merger. We might hypothesize that initial stellar rotation that is non-synchronous would enhance this effect and lead to more rapid orbital evolution. To study the magnitude of this effect, we ran simulations otherwise identical to our fiducial case with , but which had envelopes with initial angular velocities of and .

As we show in Figure 18, initial rotation does have an effect on the rate of orbital decay, with more rapid coalescence in the non-synchronous cases. The effect is most pronounced at large separations – at which the synchronous case has not yet fallen far from co-rotation. By the time , all cases exhibit similar behavior because even the synchronous case has desynchronized.

Overall, the magnitude of the effect of the initial spin of the star is relatively mild. This indicates that the changing flow morphology due to synchronous or non-synchronous rotation is a secondary aspect of runaway orbital coalescence. Instead, the the dominant aspect of the runway coalescence is the increasingly mass-rich outflow as the secondary carves deeper into the envelope of the primary star (Figures 7 and 11).

5 Connection to Astronomical Transients

The key question that motivated this study is whether impulsive, luminous transients can originate from binary systems at the onset of common envelope episodes. We have shown that runaway orbital decay leads to increasing mass ejection rate from the binary during the final plunge toward coalescence. Here we argue that this runaway orbital decay and mass ejection are the mechanisms that lead to the production of “luminous red novae” transients at the onset of common envelope phases.

5.1 Impulsive Outbursts in Stellar Coalescence Transients

Luminous red novae transients are thought to originate in the coalescence of binary systems. The strongest evidence for this identification is the galactic transient V1309 Sco, which was an eclipsing binary with decreasing orbital period (Tylenda et al., 2011) that exhibited a 10 magnitude outburst (Mason et al., 2010), after which no periodicity remained (Tylenda et al., 2011). Much of the rest of the luminous red novae class of transients are assembled on spectroscopic and photometric similarity to this transient.

Several recent extragalactic transients have offered new insights through particularly well-sampled, multicolor light curves. M31 LRN 2015, which reached peak brightness in January 2015 (Kurtenkov et al., 2015), has high-cadence multicolor photometry during the rise to peak brightness(Williams et al., 2015; Kurtenkov et al., 2015). A progenitor source is detected in Hubble Space Telescope data a decade prior to outburst (Dong et al., 2015; Williams et al., 2015) and in other archival datasets (Dong et al., 2015). Another transient, M101OT 2015-1 (Blagorodnova et al., 2017) was detected in outburst just a month later, in February 2015. The rise to peak occurred while the object was behind the sun, but there is a decade of exquisite multi-color photometry of the luminous pre-outburst source (Blagorodnova et al., 2017)

A surprising lesson from the M31 LRN 2015 transient was that the rise to peak brightness occurs over nearly the same timescale as the orbital period around the surface of the giant star observed in pre-outburst data (MacLeod et al., 2017). Another universal feature of these transients is precursor emission and a gradual brightening over tens to hundreds of orbital periods prior to coalescence, which was first commented on in the light curve of V1309 Sco (Tylenda et al., 2011). This precursor brightening is seen to be accompanied by slightly declining photosphere temperature in the multicolor photometry of M101OT 2015-1(Blagorodnova et al., 2017).

Long-term precursor emission suggests an interaction that occurs steadily over a number of orbits prior to the outburst. The progressive brightening indicates that the interaction intensifies over this time. These precursor signatures rule out a stellar collision or very eccentric encounter as an explanation for the impulsive nature of the outbursts. A relatively-circular binary system is the simplest explanation – and yet, how can such a system be stable over a stellar evolution timescale, perhaps billions of orbits, and then power an outburst over the course of a single orbit?

5.2 Connection to Simulated System

Here we argue that precursor emission and outburst share a common origin in the runaway coalescence of a circular binary.

Precursor emission: Mass transfer in our simulated system soon leads to outflows from the inner binary pair as material is unable to pile-up in the vicinity of the accretor. Pressure gradients divert flow outward from the binary, and it is trailed off from the inner pair near the outer, and Lagrange points (as described by Shu et al. (1979) and seen, for example, in the earliest snapshot of Figure 5).

We have extensively discussed the role of outflows in carrying angular momentum from the binary, but this outflow also leads to an expanding photosphere radius. Pejcha et al. (2016b, a) have modeled the gas dynamics and thermodynamics of outflows from the point of a binary system, and have applied this modeling to very convincingly explain the pre-outburst evolution of the V1309 Sco transient (Pejcha, 2014; Pejcha et al., 2017). In particular, Pejcha et al. (2017) showed that the evolving eclipse morphology in the decade prior to outburst can be traced to increasing (orbital-phase-dependent) obscuration due to increasing outflow. Pejcha (2014) and Pejcha et al. (2017) show that periodicity vanishes as V1309 Sco brightens because the photosphere exceeds the radius of the binary.

Here we have shown that these outflows arise naturally as a consequence of mass exchange, that they intensify significantly as the binary separation shrinks, and that the torques associated with generating the outflow drive the evolution of the binary orbit. We therefore argue that heated ejecta and an expanding photosphere power the observed precursor emission.

Impulsive outburst: Two features of our simulation explain the ability of coalescing binaries to produce impulsive outbursts. We observe a runaway decay of binary orbital separation (Figure 3) that becomes so rapid that the decay timescale is similar to the orbital period, , when (Figure 4). Associated with (and directly responsible for) the runway orbital decay, mass ejection from the binary intensifies rapidly in the final moments of the system’s coalescence.

Figure 19: Zoom-in on orbital decay and donor mass-loss rate in the final dynamical times prior to binary coalescence. As the orbit plunges to smaller separations within a few dynamical times of , the mass ejection rate from the donor increases by an order of magnitude. We suggest that the origin of impulsive outbursts in luminous red novae transients is the rapid escalation of the mass ejection rate when the binary separation is similar to .

We summarize simultaneous orbital decay and mass loss with a zoom-in to the final dynamical times pre-coalescence in Figure 19. Over several donor-star dynamical times, the orbital separation goes from to less than and the mass loss rate from the donor star increases by an order of magnitude. Impulsive mass ejection accompanies impulsive orbital decay as the accretor star plunges within the envelope of the donor. We argue that the outburst phase of luminous red novae transients originates from this phase of runaway coalescence leading to the engulfment of the accretor and the onset common envelope phase.

Finally, we note a changing outflow morphology as the binary coalesces – from a thin stream to a broadening fan of essentially ballistic ejecta (Figure 17). The detailed dynamics and thermodynamics of this outflow surely imprint themselves on outbursts of luminous red novae transients.

6 Conclusion

Many open questions remain – we have considered just a single binary configuration – and many physical processes have been neglected in our model, which only treats the hydrodynamics of gravitational interaction between donor-star gas and the stellar cores. Nonetheless, from the examination of this simplified system, we draw some broad conclusions that carry implications for the coalescence of binary systems at the onset of stellar mergers and common envelope episodes.

  1. Our simulation of unstable mass transfer in a binary system begins at the Roche limit separation from a primary-star donor of mass onto an accretor of mass . The binary separation decreases, at first slowly, then increasingly rapidly. The orbit eventually enters into a phase of runaway decay, plunging from separation of to less than in roughly an orbital period (Figures 3 and 4).

  2. Angular momentum lost from stellar cores transfers to the gas pulled from the donor’s envelope (Figure 6) and is ejected with average specific angular momentum between that of the accretor and the point (Figure 9).

  3. The mass loss rate (Figure 7) and specific angular momentum of material lost (Figure 9) allow one to analytically reconstruct orbital evolution, using equation (2), which produces a remarkable fit to the simulated result (Figure 10). This reconstruction indicates that escalating mass loss drives much of the runaway orbital decay.

  4. The gravitational torques (Figure 11) that exchange angular momentum between the particles and gas are mostly applied by the accretor particle (Figure 12) in its immediate surroundings (Figure 13).

  5. The morphology of outflow from the binary changes as the orbital separation tightens, from a thin, high-entropy stream from the and outer Lagrange points to a broad fan of low-entropy material on largely ballistic trajectories (Figures 5, 14, and 17). Desynchronization of the orbital motion with donor-star rotation (Figure 16) is the origin of this transition in outflow properties (Figure 17).

  6. Characteristic features of luminous red novae transients – precursor emission and impulsive outbursts – find common origin in runaway orbital decay and mass ejection leading into a common envelope phase. The orbit shrinks from well-separated to engulfed and the mass ejection rate increases by an order of magnitude in a single orbital period at separations similar to the radius of the donor (Figure 19).

Our study of the dynamics of coalescing binaries was motivated by surprises in new time-domain observational data of luminous red novae transients. Future efforts to expand these samples are underway in the form of searches for low-temperature extragalactic optical/infrared transients (e.g. Kasliwal et al., 2017; Adams et al., 2018). Future simulation work can tighten the connection to the observable properties of luminous red novae transients through the progressive inclusion of currently-unmodelled physical effects like magnetohydrodynamics, more realistic equations of state with partial ionization, and radiative transfer. More broadly, our conclusions suggest the powerful role new empirical constraints are playing in opening a window into long-uncertain physics of common envelope events – especially through comparison to simulated counterparts of these systems.

We gratefully acknowledge many colleagues for helpful conversations that shaped this work, especially Andrea Antoni, Eric Blackman, Nadia Blagorodnova, Matteo Cantiello, Orsola De Marco, Rosanne Di Stefano, Paul Duffell, Adam Frank, Jonathan Grindlay, Tomasz Kaminski, Mansi Kasliwal, Abraham Loeb, Phillip Macias, Ariadna Murguia-Berthier, Priya Natarajan, Enrico Ramirez-Ruiz, Melinda Soares-Furtado, Jennifer Sokoloski, Rich Townsend, Scott Tremaine, Matias Zaldarriaga, and the “Evolved stars and binaries” discussion group at the CfA. This research made use of astropy, a community developed core Python package for Astronomy (Astropy Collaboration et al., 2013). M.M. is grateful for support for this work provided by NASA through Einstein Postdoctoral Fellowship grant number PF6-170169 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. Support for program #14574 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The work of E.C.O. was supported by a grant from the Simons Foundation (grant no. 510940). Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.


  • Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. D, 93, 122003
  • Abbott et al. (2016b) —. 2016b, Physical Review Letters, 116, 241103
  • Abbott et al. (2017a) —. 2017a, Physical Review Letters, 118, 221101
  • Abbott et al. (2017b) —. 2017b, ApJ, 851, L35
  • Abbott et al. (2017c) —. 2017c, Physical Review Letters, 119, 161101
  • Adams et al. (2018) Adams, S. M., Blagorodnova, N., Kasliwal, M. M., et al. 2018, PASP, 130, 034202
  • Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  • Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407
  • Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223
  • Blagorodnova et al. (2017) Blagorodnova, N., Kotak, R., Polshaw, J., et al. 2017, ApJ, 834, 107
  • Bodenheimer & Taam (1984) Bodenheimer, P., & Taam, R. E. 1984, ApJ, 280, 771
  • Chen et al. (2017) Chen, Z., Frank, A., Blackman, E. G., Nordhaus, J., & Carroll-Nellenback, J. 2017, MNRAS, 468, 4465
  • Dan et al. (2011) Dan, M., Rosswog, S., Guillochon, J., & Ramirez-Ruiz, E. 2011, ApJ, 737, 89
  • De Marco & Izzard (2017) De Marco, O., & Izzard, R. G. 2017, PASA, 34, e001
  • Dong et al. (2015) Dong, S., Kochanek, C. S., Adams, S., & Prieto, J.-L. 2015, The Astronomer’s Telegram, 7173
  • D’Souza et al. (2006) D’Souza, M. C. R., Motl, P. M., Tohline, J. E., & Frank, J. 2006, ApJ, 643, 381
  • Edwards & Pringle (1987) Edwards, D. A., & Pringle, J. E. 1987, MNRAS, 229, 383
  • Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  • Galaviz et al. (2017) Galaviz, P., De Marco, O., Passy, J.-C., Staff, J. E., & Iaconi, R. 2017, ApJS, 229, 36
  • Hernquist & Katz (1989) Hernquist, L., & Katz, N. 1989, ApJS, 70, 419
  • Huang (1956) Huang, S. S. 1956, AJ, 61, 49
  • Huang (1963) Huang, S.-S. 1963, ApJ, 138, 471
  • Iaconi et al. (2017) Iaconi, R., Reichardt, T., Staff, J., et al. 2017, MNRAS, 464, 4028
  • Iben & Livio (1993) Iben, Jr., I., & Livio, M. 1993, PASP, 105, 1373
  • Ivanova et al. (2013a) Ivanova, N., Justham, S., Avendano Nandez, J. L., & Lombardi, J. C. 2013a, Science, 339, 433
  • Ivanova et al. (2013b) Ivanova, N., Justham, S., Chen, X., et al. 2013b, A&A Rev., 21, 59
  • Kalogera et al. (2007) Kalogera, V., Belczynski, K., Kim, C., O’Shaughnessy, R., & Willems, B. 2007, Phys. Rep., 442, 75
  • Kasliwal et al. (2017) Kasliwal, M. M., Bally, J., Masci, F., et al. 2017, ApJ, 839, 88
  • Katz et al. (2016) Katz, M. P., Zingale, M., Calder, A. C., et al. 2016, ApJ, 819, 94
  • Kurtenkov et al. (2015) Kurtenkov, A. A., Pessev, P., Tomov, T., et al. 2015, A&A, 578, L10
  • Lombardi et al. (2011) Lombardi, Jr., J. C., Holtzman, W., Dooley, K. L., et al. 2011, ApJ, 737, 49
  • Lubow & Shu (1975) Lubow, S. H., & Shu, F. H. 1975, ApJ, 198, 383
  • MacLeod et al. (2017) MacLeod, M., Macias, P., Ramirez-Ruiz, E., et al. 2017, ApJ, 835, 282
  • Mason et al. (2010) Mason, E., Diaz, M., Williams, R. E., Preston, G., & Bensby, T. 2010, A&A, 516, A108
  • Metzger & Pejcha (2017) Metzger, B. D., & Pejcha, O. 2017, MNRAS, 471, 3200
  • Meyer & Meyer-Hofmeister (1979) Meyer, F., & Meyer-Hofmeister, E. 1979, A&A, 78, 167
  • Motl et al. (2017) Motl, P. M., Frank, J., Staff, J., et al. 2017, ApJS, 229, 27
  • Motl et al. (2007) Motl, P. M., Frank, J., Tohline, J. E., & D’Souza, M. C. R. 2007, ApJ, 670, 1314
  • Motl et al. (2002) Motl, P. M., Tohline, J. E., & Frank, J. 2002, ApJS, 138, 121
  • Nandez et al. (2015) Nandez, J. L. A., Ivanova, N., & Lombardi, J. C. 2015, MNRAS, 450, L39
  • Nandez et al. (2014) Nandez, J. L. A., Ivanova, N., & Lombardi, Jr., J. C. 2014, ApJ, 786, 39
  • Nicholls et al. (2013) Nicholls, C. P., Melis, C., Soszyński, I., et al. 2013, MNRAS, 431, L33
  • Ohlmann et al. (2016a) Ohlmann, S. T., Röpke, F. K., Pakmor, R., & Springel, V. 2016a, ApJ, 816, L9
  • Ohlmann et al. (2017) —. 2017, A&A, 599, A5
  • Ohlmann et al. (2016b) Ohlmann, S. T., Röpke, F. K., Pakmor, R., Springel, V., & Müller, E. 2016b, MNRAS, 462, L121
  • Paczynski (1976) Paczynski, B. 1976, in IAU Symposium, Vol. 73, Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, 75
  • Paczyński & Sienkiewicz (1972) Paczyński, B., & Sienkiewicz, R. 1972, Acta Astron., 22, 73
  • Pakmor et al. (2012) Pakmor, R., Edelmann, P., Röpke, F. K., & Hillebrandt, W. 2012, MNRAS, 424, 2222
  • Pejcha (2014) Pejcha, O. 2014, ApJ, 788, 22
  • Pejcha et al. (2016a) Pejcha, O., Metzger, B. D., & Tomida, K. 2016a, MNRAS, 461, 2527
  • Pejcha et al. (2016b) —. 2016b, MNRAS, 455, 4351
  • Pejcha et al. (2017) Pejcha, O., Metzger, B. D., Tyles, J. G., & Tomida, K. 2017, ApJ, 850, 59
  • Podsiadlowski (2001) Podsiadlowski, P. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 229, Evolution of Binary and Multiple Star Systems, ed. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, 239
  • Postnov & Yungelson (2014) Postnov, K. A., & Yungelson, L. R. 2014, Living Reviews in Relativity, 17, 3
  • Pribulla (1998) Pribulla, T. 1998, Contributions of the Astronomical Observatory Skalnate Pleso, 28, 101
  • Rasio & Shapiro (1995) Rasio, F. A., & Shapiro, S. L. 1995, ApJ, 438, 887
  • Ricker & Taam (2008) Ricker, P. M., & Taam, R. E. 2008, ApJ, 672, L41
  • Ricker & Taam (2012) —. 2012, ApJ, 746, 74
  • Sandquist et al. (1998) Sandquist, E. L., Taam, R. E., Chen, X., Bodenheimer, P., & Burkert, A. 1998, ApJ, 500, 909
  • Shiber et al. (2017) Shiber, S., Kashi, A., & Soker, N. 2017, MNRAS, 465, L54
  • Shiber & Soker (2017) Shiber, S., & Soker, N. 2017, ArXiv e-prints, arXiv:1706.00398
  • Shore et al. (1994) Shore, S. N., Livio, M., van den Heuvel, E. P. J., Nussbaumer, H., & Orr, A., eds. 1994, Interacting binaries
  • Shu et al. (1979) Shu, F. H., Anderson, L., & Lubow, S. H. 1979, ApJ, 229, 223
  • Shu et al. (1976) Shu, F. H., Lubow, S. H., & Anderson, L. 1976, ApJ, 209, 536
  • Smith (2011) Smith, N. 2011, MNRAS, 415, 2020
  • Smith et al. (2016) Smith, N., Andrews, J. E., Van Dyk, S. D., et al. 2016, MNRAS, 458, 950
  • Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  • Taam & Bodenheimer (1989) Taam, R. E., & Bodenheimer, P. 1989, ApJ, 337, 849
  • Taam & Bodenheimer (1991) —. 1991, ApJ, 373, 246
  • Taam et al. (1978) Taam, R. E., Bodenheimer, P., & Ostriker, J. P. 1978, ApJ, 222, 269
  • Taam et al. (1994) Taam, R. E., Bodenheimer, P., & Rozyczka, M. 1994, ApJ, 431, 247
  • Taam & Sandquist (2000) Taam, R. E., & Sandquist, E. L. 2000, ARA&A, 38, 113
  • Tang et al. (2017) Tang, Y., MacFadyen, A., & Haiman, Z. 2017, MNRAS, 469, 4258
  • Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170
  • Terman & Taam (1996) Terman, J. L., & Taam, R. E. 1996, ApJ, 458, 692
  • Terman et al. (1994) Terman, J. L., Taam, R. E., & Hernquist, L. 1994, ApJ, 422, 729
  • Terman et al. (1995) —. 1995, ApJ, 445, 367
  • Toonen & Nelemans (2013) Toonen, S., & Nelemans, G. 2013, A&A, 557, A87
  • Tylenda & Kamiński (2016) Tylenda, R., & Kamiński, T. 2016, A&A, 592, A134
  • Tylenda et al. (2011) Tylenda, R., Hajduk, M., Kamiński, T., et al. 2011, A&A, 528, A114
  • Webbink (1977) Webbink, R. F. 1977, ApJ, 211, 881
  • Williams et al. (2015) Williams, S. C., Darnley, M. J., Bode, M. F., & Steele, I. A. 2015, ApJ, 805, L18
  • Yorke et al. (1995) Yorke, H. W., Bodenheimer, P., & Taam, R. E. 1995, ApJ, 451, 308

Appendix A Numerical Convergence Studies

In this section we present a series of tests validating the convergence and robustness our simulation results with respect to numerical parameters.

a.1 Donor Hydrostatic Equilibrium

Figure 20: One-dimensional tests of hydrostatic equilibrium at varying spatial resolution. Zones are logarithmically-spaced in radius, and therefore have constant .
Figure 21: Three-dimensional test of hydrostatic equilibrium at the same resolution as our fiducial simulation. The left panels show polar slices through the simulation domain 100 dynamical times after the simulation is initialized. In the velocity slice, the black contour shows a density of . The right-hand panels compare profiles near the pole to those near the equator with the highest-resolution 1D profile.

The ability to preserve a structure in hydrostatic equilibrium with the hydrodynamic algorithm is a key aspect of our simulation method. Stars feature very small scale heights of density and pressure near their limbs, which pose significant challenges for any multi-dimensional approach. Approaches with too low resolution near the limb artificially extend the star’s atmosphere with dramatic effects on the phase of early, low-level mass transfer that we simulate here. Lagrangian and Lagrangian-Eulerian approaches tend to have very poor resolution near the stellar surface, leading the radius of the numerical star to be larger by a significant factor than the original model (see, for example, Figures 7-10 of Ohlmann et al., 2017, which show doubling of the radius of low-density material near the stellar limb). Approaches in Eulerian codes with cartesian meshes struggle to capture the steepest gradients along non-cardinal directions of the mesh.

For the above reasons, we have adopted a spherical polar coordinate system centered on the donor star in our simulation. Figure 20 shows 1D tests of hydrostatic equilibrium of an envelope with no perturber object. The domain is located along the midplane at and the envelope is initialized with spin frequency equal to that of a test mass at a separation of . Zones are logarithmically spaced in the direction such that is constant (with total extent in from 0.3 to 30). We vary the total number of radial zones from 96 () to 1536 (). Profiles are initialized with a relaxation phase of 15 dynamical times, and the simulations run to 100 dynamical times. This test examines the resolution needed to preserve donor hydrostatic equilibrium in the radial direction.

The left panels of Figure 20 show the full stellar envelope profile, the right panels zoom in near the stellar limb, where differences are most apparent. From top to bottom, panels show logarithm of density, logarithm of pressure, and velocity. Profiles are plotted 100 envelope dynamical times after initialization. The radius of the relaxed stellar profile is larger than because of centrifugal forces not accounted for in the initial, spherical model. The relaxation method allows a smooth transition to the hydrostatic rotating solution. All resolutions above the mid-level show reasonable agreement, hence we conclude that the converged hydrostatic solution is attained.

Next, we test the hydrostatic equilibrium of an identical envelope but in the full 3D domain. We apply our fiducial resolution used in the science runs, with 12 () x 8 () x 16 () mesh blocks of zones and two levels of nested static mesh in the equatorial plane. This corresponds to at the base level (near the poles) and at the maximum level (near the equator). Figure 21 displays the results of this calculation, 100 dynamical times after initialization. Slices are plotted in the vertical axis, such that we see the oblate stellar profile due to rotation. The radial velocity map (with contour showing a density of ) shows some residual flow inside (and especially outside) the star, where the extra degrees of freedom relative to the 1D simulation result in more fluid motion. We also note that this profile, with is marginally convectively unstable under the Schwarzschild criterion (that states that convective stability is realized when ).

The right panels of Figure 21 compare the highest-resolution 1D calculation to profiles from the 3D run along the pole and equator. The 3D star is somewhat more oblate than its 1D counterpart due to the 3D effect of additional fluid pressing toward the equatorial belt from latitudes above and below the equator.

The overall result of the tests presented in this section is that, at sufficiently high spatial resolution, the hydrostatic profile of the donor star can be preserved with a high degree of stability. This ensures that any evolution of the binary system is a direct result of the binary’s gravitational interaction, not the artificial erosion of the original profile of the stellar model.

a.2 Numerical Parameters and Spatial Resolution

Figure 22: Varying the inner boundary condition radius within the donor star. Left panel shows orbital separation as a function of time (relative to , when the separation is ). The right panel shows total (gas plus particle) system angular momentum as a function of binary separation. Ideally, total angular momentum should be conserved, the amount of departure from conservation is a useful benchmark of the simulations (detailed discussion of reasons for imperfect conservation in the text).
Figure 23: Same as Figure 22, but varying the softening radius around the accretor particle.
Figure 24: Same as Figure 22, but varying spatial resolution. We preserve the same 12 () x 8 () x 16 () base-level mesh block structure and two levels of nested static mesh refinement, but change the number of zones in each block as labeled in the legend.

In this section we test the robustness of our numerical results against several model parameters. These are

  1. the inner boundary radius, controlling how much of the star is excised, ,

  2. the softening radius of the secondary star’s gravitational point mass, ,

  3. and the spatial resolution of our mesh.

In what follows we describe tests against each of these parameters that justify the choices adopted for science runs presented in the main text of this paper. Our fiducial model for this section has an initial separation of , an inner boundary of , a softening radius of , and 12 () x 8 () x 16 () base-level mesh blocks of zones, with additional two levels of static mesh refinement about the orbital plane.

In Figure 22 we test the sensitivity of the model to the extent of stellar envelope excised from the center. The left panel of Figure LABEL:fig:conv_in examines the orbital separation evolution, the right the total angular momentum as a function of separation of the two point masses. In theory this should be conserved, in practice, it is not perfectly – for more discussion see Section LABEL:sec:Lconvergence. We see that when the excision is less than , the orbital evolution is converged with respect to smaller excision regions. This suggests that material at those smallest radii within the envelope is not significantly distorted from spherical within the extent of the calculation, so its replacement with a point mass and boundary has little effect on the outcome.

We test the softening radius of the secondary, accretor particle, adopting . We test values from to . These range of softening lengths range from 4 zones across the softening radius at widest binary separation for the smallest softening length to 16 zones across the largest softening length. Figure 23 demonstrates that values of the softening length less than have little effect on the simulation outcome.

Finally, we test for convergence of results against spatial resolution of the computational mesh. We do this by fixing the other parameters, and the base mesh and statically-refined mesh block structure, while changing only the number of zones per mesh block. This has the effect of holding everything fixed except the spatial resolution of the calculation. We show a range of resolutions from zones per mesh block to zones per mesh block in Figure 24. This lowest resolution barely preserves the hydrostatic equilibrium of our chosen stellar profile, while higher resolutions become prohibitively expensive to run. Across this range, we show that neither the runaway orbital behavior nor the overall conservation of angular momentum are sensitive to spatial resolution.

a.3 Sensitivity to Initial Binary Separation

Figure 25: Test of varying initial binary separation, otherwise similar to Figure 22.

Finally we consider how the initial binary separation may alter the results of binary coalescence. To test sensitivity to our adopted initial separation ( in the tests above and in our science runs), we consider a range from to . The deeper into the phase of runaway mass transfer the simulation is initialized, the less time the simulation needs to run leading into the coalescence of the pair of objects. We run calculations with identical initial stellar spin frequency (equal to the orbital frequency at the analytic Roche limit separation of ), but orbital separations of , 2.0, 1.9, 1.85 . The left panel of Figure 25 demonstrates that the overlapping stages of runaway orbital evolution proceed very similarly in each of these cases. We note that the different runs of Figure 25 have different initial angular momenta because of their differing initial binary separations.

a.4 Conservation of Total Angular Momentum

Finally, we comment on the degree of non-conservation of total angular momentum, equation (27), observed in the models in Figure 6 and Figures 22, 23, 24, and 25. This is due to contributions of several unmodeled physical effects. These effects include boundary interactions and the treatment of self-gravity. The outer boundary allows outflow of gas (and angular momentum) from the domain, but we keep track of this, as described in equation (26). Unaccounted for effects include the back-reaction of the inner, reflecting boundary condition surrounding the donor star core pushing on the gas. One might imagine this adds angular momentum to the system in the case of orbiting particles in gas. Indeed, with the largest central boundary in Figure 22, the total angular momentum increases slightly. The other physical effect, which we suggest is most significant, is the unmodeled gas self-gravitational interaction. We treat gas self-gravity statically, by applying the potential of the initial, spherical, undisturbed star, as described by equation (10). This becomes less accurate as mass is removed from the donor and distributes in increasingly complex configurations around the binary. That said, the overall magnitude of the non-conservation (about 2% given a factor of similar to 2 change in component momenta) justifies this approximation for the purposes of our analysis.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description