Small-N Collisional Dynamics: Pushing Into the Realm of Not-So-Small-N

In this paper, we study small-N gravitational dynamics involving up to six objects. We perform a large suite of numerical scattering experiments involving single, binary, and triple stars. This is done using the FEWBODY numerical scattering code, which we have upgraded to treat encounters involving triple stars. We focus on outcomes that result in direct physical collisions between stars, within the low angular momentum and high absolute orbital energy regime. The dependence of the collision probability on the number of objects involved in the interaction, N, is found for fixed total energy and angular momentum. Our results are consistent with a collision probability that increases approximately as N^2. Interestingly, this is also what is expected from the mean free path approximation in the limit of very large N. A more thorough exploration of parameter space will be required in future studies to fully explore this potentially intriguing connection. This study is meant as a first step in an on-going effort to extend our understanding of small-N collisional dynamics beyond the three- and four-body problems and into the realm of larger-N.


INTRODUCTION
The gravitational three-body problem was first studied by Sir Isaac Newton in his Principia (Newton 1686). After solving the two-body problem, Newton boldly added a third body into the mix and attempted to create a similar mathematical formalism to describe the motion of three celestial bodies under their mutual gravitational attraction. The solution evaded Newton until his end, leaving the problem at the mercy of his disciples. It was later taken up by Euler (Euler 1772), Lagrange (Lagrange 1811), Jacobi (Jacobi 1836), Poincare (Poincare 1892), Hill (e.g. Hill 1878), Henon (e.g. Henon 1969), and a host of others, often with a focus on the Earth-Moon-Sun system (e.g. Valtonen & Karttunen 2006).
Despite the considerable efforts of numerous researchers, a simple analytic solution to the general threebody problem has never been found. Sundman (1912) discovered a uniformly convergent infinite series involving familiar functions that technically solves the three-body problem. However, in order to attain a reasonable level of accu-⋆ E-mail: leighn@mcmaster.ca (NL); a-geller@northwestern.edu (AG) racy, the solution must contain on the order of 10 8,000,000 terms (Valtonen & Karttunen 2006). More practical solutions require a number of simplifying assumptions to make the general three-body problem tractable (e.g. Henon 1969). As a result, the most useful analytic solutions tend to be solely applicable to a very narrow subset of the total allowed parameter space.
The introduction of computers within the last few decades revolutionalized our understanding of the threebody problem. This allowed for the direct integration of the equations of motion for each particle over small time steps, incrementally moving each body forward in an iterative fashion. Despite the fact that this approach completely transformed our tool-set for studying the three-body problem, it too has its short-comings. For example, the times required to run the simulations to completion are often quite long. The precise definition of a "complete encounter" can often be ambiguous as well. Quasi-stable configurations can form that remain bound for millions or even billions of years before eventually breaking apart (e.g. Mikkola & Valtonen 1986). There is also the issue of errors in computing the trajectories of the particles in position-space which are introduced at each time-step. These arise as a result of moving one body forward at a time, instead of moving all bod-ies simultaneously. Such errors can be minimized by taking suitably short time-steps. However, this comes at the often considerable cost of simulation run-time.
When only three bodies are involved, qualitative descriptions of the outcomes of dynamical interactions are often relatively straight-forward. They include ionizations, exchanges, fly-bys, and collisions (e.g. Hut & Bahcall 1983;Mikkola 1983Hut 1984;McMillan & Hut 1996;Fregeau et al. 2004). However, the number of possible outcomes quickly increases with increasing N , where N is the total number of objects involved in the interaction (Leigh & Sills 2011). This complicates descriptions of the interactions, and introduces a considerable challenge in developing a physical understanding of the evolution of the system. For example, nearly 100 generic outcomes are possible for encounters involving six objects. Additional bodies not only increase the parameter space to be explored, they also increase integration times for computer simulations. Consequently, most previous numerical scattering studies considered only single-binary and, to a lesser extent, binary-binary encounters (e.g. Heggie 1975;Hut 1983;McMillan 1986;Sigurdsson & Phinney 1993;Davies 1995;Bacon, Sigurdsson & Davies 1996;Fregeau et al. 2004).
Many of the numerical scattering studies conducted to date considered equal-mass particles, and devoted their attention to studying the effects of varying the initial relative velocity or impact parameter (e.g. Hut & Bahcall 1983). For example, Hut (1983) found analytic approximations for high-velocity encounters, and provided simple formulae for the corresponding cross-sections for exchanges and ionizations to occur. Similar analytic formulae were derived by  for encounters involving two binaries having unequal orbital energies but equal masses.
An extensive analysis that considered unequal mass particles was conducted by Sigurdsson & Phinney (1993), who studied the effects of dynamics on the stellar and remnant populations in a globular cluster. A number of other scattering experiments were designed purely to study the formation of different types of stellar exotica in globular clusters, including blue stragglers (e.g. Leonard 1989), low-mass xray binaries (e.g. Sigurdsson & Phinney 1995), cataclysmic variables (e.g. Ivanova et al. 2006), and millisecond pulsars (e.g. Ivanova et al. 2008). Many of these also considered a range of particle masses. The effects of general relativity have also been studied in the context of the three-and fourbody problems. For example, Mikkola & Valtonen (1990) and Valtonen et al. (1994) considered interactions involving binary super-massive black holes during the mergers of galaxies, and identified several important trends arising from energy losses due to gravitational radiation.
In this paper, we study gravitational interactions involving up to six objects. Our goal is to better understand how the outcome of an encounter depends on the number of interacting objects. This is a first step toward extending techniques developed for the three-body problem, where the vast majority of research efforts have thus far been concentrated, to treat larger-N encounters. To this end, we perform > 10 7 numerical scattering experiments involving single, binary and triple star systems. The number of possible encounter outcomes is a steeply increasing function of N . This presents a considerable challenge when trying to draw parallels between encounters involving different numbers of objects. To minimize this issue, we focus only on outcomes resulting in direct physical collisions, which we consider to have occurred when the radii of any two stars overlap.
In Section 2, we describe the set-up for our numerical scattering experiments, including the range of initial conditions considered. Additionally, we develop an analytic formula for the collision probability as a function of N , that is based on a simple analytic model originally derived for 1+2 encounters. In Section 3, we present the results of our analysis of this large suite of single-binary (1+2), binarybinary (2+2), single-triple (1+3), binary-triple (2+3), and triple-triple (3+3) scattering experiments. Here, we fit the analytic formula to the results from these numerical scattering experiments, thereby deriving the N -dependence of the collision probability. The implications of our analysis for small-N collisional dynamics are discussed in Section 4. Concluding remarks are given in Section 5. Finally, in an appendix, we present a formalism for creating schematic diagrams that quantitatively depict the evolution of an interaction in energy-space, and briefly discuss some of their possible applications.

