The chaotic four-body problem in Newtonian gravity I: Identical point-particles
Nathan W. C. Leigh, Nicholas C. Stone, Aaron M. Geller, Michael M. Shara, Harsha Muddu, Diana Solano-Oropeza, Yancey Thomas
aa r X i v : . [ a s t r o - ph . S R ] A ug Mon. Not. R. Astron. Soc. , 1–18 (2016) Printed 29 August 2016 (MN L A TEX style file v2.2)
The chaotic four-body problem in Newtonian gravity I:Identical point-particles
Nathan W. C. Leigh , Nicholas C. Stone , , Aaron M. Geller , , Michael M. Shara ,Harsha Muddu , Diana Solano-Oropeza , Yancey Thomas ⋆ Department of Astrophysics, American Museum of Natural History, Central Park West and 79th Street, New York, NY 10024 Columbia Astrophysics Laboratory, Columbia University, New York, NY, 10027, USA Einstein Fellow Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy,Northwestern University, 2145 Sheridan Rd, Evanston, IL 60208, USA Adler Planetarium, Dept. of Astronomy, 1300 S. Lake Shore Drive, Chicago, IL 60605, USA Student Research and Mentoring Program, Department of Education, American Museum of Natural History,Central Park West and 79th Street, New York, NY 10024
29 August 2016
ABSTRACT
In this paper, we study the chaotic four-body problem in Newtonian gravity. Assum-ing point particles and total encounter energies
0, the problem has three possibleoutcomes. We describe each outcome as a series of discrete transformations in energyspace, using the diagrams first presented in Leigh & Geller (2012; see the Appendix).Furthermore, we develop a formalism for calculating probabilities for these outcomesto occur, expressed using the density of escape configurations per unit energy, andbased on the Monaghan description originally developed for the three-body problem.We compare this analytic formalism to results from a series of binary-binary encoun-ters with identical point particles, simulated using the
FEWBODY code. Each of ourthree encounter outcomes produces a unique velocity distribution for the escapingstar(s). Thus, these distributions can potentially be used to constrain the origins ofdynamically-formed populations, via a direct comparison between the predicted andobserved velocity distributions. Finally, we show that, for encounters that form sta-ble triples, the simulated single star escape velocity distributions are the same as forthe three-body problem. This is also the case for the other two encounter outcomes,but only at low virial ratios. This suggests that single and binary stars processed viasingle-binary and binary-binary encounters in dense star clusters should have a uniquevelocity distribution relative to the underlying Maxwellian distribution (provided therelaxation time is sufficiently long), which can be calculated analytically.
Key words: gravitation – binaries (including multiple): close – globular clusters:general – stars: kinematics and dynamics – scattering – methods: analytical.
The N -body problem is a longstanding issue hailing fromSir Isaac Newton’s day (Newton 1686), as discussed ina previous paper in this series (Leigh & Geller 2015;Valtonen & Karttunen 2006). It has been the subject ofintense study for centuries, with a flurry of rapid progressover the last few decades due to the introduction of com-puters. And yet, small- N chaos is an essential puzzle relatedto dynamical phenomena such as direct stellar and binaryinteractions in globular (e.g., Heggie 1975; Leigh & Geller ⋆ E-mail: [email protected] (NWCL) N -body, problem has never been fully solved analyti-cally; instead, approximate methods are used along withsimplifying assumptions that make the problem tractable.Examples include simple three-body configurations suchas Burrau’s Pythagorean problem (Burrau 1913), therestricted three-body problem often applied to the Earth-Sun-Moon system, periodic solutions to the non-hierarchicalthree-body problem (e.g. Suvakov & Dmitrasinovic 2013)and finally complex numerical simulations of very large- N Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y. systems. These N -body simulations can be prohibitivelycomputationally expensive, with integration times scalingwith the number of particles as N . Furthermore, starcluster simulations that include both stable and chaotically-interacting multiple star systems can significantly slow theintegration times down, since the need for small or shorttime-steps during these small- N interactions occurringwithin the larger framework of the cluster simulationcan contribute to a significant increase in the overallcomputer run-times of the simulations (Hurley et al. 2005;Valtonen & Karttunen 2006; Geller, Hurley & Mathieu2013). This is particularly problematic, since simula-tions have shown that encounters involving binaries andtriples can be crucial for not only the overall clusterevolution (e.g. Hut 1983a; Hut, McMillan & Romani1992), but also for the formation of exotic populationssuch as blue stragglers (Perets 2009; Leigh & Sills 2011;Geller, Hurley & Mathieu 2013; Naoz & Fabrycky 2014)and even accreting (Mapelli & Zampieri 2014) or massive(Stone, Kuepper & Ostriker 2016) black holes. What’smore, observations have now revealed that higher-ordermultiple star systems are present in young star clusters innon-negligible numbers (e.g. Leigh & Geller 2013).We argue here that N = 4 is an optimal number forstudying the chaos of gravitationally-interacting particles,with an emphasis on going beyond the already well-studiedthree-body problem. This is because N = 4 offers a rea-sonable balance between the computer run-times for thesimulations, and the statistical significance of our analysis,which depends directly on the total number of simulationsperformed for a given set of initial conditions. Theoreti-cally, binary-binary encounters should dominate over single-binary encounters in any star cluster with a binary fraction f b &
10% (Sigurdsson & Phinney 1993; Leigh & Sills 2011).For these reasons, we focus on the chaotic N = 4 prob-lem in this paper. We further narrow our focus to consideronly identical, equal-mass point particles. Regardless, in Sec-tion 4, we use our results to predict the expected behavior ofinteractions involving non-identical particles with differentmasses, which we intend to test directly in future work.The three-body system with N = 3 typically evolves viaa series of close triple encounters (e.g. Agekyan & Anosova1967; Anosova & Orlov 1994). Between each such event, oneof the objects is temporarily ejected but remains boundto the three-body system. This object recoils some dis-tance from the remaining binary before returning to initiateanother triple encounter. Eventually, one of the bodies isejected with a sufficiently high velocity to become unbound,and it escapes to infinity. This chaotic progression or evolu-tion can be simplified as follows.The time evolution of the chaotic N = 3 problem canbe broken down into a single discrete transformation, whichoccurs in energy and angular momentum space. The initialconditions for the encounter define the relative energies andangular momenta of the single star and the binary. Aftera transformation is applied in energy and angular momen- We will use the term ”ejection” to refer to these types of events,where the ejected object actually remains bound to the system.We will use the term ”escape” to refer to events where particlesbecome unbound from the remaining system. tum space, the final state of the system is qualitatively thesame as the initial state (i.e., a single star and a binary areleft over), but quantitatively the relative energies and angu-lar momenta are distributed differently among the particlesin the final state. Hence, the time evolution of the systemconnecting the initial and final states can, to first-order, bedescribed as a single transformation in energy and angularmomentum space, where the intermediary chaos is neglectedsuch that statistical ensembles of outcomes are more easilyconsidered, rather than individual choices of scattering pa-rameters.In this paper, we extend this approach to the four-bodyproblem which, as we will show, can similarly be brokendown into a series of discrete transformations in energy space(and angular momentum space). In the case of the four-bodyproblem with point particles, the transformations lead to oneof three final states for the system (instead of one final statein the three-body problem). We use the numerical scatteringcode
FEWBODY to simulate a series of binary-binary encoun-ters involving identical point particles, for different valuesof the virial ratio. In Section 2, we describe the simulationsused in this study, and present the resulting distributionsof final virial ratios, escape velocities and encounter times.We go on to apply the time-averaged virial approximationto better understand the response of the interacting systemto particle escapes in Section 3, and adapt the Monaghanformalism, originally derived by Monaghan (1976a) for thethree-body problem, to describe the distributions of singlestar escape velocities. The significance of our results for as-trophysical observations are discussed in Section 4, and ourkey results are summarized in Section 5.
In this section, we present the numerical scattering experi-ments used to study the time evolution of the chaotic 4-bodyproblem as a function of the initial virial ratio.
We calculate the outcomes of a series of binary-binary(2+2) encounters using the
FEWBODY numerical scatteringcode . The code integrates the usual N -body equations inconfiguration- (i.e., position-) space in order to advance thesystem forward in time, using the eighth-order Runge-KuttaPrince-Dormand integration method with ninth-order errorestimate and adaptive time-step. For more details about the FEWBODY code, we refer the reader to Fregeau et al. (2004).The outcomes of these 2+2 encounters are studied as afunction of the initial virial ratio k , defined as: k = T + T E b , + E b , , (1)where the indexes 1 and 2 correspond to the two initial bina-ries. The initial kinetic energy corresponding to the centreof mass motion of binary i is: T i = 12 m i v , i , (2) For the source code, see http://fewbody.sourceforge.net. he 4-Body Problem in Newtonian gravity I where m i = m i , a + m i , b is the total binary mass and v inf , i is the initial centre of mass velocity for binary i . The initialorbital energy of binary i is: E b , i = − Gm i , a m i , b a i , (3)where m i , a and m i , b are the masses of the binary componentsand a i is the initial orbital separation. Given this definitionfor the virial ratio, k = 0 corresponds to the binaries startingfrom rest and k = 1 corresponds to a relative velocity equalto the critical velocity v crit , defined as the relative velocityat infinity needed for a total encounter energy of zero. Thatis, Equation 1 can be re-written as: k = (cid:16) v rel v crit (cid:17) , (4)where v rel is the initial relative velocity at infinity betweenthe binaries and v crit the critical velocity. We consider initialvirial ratios of k = 0.00, 0.04, 0.16, 0.36, 0.64 and 1.00.All objects are point particles with masses of 1 M ⊙ . Allbinaries have a i = 1 AU initially, and eccentricities e i = 0.We fix the impact parameter at b = 0 for all simulations.The angles defining the initial relative configurations of thebinary orbital planes and phases are chosen at random. Weperform 10 numerical scattering experiments for every ini-tial virial ratio. For comparison purposes, we also run oneset of 10 simulations at k = 0 .
04 with identical initial con-ditions, but assuming initial binary orbital separations of 1AU and 10 AU. These simulations are shown by the solidtriangles in Figure 4 and are compared to our fiducial set ofsimulations for two initially 1 AU binaries in Figure 9.We use the same criteria as Fregeau et al. (2004) to de-cide when a given encounter is complete. To first order, thisis defined as the point at which the separately bound hier-archies that make up the system are no longer interactingwith each other or evolving internally. More specifically, theintegration is terminated when the top-level hierarchies havepositive relative velocity and the corresponding top-level N -body system has positive total energy. Each hierarchy mustalso be dynamically stable and experience a tidal perturba-tion from other nodes within the same hierarchy that is lessthan the critical value adopted by FEWBODY , called the tidaltolerance parameter. For this study, we adopt the a tidaltolerance parameter δ = 10 − for all simulations. Thischoice for δ , while computationally expensive, is needed tomaximize the accuracy of our simulations, and ensure thatwe have converged on the correct encounter outcome (seeGeller & Leigh 2015 for more details). The chaotic four-body problem involving point particles hasthree possible outcomes, provided the total encounter energysatisfies E That is, when an encounter is over, the The more stringent the tidal tolerance parameter is chosen tobe, the closer to a ”pure” N -body code the simulation becomes. We ignore encounters producing four single stars, since theserequire positive total encounter energies and hence very large rel-ative velocities at infinity. As explained in Appendix A, such en-counters are unlikely to occur in dense stellar systems, since dy-namically “soft binaries dissociate rapidly from the cumulative remaining configuration is described by one of the followingoutcomes:* two binaries (2+2)* a triple and a single star (3+1)* a binary and two single stars (2+1+1)Although complete ionizations, or 1+1+1+1 outcomes, aretechnically possible in some scenarios, they are extremelyrare in any realistic star cluster environment. We defer amore in depth discussion of this to Appendix A.The (qualitative) time evolution of three example sim-ulations, each ending in one of the three outcomes listedabove, are depicted in the top insets of Figures 1, 2 and 3 us-ing the schematic diagrams first introduced in Hut (1983b).The three outcomes are also depicted schematically in thebottom insets of Figures 1, 2 and 3 using the energy dia-grams first presented in Leigh & Geller (2012). The latterdiagrams quantify the exchange of energy between particlesduring an encounter. Briefly, each angle of the polygon cor-responds to the fraction of the total encounter energy con-tained in star i . The solid and dashed lines connect particlesthat are gravitationally bound and unbound, respectively.For more details regarding the details of these energy dia-grams, see the Appendix in Leigh & Geller (2012).The energy diagrams shown in the bottom insets of Fig-ures 1, 2 and 3 correspond to nearly instantaneous, or ”dis-crete” events (i.e., close or strong encounters between indi-vidual particles); to first order, each time an object passesnear the centre of mass of the system, its motion is signifi-cantly perturbed and it recedes from the centre of mass ona new orbit, before returning to repeat the process. Thus,in energy space, the time evolution of the four-body systemcan be described by a series of discrete transformations inenergy space, as depicted schematically in the bottom panelsof Figures 1- 3. In principle, this has the potential to sim-plify the problem, since we are not concerned with the evo-lution in position and velocity space, where long excursionscan occur with little to no exchange of energy (or angularmomentum) between particles. This same procedure can beapplied analogously in angular momentum space.Importantly, it is the first escape event that decides theoutcome of a chaotic 4-body encounter with non-positivetotal energy; that is, the first object to become gravitation-ally unbound and leave the system with positive total en-ergy (ignoring the binding energies of any left-over binaryorbits) decides the final outcome. This can be understoodas follows. If a binary is the first object to escape, thenthe encounter is over and the outcome is 2+2. If, on theother hand, a single star is the first object to escape, thenthe encounter either terminates immediately if the remain-ing triple is dynamically stable (i.e., the outcome is 3+1),or continues to evolve chaotically until the temporary tripleeventually disrupts (i.e., the outcome is 2+1+1). It is notpossible for an initially unstable isolated triple to becomestable (Littlewood 1952). Thus, the stability of the tripleconfiguration immediately after the escape of the first sin- effects of many weak three-body encounters. Only the far tail ofthe Maxwellian will be capable of completely ionizing two hardbinaries, and the rate of this will be exponentially suppressedunless both binaries are only marginally hard. Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Figure 1.
The time evolution of an example 2+2 simulation pro-ducing a 2+2 outcome is shown schematically. In the top diagram,time increases from left to right along the x-axis, whereas the y-axis shows the progression of configurations. In the bottom dia-gram, with two distinct insets, shows the initial and final states ofthe encounter in energy space, using the energy diagrams first in-troduced in Leigh & Geller (2012). Note that most 2+2 outcomesin our simulations produce binaries with a very large differencein their orbital energies, even more exaggerated than shown here. gle star, and, more generally, the first escape event, entirelydecides the encounter outcome.In Figure 4, we show the fraction of simulations thatend in each of our three encounter outcomes, as a function ofthe initial virial ratio. Note that these results are sensitive toour choice of particle masses and initial binary energies. Theblue, red and green points correspond to the 2+1+1, 2+2and 3+1 outcomes, respectively. Encounters producing twosingle stars dominate for our choice of the initial conditions,giving outcome fractions of ∼ k ∼ ∼
10% of the outcomes are 2+2, and ∼
10% are 3+1). Asthe virial ratio increases, however, the fraction of 3+1 out-comes slowly decreases reaching ∼
0% at k = 1, while thefractions of 2+2 and 2+1+1 outcomes both increase by ∼ afew percent. These results are in good agreement with thosepresented by Mikkola (1983), who performed numerical scat-tering experiments using very similar initial conditions (i.e.,all identical particles, and identical binaries).In Figure 5 we show the distributions of initial and fi-nal virial ratios, and in Figure 6 we show the distributionsof initial and final velocities (in km s − ). Figure 7 shows thedistributions of the fractional change in energy, for all es-cape events, provided they satisfy ∆ E / | E | >
0. That is, thefractional change in energy is plotted for the escaping sin-gle star during 3+1 outcomes (green), and for both escapingsingle stars during 2+1+1 outcomes (blue). Note that thequantity ∆ E / | E | is positive if the escaping object is a single Figure 2.
The time evolution of an example 2+2 simulation pro-ducing a 3+1 outcome is shown schematically. In the top diagram,time increases from left to right along the x-axis, whereas they-axis shows the progression of configurations. The bottom dia-gram, with two distinct insets, shows the initial and final states ofthe encounter in energy space, using the energy diagrams first in-troduced in Leigh & Geller (2012). Note that most 3+1 outcomesin our simulations produce a relatively low-velocity escaping sin-gle star and a triple with a very large difference in their inner andouter orbital energies, even more exaggerated than shown here. star, or if the escaping object is a binary and the absolutevalue of its orbital energy is less than its (translational) ki-netic energy. Hence, to avoid confusion, we omit the 2+2outcome in Figure 7, and note that in Figures 2 and 3 thequantity ∆ E / | E | corresponds to the angle(s) of the polygonin the final state (i.e., the lower right panels) with dashedsides. We caution that at high virial ratios, very few triplesform, and small-number statistics become a concern in thegreen histograms. Note that the distributions of single starescape velocities are noticeably different for encounters thatproduce stable triples, relative to those that produce twosingle stars and a binary.In Figure 8, we show the distributions of total encounterdurations (in years) for each of our three encounter out-comes. Note that the initial displacement toward longer en-counter durations, which is more severe for low virial ratios,is due to the initial infall time for the two binaries, whichis determined by our choice for the tidal tolerance parame-ter (and is hence not physical). The encounter durations areshortest for the the 2+1+1 outcome, and longest for the 3+1outcome. This is at least due in part to the fact that it is theouter orbit of the top-level node in FEWBODY that decides thefinal stability of the system, and the encounter cannot becompleted until one orbital period has occurred. Notwith-standing, these results serve to further strengthen the con-clusion that each outcome can be regarded, to first order, asa unique transformation in energy and angular momentum he 4-Body Problem in Newtonian gravity I Figure 3.
The time evolution of an example 2+2 simulation pro-ducing a 2+1+1 outcome is shown schematically. In the top dia-gram, time increases from left to right along the x-axis, whereasthe y-axis shows the progression of configurations. The bottomdiagram, with four distinct insets, we show the time evolution ofthe encounter in energy space, using the energy diagrams first in-troduced in Leigh & Geller (2012). This figure is reproduced fromFigures A2 and A3 in Leigh & Geller (2012). space, since each outcome produces distinct distributions ofvirial ratios, escape velocities and encounter durations.Finally, in Figure 9, we compare our fiducial simulationsinvolving two identical binaries with initial separations of a = a = 1 AU (left panels) to what we obtain assuminginstead that the binaries have initial separations of a =1 AU and a =10 AU (right panels). This is done for onechoice of the initial virial ratio, namely k = 0 .
04, since thecomputational expense for these additional simulations iseven higher than for our fiducial runs. In the top, middleand bottom insets we show, respectively, the distributionsof final virial ratios, escape velocities (in km s − ) and totalencounter durations (in years). Where comparisons are pos-sible, our results are in overall good agreement with those ofMikkola (1984), who performed a similar study of encoun-ters involving binaries with unequal energies.A few trends are immediately clear from the compar-isons shown in Figure 9. First, the fraction of simulationsproducing stable triples is higher by a factor &
2. Corre-spondingly, the distributions of escaper velocities have beenshifted to lower velocities for the 1 AU + 10 AU case which,as previously illustrated, are more conducive to triple for-mation. In the bottom panels, we see that the encounter du-rations corresponding to triple formation are much shorterfor the 1 AU + 10 AU case, relative to our fiducial simu-lations, whereas encounters producing two binaries tend tohave longer encounter durations. We speculate that this isdue to the higher angular momentum characteristic of theseadditional simulations with 1 AU + 10 AU binaries; in or-der to conserve angular momentum any ejected binary will
Figure 4.
The fraction of simulations that end in each of ourthree encounter outcomes are shown as a function of the initialvirial ratio. The blue, red and green points correspond to the2+1+1, 2+2 and 3+1 outcomes, respectively. For comparison,we also show by the solid triangles these fractions for k = 0 . a =1 AU and a = 10 AU. be most likely to have a large orbital separation and hencea small absolute orbital energy. This puts it very close tothe dissociation border, such that even a small additionalamount of positive energy will dissociate the binary, whichlikely results from tidal effects imparted by the remaining(more compact) binary. More work needs to be done to bet-ter understand the underlying physics responsible for thedifferences highlighted in Figure 9, which arise due to the dif-ferent ratios between the initial binary orbital separations.The take-away message from this figure is that larger ratiosbetween the initial binary orbital separations correspond toa higher probability of triple formation at a given virial ra-tio, relative to what is shown for our fiducial simulations. Consider a small star cluster of N identical point-particleseach with mass m . If the system is in dynamical equilibrium,then the mean radius R of the cluster is determined by thevirial radius (Valtonen & Karttunen 2006): R ∼ GMv = GM | E | , (5)where M = Nm is the total cluster mass, v rms = ( | E | /M ) / is the the root-mean-square velocity of the virialized systemand E is its total energy. Let us assume that a single starleaves the cluster, escaping to spatial infinity with a positivetotal energy ∆ E . After the system re-achieves dynamical Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Figure 5.
The distributions of final virial ratios are shown. The blue, red and green histograms correspond to the 2+1+1, 2+2 and3+1 outcomes, respectively. Each panel shows the distributions for a different value of the initial virial ratio. All histograms have beennormalized by the total number of simulations that resulted in the corresponding outcome. equilibrium, the new virial radius R ′ is: R ′ ∼ G ( M − ∆ M ) | E | + ∆ E , (6)where ∆
M > E is positive if the escaping object is asingle star, or if the escaping object is a binary and the abso-lute value of its orbital energy is less than its (translational) kinetic energy. Equation 6 can be re-written: R ′ = GM E (1 − ∆ M/M ) (1 + ∆ E/ | E | )= R (1 − ∆ M/M ) (1 + ∆ E/ | E | )= αR. (7)If α > < α > E/ | E | <
0, and hence can only occur when a binary is the he 4-Body Problem in Newtonian gravity I Figure 6.
The distributions of initial and final escape velocities are shown in km s − . The blue, red and green histograms correspond tothe 2+1+1, 2+2 and 3+1 outcomes, respectively. We use different line types to show the distributions for the different objects producedfor each outcome; solid lines correspond to multiple star systems (i.e., binaries and triples) and dashed lines correspond to single stars(with the exception of the 2+2 case, for which the dashed red lines simply correspond to the other binary). Each panel shows thedistributions for a different value of the initial virial ratio. All histograms have been normalized by the total number of simulations thatresulted in the corresponding outcome. escaping object (i.e., for a 2+2 outcome). In encounters fea-turing very close pericenter passages, orbital expansion canalso in principle occur through dissipative processes such asgravitational wave emission (Peters 1964) or tidal excitationof stellar modes (Press & Teukolsky 1977), but we neglectfinite-size and general relativistic effects in this work. In gen-eral, from Equation 7, we see that for ejections correspond- ing to very high escape velocities but low particle massescompared to M = Nm (i.e., | ∆ E/E | ≫ | ∆ M/M | ) the sys-tem should expand. Conversely, if | ∆ E/E | ∼ | ∆ M/M | , thesystem will typically contract.Returning to the four-body problem, if ∆ E >
0, whichis always the case for the escape of a single star (since wehave taken the absolute value of the system binding energy
Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Figure 7.
The distributions of the fractional change in energy are shown for all escape events. The different panels show the resultsfor different values of the initial virial ratio. The colour coding is the same as in Figures 4-6. We plot the fractional change in energyonly if ∆ E / | E | > | E | in Equation 7), then α < E / | E | value. However, the converse is notalways true, due to angular momentum conservation and thechaotic nature of the system evolution. In other words, if asingle star escapes with a low ∆ E / | E | value, it is not guar-anteed that a stable triple will form. As shown in Figure 7,at low virial ratios, stable triples are only more likely to be he 4-Body Problem in Newtonian gravity I Figure 8.
The distributions of total encounter durations are shown in years. The blue, red and green histograms correspond to the2+1+1, 2+2 and 3+1 outcomes, respectively. Each panel shows the distributions for a different value of the initial virial ratio. Allhistograms have been normalized by the total number of simulations that resulted in the corresponding outcome. associated with low ∆ E / | E | values by a factor of ∼ a few,whereas this factor increases markedly with increasing virialratio. Finally, as shown in Figure 7, due to conservation ofenergy and linear momentum, the vast majority of encoun-ters leading to stable triple formation correspond to escapeevents satisfying ∆ E / | E | ≪ E / | E | (as shown in Figure 7) since these tend toresult in the greatest contraction of the system post-escape(but not always; see Figure 7). More compact three-bodyconfigurations have a lower probability of being stable, dueto the smaller parameter space corresponding to stabilitythat is available to them. Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Figure 9.
The distributions of final virial ratios (top panels), escape velocities (in km s − ; middle panels) and total encounter durations(in years; bottom panels) for an initial virial ratio of k = 0 .
04. The left panels show the results for simulations involving two binarieswith identical initial separations of a = a =1 AU, whereas the right panels show the results for simulations involving binaries withinitial separations of a = 1 AU and a =10 AU. The blue, red and green histograms correspond to the 2+1+1, 2+2 and 3+1 outcomes,respectively. All histograms have been normalized by the total number of simulations that resulted in the corresponding outcome. Next, we adapt the Monaghan formalism to the four-body problem. This was originally developed by Monaghan(1976a) for the three-body problem to estimate the statis-tical distribution of outcomes from three-body scatterings.We make a historical note here for completeness, namelythat the Monaghan distribution derived a power-law in- dex of -5/2 for the distribution of binary binding ener-gies (originally derived by Jeans (1928)). This was shownto provide a poor fit to numerical data, but was latercorrected by Valtonen & Karttunen (2006) using the loss-cone approximation, changing the power-law index to -9/2 (originally derived by Heggie (1975) and confirmed inSaslaw, Valtonen & Aarseth (1974)).The Monaghan formalism relies on the density of es- he 4-Body Problem in Newtonian gravity I cape configurations per unit energy σ , obtained by integrat-ing over the phase space volume. As we will show, this canbe done relatively precisely for the 3+1 outcome (i.e. theescaping object is a single star and the left-over system isa stable hierarchical triple), whereas this is not the case forthe other two encounter outcomes. This is due to a strongparallel that can be drawn between the 3+1 outcome andthe escape of a single star during a three-body interaction, asthe outer orbit of the final stable triple contains a negligibletotal energy. We will return to this point in the subsequentsub-section.Given that there are three possible outcomes for the E < P + P + P , (8)where P , P and P are the probabilities corre-sponding to the encounter outcomes, respectively, 2+2, 3+1and 2+1+1. Equation 8 can also be re-written as:1 = Z (cid:16) σ ( | E | ) + σ ( | E | ) + σ ( | E | ) (cid:17) d | E | , (9)where each σ value denotes the (normalized) density of es-cape configurations per unit energy for each of the indicatedoutcomes. We will return to Equation 9 in Section 3.2. Consider a chaotic four-body interaction that produces astable hierarchical triple system. The total energy of thefour-body system is: E = E e + E b , (10)where the total encounter energy is divided between the es-caper energy E e , and the energy of the left-over (bound)system E b which is in this case a stable triple. Hence, thefinal system energy can be decomposed into the inner andouter orbital energies: E b ≈ E b , a + E b , b ∼ E b , b , (11)where the total system energy can be decomposed into theinner and outer orbital energies: E b , a = − Gm a , m a , a b , a (12)and E b , b = − Gm b , m b , a b , b , (13)where m b = m b , + m b , . We assume that a b , a = βa b , b forsome positive constant β &
10 (i.e. for a stable hierarchicaltriple; Mardling 2001). This gives, to first order, E b ∼ E b , b .For now, we stick with the variable E b for the rest of thederivation of the distribution of single star escape velocities,to highlight the analogy with the three-body problem.The density of escape configurations per unit energy is: σ = Z ... Z δ (cid:16) p m + V e + E b − E (cid:17) d r e d p e d r d p , (14) Note that we have replaced the variables E s and E B in theoriginal Monaghan derivation with, respectively, E e and E b . where m = m b m e M (15)and M = m b + m e (16)Now, assuming (49/4)( a b , a / r e ) for the loss-cone factor asin Valtonen & Karttunen (2006), the integrations over p e and r e are carried out analogously and we obtain: σ = 98 √ π ( GMR ) / m ( Gm b , m b , ) Z ... Z d r d p | E b | , (17)where the upper limit R of the integration over the r e param-eter is a free parameter, but must be chosen to be relativelysmall (i.e., taken to be R = 3 a in Valtonen & Karttunen2006, where a is the initial semi-major axis of the binariesgoing into the four-body interaction, which we assume forsimplicity are equal).Here, we recall the approximation E b ∼ E b , b , whichimplies that the contribution to the total energy from theouter orbit of the stable hierarchical triple can be ignoredin the derivation for σ . In other words, all of the totalencounter energy is distributed among the escaping star andthe inner binary of the triple in the final (stable) configu-ration. Only the angular momentum of the outer orbit issignificant. Hence, as we show in the next section, we canmodel a 3+1 outcome as a single escape event, with m b = m b , b = m b , + m b , and m e = M − m b , or m e = m b = M /2 for all identical particles. In other words, we model a3+1 outcome using the same formulation as originally de-veloped by Monaghan (1976a) for the three-body problem,but modified slightly by assuming that two stars are effec-tively ”ejected” together. One of these stars escapes to in-finity and the other ends up loosely bound to the remainingbinary. Only the escaping single star and the inner binaryof the triple exchange significant energy, while the remain-ing star becomes the outer triple companion to the binaryand carries only a negligible amount of the total encounterenergy.The rest of the derivation for the 3+1 outcome is effec-tively identical to the original Monaghan formalism, givenour assumption that only the inner binary of the triple andthe escaping single star exchange significant energy. Hence,we do not repeat the full derivation here, and only describethe essential steps (see Valtonen & Karttunen 2006 for moredetails). Adopting spherical polar coordinates (r, θ , Φ), the(inner) binary energy can then be written: E b = 12 p µ − Gm b , m b , r , (18)where µ = m b , m b , /( m b , + m b , ). Following We note that this loss-cone factor was determined from 1+2scattering experiments, and should be tested specifically for thefour-body problem in future studies concerned with the normal-ization factor in front of the integral in Equation 20. This nor-malization factor enters into the ”branching ratio”, or the proba-bility of the encounter ending in a given outcome. Crucially, ourassumption for this loss-cone factor does not affect our derived ve-locity distributions, only the normalization factor or ”branchingratio”. Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Valtonen & Karttunen (2006), several integrations canbe carried out, leaving only the integrals over the binaryenergy E b and angular momentum L, or: σ = 2 × π ( Gm b , m b , ) / R / m / M − / m Z Z d | E b || E b | / LdL. (19)Re-writing Equation 19 in terms of the binary eccentricitye gives, finally: σ = 98 π ( Gm b , m b , ) / R / m / M − / m µ Z Z dE b | E b | / ede. (20)The quantities following the integral signs in Equa-tion 20 correspond to the distributions over which one mustintegrate in order to obtain the total phase space volumecorresponding to the 3+1 outcome. Thus, the distributionof binary energies | E b | , normalized to unity, is: f ( | E b | ) d | E b | = 3 . | E | / | E b | − / d | E b | . (21)This corresponds to the low total angular momentumcase. More generally, numerical scattering experimentshave shown that, in order to cover the full range intotal angular momenta, Equation 21 can be re-written(Valtonen & Karttunen 2006): f ( | E b | ) d | E b | = ( n − | E | n − | E b | − n d | E b | , (22)where, for the three-body problem, the power-law index n ranges from n = 3 at L = 0 to n = 14.5 at L =0.8 L max , and L is the total encounter angular momentum(Valtonen & Karttunen 2006). Hence, n is a positive con-stant that depends on the total angular momentum via thesubstitution L = L / L max , or: n − L , (23)as found from numerical scattering experiments for thethree-body problem (Valtonen et al. 2003). Importantly, wehave not verified in this paper that Equation 23 holds forthe chaotic four-body problem, and this should be tested infuture work. To do this, significantly more simulations mustbe performed, in order to sample the full range of angularmomentum space. In this paper, we focus on the minimumangular momentum case (i.e., zero impact parameter and n = 3 in Equation 23), for which Equation 23 does indeedprovide a good match to the simulations (see Figures 10and 11).Equation 22 can ultimately be used to derive the distri-bution of escape velocities for the escaping single star pro-duced during a 3+1 outcome, as given by Equation 24 with m e = m b = M /2 for all identical particles. That is, theintegration in Equation 20 can instead be re-written to becarried out over the single star escape velocity v e to ob-tain the escape velocity distribution f ( v e ) dv e . This givesthe following functional form (Equations 7.19 and 7.26 inValtonen & Karttunen 2006): f ( v e ) dv e = ( n | E | n − ( m e M/m b )) v e dv e ( | E | + ( m e M/m b ) v ) n . (24)As we will show, the above derivation works only toreproduce the distribution of escape velocities for 3+1 out-comes. For the 2+2 and 2+1+1 outcomes, the escape ve-locity distributions are reproduced by Equation 24 at lowvirial ratios only. This is because the assumption E b ∼ E b , b can only be applied for the 3+1 outcome. In the 2+2 and2+1+1 cases, this approximation is not valid, since signif-icant additional energy can be stored in either of the tworemaining orbits (after the first ejection event, whether itbe a single star or a binary). This seems to account for theadditional positive energy imparted at large virial ratios tothe escaping object in the 2+2 and 2+1+1 outcomes, rel-ative to the Monaghan formalism (see Figures 9 and 10).We will return to this important and interesting result inSection 4. For now, we note that we intend to apply thistheory to the other two encounter outcomes in a forthcom-ing paper. This will ultimately require integrating over theadditional degrees of freedom brought in by the additionalbinary orbit (i.e., relative to the more familiar three-bodyproblem). In this section, we compare the theoretical distributions ofescaper velocities derived in the previous section for the 3+1outcome to the results of our numerical scattering exper-iments. For comparison, we also show the distributions ofescaper velocities for the 2+1+1 (see Figure 10) and 2+2(see Figure 11) outcomes, beginning with the former.In Figure 10 we show, for different initial virial ratios,the distribution of escape velocities for the single star for a3+1 outcome (green lines), as well as for both single starsfor a 2+1+1 outcome (solid and dashed blue lines). The keypoint to take away from Figure 10 is that escaping singlestars that leave behind a stable triple have a unique velocitydistribution relative to the corresponding distributions forthe escaping single stars leaving behind a binary. Note thatthe increased scatter in the green histograms seen at largevirial ratios is due to the fact that the number of simulationsresulting in a 3+1 outcome is very small here (i.e., σ → E → v e > df/dv e = 0. The result can be written inthe general form: v e , peak = ǫ r ( M − m e ) m e M p | E | , (25)where ǫ is a positive constant of order unity equal to: ǫ = 1 q n − . (26)Equation 24 provides a reasonable fit to the simulated he 4-Body Problem in Newtonian gravity I escape velocity distributions for the 2+1+1 outcome onlyat low virial ratios, as shown by the solid black lines inFigure 10. The agreement gets worse as we go to largervirial ratios, and begins to better approximate the simu-lated escape velocity distribution for the 3+1 outcome. Atthe same time, the distribution of single star escape veloci-ties for the 2+1+1 outcome is remarkably insensitive to theinitial virial ratio (note that the distributions for both es-caping single stars overlap and are almost identical for allvirial ratios), as originally found by Mikkola (1983). This isa curious result worthy of further investigation. As pointedout by Mikkola (1983), this trend can perhaps be understoodby noting that the mean change in the total binding energydecreases roughly linearly with increasing impact energy, orvirial ratio.We see a similar scenario in the distributions of binaryescape velocities for the 2+2 outcome in Figure 11, as shownby the solid and dashed red lines. That is, Equation 24 pro-vides a good agreement to the simulated data at low virialratios only, as shown by the solid and dotted black lines.At larger virial ratios, the agreement becomes significantlyworse. Unlike the blue histograms in Figure 10, however,the red histograms in Figure 11 do change with increasingvirial ratio, shifting to slightly higher binary escape veloci-ties. Hence, at large virial ratios, Equation 24 under-predictsthe peak escaper velocity.Only the distribution of single star escape velocities forthe 3+1 outcome, shown by the green histograms, matchesEquation 24 for all virial ratios. This is illustrated via thesolid and dashed black lines in, respectively, Figures 10and 11. We will return to this intriguing result in Section 4. The key point to take away from the preceding sections isthat the different outcomes of the chaotic four-body problemgenerally correspond to different discrete transitions in en-ergy and angular momentum space. Consequently, the dis-tributions of escape velocities are unique for the differentencounter outcomes. Thus, these distributions can poten-tially be used to constrain the origins of dynamically-formedpopulations, whether they are single, binary or triple stars.This is done via a direct comparison between the predictedand observed velocity distributions. For example, compar-isons can be made to populations of single stars thoughtto be connected to triple formation. These types of four-body encounters have also been proposed to produce hy-pervelocity stars via the capture of binaries around mas-sive black holes when they interact with interloping triplestars (Perets 2009). Conversely, comparisons can be madeto the observed velocity distributions of triples thought tobe formed dynamically, preferably in massive clusters withvery long relaxation times.Interestingly, the distributions of escape velocities forthe chaotic four-body problem are surprisingly similar to thedistributions originally derived by Monaghan (1976a) for thethree-body problem.
The single star escape velocity distri-butions during stable triple formation (i.e., 3+1 outcomes)are well described by the Monaghan formalism at all virialratios and arbitrary total angular momentum, provided weassume that both the escaping single star and the outer third companion of the resulting triple are effectively ”ejected” to-gether , analogous to a 2+2 outcome. The escape velocitydistributions for the 2+1+1 and 2+2 outcomes are also welldescribed by the Monaghan formalism, but only at low virialratios (i.e., k . &
80% of all encounters between identical binaries result inthe 2+1+1 outcome, as shown in Figure 4. As shown in Fig-ure 8, the resulting distribution of single star escape veloci-ties is well described by the same escape velocity distributionas originally derived by Monaghan (1976a) for the three-body problem, provided the virial ratio is small. Thus, anypopulation of single and binary stars in a cluster that havegone through either a single-binary (i.e., 1+2) or binary-binary (i.e., 2+2) encounter should have post-encounter ve-locities that are roughly described by Equation 24, assumingthat at least one relatively hard binary is included (to en-sure a small virial ratio) in the encounter and no other fac-tors influenced the final escape velocity distribution (e.g.,relaxation in a star cluster). In a dense globular cluster, forexample, the post-encounter single star velocity distributionshould be described by averaging | E | in Equation 24 overa Maxwellian velocity distribution. This offers a potentialmeans of identifying evidence for dynamical processing indense, massive star clusters with very long relaxation times.That is, we hypothesize that if the relaxation time is verylong compared to the time-scale for 1+2 and 2+2 encoun-ters to occur, then the velocity distributions for both thebinaries and single stars in the cluster core should deviatefrom Maxwellian distributions, approaching more and morethe distributions given by Equation 24 as time goes on (untilsome steady-state balance is achieved).As stated above, the simulated distribution of singlestar escape velocities for the 2+1+1 outcome agrees wellwith the Monaghan formalism, but only at low virial ra-tios. At larger virial ratios (i.e., near unity), the Monaghanformalism over-predicts the fraction of low-velocity escapersrelative to the simulated data. Interestingly, the resultingescaper velocity distributions in Figure 9 (especially for the2+1+1 outcome) appear roughly insensitive to the total en-counter energy | E | (or virial ratio), and only the relativeprobabilities for the 3+1 and 2+1+1 outcomes change. Thissuggests that the Monaghan formula for the escape velocitydistributions at k . | E | ∼
0) should also fit thesimulated distributions for the 3+1 outcome, nearly inde-pendent of | E | (since there is only a weak dependence onvirial ratio). Similarly, the Monaghan formula for the es-cape velocity distributions at k ∼ | E | . To summarize, ignoring the normalization constants(and hence the ”branching ratios” for the three outcomes ofthe four-body problem), we have in the limit | E | → →
0) that ˜ σ ∼ ˜ σ and ˜ σ ∼ ˜ σ ,where ˜ σ corresponds to the density of states in the orig-inal three-body body formulation of Monaghan, divided bythe normalization constant out in front of the integral. Thisneeds to be confirmed in future studies that consider otherinitial conditions than in this paper.In a forthcoming paper, we will consider the implica-tions of different particle masses for the results presented Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.
Figure 10.
The distributions of single star escape velocities are shown in km s − for the 3+1 (green) and 2+1+1 (blue) outcomes.The black lines show the distribution of escape velocities calculated using Equation 24 for a 3+1 outcome and assuming n = 3 (whichcorresponds to the minimum angular momentum case in Equation 23). For the solid and dashed black lines we assume, respectively, m e = m b /3 = M /4 and m e = m b = M /2. The different insets show the distributions for different virial ratios, as indicated. in this paper. We can nonetheless use Equation 24 alongwith the results presented here to extrapolate to extrememass ratios. For instance, consider a scenario where two ofthe four interacting stars are much less massive than theother two. In this case, the low-mass particles have the high-est probability of either escaping or ending up as the outercompanion to the remaining stable triple (i.e., the two mostmassive stars constitute the inner binary of the triple). Ifwe take the limit m e → E becomes (∆ E ) peak = m e v , peak = ǫ | E | /2. In this paper, however, all particleshave the same mass. This yields a peak escaper kinetic en-ergy of (∆ E ) peak = 3 ǫ | E | /8, which is larger than we foundin the limit m e →
0. As shown in Figure 7, encountersthat result in stable triple formation prefer smaller valuesof ∆
E/E . Thus, we predict that four-body encounters in-volving two very low mass particles should have a higher he 4-Body Problem in Newtonian gravity I Figure 11.
The distributions of binary star escape velocities are shown in km s − for the 2+2 outcomes are shown by the solid anddashed red lines. The solid black lines show the distribution of escape velocities calculated using Equation 24 and assuming n = 4.5(this corresponds to small but non-zero total angular momentum, via Equation 23). For comparison, we also show the single star escapevelocity distributions for a 3+1 outcome and assuming as before n = 3, with both m e = m b /3 = M /4 (dashed black lines) and m e = m b = M /2 (dotted black lines). The different insets show the distributions for different virial ratios, as indicated. probability of triple formation than found in this paper,since we assume all identical particles. We predict a sim-ilar phenomenology for encounters with a wider range ofparticle masses, since the most probable stars to escape al-ways have the lowest masses. Hence, typically, encounterswith a range of particle masses should produce outcomesthat satisfy m e ≪ m b ∼ M . These predictions, derivedfrom our results, are in good qualitative agreement with previous studies (e.g. Saslaw, Valtonen & Aarseth 1974;Mikkola & Valtonen 1990; Valtonen et al. 1994).A similar argument can perhaps be applied to encoun-ters with higher total angular momentum since, as shown inFigure 7.11 of Valtonen & Karttunen (2006) for the three-body problem, a higher angular momentum translates into alower peak escape velocity (which is preferred during tripleformation). This should be confirmed in future studies. Like-wise, we have already seen that the probability of stable Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y. triple formation is enhanced when considering encounters ofbinaries with varying internal energies (see Figure 4).Putting this all together, the probability of formingstable triples during chaotic 4-body encounters should in-crease relative to the results shown in this paper, whichfocused primarily on equal mass stars in identical binarieswith zero angular momentum. Varying any of these ideal-ized assumptions (mass ratio, binary binding energy, totalangular momentum) seems likely to increase the probabil-ity of the 3+1 outcome. This is of great interest for for-mation of exotic stars or transient phenomena that canbe catalyzed through the Kozai mechanism in hierarchicaltriples; examples include blue stragglers (Perets & Fabrycky2009), direct-collision Type Ia supernovae (Thompson 2011;Kushnir et al. 2013), or gravitational waves from compactobject inspirals, particularly eccentric compact mergers(Antonini et al. 2016), which may be impossible to produceat interesting rates through non-Kozai channels. Althoughdetailed rate calculations are beyond the scope of this pa-per, the greater understanding of dynamical triple formationrates developed here is a necessary step towards better quan-tifying the production of interesting Kozai-induced exoticain dense stellar systems.
In this paper, we study the chaotic four-body problem inNewtonian gravity, assuming point particles and total en-counter energies
0. The problem has three possible out-comes, and we describe each outcome as a series of dis-crete transformations in energy space using the diagramsfirst presented in Leigh & Geller (2012; see the Appendix).We further adapt the original Monaghan formalism for thethree-body problem to treat the chaotic four-body problem,based on the density of escape configurations per unit en-ergy. This gives theoretical predictions for the probability ofan encounter producing a given set of outcome parameters.We focus on encounters that produce stable triples for thisanalytic derivation, since the outer orbit of the triple car-ries negligible total energy. Hence, for this outcome alone,the addition of a fourth particle to the original Monaghanformalism does not add an additional degree of freedom.We simulate a series of binary-binary encounters withidentical point particles using the
FEWBODY code. We com-pare the resulting single star velocity distributions with ouranalytic predictions, and show that each of the three en-counter outcomes produces a unique escape velocity distri-bution. Thus, these distributions can potentially be usedto constrain the origins of dynamically-formed populations,via a direct comparison between the predicted and observedvelocity (and binary binding energy) distributions. Finally,we show that, for encounters that form stable triples, thesimulated single star escape velocity distributions are wellreproduced by our theoretical predictions, in the low an-gular momentum regime. Interestingly, this is the case forthe other two encounter outcomes as well, but only at lowvirial ratios; as encounters become more energetic, the ana-lytic three-body formalism of Monaghan breaks down. Thispredicts that single (and hard binary) stars processed viasingle-binary and binary-binary encounters in dense starclusters should have a unique and computable velocity dis- tribution (that can be calculated) relative to the underlyingMaxwellian distribution, provided the relaxation time is suf-ficiently long.
APPENDIX A: THE RARITY OF COMPLETEIONIZATIONS DURING 2+2 ENCOUNTERS INSTAR CLUSTERS
If the total energy of the binary-binary encounter is positive,a fourth outcome becomes possible: a 1+1+1+1 total disso-ciation of the two scattered binaries. However, this outcomeis inherently unlikely in any old, collisional stellar system,because all binaries on the soft side of the hard-soft bound-ary will have long ago dissociated, and the surviving systemsthat engage in binary-binary scattering will be dynamicallyhard. For their encounters to have
E > m = M ⊙ .Next, we assume that all binaries follow Opik’s law, with asemimajor axis ( a ) distribution given by A ( a ) = ( a ln( a max /a min )) − . (A1)Here the minimum binary separation is set to a min = 3 R ⊙ ,and the maximum separation a max is set to the hard-softboundary. If the predominant scattering events in the clusterare binary-single (binary-binary) interactions, then a max = Gm/v ( a max = Gm/v ). The critical velocity requiredto render total dissociation energetically possible is v = r Gma + Gma , (A2)where a and a are the binary semimajor axes, drawn from A ( a ).Next, we assume that the relative velocity of the binary-binary encounter is drawn from a Maxwellian distribution, f m ( v ) = 6 r π v v exp (cid:18) − v v (cid:19) . (A3)The fraction of this distribution with velocities above a ve-locity v is F v ( v ) = erfc r vv rms ! + r π vv rms exp (cid:18) − v v (cid:19) . (A4)Here erfc is the complementary error function.The fraction of binary-binary scatterings with E > χ = Z a max a min Z a max a min F m ( v ( a , a )) A ( a ) A ( a )d a d a . (A5)This integral cannot be evaluated analytically, so we plotnumerical results in Figure A1. As the velocity dispersionof a star cluster increases, surviving binaries become harderto dissociate (decreasing χ ), but larger encounter velocitiesbecome accessible in the Maxwellian (increasing χ ). Thesecond effect dominates; χ increases with increasing v rms .Nonetheless, χ ≪ he 4-Body Problem in Newtonian gravity I × - [ km / s ] Figure A1.
The fraction χ of all binary-binary scatterings thatmake total dissociation (1 + 1 + 1 + 1) energetically possible,plotted against cluster root-mean-square velocity v rms . In general χ ≪
1, rendering the fourth outcome of binary-binary scatteringquite unlikely in older collisional star systems where the hard-softboundary sets the upper limit of the binary semimajor axis dis-tribution. The black line shows χ when the hard-soft boundary isset by binary-single scatterings; the blue line shows χ when theboundary is instead set by binary-binary scatterings. focus on the three E < a max in the binary semima-jor axis distribution, 1 + 1 + 1 + 1 outcomes may be muchmore common. ACKNOWLEDGMENTS
A. M. G. is funded by a National Science Foundation As-tronomy and Astrophysics Postdoctoral Fellowship underAward No. AST-1302765. Financial support was providedto N. C. S. by NASA through Einstein Postdoctoral Fellow-ship Award Number PF5-160145. N. C. S. would also like tothank the Aspen Center for Physics for its hospitality duringthe completion of this research.
REFERENCES
Agekyan T. A., Anosova J. P. 1967, Astronomicheskii Zhur-nal, 44, 1261Agekyan T. A., Anosova J. P. 1967, Astrophysics, 19, 66 Anosova J. P. 1969, Astrophysics, 5, 81Anosova J. P., Orlov V. V. 1983, Trudy Astronomical Ob-servatory of Leningrad, 38, 142Anosova J. P., Orlov V. V. 1986, Soviet Astronomy, 30, 380Anosova J. P., Orlov V. V. 1994, Celestial Mechanics andDynamical Astronomy, 59, 327Antognini J. M. O., Thompson T. A. 2015, MNRAS, ac-cepted (arXiv:1507.03593)Antonini F., Chatterjee S., Rodriguez C. L., Morscher M.,Pattabiraman B., Kalogera V., Rasio F. A. 2016, ApJ,816, 65Bahramian A., Heinke C. O., Sivakoff G. R., Gladstone J.C. 2013, ApJ, 766, 136Binney J., Tremaine S., 1987, Galactic Dynamics (Prince-ton: Princeton University Press)Brockamp M., Baumgardt H., Kroupa P. 2011, MNRAS,418, 1308Brown W. R., Geller M. J., Kenyon S. J., Kurtz M. J. 2005,ApJL, 622, L33Brown W.R., Cohen J. G., Geller M. J., Kenyon S. J. 2012,ApJL, 754, L2Burrau C. 1913, Astronomische Nachrichten, 195, 113Cohn H., Kulsrud R. M. 1978, ApJ, 226, 1087Davies M. B., Blackwell R., Bailey V. C., Sigurdsson S.1998, MNRAS, 301, 745Demarque P., Virani S. 2007, A&A, 461, 651Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F.A. 2004, MNRAS, 352, 1Fregeau J. M., Ivanova N., Rasio F. A. 2009, ApJ, 707,1533Geller A. M., Mathieu R. D., Harris H. C., McClure R. D.2008, AJ, 135, 2264Geller A. M., Hurley J. R., Mathieu R. D. 2013, AJ, 145,8Geller A. M., Leigh N. W. C. 2015, ApJL, 808, 25Goldreich P., Lithwick Y., Sari R. 2004, ARA&A, 42, 549Heggie D. C., Hut P. 2003, The Gravitational Million-BodyProblem: A Multidisciplinary Approach to Star ClusterDynamics (Cambridge: Cambridge University Press)Heggie D. C. 1975, MNRAS, 173, 729Henon M. 1969, A&A, 1, 223Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A. 2005,MNRAS, 363, 293Hut P. 1983, ApJL, 272, L29Hut P., Bahcall J. N. 1983, ApJ, 268, 319Hut P., Verbunt F. 1983, Nature, 301, 587Hut P., Murphy B. W., Verbunt F. 1991, A&A, 241, 137Hut P., McMillan S., Romani R. W. 1992, ApJ, 389, 527Jeans J. H. 1928, Astronomy and Cosmogony. Cambridge:Cambridge University PressKushnir D., Katz B., Dong S., Livne E., Fernandez R. 2013,ApJ, 778, 37Leigh N., Knigge C., Sills A. 2007, ApJ, 661, 210Leigh N., Sills A. 2011, MNRAS, 410, 2370Leigh N., Geller A. M. 2012, MNRAS, 425, 2369Leigh N., Geller A. M. 2013, MNRAS, 432, 2474Leigh N., Giersz M., Webb J. J., Hypki A., De Marchi G.,Kroupa P., Sills A. 2013, MNRAS, 436, 3399Leigh N., Geller A. M. 2015, MNRAS, 450, 1724Leigh N. W. C., Shara M. M., Geller A. M. 2016, MNRAS,459, 1242Leonard P. J. T. 1989, AJ, 98, 217 Leigh N. W. C., Stone, N. C., Geller A. M., Shara M. M., Muddu H., Solano-Oropeza D., Thomas Y.