Twin Binaries: Studies of Stability, Mass Transfer, and Coalescence
James C. Lombardi Jr., William Holtzman, Katherine L. Dooley, Kyle Gearity, Vassiliki Kalogera, Frederic A. Rasio
aa r X i v : . [ a s t r o - ph . S R ] J un Draft version September 4, 2018
Preprint typeset using L A TEX style emulateapj v. 04/20/08
TWIN BINARIES: STUDIES OF STABILITY, MASS TRANSFER, AND COALESCENCE
J. C. Lombardi, Jr. , W. Holtzman , K. L. Dooley , K. Gearity , V. Kalogera , , F. A. Rasio , Draft version September 4, 2018
ABSTRACTMotivated by suggestions that binaries with almost equal-mass components (“twins”) play an im-portant role in the formation of double neutron stars and may be rather abundant among binaries,we study the stability of synchronized close and contact binaries with identical components in circularorbits. In particular, we investigate the dependency of the innermost stable circular orbit on the coremass, and we study the coalescence of the binary that occurs at smaller separations. For twin binariescomposed of convective main-sequence stars, subgiants, or giants with low mass cores ( M c . . M ,where M is the mass of a component), a secular instability is reached during the contact phase, accom-panied by a dynamical mass transfer instability at the same or at a slightly smaller orbital separation.Binaries that come inside this instability limit transfer mass gradually from one component to theother and then coalesce quickly as mass is lost through the outer Lagrangian points. For twin giantbinaries with moderate to massive cores ( M c & . M ), we find that stable contact configurationsexist at all separations down to the Roche limit, when mass shedding through the outer Lagrangianpoints triggers a coalescence of the envelopes and leaves the cores orbiting in a central tight binary.In addition to the formation of binary neutron stars, we also discuss the implications of our resultsfor the production of planetary nebulae with double degenerate central binaries. Subject headings: binaries: close — binaries: general — hydrodynamics — instabilities — methods:numerical — stars: general INTRODUCTION AND MOTIVATION1.1.
Formation of Binary Neutron Stars
The evolutionary history and formation of close bina-ries with two neutron stars similar to the Hulse-Taylorpulsar B1913+16 (Hulse & Taylor 1975) and the dou-ble pulsar J0737-3039 (Burgay et al. 2003) is a topicof intense current interest. Most recent studies of theknown double neutron stars focus on the stages goingback to the time of the second supernova explosion andthe formation of the youngest of the two neutron stars(Dewi & van den Heuvel 2004; Willems & Kalogera 2004;Willems et al. 2004; Piran & Shaviv 2005; Stairs et al.2006; Wang, Lai, & Han 2006; Wong et al. 2010, andreferences therein). Although these studies provide veryinteresting constraints on the properties of the stellarprogenitor of the second neutron star, they do not probethe earlier evolutionary history. That part remains un-certain and more difficult to constrain empirically basedon the measured properties of observed systems.Since the discovery of the Hulse-Taylor binary, the ori-gin of double neutron stars has been naturally connectedto the evolution of massive binaries, with stellar com-ponents that are massive enough to form two neutronstars at the end of their lifetime. Over the years, a qual-itative consensus of understanding for the formation ofdouble neutron stars developed (see, e.g., Bhattacharya& van den Heuvel 1991): massive binaries experience a Department of Physics, Allegheny College, Meadville, PA,16335 Department of Physics, University of Florida, Gainesville, FL32611 Department of Physics and Astronomy, Northwestern Univer-sity, 2145 Sheridan Rd., Evanston, IL 60208 Center for Interdisciplinary Exploration and Research in As-trophysics (CIERA), Northwestern University. phase of stable mass transfer when the primary over-flows its Roche lobe revealing its helium core; this coreends its lifetime in a supernova forming the first neu-tron star in the system; the binary becomes a high-massX-ray binary until the massive secondary fills its Rochelobe and the binary enters a dynamically unstable phaseof mass transfer leading to inspiral in a common enve-lope that engulfs the neutron star; during this phase theneutron star is thought to be spun up through recyclingand, if the binary avoids a merger in the inspiral, thehelium core of the secondary is revealed after the enve-lope ejection; this core explodes in a supernova, formingthe second neutron star in the system, and the doubleneutron star further evolves through orbital contractionand gravitational-wave emission. Variations of this evo-lutionary sequence have been shown to be realized bytheoretical binary population studies (e.g., Belczynski etal. 2002).Brown (1995) argued that the inspiral of the neutronstar during the common envelope phase in the standardmodel is problematic: the neutron star is expected to ex-perience hypercritical accretion (Chevalier 1993) at ratesmany orders of magnitude above the photon Eddingtonlimit through a neutrino-cooled accretion flow (see alsoFryer et al. 1996). Such a rapid accretion phase alongwith the adoption of a low maximum neutron star mass( ∼ ⊙ derived for a soft equation of state for neutronstar matter with kaon condensation) led Brown to con-clude that all neutron stars in common envelope phaseswill accrete enough matter to collapse into a black hole.Consequently, he argued that the “standard” evolution-ary sequence described above aborts the formation ofdouble neutron stars and instead leads to the formationof binaries with low-mass black holes and neutron starcompanions. We note however that, even with the same Lombardi et al.treatment of hypercritical accretion, a maximum neu-tron star mass of ∼ M ⊙ (corresponding to a more reg-ular stiff equation of state) does prevent a good fractionof neutron stars from being transformed into low-massblack holes (Belczynski et al. 2002).Brown (1995) also noted how the masses measured inknown double neutron stars are very close to being equal(Nice et al. 1996; Thorsett & Chakrabarty 1999; Stairset al. 2002; Weisberg & Taylor 2005; Jacoby et al. 2006;Kramer et al. 2006, and references therein). Motivatedby these two points, he proposed a different formationchannel for double neutron stars. Brown suggested thatdouble neutron stars form from massive binaries withcomponent masses that are within ∼
4% of one another.Consequently the red giant phases of the two compo-nents coincide in time and when mass transfer ensuesfrom the primary, both components have deep convectiveenvelopes and well developed helium cores, so that a dou-ble core phase develops, where the two helium cores orbitwithin the combined envelopes of the two massive stars.Provided that there is enough orbital energy, the com-mon envelope is ejected before the two cores merge, anda tight binary with two helium cores is formed. Thesetwo cores differ very little in mass and reach core col-lapse one very soon after the other ( ∼ yr based onhelium-star models, Habets 1986; Pols 1994), forming aclose double neutron star.The advantages of this hypothesized evolutionarychannel are (i) a neutron star never experiences com-mon envelope spiral-in and hypercritical accretion, and(ii) the two stellar components have so similar massesthat they naturally form neutron stars of almost equalmass as observed (Bethe & Brown 1998; Bethe, Brown,& Lee 2007). On the other hand, this channel (i) re-quires that mass transfer between the red giant progeni-tors will indeed lead to the inspiral of the two cores in acommon envelope and (ii) requires fine-tuning the condi-tions for recycling the first neutron star, as this spin upmust occur during the very short interval between thetwo supernovae through the stellar wind or possibly dur-ing the brief Roche-lobe overflow from the lower-masshelium star (Dewi et al. 2006). A potential additionaldisadvantage is that this channel is very restrictive inthat it requires progenitors that are at most only ∼ . .
5, Hogeveen 1992a,b,c); therefore the con-tribution of twins may not be as significant as implied bythe most recent observations. Quantitative modeling ofthe associated selection effects is required to derive more reliable statistical conclusions.Apart from uncertainties with the initial properties ofthe binary population, a number of physical processesrelated to these formation channels make it hard to as-sess their relative contribution to double neutron starformation. The physics of common envelopes and hyper-critical accretion, as well as of neutron star equations ofstate and their maximum mass, is not well understood(but see Lee, Park, & Brown 2007). Also, the develop-ment of a dynamical instability and subsequent commonenvelope phase with the inspiral of the two cores is as-sumed, but has not been investigated before in any de-tail. In this study, we attempt to understand better oneof the aspects related to the formation channel suggestedby Brown (1995): the fate of mass transfer between twostars of almost equal mass.1.2.
Binaries and Planetary Nebulae
Although our primary motivation is to investigatea formation channel for binary neutron stars, our re-sults are relevant to other scenarios as well. For exam-ple, through common envelope evolution, binaries couldtransform quite naturally into planetary nebulae (PNe)with a central close binary. The influence of binaries onthe production and morphology of PNe has received in-creased attention in recent years (e.g., Miszalski et al.2009; Jones et al. 2010), and observations to date haveidentified approximately 40 close binaries in the centersof PNe: see De Marco (2009) for a summary of both therelevant theory and observations. As lifetimes of PNeare less than only 10 years, an embedded central binaryhas undergone no significant evolution since the commonenvelope phase that presumably formed it.About a quarter of the known close binary systemswithin PNe are thought to be double degenerates, thatis, binaries in which the components are pre-white dwarfs(also known as extreme horizontal branch stars, hot sub-dwarfs, or subdwarf B and O stars) or white dwarfs.Such double degenerate binaries, with or without plan-etary nebulae, are particularly interesting. As possibleprogenitors for Type Ia supernovae, they are the desiredtargets of the ESO SN Ia Progenitor Survey (see Geieret al. 2011, and references therein).Interestingly, three of the five well studied double de-generate binaries within PNe (see Hillwig 2011, and ref-erences therein) may contain nearly equal mass compo-nents and therefore are possibly the progeny of twin bi-naries. In particular, models of a double hot subdwarfbinary with q ≈ q = 1. As we investigate in this paper, the coales-cence of a twin giant binary could lead to the formationof double degenerate core surrounded by a circumbinaryenvelope, a type of proto-PN.1.3. Theoretical Work on Close Binaries
Most of the classical theoretical work on close bina-ries was done in the limit of a self-gravitating incom-win Binaries 3pressible fluid (see Chandrasekhar 1969, and referencestherein). An essential result found in the incompress-ible case is that the hydrostatic equilibrium solutions forsufficiently close binaries can become globally unstable(Chandrasekhar 1975; Tassoul 1975). The classical an-alytic studies for binaries were extended to polytropesin the work of Lai et al. (1993a,b, 1994a,b,c). In theirapproach, the stars are modeled as self-gravitating com-pressible ellipsoids, and an energy variational principleis used to construct approximate equilibrium configura-tions and study their stability. These treatments, alongwith complementary numerical hydrodynamic calcula-tions (Rasio & Shapiro 1992, 1994, 1995), demonstratethat dynamical instabilities persist in the compressibleregime and can cause a binary to coalesce to form arapidly rotating spheroidal object. Such a dynamical in-stability can trigger a merger in just a few orbital periods(Rasio & Shapiro 1992) or an episode of mass transferthat lasts many orbits (Motl et al. 2002; D’Souza et al.2006; Frank 2008; Dan et al. 2009; Lor´en-Aguilar et al.2009).The evolution of a close binary system can also be af-fected by another type of global instability. It has beenreferred to by various names, such as the secular instabil-ity (Lai et al. 1993a,b, 1994a,b,c), tidal instability (Coun-selman 1973; Hut 1980), gravogyro instability (Hachisu& Eriguchi 1984), and Darwin instability (Levine et al.1993). Its physical origin is easy to understand (Lai et al.1993a, 1994b; Rasio 1994; Webbink 2006). There existsa minimum value of the total angular momentum J fora synchronized close binary. This is simply because thespin angular momentum, which increases as the separa-tion r decreases for a synchronized system, can becomecomparable to the orbital angular momentum for suffi-ciently small r . A system that reaches the minimum of J cannot evolve further by angular momentum loss andremain synchronized. Instead, the combined action oftidal forces and viscous dissipation will drive the systemout of synchronization and cause rapid orbital decay asangular momentum is continually transferred from theorbit to the spins. The orbital decay then proceeds ona timescale comparable to the synchronization time of astable binary.In this paper, we pay particular attention to the onsetof orbital instabilities, including those characterized bymass transfer, as well as the subsequent inspiral of thecores. We extend previous hydrodynamic studies of closebinary systems to cases that involve identical giants, thatis, to twin stars with dense stellar cores and extendedenvelopes. So that our results can be scaled to stars ofany size or mass, we approximate each star as a con-densed polytrope, namely, as a point mass surroundedby a uniform specific entropy fluid of adiabatic indexΓ = 5 /
3. Such a model is appropriate for a fully convec-tive monatomic ideal gas surrounding a compact core.We consider fractional core masses m c that cover theentire range of theoretical possibilities. Zero core massmodels correspond to normal (or “complete”) n = 1 . m c = 1, there is no massin the gaseous envelope, and the binary is simply thatof two point masses. By varying the core mass betweenthese extremes, we are able to study in a systematic way the full parameter space of twin binaries and determineunder what conditions mass transfer develops or an in-nermost stable circular orbit exists.For typical compositions, the subgiant phase beginswhen the core mass grows to ∼
10% of the total mass,the so-called Sch¨onberg-Chandrasekhar limit (Sch¨onberg& Chandrasekhar 1942). By the time the star reaches thebase of the red giant branch, the core mass has increasedby a few more percent of the total mass. Roughly speak-ing then, our models with m c . .
1, 0 . . m c . . m c & .
13 correspond to main sequence, subgiant,and giant stars, respectively. The assumption that thestellar envelope has constant specific entropy makes ourmodels most relevant to stars that are fully convective(as with low mass main sequence stars) or that have deepconvective envelopes (as in many red giants).The use of condensed polytropes as a model of redgiants has a rich history, including seminal work byChandrasekhar (1939), Osterbrock (1953), and H¨arm &Schwarzschild (1955). Mass transfer in close binary sys-tems has been modeled with the help of condensed poly-tropes as well: the response of the donor due to mass lossis considered, with its resulting contraction, or expan-sion, compared against that of its Roche lobe (Paczyi´nski& Sienkiewicz 1972; Hjellming & Webbink 1987). If atthe onset of mass transfer the star contracts less rapidlythan its Roche lobe (or expands more rapidly than it),then the mass exchange is dynamically unstable, andthe binary will evolve rapidly toward a new, often qual-itatively different, equilibrium. Hjellming & Webbink(1987) find that equal mass condensed polytrope con-tact binaries with fractional core masses m c . .
46 ex-perience stable mass transfer. More recently, Krumholz& Thompson (2007) have used condensed polytropes tostudy the formation of twin star systems.Although essential for a qualitative understanding ofmass transfer, such treatments of close binaries do, how-ever, make several simplifying approximations: most im-portantly, (i) the dynamics of the orbit and size of theRoche lobe are treated in the point mass approximation,(ii) the response of the binary components to mass lossor gain is modeled as if each star were spherical and inisolation, and (iii) mass that overflows a Roche lobe isconsidered to leave that star. These approximations arequite reasonable for semidetached binaries, but their va-lidity can be questioned for contact binaries. In suchcases, a common envelope persists in equilibrium out-side of the Roche lobes (the inner Lagrangian surface)so that the pressure and density on the Roche lobes arenon-zero: mass that overflows a Roche lobe is not neces-sarily transferred to the other star but rather can persistin equilibrium inside the outer Lagrangian surface. Aprimary goal of this paper is therefore to relax the ap-proximations of previous works by using accurate hydro-dynamical calculations to study contact binary systems.Our paper is organized as follows. In § § § NUMERICAL METHODS AND ASSUMPTIONS Lombardi et al.2.1.
The SPH code
To generate our models, we use a modified version ofthe SPH code originally developed by Rasio (1991) thathas been updated to include the variational equationsof motion derived in Gaburov et al. (2010). SPH is aLagrangian particle method, meaning that the fluid isrepresented by a finite number of fluid elements or “par-ticles.” Associated with each particle i are, for example,its position r i , velocity v i , and mass m i . Each particlealso carries a purely numerical smoothing length h i thatdetermines the local spatial resolution and is used in thecalculation of fluid properties such as acceleration anddensity. Details of our SPH code, such as the particularform of the artificial viscosity Π ij and smoothing kernel W ij implemented, are described in Gaburov et al. (2010).See Rasio & Lombardi (1999) and Rosswog (2009) for re-views of SPH.Because the gas in our initial stellar models is of con-stant specific entropy, we find it convenient to integratethe so-called entropic variable A i of each particle i . Theentropic variable is simply the proportionality constantin the polytropic equation of state p = Aρ Γ , where p is pressure and ρ is density. The entropic variable isso named because of its close connection to entropy:both quantities are conserved in reversible processes andstrictly increase otherwise. We therefore use dA i /dt = 0in the relaxations of our single star models and in thecalculations of our binary equilibrium sequences. Forour dynamical calculations of merger scenarios (see § § A i according to the discretized SPHversion of the first law of thermodynamics: dA i dt = Γ − ρ Γ − i X j m j Π ij ( v i − v j ) · ∇ i W ij ( h i ) . (1)To calculate the gravitational accelerations and po-tentials, we use direct summation on NVIDIA graphicscards, softening with the usual SPH kernel as in Hern-quist & Katz (1989). The use of such a softening withfinite extent (as opposed, for example, to Plummer soft-ening) increases the accuracy and stability of our SPHmodels, consistent with the studies of Athanassoula etal. (2000) and Dehnen (2001). The gravity of core pointsin our models is similarly softened, applying a constantsmoothing length comparable to the minimum smooth-ing length in the system.2.2. Single Star Models
In this section we present our procedure for modelingthe stars that are used in the binary simulations of § / m c . Each of the condensed poly-tropes in our family of models has total mass M = 1and radius R = 1. The unit system is completed bychoosing Newton’s gravitational constant G = 1. Con-densed polytrope models have not only the advantage ofreproducibility but also of scalability to any stellar massand radius: although our focus here is on stars massiveenough ultimately to yield a neutron star, the same cal-culations are also valid for low mass systems. In thelimiting scenarios of m c = 0 and m c = 1, we recover thewell-studied cases of a n = 3 / (Portegies Zwart et al. 2009): we evolve 10 and25 M ⊙ stars with initial helium abundance Y = 0 . Z = 0 .
02 until they obtain core massesof approximately m c = 0 .
16 and 0.27, respectively. Thiscorresponds to our initially 10 M ⊙ star being at the verytip of the red giant branch and the initially 25 M ⊙ starbeing on the upper part of the branch. The profiles of thesimple condensed polytrope models follow the same gen-eral trend as those of the red giants, indicating that oursimple models can indeed help provide an understandingof real red giant binaries.Table 1 presents the models used in this paper. Col-umn (1) gives the core mass m c , while column (2) lists thecorresponding E O value: this is the parameter E in theoriginal notation of Osterbrock (1953), as well as in thenotation of H¨arm & Schwarzschild (1955) and Hjellming& Webbink (1987). (We prefer to reserve the variable E for energy.) The value of E O controls the shape of thedensity profile: to lowest order near the surface ( r ≈ R ), ρ ( r ) ≈ (2 / / E O (1 − r/R ) / M/ (4 πR ). In practice,we adjust E O to achieve the desired core mass m c . TABLE 1Parent Star Characteristics m c E O (1) (2)0 45.48080.05 39.92500.1 35.24030.125 33.16610.15 31.24070.175 29.44610.2 27.76730.25 24.70720.3 21.97820.4 17.28610.5 13.35820.6 9.99040.7 7.04960.8 4.44430.9 2.10920.99 0.1984 We begin by making an SPH model of a single star http://muse.li win Binaries 5 Fig. 1.—
The solid curves show the quasi-analytic radial profilesof the pressure p and density ρ for m c = 0, 0.125, and 0.5 condensedpolytropes; for comparison, the dashed and dotted curves representred giants, as modeled by the TWIN stellar evolution code, withinitial masses respectively of 10 and 25 M ⊙ . For the condensedpolytropes, curves associated with larger m c are higher on the leftedge of the figure and lower on the right edge. Units are such that G = M = R = 1. in isolation. Unless stated otherwise, we use N =19938 SPH particles initially placed on a hexagonal closepacked lattice with a lattice spacing constant a =0 . h = 0 . M = 1 is precisely achieved. Theentropic variable A i of each SPH particle is set to thedesired polytropic constant K , which is determined fromthe mass-radius relation for n = 3 / R = (4 π ) − / G − KE / M − / . (2)After the initial parameters of the particles have beenassigned, we relax the SPH fluid into hydrostatic equilib-rium. This relaxation is effected by including an artificialviscosity contribution in the acceleration equation (with α = 1 and β = 2 in equation (A19) of Gaburov et al.2010), while still keeping A i constant for all particles. Inthis way, the entropy of the system is preserved while thesystem approaches equilibrium. The total energy typi-cally decreases by less than a percent in the process, in-dicating that our initial assignment of particle propertieswas indeed very close to an equilibrium state.This approach allows us to model the desired profiles Fig. 2.—
Properties of the model with core mass m c = 0 .
25 as afunction of radius, after relaxation for 500 time units. The framesin the left column compare calculated pressure p and density ρ pro-files of the star (dashed curves) against particle data from our SPHmodel (dots). The right column provides additional SPH particledata: individual SPH particle mass m , smoothing length h , num-ber of neighbors N N , and radial component of the hydrodynamicacceleration a hydro (upper data) and gravitational acceleration g (lower data). very accurately. An example is presented in Figure 2,where we plot desired profiles and SPH particle data forour relaxed m c = 0 .
25 star. Although the core and sur-face of the star of course cannot be resolved on the lengthscale of a smoothing length (typically 0.05 to 0.08 lengthunits), the thermodynamic profiles of the SPH modelnicely reproduce the quasi-analytic curves. Indeed, theSPH data in the left column of Figure 2 are difficult todistinguish from the desired pressure and density pro-files throughout most of the star. We also note that thehydrodynamic and gravitational accelerations are verynearly equal in magnitude and opposite in direction, asnecessary for hydrostatic equilibrium.2.3.
Binary Equilibrium Configurations
SPH Calculations
The ability of our code to model close and contact bi-nary systems for thousands of orbits or longer is pre-sented in Gaburov et al. (2010). Here we present ourmethods for modeling equilibrium sequences of twin bi-naries, that is, binaries that consist of two identical starsin synchronized orbit. First, we place identical relaxedstellar models along the x axis with their centers of massseparated by r . While the relaxation of the binary takesplace, the entropic variable particle values are held con-stant. The center of mass of the entire system is fixed atthe origin. In addition, the positions of the particles arecontinuously adjusted (by a simple uniform translationalong the binary axis), so that the separation betweenthe centers of mass equals the desired separation r .The orbit is chosen to occur in the xy plane. The an- Lombardi et al.gular velocity Ω orb defining the corotating frame is up-dated at every timestep so that the centrifugal and iner-tial accelerations acting on the fluid cancel. In particu-lar, we wish to find configurations in which Ω ( x i ˆx + y i ˆy) = − ( ˙ v x,i ˆx + ˙ v y,i ˆy), where the velocity derivativeson the right hand side are components of accelerationin the inertial frame. By taking the dot product with m i ( x i ˆx + y i ˆy) and then summing over all particles, weobtain Ω = − P i m i ( x i ˙ v x,i + y i ˙ v y,i ) P i m i ( x i + y i ) . (3)We also include a drag force that opposes the veloc-ity and provides a contribution to the acceleration of − v i /t relax . We use t relax = 3, approximately the funda-mental period of oscillation for our parent models. Wedo not include any artificial viscosity contribution whenfinding binary equilibrium configurations.In order to find equilibrium configurations for a pre-cisely equal mass binary, even for configurations unsta-ble to mass transfer, we enforce a symmetry in particleproperties: for each particle i in star 1, there is a part-ner particle j in star 2 at x j = − x i and y j = − y i withvelocity components v x,j = − v x,i , v y,j = − v y,i and withacceleration components ˙ v x,j = − ˙ v x,i , ˙ v y,j = − ˙ v y,i . Allother properties are identical for any such pair of parti-cles.The separation r between the centers of mass can beallowed to drift slowly so that an equilibrium sequenceis constructed: a so-called “scanning run.” In practice,we start runs that will scan over separations by hold-ing the centers of mass fixed at an initial separation r (0) for 40 time units, allowing the system to approacha tidally bulged equilibrium configuration. At an addi-tional amount of time t , the separation is set accordingto r ( t ) = r (0) [ r ( t scan ) /r (0)] t/t scan . This form for r ( t ) al-lows the change in r to occur at a decreasing rate as thestars approach and interact more strongly, although theexact form is not critical to our results. We typically use t scan = 300, r (0) = 3 .
3, and r ( t scan ) = 2 . Data Reduction Methods
Once an SPH binary calculation has completed, weanalyze the system at various separations along the se-quence. To this end, a useful quantity to consider is theeffective potential, calculated asΦ e ( x, y, z ) = Φ( x, y, z ) −
12 Ω ( x + y ) , (4)where Φ is the gravitational potential, the coordinate y measures perpendicular to the binary axis in the or-bital plane, and z measures parallel to the rotation axis.Along the binary axis ( y = z = 0), the effective potentialhas a local maximum Φ ( i ) e at x = 0 (the inner Lagrangianpoint) and global maxima Φ ( o ) e at | x | = x o (the outer La-grangian points). There are two minima at | x | = x c , cor-responding to the cores of the two components. In equi-librium, the fluid will fill up the effective potential wellto some maximum, constant level Φ ( s ) e . Borrowing the Note that in general r = 2 x c , because we define the binaryseparation r as the distance between the centers of mass of thetwo components. terminology from models of W UMa binaries (Rucinski1992, and references therein), we follow Rasio & Shapiro(1995) and define the degree of contact η as η ≡ Φ ( s ) e − Φ ( i ) e Φ ( o ) e − Φ ( i ) e . (5)Clearly, we have η < < η <
1, the effec-tive potential of the fluid near x = 0 does exceed thebarrier, and the system is classified as a contact binary.For η >
1, the envelopes overflow beyond the outer La-grangian surface, and no dynamical equilibrium configu-ration can be achieved; that is, the system has exceededthe Roche limit. Calculations in which we slowly scanto smaller separations can therefore determine positionof first contact ( η = 0), the secular stability limit (atthe minimum energy and angular momentum along thesequence), and the Roche limit ( η = 1).It is important to realize that the equilibrium sequenceof a twin binary passes smoothly from detached to con-tact configurations as the separation r decreases. Thisis in contrast to all binary equilibrium sequences withmass ratio q = 1, which terminate at a Roche limit cor-responding to the onset of mass transfer through the in-ner Lagrangian point (that is, once one of the binarycomponents overflows its Roche lobe). For twin binaries,however, the Roche limit, which we still define as the lastequilibrium configuration along a sequence with decreas-ing r , corresponds to the onset of mass shedding throughthe outer Lagrangian points: as an example, note thatseveral particles have been shed to the far left and farright of the r = 2 .
22 frame in Fig. 3.We estimate Φ ( s ) e from our SPH models by finding themaximum effective potential of the points along the x -axis that are within one smoothing length of the cen-ter of an SPH particle. Thus, even if the centers of allSPH particles are within the outer Lagrangian surface,we may still consider the system as having reached theRoche limit when some smoothing kernels extend sub-stantially beyond the outer Lagrangian surface. Such anestimate accounts for the fact that an SPH particle isnot a point mass but instead represents a parcel of fluidwith a density profile described by the smoothing kernel.Our means of estimating Φ ( s ) e allow our critical separa-tion results to converge quickly to a steady value as theresolution is increased up to the resolution presented inthis paper (see § Dynamical Calculations
We generate initial conditions for our dynamical runsby taking a configuration at the desired separation r froma scanning run and then relaxing for an additional 200time units. If a particle escapes past an outer Lagrangianpoint during this time interval, then we end the relax-ation stage and begin following the dynamics immedi-ately. During dynamical calculations, we include no dragforce and move the particles according to their velocitiesin the usual way (for details, see Gaburov et al. 2010).Particles are again treated independently so that masstransfer events can be followed; that is, unlike the scanswin Binaries 7described in § orb calculated once at the beginning of the dynamical evolu-tion and thereafter held constant when applying Coriolisand centrifugal forces. RESULTSUsing the methods described above, we construct twinbinary sequences for 16 different core masses m c listed inTable 1 and covering the range from 0 to 0.99. In § § r then allow us to test the dynamical stability of thecontact configurations and to follow any mass transferand the merger in unstable systems.3.1. Equilibrium Sequences
Representative snapshots along the m c = 0 . xy plane) and in terms of the effectivepotential Φ e . The thick solid curves in Figure 3 are thesurfaces of constant effective potential Φ e that mark theinner and outer Lagrangian surfaces. For fixed x , Φ e isminimum on the binary axis ( y = z = 0), and this min-imum value is given as a dashed curve in Figure 3. Inhydrostatic equilibrium, the fluid fills up to a constantlevel Φ ( s ) e that is independent of x .Referring to Figure 3 we see that at the initial separa-tion in our scan, r = 3 .
30, the system is tidally bulgedand the binary is detached: the fluid does not extend outto the inner Lagrangian surface (also known as the Rochelobe). At a separation r = 2 .
71, the binary stars fill theinner Lagrangian surface and make first contact throughthe L1 point, located at the origin. Once the separationdecreases to r = 2 .
42, the binary reaches the secular in-stability limit, as marked by a minimum in the energy E and angular momentum J along this equilibrium se-quence (see below). The critical separation r = 2 . § § r = 2 .
24, thefluid extends out to the outer Lagrangian surface, mark-ing the Roche limit. At even smaller separations, forexample at the separation r = 2 .
22 shown in the figure,no equilibrium configurations exist and the stars shedmass through the outer Lagrangian surface near the L2point. The variation of critical effective potentials andthe degree of contact η along this m c = 0 . η on theseparation r .From Figure 5, we note that as m c increases toward 1,our results approach the Keplerian solution of two orbit-ing point masses. As expected, deviations from the point Fig. 3.—
Sequence of binary equilibrium configurations for twoidentical condensed polytropes of core mass m c = 0 .
1. Projectionsonto the orbital plane (the xy plane) are shown at six differentbinary separations for those SPH particles with | z | < .
06. Thethick solid curves represent two surfaces of constant effective po-tential Φ e (see eq. [4]): namely, the inner and outer Lagrangiansurfaces passing through the points L1 and L2. Shown beneatheach configuration are corresponding projections onto the ( x, Φ e )plane for the same particles. The dashed curves give the variationof Φ e along the binary axis ( y = z = 0). Contact configurationsare obtained when the binary separation r . .
71 (in units wherean isolated binary component has radius R = 1). For r . . mass result increase as a given binary becomes moredeeply in contact or as we consider a binary associatedwith a smaller core mass. From the bottom frame of Fig-ure 5, we note that smaller core masses (correspondingto stars with less centrally concentrated density profiles)have a smaller orbital period at any given separation. For m c <
1, tidal interactions between the two stars makethe effective potential stronger than 1 /r and shorten therotation period compared to a point mass system. For m c = 0 .
1, for example, the deviation of the orbital pe-riod from the point mass result is approximately 1% atfirst contact, 2% at the secular instability limit, and 3%at the Roche limit.As in Rasio & Shapiro (1995), we determine the secu-lar stability limit along the equilibrium sequence by lo-cating where both the total energy E and total angularmomentum J are minimum in curves such as those ofFigure 5. Our numerical results provide an accurate de-termination of this point for a close binary system, as theseparate minima in E and J coincide to high numericalaccuracy. This is in accord with the general propertythat dE = Ω orb dJ along any sequence of uniformly ro-tating fluid equilibria (Ostriker & Gunn 1969).For the m c = 0 binary, secular instability occurs soonafter contact along this sequence and therefore stable,long-lived equilibrium configurations can exist only inshallow contact, η . .
2. In contrast, the sequences withnon-zero core masses permit much deeper contact before Lombardi et al.
Fig. 4.—
Variation of critical effective potentials along the equi-librium sequence presented in Fig. 3. Values of the effective poten-tial at the outer Lagrangian surface (solid curve o), the inner La-grangian surface (solid curve i), and the fluid surface (long dashedcurve s) are shown as a function of binary separation r in the topframe. The degree of contact η (eq. [5]) is shown in the bottomframe as the long dashed curve. The short dashed curves give thepositions of first contact ( η = 0) and of the Roche limit ( η = 1). the secular instability is reached. For example, a binarywith core masses of 0 .
125 does not reach the secular in-stability until nearly η = 0 . . For core masses m c & . η = 0), the secularinstability limit (where E and J are minima), and at theRoche limit ( η = 1) for our sequences of various coremass m c . Additional runs at varying resolution indicatethat the results in our tables have converged to within ∼
1% (e.g. see Fig. 6).These results for critical separations are summarizedin Figure 7. Due to tidal effects, the volume of each staris typically 1-2% larger at first contact than it is for thatstar in isolation. The separation r fc at first contact isonly weakly dependent on m c , being within 2% of 2.7for any core mass. For comparison, we note that thestandard simple treatment of twin binaries would implya first contact separation of 1 / . .
63, where thefactor 0.3799 comes from numerical integration of Rochelobe volumes around identical point masses (e.g., Eggle-ton 1983) and any change in the volumes of the stars dueto tidal effects is neglected. We see therefore that finitesize effects act to increase the separation of first contact.All three of the critical separations considered (firstcontact, secular instability, and Roche limit) tend to de-crease as the core mass increases. It is straightforward tofind fitted formulas for the critical separations that areaccurate to within ∼
1% for any core mass m c : for firstcontact r fc ≈ .
66 + 0 . − m c ) (6) Fig. 5.—
Variation of the system energy E (relative to the totalself-energy E ∞ of the binary components at infinity), angular mo-mentum J , and orbital period P along the equilibrium sequencefor twin binaries with m c = 0, 0.1, 0.25, 0.5, and 1. In the E and J frames, higher curves correspond to smaller core mass, while inthe P frame higher curves correspond to larger core mass. Theorbital period P is normalized to the analytic point mass result P Kepler = 2 / πr / . The curves are from SPH scans of the equi-librium sequence, except for the m c = 1 curve which is the analyticresult for two point masses. The m c = 0 and 0 . E and J , marking the position of the secular instabil-ity limit. The curves from the SPH scans terminate at the Rochelimit, where mass shedding through the outer Lagrangian pointcommences. Critical points along our SPH scans are presented inTables 2, 3, and 4. The individual data points (open circles) inthis figure result from relaxing for an additional 200 time units atthe given separation r . The agreement between these points andtheir corresponding scan helps to confirm that our scans are indeedproducing equilibrium sequences. TABLE 2FIRST CONTACT ALONG THE EQUILIBRIUMSEQUENCES OF TWIN BINARIES a m c r P E − E ∞ J a Units are defined such that G = M = R = 1, m c is the coremass of each component, r is the binary separation, η is the degreeof contact (eq. [5]), P is the orbital period, and E − E ∞ and J are the orbital energy and angular momentum, respectively; theenergy E ∞ is the total equilibrium energy at infinite separation(that is, twice the energy of a single component in isolation). win Binaries 9 TABLE 3SECULAR INSTABILITY ALONG THE EQUILIBRIUMSEQUENCES OF TWIN BINARIES b m c r η P E − E ∞ J b In twin binary sequences with core masses m c & . , the Rochelimit is reached before the secular limit. Units and column headingsare as in Table 2, footnote a; the degree of contact η is defined byeq. 5. TABLE 4ROCHE LIMIT c ALONG THE EQUILIBRIUM SEQUENCESOF TWIN BINARIES d m c r P E − E ∞ J c Defined as the equilibrium configuration with the minimum bi-nary separation. d Unit and column definitions are identical to those in Table 2,footnote a.
Fig. 6.—
Critical separation r , energy E − E ∞ , angular mo-mentum J , and orbital period P at first contact (squares), secularinstability (triangles), and the Roche limit (circles) versus totalparticle number N for several m c = 0 . N & . Most of the binary calcu-lations in this paper employ N ≈ × particles. Fig. 7.—
Separation at first contact (squares), the onset of secularinstability (stars), and the Roche limit (circles) versus core mass m c . The data points are determined from SPH scans of equilibriumsequences, as reported in Tables 2, 3, and 4. The lines are simplyto help guide the eye. and for the secular instability limit r sec ≈ . − m c . (7)We note that the 1 − m c in equation (6) equals the enve-lope mass (in units where the total stellar mass M = 1).We give our fit for the Roche limit separations in thenext subsection, where we can determine these data withslightly better accuracy. We have not fit for the slightincrease in the first contact data as the core mass is in-creased from m c = 0 . m c = 0 .
99 single star equilibriummodel settling to a radius a few percent larger than 1.The critical core mass m c ≈ .
15, for which the secularinstability and Roche limits coincide, can be determinedgraphically from Figure 7 by extrapolating the line con-necting the secular instability data down to the Rochelimit curve. This intersection is important because it im-plies that all of our equal mass binaries with m c & . m c . .
15, corresponding primarilyto main-sequence stars and subgiants, reach the secularinstability limit at a larger separation than that of theRoche limit. As we will see in the next subsection, thesecular instability limit is usually accompanied by dy-namical mass transfer at the same or a slightly smallerseparation. Therefore, main sequence and subgiant twinbinaries, as contrasted to most red giant twins, will startcoalescing (a) when in more shallow contact and (b)through mass transfer across the inner Lagrangian point.3.2.
Dynamical Integrations
Fig. 8.—
Results of dynamical integrations: stable binaries (opensquares), dynamically unstable binaries (asterisks), and configura-tions with no equilibrium (filled squares). The curves represent thecritical separations as determined by scanning runs, as in Fig. 7.The agreement between the results of the equilibrium sequencesand dynamical calculations is excellent.
We now study the stability of binary configurationswith fully dynamical SPH integrations (see § m c andinitial r values. We find the Roche limit to be verynearly at the separations determined from the equilib-rium scans, and the additional relaxation that we per-form before beginning a dynamical calculation allows usto determine these separations even more accurately. Inaddition, we find that most systems that are secularlyunstable are also dynamically unstable to mass transferand then merger, as discussed below.The time evolution of the separation of m c = 0 . r sec = 2 . § r & r sec are clearly dynamically stable,while those with r . r sec are dynamically unstable. Thatis, the secular and dynamical stability limits coincide orvery nearly coincide, at least at this core mass.The bottom plot of Figure 9 shows the dynamical evo-lution of two cases that straddle the instability limit. Forthe r = 2 .
56 case, an epicyclic period of 350 time unitsis clearly evident, much larger than the orbital periodof 17.9 time units. The large difference in these periodsis an indication of how close the system is to an insta-bility limit (Rasio & Shapiro 1994). Indeed, if r wereprecisely at the dynamical stability limit, the period ofsmall epicyclic oscillations would formally be infinite. Visualizations of selected simulations are available athttp://webpub.allegheny.edu/employee/j/jalombar/movies/.
Fig. 9.—
Top plot : Separation r versus time in several differentdynamical SPH calculations for the core mass m c = 0 .
05 twinbinaries. The initial separations are r = 2 .
64, 2.62, 2.56, 2.54,2.40, 2.30 and 2.28.
Bottom plot:
Zoomed in view of the r = 2 . r = 2 . Figure 10 presents projected particle positions (top)and column density plots (bottom) at six different timesin the m c = 0 . r = 2 .
54 dynamical calculation, acase just inside the instability limit. Colors in the parti-cle plot are used to indicate from which component theparticles originated. The coordinate system used here ro-tates counterclockwise with a period of 17.64 time units,which equals the orbital period of this binary at earlytimes, before significant mass transfer has occurred. Theinstability initially manifests itself in the form of a nar-win Binaries 11row arm of gas that begins in the outer layers of one star,gradually flows across the neck surrounding the inner La-grangian point, and then creeps around the other star(see the t = 746 and 771 particle plot frames). The masstransfer drives the binary components closer, trigger-ing the excretion of mass through the outer Lagrangianpoints ( t = 771) and accelerating the inspiral of the twocores. By t = 784, the mass transfered from one starhas completely engulfed the other. At later times, themerger product approaches a rapidly rotating, axisym-metric configuration centered on the two cores orbitingin a tight binary (see the t = 792 and 961 frames).Figure 11 shows the evolution of binary separation for m c = 0 . § r = 2 .
22 (or larger)orbits stably. In contrast, a binary at r = 2 .
20 grad-ually loses mass through the outer Lagrangian points,triggering a stage of rapid coalescence. The period ofmass loss persists for several orbital periods: each of theoscillations superposed on the decreasing r = 2 .
20 curvein the bottom plot of Fig.11 corresponds to one orbitalperiod. These dynamical calculations indicate that theRoche limit for a m c = 0 . r = 2 .
2, in excellent agreement with the r = 2 . m c = 0 . r = 2 .
20 dynamical calculation, a case just inside theRoche limit. The coordinate system used here rotatescounterclockwise with a period of 14.22 time units, whichequals the orbital period of this binary at early times.Gas is excreted almost immediately, with each parcelof gas carrying a specific angular momentum essentiallyequal to that of the outer Lagrangian points. In contrastto mergers with m c . .
15, the excreted gas originatesequally from both binary components and flows past theouter Lagrangian points symmetrically. As the outer La-grangian points are the outermost positions at which gascan be in static equilibrium, they are also the positions oflargest possible specific angular momentum in rigidly ro-tating equilibrium twin binaries. Consequently, the massloss necessarily decreases the average angular momentumper unit mass of the gas remaining within the outer La-grangian surface, causing the binary components, alongwith their cores, to inspiral: see the appendix of Webbink(1976) for a rigorous analysis of the angular momentumbudget during mass excretion. As the components getcloser, the excretion rate increases, and in addition theresulting arms become more tightly wound (compare the t = 202 frame to later ones). By t = 233, the central re-gions of the binary components have effectively merged.At later times, the merger product approaches a rapidlyrotating, axisymmetric configuration (see Fig. 13).Figure 14 shows energies versus time for the m c = 0 . r = 2 .
20 calculation. The kinetic energy T graduallyincreases as the binary components inspiral, until thecores approach closely at t ≈ U and increase in the gravitationalpotential energy W . The rapid variations in T and W at late times are due to the eccentric orbit of the centraldouble core. The total energy E is well conserved in this simulation, varying by only 0.4% from its minimum tomaximum values over the entire time interval shown inFigure 14. Energy conservation in other runs is typicallyat least this good and often even much better. In oursimulations, most of the small non-conservation in energyoccurs at late times once the core particles have entereda tight orbit.In all our merger simulations of twin binaries with frac-tional core masses of 0.15 or larger, each star loses massthrough an outer Lagrangian point. Most of this massultimately ends up in a circumbinary envelope gravita-tionally bound to the central binary of cores. The massejected varies from ∼ .
5% (for m c = 0 .
15) up to ∼ m c = 0 .
8) of the total system mass. A trend ev-ident from the simulations is that more massive coresultimately remain more widely separated, after inspiral-ing, than less massive cores. The top frame of Figure 15provides a closer look at how the separation of the coresevolves in several dynamical simulations that begin justinside the Roche limit. For very large fractional coremasses ( m c & . . . m c . . . . m c . .
5, which cor-responds to most red giants in nature, the cores rapidlyinspiral to separations less than 0.1. Although our codehas no mechanism for merging the cores, we do not ex-pect such an effect to be relevant here: the size of a corerelative to the stellar radius is typically only ∼ − orless for a giant, so that a tremendous amount of angularmomentum would have to be removed from the doublecore before they could merge.The bottom frame of Figure 15 concerns fractional coremasses less than 0.15, namely those cases that reach thesecular instability limit before the Roche limit. Thecores in these merger simulations inspiral to a separa-tion . .
01, considerably less than the spatial resolution.Whether or not the cores would merge in such circum-stances will likely depend on the details of the parent starstructure, with simulations that resolve the cores neces-sary to address the issue fully. We find that less than0.5% of the system mass is ejected whenever the mergeris triggered by mass transfer.Figure 16 helps to quantify the entropy evolution dur-ing mergers by plotting the mass-average < ln A > overtime for several cases. We note that, for our polytropicequation of state and gas of uniform composition, thespecific entropy of a parcel of gas is proportional to ln A plus a constant. From the curves of Figure 16, we seethat the change in entropy per unit mass tends to belarger, and develops on a longer timescale, for cases in-volving larger core masses. We also note that the entropyis still gradually increasing even at the end of our simula-tions, due to the influence of the central binary embeddedwithin the system.Like the critical separations for first contact and sec-ular instability, the Roche limit separation tends to de-crease as the core mass increases. A fitted formula con-sistent with our dynamical integrations to within ∼ Fig. 10.—
Top frames:
Projection of SPH particles onto the orbital plane at six different times in the merger of a m c = 0 .
05 binary withinitial separation r = 2 .
54. Particles are colored according to the star in which they originated.
Bottom frames:
Column density plots atthe same six moments, with the asterisks representing the positions of the compact cores. for any core mass m c is r Roche ≈ .
11 + 0 . − m c ) . (8)Because the degree of contact η varies nearly linearlywith separation r (for two examples, see Fig. 5 of this pa-per and Fig. 2 of Rasio & Shapiro (1995)), this formula,along with others from § η at the secular instability limit: η sec ≈ ( r fc − r sec ) / ( r fc − r Roche ). DISCUSSION AND FUTURE WORK4.1.
Summary of Main Results
We have determined equilibrium sequences and per-formed dynamical calculations of twin binaries, focusingprimarily on configurations in which the stars are in con-tact. Our equilibrium sequences of § m c . For m c . .
15, the innermost stableorbit occurs at the secular instability limit (marked bya minimum in energy and angular momentum along theequilibrium sequence); for m c & .
15, the innermost sta-ble orbit occurs at the Roche limit (defined as the min-imum separation for which an equilibrium configurationexists). Our dynamical calculations of § r on the vertical axisis scaled to the unperturbed stellar radius R , while thecore mass m c on the horizontal mass is scaled to the stel-lar mass M . Thus, as the components in a twin binaryexpand and gradually increase their core masses due tostellar evolution, the corresponding position in the pa-rameter space of Figure 17 will shift down and slightlyto the right. When this position drops below the topcurve, which marks first contact, the binary enters thewin Binaries 13 Fig. 11.—
Like Fig. 9, but for twin binaries with core mass m c = 0 .
2. At this core mass, there is no secular instability limit:instead, the system become unstable only once it reaches the Rochelimit at r ≈ .
2. The initial separations shown in the top plot are r = 2 .
60, 2.40, 2.22, 2.20, while a zoomed in view of the r = 2 . stable contact phase. When the position drops into ei-ther the unstable or no equilibrium portions of parameterspace, the components merge.The curves shown in Figure 17 are given by equations(6), (7), and (8); population modelers can use these fittedformulae in treatments of twin binaries. Consider, forexample, such a binary with a given orbital separation a . The stellar evolution of each component gives thetime dependence of the stellar radius R , the stellar mass M , and the core mass M c . The dimensionless separation r = a/R and the fractional core mass m c = M c /M are thus known functions of time, and the evolutionary trackcan be placed in the theoretical r versus m c diagram(Figure 17) to determine the final fate of the system.We note that the volume of parameter space where realbinaries would ultimately end up in the stable contactregion, without crossing the instability limit or the Rochelimit, is small. Such a situation would require a finetuning of the initial semimajor axis a . For example, astar with an initial mass of 8 M ⊙ will expand to R ≈ R ⊙ and reach a fractional core mass of m c ≈ . r fc ≈ . r Roche ≈ .
2. Thus, a twinbinary composed of such stars will remain detached if a & r fc R ≈ R ⊙ and will ultimately surpass the Rochelimit if a . r Roche R ≈ R ⊙ . Only if 800 R ⊙ . a . R ⊙ will the binary reach the contact phase withoutthe cores also inspiraling to form a tight binary.The dynamic simulations of § m c = 0 .
125 andless, we find a dynamical instability to exist at or slightlyinside the secular instability limit. This dynamical in-stability is a global instability of the equilibrium state,triggered by small numerical noise and characterized bya growing asymmetric mode. Binaries that come in-side this instability limit first transfer mass graduallyfrom one component to the other and eventually coa-lesce quickly as mass is lost through the outer Lagrangianpoints. The cores are left in a tight binary surroundedby a circumbinary envelope.4.2.
Comparison with Other Works
The merger of our small core mass twin binaries pro-ceeds in a fashion qualitatively similar to that of theQ1.3 model of D’Souza et al. (2006). In that model,the binary is composed of purely polytropic components( m c = 0) with mass ratio (donor to accretor) q = 1 .
3. Inboth our simulations and theirs, the dynamical instabil-ity manifests itself as a gradually developing mass trans-fer flow, followed by excretion of gas through the outerLagrangian points and merger of the stellar envelopes.One important difference, however, is that the instabil-ity in our m c = 0 case does not develop until the starshave reached a degree of contact η ≈ .
17, whereas theinstability is present in the Q1.3 model while the binaryis still semidetached. This difference highlights the stabi-lizing influence of the common envelope in twin binaries,even though the instability still exists even for q = 1.We can directly compare our results for the limitingcase m c = 0 with those of Rasio & Shapiro (1995).The agreement in the first contact, secular instability,and Roche limit separations is excellent (better than1%). The computational resources of the time, how-ever, limited Rasio & Shapiro (1995) to follow up toonly ∼ r ≈ .
45, well insidethe secular instability limit. In contrast, by following thedynamical evolution of m c = 0 twin binaries for up to ∼
100 orbits, we find that the dynamical instability limit4 Lombardi et al.
Fig. 12.—
Like Fig. 10, but for the m c = 0 . r = 2 .
20, just inside the Roche limit. actually coincides, or nearly coincides, with the secularinstability limit at r ≈ . Relevance to Binary Neutron Stars
We find that twin binaries with m c & .
15 exist stablyat separations all the way down to the Roche limit, wheremass is then excreted symmetrically through the outerLagrangian points. This excretion carries away angularmomentum and causes the stars, along with their cores,to inspiral on a dynamical timescale. For core masses m c . .
9, the cores inspiral to a final separation thatis a fraction of the stellar radius. Indeed for m c . . ∼ .
1, our dynamic simulations that lead to cores in a tight binary can provide only an upperlimit on the separation at the end of their inspiral. Oursimulations of the m c = 0 .
15 case, for example, placethis upper limit at ∼ .
03 times the stellar radius (seeFig. 15). Thus, a binary composed of twin M = 10 M ⊙ , R = 200 R ⊙ , M c = 1 . M ⊙ red giants would have theircores spiral to a separation of less than 0 . R ≈ R ⊙ .The gradual transfer of energy to the circumbinary en-velope could easily decrease the separation of the coresby a factor of ∼ . M ⊙ point masses separated byup to ∼ R ⊙ to contact in less than a Hubble time, suchsystems could have evolved to become arbitrarily tightby the present time. We therefore believe that such dou-ble cores are indeed excellent candidates for binary NSprogenitors, as proposed in the Brown scenario.4.4. Relevance to Planetary Nebulae
Given the relatively short timescale covered by hydro-dynamic simulations such as ours, the circumbinary en-win Binaries 15
Fig. 13.—
The state of the simulation presented in Fig. 12 at alate time ( t = 422), with the column density projected in both the xy and xz planes. Fig. 14.—
Energies versus time t for the simulation presented inFig. 12. From the bottom curve to the top one, we show gravi-tational energy W , total energy E , kinetic energy T , and internalenergy U . velopes at the end of our dynamical simulations are stilloptically thick. Nevertheless, it is natural to think ofthe final states of our merger calculations as a type ofproto-PNe, specifically, as the immediate precursors ofPNe with equal mass central binaries. Future work coulduse the end state of hydrodynamic calculations as initialconditions in models of PN evolution, with particular at-tention paid to any transition to the optically thin regime(revealing the central binary) and to the morphology ofthe gas distribution. Fig. 15.—
The separation r c between the cores as a functionof time t for binaries that begin just inside the Roche limit (topframe) and those that begin just inside the secular limit (bottomframe). Each curve is labeled by the factional core mass m c . Thehorizontal dashed line indicates the minimum separation for whichthe gravitational interaction of the two cores is unsoftened in oursimulations. To examine whether the dynamical calculations of thispaper yield final states that are consistent with the obser-vations of the q ≈ M ⊙ components. We use theTWIN stellar evolution code to evolve each star in isola-tion. The calculation shows that the star expands to only R ≈ R ⊙ as it ascends the giant branch and to nearly500 R ⊙ on the asymptotic giant branch. Therefore, forany reasonable distribution of initial orbital separations,twin binaries made of such stars would much more often6 Lombardi et al. Fig. 16.—
A measure of entropy evolution due to shock heatingduring several mergers. The natural logarithm of the entropic vari-able A is averaged by mass over the gas of the system and plottedversus a shifted time coordinate ∆ t , chosen such that the inflectionpoint in each curve (when the entropy is increasing most rapidly)occurs at ∆ t = 0. Each curve is labeled by the factional core mass m c : as m c increases, so do both the initial and final values of < ln A > . For each core mass, we present our simulation of largestinitial separation that results in a merger, corresponding in Fig. 8to the uppermost asterisk (for a given m c < .
15) or filled square(for a given m c ≥ . A equals the K in equation (2), as calculated from theappropriate E O in Table 1. The entropic variable A is in units of GM / R . come into contact while on the asymptotic giant branchthan when on the red giant branch, consistent with thefact that most or all of the well-observed double degen-erate stars in PNe seem to have masses too large to beHe white dwarfs.For the sake of discussion, consider such a twin binarythat reaches the Roche limit when R = 200 R ⊙ . Ac-cording to the stellar evolution calculation, an isolatedstar at that time would have a carbon-oxygen core ofmass M c = 0 . M ⊙ and, accounting for stellar winds,a total mass M = 2 . M ⊙ . At this fractional core mass m c = 0 . / . .
20, the dimensionless Roche separa-tion (from the results of this paper) is r Roche ≈ .
2, andthus the initial semimajor axis of this binary would havebeen a = r Roche R ≈ R ⊙ . From the m c = 0 . ∼ . R = 10 R ⊙ on the final separation of thetwo cores, which, from Kepler’s third law, corresponds toan upper limit on the orbital period of 3 days. If instead R = 100 R ⊙ at the time of merger, then a similar calcula-tion gives a core mass M c = 0 . M ⊙ and an upper limiton the orbital period of 0.9 days.For comparison, the central binaries in NGC 6026,Abell 41, and Hen 2-428 have comparable masses to thesecore masses and orbital periods in roughly the 0.2 to 0.5day range. We conclude that the results of our dynami-cal simulations, when applied to intermediate mass twinbinaries, are consistent with observed characteristics of Fig. 17.—
The parameter space of twin binaries. Here andthroughout the paper, the core mass m c and binary separation r are given in units of the total mass and radius, respectively, ofan isolated binary component. The top curve, separating detachedand contact binaries, marks configurations of first contact. Themiddle curve, spanning 0 ≤ m c . .
15, is the secular instabilitylimit that separates stable and unstable contact binaries: in thiswork, we find that most secularly unstable systems are also dy-namically unstable to mass transfer across the inner Lagrangianpoint. The bottom curve, separating contact binaries from config-urations with no equilibrium, represents the Roche limit. Binariesthat are unstable or that cannot exist in equilibrium have theircomponents inspiral and merge, a process we follow with dynami-cal calculations. double degenerate central binaries in PNe. In future sim-ulations, a reduced gravitational softening between thetwo cores could allow for an even more precise determi-nation of their final orbital separation. In this way, ob-served orbital parameters of degenerate binaries in PNecould be more readily connected to the properties of theparent stars from which they may have originated.4.5.
Additional Comparisons and Future Work
For core masses m c & .
15, we find that a twin bi-nary can exist stably in deep contact, at separationsall the way down to the Roche limit. In contrast, thesemi-analytic condensed polytrope models of Hjellming& Webbink (1987) predict instead that twin binaries willexperience sustained mass transfer once the componentscome in contact, provided only that m c < .
458 (a rangethat includes the vast majority of giant stars). The pri-mary oversimplification in the semi-analytic treatmentappears to be the approximation that mass outside ofthe Roche lobe cannot help to contain the star within it.Our numerical calculations, however, model the commonenvelope that exists outside of the Roche lobe and thatacts to suppress mass transfer. In addition, our fullythree-dimensional calculations remove the point massand spherical structure approximations implicit to thesemi-analytic method.The models of Hjellming & Webbink (1987) seem bestsuited to semidetached binaries, where there is no com-win Binaries 17mon envelope to complicate the dynamics of the massflow. A comparison of such cases with our results is notpossible, as our work is limited to binaries with identicalcomponents. Natural future work would include relax-ing this constraint so that binaries with mass ratio q = 1can be studied and compared with semi-analytic models.Of particular interest to the binary neutron star problemwould be cases in which the mass ratio deviated from oneby only a few percent or less.The modeling of giants as Γ = 5 / M & M ⊙ , according tocalculations with the TWIN stellar evolution code). Forsuch cases, our treatment of the envelope as a constantentropy, Γ = 5 /3 gas can be legitimately questioned.While the effects of employing more realistic stellar mod- els would be worthwhile to study in future work, we donot expect our results to change qualitatively. Regard-less of the equation of state, gas that flows past the outerLagrangian points will still necessarily carry away a spe-cific angular momentum larger than the system average,forcing the remaining gas to configurations of smaller an-gular momentum per unit mass. We conclude that theinspiral of cores should be a common outcome whenevera real twin binary exceeds the Roche limit.We thank Evghenii Gaburov, Zachary Proulx, AdamSimbeck, Eric Theriault, and the anonymous referee forhelpful input. This material is based upon work sup-ported by the National Science Foundation under grantno. 0703545 and has made use of both the SPLASH visu-alization software (Price 2007) and NASA’s AstrophysicsData System. F.A.R. acknowledges support from NSFgrant PHY–0855592. V.K. acknowledges support fromNSF grant PHY–0969820.