METHOD
In this section, we describe the details of our numerical scattering experiments, and present the general form for the collision probability to which our results will be compared in Section 3.

Scattering Experiments
We calculate the outcomes of a series of single-binary, binary-binary, single-triple, binary-triple, and triple-triple encounters using the FEWBODY 1 numerical scattering code. The code integrates the usual N -body equations in configuration-(i.e. position-) space in order to advance the system forward in time. For details concerning the adaptive integration, classification techniques, etc. used by FEW-BODY, we refer the reader to Fregeau et al. (2004).
We adapted the FEWBODY code to handle all encounters involving singles, binaries, and triples. (Specifically we created additional subroutines to simulate 1+3 and 3+3 encounters 2 ; codes to simulate encounters between binaries and singles only, as well as a 2+3 encounter code, were previously available in the FEWBODY package.) We use the same criteria as Fregeau et al. (2004) to decide when a given encounter is complete, defined as the point at which the separately bound hierarchies that make up the system are no longer interacting with each other or evolving internally.
To perform physical collisions between stars, FEW-BODY uses the "sticky-star" approximation. This treats stars as rigid spheres with radii equal to their stellar radii. A physical collision is assumed to occur when the radii of the stars overlap. When this happens, the stars are merged assuming conservation of linear momentum and no mass loss. This does not account for tidal effects, which could significantly increase the collision probability (e.g. McMillan 1986), but are beyond the scope of this work. For this reason, we consider the collision probabilities presented in this paper to be lower limits for the true values.
Previous scattering experiments have shown that the total energy and angular momentum are the most important parameters in deciding the outcomes of 1+2 interactions (e.g. Valtonen & Karttunen 2006). Therefore, we will fix the total energy and angular momentum when comparing between encounters involving different numbers of objects. By fixing these quantities, we aim to remove the dependences of the encounter outcomes on energy and angular momentum (when comparing between different encounter types), thereby normalizing the comparisons to reveal the N -dependence of the collision probability. We consider several different combinations of the total energy and angular momentum, which we henceforth refer to as Runs (see below). For each of these combinations or Runs, both the total energy and angular momentum are always chosen to be the same to within a factor of ∼ 2 for every type of encounter.
We focus on three primary Runs (labeled Runs 1, 2, and 3 in Table 1) in Section 3. However, in order to check the robustness of our results, we also perform eight additional Runs with similar total energy and angular momentum as adopted in Runs 1, 2, and/or 3. We selected the following additional combinations of semi-major axes for the two orbits used in each of these additional Runs: 0.5 AU and 3.5 AU, 0.5 AU and 6.5 AU, 0.5 AU and 8 AU, 0.5 AU and 10 AU, 1 AU and 8.5 AU, 1 AU and 11.5 AU, 1 AU and 15 AU, 2 AU and 15 AU. For all three primary Runs, we perform a total of 800,000 numerical scattering experiments for each of the different encounter types, whereas this number is reduced to 200,000 for the other Runs. In total, this amounts to 12 million simulations. This number is chosen to be a balance between statistical significance and simulation run-time, which can be quite long for simulations involving five or more objects.
As we design our experiment such that the total angular momentum is roughly the same for every encounter type within a given Run, we choose the sum of the semi-major axes of the two objects involved in the interaction (which we will call the geometric cross-section) to be equivalent to within a factor of ∼ 2 for every encounter type. These cross-sections are determined by the initial semi-major axes of any binaries and/or triples involved in the interaction, since the physical radii of the stars are small in comparison. (For triples, the cross section is primarily determined by the semi-major axis of the outer orbit.) Specifically, we choose the semi-major axes shown in Table 1. A semi-colon is used to separate different objects, whereas a comma is used to separate the orbits within triples. Parantheses are used to enclose the semi-major axes of triples, with the smaller of the two semi-major axes always corresponding to the inner binary.
We assume equal masses of 1 M⊙ and stellar radii of 1 R⊙ for all stars, and circular orbits for all binaries and triples. This is done to make our exploration of the relevant parameter space more tractable. However we expect the collision probability to in general depend on the mass ratios and orbital eccentricties (e.g. Sigurdsson & Phinney 1993). Unequal masses and non-zero eccentricities are be-yond the scope of the present paper. However we intend to address these issues in future work. The ratio between the inner and outer semi-major axes is chosen to be 7 for all triples, since this roughly corresponds to the stability boundary for the internal evolution of triples composed of equal mass-particles with initially zero-eccentricity inner orbits and moderate-eccentricity outer orbits (Mardling 2001). The initial angle of inclination between the inner and outer orbits of all triples, in addition to their initial relative phases, are chosen at random.
For each Run, we populate a grid of scattering experiments varying only the relative velocity at infinity and the impact parameter. Specifically, we select values for the relative velocity at infinity from 0 to 1.1 in equally-spaced intervals of 0.004, in units of the critical velocity vcrit (defined as the relative velocity that gives a total energy of zero for the encounter). We select values for the impact parameter from 0 to 7 in equally-spaced intervals of 1, in units of the geometric cross-section for collision (e.g. the semi-major axis of the binary for a 1+2 encounter, the sum of the semi-major axis of the outer orbit of the triple and that of the binary in a 2+3 encounter, etc.). In this way, for a given impact parameter, we ensure that the range in both the total angular momentum and the total energy are equal for every encounter type to within a factor of ∼ 2. Finally, at each point in this grid, we perform multiple scattering experiments randomizing all angles in the encounter, including the angles at impact between the orbital planes of any binaries and triples involved in the encounters.

