Small-N collisional dynamics IV: Order in the realm of not-so-small-N
Nathan W. C. Leigh Aaron M. Geller Michael M. Shara Lukas Baugher Vianny Hierro De'Andre Ferreira Elizabeth Teperino


firstpage–lastpage 2016

Combined Mutiplicative-Heston Model for Stochastic Volatility
    [
    
=============================================================






In this paper, the fourth in the series, we continue our study of combinatorics in chaotic Newtonian dynamics.  We focus once again on the chaotic four-body problem in Newtonian gravity assuming finite-sized particles, and interactions that produce direct collisions between any two particles.  Our long-term goal is to predict the probability of a given collision event occurring over the course of an interaction, as a function of the numbers and properties of the particles.  In previous papers, we varied the number of interacting particles, as well as the distributions of particle radii and masses.  Here, we refine the methods developed in these preceding studies, and arrive at a final and robust methodology that can be used to study collisional dynamics in a variety of astrophysical contexts, ranging from stars in star clusters, galaxies in galaxy groups and clusters and even the collisional growth of planetesimals in protoplanetary disks.  We further present and refine the concept of a Collision Rate Diagram (CRD), the primary tool we use to quantify the relative rates for different collision scenarios to occur.  The agreement between our final theoretical model and the results of numerical scattering simulations is excellent.








gravitation – binaries (including multiple): close – globular clusters: general – stars: kinematics and dynamics – scattering – methods: analytical.




§ INTRODUCTION
 















Particle-particle collisions during small-number fewbody interactions are the cause for several ubiquitous astrophysical phenomena.  These include, but are not limited to, blue straggler formation in globular and open star clusters due to stellar collisions <cit.>, the production of anomalously blue stars in galactic nuclei <cit.>, the formation of intermediate-mass black holes via runaway stellar collisions that could also serve as the seeds for the formation of supermassive black holes <cit.>, the collisional growth of protoplanetary disks <cit.>, the production of runaway stars from young star-forming regions <cit.>, the origins of elliptical galaxies <cit.>, etc. 




Here, we continue our study of direct collisions between particles during chaotic few-body interactions.  We paraphrase the results of our previous works in this series here for completeness.  In Paper I <cit.>, we studied how the collision probability depends on the number of interacting particles.  We found a connection between the mean free path approximation and the binomial theorem. We showed that, for identical particles and a given total encounter energy and angular momentum, the collision probability scales roughly as N^2, where N is the number of interacting particles.  The physical origin of this N-dependence comes from the binomial theorem; the number of ways of selecting any pair of particles from a set of N identical particles is N  2 = N(N-1)/2.

In Paper II <cit.>, we found that, for (near-)identical mass particles, the collision probability is directly proportional to the collisional cross-section for the types of small-number interactions expected to occur in actual star clusters.  The dynamics of such gravitationally-bound systems of chaotically-interacting finite-sized particles are analogous to a system of pendulums; the particles oscillate semi-periodically about the system centre of mass.  Here, the cross-section for any two particles to collide directly is, to first order for particles with similar masses and large radii, proportional to the square of the sum of their radii.  By means of a combinatorics-based approach, it follows that the collision probability can be expressed analytically for any number of particles and any combination of particle radii.

In Paper III <cit.>, we derived analytic formulae for the time-scales for different collision scenarios to occur, and compared the results to numerical scattering simulations of binary-binary interactions.  We showed that the simulated relative probabilities for the different collision scenarios are bounded by the corresponding analytic predictions, assuming either purely radial or purely tangential motions for the particles.   We further showed that, in the purely radial limit, our analytic time-scales provide good order-of-magnitude estimates for the mean time-scales for direct collisions to occur in our simulations.

In this paper, the fourth in the series, we study the probabilities for different collision scenarios to occur, while simultaneously varying the distribution of particle masses and radii.  We first describe the framework underlying our model in the Newtonian limit, which is founded on a combinatorics-based backbone and is designed to calculate the time-scales or rates for direct collisions to occur during chaotic gravitational interactions involving finite-sized particles with different (but comparable) masses.  



In Section <ref>, we apply the mean free path approximation to derive theoretical collision time-scales and relative collision probabilities.  We further introduce the Collision Rate Diagram (CRD), which illustratively quantifies how well our derived collision time-scales are able to reproduce the simulated data.  We also describe the simulations used in this study to test our model.  In Section <ref>, we present and compare the resulting simulated and theoretical collision probabilities and rates, using Collision Rate Diagrams as a guide to constrain to constrain the dominant physics needed to be reproduced via our analytic model.  The assumptions and limitations underlying our model are discussed along with their applicability to astrophysical systems in Section <ref>.  Our key results are summarized in Section <ref>.
 
 

§ METHOD
 

In this section, we re-visit the concept of a Collision Rate Diagram (CRD), first presented in Paper III of this series. In this previous study, we derived different collision time-scales, and from these the various rates and probabilities for different collision scenarios to occur.  This will ultimately facilitate our ability to quantify in this paper, via the CRD, the effects of incorporating different assumptions in deriving these time-scales and rates.  We go on to present the numerical scattering experiments of binary-binary encounters involving finite-sized particles with different combinations of particle masses and radii.  We compare the results of these simulations to our analytic predictions in Section <ref>, for different input assumptions to our model (e.g., setting the collisional cross-section equal to the geometric cross-section, setting it equal to the gravitationally-focused cross-section, with/without the assumption of time-averaged virial equilibrium, etc.), in order to identify the dominant physical processes deciding the relative collision rates and probabilities.  

Throughout this paper, we define a direct collision as occurring when the particle radii overlap directly, following the "sticky-star" approximation.  
 


 §.§ Model
 

In this section, we first re-visit the concept of a Collision Rate Diagram, before going on to present the numerical scattering experiments performed in this paper.  Later, these will be used, in conjunction with the CRDs presented in this section, to identify the dominant physical processes deciding the relative collision rates as a function of the the distribution of particle masses and radii, etc.



  §.§.§ Collision Rate Diagram
 

In this section, we re-introduce the concept of a Collision Rate Diagram (CRD), first presented in Paper III of this series.  This diagram provides an immediate and visual comparison between the predictions of our analytic derivations for the relative rates of collisions between different particle types (and their underlying assumptions; see Paper III) and the results of numerical scattering experiments.  This is because the area corresponding to a particular collision event is directly proportional to the probability of that outcome occurring.  Hence, as we will show, it provides a fast and efficient means of comparing theoretical predictions to simulated data.  By changing the assumptions underlying a given model, the parameter space indicating the rate dominance for each collision scenario will change.  Hence, in this way, the CRD is a potentially efficient tool for isolating the dominant physics deciding the relative rates for different collision scenarios to occur, and can be robustly applied to any astrophysical problem that touches upon the collisional regime of gravitational dynamics.

We begin by describing the original CRD from <cit.>, and repeat it here for completeness.  First, consider interactions involving three different types of particles, labeled A, B, and C.  In this case, we can use our derived time-scales to construct a Collision Rate Diagram, using a similar formalism as outlined in Paper III, and earlier in <cit.> and <cit.>.  A CRD is a diagram that illustrates the parameter space for which the rates of the different types of collisions (e.g., A+A, A+C, B+C, etc.) each dominate over all others.  In Paper III of this series, we considered only 3-dimensional CRDs, which have only a single quadrant, facilitated by writing the fraction of single stars as the sum of the fractions of binaries and triples, combined with the critical assumption of mass conservation for the entire system.

The procedure for producing such a CRD is as follows.  First, we note that, for three particle types, we can write the total number of particles N involved in an interaction as:

    N = N_ A + N_ B + N_ C,

where N_ A, N_ B and N_ C denote, respectively, the total number of particles of type A, B and C.  Then, the fraction of objects of a given particle type i can be written:

    f_ i = N_ i/N,

and the sum of their total must of course satisfy the relation:

    1 = f_ A + f_ B + f_ C


Now, to produce a CRD, every pair of collision rates (e.g., Γ_ A+A, Γ_ A+B, Γ_ A+C, etc.) should be equated, and the resulting relation plotted in f_ B-f_ C-space.  The region of parameter space in the f_ B-f_ C-plane for which each type of collision occurs at the highest rate can then be identified, and a corresponding boundary can be drawn in the CRD.  This produces a diagram that identifies the parameter space in the f_ B-f_ C-plane for which the rates of the different types of collisions (e.g., A+A, A+C, B+C, etc.) each dominate.  

Figure <ref> shows an example of a 3-D Collision Rate Diagram.  To construct this figure, we assume for simplicity that the rate of collisions between particles of type i and j can be written:

    Γ_ i+j = f_ if_ jN_ in_ jσ_ i+jv_ i+j,

where we do not yet include gravitational-focusing in our estimate for the collisional cross-section σ_ i+j, and instead set it equal to the geometric cross-section for collision.  We emphasize that this is the rate for any particle of type i to collide with any particle of type j.  This is in contrast to the rate for a given or specified object of type i to collide with any particle of type j, which is smaller than the previous rate by a factor N_ i.


Equation <ref> can be modified to generate the simplest possible form for this CRD.[We note that Equation <ref> is itself not strictly correct, since the relative rates should include a combinatorial correction to avoid over-counting collisions between identical particles.  We will return to this later, and properly include this correction in our formulation.]  Specifically, we can assume the same particle number density n = n_ j and relative velocity at infinity v = v_ i+j for each type of collision (i.e., for all i and j).  This simplifying assumption neglects the more complicated geometry considered in Paper III, but probes the simplest physical assumptions possible for the problem at hand (i.e., a suitable starting point for the method).  Hence, the simplified CRD shown in Figure <ref> quantities only the importance of the particle number and the geometric collisional cross-section in determining the relative collision rates.  




Now, let us add a fourth particle type in to the mix, with label D.  Here, we can generate a 4-dimensional CRD, which has four quadrants and ultimately represents a 2-dimensional slice of the over-arching 4-dimensional parameter space.  In each of the four quadrants, we set the number of one particle type to be zero.  Then, each quadrant is effectively analogous to the CRD shown in Figure <ref>.  Combining all four quadrants allows for a more thorough comparison between theoretical predictions and the simulations for a larger subset of the total possible parameter space.  An example of a 4-D CRD is shown in Figure <ref>, adopting the same assumptions as in Figure <ref> in order to generate the simplest form of the CRD.



In the subsequent sections, we will compare the results of numerical scattering simulations to our analytic predictions.  This will be done using various forms of the CRD as our guide toward isolating the dominant physics deciding the relative collision rates or probabilities.  In this section, we have presented an over-simplified form of the CRD.  By introducing additional physics in to our model in the subsequent sections, specifically the gravitationally-focussed cross-section and a combinatorial correction, we will illustrate and quantify the effects of each of these physical components on our over-arching model.  



 §.§ Numerical scattering experiments
 

As in Paper III of this series, we calculate the outcomes of a series of binary-binary (2+2) encounters using the  numerical 
scattering code[The source code can be found at http://fewbody.sourceforge.net.].  As discussed in more detail in <cit.>, the code 
integrates the usual N-body equations in position-space in order to advance the system forward in time.  This is done using the 
eighth-order Runge-Kutta Prince-Dormand integration method with adaptive time-stepping and ninth-order error estimation.  


In this paper, we set m_ A = 1 M_⊙, m_ B = 2 M_⊙, m_ C = 3 M_⊙ and m_ D = 4 M_⊙, with R_ A = 1 R_⊙, R_ B = 2 R_⊙, R_ C = 3 R_⊙ and R_ D = 4 R_⊙.  We then consider different combinations of these four over-arching particle types.  

For these simulation sets, all particles are assumed to have finite radii (i.e., spherical) and we adopt the indicated combinations of masses and radii in Table <ref>.  In all simulations all binaries have a_ A = a_ B = 5 AU initially, and eccentricities e_ A = e_ B = 0.  We set the impact parameter to zero and the initial relative velocity at infinity v_ rel to 0.3v_ crit, where v_ crit is the critical velocity.  It is defined as the relative velocity at infinity required for a total encounter energy of zero.[Note that this choice of relative velocity is typical for dense star clusters.]  As found in previous studies <cit.>, such low relative velocities at infinity and small impact parameters maximize the probability of long-lived resonant interactions occurring, for which the assumption of ergodicity is upheld.  All angles defining the relative configurations of the binary orbital planes and phases are chosen randomly.  We perform 4 × 10^4 numerical scattering experiments for every combination of particle masses.  


As in previous papers in this series, all simulations are terminated at the instant the first collision occurs.  If no collisions occur, we use the same criteria as in <cit.> to determine when an encounter is complete.  We refer the reader to <cit.> and the previous papers in this series for the precise implementation of the stopping criteria used in this paper.  






As in previous studies, we adopt a tidal tolerance parameter 
δ = 10^-7 for all simulations.  The reader can refer to previous papers in this series for the full justification underlying this choice for δ (see also  and ).



§ RESULTS
 

In this section, we present the results of our numerical scattering experiments and compare them to our theoretical predictions.  This is done by plotting the results of our numerical experiments in different manifestations of the Collision Rate Diagram, each time incorporating additional physics in to the model that generates the CRD.  These results are summarized below in Table <ref>, and illustrated in Figures <ref>,  <ref> and <ref>.  



 §.§ Confronting the analytic rates with simulated data
 

In order to compare the results of the simulated data with the predictions of the CRD, we can include points in Figure <ref> for those combinations of f_ A, f_ B, f_ C and f_ D included in Table <ref> via our simulations.  The points are assigned shapes according to the outcome specified in each figure.  If the relative collision probabilities in the simulations agree exactly with the analytic predictions shown in the CRD, we plot the points as filled (i.e., the points fall in a region of the CRD corresponding to the dominant outcome found in our simulations).  If the points fall in an incorrect region of the CRD, we leave them unfilled.  




Figure <ref> clearly shows that the agreement between the simulations and the predictions of the CRD is poor without including additional physics, such as gravitational-focusing and a combinatorial correction.  If we include these additional effects in our model, can we improve upon the reported disagreement between theory and simulations?  This is directly addressed below.

To address this, we re-construct Figure <ref>, but adopting the gravitationally-focused cross-section for collision instead of the geometric cross-section.  This is motivated by the derived collision time-scales and rates presented in Paper III of this series.  As shown in Figure <ref>, the net effect of these alternations to the CRD is to increase the probability of collisions involving heavier particles.  The inclusion of this additional physics in our model only slightly improves the agreement between theory and simulations.






Finally, we include a combinatorial correction in our estimates for the relative collision rates in Equation <ref>.  Specifically, the number of ways of selecting two identical particles from a sample of N identical particles is N  2.  To account for this, we directly correct our assumed number fractions in Equation <ref> if i = j, and leave the rates unchanged otherwise.  This corrects for over-estimating the rate of collisions between identical particles in our simplified formulation presented in Section <ref>.  

As shown in Figure <ref>, this additional correction to our base model further improves the agreement between theory and simulations.  Now, all data points agree with our theoretical model.  It follows that our theoretical model successfully reproduces the simulated data for all of our simulations.





 §.§ The validity of the CRD for unexplored regions of parameter space
 

We consider in total 12 combinations of particle types (A, B, C or D), but many more are also possible.  This choice was done for simplicity and to ensure that each simulated data point would appear only once in the CRD, in one of the four quadrants shown in Figure <ref>.  We do not include combinations of only two particle types, for example, since these simulated data points would be degenerate, and appear in more than one place in the CRD (i.e., two different quadrants).

With that said, we did perform a few additional scattering experiments assuming only two types of particles in a given four-body interaction.  Upon comparing the results of these simulations to our model and plotting them in the CRD, we consistently find excellent agreement between the theoretical predictions (i.e., CRD) and the simulations.  For example, we performed analogous simulations to those described in Section <ref> but assuming two particles of type B and two particles of type D.  The simulated data predict that B+D collisions should be dominant for f_ B = f_ D = 0.5.  A simple comparison to Figure <ref>, and specifically the lower right and left panels, reveals that the simulation outcomes are once again successfully described by the CRD.   Upon assuming two particles of type B and two particles of type C, the simulated data predict that B+C collisions should be dominant for f_ B = f_ C = 0.5.  Once again, a simple comparison to Figure <ref>, and specifically the upper and lower right panels, confirms that the simulation outcomes agree with the CRD.  We also point that, for any set of simulations assuming all identical particles, the resulting data points will always appear at the extremums of the various axes in Figure <ref>.  Consequently, by definition, these simulated data points will also always agree with the predictions of the CRD.  All of these examples further supports the validity the CRD developed and tested in this paper.






 










§ DISCUSSION
 

In this paper, we present a method for directly comparing to simulated data the analytic rates for different collision scenarios to occur during chaotic gravitational interactions involving arbitrary numbers of finite-sized particles, with any distributions of particle masses and radii.  The method is flexible in the sense that the assumptions underlying the derivations of the analytic rates can be freely modified, and the subsequent impact on the agreement with the simulated data is easily quantified.  This facilitates a quick identification of the dominant physics deciding the relative rates of collisions, for the chosen initial conditions and particle properties.  To illustrate this, we compare our analytic predictions to the results of numerical scattering experiments of four-body interactions involving finite-sized particles.  Overall, the agreement between our analytic predictions and the simulations is excellent.  

This paper is meant as a more in-depth introduction to the CRD than given in Paper III of this series, while clearly articulating how to properly use it to construct a robust analytic framework for predicting collision rates, depending on the particle types and important physics (sticky-star approximation for collisions?  dissipation?  GW emission? etc.) for the problem to which the user wishes to apply their version of the CRD.  In other words, the validity of the CRD depends on the astrophysical environment of interest, since this in turn decides the types and properties of the particles (e.g., stars, comets, galaxies, etc.) in addition to the dominant physics affecting the collision rates (e.g., gravitational-focusing, combinatorics, various forms of dissipation such as from tides or GWs, etc.).





 §.§ Caveats
 

A few cautionary notes should be kept in mind when applying the Collision Rate Diagrams presented in this paper.  First, we expect our base model to be valid only for mass ratios ≲ 10 <cit.>.  Above this, we expect the most massive objects to sufficiently dominate the gravitational potential, such that the collision rate begins to enter the loss-cone regime (see <cit.> for more details) and/or the lowest mass particles are quickly ejected from the system without entering a resonant interaction state.  In general, a prolonged resonant state and the assumption of ergodicity being upheld are key requirements for our base model to be applicable <cit.>.  More work will be needed to better understand how to adapt our model to smoothly transition between these two regimes, where these assumptions start to break down.  Provided the assumption of ergodicity is upheld, we expect some variation of the collision rate estimates presented in Paper III to accurately capture the physics, in all but the most extreme cases (see below).

Finally, we note that more work is still needed to better understand the dependence of the collision probability on the impact parameter of the encounter, which was consistently set to equal zero throughout this paper.  An increase in the impact parameter should mostly act to increase the total angular momentum of the interaction (for a fixed relative velocity at infinity).  In future work, we intend to vary this parameter to better understand what adjustments to our base model might be needed to properly accommodate this additional free parameter.  The analytic model presented in Paper III of this series for the high-angular momentum regime, which assumes purely tangential motions relative to the system centre of mass, could provide the needed adjustments.




 §.§ Implications for dense stellar environments
 

In this section, we describe the implications and applicability of our results to real astrophysical environments, with a focus on dense stellar clusters.



  §.§.§ Old globular clusters
 

Old globular clusters are the ideal environments for applying the methods presented in this paper.  There are three main reasons for this: 1) direct single-binary and binary-binary encounters (and possibly interactions involving triples as well) occur commonly in the dense cores of old GCs <cit.>; 2) direct collisions between stars occur frequently during 3-body, 4-body, 5-body, etc. interactions in GCs <cit.>; and 3) the range of stellar masses among such old stellar populations are typically ∼ 0.08 - 0.8 M_⊙, which limits the mass ratios of the fewbody interactions to q ≳ 0.1 and consequently ensures that the assumption of ergodicity (which is needed to apply our model and construct a CRD) should be approximately upheld for most chaotic fewbody interactions in old GCs.

