An Analytical, Statistical Approximate Solution for Dissipative and non-Dissipative Binary-Single Stellar Encounters

We present a statistical approximate solution of the bound, non-hierarchical three-body problem, and extend it to a general analysis of encounters between hard binary systems and single stars. Any such encounter terminates when one of the three stars is ejected to infinity, leaving behind a remnant binary; the problem of binary-single star-scattering consists of finding the probability distribution of the orbital parameters of the remnant binary, as a function of the total energy and the total angular momentum. Here, we model the encounter as a series of close, non-hierarchical, triple approaches, interspersed with hierarchical phases, in which the system consists of an inner binary and a star that orbits it -- this turns the evolution of the entire encounter to a random walk between consecutive hierarchical phases. We use the solution of the bound, non-hierarchical three-body problem to find the walker's transition probabilities, which we generalise to situations in which tidal interactions are important. Besides tides, any dissipative process may be incorporated into the random walk model, as it is completely general. Our approximate solution can reproduce the results of the extensive body of past numerical simulations, and can account for different environments and different dissipative effects. Therefore, this model can effectively replace the need for direct few-body integrations for the study of binary-single encounters in any environment. Furthermore, it allows for a simply inclusion of dissipative forces typically not accounted for in full N-body integration schemes.


I. INTRODUCTION
The three-body problem -that of the interaction of three gravitating objects -is one of the oldest physical and astrophysical problems, and has been studied for centuries; solving it requires an understanding of the outcomes of encounters between binary systems and single stars. Such encounters play a key rôle in the evolution of binary systems, both in dense environments such as globular or nuclear clusters [1,2], and in the field, too [3,4]. These encounters can result in the formation of compact binaries, and consequently engender a plethora of physical phenomena, ranging from mergers of stars and compact objects, to the production of electromagnetic and gravitational-wave transients, such as type Ia supernovae, short gamma-ray bursts, et cetera. Currently, full numerical integrations of few-body systems are required in order to characterise the results of such interactions. However, the large number of such few-body simulations needed for one cluster comes with a significant computational cost; any change in the basic aspects of the problem (e.g. different masses, binary configuration etc.) necessitates a new set of simulations; and likewise different environments (e.g. external, outer potentials in clusters) require different sets of simulations (and in many cases are not even consistently accounted for). Moreover, realistic physical processes that may affect the outcomes, such as dissipative forces (e.g. tidal interactions, or gravitational-wave dissipation) acting on the interacting stars are rarely incorporated into simulations of the evolution of binary systems through binary-single encounters. Finally, such simulations, by their very nature, provide neither any explanation for their outcomes, * ginat@campus.technion.ac.il; hperets@physics.technion.ac.il nor any direct insight into them, and are used, to some extent, as black-box ingredients in simulations of largescale systems. Despite much progress in the analytical understanding of the problem, the full solution for the outcomes of non-hierarchical three-body systems, and in particular the general understanding of binary-single encounters and their long-term end-states under realistic conditions remains unsolved. Furthermore, accounting analytically for dissipative aspects of this problem was not done previously, to the best of our knowledge.
Here we provide a full analytical, statistical model that solves the non-hierarchical three-body problem in the limit of encounters with hard binaries. We show how our approach can account for the long-term evolution of a non-hierarchical triple system and characterise it as it goes through consecutive close binary-single encounters until the system is destroyed. Our solution provides the detailed cross-section for the final outcomes, including both cases of ejection of one of the components or the collision/merger of two of the components, as well as, most importantly, the characterisation of the final remnant binary properties. Our approach can also account for different spatial cut-offs due to external perturbations. Last, but not least, it can generally account for dissipative forces -we exemplify its application for dissipative tidal forces during binary-single encounters. Thus, the use of our model, which is based on a randomwalk approach, can potentially replace the need for direct integration of binary-single encounters and provide a general, robust tool for the understanding the outcomes of binary-single encounters.
Binary-single encounters can generally be divided between the cases of encounters with hard or soft binaries, i.e. according to the ratio between the binary binding energy and the kinetic energy of the incoming third star. When the binary star is wide -when its binding energy is considerably smaller than the typical kinetic energy of a star in the host cluster, k B T (where T is proportional to the squared velocity dispersion of the cluster), the encounter is well-described by a composition of two two-body interactions [1,5,6]. In the case where the binary is 'hard', i.e. when its binding energy is much larger than k B T , an analytical description of the encounter is more difficult. If the encounter stays hierarchical, it may be treated by means of perturbation theory [5], but if it doesn't, it becomes a so-called 'resonant encounter', where there are phases when all three bodies are close to each other, and the evolution is inherently chaotic. One could still find some analytical insight by investigating the phase-space evolution of a set of various initial conditions during such a close encounter. For example, it has been shown at the early 20-th century that if the system starts as a binary and an unbound third star, then, eventually, it will end up with at least one unbound star (apart from a set of measure zero); see Arnold et al. [7, § §2.4, 2.6] for a review of this theorem and related ones.
Here we wish to present an analytical, statistical model for the evolution of such encounters with hard binaries. We will endeavour to derive a probability distribution function for the possible end-states. The initial binary is hard, so the total energy of the system, E, must be negative, whence an end-state with three stars unbound is forbidden -the set of all possible end-states is therefore the set of all possible final binary configurations, multiplied by the set of configurations of the ejected star (which, of course, may be a different star from the initial unbound one -see §III below). We derive a distribution function, f bin (E bin , S|E, J), for the resultant binary to have energy E bin and angular momentum S, given that the total energy is E and that the total angular momentum is J.
Let m 1 , m 2 , m 3 denote the masses of the three stars, which we take to be of similar magnitude, and let M = m 1 +m 2 +m 3 be the total mass. Further let the subscript 'bin' denote any quantity related to the binary, such as its total mass m bin , the reduced mass of its two components µ bin , et cetera. Likewise, let a subscript s denote quantities pertaining to the ejected star, like µ s = m s m bin /M . We denote the total energy of the triple by E or sometimes by E 0 , and its total angular momentum by J. More explicitly then, the problem we investigate is the following one: we consider a single star with velocity v 0 at infinity, which is scattered on a hard binary, such that the total energy is E and the total angular momentum is J. As the binary is hard, E ≈ −Gm a m b /(2a 0 ), where a 0 is the initial semi-major axis of the binary, consisting, initially, of stars a and b. For a given impact parameter b of the single star, relative to the binary centre-of-mass (as well as all the initial anomalies of the binary), it is possible, in principle, to predict the outcome of the encounter, i.e., the final values of E bin , S and E s . However, due to the chaotic nature of the close three-body interaction, doing so is impracticable. One could study the problem either by performing numerical simulations (see, e.g., refs. [8][9][10][11][12][13][14][15][16][17][18][19][20][21]), or, statistically, using various analytical approximations. Such analytical treatments presuppose that the results of such numerical scattering experiments may be treated as random variables, drawn from some distribution; in effect, it reduces the problem to finding that distribution. Numerical simulations have shown that the encounter proceeds as a sequence of many close triple approaches, after each of which a single star is ejected [11,12,19]: if the star is still bound to the other binary, it returns eventually, whereupon a new close triple approach ensues, and so on, until the single star is ejected to infinity with positive energy.
There are two dominant analytical approaches in the literature (to our knowledge): one, initiated by Monaghan in [22], relies on chaotic mixing during the close triple approach to argue that the distribution is ergodic in the relevant part of the system's phase space. This approach was refined later by various authors, e.g., refs. [6,[23][24][25], to account more accurately for angular momentum conservation, until Stone and Leigh [26] succeeded in doing so fully, for the unbound case.
The other approach, which dates back to Heggie [5], was to draw upon the principle of detailed balance to deduce what the outcome distribution, f bin (E bin , S|E, J), must be in order for the number of bound triples formed by such encounters to be fixed, in a cluster in thermal equilibrium (see, e.g., refs. [2,5,15]). Various authors also attempted to account, in simulations, not only for the Newtonian interaction of three point-particles, but rather include other physical effects, such as collisions [17,27], and more recently also gravitational-wave emission and tidal dissipation [19,28]. The reader is referred to chapters 7-8 of ref. [6] for a review of some of the work on binary-single scattering.
The contribution of this work to the understanding of binary-single encounters is threefold: first, we present the first (to our knowledge) closed-form statistical approximate solution of the bound, non-hierarchical three-body problem that takes both energy and angular momentum conservation into account, which complements the solution of the unbound case of ref. [26]; using the results presented here, one can compute the outcome distribution of each intermediate close approach, not only of the final one. 1 Secondly, we show how to derive the solution both using detailed balance arguments and using ergodic arguments, thereby unifying the two approaches. That they give the same result is hardly surprising, since the principle of detailed balance is essentially a statement about phase-space volumes. Thirdly, and most significantly, we elevate our solution of the bound case to a random-walk model of the entire encounter, in which the binary actions and the three-body conserved quantities perform a random walk, whose transition probabilities are related to f bin (E bin , S|E, J). This random walk model is extremely general, and can incorporate many physical processes beyond the Newtonian gravitational three-body interaction and in addition to it. Here we provide as an example the important case of dissipation due to tidal forces, but our approach can be generalized to any other type of additional perturbations and physical processes, which will be done in future work.
The structure of this paper is as follows: we start with an explicit calculation of f bin (E bin , S|E, J) using ergodic arguments, in §II, while §IV presents a computation of the same function using the principle of detailed balance. In §III and §V we discuss two marginal distributions: the probability function of the ejected star's mass, and the marginal energy distribution. Our paper culminates in §VI, where we move on to present the random walk model in its full generality. In §VII we compare the results of § §III and V to past numerical simulations, and in §VIII we apply the random walk model to tidal dissipation, while also comparing it with relevant simulations. All of our results mesh well with the simulations.