Collision Probability
In this section, we present a general functional form for the collision probability. We begin by comparison to the 1+2 encounters studied in Section 8.6 of Valtonen & Karttunen (2006), and then generalize this formula to N > 3. The bestfitting values for all free parameters in this analytic collision probability formula are then found in Section 3 through fits to the results of our numerical scattering experiments.

Collision Probability for N = 3
The simple analytic model we use to guide our choice for the collision probability was originally found for resonant 1+2 interactions. The model is described in detail in Section 8.6 of Valtonen & Karttunen (2006), to which we refer the reader for its full derivation. It is based on the assumption that the encounter will evolve via a series of successive ejections, eventually leading to the escape of one of the stars from the three-body system. We use the term ejection to describe a scenario in which one object ends up at a great distance from the others before falling back to resume the interaction. Some of these ejections are of considerably longer duration than others. Every ejection counts, however, since they all produce a temporary binary independent of the duration of the ejection. The temporary binary experiences close two-body encounters at the pericentre q of every orbit (e.g. Szebehely 1967; Saslaw 1974; Valtonen & Karttunen 2006).
For a 1+2 interaction, the probability that two stars  will pass within a distance q from each other can be approximated by: where a0 is the initial semi-major axis of the binary, and we require q ≪ a0. The factor C(L) is needed to account for the fact that the lifetime of the system depends on the total angular momentum L (e.g. Valtonen 1974;Anosova & Orlov 1986). If we take q to be the physical sizes of the objects involved in the encounter, then Equation 1 gives the collision probability for a resonant 1+2 interaction. It has been shown to give good agreement with the results of numerical scattering experiments (Valtonen & Karttunen 2006).

