Bounds on energy and angular momentum loss in the full n-body problem

A new bound on the amended potential in an n-body system is derived and applied to the partitioning of energy and angular momentum in a disrupting gravitational aggregate. This result provides analytic insight into how energy and angular momentum can be lost or partitioned between different collections of bodies as they escape from each other. To better understand the possible outcomes in such a situation, some specific numerical tests are also performed for systems of N=3,4,5,6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 3, 4, 5, 6$$\end{document}. The results confirm that disrupting systems always lose energy, with a characteristic stochastic distribution of the lost energy. We also find that the system loses most of its angular momentum, although individual escaping components can retain significant nonzero angular momentum vectors that can be uniformly distributed in space as the number of bodies in the system increases.


Introduction and problem statement
The study of escape trajectories in the n-body problem has a long history and many significant results and discoveries. While there is a deep general understanding of the behavior of solutions asymptotic in time (Saari 1976), especially when we remove collisions between the point mass bodies, there are relatively few sharp results delineating necessary or sufficient conditions for the disruption of the n-body problem. The approach to understanding and constraining such final motions generally involves specifying the total energy and angular momentum of a system and then working out the limits for these. For example, for a system with nonzero angular momentum and positive energy, it can be shown that this is a sufficient condition for at least one component to escape. Going beyond this, if the system satisfies some specific conditions on relative placement and relative velocities, sufficient conditions can also be found for when additional bodies must escape (Patnaik 1975). For these disrupting systems, it has also been shown that the escaping sub-components may asymptotically become distinct n-body problems with their own conserved angular momentum and energy (Marchal and Saari 1976). Necessary conditions for the escape of system components seem to be more rare for the point-mass problem. The main results here exist for the three-body problem, where conditions can be found under which two of the bodies are guaranteed to be trapped in a bound orbit, in the presence of a third body (Marchal and Saari 1975).
For motion in the traditional point-mass n-body problem, the singular nature of the bodies leads to mathematical difficulties and physically unreal conditions, such as two-body collisions (which can be addressed) or three-or more-body collisions (which remain a fundamental difficulty in the problem). When the n-body problem is reformulated such that finite densities and rigid bodies are introduced, these singularities are removed and replaced by fundamental minimum distances between the bodies (Scheeres 2002). This allows for significant changes to the n-body problem, such as the existence for minimum energy configurations for any value of n, a result that does not exist in the point-mass problem for n ≥ 3 (Moeckel 1990). It also introduces new necessary conditions for the escape of components (Scheeres 2020), which again do not exist in the point-mass n-body problem. In the current paper, some new results that are relevant for these necessary escape conditions for the finite density (or full) n-body problem are derived. While these new results also apply to the point-mass n-body problem, the main motivation and focus of the paper are on the finite density problem.
Previous research has established that a key function for the dynamical analysis of a collection of rigid mass distributions interacting through gravity and possible contact forces can be captured by the amended potential, also called the minimum energy function (Scheeres 2012): where H is the total angular momentum of the system, E is the total energy, I H is the total (instantaneous) moment of inertia of the system about the angular momentum vector, and U is the instantaneous gravitational potential energy of the system. We note that both the moment of inertia and the gravitational potential are only functions of the system configuration space consisting of the relative positions of all of the bodies and their relative orientations with respect to each other. Given the configuration of a set of known mass distributions and their total energy and angular momentum, analysis of the amended potential enables one to find the relative equilibria of the system and their stability, a rigorous lower bound on the system energy at any time, conditions under which components of the system can escape due to gravitational interactions, and conditions for when they are bound to each other (see Scheeres 2002Scheeres , 2016aScheeres , b, 2020. In a previous paper sharp necessary conditions for the escape of collections of bodies were derived using this amended potential (Scheeres 2020). However, when using the amended potential to constrain whether a system can have components escape from each other, one runs into a practical issue once such an escape has occurred. As one body of the system asymptotically departs the remaining components, no matter how small it is, we generically find that I H → ∞. Thus, once escape of any mass occurs the amended potential function for the system collapses to where U ∞ consists of the gravitational potential energy of all bodies at a finite distance to each other, but whose terms that mutually escape all go to zero. Thus, once mutual escape of any collection of bodies occurs the amended potential loses all of its information encoding the system total angular momentum. In this paper we will prove that for an n-body system with m groups of component collections the amended potential can be bounded by a disaggregated version of the amended potential that takes the form: where E is as defined above and E i and E i define a disaggregated version of the amended potential and energy. Motivated by this result we discuss several implications for N -body systems and provide some numerical examples that shed light on the disaggregation process.
The significance of this result, as will be shown, is that an individual or multiple clusters can escape with the separate amended potentials still being well defined, and that the total energy remaining in the collection of bodies can be precisely tracked. Further, through numerical examples we provide initial insight into the range of energy and angular momentum that can be lost by systems as they disaggregate.
The paper is laid out as follows. In Sect. 2 we define the overall model and in Sect. 3 we introduce and apply the Jacobi transformation to the system. Then in Sect. 4 we present and prove our main theorem on the disaggregation of the amended potential and demonstrate the sharpness of the constraint. Following this in Sect. 5, we explore some implications of this theorem for systems that suffer a mutual escape of components. These results raise some simple questions about how an escaping component of the n-body problem affects the total energy and angular momentum of the remaining components in the system. To establish a better empirical understanding of these effects, we study a series of unstable few-body pointmass problems as they lose energy and angular momentum due to component ejections in Sect. 6. Finally, in Sect. 7 we summarize our main findings.

Model
We consider the finite-density, multiple-body problem where it is assumed that all of the component bodies are rigid and can thus rest on each other to form new varieties of equilibrium (Scheeres 2002(Scheeres , 2016a. We consider an arbitrary number of such bodies and borrow some key definitions from Scheeres (2020) to properly pose the system. We also choose a numbering system that leads to a more convenient notation for the main results.
We are given N + 1 finite density bodies P i , each of which has a mass distribution defined as B i , where i = 0, 1, 2, . . . , N . Each body has a mass m i , a rigid-body inertia dyadicĪ R i , a position r i , an orientation dyadicT i , a velocity v i , and an angular velocity i , all generally specified with respect to an inertial frame, except that the orientation dyadic will orient the ith body frame in inertial space and the ith rigid-body inertia dyadic is specified in the ith body-fixed frame. Each body will contribute its energy and angular momentum to the system, allowing us to define the total angular momentum, energy and inertia dyadic as where U i j are the mutual gravitational potentials between two rigid bodies and are a function of the relative position and orientation between the bodies, r i j = r j − r i andT i j =T T j ·T i . The term U ii is the self-potential of a given body,Ū is the unity dyadic, the product rr is a dyad, and the rigid-body angular momentum and rotational inertia are mapped from a body-fixed frame into the inertial frame with the orientation dyadicT i .
One important implication of the finite density rigid-body assumption is that any two bodies have a constraint on how close their centers of mass can come to each other, denoted as r i j ≥ d(r i j , T i j ), where equality occurs when the two bodies touch. This form of the constraint implicitly implies that the bodies are mutually convex, although this assumption can be relaxed. In the simplest model, the individual bodies are spheres with a given size and density, however the problem is easily generalized to arbitrarily shaped components that exert gravitational torques on each other and whose rotational dynamics must also be accounted for Scheeres (2002).
A single collection of bodies can be uniquely specified by a set of integers that specify the individual bodies. For example, I = {i 1 , i 2 , i 3 , . . . , i n } is a list of the bodies P i where i ∈ I. Then this collection of bodies is specified as P(I) = {P i |i ∈ I}.
Given this notation, we can define the global partition of our N + 1-body problem into The M subscript will be used when we wish to emphasize the number of components, but will be suppressed in other situations. Such a partition P(I) is a global partition if its individual elements P(I i ) are ordered, non-empty, disjoint and if their union consists of all the mass in the system. This is discussed in more generality in Scheeres (2020).

Jacobi coordinates
A crucial component to establishing our result is to transform the positions and velocities of the bodies using Jacobi Coordinates, which makes them all relative to each other in a recursive way. We note that the rigid-body components are not affected by this transformation, which only operates on the centers of mass of each body. Thus, we will generally only discuss the rigid-body terms when necessary and focus on the locations and velocities of the centers of mass of the bodies.

The Jacobi transformation
The Jacobi formulation defines a sequence of transformations where the position and velocity of each body P i is measured from the collective center of mass of the bodies P j , j = 0, 1, . . . , i − 1. The numbering of the bodies is arbitrary, a fact we use later, so without loss of generality we keep the current numbering of the system. To start the sequence, we define: Then the remaining Jacobi coordinates are defined as all for i = 1, 2, . . . , N . We can also find the convenient relationships The center of mass position and velocity vectors R C i and V C i are computed accounting for all bodies with that index and lower. The relative position and velocity vectors R i and V i are of body i relative to the center of mass of all bodies with index i − 1 and lower.
A special point to make is that the vector R C N is the total center of mass of the entire system. Keeping with usual practice, we can arbitrarily set this position to zero and have it remain fixed by the choice R C N = V C N = 0. We note that this assumption need not be made, and then the system center of mass would instead translate uniformly in space, with all other positions and velocities only providing relative information between each other (Fig. 1).

Angular momentum, energy and moments of inertia
The main advantage of the Jacobi coordinates are that they explicitly decouple the key quantities of angular momentum, moment of inertia and kinetic energy at a given order from Then the system total angular momentum, kinetic energy and inertia dyadic (ignoring the rigid-body contributions for the moment) are: Then the total system energy can be defined as where U i = i−1 j=0 U ji . All of these results can also be generalized to the case where each body is a rigid body. To do this at each order i, we just add to the definition of the quantities Fig. 2 Splitting a six-body system into two-to three-body systems and their interaction whereT i is the attitude dyadic that takes the body frame into an inertial frame, i is the angular velocity vector of the ith rigid body,Ī R i is the inertia dyadic of the rigid body in a body-fixed frame and U ii is the self-potential of the rigid body, which is in general nonzero for finite density mass distributions.

Jacobi transformation for multiple clusters
With the above formulation, we can see that a full set of N + 1 bodies can be represented by its total angular momentum and energy broken into two main terms, an internal term N i=1 H i and N i=1 T i and an external term for the angular momentum and kinetic energy of the center of mass The internal terms at a given index i capture all interactions within this group of masses and has no involvement with any bodies outside of it.
Given this observation, it is simple to generalize the Jacobi transformation to a set of M collections of bodies, P(I α ) for α = a 1 , a 2 , . . . , a M with each collection P α having n α bodies. Then the total angular momentum, energy and inertia dyadic for each cluster P α can be defined as the sum of the individual bodies in each component H i α , E i α andĪ i α from i α = 0, 1, 2, . . . , N α and yielding a total H α , E α andĪ α for each cluster.
Having defined the clusters, we describe how they are combined. We will just consider the case when we can separate our problem into two clusters. The mass and center of mass of each cluster are M α and R C α , respectively, where α = a, b. These two clusters can then be combined as described in the following and illustrated in Fig. 2.
The total mass of the system is M a + M b = M. The total center of mass of the system is and similar for the velocity. The final relative state between the collections is then and the associated relative angular momentum, energy and inertia between these collections are where we note that U ab = n a i a =1 This means that the gravitational potential contains all of the mutual potentials between the bodies in P a with the bodies in P b .
Then the total angular momentum, energy and inertia of the entire system can be expressed as This can be easily generalized to the combination of multiple clusters.

Disaggregated amended potential
Now recall the total angular momentum of the system H = HĤ, where the magnitude and the unit vector direction of this constant are important to define. Define the projection of an angular momentum H i and inertia dyadicĪ i onto the total direction of the angular momentum, H.
We note the useful identities which are easily established.
We first prove Lemma The following inequality can be established: Proof We refer to the proof of Theorem 1 in Scheeres (2002) for the approach. We do not explicitly add in the rigid-body rotation terms, however they can be accommodated as shown in Scheeres (2002). First note that It is important to note that the energy and angular momentum defined here, E i and H H i , are not constants of motion and that the bound involves all of the relative positions with respect to all of the bodies with an index j ≤ i, although they are independent of all bodies j > i. With this we can prove our main result Theorem The amended potential can be bounded by a disaggregated summation of amended potential-like terms, while still retaining its major inequalities.
Further, this inequality is sharp, meaning that both sides of the inequality can be equal at an appropriate state.
Proof First note that N i=1 E i = E is established by the property of the Jacobi Transformation, under the assumption that the total center of mass is stationary at the origin. Also from the Lemma we see that the last inequality on the right holds. In addition, we note that U = N i=1 U i ; thus, the gravitational potential can be removed from the remaining inequality. Thus, it is only necessary to prove the following inequality: On the left-hand side, we have the identity: On the right-hand side, we pull out the j = i term of the I H j sum to get ⎛ The sums of the H 2 H i cancel on each side of the inequality, and reordering the inequality, we find Finally, this is easily factored, leading to which is trivially true, establishing the result.
To establish that these inequalities are sharp, we need to evaluate the system at a condition that yields H H i for all i and j. We must also find a condition where Such conditions can be conveniently found when the system is in a relative equilibrium, which as established in Scheeres (2002) can occur for a finite density system at any given value of angular momentum.
When in a relative equilibrium, the system will act as a single body with all components rotating at a uniform rate about a single axisĤ. Thus, the individual terms H H i =Ĥ · H i = H ·Ī i ·Ĥ = I H i . Substituting this into the sharpness condition then yields I H i I H j = I H j I H i , trivially true and establishing the sharpness of the lower bound, or E = N i=1 E i when in a relative equilibrium. For the energy terms, the gravitational potentials are trivially equal and we just need to note the following equalities: ·Ī i · , as the shared angular velocity rotates about the angular momentum direction. Thus, the equality E i = E i is established, and the upper sharpness condition.
For an even shorter version of the proof, we can use the fact that when a system is in a relative equilibrium the equality E = E must hold.

Implications
For convenience we define the individual amended potential E i as It is tempting to give this function some dynamical meaning; however, while this expression is bounded by E i , this individual "amended potential" does not have all of the same properties as the full system amended potential E. First, the quantities H i and E i are not constant and vary with the entire system, with only the total H and E remaining constant. Further, the function E i does not only depend on body i alone, but also has terms that involve the positions and attitudes for j ≤ i. What can be asserted, to benefit, is that E i is independent of all terms j > i, and thus these higher-order terms do not directly affect the current E i term.
Where this result becomes more powerful is if some collection of bodies escape from the system. To start with we will explore what constraints we retain when a single body escapes. Then we will give the generalizations for when rigid bodies are involved. Finally, we will consider what happens when clusters mutually escape.

Escape of a single body
As explained in Scheeres (2020), if the system energy is sufficiently high it is possible for a single body to escape, separating the N + 1-body system into two collections escaping from each other, a system with N bodies and a single body, with their relative distance becoming arbitrarily large and asymptotically approaching ∞ in the limit. Conditions for this to occur have been worked out in detail for the point-mass n-body problem (Patnaik 1975). It is important to note that such escapes can occur even if the total system energy is negative.
As our ordering of the Jacobi coordinates was arbitrary, and the coordinates can be reordered at will, we will do so making the body i = N mutually separating from the other N bodies, which for the moment we assume remain at finite distance from each other. We can then split the system into two components (Marchal and Saari 1976) The position and velocity R N and V N are, by definition, the mutual position and velocity of body N with respect to the center of mass of the other N bodies. Thus by definition we see R C Then the position of the remaining N -body center of mass and the relative position of the single body are It is significant to note that neither R C N −1 or V C N −1 are included in the inequality, and thus they don't impact it at all.
When mutual escape occurs between body N and the remainder N bodies we have |R N | → ∞ and |V N | → V ∞ , but that in general all other |R i<N | < ∞. Thus, the mutual potentials between body N and the other collection all collapse to zero in the limit. Thus, U j N → 0 for all j < N . As I H N → ∞ in addition, we see that E N → 0.
We note that the total potential energy of the system then approaches U → 0≤ j<i≤N −1 U ji = U 0≤ j<i≤N −1 , with the potential energy relative to body N disappearing. Also, since I H N → ∞, we also have I H → ∞ and the amended potential function reaches the limiting form E → U 0≤ j<i≤N −1 .
The energy and angular momentum of the mutual system have a well defined limit, however, which we find from the standard two-body problem: where V ∞ is the hyperbolic excess speed and b ∞ is the hyperbolic semi-latus rectum, also called the "impact parameter." Recall that M N −1 = N −1 i=0 m i and M N = M N −1 +m N ; thus, the mass term is the reduced mass of the system computed between the two collections of bodies. Here we are also using an "osculating" two-body hyperbola to describe the magnitude of the angular momentum, using the parameters V ∞ and b ∞ which are only defined when the distance between the components is in the limit. We note that the mutually escaping bodies will also have a direction to the angular momentum,Ĥ N . All of the values V ∞ , b ∞ andĤ N cannot be predicted so long as escape was possible; however, it is possible to place constraints on the energy and angular momentum that can be lost using the limits defined by the amended potential. Given the chaotic nature of the n-body problem, the values within these bounds can at best be described in a probabilistic sense, or found from specific numerical simulations as shown later in the paper.
It is instructive to rewrite the theorem in the asymptotic escape limit. This gives us We explore each inequality individually. First is , meaning that the escape energy term is extraneous. Thus, we can rewrite this inequality into a more useful form: This explicitly shows that the total energy of the remaining N -body cluster is reduced by the mutual escape energy. If this process repeats itself, one ejected mass at a time, it is possible for the remaining bodies to eventually have a low enough energy so that they cannot escape (if they are of finite density), see Scheeres (2020). The argument is not as simple for the angular momentum. The mutual angular momentum of the escaping system approaches a limiting value which represents a "loss" of angular momentum from the total H. If the single body carries no angular momentum with itself (which will not be the case if it is a rigid body or if it is a collection of escaping bodies), then the remaining bodies will carry the remaining angular momentum H = H−H N . As the angular momentum is a vector, this new angular momentum can actually increase, decrease or stay the same in magnitude. Also, its direction in inertial space can also shift, depending on what the angular momentum of the escaping hyperbola is. There will be limits on what the final angular momentum and energy can be, shown with examples later. One issue is that the new "remaining" amended potential E 1→(N −1) is computed with the original angular momentum magnitude and direction, and thus is not the proper value to be used for the amended potential to find the sharpest limits. Nonetheless, from the Lemma the inequality E i ≤ E i still holds.
Still, to restart the analysis we can restate the problem with the new total angular momentum. Given the new n-body system (for this example with N bodies) with a total energy and angular momentum E and H , the new amended potential can be stated as: Starting from this, the system can be analyzed as before, the amended potential can be split again, following the Theorem, and the analysis can continue until the next escape occurs. This whole argument is justified by the results shown in Marchal and Saari (1976). The point is clear, however, in that the energy and angular momentum between the escaping components is lost to the system, especially to the remaining collection of N bodies. And that the overall energy for the remaining system to evolve has been reduced, independent of what its change in angular momentum has been.

Consideration of rigid bodies
If we extend this scenario to involve finite densities and hence rigid-body rotations, then nothing is essentially changed except that the remaining collection of N bodies will also have internal dynamics involving rotational motion, and that the escaping body can now retain some energy and angular moment in and of itself. In the limit when the bodies undergo their mutual escape, body N will also carry with it an energy E N = 1 2 N ·Ī N · N + U N N and will also carry with it some angular momentum H N =T N ·Ī N · N . In the limit the single rigid body will conserve these quantities as it will be isolated, and thus the energy retained by the remaining collection will be Each of the terms in these equations will be constant in the limit, discounting energy dissipation effects in the rotating rigid body and the remaining collection of masses. We note that for asteroid pairs, this sort of constraint can lead to a body rotating slowly after being ejected from an asteroid system (Pravec et al. 2010(Pravec et al. , 2018.

Escape of a collection of bodies
The above analysis holds for a single particle ejection, or for multiple particle ejections, so long as the bodies all escape from each other and from the remaining mass collection. However, it is also possible for collections of bodies P α to escape from each other yet the bodies remaining in each collection remain bound to each other as a group, as discussed in Marchal and Saari (1976).
We analyze the case if two collections of bodies escape from each other, but themselves stay relatively bound. First, both collections are aggregated into two disjoint groups we call P a and P b , and each of these collections have their total masses, M a and M b , and their centers of mass and velocity, R C a , R C b , V C a , and V C b . Forming the system center of mass, which is chosen to lie at zero, gives us the constraint M a R C a + M b R C b = 0, and similarly for the velocities. The relative position of the two collections is then called R ab = R C b − R C a and similarly for the velocities. The final angular momentum term becomes The gravitational potential U ab are the mutual potentials between all of the P a bodies with the P b bodies. Finally, the additional amended When escape occurs we have |R ab | → ∞ leading to U ab → 0 and I H ,ab → ∞. Then the inequality will become: In this situation the three energy terms on the right are all decoupled from each other and independently constant. Even though we cannot predict what they will in general be, they must sum to the total energy. Thus, the system can be decoupled into two separate terms, leading to Each of the collections P a and P b will have complex internal dynamics that can consist of coupled translational and rotational motions. In this initial form, the new inequalities will not be sharp. They can be modified after defining the new decoupled angular momenta H a and H b , defined such that

Evaluations of the energy and angular momentum loss
While the inequalities defined and applied above provide real constraints on the systems, there are elements of the final terms that cannot be predicted analytically. This because the escape process is inherently chaotic and in general cannot be predicted. This is an extremely important point and raises questions about what the likely statistical properties of the energy and angular momentum that are "lost" when escape occurs.
In the following we present some definite scenarios that allow us to accumulate some statistics on this question. To simplify the process, we only look at the interactions of a traditional point-mass n-body problem with equal masses. However, we apply many of the observations defined above in our simulations, including the use of Jacobi coordinates that are reordered once escape occurs, and explicitly tracing the energy and angular momentum that is removed from the remaining bodies. The simulations provide useful insight into the likely amount of energy and angular momentum that may be shed from a system when components are ejected.

Initial system configuration
To simulate an n-body system, an initial condition for the system is needed. In an attempt to obtain consistent results for analysis between n-body systems with different values of n, the same approach was used to define the initial system configuration used in the simulations. To avoid confusion where systems of N + 1 masses were considered earlier, a n-body system here will have N + 1 bodies.

Symmetric central configuration
For each n-body system, the positions of each body were generated by placing the bodies at the vertices of a symmetric n-polygon where each body was separated from the two bodies next to it by a length of l cc = L. This choice was made as this is an inherently unstable central configuration, implying that random deviations from it will also be unstable, allowing us to explore different possible evolutions for the same system. The value of L was selected to be the same value for all systems. The polygon was generated so that all the vertices were located in the x y-plane of an arbitrary coordinate system with unit vectorsî in the x-direction,ĵ in the y-direction, andk in the z-direction. Based on the specified value of L and n, the positions of each body at some initial time t 0 were calculated using Eq. 68. A diagram showing the locations of the vertices for the symmetric n-polygons generated using Eq. 68 is provided in Fig. 3.
The system was given a central configuration spin rate k which can be calculated using Eq. 70 (Scheeres and Vinh 1991). The initial velocities for each body were then calculated

Nondimensionalization
The dimensional quantities of length, time, and mass were nondimensionalized by some constant, i.e., q = Cq andq = C −1 q where q is the dimensional quantity,q is the nondimensional quantity, and C is the constant. Length was nondimensionalized by the initial distance between the masses, i.e., d = C dd where C d = L is the initial distance between the masses before the perturbation. Time was nondimensionalized by the initial central configuration spin rate, i.e., t = C tt where C t = −1 is the initial central configuration spin rate. Mass was nondimensionalized by the total system mass, i.e., m = C mm where C m = M is the total system mass. Based on these three conversions, some additional conversions are: C R = L for position, C V = L for velocity, C A = L 2 for acceleration, C E = M L 2 2 for energy, C SE = L 2 2 for specific energy, C H = M L 2 for angular momentum, and C S H = L 2 for specific angular momentum. Up to this point, all quantities have been presented in their dimensional form. For the remainder of this paper the majority of quantities will be presented in their nondimensional form. The relevant dimensional form can be obtained by using the nondimensionalization constants presented in the previous paragraph.

Initial condition variations
Each system was then given a unique initial condition by shifting the position of one of the bodies. No change was made to the initial velocities. As a result, for all simulations in this analysis V C N = 0. The perturbation was implemented so that the position of P N was shifted to a new position within a uniformly distributed sphere with a radius of L/2 centered at the original position of P N . The positions of the bodies were then translated so that the COM of the perturbed initial system was at the origin, i.e., R C N = 0. The resulting initial condition was then used to generate a simulation for the motion of the perturbed system. Aggregate statistics for each n-body system were generated by repeating this process multiple times for each n-body system.
While a random perturbation was applied in this analysis, there is not expected to be any significant changes between these results and the results that would be obtained from simulations where a systematic perturbation is applied. This expectation is based on the randomness of n-body motion and was confirmed based on preliminary results obtained from analysis on three-body systems where a systematic perturbation was applied to one of the masses first in the y-direction in one set of simulations and then in the z-direction in another set of simulations. The results for these two sets of simulations were similar to the results for the random sets for 3 ≤ n ≤ 6.
In the "Appendix" we describe the integration routine used for the simulations. We also detail the conditions used to determine when a single or multiple set of bodies have been mutually ejected.

Final system configuration
As mentioned previously, a n-body system with n > 2 can have bodies that are ejected from the system. As the total system energy was negative for all simulations based on the initial system configuration used, bodies will be ejected, but at least one pair of bodies must remain together as complete disaggregation requires a positive or zero total energy (Marchal and Saari 1976). For this reason, a simulation was considered to have reached a "complete" end condition when n − 2 bodies had been ejected from the system. If less than n − 2 bodies were ejected up to the point the maximum integration time was reached, the simulation was considered to have reached an "incomplete" end state. It is possible to have cases where n ejc < n − 2 but the "remaining" system is a long-period system with n rem > 2 bodies. In some cases, these remaining systems with n rem > 2 appear somewhat "stable" over the time scale considered in this analysis. For that reason, any simulation where at least one body was ejected from the system were included in the results presented in this paper. However, often times the results from the "incomplete" simulations and the results from the "complete" simulations are distinguished from each other.
There were also a few cases where the integration tolerance could not be met over the time range considered, due to close approaches between the bodies. The simulations where this occurred were considered to have reached a "failed" end state. In these cases, if any bodies had been ejected before the integration tolerance could not be met, the ejection was included in the results as part of the "incomplete" simulations. However, the characteristics of the "remaining" system in these cases were not included.
A subscript of C will be given to terms related to the separation of the two clusters, and a subscript of S will be given to terms related to the internal interaction of a cluster. Please refer to the "Appendix" for more detail on this notation. If there were multiple system ejections detected at the same time, the bodies were ordered so that the bodies in the ejected system with the highest energy E S were given the highest indexes when reordering the bodies. While any "ejected system" could technically be considered the "remaining system" by reordering the bodies so that the indices of the bodies in the ejected system were 1 and 2, the order where the bodies in earlier detected system ejections are given the highest indexes was used. In terms of implementation, if n rem = 3 only single ejection events were allowed to occur so that there was always a remaining binary system. As a result, system ejections were not considered for the n = 3 systems because if a system ejection event was detected, the bodies in the ejected system were reordered as P 1 and P 2 and therefore were viewed as the "remaining" system.

Computational performance and end conditions identified
In this analysis n-body systems for n = 3 to 6 were considered, with 1000 simulations for n = 3 to 5 and 400 simulations for n = 6. The aggregate statistics resulting from these simulations were analyzed to determine if there were any similarities in the results between n-body systems with different values of n.
Several MATLAB scripts were written to generate and analyze the simulation results. These scripts can be classified as those related to integrating the system and those related to processing the results from the integration.

Initial system states
The procedure discussed in Sect. 6.1 was used to generate the initial conditions used in the simulations. The position space representation of the initial conditions used in the simulations is provided in Fig. 4 . Plots showing the total energy and total angular momentum for the systems used in each simulation are provided in Fig. 5.

Program execution time
Data related to the execution time of both the integration and processing routines used to generate and analyze the results, as well as the number of simulations that achieved each end condition, are provided in Table 1. As expected, the execution time for the integration routine was much longer than the execution time for the processing routine. It should be noted that even though parallel loops were used in both the integration and processing routines, a small number of the simulations accounted for a considerable portion of the total integration time. In the n = 5 simulations, six of them had an integration execution time that were at least 5% of the total execution time for all 1000 of the n = 5 simulation. In the n = 6 simulations, this was true for 16 of the 400 simulations.
Regarding the number of simulations that achieved each end condition, as the number of bodies in the system increased, the percentage of simulations where a complete end condition was reached decreased, as did the percentage of simulations where no bodies were ejected. Also, no simulations reached a failed end state for n = 3. For the other values of n, around 1% of the simulations reached a failed end state. One possible explanation for these simulations reaching a failed end state is a close approach between two of the bodies that is too small to meet the predefined integration tolerance. While this could be remedied by using a regularization approach, the numbers involved were small enough to not do this.

Timing of ejection events
As anywhere from 13% to 52% of the simulations reached an incomplete end condition, it may be natural to wonder if this number could be decreased by simply increasing the maximum allowed propagation time so more ejections can be detected. To evaluate whether extending the integration time would yield a useful number of additional results, an exponential curve was fit to the number of cumulative ejections detected at each time across all simulations for each value of n with a model of the form f (τ ) = a 1 − e −bτ . The data used in the  The fit function is of the form f (τ ) = a 1 − e −bτ and represents the number of ejections at some time τ normalized by the total number of theoretically expected ejections. Only the data corresponding to the second half of expected ejections were used in the fit. The 95% confidence interval for the values of the coefficients are given in parentheses fit function was normalized by the total number of bodies that should be ejected across all simulations for each value of n, i.e., the number of cumulative ejections for each collection of n-body simulations was divided by s(n − 2) where s is the number of simulations for that particular value of n. Using this form, the coefficient a should be a number between 0 and 1 that represents the fraction of bodies that were ejected compared to the number of bodies that were expected to be ejected. As the number of ejections later in time were of more interest than when ejections occurred on a shorter time-scale, another fit of the same form was determined using only the data corresponding to the second half of all expected ejections. The results from the second fit were selected and will be referred to as the "final fit". The coefficient b acts as a time scaling parameter. While the fit analysis did yield slightly different coefficients for different values of n, the values of the coefficients were in similar ranges. Also, the R 2 value for the final fits were at least 0.99 for all values of n considered in this analysis. The fits are provided in Fig. 6, and the coefficients for the final fit are provided in Table 2. The results indicate that increasing the simulation time may not result in significantly more ejections, as the number of ejections is predicted to slow exponentially. Thus, we consider the stopping conditions we've chosen to be a reasonable compromise. By manipulating the fit models for an n-body system with fit coefficients a and b, the additional time required for n ejc more bodies to be ejected in s simulations can be determined based on time using Eq. 71a or Eq. 71b if the number of bodies that have been ejected n ejc across the s simulations is already known. Based on the form of these equations, as the number of ejected bodies increases, the expected time for the next body to be ejected increases as well. For example, using the coefficients for the n = 6 systems, once one body has been ejected, the expected additional time until the next body is ejected is τ = 326. However, for the same system, once three bodies have been ejected, the expected additional time until the fourth body is ejected is τ = 3624, more than ten times the original expected time and on a similar scale as the original total propagation time used in this analysis. It should be noted that these calculations are based on the aggregate statistics used for the fit data and as such there is no guarantee that these results directly apply to any single simulation or accurately represent the behavior of ejections that occur relatively early on in the simulations. Also note that these equations do not extend to situations where n ejc > as(n − 2) or where n ejc + n ejc > n − 2.
Several simulations were run for systems with n = 8. The computation time required for these simulations was significantly longer than the simulations for n ≤ 6. Not enough of these simulations were completed to present aggregate statistics, but of the five simulations for n = 8, only one reached a complete end state and that simulation had only single-body ejection events. We note that these exponential fits are only to guide our choice of total integration time, and do not have any implications about the long-term behavior of such chaotic dynamical systems.

Characteristics of ejected bodies and systems
The characteristics of the orbits related to the ejected bodies and the remaining systems will now be analyzed. Note that for all histograms in this subsection where the y-axis is the percent of simulations, that percentage is with respect to the number of simulations that met the conditions to be plotted, not the total number of simulations.

Hyperbolic orbital parameters of ejected bodies and systems
Whether a single body or a system are ejected, there will be a part of the total system energy that can be used to describe a two-body orbit between the ejected body(ies) and the remaining bodies. If the ejected body(ies) are not to return to the remaining system, then this orbit must be parabolic or hyperbolic. For this reason, the hyperbolic excess velocity and impact parameter of these orbits were calculated for every detected ejection and the aggregate results are provided in Figs. 7 and 8 . Note that these histograms are normalized based on the total number of ejections detected for each value of n. As can be seen in Fig. 7, the distributions of V ∞,C and b ∞,C are similar for simulations with different values of n. The peaks in the V ∞,C distributions were around approximately 0.2 to 0.5 for n = 3 and 0.4 to 0.6 for 4 ≤ n ≤ 6 and the peaks in the b ∞,C distributions were around approximately 2 to 6 for 3 ≤ n ≤ 6. Also, the b ∞,C distributions tended to have longer tails than the V ∞,C distributions. When considering the ejections across all simulations for a single value of n, the distributions of the parameters for single-body ejections and system ejections were similar to the overall distribution depicted in Fig. 7. This similarity was also identified when analyzing the parameters based on the order of the ejection event. 2 V 2 ∞,C (bottom) for 3 ≤ n ≤ 6 Simulations. Note that the x-axis ranges for the n = 3 and n = 6 plots are different than for the n = 4 and n = 5 plots In Fig. 8 the b ∞,C and V ∞,C terms are combined into the hyperbolic angular momentum b ∞,C V ∞,C and escape energy 1 2 V 2 ∞,C to track these quantities in the ejected bodies. We see that the n = 3 distributions are unique relative to n > 3, and that the distributions look similar for the larger numbers of bodies. The prevalence of low angular momentum and low energy in these ejections is clear and will show up later when the total system angular momentum and energy is considered.

Characteristics of ejected and remaining systems
When a system is ejected, there will be another part of the total system energy that can be used to describe a two-body orbit between the ejected body(ies). If the ejected body(ies) Fig. 9 Semi-major axis a S (top) and eccentricity e S (bottom) for 3 ≤ n ≤ 6 simulations. The distributions for the ejected systems are in green. The distributions for the remaining systems are in blue for the complete simulations and light blue for the incomplete simulations. Note that system ejections were not possible for n = 3 simulations, and that the x-axis ranges may be different are not to be ejected from each other, then this orbit must be circular or elliptical. As the ordering of the masses is arbitrary, the orbits of a remaining system consisting of just P 1 and P 2 in the updated order of bodies were included in these results. These two-body remaining systems were included for all simulations except for those that had failed end conditions, or the orbit of those two bodies was parabolic or hyperbolic. The subscript of S in these results reference the collection of these two-body remaining systems and the ejected systems, not just the ejected systems. Note that these histograms are normalized based on the sum of the number of system ejections detected and the number of these two-body remaining systems considered for each value of n. The semi-major axis and eccentricity of these orbits were calculated for every detected system ejection and the aggregate results are provided in Fig. 9 .
The distributions of a S and e S were similar for simulations with different values of n. The distributions of a S had a longer tail on the right-hand side of the peaks than the left-hand side, most likely in part due to the minimum allowed value of zero. The e S distributions were in the range of 0 ≤ e S < 1 as expected, and as the value of e S increased the relative rate of occurrence increased. However, even though the distributions of the semi-major axis values and the distributions of the eccentricity values had similar shapes for different values of n, only the eccentricity distributions covered the same range of values for all three types of systems described in the figures. Looking at Fig. 9, it is apparent that the distributions of the remaining systems in both the complete and incomplete simulations span similar ranges while the distributions of the ejected systems spans a different range. For this reason, and because system ejections were not possible in n = 3 simulations, the semi-major axis values in Fig. 9 were separated, with the values for both the complete and incomplete simulations and the values for the ejected systems provided in Fig. 10.
As before, in Fig. 11we combine the semi-major axis and eccentricity of the binary systems to compute the energy and angular momentum of the binary systems. Again, the case of n = 3 and n > 3 shows distinct distributions, with the larger number of bodies showing similarities between their results. Fig. 10 Semi-major axis a S for the Remaining Systems (top) and for the Ejected Systems (bottom) in the 4 ≤ n ≤ 6 Simulations. The distributions for the remaining systems are in blue for the complete simulations and light blue for the incomplete simulations. Note that the x-axis range for the n = 4 plot is different than for the n = 5 and n = 6 plots (top) and angular momentum μa S 1 − e 2 S (bottom) for 3 ≤ n ≤ 6 Simulations. Note that the x-axis ranges for the n = 3 plot is different than for the n = 4, n = 5, and n = 6 plots Fig. 12 Graphic showing the vector geometry between the total, ejected and remaining angular momentum (left) and the sum of the total, ejected and remaining system energy (right) Fig. 13 Taken energy for 3 ≤ n ≤ 6 simulations. The distributions for the complete simulations are in red and the distributions for the incomplete and failed simulations with at least one ejection are in pink. Note that the x-axis range for the n = 3 plot is different than for the n = 4, n = 5, and n = 6 plots

Partitioning of conserved quantities
The partitioning of the total energy and total angular momentum will now be discussed. Figure 12 shows the basic geometry and constraints between the total energy and angular momentum, the energy and angular momentum lost to the ejection processes, and the energy and angular momentum that remains in the final orbiting systems.

Partitioning of energy
With respect to energy, looking at the total amount of energy in the E C terms for each ejection corresponds to the energy lost due to ejections. The total of the E S terms related to the remaining system and the ejected systems (if applicable) is the amount of energy that is kept active in the system. Note that the sum of these two energies should equal to the total system energy. The percent of the total energy taken by the ejected bodies in each simulation is shown in Fig. 13 . This is computed by dividing the energy E C by the total initial energy E and multiplying by 100. Note that the values are negative as the total system energy was negative for all simulations and the taken energy was always positive. In this figure red represents the results for the completed simulations. For both the incomplete and failed simulations where at least one body was ejected, there is an amount of energy that Fig. 14 Taken angular momentum for 3 ≤ n ≤ 6 simulations. The distributions for the complete simulations are in red and the distributions for the incomplete and failed simulations where at least one body was ejected are in pink. Note that the x-axis range for the n = 3 plot is different than for the n = 4, n = 5, and n = 6 plots is taken, but it is not the final amount of energy that would be taken if the other n rem − 2 bodies are eventually ejected. For this reason the value of taken energy in these simulations are included in the results, and are represented by the pink bars in the figure to distinguish them from the values for the complete simulations.
The distribution of the percent of total energy taken by the ejected bodies is similar between systems with different values of n. However, as the number of bodies in the system increased, the percent of total energy taken by the ejected bodies appeared to increase slightly. For all values of n examined here, as the magnitude of the energy taken increased (i.e., the percent of energy taken got more negative), the number of simulations that produced that result decreased. As expected, the magnitude of energy taken by the ejected bodies was generally smaller for the incomplete simulations than for the complete simulations. Also, for the n > 3 systems, there were a few simulations where the taken energy was larger in magnitude than the magnitude of the total system energy.

Partitioning of angular momentum
For the angular momentum, the sum of the H C vectors for each ejection was computed, and the magnitude of that resulting vector was the amount of angular momentum taken by the ejected bodies. The vector sum of the total taken angular momentum, the angular momentum in the remaining system, and the angular momentum stored in ejected systems should equal the total system angular momentum vector. Note this does not mean the sum of the scalar magnitudes of those vectors will equal the magnitude of the total system angular momentum. With that being said, the magnitude of the taken angular momentum vector was compared to the magnitude of the total system angular momentum vector when computing the percent of the total angular momentum taken. The total angular momentum taken in each simulation is shown in Fig. 14 . In this figure red represents the results for the completed simulations and pink represented both the incomplete and failed simulations where at least one body was ejected.
As was the case for the taken energy, the distribution of the percent of total angular momentum taken by the ejected bodies is similar between systems with different values of n. The distributions of the taken angular momentum appear to be more Gaussian in nature with peaks between 80 and 105% of the total angular momentum. It is interesting to note that for a considerable number of the simulations, the magnitude of the taken angular momentum was larger than the magnitude of the total angular momentum. These distributions had a longer tail on the left-hand side of the peaks due to the incomplete and failed simulations. This behavior was expected for these simulations as the taken angular momentum could still increase as more bodies are ejected at later times. We also note that for n > 3, as the systems Fig. 15 Angle between the remaining h and total H for 3 ≤ n ≤ 6 Simulations. The distributions for the complete simulations are in blue, and the distributions for the incomplete simulations are in light blue eject more bodies their distribution becomes more focused about total (100%) loss of angular momentum. We note that this can occur while the ejected orbital systems still have nonzero angular momentum due to the vectorial nature of angular momentum.
The direction of the angular momentum is also of interest. H 1 is the angular momentum corresponding to the interaction between just P 1 and P 2 in the updated order of the bodies if all other bodies are ejected. For that reason, the angle between H 1 and H (represented by θ ) was calculated for each completed simulation in order to analyze how the ejections affect the direction of the angular momentum in the remaining system. If the system ejected some bodies but did not reach completion and had m bodies left, the angle θ is computed as the angle between H m−1 and H. The results for both sets of angles are provided in Fig. 15. The angles resulting from complete simulations are in blue. The angles resulting from incomplete simulations where at least one body was ejected are shown in light blue. It's important to note that the manner in which this quantity is tracked means that the orientation of any ejected binary system is not represented.
Unlike for the taken energy and angular momentum, the distribution of θ was not the same between systems with a different value of n. For the n = 3 systems, almost none of the simulations produced a result where θ > 90 • , and the majority of the simulations produced a result where θ ≤ 30 • which was the most common result. This is likely related to the well-known restrictions on inclination for the three-body problem (Marchal and Saari 1975). As n increased, the distribution of θ shifted rightward and flattened out. For n = 4, the most common result was θ ≤ 30 • , but this only occurred in approximately a quarter of the simulations considered where for n = 3, θ ≤ 30 • in approximately two-thirds of the simulations considered. The most common result for 40 • ≤ θ ≤ 70 • for n = 5 and 70 • ≤ θ ≤ 100 • for n = 6. Also, there were a considerable number of simulations where θ > 90 • for systems where n > 3. As a result of the initial conditions used in the simulations, the total angular momentum vector had a large +z-component and smaller x-and y-components. While the z-component was not the only non-zero component, as the total angular momentum was close to being perpendicular to the plane where the bodies were initially located, if θ was greater than approximately 90 • , then H 1 had a negative z-component, implying that it was rotating in the opposite sense of the initial configuration.

Summary of results for systems with different n
For the majority of the quantities examined in this analysis, there was a consistent pattern in the distributions between systems with different values of n for 3 ≤ n ≤ 6. First, the energy always decreased when bodies or systems were ejected. This is explicitly predicted by the bounds derived earlier and indicates that ejections will systematically decrease the level of energy in the remaining systems, thus potentially allowing them to stabilize by shedding particles. While the total percentage of initial energy lost is not a meaningful scientific term, it does track the magnitude of energy that can be lost in this way. While there were many systems that reduced their energy by a significant amount, the most common outcome is clearly for the systems to only decrease their energy by a relatively minor amount, with the taken energy near zero but negative. This can provide a potentially useful analytical approximation when considering ideal cases of ejection. When rigid bodies effects are included in these systems it is significant to note that many ejected systems will take rotational energy from the rotating bodies to enable escape. This has been predicted analytically and has been seen in clusters of asteroids that have undergone mutual escape, and which have the signature of a slower rotation rate (Pravec et al. 2018).
While the energy had a systematic change due to body ejections, it is notable that the angular momentum has no such constraints. The lack of a constraint on angular momentum "loss" was noted analytically, however the numerical results clearly demonstrate this fact. The range of angular momentum that was taken by the system spanned the entire theoretical possibility for n > 3, with the extremes being no angular momentum lost up to 200% being lost (meaning that the lost angular momentum is completely reversed in its direction and magnitude from the initial system). For the case of n = 3 there seem to be stricter constraints on the angular momentum loss, with the most likely loss value being a bit less than 90% (see Marchal and Saari 1975). For the n > 3 case the most likely case was to lose 100% of the angular momentum. For both cases, it is important to note that the direction of the angular momentum loss was also quite variable, meaning that the individual systems could have a final angular momentum with a direction randomly distributed relative to the initial, total angular momentum.
Here it is relevant to talk through the angular momentum results. First consider a n = 3 case, where one body is ejected. The ejected body in our computation carries no angular momentum with it, so the total angular momentum equals the sum of the final angular momentum in the remaining binary system and the "lost" angular momentum of the mutual orbit of the escaping components. For the most typical n = 3 body case, the remaining binary contains around 10% of the initial angular momentum magnitude, and 90% of the angular momentum is taken by the mutual hyperbolic orbit, with these two being aligned (the case when θ = 0). If instead, we have a case where θ = 60 • , with the taken magnitude being 90% again, then the angular momentum remaining in the binary would have a magnitude of 74.5% of the initial value. As we go to n > 3, we see the most likely value approach 100% and the final direction becoming more uniformly distributed, with a notable excess emerging at 90 • as n increases to 6. Here, for an n = 6 case, consider that 4 single ejection events occur, leaving a final binary asteroid and four separate loss angular momentum vectors that must be summed to equal a net lost angular momentum of 110% and with the final binary system's angular momentum pointing perpendicular to the original angular momentum (θ = 90 • ). Then the angular momentum in the final system must have a magnitude 45.8% of the initial magnitude. For the same system, if the lost angular momentum magnitude equals the total angular momentum magnitude, and the final angle of the bound system is 60 • , then the diagram forms an equilateral triangle and the angular momentum H S will be equal to the original angular momentum magnitude. What both of these examples show is that a large loss of angular momentum does not mean that the final angular momentum in the ejected or binary systems need be small. Thus ejection does not have a systematic effect on the final angular momenta of the bodies, even though energy does.
As a final comment, we recall that in the two-body problem there is a simple constraint between energy and angular momentum. For a given binary system defined by an energy E S , an angular momentum magnitude H S and its total mass parameter μ we have From our above examples and discussion we note that the energy E S is consistently decreased, potentially by a large amount. On the other hand, we note that the angular momentum H S is not necessarily systematically driven to a small value, and can retain an appreciable fraction of the initial system angular momentum. The combination of these effects should push these systems closer to the circular orbit condition (consistent with the equality), or for finite density bodies potentially to a resting configuration. The manifestation of these constraints for finite density systems is an interesting question and should be investigated in future research. A question of specific interest is how these mechanics influence the creation of rubble pile bodies in the aftermath of cataclysmic collisions, which is the currently accepted model for how these bodies are created (Fujiwara et al. 2006).

Conclusions
In this paper constraints on the energy and angular momentum of an n-body system undergoing mutual escape are derived and studied. The paper presents and proves rigorous analytical constraints between the energy and angular momentum of the total system and sub-sets of the system. Then it provides some specific examples related to these results, filling in details where analytical constraints are not possible due to the chaotic nature of n-body dynamics.
The results show that energy is systematically lost during the ejection process, thus decreasing the overall energy available in any remaining n-body clusters that are not disrupted. In contrast, the examples and theory indicate that there are not such strong constraints on system angular momentum, due in part to its vectorial nature. Arising from this balance we show that ejections should in general stabilize remaining orbital systems.

Mass ejections
If a system consists of n > 2 bodies, then the system is in general unstable and body(ies) can be ejected until only single or pairs of bodies remain. With that in mind, a method to detect when a single body was ejected and when pairs of bodies were ejected was developed. A single body being ejected will be referred to as a single-body ejection event and a pair of bodies being ejected will be referred to as a system ejection event. While a much simpler evaluation is conducted in this analysis, conditions for detecting ejections of a body in the three-body problem and different system end states are discussed in Marchal (1974).

Single-body ejection events
For a n-body system, if one of the bodies is sufficiently far away from the other bodies, then the remaining n − 1 bodies can be approximated as a single body with the combined mass and averaged position and velocity of those n − 1 bodies. In this discussion of single-body ejections, "System A" and "System B" will refer to B A comprising of P 0 , . . . , P N −1 and B B comprising of P N , respectively. If the two-body system of B A and B B has a specific energy E ab(2) ≥ 0, then it can be concluded that P N has been ejected from the system and will not return to influence the motion of the other bodies. Granted, this is assuming that the ejected body P N is sufficiently far away from the remaining masses so that the two-body approximation is valid. The Jacobi formulation vectors are particularly useful in this ejection analysis as R N and V N are the position and velocity of P N relative to the COM of the other bodies. Therefore, R N and V N can be used as the position and velocity vectors in the two-body specific energy equation in Eq. 73.
To determine if P N had been ejected, three conditions were evaluated: Condition 1 is the most important criteria to determine if a single body has been ejected because when this condition is true, P N is on a parabolic or hyperbolic orbit with respect to the COM of the remaining n − 1 bodies. Condition 2 is used to ensure that the two-body assumption used in Condition 1 is valid. U i is a function of only the Jacobi Formulation vectors R k for k = 1 to i as can be seen in Eq. 74. Please note that the form of Eq. 74 assumes i > j.
If none of the masses P 1 through P i−1 have been ejected, then as P i is ejected, R i will be much larger than R k for k = 1 to i − 1. Considering Eq. 74, in this case R i − R j can be approximated as R i , and U i → − G M i−1 m i |R i | which is similar to the form of the potential energy for a two-body system. Condition 3 provided an additional check that was similar to Condition 1. Note that under this ejection assumption with i = N , M N M N −1 m N E N → E ab (2) . If all three of these conditions were met, then P N was considered to have been ejected from the system. If P N is ejected, E N and H N at that time step were treated as the energy and angular momentum taken by the ejected body. The parameters V ∞,C and b ∞,C were calculated for the orbit of the two-body system of B A and B B . These parameters are given a subscript C in the results to indicate they are associated with ejection orbits.
Once P N is ejected the "remaining" system consists of the remaining n rem = n −1 bodies. To determine if any of the remaining bodies are ejected at a later time, the same process can be repeated. However, instead of considering bodies P 0 through P N , the remaining system of P 0 through P N −1 can be analyzed as once P N has been ejected, the gravitational potential terms associated with the interaction between that ejected mass and the remaining n rem masses in the system go to zero. Once P N −1 is ejected, n rem = n rem − 1 and the "remaining" system consists of P 0 through P N −2 . This process can be repeated until n rem ≤ 2, i.e., the number of bodies ejected is n ejc = n − 2. A figure depicting this process is shown in Fig. 16.
As the initial ordering of the bodies is arbitrary, there is no guarantee that the body initially ordered as P N will be the first body ejected which adds an additional level of complexity to the problem. However, because the order of the bodies is arbitrary, at a particular time the bodies can be reordered in any way without loss of generality or accuracy. The ejection analysis can then be conducted on this reordered system and then parts of the new order can be kept, or the old order can be reinstated. As all integration was performed before the ejection analysis was conducted, this reordering process does not affect the results obtained from the integration of the system. As the simple relationship between the Jacobi formulation Fig. 16 Single-body ejection event for a n-body system with N = 6 Fig. 17 Reordering process in single ejection analysis for a n-body system with N = 6. The old order 012345 is rearranged as 401253 vectors and the two-body position and velocity in the specific energy equation are only valid for an n-body system if the body being considered is P N , the system can be reordered n times (including the original order) such that each of the n bodies is considered as the P N with respect to the other n − 1 bodies in the reordered system. In order to determine if a body has been ejected, the bodies are reordered at each time step so that each body P i is treated as P N , i.e., P 1 was treated as P N , then P 2 was treated as P N ,…then P N was treated as P N . A figure depicting this process for an initial system where n = 6 and i = 3 is shown in Fig. 17. Note that the order of the other bodies P j for 0 ≤ j < N in the reordered system do not matter.
If P i is ejected, the order used in future analysis is updated so that P i is reordered as P N . This process is repeated at future time steps and each time a mass is ejected the order is updated so that the newly ejected mass is ordered as P N rem −1 and the remaining masses are considered as a system with n rem = n rem − 1. This process is repeated until only two bodies remain or a maximum integration time limit is reached. So, in a system where only single-body ejection events occurred (i.e., no system ejection events occurred) the i th ejected body would be ordered as P N +1−i in the final order of the bodies.

Multiple-body ejection events
The case for a system ejection event is similar to the case for a single-body ejection event. For a n-body system, if two of the bodies are close to each other but are sufficiently far away from the other bodies, then the remaining n − 2 bodies can be approximated as a single mass with the combined mass and averaged position and velocity of those n − 2 bodies. In this discussion of single-body ejections, "System A" and "System B" will refer to B A comprising of P 0 , . . . , P N −2 and B B comprising of P N −1 , P N , respectively. If the two-body system of B A and B B has a specific energy E ab(2) ≥ 0, then it can be concluded that both P N and P N −1 have been ejected from the system and will not return to influence the remaining bodies. A process similar to the one used in the single-body ejection analysis was performed when analyzing system ejection events. However, there were several key differences. B B in this analysis contains two masses while earlier B B contained only one mass. As a result, there is a certain amount of energy and angular momentum stored in the orbit between P N −1 and P N . Also, the third evaluation criterion was changed to ensure that P N −1 and P N would not be ejected from each other.
Condition 1 is the main criteria used to determine if a system has been ejected, and is functionally the same as Condition 1 in the single-body ejection analysis. Condition 2 is used to ensure that the two-body assumption used in Condition 1 is valid. The first and second summation terms are to ensure that P N and P N −1 are both sufficiently far away from P 0 through P N −2 . Condition 3 serves a different purpose than the third condition used in the single-body ejection analysis. This condition is used to ensure that P N and P N −1 will remain in an elliptical or circular two-body orbit about each other, i.e., this condition ensures that P N will not be ejected from P N −1 . If all three conditions are satisfied, the resulting values for E ab and H ab were considered the energy and angular momentum taken by those two ejected bodies in System B. E S and H S (the quantities associated with the internal interactions in System B) were considered the energy and angular momentum stored by those two ejected bodies. The parameters V ∞,C and b ∞,C were calculated for the orbit of the two-body system of B A and B B . Regarding the orbit of the two-body system comprising of P N −1 and P N , the semi-major axis and eccentricity will be represented as a S and e S , respectively. Note that the S subscript is also used for quantities related to the "remaining" system, not just "ejected" systems. A figure depicting this process is shown in Fig. 18.
While there is no guarantee that the pair of bodies initially ordered as P N −1 and P N will be the first bodies ejected, the bodies can be reordered arbitrarily. The system can be reordered multiple times such that every pair P i and P j of the n bodies have the two highest indexes in the system (i.e., N − 1 and N ) in the reordered system. A figure depicting this process for an initial system where n = 6, and analyzing if P 5 and P 3 are ejected as a pair, is shown in Fig. 19. Note that the order of the other bodies P k for 0 ≤ k < N − 1 in the reordered system do not matter. Also, whether P i or P j is ordered as P N (and the other as P N −1 ) does not matter. As a result, for a n-body system there are 1 2 n(n − 1) orderings that need to be Fig. 18 System ejection event for a n-body system with n = 6

Fig. 19
Reordering process in system ejection analysis for a n-body system with n = 6. The old order 012345 is rearranged as 4012 (35) considered (including the original order). With that in mind, each pair is analyzed where P i is reordered as P N for i = 1 to N and the corresponding P j is reordered as P N −1 for j = 0 to i − 1.