To properly apply the CRD and methods presented in this paper to old GCs, a mass-radius relation for main-sequence stars should be adopted.  Typically, for such an old stellar population, the relation R/R_⊙ = (m/M_⊙)^0.75 is used, where m is the main-sequence mass and R is the corresponding stellar radius <cit.>.  This ensures that the particle masses and radii are chosen appropriately when constructing a CRD designed for chaotic fewbody interactions in old GCs.  Such a change would likely have the affect of reducing the importance of the particle radius in the CRD, such that the difference in area between the largest and smallest zones (e.g., D+D and A+A) would be reduced.



  §.§.§ Young massive clusters
 

Young massive clusters could also be well suited to applying the methods presented in this paper.  For example, <cit.> showed using N-body simulations that runaway collisions of massive stars can occur in the dense cores of primordial GCs, possibly forming a supra-massive star or even an intermediate-mass black hole.  We emphasize that, since we are only interested in the relative collision rates, the model is scale-free for our purposes if we fix the ratio of particle mass to particle radius such that it remains unity (as assumed in Figure <ref>, for example).  By extension, the CRD shown in Figure <ref> is equally well-suited to larger, more massive particles, provided the masses are all directly proportional to the masses assumed in Figure <ref>.  For example, if we replace the particles with masses of 1, 2, 3 and 4 M_⊙ with more massive particles with masses of 10, 20, 30 and 40 M_⊙, then Figure <ref> still applies provided the stellar radii are, respectively, 10, 20, 30 and 40 R_⊙.  Thus, our method is just as applicable to treating chaotic interactions involving more massive stars, provided the minimum mass ratio in each interaction satisfies q ≳ 0.1, as previously discussed.