Collision Probability for N > 3
We write the collision probability for an encounter involving N > 3 objects as: where and α, β, δ, and γ are all constants. We take a0 to be the semi-major axis of the shortest-period orbit involved in the interaction, and C(N, L) is a function of both the total angular momentum L and the number of objects N . As we will show in Section 3, Equation 2 gives excellent agreement to the results of our numerical scattering experiments. Below, we justify our choice for this functional form for the collision probability. First, we are interested in quantifying the dependence of the collision probability on the number of objects N involved in the interaction. This is accounted for in Equation 2 by a general power-law dependence on N (i.e. β). Second, based on the results of previous numerical scattering experiments, we expect that an encounter outcome will depend both on the total energy and the total angular momentum. Therefore, we expect that these two quantities should appear somewhere in the equation for the collision probability. These dependences are included in Equation 2 via the terms C(N, L) and q/a0. The first term C(N, L) is a function of the total angular momentum L, whereas the second term is roughly proportional to the total energy of the encounter. This last point follows from the fact that we have defined a0 to be the semi-major axis of the shortest-period orbit, and it is this orbit that has the largest absolute orbital energy. Moreover, previous numerical scattering experiments of 1+2 and 2+2 encounters have shown that the components of the shortest-period binary involved in the interaction are the most likely to experience a collision during an encounter (e.g. Fregeau et al. 2004).
We allow for a possible N -dependence in the function C(N, L), since previous numerical scattering experiments performed to constrain this function considered only 1+2 interactions. Therefore, we do not know whether the total number of stars involved in the interaction will affect its lifetime, and play a role in determining the function C(N, L). As we will show in Section 3, the specific functional form we have adopted for C(N, L) in Equation 3 is needed to ensure that the correct agreement with the results of our numerical scattering experiments persists as we move to larger total angular momentum.
It is important to note that Equation 2 does not apply to 1+2 encounters. Therefore, we do not include them when finding the best-fit parameters. This is because the conditions imposed by our assumptions (e.g. equal masses for all stars, only circular orbits, etc.) are such that 1+2 interactions cannot always be made to fit into our normalization for comparing between the different encounter types without significantly modifying the initial parameters of the encounter. Specifically, 1+2 encounters initially involve only a single bound orbit (via the binary). However, all other encounter types initially involve multiple orbits, which affords us additional free parameters. These can be adjusted to get the total energy and angular momentum to within our required factor of 2, and therefore define a suitable normalization between encounter types. This is not always possible within the confines of our assumptions when 1+2 encounters are also included. We will come back to this issue in Section 4.

RESULTS
The collision probability is calculated from the output of our numerical scattering experiments for a given encounter type and a given Run as: where N coll is the total number of encounters that result in a direct physical collision, and N is the total number of encounters performed. The uncertainty for the collision probability is calculated using Poisson statistics according to: Technically, this uncertainty should include scattering experiments that result in unresolved outcomes (Hut 1983). However, we find that the number of unresolved outcomes is sufficiently small that it does not significantly contribute to ∆P coll , and we do not include it in its calculation.
We show these collision probabilities as a function of the total angular momentum for 1+2, 2+2, 1+3, 2+3, and 3+3 encounters in Figures 1, 2 and 3 for Runs 1, 2, and 3, respectively. The total angular momentum, shown on the x-axis, is provided in units of M⊙σgeovcrit, where σgeo is the geometrical cross-section for collision in AU (which is determined by the semi-major axes of the orbits going into the interaction).
We then fit Equation 2 to each of our individual Runs in order to derive the best-fitting values for the free parameters α, β, δ, and γ, that correctly predict the collision probability simultaneously for the 2+2, 1+3, 2+3, and 3+3 encounters. We use the Levenberg-Marquardt least-squares fitting technique implemented in the MPFIT code (Markwardt 2008;Moré 1978) to derive these best-fitting values. In order to obtain realistic uncertainties on these fit parameters, we scale the uncertainties on our measurements by a constant factor to make the reduced χ 2 value equal to unity.
The solid lines in Figures 1, 2 and 3 show our fits to Equation 2. The best-fitting parameters are shown in Table 2 for all three Runs, along with their corresponding uncertainties.
The mean power-law index for all 11 Runs isβ = 2.40 ± 0.12. We see weak evidence for an increase in our β values with increasing geometric cross-section, with a Pearson's linear correlation coefficient of 0.73. We suspect that this weak trend is due to our normalization beginning to break down at large impact parameters, which can amplify small differences in the geometric cross-sections for different encounters types. The best-fitting power-law index on N is consistent with a value of β = 2 to within roughly three standard deviations for all Runs. Moreover, a chi-squared test between all 11 of our β values and a constant powerlaw index of β = 2 returns no significant distinction, with a reduced chi-squared value of 1.73. Therefore, we conclude that our results are consistent with a collision probability that increases approximately as N 2 .
Additionally, we find an exponential drop-off in the collision probability with increasing angular momentum that is steeper for larger N . Once again, this is the case for all Runs. We will return to these two features of the collision probability relation in Section 4.
We find some small disagreement between the collision probabilities for 2+2 and 1+3 encounters for Run 3. This can be understood as follows. For 1+3 encounters, the geometric cross-section for collision is equal to the semi-major axis of the outer orbit of the triple (since the radius of the single star is negligible in comparison). For 2+2 encounters, however, the geometric cross-section is equal to the sum of the semi-major axes of the two binaries. Therefore, within the framework of our experimental set-up, the approximation that the geometric cross-sections for 1+3 and 2+2 encounters will be equal only holds provided one of the binaries in the 2+2 case has a much greater orbital separation than the other binary. This assumption is a good approximation for Runs 1 and 2, however it begins to break down for Run 3 with increasing impact parameter (i.e. increasing angular momentum in Figure 3).  Figure 1. The probability of a single collision occurring as a function of the total angular momentum is shown for Run 1 (see Table 1) for every type of encounter. The open stars correspond to 3+3 encounters, the solid triangles to 2+3 encounters, the open circles to 1+3 encounters, the solid circles to 2+2 encounters, and the crosses to 1+2 encounters. The solid lines show our best-fits to the data using Equation 2. The top line corresponds to the case N = 6, the line below to N = 5, and the bottom line to N = 4. We show the scaled uncertainties here (and in Figures 2  and 3), as discussed in Section 3.

