# Illustrating chaos: a schematic discretization of the general three-body problem in Newtonian gravity

Illustrating chaos: a schematic discretization of the general three-body problem in Newtonian... Abstract We present a formalism for constructing schematic diagrams to depict chaotic three-body interactions in Newtonian gravity. This is done by decomposing each interaction into a series of discrete transformations in energy- and angular momentum-space. Each time a transformation is applied, the system changes state as the particles re-distribute their energy and angular momenta. These diagrams have the virtue of containing all of the quantitative information needed to fully characterize most bound or unbound interactions through time and space, including the total duration of the interaction, the initial and final stable states in addition to every intervening temporary meta-stable state. As shown via an illustrative example for the bound case, prolonged excursions of one of the particles, which by far dominates the computational cost of the simulations, are reduced to a single discrete transformation in energy- and angular momentum-space, thereby potentially mitigating any computational expense. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening. Finally, we provide a method for dynamically evolving entire populations of binaries via three-body scattering interactions, using a purely analytic formalism. In principle, the techniques presented here are adaptable to other three-body problems that conserve energy and angular momentum. gravitation, scattering, methods: analytical, binaries: close, stars: kinematics and dynamics, globular clusters: general 1 INTRODUCTION After hundreds of years of study (e.g. Newton 1686), the chaotic three-body problem remains pertinent to the field of stellar dynamics, as it relates to our understanding of interactions within globular (e.g. Heggie 1975; Leigh & Geller 2012; Leigh et al. 2013; Antonini et al. 2016), open (e.g. Leonard 1989; Hurley et al. 2005; Leigh & Sills 2011; Geller, Hurley & Mathieu 2013; Leigh & Geller 2013; Mapelli & Zampieri 2014; Geller & Leigh 2015; Leigh et al. 2016), and even nuclear (e.g. Davies et al. 1998; Merritt 2013) star clusters. The problem remains analytically unsolved. Historically, any solutions involved approximate methods or the application of constraining assumptions used to simplify the problem. The best known such example of this approach is the restricted three-body problem, as applied to the Earth–Sun–Moon system. But modern computer technologies have come a long way and, in so doing, sparked a new wave of intensive study of such systems over the last few decades. Most simulations for the chaotic three-body problem run very quickly on modern computers. Every now and then, however, you will encounter one simulation that gets hung up, and takes many orders of magnitude more computer time to run to completion. Thus, upon performing a suite of numerical experiments of single–binary scattering, for example, it is a handful of individual runs that can consume the vast majority of the total computational cost. For small tidal tolerance parameters (see Section 2.2 for more details), this is due to prolonged excursions, where one of the stars ends up on a near-zero energy orbit. In fact, some authors have argued that the mean disruption time should formally be infinite for resonant three-body interactions (e.g. Valtonen 1975; Valtonen & Karttunen 2006; Orlov, Rubinov & Shevchenko 2010; Shevchenko 2010; Leigh et al. 2016); if an infinite number of simulations could be performed, the probability of obtaining one or two simulations with near infinite excursion time events becomes finite. These individual events can be bottlenecks in studies of the chaotic three-body problem, and compound to become even worse when additional particles are included into the mix. Chaotic three-body systems interacting resonantly display a characteristic evolution in time that begins when one particle is sent on a temporary excursion while remaining bound to the system. Upon returning to periastron, the particle initiates another close triple encounter, interacts strongly with the other two particles, and the process repeats (e.g. Agekyan & Anosova 1967; Anosova & Orlov 1994). To refer to events where the ejected particle is sent on a bound excursion we will use the term ‘ejection’. Conversely, to refer to events where the particle becomes unbound from the system and does not return, we will use the term ‘escape’. The end result of this chaotic dance is always an escape event where one of the bodies is ejected with a sufficiently high velocity to become unbound, and escapes to spatial infinity. The chaotic progression or evolution described above can be simplified or reduced to a single discrete transformation in energy and angular momentum space. The initial and final states of the system are qualitatively the same, since the unbound system can be decomposed into a single star and a binary. Once the initial relative energies and angular momenta of the single star and the binary are decided (via the input parameters that define a given scattering experiment), a transformation is applied. This changes the relative energies and angular momenta such that they are distributed differently among the particles in the initial and final states. By neglecting the intermediary chaos that occurs in the bound state, statistical ensembles of outcomes can be considered directly, rather than relying on specific choices of scattering parameters. In this paper, we present a schematic formalism to discretize the general three-body problem in Newtonian gravity, in both the bound and unbound states. This is done in such a way that, for most cases of interest, our diagrams fully and quantitatively characterize the time evolution of a given chaotic single–binary interaction. We use the numerical scattering code fewbody to simulate a series of binary-single encounters involving identical point particles in Section 2. We choose one particularly interesting example, that features a very prolonged excursion of one of the stars, and generate our schematic diagrams to depict the time evolution of this interaction. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening, and even describe a purely analytic method for dynamically evolving entire populations of binaries via three-body scattering interactions. In Section 4, we discuss the utility of our schematic formalism, address possible shortcomings and applications of our method and explain how our formalism might eventually replace the need for computer simulations of single–binary scatterings. We summarize in Section 5. 2 METHOD In this section, we introduce our formalism for constructing schematic diagrams of the general three-body problem, adapted from Leigh & Geller (2012) and Leigh et al. (2016). We further present the numerical scattering experiments used to study the time evolution of the chaotic three-body problem, from which we construct an illustrative example of our schematic diagram formalism. 2.1 Constructing the schematic diagrams In this section, we begin by adapting the schematic diagrams first presented in Leigh & Geller (2012) to depict individual snapshots of the system in time for application to the general three-body problem, including arbitrarily long excursions of one of the particles, followed by a description of our method for discretizing the interaction. First, the total energy of the three-body system can be written as   $$E = \sum _{\rm i=1}^{\rm N=3} E_{\rm i}^{\prime } = \sum _{\rm i=1}^{\rm N=3} \Big ( \frac{1}{2}m_{\rm i}v_{\rm i}^2 - \frac{1}{2}\sum _{\rm j=1, j \ne i}^{N=3} \frac{Gm_{\rm i}m_{\rm j}}{r_{\rm ij}} \Big ),$$ (1)where mi is the mass of the ith particle, vi is its velocity relative to the system centre of mass, ri is the magnitude of its distance from the system centre of mass, and $$r_{\rm ij} = |\bar{r}_{\rm i} - \bar{r}_{\rm j}|$$. Re-writing equation (1), we get   $$E - E_{\rm 0} = \sum _{\rm i=1}^{\rm N=3} \Big (E_{\rm i}^{\prime } - E_{\rm 0}/3 \Big ) = \sum _{\rm i=1}^{\rm N=3} E_{\rm i},$$ (2)where E0 is a negative constant with units of energy. The addition of the constant E0 serves only to shift the zero-point of the total encounter energy. This is needed to avoid the appearance of negative angles in our formalism. We adopt this approach for the bound case in Section 3.1. In Section 3.2 we consider a series of unbound three-body scattering interactions, and explore a different procedure for handling negative angles. Here, negative angles instead correspond to unbound particles. We will return to this in Section 3.2. Next, we point out that the sum of the angles of any polygon always add to a constant. For example, the sum of the angles of a triangle sum to 180°. Similarly, the energy contained in each particle Ei must always add up to equal the (shifted) total encounter energy E − E0. Hence, we can write   $$\frac{E - E_{\rm 0}}{180^{\circ }} = \frac{E_{\rm i}}{\theta _{\rm i}},$$ (3)where θi denotes the angle of a triangle corresponding to particle i with total energy Ei, as given by equation (2). For each triangle, we fix the length of the side of the triangle oriented along the x-axis (i.e. parallel to the direction of increasing time) to be equal to unity. This ensures that, in energy-space, the symmetry of each triangle is positively correlated with the duration of the corresponding excursion event. More equilateral triangles correspond to longer excursion events (see Section 3 for more details). A similar exercise can be performed to generate triangles for each component of the total angular momentum, namely Lx, Ly, and Lz. Using the x-component of the total angular momentum as an example we can write   $$L_{\rm x} - L_{\rm x,0} = \sum _{\rm i=1}^{\rm N=3} \Big (L_{\rm x,i}^{\prime } - L_{\rm x,0}/3 \Big ) = \sum _{\rm i=1}^{\rm N=3} L_{\rm x,i},$$ (4)and hence:   $$\frac{L_{\rm x} - L_{\rm x,0}}{180^{\circ }} = \frac{L_{\rm x,i}}{\theta _{\rm x,i}},$$ (5)where θx, i denotes the angle of a triangle corresponding to particle i with an x-component of the total angular momentum Lx, i, as given by equation (4). Next, we describe our formalism for discretizing the time evolution of the interaction. As described in Section 1, the evolving three-body system can usually be decomposed into a temporary binary, with the other star acting as a temporary single star on an excursion. This single–binary pair forms a closed (bound) temporary orbit. In order to formally define this criterion, we output from the simulations described in the subsequent section the critical time-steps at which the radial component of the centre of mass velocity of the temporary binary changes sign from positive to negative. This always coincides with the time-step at which the radial component of the centre of mass velocity of the temporary single star also changes sign. Every time this transition occurs, we record the energies of the particles and all three components of their total angular momentum. We emphasize that, since there are three particles, in this scenario there will always be one particle that is furthest from the system centre of mass, and is closest to the E ∼ 0 boundary. Using the formalism described above for generating triangles that quantitatively describe the temporary system in energy- and angular momentum-space, we then construct four triangles for each time-step corresponding to this change in sign of the radial velocities, one for energy and three for the total angular momentum (one for each component). The angle corresponding to each vertex of a given triangle is directly proportional to the energy or angular momentum content of the indicated particle, as described above. The angle corresponding to the top vertex is defined to always correspond to the temporary single star. The triangles are arranged sequentially from left to right, to indicate the direction of increasing time. The number of triangles is one greater than the number of applied transformations, to include the initial state of the system. This number indicates the total number of temporary excursion events, and clearly reveals any very long-lived computationally expensive excursions (via the relative symmetry of the triangles; see above and Section 3). 2.2 Numerical scattering experiments We calculate the outcomes of a series of single–binary (1+2) encounters using the fewbody numerical scattering code.1 The code integrates the usual N-body equations in configuration- (i.e. position-) space in order to advance the system forward in time, using the eighth-order Runge–Kutta Prince–Dormand integration method with ninth-order error estimate and adaptive time-step. For more details about the fewbody code, we refer the reader to Fregeau et al. (2004). All objects are point particles with masses of 1 M⊙. The initial binary has ai = 1 au, and eccentricity ei = 0. We fix the impact parameter at b = 0 for all simulations. The angles defining the initial relative configurations of the binary orbital plane and phases relative to the velocity vector of the incoming single star are chosen at random. We perform 1000 numerical scattering experiments to find a suitable illustrative example, including a prolonged excursion of one of the particles. This is illustrated in Fig. 1. Figure 1. View largeDownload slide The time evolution of an example single–binary simulation in position-space, projected into two dimensions. The unit of distance is 1 au on both axes. Notice the prolonged excursion of one of the particles, which dominates the total integration time of the simulation. Figure 1. View largeDownload slide The time evolution of an example single–binary simulation in position-space, projected into two dimensions. The unit of distance is 1 au on both axes. Notice the prolonged excursion of one of the particles, which dominates the total integration time of the simulation. To decide when an encounter is complete, we adopt the same criteria as Fregeau et al. (2004). To first order, we define this as the time when the separately bound hierarchies that constitute the system are no longer interacting with each other, or evolving internally. Formally, the integration is terminated when the top-level hierarchies have positive relative velocity and the corresponding top-level N-body system has positive total energy. We set the tidal tolerance parameter to δ = 10−7 for all simulations, which has been previously shown to be sufficient for the virial ratios (defined as the ratio between the total system kinetic energy and the total potential energy) considered here (see Geller & Leigh 2015 for more details). 3 RESULTS In this section, we generate and present our schematic diagrams for the discretization of the general three-body problem. We begin with a single three-body interaction that remains bound and in a resonant state for several system crossing times, and apply our method to the illustrative example shown in Fig. 1. We then present a modified version of our formalism to treat sequences of unbound single–binary scattering interactions, as occur during binary hardening in dense stellar environments, and show how (in some cases) our schematic diagrams can be generated purely analytically without any numerical simulations. Using this basic approach, we develop an analytic formalism for dynamically evolving entire populations of binaries in dense stellar environments due to three-body scattering interactions. 3.1 Bound case Schematic diagrams for the single–binary interaction depicted in Fig. 1 are shown in Fig. 2. Four separate triangle series are shown, one for energy-space (top) and three for angular momentum (bottom) with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. For each triangle, we fix the length of the side of the triangle oriented along the x-axis (i.e. parallel to the direction of increasing time). Combined with the imposed offset in energy E0, this ensures that equilateral configurations of our triangles in energy-space are positively correlated with the duration of the corresponding excursion event. The initial and final states correspond to, respectively, the first and last triangles. Hence, the progression of time is from left to right. Figure 2. View largeDownload slide The time evolution of the example single–binary simulation shown in Fig. 1 in both energy- and angular momentum-space, as described in the text. Four separate triangle series are shown, one for energy (top) and three for angular momentum (bottom), with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. At the vertex of the apex of each triangle opposing the horizontal axis, we indicate which particle corresponds to the temporary single star, using the same choice of labelling as in Fig. 1. Left to right denotes the direction of increasing time. We set E0 = 6.6 × 1038 J, Lx0 = 4.0 × 1048 kg m2 s−1, Ly0 = 4.0 × 1048 kg m2 s−1, and Lz0 = 2.0 × 1048 kg m2 s−1 for the zero-points, or equivalently E0 = 10|E|, Lx0 = 2300Lx, Ly0 = 510Ly, and Lz0 = 1300Lz. Figure 2. View largeDownload slide The time evolution of the example single–binary simulation shown in Fig. 1 in both energy- and angular momentum-space, as described in the text. Four separate triangle series are shown, one for energy (top) and three for angular momentum (bottom), with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. At the vertex of the apex of each triangle opposing the horizontal axis, we indicate which particle corresponds to the temporary single star, using the same choice of labelling as in Fig. 1. Left to right denotes the direction of increasing time. We set E0 = 6.6 × 1038 J, Lx0 = 4.0 × 1048 kg m2 s−1, Ly0 = 4.0 × 1048 kg m2 s−1, and Lz0 = 2.0 × 1048 kg m2 s−1 for the zero-points, or equivalently E0 = 10|E|, Lx0 = 2300Lx, Ly0 = 510Ly, and Lz0 = 1300Lz. The diagrams shown in Fig. 2 fully and quantitatively describe the time evolution of the example single–binary interaction depicted in Fig. 1. In principle, using these diagrams, the time evolution of the interaction in position- and velocity-space can be re-constructed. This is because the energy and angular momentum of each particle is quantitatively encoded in the diagrams at each discrete exchange event via equations (3) and (5) combined with equations (1), (2), and (4), and the duration of each corresponding excursion is quantitatively encoded by the shape of each triangle (see below). Specifically, these diagrams encode 12 equations with 18 unknowns. The remaining six equations are found from the symmetry properties of the initial, final, and temporary binaries (i.e. given the provided information in our diagrams, only the positions and velocities of one of the binary components are needed to compute those of the other). This is done as follows. Consider a binary with component masses m1 and m2, oriented in the xy-plane. We transform our coordinate system such that it is stationary in the centre-of-mass frame of reference of the binary. Then, if the (instantaneous) positions and velocities of the component with mass m1 are (x1, y1, z1) and (vx, 1, vy, 1, vz, 1), then the positions and velocities of the second component are   \begin{eqnarray} x_{\rm 2} = -\frac{m_{\rm 1}}{m_{\rm 2}}x_{\rm 1}, \end{eqnarray} (6)  \begin{eqnarray} y_{\rm 2} = -\frac{m_{\rm 1}}{m_{\rm 2}}y_{\rm 1}, \end{eqnarray} (7)  \begin{eqnarray} z_{\rm 2} = z_{\rm 1} = 0, \end{eqnarray} (8)  \begin{eqnarray} v_{\rm x,2} = -\frac{m_{\rm 1}}{m_{\rm 2}}v_{\rm x,1}, \end{eqnarray} (9)  \begin{eqnarray} v_{\rm y,2} = -\frac{m_{\rm 1}}{m_{\rm 2}}v_{\rm y,1}, \end{eqnarray} (10)  \begin{eqnarray} v_{\rm z,2} = v_{\rm z,1} = 0. \end{eqnarray} (11)This follows from the definition of the system centre of mass, which defines the equation:   $$m_{\rm 1}\boldsymbol {{R}}_{\rm 1} + m_{\rm 2}\boldsymbol {{R}}_{\rm 2} = 0,$$ (12)where $$\boldsymbol {{R}}_{\rm 1}$$ and $$\boldsymbol {{R}}_{\rm 2}$$ are the radius vectors of the bodies with respect to the binary centre of mass, and   $$|\boldsymbol {{R}}_{i}| = \sqrt{x_{\rm i}^2 + y_{\rm i}^2 + z_{\rm i}^2},$$ (13)and i = 1 or 2. Now, consider the very prolonged excursion event indicated in Fig. 1. This corresponds to the third triangle (from the left) in Fig. 2. This is very close to a completely equilateral triangle in energy-space, to within fractions of a degree, whereas triangles corresponding to shorter excursion events are much less symmetric in energy-space. This behaviour is also evident by the corresponding triangles in Lx-, Ly-, and Lz-space, in which one particle clearly dominates the angular momentum budget during the excursion. We define the following parameter:   $$\epsilon = \sum _{\rm i=1}^{3} \left|\theta _{\rm i} - \frac{\pi }{3}\right|.$$ (14)The longest excursion events correspond to ε ∼ 0, whereas shorter excursion events satisfy ε > 0. The reason for ε ∼ 0 corresponding to the longest excursions is that, here, the temporary single star is nearly unbound such that its energy is very close to zero. Meanwhile, the temporary binary is approximately isolated, such that it (approximately) obeys the virial theorem. Hence, twice the total kinetic energy is roughly equal to the absolute value of the total potential energy in the temporary binary. This translates into E΄ ∼ 0 in equation (1), and hence all θi are approximately equal to π/3 via equations (2) and (3). 3.2 Unbound case: one binary In this section, we briefly reproduce the analytic formalism for unbound scattering presented in Valtonen & Karttunen (2006), in particular the key formulae and distribution functions to be used and discussed in the subsequent sections. These are presented here via an example application of our schematic formalism to the unbound case. In dense stellar environments, as occur in the cores of star clusters, a single–binary star can undergo repeated interactions with many different single stars. This sequence of single–binary scatterings typically acts to harden the binary (i.e. reduce its orbital separation). Our schematic diagrams can be adapted to quantitatively describe this process. The required distribution functions are already available in analytic form for the initial and final unbound states (i.e. the first and last triangles) (Valtonen & Karttunen 2006). This is thanks to pioneering efforts by Monaghan (1976a) and Monaghan (1976b), who re-worked the expression for the density of configuration states per unit energy in order to derive the distributions of binary orbital energies, eccentricities, and escaper velocities in the low- and high-angular momentum limits. This was later parametrized to be applicable to any total angular momentum (Valtonen & Karttunen 2006). Thus, we can already generate the initial and final states from these distribution functions. One way to do this is as follows. We write for the total encounter energy of a given single–binary interaction:   $$E = E_{\rm s} + E_{\rm B},$$ (15)where EB is the binary orbital energy and $$E_{\rm s} = \frac{1}{2}m_{\rm s}v_{\rm s}^2$$ is the kinetic energy of the single star in the centre of mass frame of the binary. From equation 7.19 in Valtonen & Karttunen (2006), the distribution of escaper velocities is   $$f(v_{\rm s}){\rm d}v_{\rm s} = \frac{(n-1)|E|^{n-1}(m_{\rm s}M/m_{\rm B}))v_{\rm s}{\rm d}v_{\rm s}}{(|E| + \frac{1}{2}(m_{\rm s}M/m_{\rm B})v_{\rm s}^2)^{n}},$$ (16)where ms is the mass of the single star, mB is the mass of the binary, and M = ms + mB is the total system mass. From equation 7.26, we have for the final distribution of binary orbital energies:   $$f(|E_{\rm B}|){{\rm d}}|E_{\rm B}| = (n-1)|E|^{n-1}|E_{\rm B}|^{-n}{{\rm d}}|E_{\rm B}|.$$ (17)Equations (16) and (17) can be sampled from directly,2 and hence final states constructed from the initial state. The parameter n corrects each distribution for the total angular momentum, via the simple equation:   $$n = 3 + 18L^2,$$ (18)and here L is a normalized version of the total encounter angular momentum $$\boldsymbol {L}$$. Equation (18) has been calibrated from numerical scattering experiments of single–binary interactions (Valtonen & Karttunen 2006). Now, to construct the schematic triangles, we can set   $$\frac{E}{180^{\circ }} = \frac{E_{\rm s}}{\theta _{\rm s}} = \frac{E_{\rm B}}{\theta _{\rm B}},$$ (19)where θs is the angle corresponding to the single star, which is found at the apex of the triangle, and θB = θB1 + θB2 are the angles for each binary component, at the base of the triangle. We assume θB1 = θB2 for this example, to keep things simple. If Es exceeds the critical energy for escape, which is the situation for the unbound case, then θs is negative. Hence, to correct for this, we set $$\theta _{\rm s}^{\prime } = -\theta _{\rm s}$$ along with $$\theta _{\rm B}^{\prime } = 360^{\circ } - \theta _{\rm B}$$. We emphasize that, unlike in the bound case in Sections 2.1 and 3.1, we allow for negative angles in this section for the unbound case, to treat unbound escaping stars appropriately (given that their total energy is positive, and we began by assigning positive angles to negative energies). An example of our schematic diagrams for the unbound case is shown in Fig. 3, which depicts five successive binary hardening events from single–binary scattering. All particles are assumed to have masses 1 M⊙. The initial binary has a semimajor axis of 5 au. We assume a velocity dispersion of 5 km s−1 for the surrounding stellar environment which provides the single stars for each scattering interaction. Hence, we take vinf = 5 km s−1 as the initial relative velocity at infinity for every single–binary interaction. We further assume that, in every final state, the escaping single star is ejected with a velocity equal to the peak of the velocity distribution in equation (16), or   $$v_{\rm s,peak} = \alpha \sqrt{|E|\frac{M - m_{\rm s}}{m_{\rm s}M}}.$$ (20)We set α = 0.5, which corresponds to typical intermediate values of the angular momentum L, and hence to isotropic scattering. We emphasize that this choice of adopting the peak velocity for the escaper [instead of drawing at random from the escaper velocity distribution in equation (16)] is arbitrary, and is adopted only for simplicity in this illustrative example. Figure 3. View largeDownload slide Schematic diagrams for the unbound case in energy-space. We show five successive binary hardening events due to single–binary scattering. The initial and final states for each scattering event are shown, respectively, at the top and bottom. The incoming and outgoing single stars always correspond to the vertex at the apex (i.e. furthest from the horizontal axis) of a given triangle. Left to right denotes the direction of increasing time. Figure 3. View largeDownload slide Schematic diagrams for the unbound case in energy-space. We show five successive binary hardening events due to single–binary scattering. The initial and final states for each scattering event are shown, respectively, at the top and bottom. The incoming and outgoing single stars always correspond to the vertex at the apex (i.e. furthest from the horizontal axis) of a given triangle. Left to right denotes the direction of increasing time. A similar formalism can potentially be applied in angular momentum-space, albeit this should be done with caution. In particular, the exact methodology for this depends on exactly which distributions have been derived, and ultimately whether or not they are joint or marginalized distributions. We have for the total angular momentum:   $$\boldsymbol {L} = \boldsymbol {L}_{\rm s} + \boldsymbol {L}_{\rm B}.$$ (21)From this, an analogous procedure can (in principle) be applied to define equivalent diagrams in angular momentum-space as done above in energy-space, for all three components of the angular momentum. We emphasize that the formalism presented here for characterizing the time evolution of a given binary star system as it experiences scattering interactions with single stars in a dense stellar environment is purely analytic. Consequently, any potential issues related to computational expense are entirely mitigated, whereas this is never really the case for either large-scale N-body simulations of star cluster evolution or suites of individual numerical scattering simulations. 3.3 Unbound case: many binaries Consider applying our formalism for the unbound case to an entire population of binary star systems. This is equivalent to dynamically evolving a population of binary stars in a live dense stellar environment. Using the analytic distribution functions from the previous section for the properties of the products of three-body scattering interactions, we can construct a purely analytic formalism to propagate an initial distribution of binary orbital parameters forward through time due to single–binary scattering interactions. This is done as follows. We begin with a population of N binaries, and assume a distribution of orbital separations f(a), a distribution of mass ratios f(q) and a primary mass distribution f(m1). These distribution functions can be convolved to obtain the distribution of binary orbital energies fB(EB), using the chain rule:   \begin{eqnarray} f_{\rm B}(E_{\rm B}) &=& \Big ( \frac{{\rm d}N}{{\rm d}a} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}a} \Big )^{-1} + \Big ( \frac{{\rm d}N}{{\rm d}q} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}q} \Big )^{-1} \nonumber\\ &&+\, \Big ( \frac{{\rm d}N}{{\rm d}m_{\rm 1}} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}m_{\rm 1}} \Big )^{-1}. \end{eqnarray} (22)The rate of change of fB(EB) is   $$\frac{{\rm d}f_{\rm B}(E_{\rm B})}{{\rm d}t} = \frac{{\rm d}f_{\rm B}}{{\rm d}E_{\rm B}}\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}N} \Big )\Big ( \frac{{\rm d}N}{{\rm d}t} \Big ) = \frac{{\rm d}f_{\rm B}}{{\rm d}E_{\rm B}}\Big ( \frac{\Gamma (E_{\rm B})}{f(|E_{\rm B}|)} \Big ),$$ (23)where Γ(EB) is the rate of single–binary interactions (Leigh & Sills 2011) and f(|EB|) is the distribution of binary orbital energies given in equation (17). The interaction rate Γ depends on the binary orbital energy via the orbital separation and total binary mass, in addition to the properties of the host stellar environment (density, velocity dispersion, etc.). The angular momentum dependence enters via the assumed impact parameter distribution f(b). For isotropic scattering in spherical star clusters, this typically takes the form f(b)db = bdb and is related to the total angular momentum via L = μvinfb, where μ is the reduced mass of the initial single–binary system. Ultimately, these relations can be used to relate the assumed impact parameter distribution to a distribution of power-law indices n in equation (17), via equation (18). The only requirement for the above analytic methodology to be accurate, is for the single–binary interactions to enter a long-lived ‘resonant’ state, in which all particles are thoroughly mixed in phase space and the assumption of ergodicity is upheld. Fortunately, the cross-section for such a resonant encounter is straight-forward to calculate analytically (Hut 1983b). 4 DISCUSSION In this paper, we present a schematic discretization of the general three-body problem in Newtonian gravity. Our diagrams fully and quantitatively describe the time evolution of chaotic three-body interactions involving point particles for many, if not most, encounters. Hence, using our diagrams, the time evolution of most (bound or unbound) interactions can in principle be reconstructed in position- and velocity-space. In this section, we address the shortcomings and utility of our ‘chaos diagrams’, discuss possible applications of our method and look to future work. 4.1 The E = 0 boundary Particles can end up asymptotically approaching the E = 0 boundary, although the probability for such an event is vanishingly small in the limit E → 0 for negative total encounter energies. The duration of these asymptotic excursions can be exceedingly long, and by far dominate the total computer run times for simulations. Nearly all of the computational expense lies in these prolonged excursions. Our schematic formalism reduces these prolonged excursion events to single discrete transformations in energy- and angular momentum-space. Thus, if numerical scattering simulations of three-body interactions were to be replaced with our chaos diagrams, this computational cost would be altogether mitigated. 4.2 How to perform gravitational scattering experiments without a computer For many, if not most, resonant three-body interactions, all of the quantitative information is encoded in our chaos diagrams. Thus, in principle, these diagrams could remove the need to perform computer simulations, at least for some regions of the total relevant parameter space. It is just a matter of generating the distributions of particle energies and angular momenta at each discrete excursion event, without relying on computer simulations. This amounts to somehow calculating distribution functions for these quantities, from which we can sample at each discrete excursion event. This has already been done semi-analytically for the unbound case, which has provided simple analytic functions to work with [see Section 3.2 and equation (17)]. No such analytic functions are known for the bound case. If such a function is derived, however, then in principle the formalism presented in this paper could completely replace numerical scattering simulations in those regions of parameter space where the assumption of ergodicity is upheld. For this to be the case, two additional criteria must also be satisfied: The approximation E ∼ Es + EB must remain valid throughout the entire (bound) interaction, such that each transformation in energy- and angular momentum-space corresponds to a single orbit of the (temporary) single star about the (temporary) binary. The state of each such temporary orbit is evaluated at apocentre. An analytic formalism for mapping the system from one region of phase space to another between successive transformations (i.e. temporary ejection events) during the bound state. For example, one way this could occur is if each transformation in energy- and angular momentum-space turns out to be independent of energy and angular momentum, and correspond to a random sampling of the corresponding distribution functions. That is, the different meta-stable states corresponding to each discrete event are independent, such that the probability of the system being in one region of parameter space does not affect the probability of it finding itself in another after a given transformation. The first criterion is not always upheld, and it can occur that all three particles have comparable energy and angular momenta, such that the approximation E ∼ Es + EB temporarily breaks down during the bound state. Thus, currently, the formalism presented in this paper cannot do this accurately, since we ignore cross-terms in the energy equation that become more important during very short excursions. However, as already stated, for many three-body interactions this assumption is strictly upheld for the entire duration of the encounter. In future work, we intend to study this interesting limit of the chaotic three-body problem in more detail. We intend to better define and quantify this transition, in order to understand how to include in our schematic formalism this more ambiguous phase in the evolution of the system. The second criterion can in principle be tested using numerical scattering experiments via the procedure described in Section 2. In an effort to quantify and illustrate this, we identified one especially prolonged single–binary interaction from our suite of simulations. This yielded a total of 120 discrete transformation events (i.e. triangles in Fig. 2), and provided the needed statistics to get an approximate idea of the distribution of binary orbital energies during the bound state of a given encounter, as shown in Fig. 4. As is clear, the simulated distribution deviates from what is expected from a smooth continuation of equation (17) from the bound state to the unbound state, and resembles more of a skewed Gaussian. Unfortunately, the reason for this is not clear, and we do not know how sensitive it is to our adopted criterion for sampling the binary orbital energies. As already stated, clearly defining an excursion event can become ambiguous during very short excursions. The data shown in Fig. 4 naively suggest that the system remains close to its initial location in phase space throughout the entire resonant phase of the interaction. More work will be needed to verify and better understand the basic results shown and discussed here. In the subsequent sections, we discuss two possible methods for better accomplishing the overarching goal of quantifying the distributions of binary orbital energies and angular momenta during the bound state. Figure 4. View largeDownload slide The distribution of binary orbital energies for the bound state, normalized by the total encounter energy as z = |E|/|EB|. We show the results of a single very prolonged three-body interaction, with 120 independent discrete transformations in energy- and angular momentum-space. We adopt a bin size of Δz = 0.05. The dotted (n = 3) and dashed (n = 3.5) lines show the binary orbital energy distribution function given by equation (17) for the unbound case in the low angular momentum limit, for two different values of the power-law index n. Figure 4. View largeDownload slide The distribution of binary orbital energies for the bound state, normalized by the total encounter energy as z = |E|/|EB|. We show the results of a single very prolonged three-body interaction, with 120 independent discrete transformations in energy- and angular momentum-space. We adopt a bin size of Δz = 0.05. The dotted (n = 3) and dashed (n = 3.5) lines show the binary orbital energy distribution function given by equation (17) for the unbound case in the low angular momentum limit, for two different values of the power-law index n. 4.2.1 Embracing the machine Motivated by our results in this paper and especially Fig. 4, one way to obtain distribution functions to describe the probabilities of each particle having a given combination of energy and angular momentum during a given excursion event is machine learning. That is, using a machine learning algorithm to calculate these distribution functions directly from millions (and preferably billions) of three-body simulations, while implementing supervised learning techniques. For example, an artificial neural network learning algorithm could be implemented, with the interconnections between artificial neurons defined by repeated fine-graining of different volume elements of the total phase space. This could be used to evaluate the degree to which these regions can (or cannot) be causally connected by a single discrete transformation, for each level or degree of fine-graining. In other words, the machine could be trained to learn how the different meta-stable states corresponding to each discrete event are inter-related, and to what extent (if any) the probability of the system being in one region of parameter space affects the probability of it finding itself in another after a given transformation. This offers an ideal method for testing the assumption of ergodicity during the bound state; that is, that the available phase space should be populated uniformly. Thus, this method will not only reveal the underlying patterns characteristic of the time evolution of the system throughout phase space, but also the underlying physical mechanism(s) driving this evolution. With the necessary distribution functions in hand, an operator could potentially be generated that propagates the system forward through time, transforming one triangle into the next. Consider the following transformation $$\hat{T}$$ applied to a three-dimensional vector $$\boldsymbol {\theta } = [\theta _{\rm 1},\theta _{\rm 2},\theta _{\rm 3}]$$:   \begin{equation*} \hat{T} {\left[\begin{array}{cccc}\theta _{\rm 1} \\ \theta _{\rm 2} \\ \theta _{\rm 3} \end{array}\right]} = {\left[\begin{array}{cccc}\theta ^{\prime }_{\rm 1} \\ \theta ^{\prime }_{\rm 2} \\ \theta ^{\prime }_{\rm 3} \end{array}\right]}. \end{equation*} The primes denote the vector components in the new transformed state. Using machine learning, an operator $$\hat{T}$$ could be generated for each of the four three-dimensional vectors relevant to our problem, namely E, Lx, Ly, and Lz. Each operator would be applied to the corresponding vector at each discrete excursion event, transforming it to a new state. Machine learning algorithms are an option potentially worth exploring more, since they offer one avenue towards achieving the overarching goal of eventually being able to generate our chaos diagrams without any intervention from computers. But the irony is rich. Below, we explore another promising option for accomplishing our goal, which could potentially take us all the way beyond the machine. 4.2.2 Beyond the machine In Section 3.2, we used analytic distribution functions from Valtonen & Karttunen (2006) for the products of single–binary interactions to calculate the energies of the particles in the final state, and presented chaos diagrams for the initial and final states using a purely analytic approach. As already discussed, to complete a purely analytic approach, this leaves only the intermediate bound meta-stable states, in which the ejected single star eventually returns to interact with the binary at least once more. If this distribution can be derived, we can reduce the creation of the chaos diagrams presented in Fig. 2 to a purely analytic problem. We emphasize that this could entirely replace the need for computer simulations, for interactions that satisfy the criteria discussed in the preceding sections. This can already be done in energy-space (and angular momentum-space), as shown in Section 3.2 for the unbound states. In a forthcoming paper (Stone & Leigh in preparation), we re-visit the original Monaghan (1976a) formalism for calculating the density of configuration states per unit energy. A more complete description of the outcome distributions can be derived by incorporating an angular momentum-dependence into the integrand in the expression for the density of states. The hope is that this will provide the required distribution functions for the angular momenta of the final binary and the ejected single star, in addition to the corresponding distribution functions for both energy and angular momentum for the bound case. These joint distributions can in turn be marginalized over to calculate specific distribution functions of interest. We note that, if the assumption of ergodicity is upheld collectively for both the bound and unbound states, then the mean number of excursions for a collection of three-body interactions should be equal to the total volume in phase space available to the system (i.e. bound and unbound states) divided by the total volume corresponding to bound states alone. This potentially provides a means of calculating an expectation value for the number of excursion events for a given (resonant) three-body interaction. The potential to bypass the need for computers in studying the chaotic three-body problem is tempting, to say the least. If the bound state distribution functions (i.e. energy, angular momentum, as well as the duration and number of excursions) can indeed be derived using a First Principles approach, then our entire formalism could become purely analytic. Numerical simulations would only be needed to better understand these key distribution functions, and verify their accuracy. This will be the focus of future work. Finally, we point out that, although the diagrams presented here were primarily used for illustrative purposes throughout this paper, they have the potential to fully depict the full time evolution of individual three-body interactions, along with all of the associated quantitative data needed to uniquely identify any one interaction and even reproduce it as a visualization of the time evolution of the three-body system in position-space (i.e. what is often done to generate a visual aid for these types of interactions). We would argue, however, that this information can be conveyed to an audience using our diagrams much more quickly and efficiently than upon using a time series visualization, such as a movie. The underlying formalism of our method also provides an extremely efficient means of storing the data needed to fully and quantitatively reproduce a given scattering interaction, compared to other output files in common use, with the positions outputted at regular time-steps. 5 SUMMARY In this paper, we present a formalism for constructing schematic diagrams to depict chaotic three-body interactions in Newtonian gravity. This is done by decomposing each interaction into a series of discrete transformations in energy- and angular momentum-space. With each transformation, the system changes state as the particles re-distribute their energy and angular momentum. These diagrams have the virtue of containing all of the quantitative information needed to fully characterize many (if not most) bound or unbound interactions through time and space, including the total duration of the interaction, the initial and final stable states in addition to every intervening temporary meta-stable state. As shown via an illustrative example for the bound case, prolonged excursions of one of the particles by far dominates the computational cost of the simulations. In our formalism, these are reduced to a single discrete transformation in energy- and angular momentum-space, thereby potentially mitigating any computational expense. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening. Finally, we provide a method for dynamically evolving entire populations of binaries via three-body scattering interactions, using a purely analytic formalism. In principle, the techniques presented here are adaptable to other three-body problems that conserves energy and angular momentum. Acknowledgements We thank an anonymous reviewer for very helpful comments and suggestions that improved the quality of our work. We further thank Johan Samsing and Nicholas Stone for an early read of our manuscript and critical feedback. NWCL acknowledges support from a Kalbfleisch Fellowship at the American Museum of Natural History and the Richard Gilder Graduate School, as well as support from National Science Foundation Award AST 11-09395. Footnotes 1 For the source code, see http://fewbody.sourceforge.net. 2 Although, in order to conserve energy and momentum, it is easiest to sample from only one of them. In our example case, we sample from the velocity distribution in equation (16), and use this to calculate a corresponding binary orbital energy using conservation of energy. REFERENCES Agekyan T. A., Anosova J. P., 1967, Astron. Zh. , 44, 1261 Anosova J. P., Orlov V. V., 1994, Celest. Mech. Dyn. Astron. , 59, 327 https://doi.org/10.1007/BF00692101 CrossRef Search ADS   Antonini F., Chatterjee S., Rodriguez C. L., Morscher M., Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ , 816, 65 https://doi.org/10.3847/0004-637X/816/2/65 CrossRef Search ADS   Davies M. B., Blackwell R., Bailey V. C., Sigurdsson S., 1998, MNRAS , 301, 745 https://doi.org/10.1111/j.1365-8711.1998.02027.x CrossRef Search ADS   Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A., 2004, MNRAS , 352, 1 https://doi.org/10.1111/j.1365-2966.2004.07914.x CrossRef Search ADS   Geller A. M., Leigh N. W. C., 2015, ApJ , 808, L25 https://doi.org/10.1088/2041-8205/808/1/L25 CrossRef Search ADS   Geller A. M., Hurley J. R., Mathieu R. D., 2013, AJ , 145, 8 https://doi.org/10.1088/0004-6256/145/1/8 CrossRef Search ADS   Heggie D. C., 1975, MNRAS , 173, 729 https://doi.org/10.1093/mnras/173.3.729 CrossRef Search ADS   Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS , 363, 293 https://doi.org/10.1111/j.1365-2966.2005.09448.x CrossRef Search ADS   Hut P., Bahcall J. N., 1983, ApJ , 268, 319 https://doi.org/10.1086/160956 CrossRef Search ADS   Leigh N., Geller A. M., 2012, MNRAS , 425, 2369 https://doi.org/10.1111/j.1365-2966.2012.21689.x CrossRef Search ADS   Leigh N., Geller A. M., 2013, MNRAS , 432, 2474 https://doi.org/10.1093/mnras/stt617 CrossRef Search ADS   Leigh N., Geller A. M., 2015, MNRAS , 450, 1724 https://doi.org/10.1093/mnras/stv685 CrossRef Search ADS   Leigh N., Sills A., 2011, MNRAS , 410, 2370 https://doi.org/10.1111/j.1365-2966.2010.17609.x CrossRef Search ADS   Leigh N., Giersz M., Webb J. J., Hypki A., De Marchi G., Kroupa P., Sills A., 2013, MNRAS , 436, 3399 https://doi.org/10.1093/mnras/stt1825 CrossRef Search ADS   Leigh N. W. C., Shara M. M., Geller A. M., 2016, MNRAS , 459, 1242 https://doi.org/10.1093/mnras/stw735 CrossRef Search ADS   Leonard P. J. T., 1989, AJ , 98, 217 https://doi.org/10.1086/115138 CrossRef Search ADS   Mapelli M., Zampieri L., 2014, ApJ , 794, 7 https://doi.org/10.1088/0004-637X/794/1/7 CrossRef Search ADS   Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei . Princeton Univ. Press, Princeton, NJ Miller M., Davies M. B., 2012, ApJ , 755, 81 https://doi.org/10.1088/0004-637X/755/1/81 CrossRef Search ADS   Monaghan J. J., 1976a, MNRAS , 176, 63 https://doi.org/10.1093/mnras/176.1.63 CrossRef Search ADS   Monaghan J. J., 1976b, MNRAS , 177, 583 https://doi.org/10.1093/mnras/177.3.583 CrossRef Search ADS   Newton I., 1760, Philosophiae Naturalis Principia Mathematica . Regalis Societatis Praesses, Trinity College Orlov V. V., Rubinov A. V., Shevchenko I. I., 2010, MNRAS , 408, 1623 https://doi.org/10.1111/j.1365-2966.2010.17239.x CrossRef Search ADS   Shevchenko I. I., 2010, Phys. Rev. E , 81, 6216 https://doi.org/10.1103/PhysRevE.81.066216 CrossRef Search ADS   Valtonen M. J., 1975, Mem. RAS , 80, 77 Valtonen M., Karttunen H., 2006, The Three-Body Problem . Cambridge Univ. Press, Cambridge Google Scholar CrossRef Search ADS   Valtonen M. J., Mylläri A. A., Orlov V. V., Rubinov A. V., 2003, Astron. Lett. , 29, 41 https://doi.org/10.1134/1.1537377 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Illustrating chaos: a schematic discretization of the general three-body problem in Newtonian gravity