§ SUMMARY
 

In this paper, the fourth in the series, we push forward in our study of chaotic Newtonian gravity involving small numbers of finite-sized particles.  Our focus remains direct collisions between pairs of particles in the "sticky-star" approximation.  Our over-arching goal in this series of papers is to develop a method to calculate the probability of any two particles colliding during a chaotic (bound) resonant gravitational interaction involving any number N of particles with any combination of particle masses and radii, as well as to directly compare its predictions to simulated data.  

In our previous papers, we showed that (1) the probability of a collision occurring during interactions involving identical particles is approximately proportional to N^2, which comes from combinatorics and, specifically, the number of ways of selecting any two particles from a larger set of N identical particles, or N  2 <cit.>; (2) for strongly bound gravitational encounters (i.e. with E ≪ 0, where E is the total encounter energy, and having small impact parameters) involving small numbers of particles, the collision probability is directly proportional to the collisional cross-section.  For identical particle masses and large particle radii, the collisional cross-section is roughly equal to the sum of the cross-sectional areas of the colliding particles <cit.>; and (3) for different particle masses but identical particle radii, the mean free path approximation can be used in conjunction with (or without) the assumptions of time-averaged virial equilibrium and energy equipartition to derive estimates for the relative collision rates that agree with the simulated data at the order-of-magnitude level <cit.>.  