DISCUSSION
Within the angular momentum and energy regime studied here, we find that the probability of a direct physical collision occurring during an interaction increases roughly as N 2 . One interpretation of this result can be understood as follows. Previous numerical scattering experiments of 1+2 interactions have shown that the collision probability is proportional to the average number of close approaches experienced by the system per crossing time, multiplied by the number of crossing times survived through (Valtonen & Karttunen 2006). Therefore, this N 2 dependence may arise physically from an N 2 dependence in the total number of close approaches experienced by the system. This can occur in one of (at least) two ways. Either the number of close approaches per crossing time scales as N 2 while the number of crossing times survived through is constant, or the number of crossing times the system survives through (until the time when the first collision occurs) scales as N 2 while the number of close approaches per crossing time is constant. We find that the average number of crossing times at the time of the first collision is the same to within roughly a factor of ≈ 2 for all encounter types, and that there is no trend with N . Therefore, it is unlikely that the latter scenario described above is responsible for producing the N 2 dependence. This suggests that perhaps the N 2 dependence 2.25 ± 0.11 δ -6.47e-2 ± 4.40e-3 -4.60e-2 ± 3.21e-3 -1.15e-1 ± 9.97e-3 γ 1.33e-1 ± 1.56e-2 1.33e-1 ± 1.56e-2 2.31e-1 ± 2.33e-2  Figure 2. The probability of a single collision occurring as a function of the total angular momentum is shown for Run 2 (see Table 1) for every type of encounter. The symbols representing each encounter type, as well as the lines showing our best-fits to the data using Equation 2 for each encounter type, are the same as in Figure 1.
in the collision probability may arise from a similar dependence in the number of close approaches per crossing time. However, in practice, this hypothesis is much more difficult to test quantitatively. In particular it is not immediately clear how to define a "close approach". We attempt to quantify this effect by examining animations of a handful of encounters in position-space, and counting the number of close approaches per crossing time by eye. It is clear that the number of close approaches per crossing time increases with increasing N . However, this method is far from adequate to determine the precise N -dependence of this relation at any statistically significant level.
In Appendix A, we present a method that will improve upon this component of our analysis in future studies. Specifically, we present a prescription for generating schematic diagrams that depict the evolution of an interaction in energy-space. As is explicitly shown in Appendix A, the advantage of this technique is to provide a straightforward means of defining the criterion of "close approach" in terms of the fraction of the total energy exchanged between two individual stars. This directly relates the definition of close approach to the total energy and therefore the  Figure 3. The probability of a single collision occurring as a function of the total angular momentum is shown for Run 3 (see Table 1) for every type of encounter. The symbols representing each encounter type, as well as the lines showing our best-fits to the data using Equation 2 for each encounter type, are the same as in Figure 1.
initial conditions of the encounter, which is not the case if this criterion is defined purely in terms of a distance. Also note that, in a sense, a collision could be interpreted as a very strict definition of a "close approach". Here multiple crossing times are required before this definition of close approach is satisfied. As we find that the number of crossing times until the time of first collision is roughly equivalent (to within a factor of ≈ 2) for all encounter types, our result can potentially also be expressed as the number of "close approaches", defined in this manner, per crossing time scaling as N 2 . We will explore a range of other criteria to define a close approach in a future study and investigate explicitly the dependence of the number of close approaches per crossing time on the total number of stars involved in the interaction.
Intriguingly, the collision rate, as derived from the mean free path approximation (e.g. Leonard 1989), also scales as N 2 . In the limit of very large N , we would expect our relation for the collision probability to agree with what is predicted from the mean free path approximation. On the other hand, in the limit of small-N , it is not clear that the standard assumptions of the mean free path approximation should still hold. For the angular momentum and energy regime studied here, the collision probability appears consistent with that of the mean free path approximation, at least for N = 4, 5, 6. Further study is required to determine, e.g., how different energy and angular momentum regimes affect this relationship, what the minimum N is for which the standard mean free path approximations are still valid, etc.
The second most noticable similarity between all Runs is an exponential decline in the collision probability with increasing angular momentum that is steeper for larger N . This effect is small, however, relative to the total drop in the collision probability as N decreases within a given Run (see Figures 1, 2, and 3). One possible explanation that is consistent with our results is that the number of close approaches per crossing time may decrease with increasing angular momentum more steeply for larger N encounters.
Our results are applicable to a regime of low angular momentum and high absolute orbital energy. As the integration time increases with increasing angular momentum, we focussed our attention on the low-angular momentum regime in this paper to maximize the number of simulations performed for each Run, and thereby increase the statistical significance of our results.
We cannot probe lower angular momentum encounters without violating our assumptions. There are two reasons for this. First, the lower boundary for the long-term stability of triple systems corresponds to a ratio between the inner and outer orbital separations 7 (Mardling 2001). Therefore, we cannot lower the total angular momentum by decreasing the semi-major axis of the outer orbit of the triple without simultaneously reducing the semi-major axis of its inner orbit. This brings us to our second requirement, namely that q ≪ a0 (where a0 is the semi-major axis of the shortestperiod orbit initially involved in the interaction and q is the physical sizes of the objects involved in the encounter). We performed two additional Runs with a0 = 0.1 AU to test the limit of the assumption q ≪ a0. In both cases, the resulting chi-squared values found using the best-fits for our free parameters in Equation 2 were significantly larger than we found for our other Runs. We interpret this as being due to a breakdown of the requirement q ≪ a0. Our results suggest that the assumption q ≪ a0 holds provided a0 100q. This is only a rough estimate for the lower limit for the ratio q/a0, and more work will be needed to better constrain its precise value.
In future work, we intend to address the N -dependence of the collision probability for higher angular momentum encounters, as well as different mass-ratios and eccentricities. A non-circular orbit provides an additional free parameter that can be changed to affect the total angular momentum, but not the total energy. Therefore, the addition of a nonzero eccentricty should in principle allow us to include 1+2 encounters in our estimates for the collision probability using the same or similar normalization method (i.e. fixing the total angular momentum and energy) as adopted in this paper to compare between the different encounter types.