II. CROSS-SECTION VIA PHASE-SPACE INTEGRATION
All the system's phase-space is divided into three parts: one where the system is hierarchical, but the outer body is unbound; one where the system is still hierarchical but the outer body is bound; and a chaotic region in which all three bodies interact closely. Let us denote them by A, B and C, respectively. The latter is taken to be the collection of points where the maximum relative separation between any one body and the centre-of-mass of the other two is no more than some value R, which is defined below. Then, if chaos is sufficiently strong to lead to phase-space mixing inside C, the cross-section of the outcome of a close approach is simply an integral over C, as assumed in previous studies [6,[22][23][24]26] where r j and p j are the position vector and the momentum of the j-th particle, H is the Hamiltonian, and E, P CM and J are the total energy, momentum and angular momentum, respectively. This equation is equivalent to the statement that the probability of the system exiting C through a point w is actually independent of w; we give a heuristic argument for this statement in appendix B.
One can calculate this integral by transforming to the co-ordinate system of a binary and a lone star orbiting its centre-of-mass, which amounts to demanding that the system becomes hierarchical when it leaves C; this requirement is what defines R below. Working in the centre-of-mass frame: where E bin = −Gm 1 m 2 /2a bin is the binary energy, and where the sign is negative/positive when the system goes into A/B respectively. The spin of the binary is S and the angular momentum of the third body about the binary is L. Below the integral over the remaining binary phasespace is shortened to bin . As there are three stars, with possibly different masses, there are three distinct outcomes, depending on which of the three bodies ends up being ejected. As the original binary was hard, the final binary must also be hard. Formally, we weigh the binding energies between each of the three pairs, and the one whose binding energy is considerably lower (i.e. more negative) is the remaining binary. This implies that the total cross-section splits into a sum of three cross-sections where σ i is the cross-section for a break-up with star i ejected as the lone star in the hierarchical system. This also gives an infra-red cut-off for E bin naturally, forcing it to be less than E/3, for if it were more, one of the other pairs would actually be harder. Below we calculate such a σ i ; because all three are symmetric, we ignore this complication for now, and return to it in §III.

A. C in Angle-Action Variables
Performing this integral in angle-action variables is the best way to proceed [26], but first it is essential to determine the integration region explicitly: in the initial encounter, the third body must come close enough to the binary, such that its pericentre distance is of the same order as R. This implies that a s (1 − e s ) R for an elliptic orbit, or a s (e s − 1) R for a hyperbolic orbit.
Equation (5) needs to be supplemented by another condition, which says that while the system becomes hierarchical, the lone star's apoapsis is sufficiently larger than R. This condition is automatically verified for a hyperbolic orbit (where the lone star escapes and the preceding close encounter is the final one), but for an elliptic orbit (where the lone star eventually returns), one must have for η ≈ 5 (other values are fine, too).
Let us discuss how R is related to the binary parameters: R is defined as the critical distance between the would-be ejected star to the centre-of-mass of the other two stars, where the problem becomes hierarchical. Conversely, for a hierarchical triple, one might write the Hamiltonian (in the centre-of-mass frame) as where m a and m a are the inner binary's masses, r bin is the distance between them, Φ is the angle between r bin and r s , and P n is the n-th Legendre polynomial. We refer the reader to refs. [6,29] for details and for more references on the hierarchical three-body problem. Now, the triple ceases to be hierarchical if the multipole series becomes as large as the leading order term in the Hamiltonian, namely, as big as the energy E = E bin + E s . The leading term in this series is the quadrupole (for non-extreme mass ratios). Approximating r 2 bin P 2 (cos Φ) ≈ a 2 bin , we find where β is a constant of order unity. If the masses are un-equal, this formula might lead to an over-estimate, since then it is possible that an exchange of a light star with a heavy star would yield R a bin , in which case the problem is still quite visibly hierarchical. To account for that, we define R as where, again, β ≥ 1 is of order unity; the cross-section σ depends on β only weakly (cf. ref. [26] for the unbound case). In future work we will explore the consequences of the existence of this threshold for the stability of hierarchical triples. If the reader is concerned about the crude approximation of the Legendre polynomial, we offer a more refined one, and test both in appendix D, but we recommend that it be read only after the next section. Both conditions (5) and (6) may be translated into conditions on the angular momentum L (they are, evidently, independent of its direction). Using energy conservation one may express E s in terms of E bin , whence conditions (5) and (6) may be written as L ≤ A(E bin ). Let us start with the apoapsis condition: First, if ηR < a s , then this condition is fulfilled automatically. If not, it simplifies to As the left-hand-side is non-negative, this implies that ηR ≤ 2a s . Therefore a s can either be more than ηR, or more than ηR/2, so in both cases more than ηR/2. Thus, one has a restriction on the energy-difference |E − E bin |, viz.
The elliptic periapsis condition is satisfied trivially if a s < R, but if η > 2 this cannot be the case. Otherwise, one has Which of conditions (11) and (13) is more stringent depends on the masses, and on a s . For the unbound, hyperbolic case, the periapsis condition is equivalent to note that as a consequence of the plus sign inside the brackets here, there is no restriction on a s in this case, and it may be as small as one pleases. The function A(E bin ) is therefore defined as The approximation made here, which includes a sepa-ration of phase space into the three regions A, B and C, and the different treatment of each, inadvertently induces some uncertainty. We attempt to gauge it in appendix C, but we urge the reader to peruse §IV, §III and §V before turning to this appendix.