In this paper, we continue our study by considering interactions involving particles with both different masses and radii.  As in Paper III, we consider the four-body problem in this paper.  This is because, for N = 4, we can run more simulations since we minimize the computational expense due to the small number of particles.  This contributes to a significant increase in the statistical significance of the analysis.  Using our previous results from Papers I, II and III, we derive Collision Rate Diagrams for our case of interest.   This is done first assuming that only the geometric cross-sections for collision affect the relative collision rates, and then again assuming gravitational focusing and including a combinatorial correction.  For these cases, we analyze the relative collision probabilities as predicted by our analytic formulae and compare them to the results of numerical scattering simulations performed with the  code <cit.>.  

By calculating different manifestations of the CRD, our method illustrates the following key result.  While the geometric cross-section captures most of the relevant physics for comparable particle masses, invoking the additional assumptions of a gravitationally-focused cross-section along with a combinatorial correction yield better agreement with the simulations.  With these additional assumptions, our analytic estimates reproduce the simulated data for all combinations of particle masses and radii considered here.  


Our method (or variations thereof) is suitable to direct stellar collisions in dense star clusters, the collisional growth of planetesimals in protoplanetary disks, the growth of super-massive black holes via runaway stellar collisions, the formation of giant elliptical galaxies in the galaxy clusters and groups, and even the growth of large-scale structure in the Universe.