SUMMARY
In this paper, we perform a large suite of numerical scattering experiments to study the probability of a direct physical collision occurring during single-binary, binary-binary, single-triple, binary-triple, and triple-triple interactions. We quantify the dependence of the collision probability on the number of objects involved in the interaction N for fixed total energy and angular momentum.
Our results suggest that the collision probability increases approximately as N 2 , for N = 4, 5, 6. This result is consistent with the hypothesis that the average number of close approaches per crossing time also scales as N 2 , and we will fully investigate this possibility in a future study. Interestingly, this same N -dependence is predicted by the mean free path approximation in the limit of very large N . This similarity is rather intriguing, and further work investigating the connection between these two regimes in N will be highly valuable for our understanding of collisional dynamics in the realm of not-so-small-N .

APPENDIX A: SCHEMATIC REPRESENTATION OF LIOUVILLE'S THEOREM
In this appendix, we describe the construction of schematic diagrams that quantitatively describe the evolution of an interaction in energy-space.

A1 Motivation
Previous numerical scattering experiments for 1+2 interactions have demonstrated that the probability for a dynamical encounter to result in a direct collision is approximately proportional to the number of close approaches per crossing time multiplied by the number of crossing times the system survives through (e.g. Valtonen & Karttunen 2006). This is the basis for the analytic model described in Section 2.2.1 to predict the collision probability during resonant 1+2 encounters, and we have extended this formalism to also describe encounters with N > 3. We briefly discuss the difficulties in determining the precise number of close approaches during an encounter from an analysis of position-space in Section 4.
Here we present an improved method for determining the number of close approaches per crossing time by analyzing encounters in energy-space. There are three clear benefits to this technique, which we will expand upon below. First, defining the condition for a close approach in terms of energy rather than position maps each close approach to a nearly discrete change in energy, which is easily calculated. Second, the energy-space diagrams presented in this paper contain quantitative information pertaining to the relative energies of the objects during each stage of the encounter, whereas position-space diagrams do not (and are often very hard to interpret visually). Third, this method will potentially allow for a more general formalism to describe the collision probability during a dynamical encounter.
Below, we outline our method and provide a few example scenarios. We conclude by describing the future work that will be enabled by this energy-based approach.