Monthly Notices of the Royal Astronomical Society, Volume 476 (1) – May 1, 2018
8 pages

/lp/ou_press/illustrating-chaos-a-schematic-discretization-of-the-general-three-7hghTsXYZs
Publisher
Oxford University Press
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty192
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present a formalism for constructing schematic diagrams to depict chaotic three-body interactions in Newtonian gravity. This is done by decomposing each interaction into a series of discrete transformations in energy- and angular momentum-space. Each time a transformation is applied, the system changes state as the particles re-distribute their energy and angular momenta. These diagrams have the virtue of containing all of the quantitative information needed to fully characterize most bound or unbound interactions through time and space, including the total duration of the interaction, the initial and final stable states in addition to every intervening temporary meta-stable state. As shown via an illustrative example for the bound case, prolonged excursions of one of the particles, which by far dominates the computational cost of the simulations, are reduced to a single discrete transformation in energy- and angular momentum-space, thereby potentially mitigating any computational expense. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening. Finally, we provide a method for dynamically evolving entire populations of binaries via three-body scattering interactions, using a purely analytic formalism. In principle, the techniques presented here are adaptable to other three-body problems that conserve energy and angular momentum. gravitation, scattering, methods: analytical, binaries: close, stars: kinematics and dynamics, globular clusters: general 1 INTRODUCTION After hundreds of years of study (e.g. Newton 1686), the chaotic three-body problem remains pertinent to the field of stellar dynamics, as it relates to our understanding of interactions within globular (e.g. Heggie 1975; Leigh & Geller 2012; Leigh et al. 2013; Antonini et al. 2016), open (e.g. Leonard 1989; Hurley et al. 2005; Leigh & Sills 2011; Geller, Hurley & Mathieu 2013; Leigh & Geller 2013; Mapelli & Zampieri 2014; Geller & Leigh 2015; Leigh et al. 2016), and even nuclear (e.g. Davies et al. 1998; Merritt 2013) star clusters. The problem remains analytically unsolved. Historically, any solutions involved approximate methods or the application of constraining assumptions used to simplify the problem. The best known such example of this approach is the restricted three-body problem, as applied to the Earth–Sun–Moon system. But modern computer technologies have come a long way and, in so doing, sparked a new wave of intensive study of such systems over the last few decades. Most simulations for the chaotic three-body problem run very quickly on modern computers. Every now and then, however, you will encounter one simulation that gets hung up, and takes many orders of magnitude more computer time to run to completion. Thus, upon performing a suite of numerical experiments of single–binary scattering, for example, it is a handful of individual runs that can consume the vast majority of the total computational cost. For small tidal tolerance parameters (see Section 2.2 for more details), this is due to prolonged excursions, where one of the stars ends up on a near-zero energy orbit. In fact, some authors have argued that the mean disruption time should formally be infinite for resonant three-body interactions (e.g. Valtonen 1975; Valtonen & Karttunen 2006; Orlov, Rubinov & Shevchenko 2010; Shevchenko 2010; Leigh et al. 2016); if an infinite number of simulations could be performed, the probability of obtaining one or two simulations with near infinite excursion time events becomes finite. These individual events can be bottlenecks in studies of the chaotic three-body problem, and compound to become even worse when additional particles are included into the mix. Chaotic three-body systems interacting resonantly display a characteristic evolution in time that begins when one particle is sent on a temporary excursion while remaining bound to the system. Upon returning to periastron, the particle initiates another close triple encounter, interacts strongly with the other two particles, and the process repeats (e.g. Agekyan & Anosova 1967; Anosova & Orlov 1994). To refer to events where the ejected particle is sent on a bound excursion we will use the term ‘ejection’. Conversely, to refer to events where the particle becomes unbound from the system and does not return, we will use the term ‘escape’. The end result of this chaotic dance is always an escape event where one of the bodies is ejected with a sufficiently high velocity to become unbound, and escapes to spatial infinity. The chaotic progression or evolution described above can be simplified or reduced to a single discrete transformation in energy and angular momentum space. The initial and final states of the system are qualitatively the same, since the unbound system can be decomposed into a single star and a binary. Once the initial relative energies and angular momenta of the single star and the binary are decided (via the input parameters that define a given scattering experiment), a transformation is applied. This changes the relative energies and angular momenta such that they are distributed differently among the particles in the initial and final states. By neglecting the intermediary chaos that occurs in the bound state, statistical ensembles of outcomes can be considered directly, rather than relying on specific choices of scattering parameters. In this paper, we present a schematic formalism to discretize the general three-body problem in Newtonian gravity, in both the bound and unbound states. This is done in such a way that, for most cases of interest, our diagrams fully and quantitatively characterize the time evolution of a given chaotic single–binary interaction. We use the numerical scattering code fewbody to simulate a series of binary-single encounters involving identical point particles in Section 2. We choose one particularly interesting example, that features a very prolonged excursion of one of the stars, and generate our schematic diagrams to depict the time evolution of this interaction. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening, and even describe a purely analytic method for dynamically evolving entire populations of binaries via three-body scattering interactions. In Section 4, we discuss the utility of our schematic formalism, address possible shortcomings and applications of our method and explain how our formalism might eventually replace the need for computer simulations of single–binary scatterings. We summarize in Section 5. 2 METHOD In this section, we introduce our formalism for constructing schematic diagrams of the general three-body problem, adapted from Leigh & Geller (2012) and Leigh et al. (2016). We further present the numerical scattering experiments used to study the time evolution of the chaotic three-body problem, from which we construct an illustrative example of our schematic diagram formalism. 2.1 Constructing the schematic diagrams In this section, we begin by adapting the schematic diagrams first presented in Leigh & Geller (2012) to depict individual snapshots of the system in time for application to the general three-body problem, including arbitrarily long excursions of one of the particles, followed by a description of our method for discretizing the interaction. First, the total energy of the three-body system can be written as   $$E = \sum _{\rm i=1}^{\rm N=3} E_{\rm i}^{\prime } = \sum _{\rm i=1}^{\rm N=3} \Big ( \frac{1}{2}m_{\rm i}v_{\rm i}^2 - \frac{1}{2}\sum _{\rm j=1, j \ne i}^{N=3} \frac{Gm_{\rm i}m_{\rm j}}{r_{\rm ij}} \Big ),$$ (1)where mi is the mass of the ith particle, vi is its velocity relative to the system centre of mass, ri is the magnitude of its distance from the system centre of mass, and $$r_{\rm ij} = |\bar{r}_{\rm i} - \bar{r}_{\rm j}|$$. Re-writing equation (1), we get   $$E - E_{\rm 0} = \sum _{\rm i=1}^{\rm N=3} \Big (E_{\rm i}^{\prime } - E_{\rm 0}/3 \Big ) = \sum _{\rm i=1}^{\rm N=3} E_{\rm i},$$ (2)where E0 is a negative constant with units of energy. The addition of the constant E0 serves only to shift the zero-point of the total encounter energy. This is needed to avoid the appearance of negative angles in our formalism. We adopt this approach for the bound case in Section 3.1. In Section 3.2 we consider a series of unbound three-body scattering interactions, and explore a different procedure for handling negative angles. Here, negative angles instead correspond to unbound particles. We will return to this in Section 3.2. Next, we point out that the sum of the angles of any polygon always add to a constant. For example, the sum of the angles of a triangle sum to 180°. Similarly, the energy contained in each particle Ei must always add up to equal the (shifted) total encounter energy E − E0. Hence, we can write   $$\frac{E - E_{\rm 0}}{180^{\circ }} = \frac{E_{\rm i}}{\theta _{\rm i}},$$ (3)where θi denotes the angle of a triangle corresponding to particle i with total energy Ei, as given by equation (2). For each triangle, we fix the length of the side of the triangle oriented along the x-axis (i.e. parallel to the direction of increasing time) to be equal to unity. This ensures that, in energy-space, the symmetry of each triangle is positively correlated with the duration of the corresponding excursion event. More equilateral triangles correspond to longer excursion events (see Section 3 for more details). A similar exercise can be performed to generate triangles for each component of the total angular momentum, namely Lx, Ly, and Lz. Using the x-component of the total angular momentum as an example we can write   $$L_{\rm x} - L_{\rm x,0} = \sum _{\rm i=1}^{\rm N=3} \Big (L_{\rm x,i}^{\prime } - L_{\rm x,0}/3 \Big ) = \sum _{\rm i=1}^{\rm N=3} L_{\rm x,i},$$ (4)and hence:   $$\frac{L_{\rm x} - L_{\rm x,0}}{180^{\circ }} = \frac{L_{\rm x,i}}{\theta _{\rm x,i}},$$ (5)where θx, i denotes the angle of a triangle corresponding to particle i with an x-component of the total angular momentum Lx, i, as given by equation (4). Next, we describe our formalism for discretizing the time evolution of the interaction. As described in Section 1, the evolving three-body system can usually be decomposed into a temporary binary, with the other star acting as a temporary single star on an excursion. This single–binary pair forms a closed (bound) temporary orbit. In order to formally define this criterion, we output from the simulations described in the subsequent section the critical time-steps at which the radial component of the centre of mass velocity of the temporary binary changes sign from positive to negative. This always coincides with the time-step at which the radial component of the centre of mass velocity of the temporary single star also changes sign. Every time this transition occurs, we record the energies of the particles and all three components of their total angular momentum. We emphasize that, since there are three particles, in this scenario there will always be one particle that is furthest from the system centre of mass, and is closest to the E ∼ 0 boundary. Using the formalism described above for generating triangles that quantitatively describe the temporary system in energy- and angular momentum-space, we then construct four triangles for each time-step corresponding to this change in sign of the radial velocities, one for energy and three for the total angular momentum (one for each component). The angle corresponding to each vertex of a given triangle is directly proportional to the energy or angular momentum content of the indicated particle, as described above. The angle corresponding to the top vertex is defined to always correspond to the temporary single star. The triangles are arranged sequentially from left to right, to indicate the direction of increasing time. The number of triangles is one greater than the number of applied transformations, to include the initial state of the system. This number indicates the total number of temporary excursion events, and clearly reveals any very long-lived computationally expensive excursions (via the relative symmetry of the triangles; see above and Section 3). 2.2 Numerical scattering experiments We calculate the outcomes of a series of single–binary (1+2) encounters using the fewbody numerical scattering code.1 The code integrates the usual N-body equations in configuration- (i.e. position-) space in order to advance the system forward in time, using the eighth-order Runge–Kutta Prince–Dormand integration method with ninth-order error estimate and adaptive time-step. For more details about the fewbody code, we refer the reader to Fregeau et al. (2004). All objects are point particles with masses of 1 M⊙. The initial binary has ai = 1 au, and eccentricity ei = 0. We fix the impact parameter at b = 0 for all simulations. The angles defining the initial relative configurations of the binary orbital plane and phases relative to the velocity vector of the incoming single star are chosen at random. We perform 1000 numerical scattering experiments to find a suitable illustrative example, including a prolonged excursion of one of the particles. This is illustrated in Fig. 1. Figure 1. View largeDownload slide The time evolution of an example single–binary simulation in position-space, projected into two dimensions. The unit of distance is 1 au on both axes. Notice the prolonged excursion of one of the particles, which dominates the total integration time of the simulation. Figure 1. View largeDownload slide The time evolution of an example single–binary simulation in position-space, projected into two dimensions. The unit of distance is 1 au on both axes. Notice the prolonged excursion of one of the particles, which dominates the total integration time of the simulation. To decide when an encounter is complete, we adopt the same criteria as Fregeau et al. (2004). To first order, we define this as the time when the separately bound hierarchies that constitute the system are no longer interacting with each other, or evolving internally. Formally, the integration is terminated when the top-level hierarchies have positive relative velocity and the corresponding top-level N-body system has positive total energy. We set the tidal tolerance parameter to δ = 10−7 for all simulations, which has been previously shown to be sufficient for the virial ratios (defined as the ratio between the total system kinetic energy and the total potential energy) considered here (see Geller & Leigh 2015 for more details). 3 RESULTS In this section, we generate and present our schematic diagrams for the discretization of the general three-body problem. We begin with a single three-body interaction that remains bound and in a resonant state for several system crossing times, and apply our method to the illustrative example shown in Fig. 1. We then present a modified version of our formalism to treat sequences of unbound single–binary scattering interactions, as occur during binary hardening in dense stellar environments, and show how (in some cases) our schematic diagrams can be generated purely analytically without any numerical simulations. Using this basic approach, we develop an analytic formalism for dynamically evolving entire populations of binaries in dense stellar environments due to three-body scattering interactions. 3.1 Bound case Schematic diagrams for the single–binary interaction depicted in Fig. 1 are shown in Fig. 2. Four separate triangle series are shown, one for energy-space (top) and three for angular momentum (bottom) with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. For each triangle, we fix the length of the side of the triangle oriented along the x-axis (i.e. parallel to the direction of increasing time). Combined with the imposed offset in energy E0, this ensures that equilateral configurations of our triangles in energy-space are positively correlated with the duration of the corresponding excursion event. The initial and final states correspond to, respectively, the first and last triangles. Hence, the progression of time is from left to right. Figure 2. View largeDownload slide The time evolution of the example single–binary simulation shown in Fig. 1 in both energy- and angular momentum-space, as described in the text. Four separate triangle series are shown, one for energy (top) and three for angular momentum (bottom), with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. At the vertex of the apex of each triangle opposing the horizontal axis, we indicate which particle corresponds to the temporary single star, using the same choice of labelling as in Fig. 1. Left to right denotes the direction of increasing time. We set E0 = 6.6 × 1038 J, Lx0 = 4.0 × 1048 kg m2 s−1, Ly0 = 4.0 × 1048 kg m2 s−1, and Lz0 = 2.0 × 1048 kg m2 s−1 for the zero-points, or equivalently E0 = 10|E|, Lx0 = 2300Lx, Ly0 = 510Ly, and Lz0 = 1300Lz. Figure 2. View largeDownload slide The time evolution of the example single–binary simulation shown in Fig. 1 in both energy- and angular momentum-space, as described in the text. Four separate triangle series are shown, one for energy (top) and three for angular momentum (bottom), with one series for each of the three components. The Lx and Lz diagrams have been reflected about the horizontal dashed line. At the vertex of the apex of each triangle opposing the horizontal axis, we indicate which particle corresponds to the temporary single star, using the same choice of labelling as in Fig. 1. Left to right denotes the direction of increasing time. We set E0 = 6.6 × 1038 J, Lx0 = 4.0 × 1048 kg m2 s−1, Ly0 = 4.0 × 1048 kg m2 s−1, and Lz0 = 2.0 × 1048 kg m2 s−1 for the zero-points, or equivalently E0 = 10|E|, Lx0 = 2300Lx, Ly0 = 510Ly, and Lz0 = 1300Lz. The diagrams shown in Fig. 2 fully and quantitatively describe the time evolution of the example single–binary interaction depicted in Fig. 1. In principle, using these diagrams, the time evolution of the interaction in position- and velocity-space can be re-constructed. This is because the energy and angular momentum of each particle is quantitatively encoded in the diagrams at each discrete exchange event via equations (3) and (5) combined with equations (1), (2), and (4), and the duration of each corresponding excursion is quantitatively encoded by the shape of each triangle (see below). Specifically, these diagrams encode 12 equations with 18 unknowns. The remaining six equations are found from the symmetry properties of the initial, final, and temporary binaries (i.e. given the provided information in our diagrams, only the positions and velocities of one of the binary components are needed to compute those of the other). This is done as follows. Consider a binary with component masses m1 and m2, oriented in the xy-plane. We transform our coordinate system such that it is stationary in the centre-of-mass frame of reference of the binary. Then, if the (instantaneous) positions and velocities of the component with mass m1 are (x1, y1, z1) and (vx, 1, vy, 1, vz, 1), then the positions and velocities of the second component are   \begin{eqnarray} x_{\rm 2} = -\frac{m_{\rm 1}}{m_{\rm 2}}x_{\rm 1}, \end{eqnarray} (6)  \begin{eqnarray} y_{\rm 2} = -\frac{m_{\rm 1}}{m_{\rm 2}}y_{\rm 1}, \end{eqnarray} (7)  \begin{eqnarray} z_{\rm 2} = z_{\rm 1} = 0, \end{eqnarray} (8)  \begin{eqnarray} v_{\rm x,2} = -\frac{m_{\rm 1}}{m_{\rm 2}}v_{\rm x,1}, \end{eqnarray} (9)  \begin{eqnarray} v_{\rm y,2} = -\frac{m_{\rm 1}}{m_{\rm 2}}v_{\rm y,1}, \end{eqnarray} (10)  \begin{eqnarray} v_{\rm z,2} = v_{\rm z,1} = 0. \end{eqnarray} (11)This follows from the definition of the system centre of mass, which defines the equation:   $$m_{\rm 1}\boldsymbol {{R}}_{\rm 1} + m_{\rm 2}\boldsymbol {{R}}_{\rm 2} = 0,$$ (12)where $$\boldsymbol {{R}}_{\rm 1}$$ and $$\boldsymbol {{R}}_{\rm 2}$$ are the radius vectors of the bodies with respect to the binary centre of mass, and   $$|\boldsymbol {{R}}_{i}| = \sqrt{x_{\rm i}^2 + y_{\rm i}^2 + z_{\rm i}^2},$$ (13)and i = 1 or 2. Now, consider the very prolonged excursion event indicated in Fig. 1. This corresponds to the third triangle (from the left) in Fig. 2. This is very close to a completely equilateral triangle in energy-space, to within fractions of a degree, whereas triangles corresponding to shorter excursion events are much less symmetric in energy-space. This behaviour is also evident by the corresponding triangles in Lx-, Ly-, and Lz-space, in which one particle clearly dominates the angular momentum budget during the excursion. We define the following parameter:   $$\epsilon = \sum _{\rm i=1}^{3} \left|\theta _{\rm i} - \frac{\pi }{3}\right|.$$ (14)The longest excursion events correspond to ε ∼ 0, whereas shorter excursion events satisfy ε > 0. The reason for ε ∼ 0 corresponding to the longest excursions is that, here, the temporary single star is nearly unbound such that its energy is very close to zero. Meanwhile, the temporary binary is approximately isolated, such that it (approximately) obeys the virial theorem. Hence, twice the total kinetic energy is roughly equal to the absolute value of the total potential energy in the temporary binary. This translates into E΄ ∼ 0 in equation (1), and hence all θi are approximately equal to π/3 via equations (2) and (3). 3.2 Unbound case: one binary In this section, we briefly reproduce the analytic formalism for unbound scattering presented in Valtonen & Karttunen (2006), in particular the key formulae and distribution functions to be used and discussed in the subsequent sections. These are presented here via an example application of our schematic formalism to the unbound case. In dense stellar environments, as occur in the cores of star clusters, a single–binary star can undergo repeated interactions with many different single stars. This sequence of single–binary scatterings typically acts to harden the binary (i.e. reduce its orbital separation). Our schematic diagrams can be adapted to quantitatively describe this process. The required distribution functions are already available in analytic form for the initial and final unbound states (i.e. the first and last triangles) (Valtonen & Karttunen 2006). This is thanks to pioneering efforts by Monaghan (1976a) and Monaghan (1976b), who re-worked the expression for the density of configuration states per unit energy in order to derive the distributions of binary orbital energies, eccentricities, and escaper velocities in the low- and high-angular momentum limits. This was later parametrized to be applicable to any total angular momentum (Valtonen & Karttunen 2006). Thus, we can already generate the initial and final states from these distribution functions. One way to do this is as follows. We write for the total encounter energy of a given single–binary interaction:   $$E = E_{\rm s} + E_{\rm B},$$ (15)where EB is the binary orbital energy and $$E_{\rm s} = \frac{1}{2}m_{\rm s}v_{\rm s}^2$$ is the kinetic energy of the single star in the centre of mass frame of the binary. From equation 7.19 in Valtonen & Karttunen (2006), the distribution of escaper velocities is   $$f(v_{\rm s}){\rm d}v_{\rm s} = \frac{(n-1)|E|^{n-1}(m_{\rm s}M/m_{\rm B}))v_{\rm s}{\rm d}v_{\rm s}}{(|E| + \frac{1}{2}(m_{\rm s}M/m_{\rm B})v_{\rm s}^2)^{n}},$$ (16)where ms is the mass of the single star, mB is the mass of the binary, and M = ms + mB is the total system mass. From equation 7.26, we have for the final distribution of binary orbital energies:   $$f(|E_{\rm B}|){{\rm d}}|E_{\rm B}| = (n-1)|E|^{n-1}|E_{\rm B}|^{-n}{{\rm d}}|E_{\rm B}|.$$ (17)Equations (16) and (17) can be sampled from directly,2 and hence final states constructed from the initial state. The parameter n corrects each distribution for the total angular momentum, via the simple equation:   $$n = 3 + 18L^2,$$ (18)and here L is a normalized version of the total encounter angular momentum $$\boldsymbol {L}$$. Equation (18) has been calibrated from numerical scattering experiments of single–binary interactions (Valtonen & Karttunen 2006). Now, to construct the schematic triangles, we can set   $$\frac{E}{180^{\circ }} = \frac{E_{\rm s}}{\theta _{\rm s}} = \frac{E_{\rm B}}{\theta _{\rm B}},$$ (19)where θs is the angle corresponding to the single star, which is found at the apex of the triangle, and θB = θB1 + θB2 are the angles for each binary component, at the base of the triangle. We assume θB1 = θB2 for this example, to keep things simple. If Es exceeds the critical energy for escape, which is the situation for the unbound case, then θs is negative. Hence, to correct for this, we set $$\theta _{\rm s}^{\prime } = -\theta _{\rm s}$$ along with $$\theta _{\rm B}^{\prime } = 360^{\circ } - \theta _{\rm B}$$. We emphasize that, unlike in the bound case in Sections 2.1 and 3.1, we allow for negative angles in this section for the unbound case, to treat unbound escaping stars appropriately (given that their total energy is positive, and we began by assigning positive angles to negative energies). An example of our schematic diagrams for the unbound case is shown in Fig. 3, which depicts five successive binary hardening events from single–binary scattering. All particles are assumed to have masses 1 M⊙. The initial binary has a semimajor axis of 5 au. We assume a velocity dispersion of 5 km s−1 for the surrounding stellar environment which provides the single stars for each scattering interaction. Hence, we take vinf = 5 km s−1 as the initial relative velocity at infinity for every single–binary interaction. We further assume that, in every final state, the escaping single star is ejected with a velocity equal to the peak of the velocity distribution in equation (16), or   $$v_{\rm s,peak} = \alpha \sqrt{|E|\frac{M - m_{\rm s}}{m_{\rm s}M}}.$$ (20)We set α = 0.5, which corresponds to typical intermediate values of the angular momentum L, and hence to isotropic scattering. We emphasize that this choice of adopting the peak velocity for the escaper [instead of drawing at random from the escaper velocity distribution in equation (16)] is arbitrary, and is adopted only for simplicity in this illustrative example. Figure 3. View largeDownload slide Schematic diagrams for the unbound case in energy-space. We show five successive binary hardening events due to single–binary scattering. The initial and final states for each scattering event are shown, respectively, at the top and bottom. The incoming and outgoing single stars always correspond to the vertex at the apex (i.e. furthest from the horizontal axis) of a given triangle. Left to right denotes the direction of increasing time. Figure 3. View largeDownload slide Schematic diagrams for the unbound case in energy-space. We show five successive binary hardening events due to single–binary scattering. The initial and final states for each scattering event are shown, respectively, at the top and bottom. The incoming and outgoing single stars always correspond to the vertex at the apex (i.e. furthest from the horizontal axis) of a given triangle. Left to right denotes the direction of increasing time. A similar formalism can potentially be applied in angular momentum-space, albeit this should be done with caution. In particular, the exact methodology for this depends on exactly which distributions have been derived, and ultimately whether or not they are joint or marginalized distributions. We have for the total angular momentum:   $$\boldsymbol {L} = \boldsymbol {L}_{\rm s} + \boldsymbol {L}_{\rm B}.$$ (21)From this, an analogous procedure can (in principle) be applied to define equivalent diagrams in angular momentum-space as done above in energy-space, for all three components of the angular momentum. We emphasize that the formalism presented here for characterizing the time evolution of a given binary star system as it experiences scattering interactions with single stars in a dense stellar environment is purely analytic. Consequently, any potential issues related to computational expense are entirely mitigated, whereas this is never really the case for either large-scale N-body simulations of star cluster evolution or suites of individual numerical scattering simulations. 3.3 Unbound case: many binaries Consider applying our formalism for the unbound case to an entire population of binary star systems. This is equivalent to dynamically evolving a population of binary stars in a live dense stellar environment. Using the analytic distribution functions from the previous section for the properties of the products of three-body scattering interactions, we can construct a purely analytic formalism to propagate an initial distribution of binary orbital parameters forward through time due to single–binary scattering interactions. This is done as follows. We begin with a population of N binaries, and assume a distribution of orbital separations f(a), a distribution of mass ratios f(q) and a primary mass distribution f(m1). These distribution functions can be convolved to obtain the distribution of binary orbital energies fB(EB), using the chain rule:   \begin{eqnarray} f_{\rm B}(E_{\rm B}) &=& \Big ( \frac{{\rm d}N}{{\rm d}a} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}a} \Big )^{-1} + \Big ( \frac{{\rm d}N}{{\rm d}q} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}q} \Big )^{-1} \nonumber\\ &&+\, \Big ( \frac{{\rm d}N}{{\rm d}m_{\rm 1}} \Big )\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}m_{\rm 1}} \Big )^{-1}. \end{eqnarray} (22)The rate of change of fB(EB) is   $$\frac{{\rm d}f_{\rm B}(E_{\rm B})}{{\rm d}t} = \frac{{\rm d}f_{\rm B}}{{\rm d}E_{\rm B}}\Big ( \frac{{\rm d}E_{\rm B}}{{\rm d}N} \Big )\Big ( \frac{{\rm d}N}{{\rm d}t} \Big ) = \frac{{\rm d}f_{\rm B}}{{\rm d}E_{\rm B}}\Big ( \frac{\Gamma (E_{\rm B})}{f(|E_{\rm B}|)} \Big ),$$ (23)where Γ(EB) is the rate of single–binary interactions (Leigh & Sills 2011) and f(|EB|) is the distribution of binary orbital energies given in equation (17). The interaction rate Γ depends on the binary orbital energy via the orbital separation and total binary mass, in addition to the properties of the host stellar environment (density, velocity dispersion, etc.). The angular momentum dependence enters via the assumed impact parameter distribution f(b). For isotropic scattering in spherical star clusters, this typically takes the form f(b)db = bdb and is related to the total angular momentum via L = μvinfb, where μ is the reduced mass of the initial single–binary system. Ultimately, these relations can be used to relate the assumed impact parameter distribution to a distribution of power-law indices n in equation (17), via equation (18). The only requirement for the above analytic methodology to be accurate, is for the single–binary interactions to enter a long-lived ‘resonant’ state, in which all particles are thoroughly mixed in phase space and the assumption of ergodicity is upheld. Fortunately, the cross-section for such a resonant encounter is straight-forward to calculate analytically (Hut 1983b). 4 DISCUSSION In this paper, we present a schematic discretization of the general three-body problem in Newtonian gravity. Our diagrams fully and quantitatively describe the time evolution of chaotic three-body interactions involving point particles for many, if not most, encounters. Hence, using our diagrams, the time evolution of most (bound or unbound) interactions can in principle be reconstructed in position- and velocity-space. In this section, we address the shortcomings and utility of our ‘chaos diagrams’, discuss possible applications of our method and look to future work. 4.1 The E = 0 boundary Particles can end up asymptotically approaching the E = 0 boundary, although the probability for such an event is vanishingly small in the limit E → 0 for negative total encounter energies. The duration of these asymptotic excursions can be exceedingly long, and by far dominate the total computer run times for simulations. Nearly all of the computational expense lies in these prolonged excursions. Our schematic formalism reduces these prolonged excursion events to single discrete transformations in energy- and angular momentum-space. Thus, if numerical scattering simulations of three-body interactions were to be replaced with our chaos diagrams, this computational cost would be altogether mitigated. 4.2 How to perform gravitational scattering experiments without a computer For many, if not most, resonant three-body interactions, all of the quantitative information is encoded in our chaos diagrams. Thus, in principle, these diagrams could remove the need to perform computer simulations, at least for some regions of the total relevant parameter space. It is just a matter of generating the distributions of particle energies and angular momenta at each discrete excursion event, without relying on computer simulations. This amounts to somehow calculating distribution functions for these quantities, from which we can sample at each discrete excursion event. This has already been done semi-analytically for the unbound case, which has provided simple analytic functions to work with [see Section 3.2 and equation (17)]. No such analytic functions are known for the bound case. If such a function is derived, however, then in principle the formalism presented in this paper could completely replace numerical scattering simulations in those regions of parameter space where the assumption of ergodicity is upheld. For this to be the case, two additional criteria must also be satisfied: The approximation E ∼ Es + EB must remain valid throughout the entire (bound) interaction, such that each transformation in energy- and angular momentum-space corresponds to a single orbit of the (temporary) single star about the (temporary) binary. The state of each such temporary orbit is evaluated at apocentre. An analytic formalism for mapping the system from one region of phase space to another between successive transformations (i.e. temporary ejection events) during the bound state. For example, one way this could occur is if each transformation in energy- and angular momentum-space turns out to be independent of energy and angular momentum, and correspond to a random sampling of the corresponding distribution functions. That is, the different meta-stable states corresponding to each discrete event are independent, such that the probability of the system being in one region of parameter space does not affect the probability of it finding itself in another after a given transformation. The first criterion is not always upheld, and it can occur that all three particles have comparable energy and angular momenta, such that the approximation E ∼ Es + EB temporarily breaks down during the bound state. Thus, currently, the formalism presented in this paper cannot do this accurately, since we ignore cross-terms in the energy equation that become more important during very short excursions. However, as already stated, for many three-body interactions this assumption is strictly upheld for the entire duration of the encounter. In future work, we intend to study this interesting limit of the chaotic three-body problem in more detail. We intend to better define and quantify this transition, in order to understand how to include in our schematic formalism this more ambiguous phase in the evolution of the system. The second criterion can in principle be tested using numerical scattering experiments via the procedure described in Section 2. In an effort to quantify and illustrate this, we identified one especially prolonged single–binary interaction from our suite of simulations. This yielded a total of 120 discrete transformation events (i.e. triangles in Fig. 2), and provided the needed statistics to get an approximate idea of the distribution of binary orbital energies during the bound state of a given encounter, as shown in Fig. 4. As is clear, the simulated distribution deviates from what is expected from a smooth continuation of equation (17) from the bound state to the unbound state, and resembles more of a skewed Gaussian. Unfortunately, the reason for this is not clear, and we do not know how sensitive it is to our adopted criterion for sampling the binary orbital energies. As already stated, clearly defining an excursion event can become ambiguous during very short excursions. The data shown in Fig. 4 naively suggest that the system remains close to its initial location in phase space throughout the entire resonant phase of the interaction. More work will be needed to verify and better understand the basic results shown and discussed here. In the subsequent sections, we discuss two possible methods for better accomplishing the overarching goal of quantifying the distributions of binary orbital energies and angular momenta during the bound state. Figure 4. View largeDownload slide The distribution of binary orbital energies for the bound state, normalized by the total encounter energy as z = |E|/|EB|. We show the results of a single very prolonged three-body interaction, with 120 independent discrete transformations in energy- and angular momentum-space. We adopt a bin size of Δz = 0.05. The dotted (n = 3) and dashed (n = 3.5) lines show the binary orbital energy distribution function given by equation (17) for the unbound case in the low angular momentum limit, for two different values of the power-law index n. Figure 4. View largeDownload slide The distribution of binary orbital energies for the bound state, normalized by the total encounter energy as z = |E|/|EB|. We show the results of a single very prolonged three-body interaction, with 120 independent discrete transformations in energy- and angular momentum-space. We adopt a bin size of Δz = 0.05. The dotted (n = 3) and dashed (n = 3.5) lines show the binary orbital energy distribution function given by equation (17) for the unbound case in the low angular momentum limit, for two different values of the power-law index n. 4.2.1 Embracing the machine Motivated by our results in this paper and especially Fig. 4, one way to obtain distribution functions to describe the probabilities of each particle having a given combination of energy and angular momentum during a given excursion event is machine learning. That is, using a machine learning algorithm to calculate these distribution functions directly from millions (and preferably billions) of three-body simulations, while implementing supervised learning techniques. For example, an artificial neural network learning algorithm could be implemented, with the interconnections between artificial neurons defined by repeated fine-graining of different volume elements of the total phase space. This could be used to evaluate the degree to which these regions can (or cannot) be causally connected by a single discrete transformation, for each level or degree of fine-graining. In other words, the machine could be trained to learn how the different meta-stable states corresponding to each discrete event are inter-related, and to what extent (if any) the probability of the system being in one region of parameter space affects the probability of it finding itself in another after a given transformation. This offers an ideal method for testing the assumption of ergodicity during the bound state; that is, that the available phase space should be populated uniformly. Thus, this method will not only reveal the underlying patterns characteristic of the time evolution of the system throughout phase space, but also the underlying physical mechanism(s) driving this evolution. With the necessary distribution functions in hand, an operator could potentially be generated that propagates the system forward through time, transforming one triangle into the next. Consider the following transformation $$\hat{T}$$ applied to a three-dimensional vector $$\boldsymbol {\theta } = [\theta _{\rm 1},\theta _{\rm 2},\theta _{\rm 3}]$$:   \begin{equation*} \hat{T} {\left[\begin{array}{cccc}\theta _{\rm 1} \\ \theta _{\rm 2} \\ \theta _{\rm 3} \end{array}\right]} = {\left[\begin{array}{cccc}\theta ^{\prime }_{\rm 1} \\ \theta ^{\prime }_{\rm 2} \\ \theta ^{\prime }_{\rm 3} \end{array}\right]}. \end{equation*} The primes denote the vector components in the new transformed state. Using machine learning, an operator $$\hat{T}$$ could be generated for each of the four three-dimensional vectors relevant to our problem, namely E, Lx, Ly, and Lz. Each operator would be applied to the corresponding vector at each discrete excursion event, transforming it to a new state. Machine learning algorithms are an option potentially worth exploring more, since they offer one avenue towards achieving the overarching goal of eventually being able to generate our chaos diagrams without any intervention from computers. But the irony is rich. Below, we explore another promising option for accomplishing our goal, which could potentially take us all the way beyond the machine. 4.2.2 Beyond the machine In Section 3.2, we used analytic distribution functions from Valtonen & Karttunen (2006) for the products of single–binary interactions to calculate the energies of the particles in the final state, and presented chaos diagrams for the initial and final states using a purely analytic approach. As already discussed, to complete a purely analytic approach, this leaves only the intermediate bound meta-stable states, in which the ejected single star eventually returns to interact with the binary at least once more. If this distribution can be derived, we can reduce the creation of the chaos diagrams presented in Fig. 2 to a purely analytic problem. We emphasize that this could entirely replace the need for computer simulations, for interactions that satisfy the criteria discussed in the preceding sections. This can already be done in energy-space (and angular momentum-space), as shown in Section 3.2 for the unbound states. In a forthcoming paper (Stone & Leigh in preparation), we re-visit the original Monaghan (1976a) formalism for calculating the density of configuration states per unit energy. A more complete description of the outcome distributions can be derived by incorporating an angular momentum-dependence into the integrand in the expression for the density of states. The hope is that this will provide the required distribution functions for the angular momenta of the final binary and the ejected single star, in addition to the corresponding distribution functions for both energy and angular momentum for the bound case. These joint distributions can in turn be marginalized over to calculate specific distribution functions of interest. We note that, if the assumption of ergodicity is upheld collectively for both the bound and unbound states, then the mean number of excursions for a collection of three-body interactions should be equal to the total volume in phase space available to the system (i.e. bound and unbound states) divided by the total volume corresponding to bound states alone. This potentially provides a means of calculating an expectation value for the number of excursion events for a given (resonant) three-body interaction. The potential to bypass the need for computers in studying the chaotic three-body problem is tempting, to say the least. If the bound state distribution functions (i.e. energy, angular momentum, as well as the duration and number of excursions) can indeed be derived using a First Principles approach, then our entire formalism could become purely analytic. Numerical simulations would only be needed to better understand these key distribution functions, and verify their accuracy. This will be the focus of future work. Finally, we point out that, although the diagrams presented here were primarily used for illustrative purposes throughout this paper, they have the potential to fully depict the full time evolution of individual three-body interactions, along with all of the associated quantitative data needed to uniquely identify any one interaction and even reproduce it as a visualization of the time evolution of the three-body system in position-space (i.e. what is often done to generate a visual aid for these types of interactions). We would argue, however, that this information can be conveyed to an audience using our diagrams much more quickly and efficiently than upon using a time series visualization, such as a movie. The underlying formalism of our method also provides an extremely efficient means of storing the data needed to fully and quantitatively reproduce a given scattering interaction, compared to other output files in common use, with the positions outputted at regular time-steps. 5 SUMMARY In this paper, we present a formalism for constructing schematic diagrams to depict chaotic three-body interactions in Newtonian gravity. This is done by decomposing each interaction into a series of discrete transformations in energy- and angular momentum-space. With each transformation, the system changes state as the particles re-distribute their energy and angular momentum. These diagrams have the virtue of containing all of the quantitative information needed to fully characterize many (if not most) bound or unbound interactions through time and space, including the total duration of the interaction, the initial and final stable states in addition to every intervening temporary meta-stable state. As shown via an illustrative example for the bound case, prolonged excursions of one of the particles by far dominates the computational cost of the simulations. In our formalism, these are reduced to a single discrete transformation in energy- and angular momentum-space, thereby potentially mitigating any computational expense. We further generalize our formalism to sequences of (unbound) three-body interactions, as occur in dense stellar environments during binary hardening. Finally, we provide a method for dynamically evolving entire populations of binaries via three-body scattering interactions, using a purely analytic formalism. In principle, the techniques presented here are adaptable to other three-body problems that conserves energy and angular momentum. Acknowledgements We thank an anonymous reviewer for very helpful comments and suggestions that improved the quality of our work. We further thank Johan Samsing and Nicholas Stone for an early read of our manuscript and critical feedback. NWCL acknowledges support from a Kalbfleisch Fellowship at the American Museum of Natural History and the Richard Gilder Graduate School, as well as support from National Science Foundation Award AST 11-09395. Footnotes 1 For the source code, see http://fewbody.sourceforge.net. 2 Although, in order to conserve energy and momentum, it is easiest to sample from only one of them. In our example case, we sample from the velocity distribution in equation (16), and use this to calculate a corresponding binary orbital energy using conservation of energy. REFERENCES Agekyan T. A., Anosova J. P., 1967, Astron. Zh. , 44, 1261 Anosova J. P., Orlov V. V., 1994, Celest. Mech. Dyn. Astron. , 59, 327 https://doi.org/10.1007/BF00692101 CrossRef Search ADS   Antonini F., Chatterjee S., Rodriguez C. L., Morscher M., Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ , 816, 65 https://doi.org/10.3847/0004-637X/816/2/65 CrossRef Search ADS   Davies M. B., Blackwell R., Bailey V. C., Sigurdsson S., 1998, MNRAS , 301, 745 https://doi.org/10.1111/j.1365-8711.1998.02027.x CrossRef Search ADS   Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A., 2004, MNRAS , 352, 1 https://doi.org/10.1111/j.1365-2966.2004.07914.x CrossRef Search ADS   Geller A. M., Leigh N. W. C., 2015, ApJ , 808, L25 https://doi.org/10.1088/2041-8205/808/1/L25 CrossRef Search ADS   Geller A. M., Hurley J. R., Mathieu R. D., 2013, AJ , 145, 8 https://doi.org/10.1088/0004-6256/145/1/8 CrossRef Search ADS   Heggie D. C., 1975, MNRAS , 173, 729 https://doi.org/10.1093/mnras/173.3.729 CrossRef Search ADS   Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS , 363, 293 https://doi.org/10.1111/j.1365-2966.2005.09448.x CrossRef Search ADS   Hut P., Bahcall J. N., 1983, ApJ , 268, 319 https://doi.org/10.1086/160956 CrossRef Search ADS   Leigh N., Geller A. M., 2012, MNRAS , 425, 2369 https://doi.org/10.1111/j.1365-2966.2012.21689.x CrossRef Search ADS   Leigh N., Geller A. M., 2013, MNRAS , 432, 2474 https://doi.org/10.1093/mnras/stt617 CrossRef Search ADS   Leigh N., Geller A. M., 2015, MNRAS , 450, 1724 https://doi.org/10.1093/mnras/stv685 CrossRef Search ADS   Leigh N., Sills A., 2011, MNRAS , 410, 2370 https://doi.org/10.1111/j.1365-2966.2010.17609.x CrossRef Search ADS   Leigh N., Giersz M., Webb J. J., Hypki A., De Marchi G., Kroupa P., Sills A., 2013, MNRAS , 436, 3399 https://doi.org/10.1093/mnras/stt1825 CrossRef Search ADS   Leigh N. W. C., Shara M. M., Geller A. M., 2016, MNRAS , 459, 1242 https://doi.org/10.1093/mnras/stw735 CrossRef Search ADS   Leonard P. J. T., 1989, AJ , 98, 217 https://doi.org/10.1086/115138 CrossRef Search ADS   Mapelli M., Zampieri L., 2014, ApJ , 794, 7 https://doi.org/10.1088/0004-637X/794/1/7 CrossRef Search ADS   Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei . Princeton Univ. Press, Princeton, NJ Miller M., Davies M. B., 2012, ApJ , 755, 81 https://doi.org/10.1088/0004-637X/755/1/81 CrossRef Search ADS   Monaghan J. J., 1976a, MNRAS , 176, 63 https://doi.org/10.1093/mnras/176.1.63 CrossRef Search ADS   Monaghan J. J., 1976b, MNRAS , 177, 583 https://doi.org/10.1093/mnras/177.3.583 CrossRef Search ADS   Newton I., 1760, Philosophiae Naturalis Principia Mathematica . Regalis Societatis Praesses, Trinity College Orlov V. V., Rubinov A. V., Shevchenko I. I., 2010, MNRAS , 408, 1623 https://doi.org/10.1111/j.1365-2966.2010.17239.x CrossRef Search ADS   Shevchenko I. I., 2010, Phys. Rev. E , 81, 6216 https://doi.org/10.1103/PhysRevE.81.066216 CrossRef Search ADS   Valtonen M. J., 1975, Mem. RAS , 80, 77 Valtonen M., Karttunen H., 2006, The Three-Body Problem . Cambridge Univ. Press, Cambridge Google Scholar CrossRef Search ADS   Valtonen M. J., Mylläri A. A., Orlov V. V., Rubinov A. V., 2003, Astron. Lett. , 29, 41 https://doi.org/10.1134/1.1537377 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create folders to

Export folders, citations