§ ACKNOWLEDGMENTS


The authors thank an anonymous reviewer for comments and suggestions that strengthened our manuscript.  N. W. C. L. gratefully acknowledges support from the American Museum of Natural History and the Richard Guilder Graduate School, specifically the Kalbfleisch Fellowship Program, as well as support from a National Science Foundation Award No. AST 11-09395.  



99


















[Bailey & Davies1999]bailey99  Bailey V. C., Davies M. B. 1999, MNRAS, 308, 257
[Balland et al.1998]balland98 Balland C., Silk J., Schaeffer R. 1998, ApJ, 497, 541
[Binney & Tremaine1987]binney87 Binney J., Tremaine S., 1987, Galactic Dynamics (Princeton: Princeton University Press)
[Blaauw & Morgan1954]blaauw54 Blaauw A., Morgan W. W. 1954, ApJ, 119, 625










[Dale et al.2009]dale09 Dale J. E., Davies M. B., Church R. P., 
Freitag M. 2009, MNRAS, 393, 1016
[Davies et al.1998]davies98 Davies M. B., Blackwell R., 
Bailey V. C., Sigurdsson S. 1998, MNRAS, 1998, 301, 745


[Fregeau et al.2004]fregeau04
  Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A. 2004,
  MNRAS, 352, 1