A2 Liouville's Theorem
Liouville's Theorem is a key postulate in classical statistical mechanics. It states that the time evolution of the distribution function corresponding to an Hamiltonian system is constant along any trajectory in phase space (Liouville 1838). Said another way, the density of points representing particles in x,p phase space is conserved as the system evolves in time, where x and p denote the 3-D position and momentum vectors, respectively, of the particles. This is a remarkable and powerful result that can be applied to any dynamically-evolving system provided the forces acting on the particles are conservative and differentiable. This last condition excludes collisions due to the sudden introduction of additional forces, such as tides, shocks, etc.
As a system evolves dynamically, energy and angular momentum are exchanged between particles. This diffusion in energy-and angular momentum-space governs the trajectories of the particles through position-and velocity-space. Therefore, it is often the case that consideration of the evolution of a system in energy-and angular momentum-space, as opposed to position-and velocity-space, will more directly reveal the underlying physical processes that determine its outcome. As we will show in the subsequent sections, it follows that patterns in the dynamical evolution of an interaction are often most apparent in energy-and angular momentum-space.

A3 Diagrams for N = 4
We will begin by describing the rules for the construction of our schematic diagrams for the case N = 4. The total energy for a four-body encounter can be written: where mi is the mass of object i, vi is the speed of object i with respect to the centre of mass of the system, and rij is the distance separating objects i and j. We note that rij = |r i − r j |, where r i is the vector separating object i from the centre of mass of the system. We re-write Equation A1 in the form: Ei.
Now, we can form a polygon by setting each of the angles equal to: (Note that the angle 360 • here is a result of the more general formula 180 • × (N − 2), when N = 4.) As the system evolves, the total energy E remains conserved but the individual energy terms Ei will change. This is manifested visually via the transformation of our schematic diagrams with successive time-steps. The rules of our diagrams are such that the individual angles of the polygon change as the system evolves, while their sum remains constant.
Bound objects are connected by solid lines, whereas objects that are unbound are connected by dashed lines. Note that, if an object is unbound, then Equation A4 suggests that the angle corresponding to its vertex will be negative. We will come back to this below.
To illustrate our methodology, we show the evolution in energy-space of a 2+2 encounter in Figure A1. The corresponding position-space diagram is shown in Figure A2. The latter figure shows the evolution of a 2+2 encounter in position-space projected onto the x-y plane. In this interaction, two identical binaries that are composed of equal-mass (1 M⊙) stars with semi-major axes of 5 AU approach from the left and right (we denote the initial time by t = t0) to meet at the origin. In Figure A1, Stars 1 and 2 comprise one of the binaries, and Stars 3 and 4 comprise the other. Both of these bound pairs are connected with solid lines, respectively, while all remaining pairs of stars are initially unbound from each other, and therefore connected via dashed lines. Shortly after the binaries meet at the origin in Figure A2, one star (Star 2) is ejected from the system, and is expelled toward the lower left of the diagram (t = t1). Star 2 is now an unbound single object, and hence is connected to all other stars via dashed lines in the top left inset of Figure A1. The angle corresponding to its vertex is also now negative, since Star 2 has a positive total energy.
The three remaining stars then undergo a resonant interaction as they drift toward the upper right of the diagram in Figure A2 due to conservation of momentum (t = t2). At roughly (x, y) = (36, 61), another star (Star 4) is ejected toward the lower right. The left-over binary (composed of Stars 1 and 3) leaves the diagram at the upper right, with a semi-major axis (and hence orbital energy) that is smaller than those of the two initial binaries (t = t3). The final state of the system is depicted in energy-space in the lower right inset of Figure A1. Only Stars 1 and 3 are connected via a solid line, since they are the only two stars that remain bound. The angles corresponding to the vertices for Stars 2 and 4 are now negative, and the angles corresponding to the vertices for Stars 1 and 3 are both positive and obtuse (i.e. > 180 • ). Note that, although described in detail here, the specific evolution of the encounter is quite difficult to follow visually in position-space without the help of the previous text.
As a further step, we depict this encounter in Figure A3 using the formalism for making scattering diagrams presented in Hut & Bahcall (1983). In this plot, time increases from left to right, whereas vertical transitions denote the formation and/or destruction of temporary hierarchies in which two or more stars are in close proximity (i.e. strongly bound) to each other relative to all other stars. All four states of the system depicted in Figure A1 are shown in Figure A3. The initial (i.e. t = t0) and final (i.e. t = t3) states correspond to those shown at the far left and far right, respectively, of Figure A3. The state shown at t = t1 corresponds to a time shortly after Star 2 has been ejected from the system. Note that, within the triangle formed from the vertices of Stars 1, 3, and 4 at t = t1, all of the angles are comparable. This means that a temporary hierarchy has formed where Stars 1, 3, and 4 all have similar energies. This is represented in Figure A3 where the first convergence of all three lines (for Stars 1, 3, and 4) occurs simultaneously. This scenario occurs once more (i.e. at t = t2) before the interaction is Figure A1. Schematic diagram showing the evolution in energyspace of a 2+2 interaction between identical binaries. The panels have been arranged in chronological order, starting at the lower left inset and rotating clockwise (i.e. t 0 < t 1 < t 2 < t 3 ) . A detailed description of this interaction has been provided in the text. The state of the system depicted in the inset corresponding to t = t 2 shows one example of the many transformations of our diagram that occur between ejection events. This is due to the occurrence of several close approaches. eventually terminated by the escape of Star 4 from the system, leaving Stars 1 and 3 bound in a binary (which has a smaller orbital separation than the initial binaries going into the encounter).
To summarize, after Star 2 is ejected, a three-body system remains that evolves via the formation of a temporary binary that mediates a series of successive ejections of the remaining single star (which can also be exchanged into and out of the temporary binary). Eventually, the single star receives a positive total energy and escapes from the system.