B. Cross-Section Calculation
Now one may turn to performing the integration in equation (2). The goal is to find the final distribution of binary spin and energy, of the remaining binary, after the system stops being chaotic. Thus, one has to integrate over all of the lone star's phase-space, as well as the angle variables of the binary. First, though, note that there are two possibilities for each encounter: either E s < 0, or E s > 0 after the encounter (E s = 0 has zero measure). So σ = σ bd + σ ubd , where each cross-section pertains to each possible sign of E s . 3 σ ubd has been calculated by Stone and Leigh [26]. They obtained The other cross-section, σ bd , may be calculated using Delaunay variables for both the emergent binary and the binary formed by the star and the inner binary. We do so below; this is meaningful given that the whole encounter proceeds as a series of consecutive close approaches, each one having a cross-section σ bd . While this is supported by numerical simulations, as mentioned in the introduction, we also give a heuristic argument for it in appendix B. The single close approaches are combined below in §VI. The Delaunay variables are denoted by (J a , J b , J c , θ a , θ b , θ c ) and are defined in, e.g., ref. [1].
Using a superscript or a subscript s to denote variables pertaining to the outer binary, we have This implies that the angular-momentum-conserving delta-function is This equation in turn implies that the angular-momentum integral is independent of J s c , modulo the integration domain boundaries, which are (Please note that for the unbound case, α is defined simply as A(E bin ).) Theẑ-axis integral gives As in ref. [26] we perform a change of variables The Jacobian of this transformation is Integrating over z 1 , z 2 yields The integral dθ s b gives 2π, while the integral dθ s c -over the mean anomaly -gives a multiplicative factor of θ max , which is the maximum mean anomaly the star may have and still stay in C. Condition (5) implies that θ s c = 0 is in C, while condition (6) implies that θ s c = π is no longer in C. Thus, θ max restricts |r s | to |r s | ≤ R, such that the integration is carried out in C. The last lone-star integration, over J s c may now be performed, to remove the last delta function, and give where One may perform the integration over the binary angles trivially, to give an additional factor of (2π) 3 . This yields a cross-section Taking theẑ-axis in this integral to be along J implies that the integration domain is where E min is some minimum cut-off on the binary energy. This form is the same as that of σ ubd of ref. [26], which implies that where now the integration domain is and θ max is defined by equation (28) for the bound case, and by equation (16) for the unbound case. Equation (31), together with equations (32), specify the probability that the resultant binary has Delaunay actions J a , J b , J c , given conserved quantities E, J -their probability density function is simply the integrand in equation (31). It might be more useful to express equation (31) in term of E bin , rather than J c ; the Jacobian for this transformation is simply ∝ E −3/2 bin , which gives a distribution function for the outcome of one binary-single close approach,