[Giersz et al.2015]giersz15 Giersz M., Leigh N. W. C., Hypki A., Luetzgendorf N., Askar A. 2015, MNRAS, 454, 3150
[Geller & Leigh2015]geller15 Geller A. M., Leigh N. W. C. 2015, ApJL, 808, 25
[Goldreich, Lithwick & Sari2004]goldreich04 Goldreich P., 
Lithwick Y., Sari R. 2004, ARA&A, 42, 549 





[Hansen, Kawaler & Trimble2004]hansen04 Hansen C. J., Kawaler S. D., Trimble V. 2004, Stellar interiors : physical principles, structure, and evolution (New York: Springer-Verlag) 
[Hypki & Giersz2016]hypki16 Hypki A., Giersz M. 2016, MNRAS, submitted (arXiv:1604.07054)
[Hypki & Giersz2017]hypki17 Hypki A., Giersz M. 2017, MNRAS, 466, 320














[Leigh, Knigge & Sills2007]leigh07 Leigh N.,
Knigge C., Sills A. 2007, ApJ, 661, 210
[Leigh & Sills2011]leigh11 Leigh N.,
Sills A. 2011, MNRAS, 410, 2370
[Leigh & Geller2012]leigh12 Leigh N., Geller A. M. 
2012, MNRAS, 425, 2369
[Leigh & Geller2013]leigh13 Leigh N., Geller A. M.
2013, MNRAS, 432, 2474