A4 Future Work
One primary application for this method is to determine the number of close approaches per crossing time for an encounter involving N stars. We outline this method here. A "close approach" can be defined using Equation A4 to convert a minimum distance (i.e. q in Equation 1) into a minimum angle. 3 If we require that the distance of "close approach" is comparable to the physical sizes of the objects (which are assumed to be small relative to the characteristic size of the interaction region and the initial semi-major axes of any orbits going into the encounter), then these events 3 Note that in the point-particle limit a "collision" is considered to occur when two objects are within some minimum distance of each other. Therefore, our energy-space diagrams can be used to depict collisions. The encounter involves an essentially parabolic collision between two identical binaries composed of equal mass stars, initially approaching from the left and right. Distance is measured in units of the initial semi-major axes a of the binaries. When the encounter is over, two single stars and a binary emerge.
correspond to times of near instantaneous and significant energy-exchange between two individual stars. This is because the stars tend to be strongly accelerated or decelerated during very close approaches, so that significant energy and/or angular momentum can be exchanged. It follows that a "close approach" will appear as a near discrete transformation of our energy-space diagrams, whereas the evolution appears more continuous and chaotic in position-space. This is demonstrated in Figure 3 of Hut & Bahcall (1983) and also in Figure A3. Therefore, individual close approaches will correspond to high values for dαi/dt, where αi is the angle corresponding to star i in an energy-space diagram. Consequently, one can define a "close approach" as occurring when dαi/dt > η, where η is a parameter that can be varied to specify the precise definition of a "close approach". We intend to study the time-evolution of the quantity dαi/dt in future work. As is evident from Figure A2, the time evolution of the system is very difficult to observe in position-space, and counting the number of close approaches is essentially impossible from this type of a diagram. In contrast, we can easily count ten (to first order) close approaches in Figure A3. To at least first order, each time a close approach occurs, one of the the lines in a scattering diagram makes a vertical transition. Figure A3. Scattering diagram for the 2+2 interaction shown in position-space in Figure A2 and in energy-space in Figure A1. The time evolution of the interaction proceeds from left to right.
This energy-space formalism is easily extended to depict the evolution of encounters involving more than 4 objects. (In principle, there is no limit to the number of objects that can be described with this method.) Furthermore, we can, in principal, use this formalism to re-write the functional form for the collision probability shown in Equation 2 in terms of only the number of objects N , the total angular momentum, and energy. This is because our method re-defines the condition for a "close approach" in terms of energy, instead of distance q. The term q/a0 in Equation 2 can thus be replaced with a new term that depends directly on energy. This will further generalize our approach, and should be an improvement over the normalization used in this paper to compare between different encounter types within a given Run.
In addition, the same method we have presented to depict the evolution of the system in energy-space can be modified to depict its evolution in angular momentum-space. The prescription for this is analogous to the formalism presented in the previous sections for energy-space. In general, our technique can be used to depict the time evolution of any conserved quantity.
We intend to use our diagrams in future studies to simulate the dynamical evolution of small-N gravitational interactions in both energy-and angular momentum-space. This will be done in (accelerated) real time, so that we may simultaneously view the evolution of a system in position-, energy-, and angular momentum-space. This new tool will be a useful addition to studies with the over-arching goal of identifying trends in the evolution of dynamical interactions.

ACKNOWLEDGMENTS
A very big thank you to John Fregeau for trouble-shooting assistance when adapting the FEWBODY code, and for a critical read of our manuscript. Thank you to Alison Sills and Steinn Sigurdsson for valuable feedback and guidance. We would like to thank the following people for useful discussions: Hagai Perets, Maureen van den Berg, Evert Glebbeek, and Bob Mathieu. NL is funded by the European Space Agency (ESA) Postdoctoral Fellowship. AMG is funded by the Lindheimer Fellowship at Northwestern University.