III. DIFFERENT MASSES
In fact, equation (33) is proportional to the probability that a star s ∈ {1, 2, 3} escapes, leaving a binary with energy E bin and spin S. What remains is the coefficient, which depends on the masses, as in equation (31). Therefore, the probability density that star s is ejected after a close interaction, and that the remaining binary has energy E bin and spin S is where N (s) is a normalisation constant, that depends on s through the integration over (32). To obtain this constant, one should change the integrands in σ to ones that are independent of the masses, i.e. from (J a , J b , J c ) in equation (31) to semi-major axes, eccentricities and inclinations. The lone star's angular momentum and energy simply contribute a factor of µ −5/2 s , and the measure contributes an additional factor of µ 3 bin m 3/2 bin . Thus, up to dimensionless quantities, where, now, the proportionality constant is independent of the identity of the ejected star (but may still depend on the total mass or on conserved quantities). One immediate prediction of equation (35) is that, if one of the masses is much smaller than the other two, then the probability that each close interaction ends with the lighter one being shot out, rather than one of the others, dominates. By dimensional analysis, the probability that mass m escapes is therefore We emphasise that equation (36) is an approximation, and for accurate results one should integrate equation (34) over remnant binary energies or angular momenta or the marginal energy distribution, equation (48) below, over the allowed energies. 4 One can also compute the exchange cross-section, given by as the total cross-section is of course where a 0 is the initial binary semi-major axis, and v 0 is the perturber's initial velocity (see, e.g. Heggie and Hut [2]).
Let us try to obtain equation (31) by an easier means. One might, for instance, assume that the triple is part of a globular cluster which is in thermal equilibrium, and which contains binaries, single stars, and triples. Of course, the desired cross-section σ does not depend on whether there exists such a cluster, and on whether this hypothetical cluster is indeed in thermal equilibrium. Making those auxiliary assumptions would simplify the calculation, because then we could use the principle of detailed balance (cf. refs. [2,5,15]).
If the cluster is in thermal equilibrium, then the numbers of binaries with energy E bin and stars with energy E s must be constant. This number density is proportional (in the canonical ensemble, neglecting collisions and stellar evolution) to exp(−β * H ), where β * = 1/(k B T ). From this, one may deduce that the rates at which encounters between stars and binaries transfer these systems from one energy to another, and those of the reverse process, have to be equal. This is the principle of detailed balance.
Let Γ(E bin , S, p → E 0 , J)dE 0 dJdE bin dJ a dJ b be the differential rate at which binaries and single stars combine to form bound triples with energy E 0 and total angular momentum J, and likewise let Γ(E 0 , J → E bin , S)dE 0 dJdE bin dJ a dJ b be the rate of disintegration of such triples. Let n triple (E 0 , J) denote the number density of such triples, and let n bin (E bin , S) and n(p 3 ) denote the phase-space densities of binaries and single stars, respectively. Detailed balance amounts to the requirement that the number of bound triples remain constants, i.e. that the rate of the forward and backward reactions, weighed by the relevant densities, cancel each other out. Symbolically, in barycentric co-ordinates where p is the momentum of the relative motion between the third star and the binary centre-of-mass. (See, e.g., ref. [15] for a derivation of a similar expression for the reaction rate.) The advantage of equation (39) is that the rate of formation of bound triples may be simpler to compute; one could approximate it as where D(v) is a disc in thex-ŷ plane at infinity, with radius a bin 1 + GM µsv 2 a bin (see, e.g. ref. [15]). 5 The direction of theẑ-axis is chosen so that it points along v, (assumed not to be away from the binary, for then there wouldn't be any encounter). The integral over b, including the angular momentum delta-function gives The remaining delta-function should not really be there -it is just a mathematical artefact, which came from the way we defined the axes, so it is safe to omit it. Alternatively, one could justify its omission by integrating over all possible orientations of the axes: this delta function then picks out the orientation where theẑ axis is perpendicular to J − S; such an integration is in turn justified by the fact that a choice of axis is meaningful only for the right-hand side of equation (39), and not for the left-hand side. The factor of |J − S| in the denominator comes from expressing the angular momentum delta function in spherical co-ordinates: this turns the 3-dimensional Dirac delta function δ (J − S − L) (recall that µ s bvφ = L) into a product of three 1D delta functions, one for the magnitudes of J − S and L, one for one angle, and one for the other angle, divided by the appropriate Jacobian |L| 2 sin θ J−S . Since the other delta function sets θ = π/2, we are left with |L| 2 in the denominator, and in the numerator: two angular delta 5 We remark that one cannot assume that µsv 2 is much smaller than E bin ; if we did, then thermal equilibrium would be precluded, as interactions between hard binaries and single stars are known to be a source of heat for globular clusters (see, e.g, ref. [2]).
functions multiplied by δ (|J − S| − |L|). We now perform the integral d 2 b in polar co-ordinates b and ϕ, by writing |L| = µ s bv, and changing variables from |L| to b. As d 2 b = bdbdϕ, one b cancels one power of |L| from the denominator, and the integral db removes the absolute value delta function, replacing the other |L| with |J − S|. The integration dϕ removes one of the angular delta functions, and the second one is removed as explained above. The density n(p 3 ) is a Maxwell-Boltzmann distribution, which is proportional to ρm −3/2 3 exp(−β * E s ). After performing the integral over p = µ s v, keeping track of the energy-conserving delta-function, one has . Heggie [5] gives a formula for the density of binaries with energy E bin , eccentricity e and inclination i [5, equation 2.12], from which the phase-space density n bin (E bin , S) is de- How is this rate related to f bin (E bin , S|E 0 , J)? By definition, it is the number of disintegrations per unit time, i.e., it is the number of triples in which star s is between θ s c and θ s c + Ω s c dt, divided by dt, where θ s c ∈ [0, θ max ]. That is, otherwise.
(44) up to a function symmetric in all the particle masses, as in equation (33) I is the integral one must evaluate. The calculation is performed in appendix A, and the outcome is that I is well-approximated by a power-law, proportional to E −1/2 bin for low J, but to E −1 bin for large values of J. A consequence of §II A (inequality (12)) is that for E bin > E 0 -i.e. in the bound case -the final binary energy is constrained to lie close to the total energy. In this neighbourhood, θ max ∼ |E − E bin | 3/2 , which cancels the existing |E − E bin | −3/2 . Outside this region, θ max is approximately constant; thus, one may remove its eccentricity dependence by approximating e s ≈ 1, writing This implies that the marginal energy distribution is

VI. A RANDOM-WALK DESCRIPTION
We have now reached the point where we may introduce a random-walk description of the evolution between consecutive close triple approaches. Suppose that during each close approach, the constants of motion (E, J) might change by some amount, according to some probability distribution, which depends on the state of the system at the beginning of the close approach (i.e. at the end of the previous one), because of some additional astrophysical process. Denote this distribution by , where E j , etc. denote quantities at the end of the j-th close approach. A definition of such an additional astrophysical process one would like to incorporate in one's study of binary-single encounters amounts, therefore, to providing f c . Other- The total energy E, the total angular momentum J, the spin S and the binary energy E bin thus perform a random walk, where the probabilities for the j-th value are dictated by the values of these quantities at the j − 1th step -this random walk has a one-step memory. This process may describe, for example, a tidal interaction between two stars during the encounter (see §VIII below).
The beauty of this description is that now one can use it to find the ultimate binary parameter distribution P (E bin , S), when the single star leaves, never to return. Let x = (E bin , E, S, J), and let h(x|x ) denote the probability of the walker (i.e. the binary + single) moving from x to x at one step -that is, the probability that it started the close approach at x and left it at x. Explicitly (49) The ultimate reason we spent so much effort above computing f bin for the bound three-body problem is precisely so that we would know what h(x|x ) looks like. The mixing hypothesis ensures that the functional form of the way h depends on the E bin , S components of x is only through f bin as calculated in §II. In particular, it also allows us to account spatial cut-offs; e.g. in a dense cluster environment an ejected, but still bound third star, which would have otherwise eventually fallen back into C, would now be met with an external perturbation by other stars, if its separation became comparable to the distance between stars in the cluster. In other words, the environment could potentially dictate an effective binding energy limit which would be different from the clean case of an isolated interacting triple. Our model can easily incorporate this aspect.
Let us also introduce the following linear, integral operators, acting on a function ϕ(x): The first describes an encounter that ends with the third star bound, and the second -the final encounter. Suppose we start with initial probability where E 0 is the initial total energy and J 0 is the initial total angular momentum. Now, This means that the only dependence on E bin , S is outside the integral, inside an f bin -just as in equation (52).
In-so-far-as the binary actions are concerned, the action of W lim does not alter the functional form of the probability distribution.
Now suppose that f c may be expanded as a sum of changes in energy and angular momentum, relative to E 0 , J 0 , whose probabilities depend on the previous round: (this is the law of total probability in disguise) if this is the case, then where we have defined (57) Equation (56) implies that one encounter, if the initial probability distribution was some constant times f bin and a total-energy-total-angular-momentum deltafunction, turns this form into a sum of terms of similar structure. This fact helps us to find the final distribution of binary energies and spins, when tidal interactions are taken into account, as we do below.
Using the particularly special form of the action of W lim on p i (the action of W unlim is quite similar), we may determine the final distribution in terms of the initial total energy and angular momentum. This is done in a perturbative manner, assuming that the probability of a non-zero change in these quantities is small. 6 Let P n (x|x 0 ) be the probability that the full close interaction 6 Here we sum up all orders, so that this assumption is innocuous in-so-far-as the general model described in this section is concerned. In §VIII Below we truncate the series at linear order.
ends after exactly n steps, at x, given that it started out initially at x 0 . By the properties of random walks (see, e.g. Hughes [30]), P n is therefore W unlim acting once after n − 1 actions of W lim , on p i . The final probability is Equation (60) is the ultimate distribution of binary parameters, after a complete binary-single encounter, incorporating the physical process described by f c , in addition to the classical three-body dynamics. We apply this formalism to include the effects of tides and collisions in §VIII. Equation (60) implies that if the conserved quantities don't change, then P (x) is just f bin , and the entire process is rendered memory-less. We now move on to compare the theoretical predictions made in this paper with results of numerical simulations.

VII. COMPARISON WITH SIMULATIONS -MARGINAL DISTRIBUTIONS
Let us start by comparing some marginal distributions of f bin to numerical simulations, before moving on to test the full random-walk model in §VIII. All numerical integrations in this paper were done using MATLAB's integral functions.
We start with testing the ejected mass probability, which is given by equation (36), to two simulations, by Saslaw et al. [8] and by Hills [14] in figure 1. Please bear in mind that when the perturbing star's mass is much larger than the initial binary members' masses, a 0 in equation (37) needs to be modified by another multiplicative factor of [M/(m a + m b )] 1/3 , due to an increased effective total cross-section -this is what we show on the right panel of figure 1.
We also test the predictions of our model by comparing them to the simulations results of ref. [31], who calculated exchange cross-sections for a wide range of masses. In figure 2, we plot the predictions obtained by integrating equation (48) over the allowed range, as well as those of the approximate equation (36), and those of appendix A of ref. [25], for the following situation: the initial binary consists of masses m 1 = 1 M , and m 2 , and the in-coming star has mass m 3 = m 1 ; we plot the branching ration, defined by BR = σ exch (1) σ exch (2) , i.e. by the ratio of the exchange cross-section for ejecting star 1, and that of ejecting star 2. 7 Ref. [31] provides an analytical fit, which is also plotted, as well as data from ref. [17], for comparison. This fit is based on the entirety of the numerical simulations in ref. [31], which cover an extensive range of mass ratios. The data from ref. [31] in this figure are for resonant cross-sections, while the fit is, to our understanding, for the total one, which is dominated by the resonant cross-section everywhere, especially for large mass-ratios.
As one is concerned with exchange cross-sections, they decay to zero at sufficiently large impact parameters. The impact parameter serves only to determine the total angular momentum, and as the exchange cross-section tends to zero as J → ∞, one can integrate equivalently over J. There is a natural cut-off J * , which is the maximum angular momentum for which I does not vanish (for any E bin ), minimised over all ejected masses. Above J * , there is no configuration in which the triple could have been in a non-hierarchical phase before separating (cf. §VIII C below). The reader should bear in mind that if m 2 is considerably larger than m 1 , then BR 1, since there is a much-higher probability to eject the light particle, and likewise, for m 2 m 1 , m 3 , BR should decay to zero. All theoretical predictions plotted in figure 2 satisfy these limits, but only the exact prediction of equation (48) meshes well with the data and with the semi-analytical fit.
Next, we move on to show the semi-major axis distribution. Ref. [26] already obtained a good agreement between numerical simulations and the unbound crosssection, which has a similar form to equation (31). To check whether the bound one is also correct, we compare f bin to the numerical results of Sigurdsson and Phinney [17], who imposed an energy cut-off on the ejected star. There, as mentioned above, even some encounters with the lone star ejected with negative energy were deemed to be concluded, if its semi-major axis was large enough. This was meant to mimic the environmental effect of the globular cluster, where the triple resides. Clearly, once a single star is sufficiently far from the binary, it feels the cluster's potential more strongly and ceases to be bound to the binary, even if its energy is negative. This environmental cut-off implies that one has to use both σ bd and σ ubd to match the numerical results of ref. [17].
We do so by modifying the marginal energy distribution in equation (48) to account for the external cut-off criterion. Explicitly, Sigurdsson and Phinney [17] took a cut-off of a s (1 + e s ) = 960a 0 , where a 0 is the initial semi-major axis of the binary. When the apoapsis was larger than this value, they considered the third body to be unbound from the binary. This may be incorporated into f bin simply by modifying the apoapsis criterion in equation (5) accordingly. The result is compared with their simulation results in figure 3.
As one can tell from figures 1, 2 and 3, the function f bin calculated above agrees well with simulation results. We now proceed to test the random walk model in more detail in the next section.

VIII. COMPARISON WITH SIMULATIONS -TIDES
To test the random-walk model described in §VI, we assume that dissipation is caused by tides, and we compare with the extensive 3-body simulations conducted by Samsing et al. [28], which include both tidal forces and relativistic corrections. They investigated many cases, and we choose the equal mass case m = 1.2M ; the initial binary is comprised of a white dwarf and a compact object (i.e. point particle), and the incoming lone star is also a point particle. The white dwarf's radius is r * = 0.006R . The initial orbital separation is a 0 , and the initial speed of the third body is v 0 = 10 km s −1 . Its origin is sampled uniformly from a disc D(v 0 ) whose Right: The exchange cross-section, in units of the geometric cross-section, πa 2 0 , for binary masses equal to 1M , and third star mass m, which arrives with initial velocity v = 0.001v orb , compared with data from ref. [14].  [25] (green, dashed), compared with the semi-analytical fit of Heggie et al. [31] (blue, dotted with crosses) and numerical simulation data from Heggie et al. [31] (black circles) and Sigurdsson and Phinney [17] (purple asterisks). The initial binary masses are m1 = 1 M , and m2, and the in-coming star's mass is m3 = m1. Its velocity is v0 = 0.1vc (where vc is defined in the caption of figure 3) to ensure that the binary is hard. The y-axis shows the branching ratio of the cross-section for ejecting m1 relative to the cross-section for ejecting m2. There is excellent agreement with equation (48). [17, figure 2b]. The initial semi-major axis is a0 = 0.1 AU and the initial eccentricity is e0 = 0; the initial binary masses are m1 = 1.4 M , m2 = 0.56 M , and the incoming third star's mass is m3 = 1.4 M . Its initial velocity v0 is uniformly distributed between 0.05vc and 0.15vc, where vc = GM µ bin /(m3a0). The impact parameter is uniformly distributed in a disc D(v0) at infinity, whose radius is bmax(v0) = a0 (4vc/v0 + 0.6(1 + e0)). All these parameters were chosen to match those of ref. [17]. We use β = 1.5 but the results are insensitive to β. The straight line and the data marked by blue bars pertains to the exchange (1, 2)+(3) → (1)+(2, 3), while the dashed-dotted line and the data marked by purple bars are for the process (1, 2) + (3) → (1, 3) + (2), where the lighter star is ejected. Error-bars correspond to 3σ statistical (Poisson) errors, arising from a total integration number of 4000 [17, figure 9]. As expected by Sigurdsson and Phinney [17], the cut-off at large values of a is faster than exponential, occurring almost instantaneously. Data for fly-bys is not shown here, as in Sigurdsson and Phinney [17] it is dominated by adiabatic fly-bys, and the effect of resonant scattering is hard to disentangle from it. a a Observe that as the masses are different, the factor of m bin in equation (31) must be incorporated as well. Thus the normalisation of equation (48) is the sum, over all three possible final states, of the integrals of f bin (E bin |E, J) over E bin . Also note that the graphs shown here are normalised such that the cross-sections (integrated over a) match the values of Sigurdsson and Phinney [17, table 3B]. This is necessary as the total cross-section we calculate only takes close encounters into account, while the total cross-section of Sigurdsson and Phinney [17] also includes weak interactions.

FIG. 3. A comparison between equation (48) and the numerical simulations of Sigurdsson and Phinney
radius is [19] Before proceeding, let us note that equation (56), together with the form of f bin , suffices to explain the main features of numerical simulations: inserting equation (56) into equation (60), implies that where the upper ∼ signs indicate that the energy-shifts and their associated probabilities may be different from those obtained in a single action of W lim , but the structure is still the same. The final energy is 'traced out', while the initial total energy is fixed.
Below we start by describing the tidal model we adopt in this paper, then we compute P (x) perturbatively, and then compute the cross-section for a collision and a tidal in-spiral. The latter event is an in-spiral of two stars into one another due to orbital energy loss to tides. We back our analysis up by comparing its results to the numerical simulations of ref. [28].

A. Tidal Model
We adopt the tidal model of ref. [32]. We further deem any energy that goes into the tidal oscillations of the white dwarf as lost from the system (that is, the timescale on which it might return to orbital energy is much larger than the relevant dynamical time-scales), and we approximate the total angular momentum as fixed. 8 We take a tidal dissipation event in a single close approach to occur with probability where 1 < y is a free parameter, which describes, roughly, the maximum separation between two bodies which would engender sizeable tidal effects. Equation (64) is strictly correct in the limit where yr * a bin [27], which is the limit we consider here; at larger yr * it has to be modified [17,27]. It originates from the following reasoning: let u < y. For ur * a bin , the cross-section for star 1 (which we choose to be the white dwarf) to come within a distance ur * of one of the other two stars is approximately described by a two-body interaction with either one of them. Gravitational focussing thus implies that the cross-section for this event is 4πGmur * /v 2 , where v 2 is the "initial" velocity of star 1. As the tidal interaction occurs during the chaotic 3body-close-interaction phase, v should be, roughly, given by the virial speed, multiplied by √ 2 because it is a rel-ative velocity, i.e. v 2 ≈ 4 3 |E| m , with E ≈ − Gm 2 2a0 , since the original binary was hard. The cross-section should be divided by the total area available for star 1, which is approximately πR 2 , and multiplied by two to account for the two possible partners 2 and 3. Thus, as in equation (63). p tide (a bin , u) is therefore the probability of star 1 coming within ur * from star 2 or star 3. The probability density function is (65) If such an event does occur (i.e. if star 1 comes to (u + du)r * from either of its companions, but not below ur * ), the energy loss is given by (see ref. [32], with T 2 (x) ∼ To gauge the error on our computations, we also use the simpler model of ref. [34], in which Needless to say, the technique described in §VI applies to more sophisticated tidal models, too.

B. Perturbative Calculation of P (x)
Let us define the following functions of the total triple energy E (and, albeit suppressed, total angular momentum J), for an initial total energy E i (please bear in mind that ∆E < 0): Suppose that r * a bin . Then, p tide 1, and most of the steps conserve the total binary energy. Then one can view this process as two random walks: a "small-scale" random walk in E bin , and a larger-scale one in E, on top of it, and derive the equation perturbatively. What this means is that, in expanding P n , one may truncate the series at some small power of ε = r * |Ei| GM µs , where E i is the initial total energy. Such a truncation corresponds to re-summing the series of P (x), as a series in powers of ε. Suppose one stops at first order; then The first line in this equation corresponds to the occurrence of no tidal interactions, the last -to a tidal interaction in the very first close approach, and the second -to a tidal interaction in another close approach.
The sum over k is a geometric sum, and may be computed analytically. Then, using equation (60) (summing from n = 1, as we assume that the first, initial close approach always happens), one may sum over n analytically, too (this sum may be exchanged with the action of W unlim by its linearity). Then one may act with W unlim to find that the probability distribution for the final binary energy reads up to an overall normalisation and O(ε 2 ) corrections.

C. In-Spiral Cross-Section
Samsing et al. [28] define the result of an encounter to be deemed 'an in-spiral' if a bin ≤ 6r * and if a collision has not occurred. The final possible outcomes are therefore: a collision, an in-spiral, an exchange or a flyby. The exchange/fly-by cross-section may be computed by including a Heaviside function in p, q, p u , q u which ensures that there is no collision and no in-spiral. The in-spiral/collision cross-section is simply the total crosssection, minus the exchange/fly-by cross-section.
There is another possible outcome: if ∆E(u) is large enough relative to E i , then for large enough J, upon losing energy to tides, there is too much angular momentum for the triple system to interact closely again -this manifests itself in none of the conditions in §II A being satisfied. Denote the minimum such J by J * . In this case, the triple becomes hierarchical automatically, and the encounter ends. The inevitable fate of this triple is a tidal in-spiral: during each pericentre approach of the inner binary more energy is lost to tides, until the two stars collide. This is just another route to a tidal inspiral, which does not require the triple to have ejected the third star.
As in both models ∆E(u) decays quite fast with u, the dependence on y is very weak, and the cross-sections converge for sufficiently large y (see figure 5).
One may compute the tidal in-spiral cross-section as follows: compute the cross-section for there being no col-lision and no tidal in-spiral, and then subtract the result from the total cross-section, which we take, in this section, to be -and not twice this value -to be consistent with ref. [28]. The former is given by first computing the differential cross-section for a final periapsis r p , by integrating equation (71) over the disc D(v 0 ) -this is effectively an integration over the allowed range of total angular momenta -as well as over the binary actions with a Dirac delta function δ(r p − a bin (1 − e bin )), all the while enforcing the no-in-spiral-and-no-collisions condition both in equation (71)  bin , a bin , and the distribution of a bin is peaked sharply about a 0 . The total exchange/fly-by cross-section is then found by integrating this differential cross-section over r p ≥ r * (to preclude collisions); by the special nature of equation (62), this automatically precludes collisions during all close approaches, not just the final one.
Explicitly, let these quantities, when evaluated at E = E i − ∆E, are only defined for J ≤ J * , so for J > J * we define them to be zero. Further define where J(b) is the magnitude of the total angular momentum as a function of the position b of the in-coming star on the disc D(v 0 ). It follows from equation (71) that if we denote  (67), while the blue line uses equation (66); the red line is the numerical result of Samsing et al. [28]. As expected, the cross-section converges to a constant for large-enough y, despite some noise due to numerical integration.
where p tide,nc (a bin ) = 12a 0 (y −1)r * /a 2 bin , to preclude collisions; p ni (E) is the analogue of p(E), but with the same Heaviside theta function inserted as in the definition of N ni ubd and with p tide,nc (a bin ) used instead of p tide ; and q ni (E) is the same as q(E), but also with the afore-mentioned Heaviside function,-if so, then the inspiral/collision cross-section is given by We find, from equation (78), that for β ≈ 1.3, σ insp+coll , agrees with the findings of ref. [28] for both a 0 = 10 −3 and 10 −2 AU. The former improves upon their analytic estimate of 0.155 AU 2 considerably. This crosssection depends on β, as shown in figure 6. The fact that β = 1.3 fits both cases, with initial semi-major axes differing by an order of magnitude, implies that indeed equation (78) is an adequate model to describe a binarysingle encounter with tides.

IX. DISCUSSION AND SUMMARY
In this paper we introduced a random walk model for binary-single encounters in globular clusters. An encounter is viewed as a sequence of chaotic, close triple approaches, interspersed with hierarchical phases. The orbital parameters of the binary and the single star, as well as the triple system's constants of motion, perform a random walk: each step of the walk corresponds to one close approach and its subsequent hierarchical phase. We calculated the transition probabilities between steps of the walk, both in the Newtonian, point-mass approximation (in which the walk is memory-less, and the final outcome distribution is simply given by the formula of ref. [26], while the transition probabilities between intermediate steps are given by equation (33)), and in the case where there is some dissipative process involved, when equation (62) holds.
We have shown that this model reproduces numerical results well, as we exemplified for aspects such as the semi-major axis distribution, the escaper's mass distribution, and the final periapsis distribution. Including tides and collisions, our predictions match the inspiral/collision cross-section measured by numerical simulations. This validates the prescription of R and the solution to the bound problem; besides, including a tides allowed us to perform a non-trivial test on the extra step involved in elevating the statistical solution of the scattering problem to the random walk model, which it passed.
In some cases the probability distribution can be computed completely analytically, while in others it only involves a relatively simple calculation of a few integrals. Please note, that the formula for the escaper's mass distribution is valid in every intermediate step of the encounter. In conjunction with the random walk model, it implies, inter alia, that if one of the stars is considerably lighter than the other two, then it will be the one to be ejected after each close approach -not just the final one. The only free parameters in our analysis, which do not influence the results of the dissipation-free problem significantly, are β and η; the former of which was found to be equal to 1.3 -a single value that agrees with both cases considered in §VIII, and the latter does not change There is some noise due to numerical integration. Top: the case a0 = 10 −3 AU; bottom: the case a0 = 10 −2 AU. In both cases v0 = 10 km s −1 , and y is chosen to be very large.
the results much in either case. The fact that one value of β fits both supports the random walk approach further. The random walk model described here may be used to address a plethora of astrophysical phenomena using the analytical, statistical model we described in this paper: apart from incorporating tidal interactions as was done here, one could also include gravitational-wave dissipation, and the effects of stellar evolution. One could also investigate external effects, for example the tidal influence of an external gravitational potential. The main change would be a modification of the largest value a s (1 + e s ) that corresponds to a bound binary. Such cut-offs are important in globular clusters, which are the main places where one expects to find significant rates of binary-single encounters. This is made possible by our statistical approximate solution of the bound, nonhierarchical, three-body problem in equation (33). This approach also allows one to calculate the distribution of the number of consecutive close approaches before the disruption of the triple, and hence the distribution of time-scales of temporary captures which could be relevant for various astrophysical capture processes by gravitating stars and planets. Other applications include -but are not limited to -binary-single encounters in nuclear star clusters around super-massive black holes at the centres of galaxies (or, equivalently, binary-single encounters of planets, dwarf planets, moons or asteroids in the solar system), where the Hill radius effectively provides a limiting separation (as mentioned above) during consecutive encounters, as well as the velocity distribution of fast, runaway stars due to ejections through binary-single encounters, which will be treated in a separate paper.
Let us present an algorithm for how this would be done for any given astrophysical process involved in threebody physics. Suppose, that in addition to Newtonian, point-particle motion one wishes to incorporate another physical phenomenon, say, gravitational-wave emission (which we will address in future work). One would have to be able to calculate how this phenomenon changes the total energy and angular momentum E, J during a single close approach and subsequent hierarchical phase. With these data, one could compute f c (E, J|E , E bin , S , J ) from §VI. Using the random walk model, it gives the transition probabilities h(x|x ) immediately, and then all that remains is to sum up the series (60) to obtain the final probability -the probability density function of final binary parameters, given initial total energy and angular momentum. If there is a small parameter, e.g. if the probability for non-zero change in the constants of motion is small, one can re-sum equation (60), and expand in the small parameter. To summarise: 1. Compute the single-step (one close approach and one hierarchical phase) probabilities for a given change in total energy and angular momentum, f c (E, J|E , E bin , S , J ); 2. compute the transition probabilities h(x|x ) as in §VI; 3. and compute equation (60) to obtain the final distribution.
In our paper we brought the endeavour of statistical modelling of binary-single encounters closer to completion. By viewing them as concatenations of close triple approaches, we were able to model, statistically, the bound non-hierarchical three-body problem, and to use the solution to model the entire encounter as a random walk. This model has the potential to facilitate simulations of globular clusters significantly. Instead of having to resolve binary-single encounters individually using a high-resolution few-body code, one could simply implement the analytical random-walk model as a probabilistic solution of these encounters, with the additional possibility of incorporating any astrophysical process one wants. By the law of large numbers, if the number of encounters per cluster is sufficiently large, the few-body resolutions of binary-single encounters will be rendered unnecessary. We hope that this gain in speed would enable astrophysicists to study many more phenomena in clusters, planetary systems and in the field in much better detail.

ACKNOWLEDGMENTS
To compute I (2) − one needs to consider the signs carefully. The result is Inserting R from equation (10), where 1 β 2, the exact value of I (α(E bin )) follows, roughly, a power-law until J = α, where it begins to fall sharply (like I + ). The dominant contribution to I − is I (1) − , whence, when |E bin | |E|, but α > J, one may approximate I by I − , which goes like |E bin | −1 for high angular momentum, but like |E bin | −1/2 for low angular momentum. In the unbound case, when |E bin | |E|, one has a s ≈ a bin m s /µ bin , and An even better approximation would be to take into account the eccentricity dependence of θ max , when computing I . The angle θ max is equal to θ ap for e s = 1, and it vanishes when L 2 saturates the periapsis bound of §II A. Let us define ξ = R/a s and so that Then θ max (b = 0) = θ ap , while θ max (b = 1) = 0. One could write θ max = θ ap θmax θap , and then approximate the fraction. While at b = 0, ξ = 0, this fraction is unity, it differs from 1 for non-zero b even at ξ = 0. So, let us define expand θmax θap in powers of ξ; keeping only the leading term we find uniformly in b. We keep only the leading term, and, since we have shown above that I ≈ I − for most values of E bin , except possibly those where f bin is very small anyway, we will only compute a correction to I (1) − here, but corrections to I + and I (2) − may be obtained in a similar manner. The integral we need to compute is Changing variables from x to L gives Fortunately, this integral may be computed analytically in terms of inverse trigonometric functions. Let φ(u) = 1 20 1 − u 2 (12 + u 2 (1 + 2u 2 )) + 3 4 u arcsin u, the Lyapunov time λ −1 , and as long as the system is in C. Suppose that one starts with a phase-space distribution that is uniform on some ε-neighbourhood of some w 0 ∈ C, namely uniform in R = B ε (w 0 ) ⊆ C; then, after a time t, this phase-space density evolves to a uniform distribution in a sphere (in the metric in which the Lyapunov exponent is calculated) of volume ∼ V (B ε (w 0 ))e λt . Suppose now that we wish to start with an initial condition in C very close to a delta-function. Then, as time goes by, f spreads over C, until some parts of it reach C's boundaries, and enter A or B. The time-scale of evolution in C is roughly the virial time-scale of the three-body system, while the time-scales for interesting evolution in A and B are set by the frequencies Ω of the outer binary. Per definitionem, these are much longer than the virial time-scale, for otherwise the system would not be hierarchical. Thus, those parts of f which have left C, are, from the point of view of the C, stuck at the boundary, so to speak. The distribution continues to spread, until virtually all of f is on ∂C. The proportion that arrives at A never returns to C, but the rest -being in B, eventually does return to C, and starts all over. This goes on ad infinitum, or until all of the probability mass is in A.
Separation of scales thus ensures that the evolution of the system proceeds as a sequence of close three-body, chaotic encounters, and between them -hierarchical phases, until one of the three bodies is ejected. This picture meshes well with simulations [11,12,19], but here we have lent it some theoretical credence.
If w 0 is a typical initial condition inside C, then its phase-space distance from the ∂C is roughly R. According to the evolution described above, the time it would take the system to arrive at ∂C is then The power of 8 is due to the conserved quantities: the three-body phase-space is 18-dimensional, but conservation of linear momentum in the centre-of-mass frame, angular momentum and energy reduce this dimension to 8. Requiring a resolution ε/R 1 implies that the time it takes the system to leave C is larger than the Lyapunov time, consolidating the assumption of efficient phase-space mixing inside it. The simulations of ref. [21] demonstrated that the half-life time of chaotic three-body systems is (for the equal mass case, in units of τ vir , see their table 4). This is indeed sufficient for chaotic mixing. The next step is to see how much of the probability mass arrives at each point in ∂C in a single close encounter. (In a close approach that is not the first one, by linearity, the solution is a superposition of solutions that start out as delta-functions around some point in C.) Due to the scale-separation, the evolution equation for the distribution function in C should have (approximately) Dirichlet boundary conditions on ∂C, whence the probability of reaching a point w ∈ ∂C -which leads immediately to the outcome probability in A and B via Liouville's equation -is given by the so-called 'eventual hitting probability' of w [37], which is simply the time-integral of the scalar product of ∇f with the unit normal to ∂C,n (the gradient ∇ is a phase-space gradient). As f effectively evolves like a top hat that expands, ∇ f dt is very large at ∂C, but its magnitude is independent of w 0 . The scalar productn · ∇f gives an additional cosine, so that f eventual (w) ∝ cos θ(w, w 0 ), where θ(w, w 0 ) is the angle betweenn and the vector pointing from w 0 to w. However, the boundary between C and the other two sets is not well-defined, so instead, it would be useful to think of ∂C as a 'fuzzy' boundary -with some width. This implies that the cosine should be removed as un-physical, and replaced by its average π/2 −π/2 cos θdθ. Thus, the eventual hitting probability is well-approximated by a function independent of both w 0 and w.
What one is interested in is the probability distribution of an outcome in A or B, and specifically, in the action distribution there. Consider the set S(J a , J b , J c ), which is the set of all points w ∈ A ∪ B such that the (Delaunay) actions of the inner binary are between (J a , J b , J c ) and (J a + dJ a , J b + dJ b , J c + dJ c ). What we are interested in is the measure of S(J a , J b , J c ). Given that these actions correspond to something the has been in a close triple system, one can use Liouville's theorem to find translate this question into a question of finding the probability distribution function of ∂C. This is independent of the initial condition, which implies that one can simply compute it assuming a uniform distribution of initial conditions in C, which we compute in §II.
The reader should also bear in mind that due to the scale separation of the evolution in the hierarchical region B, there is an additional source of phase-space mixing: while the outer body moves along its two-body orbit about the inner binary, the mean anomaly of the latter evolves much more rapidly, and its value when the third star returns depends on its initial value sensitively, which implies that it is effectively quite random. Therefore, as the number of close approaches increases, σ should resemble equation (2) more closely -this fact is attested to by the simulations of ref. [26]. Indeed, if the above arguments apply only approximately, so that in a single close approach mixing is only efficient in a fraction ν of the total volume of C, then the phase-space distribution still tends to a fully-mixed one, as ν n -exponentially in the number n of close approaches. It may be possible that a close triple approach following a hierarchical phase end more quickly than is required for the system to mix chaotically in C [25]. In that case, this 'phase-mixing' implies that the cross-section is still mixed.

Appendix C: Approximation Error Induced By Model
As the solution described in the previous sections is a statistical approximation to a deterministic, albeit chaotic, problem, it would be useful to be able to constrain the uncertainty due to the approximation above. We do so by evaluating f bin for different values of β, and comparing the resulting plots. An analogous check may be performed for η. 10 While this test is not an ideal uncertainty estimate, it is at least a test of the robustness of the model; we rely on the fact that in many situations, the error induced by varying the parameters of a model is of a similar size to the one made by using that model as an approximation. Figure 8 shows the marginal energy distribution, and figure 9 shows contour plots of the joint distribution of E bin and S = |S|; both figures indicate a weak dependence on β, and therefore a high accuracy of the approximation model made here.

FIG. 10.
A comparison between f bin , calculated with R defined by equation (10), and by equation (D3). The first two rows show the marginal energy distribution, and the last line shows that joint energy-spin distribution (left: equation (10); right: equation (D3)). The top row displays f bin for a large value of J, while the second -for a low value. For the bottom row we used the same value of J as for the top row, and equal masses. All rows have β = 1.5.