[Leigh & Geller2015]leigh15 Leigh N., Geller A. M.
2015, MNRAS, 450, 1724
[Leigh, Shara & Geller2016]leigh16 Leigh N. W. C., Shara M. M., 
Geller A. M. 2016, MNRAS, 459, 1242
[Leigh et al.2016]leigh16b Leigh N. W. C., Antonini F., Stone N. C., Shara M. M., Merritt D. 2016, MNRAS, 463, 1605
[Leigh et al.2017]leigh17 Leigh N. W. C., Geller A. M., Shara M. M., Garland J., Clees-Baron H., Ahmed A. 2017, MNRAS, 471, 1830
[Leigh & Wegsman2018]leigh18 Leigh N. W. C., Wegsman S. 2018, MNRAS, accepted (arXiv:1801.07257)  
[Leonard1989]leonard89 Leonard P. J. T. 1989, AJ, 98, 217


[Lithwick & Chiang2007]lithwick07 Lithwick Y., 
Chiang E. 2007, ApJ, 656, 524 












[Merritt2013]merritt13 Merritt D. 2013, Dynamics and Evolution

[Mikkola1983]mikkola83 Mikkola S. 1983, 
MNRAS, 203, 1107
[Mikkola1984]mikkola84 Mikkola S. 1984, 
MNRAS, 207, 115
[Mikkola & Valtonen1990]mikkola90 Mikkola S., 
Valtonen M. J. 1990, ApJ, 348, 412 

























[Oh et al.2015]oh15 Oh S., Kroupa P., Pflamm-Altenburg J. 2015, ApJ, 805, 92
[Perets & Subr2012]perets12 Perets H. B., Subr L. 2012, ApJ, 751, 133
[Portegies Zwart et al.2004]portegieszwart04 Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W. 2004, Nature, 428, 724 
[Ryu, Leigh & Perna2017a]ryu17a Ryu T., Leigh N. W. C., Perna R. 2017, MNRAS, 467, 4447
[Ryu, Leigh & Perna2017b]ryu17b Ryu T., Leigh N. W. C., Perna R. 2017, MNRAS, accepted (arXiv:1703.08538)
[Ryu, Leigh & Perna2017c]ryu17c Ryu T., Leigh N. W. C., Perna R. 2017, MNRAS, accepted (arXiv:1703.08551)
[Sigurdsson & Phinney1993]sigurdsson93 Sigurdsson S., Phinney E. S. 1993, ApJ, 415, 631



[Shara & Shaviv1974]shara74 Shara M. M., Shaviv G. 1974, Nature, 248, 398




[Stone, Kuepper & Ostriker2017]stone17 Stone N. C., 
Kuepper A. H. W., Ostriker J. P. 2017, MNRAS, 467, 4180








[Trinchieri et al.2003]trinchieri03 Trinchieri G., Sulentic J., Breitschwerdt D., Pietsch W. 2003, A&A, 401, 173
[Valtonen &
    Karttunen2006]valtonen06 Valtonen M., Karttunen H. 2006, The
  Three-Body Problem (Cambridge: Cambridge University Press)
[Yu2003]yu03 Yu Q. 2003, MNRAS, 339, 189












