{'fraction_non_alphanumeric': 0.041997396161437994, 'fraction_numerical': 0.020557725421023895, 'mean_word_length': 4.235598065083553, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability and chaos of hierarchical three black hole configurations

We study the stability and chaos of three compact objects using post-Newtonian (PN) equations of motion derived from the Arnowitt-Deser-Misner-Hamiltonian formulation. We include terms up to 2.5 PN order in the orbital part and the leading order in spin corrections. We performed numerical simulations of a hierarchical configuration of three compact bodies in which a binary system is perturbed by a third, lighter body initially positioned far away from the binary. The relative importance of the different PN orders is examined. The basin boundary method and the computation of Lyapunov exponent were employed to analyze the stability and chaotic properties of the system. The 1 PN terms produced a small but noticeable change in the stability regions of the parameters considered. The inclusion of spin or gravitational radiation does not produced a significant change with respect to the inclusion of the 1 PN terms.


I. INTRODUCTION
The dynamics of compact objects play an important role in the evolution of galaxies and other stellar systems. Post-Newtonian (PN) techniques are a useful tool for modeling the dynamics of multiple compact objects. In astrophysical models, the gravitational radiation is included via an effective force which considers 1 PN and 2.5 PN corrections, and in some cases stellar dynamical friction. Hierarchical three black hole configurations interacting in a galactic core have been the focus of recent studies by several authors. For example, intermediatemass black holes with different mass ratios are considered in [1][2][3]; simulations of dynamical evolution of triple equal-mass supermassive black holes in a galactic nuclei were performed in [4,5]. Other astrophysical applications of multiple black hole simulations include three-body kicks [6,7] and binary-binary encounters (see e.g. [8][9][10][11][12]). In the context of the final parsec problem [13,14], threebody interactions are considered as a mechanism that can drive a binary black holes system to a separation below one parsec.
Since the 1990s there has been an increasing effort to detect and study compact objects. The most likely source of gravitational waves are binary compact objects. Recently, it was shown that the probability of more than two black holes to interact in the strongly relativistic regime is, not surprisingly, very small [15]. For practical purposes, the creation of gravitational waveform templates for gravitational wave detectors is naturally focused on binary systems. Binary systems can produced complicate waveforms when taking into account spinning black holes and eccentric orbits, e.g. [16]. On the other hand, triple systems consisting of a binary black hole and a star (e.g., a white dwarf) are potential sources of electromagnetic counterparts associated with the mergers [17][18][19][20]. * Pablo.Galaviz@monash.edu The chaotic behavior of triple systems is well known in the Newtonian case (see [21] and references therein). For binaries, it is known that chaos appears when using certain post-Newtonian approximations for spinning binary systems (see [22][23][24][25][26][27][28]). In the present work, we studied three-body systems with PN methods, where the main technical novelty is the inclusion of the 2.5 PN terms in the orbital dynamics and leading order of spin-orbit and spin-spin terms [29][30][31].
Recently, similar PN techniques were applied to the general relativistic three-body problem. Periodic solutions were studied using the 1 PN and 2 PN approximations in [32][33][34]. Examples of three compact bodies in a collinear configuration were considered in [35,36], and Lagrange's equilateral triangular solution was studied including 1 PN effects in [37]. In [38], the stability of the Lagrangian points in a black hole binary system was studied in the test particle limit, where the gravitational radiation effects were modeled by a drag force. The waveform characterization of hierarchical non-spinning three-body configurations using up to 2.5 PN terms was presented in [39].
Close interaction and merger of black holes require numerical relativistic simulations. The first complete simulations, using general-relativistic numerical evolutions of three black holes, were presented in [40,41]. These recent simulations showed that three compact object dynamics display a qualitatively different behavior than Newtonian dynamics. In [42], the sensitivity of fully relativistic evolutions of three and four black holes to changes in the initial data was examined, where the examples for three black holes are some of the simpler cases already discussed in [40,41]. The apparent horizon and the event horizon of multiple black holes have been studied in [43][44][45]. Although fully general-relativistic simulations are available, they are constrained by the number of orbits and separations between black holes.
The paper is organized as follows: in Sec. II, we summarize the equation of motion up to 2.5 post-Newtonian arXiv:1108.4485v2 [gr-qc] 3 Feb 2012 approximation for three spinning bodies. This is followed by a discussion of the chaos indicators that we used to characterize the triple system. In Sec. III A, we describe the numerical techniques used to solve the equation of motion, and we present some results for test cases. The perturbation of a binary system by a third object is presented in Sec. III B, where we performed numerical experiments in order to study the stability and chaos of the triple system. The conclusions are presented in Sec. IV.

A. Notation and units
We employed the following notation: x = (x i ) denotes a point in the three-dimensional Euclidean space R 3 and letters a, b, . . . are the particle labels. We defined r a := x − x a , r a := | r a |,n a := r a /r a ; for a = b, r ab := x a − x b , r ab := | r ab | andn ab := r ab /r ab ; here, | · | denotes the length of a vector. The mass parameter of the a-th particle is denoted by m a with M = a m a . Summation runs from 1 to 3. The linear momentum vector is denoted by p a . A dot over a symbol,˙ x, means the total time derivative, and the partial differentiation with respect to x i is denoted by ∂ i .
In order to simplify the calculations, it is useful to define dimensionless variables (see e.g. [46]). We used as basis quantities for the Newtonian and post-Newtonian calculation the gravitational constant G, the speed of light c and the total mass of the system M . Using derived constants for time τ = M G/c 3 , length l = M G/c 2 , linear momentum P = M c, spin S = M 2 G/c and energy E = M c 2 , we construct dimensionless variables. The physical variables are related to the dimensionless variables by means of scaling. Denoting with capital letters the physical variables with the standard dimensions and with lowercase the dimensionless variables. We defined for a particle a its position x a := X a /l , linear momentum p a := P a /P and mass m a = M a /M (notice that m a < 1, ∀a).

A. Equations of motion
In the ADM post-Newtonian approach, it is possible to split the orbital and spin contribution to the Hamiltonian (see e.g. [47,48]) where each component forms a series with coefficients which are inverse powers of the speed of light. We included terms up to 2.5 PN contributions where the orbital part for three bodies is given by [34,39,49] Here each term of the Hamiltonian is labeled by a superscript n that denotes the PN order (powers of c −2n ) and a subscript which distinguished between the Newtonian terms and the PN components. The spin contribution forms a power series where the leading order is given by [50][51][52][53] (see [31] for the next-to-leading order terms) In this case, the superscripts denote the leading order terms (LO) while the subscript distinguishes the kind of interaction: spin-orbit (SO), spin(a)-spin(b) (SS) or spin(a)-spin(a) (S 2 ). We specified the initial spin of the particles by using a dimensionless parameter χ a ∈ [0, 1] and the two spherical angles θ a and φ a . The initial spin of the particles is given by s a (0) = χ a m 2 a (cos φ a sin θ ax + sin φ a sin θ aŷ + cos θ aẑ ), (4) wherex,ŷ andẑ are the unitary basis vectors in Cartesian coordinates. Using the Hamiltonian (1), the equations of motion areẋ The first term in (2) is the Hamiltonian for n-particles interacting under Newtonian gravity The explicit form of the PN terms in (2) and the equation of motion for three compact objects can be found in [39] (see also [34,49]). The spinning part terms (3) for n compact objects are given in [31]. We will refer as Newtonian, 1 PN, 2 PN and 2.5 PN to the equations of motion derived from H

B. Chaos indicators
There are a number of methods for diagnosing chaos in a dynamical system: the Poincaré surface (for system with less than 3 degrees of freedom), the Kolmogorov-Sinai entropy, the characteristic Lyapunov exponent and the fractal basin boundary method, among others (see e.g. [54][55][56]). For our analysis, we employed the fractal basin boundary method and the Lyapunov exponent.
In this section, we describe the implementation of both methods.
The basin boundary method and the Lyapunov exponents are able to characterize different types of chaos. The basin boundary characterizes the sensitive dependence on initial conditions for nearby orbits with different attractors on the phase-space 1 . However, it does not provide information about the orbit itself. Alternatively, the Lyapunov indicator characterizes the regularity or chaos of an orbit in a way that is independent from the attractor.

Fractal basin boundary
A dynamical system usually exhibits a transient behavior followed by an asymptotic regime. The timeasymptotic behavior defines attracting sets in the phasespace which are called attractors. The attractor can be periodic, quasi-periodic or chaotic. The closure of the set of initial conditions, which has a particular attractor, is called basin of attraction. Basin boundaries for a typical dynamical system can be either smooth or fractal [56,57]. A system with fractal basin boundaries is sensitive to initial uncertainty. Given > 0, we call the initial condition point P * unsafe if there is another initial condition point P inside a neighborhood of size centered at P * which converges to a different attractor. The fraction f of unsafe points is related to by where D is the dimension of the phase-space and d f is the dimension of the basin boundary (see [56,58] for a detailed discussion). It is useful to define an uncertainty exponent α = D − d f , which satisfies the relation 0 < α ≤ 1. The value of α quantifies the presence of chaos in a basin boundary, α ≈ 1 for regular boundaries and α ≈ 0 for a chaotic region. There are many ways of measuring the dimension of a fractal [57,59]. We employed the box-dimension which is defined as where A is a nonempty bounded set in the d Edimensional Euclidean space and N δ (A) is the smallest number of sets in a δ-cover of A. A δ-cover of A is a collection of sets of diameter δ whose union contains A. For a physical or numerical experiment, the fractal basin boundary has finite resolution. A graphical representation of the basin boundary is given by an image formed by pixels of size l res . In order to estimate the dimension of the basin boundary, we compute the quotient 1 An attractor is a set towards which a dynamical system asymptotically approaches in the course of time evolution (see e.g. [54]). ln(N δn (A))/ ln(1/δ n ) for a set of δ n ≥ l res . The typical behavior of the resulting data is well approximated by a straight line. A linear regression method is used to fit a linear function to the resulting data (where δ n = nl res and n ∈ {1, 2, . . . , 12}). Therefore, the value of the fitted function for δ = 0 gives an estimate value of the dimension. We tested the method by computing the dimension of non-fractal one-dimensional sets (i.e., a regular basin boundary) and Sierpinski's gasket and carpet.

Lyapunov exponent
Lyapunov exponents are another way for characterizing chaos in a dynamical system (see e.g. [54]). For an autonomous n-dimensional systeṁ a given reference solution X * and a nearby solution X = X * + δ X. While the evolution of the difference δ X is given by the linearized equation where J( X * ) i,j := ∂ j F i ( X * ) is the Jacobian matrix of F evaluated at the reference solution X * . The solution of (12) is given by [60] δ X(t) = e Jt · δ X(0).
The Lyapunov exponents are defined by the eigenvalues Λ i (t) of the distortion matrix Λ := e (J+J T )t Sensitive dependence on initial conditions exists when λ i > 0. Therefore, in order to distinguish between regular and chaotic orbits, it is sufficient to compute the principal Lyapunov exponent λ p := max(λ i ). Calculating λ p directly is computationally expensive. However, by using the evolution of the difference, δ X, it is possible to obtain an estimation of λ p where, and δX(t) = |δ X|. We will refer as a Lyapunov indicator to λ * p in order to distinguish it from the principal Lyapunov exponent and the Lyapunov function, λ * p (t, δX(0)). The norm, |δ X|, plays an important role in the computation of the Lyapunov indicator. Particularly, in the case of general relativity, additional gauge effects can affect the estimation of the Lyapunov indicator [23,55,61]. Results from literature were applied to general relativistic dynamic of test particles and to binary systems in the center of mass frame. Additionally, the estimation of λ p is often computed using the twoparticle method instead of the evolution of the difference (see e.g. [62]). However, we did not noticed significant differences in the computation of λ * p given the gauge effects. We present the results in terms of the Cartesian length of δ X.
In practice, the Lyapunov indicator is computed approximating the limits, using a sufficient small δ X(0) and sufficient large integration time. In our simulations, we used |δ X(0)| = 10 −9 . The integration time depends on the system (see Sec. III). The evolution is performed by solving numerically the system of equations (5)- (7), therefore X = ( x a , p a , s a ) T is a vector with 27 components. The right hand side of (11) is the vector and the equation of motion for the difference δ X = (δ x a , δ p a , δ s a ) is where we define the auxiliary operator

A. Numerical methods
We solved the equations of motion numerically using the GNU Scientific Library (GSL) [63] and the Olliptic code infrastructure [42]. We generated the right hand side (RHS) of the equations of motion (5)-(7) with The simulations were done using the embedded Runge-Kutta Prince-Dormand (8,9) method provided by the GSL. We used a scaled adaptive step-size control (see [63] for details). In our simulations, the error control for the dynamical variables ( x a , p a , s a ) was set to 10 −11 and to 10 −6 for the variation variables (δ x a , δ p a , δ s a ). We used as an indicator of accuracy the estimate of the local error provided by the GSL routines and the conservation of the Hamiltonian for the simulations which do not included the radiation term.
For our stability analysis, we solved the system of equations (5)- (7) and (17)- (19) simultaneously. In order to perform the simulations in an efficient way, we adapted the number of variables depending on the problem. The parameters to consider are the number of bodies nb, the magnitude of the spin s a and whether we are solving a planar or non-planar problem. As we mentioned before, the RHS of Eqs. (5)-(7) were generated by a Mathematica script. For the Newtonian case and 1 PN cases, it is convenient to compute the analytical expression for the RHS of Eqs. (18)- (20). The analytical expression for the 2 PN and the spinning cases produces large computational source files. The source file for the 2 PN is of 25 megabytes while the inclusion of the spin generates files over a hundred megabytes. The compilation and optimization of large files is not practical, since compilers run out of memory even in the 16 GB of RAM server that we tried. Therefore, for the 2 PN order and the spinning case we computed the RHS of Eqs. (18)- (20) numerically. In order to estimate the errors in the numerical computation of the Jacobian, we compared the results obtained for the Newtonian and 1 PN cases employing both methods (see Sec. III A 1). The system (18)-(20) is the Jacobian matrix of (17). The problem reduces to the computation of A simple 2nd, 4th and even 6th order finite difference method with a fix step size h does not produce accurate solutions. The main issue is the nature of the components, for example, the position variables require a different optimal step-size than the momentum or the spin variables. We implemented an adaptive method based on the GSL differentiation routine. The method produces a solution with an estimated error below 10 −8 (see Sec. III A 1). An important issue in the numerical integration of a three-body system arises when two of the bodies are very close to each other. In the case of adaptive step size methods, it is necessary to reduce the step size in order to resolve properly the orbits in the close interaction phase. A usual approach to deal with this problem is to perform a regularization of the equations of motion, see e.g., [65][66][67][68][69] and references therein. However, in our simulations, we included a different criteria. Regardless of the equation of motion employed, we monitored the absolute value of each conservative part of the Hamiltonian (1) relative to the sum of the absolute values .
(23) The simulations stop when the contribution of the first post-Newtonian correction is larger than 10%. We will refer to H % 1 > 10% as strong interaction. Empirical results showed that, for a similar configuration, the resulted waveform exhibits a growth characteristic of the merger phase [39]. The dynamics after a strong interaction may lead to the merger of two of the bodies or to the escape of the lighter body. In order to determine safely whether there is a merger or not, it is necessary to employ full numerical relativistic simulation.
Additionally, we considered the Newtonian escape energy of each body in respect of the other two components where The simulation stops when E a esc > 0 for one of the bodies and the size of the system S := a | x a | is 100 times the initial size. During a close encounter of two of the bodies it is common to have a positive value of E a esc for a short time. However, imposing the second condition we exclude those cases. We performed several tests to estimate the numerical errors. Here we summarize the results of these tests.

Numerical tests
Our first test was done with Lagrange's triangle solution using Newtonian equations of motion. In Lagrange's solution each body is sitting in one corner of an equilateral triangle (see e.g. [70]). We set the sides of such triangle to r 12 = r 23 = r 31 = 100, the mass ratio to 1:10:100, and the eccentricity to zero. Then each body follows a circular orbit (with different radii) around the center of mass. This configuration is unstable and after 14 periods one of the bodies escapes from the system. The unstable property is reflected in the Lyapunov indicator λ * p . Fig. 1 (a) shows log 10 (λ * p ) computed using two methods for the prescription of the RHS of Eqs. (17)- (19). The first method is the analytic expression (labeled J a ) and the second, the numerical solution (labeled J n ). The two methods produce solutions that cannot be distinguishable within the plot. The maximum relative difference is 4 × 10 −6 . The Lyapunov indicator exhibits the characteristic behavior of a chaotic system. After t = 10 4 , λ * p (t) oscillates around a constant value. Fig. 1 (b) shows the relative variation of the Hamiltonian for the same solution. The conservation of the Hamiltonian has a similar behavior either using J a or J n . The Hamiltonian exhibits a variation of around 10 −12 before the evolution stops. Fig. 1 (c) illustrates the L 2 norm of the estimated errors for Lagrange's triangle solution. By looking at the errors, the difference between the solution produced by J n and J a becomes clear. The numerical solution generates errors which are two orders of magnitude larger than the analytical Jacobian. However, even when using J a , the accuracy is good, with errors below 1.5 × 10 −8 . A second test was done for a stable system. Using the 1 PN equations of motion, we computed the Lyapunov indicator λ * p , the relative variation of the Hamiltonian ∆H and the L 2 norm of the estimated errors for the Hénon's criss-cross solution [32,71,72]. We used initial parameters given by wherex,ŷ andẑ are the unitary basis vectors in Cartesian coordinates, and λ is a scaling factor (for our simulation λ = 10). Notice that for this test we used the parameters given in [73] with the scaling factor λ, and doing a change of variables from initial velocity to initial momentum. Therefore, we are not including post-Newtonian corrections to the initial parameters. As in the previous case, λ * p does not show significant differences when computed using J a or J n (see Fig. 2). The value of λ * p decreases monotonically almost like a straight line; this is the characteristic behavior of a regular solution. The maximum relative difference in λ * p using J a and J n is 8×10 −6 . The Hamiltonian is conserved with a maximum variation of 5 × 10 −11 . The L 2 norm of the estimated errors shows noisy errors using J n while J a exhibits a flat trend. A third test was done with a known chaotic binary system. We used a chaotic eccentric configuration with mass ratio 3:2 and initial parameters: (26) A similar system, with a different unit convention, was studied in [74] (Sec. IV-C1). Following [74], we solved the equation of motion using up to 2 PN correction for the orbital part and leading order in the spin. As was emphasized in [74], in this configuration, the two compact bodies are rather close to each other. In order to obtain an accurate physical description, it is necessary to include, for the orbital part, corrections higher than 2 PN and, for the spinning part, higher than the leading order. The 1 PN contribution to the Hamiltonian (22) during a close encounter of the bodies can reach 50%. At this point it is important to recall one of the conclusions of [74]. When using up to 2 PN Hamiltonian with leading order spin contributions, there is no evidence of chaotic systems for configurations which are physically valid to the PN order considered. Therefore, for comparison, this test was performed without using the stopping criteria H % 1 > 10%. The result of the third test is presented in Fig. 3. The Spinning binaries with an eccentric orbit. Lyapunov function λ * p (t, δX(0)) (a), conservation of the Hamiltonian ∆H (b) and the L2 norm of the estimated errors, for the system with initial parameters given by (26) (solid line) and the same configuration except for s1(0) = 0 (dashed line).
system was evolved using the parameters given in (26), and for reference, a non-chaotic orbit with the same initial condition but with s 1 (0) = 0 was used. Fig. 3-(a) shows the Lyapunov indicator for the chaotic configuration (solid line) and the arithmetic mean computed for t > 10 6 . The estimated value of the Lyapunov indicator is λ * p = 10 −5 (dotted line) with a standard deviation of 1.3 × 10 −6 , corresponding to a relative variation of 13%. In contrast, the regular solution shows that λ * p decreases almost monotonically after t = 10 4 . In both simulations, the Hamiltonian is numerically conserved below 10 −8 (see Fig. 3-(b)), and the estimated errors are of order 10 −6 for the chaotic solution and 10 −10 for the regular solution. The main source for the errors, in the case of the chaotic solution, is the numerical calculation of the Jacobian matrix.
We performed, additional tests already presented in Sec. III-A of [39]. A direct comparison with the results in [39] does not showed significant differences in the results. For conciseness, we do not display the results of the tests.

B. Stability of hierarchical systems
Here, we consider the strong perturbation of a binary compact object system due to a third smaller compact object. As a basic configuration, we studied a Jacobian system with mass ratio 10:20:1. The inner binary system had initial apo-apsis ( i.e. maximum separation of a Hierarchical system (see [39]). Initial configuration of the inner and external binaries. The initial momentum of the third body is given by considering the external binary as a Newtonian binary. Shown are the osculating orbital planes Πin and Πext for inner and external binary orbits. The two planes are inclined by an angle ι.
Keplerian orbit) r b (0) = 130 and eccentricity e b (0) = 0. Although there are methods in post-Newtonian dynamics to specify the initial parameters of a binary system with a given eccentricity (see e.g. [74][75][76][77][78][79]), in this study the inner binary is strongly-perturbed by a third body. Therefore, it is necessary to account for additional effects. For simplicity, we set the initial parameters considering only the Newtonian dynamics of a non-perturbed binary, where the eccentricity refers to the Newtonian case. In this approach, we view the third compact body and the center of mass of the inner binary as a new binary (we will refer to it as the external binary). The bodies start from a configuration where the initial radial vector r 12 is perpendicular to initial vector position x 3 of the external body (see Fig. 4). We denote the inclination angle between the osculating orbital planes Π in and Π ext by ι (see Fig. 4). To be more specific, the initial position and momentum of each body is given by Notice that we set the momentum in terms of the eccentricity and the apo-apsis. Therefore, a higher eccentricity given to the external binary implies a stronger perturbation of the inner binary.
The main goal of this study was to characterize the stability and chaos of a hierarchical system as function of the initial apo-apsis r 3 , the eccentricity e 3 for the external binary and the inclination angle ι between the osculating orbital planes. We explored the influence of the post-Newtonian corrections.

Final state survey
An analysis of the asymptotic behavior of the system as function of three parameters r 3 , e 3 and ι was performed. For a given osculating angle ι, we produced a map which is a subset of the space of configurations (r 3 , e 3 ). We define three possible outcomes which characterize the asymptotic behavior of the system: the escape of one of the bodies, a strong interaction H % 1 > 10% and a system which remains stable up to t f = 2 × 10 7 . We assign a color to each outcome in the map of initial configurations (r 3 , e 3 ): dark gray (color online) for an escape, white for a strong interaction and black for a stable configuration. Strictly speaking, the result is not a basin boundary map since the outcomes are not attractors in the phase space and the parameters (r 3 , e 3 ) are not coordinates of the phase space. However, Eqs. (27) gives the link between (r 3 , e 3 ) and the coordinates of the phase space. On the other hand, a system where one of the bodies escapes will be attracted to a hyper-plane in the phase space with one of the spatial coordinates at infinity; a strong interaction may become an escape or the merger of two of the bodies. A merger can be represented by a hyper-plane where two of the bodies share the same values of position and momentum. For a stable system, it is hard to define a single attractor because it can exhibit a different asymptotic behavior for t > t f . For simplicity, we refer to the boundary between the various color in the map as basin boundary since it can capture some of the features of a basin boundary.
The choice of t f = 2 × 10 7 is motivated by the dynamics of a similar Jacobian system presented in [39]. Including up to 2.5 PN terms, the inner binary arrives to the merger phase after t mgr ≈ 2.7 × 10 7 . The assumption is that if a conservative system remains stable until that time, then the corresponding radiative system will remain stable until the merger of the inner binary. We selected a shorter time (i.e., t f < t mgr ) due to computational costs. However, a set of numerical experiments for different values of t f suggest that a non-stable orbit is manifest before t f .
We explored the parameter space r 3 , e 3 using 4 regions defined by The domain of each region is the Cartesian product and its resolution is denoted by ∆r 3 , ∆e 3 . Each map was done computing 300 × 400 orbits except for region R 3 that contains 200 × 400 orbits. Fig. 5 shows the result for a Newtonian simulation and osculating angle i = 0.  Notice that R 3 ⊂ R 2 ⊂ R 1 ⊂ R 0 , i.e. each R i is a magnification of the region R i−1 (in the figure, the sub-regions are indicated by rectangles). Fig. 5 shows the characteristic behavior of a fractal set; every magnification reveals a more complex structure from which it is possible to distinguish some degree of self-similarity. Table I summarizes the quantitative results for the basin boundary. Notice that the box-dimension reflects the fact that the sets are not perfectly self-similar. We choose R 1 and R 2 to perform additional simulations. A comparison using Newtonian and 1 PN orbits for regions R 1 and R 2 is presented in Fig. 6. There is a clear difference in the set of stable orbits. In the Newtonian case, there are two disjoint sets where the 1 PN system exhibits a single set (compare the black point in panels (a) and (b) of Fig. 6). The magnification R 2 shows differences in the escape and strong interaction zones. The  Fig. 6). Our comparison includes a set of simulations for leading order spinning 1 PN particles. We employed maximally spinning particles χ a = 1 and θ a = 0, φ a = 0 except by θ 1 = π/2, θ 2 = 2π/2 (see Eq. (4)). Therefore, the initial spin of the particles (configuration-a) is given by

PN case exhibits slightly different substructures (see panels (c) and (d) of
Opposite to the non-spinning case, the orbits for our configuration of spinning particles are not coplanar. However, the change in the final states is negligible, the resulting map is almost indistinguishable to the 1 PN case (for brevity, we do not display the basin boundaries which are similar to panels (b) and (d) of Fig. 6). The quantitative analysis of the six basin boundaries under comparison is displayed in Table II. The difference on the substructures between the Newtonian case and the 1 PN case is reflected by the box-dimension (smaller value for the 1 PN case). In the region R 1 , there is a small change in the distribution of 'Regular', 'Escape' and strong interactions between Newtonian and 1 PN. Nevertheless, in region R 2 there are noticeable differences in the percentages of 'Escape' and strong interactions. The similarity between the spinning and the non-spinning 1 PN case appears for the dimension and the distribution of the three outcomes (with minor differences in the percentages). In order to investigate the influence of the gravitational radiation, we performed an evolution which contains the 2.5 PN terms. The inclusion of PN terms is computationally expensive. Therefore, we used a radiative 1 PN Hamiltonian (i.e., we include the Newtonian, 1 PN and 2.5 PN terms in (2)). For similar configurations, the difference in the dynamics between the full 2.5 PN system and the radiative 1 PN is relatively small [39]. Fig. 7 shows the fractal basin boundary and Table II the corresponding quantitative measures. The radiative term does not seem to produce a significant change in respect of the inclusion of 1 PN terms (compare Fig. 7 with Fig. 6-(b)). Particularly, the box-dimension of the boundary is the same, and there is only a small difference in the distribution of the percentages of escape and strong interactions.
The osculating angle ι has a strong influence in the region of the initial parameters R 1 . For the Newtonian case, we produced a set of eight additional maps for region R 1 where the osculating angle takes the values FIG. 8. Fractal basin boundary of region R1 as function of the osculating angle ι (Newtonian). The meaning of the color is the same as in previous plots (see e.g. Fig. 6). From the top to the bottom and from the left to the right ι takes the values 0, π/16, π/8, 3π/16, π/4, 5π/16, 3π/8, 7π/16, π/2.  ι (1 PN). The meaning of the color is the same as in previous plots (see e.g. Fig. 6). From the top to the bottom and from the left to the right ι takes the values 0, π/16, π/8, 3π/16, π/4, 5π/16, 3π/8, 7π/16, π/2. ι = nπ/16, n ∈ {1, 2, . . . , 8}. The results, including the case ι = 0, are presented in Fig. 8. An analogous set of simulations was produced using the 1 PN equation of motion. Fig. 9 shows the corresponding basin boundaries. Table III shows the quantitative results for both sets. Notice that for i > π/4 the percentage of stable points considerably increases for the 1 PN case in respect of the Newtonian. The box-dimension in both cases has a growing-oscillatory behavior.
The slope of the linear part of the 1 PN is 1.6 times larger than the Newtonian one. However, the oscillatory part is nearly 0.7 times the 1 PN value for the Newtonian case. This result shows that the chaotic properties of this hierarchical configuration increases in both cases, reaching the maximum at ι = π/2. From the basin boundary figures, we can notice an increasing number of unsafe points. There is an evident difference for ι = π/2 between the Newtonian and 1 PN basin boundary (comparing the bottom-right panel of Figs. 8 and 9). We analyze this particular case in the next section.

Analysis of orbits
We used the basin boundaries as a guide for a detailed study of specific orbits. Given the region R 1 , an osculating angle ι, two different basin boundaries Ω(r 3 , e 3 ) and Ω * (r 3 , e 3 ), it is possible to produce a new map D that shows the differences between Ω(r 3 , e 3 ) and Ω * (r 3 , e 3 ). We constructed the map as follows; D(r 3 , e 3 ) takes value 1 if Ω(r 3 , e 3 ) = Ω * (r * 3 , e * 3 ) in a neighborhood of (r 3 , e 3 ) of size at least 6dr 3 × 6de 3 , i.e., for r * 3 = r 3 + idr 3 , e * 3 = e 3 + ide 3 and i ∈ {−3, −2, . . . , 3}. If the previous condition is not satisfied, then the D(r 3 , e 3 ) is set to 0. D(r 3 , e 3 ) provides a way for identifying zones where the two basin boundaries differed in a relatively large region. From the analysis of this map and the basin boundaries, we explored several initial parameters where we compared the orbits for different configurations of the PN equations of motion. In the following, we describe 5 representative comparisons. For the region R 1 and ι = π/2 there is a clear difference between the Newtonian and 1 PN evolution (i.e., the difference between the bottom-right panel of Figs. 8 and 9). The bottom-right zone of parameters leads to a different outcome. In the Newtonian case, the simulation with parameters e 3 < (r 3 − 310)/320 ends when H 1PN > 10%. However, for the same set of parameters the 1 PN case remains stable. For the Newtonian case, there is a fast growth on the eccentricity of the inner binary induced by the external body (see Fig. 11-(a)). This effect is produced by a resonance between the orbital dynamic of the inner and external binary (see [80] for a detailed description). For the 1 PN case the orbits remain stable during the simulation (see Fig. 11-(b)). In both cases, the evolution of the Lyapunov indicator does not showed evidence of exponential growth behavior (see Fig. 11-(c)). Notice that in the Newtonian case, the strong interaction is between the components of the inner binary. The behavior of the Lyapunov function suggests that the strong interaction between the heavy components of the system does not follow chaotic trajectories.
For the region R 1 and ι = 0, there is a significant difference between the Newtonian and the 1 PN basin boundaries for stable orbits (i.e., the black dots in panels (a) and (b) of Fig. 6). The Newtonian configuration remains stable without noticeable instability. On the other hand, the coupling between the orbits in the 1 PN case results in the escape of the third body (see panels (a) and (b) of Fig. 12). The Lyapunov function shows for 1 PN dynamics the characteristic behavior of escape orbits. During a close encounter between the lighter body and one of the heavy components of the inner binary, the difference vector has an exponential growth which is followed by a regular behavior. The regular behavior is expected for an uncoupled binary-single body system (see the kink in the Lyapunov function of Fig. 12-(c)). The exponential growth for the difference vector is an indication of sensitivity to initial conditions. In general, the resulting escape orbit is chaotic; a small change in the initial condition produces a significant change in the escape direction.
A comparison between 1 PN and 2 PN dynamics showed that almost identical orbits are produced by an evolution where the lighter body has a quick encounter with the inner binary which is followed by an escape or a strong interaction. Fig. 13 displays a comparison between the Newtonian, 1 PN and 2 PN evolutions for 11 nearby configurations. The initial parameters are ι = 0 and e 3 = 0.1 with r 3 ∈ {269.95, 269.96, . . . , 270.05}. The relative change in the initial separation r 3 between the configurations is ∼ 0.004%. The Newtonian evolution (Fig. 13-(a)), exhibits a chaotic behavior. Every simulation ends with the escape of the lighter body. There are significant differences in the trajectories. On the other hand, the 1 PN and 2 PN evolutions (panels (b) and (c) of Fig. 13) produce a quick escape which are similar for each initial condition. The maximum final difference between the orbits is 1.3% for the 1 PN case and 1.5% for the 2 PN simulations. Moreover, the final difference  between the 1 PN and 2 PN reference evolution r 3 = 270 is 6%. However, for successive encounters the evolution is significantly different. Fig. 14 shows the resulted orbits for initial parameters ι = 0, r 3 = 240, e 3 = 0. For these initial parameters, the resulting outcome is the escape of the lighter body. Nevertheless, in the 1 PN evolution there is an ejection 2 of the lighter body before the final escape. In contrast, the ejection is not present in the 2 PN orbit (compare panels (a) and (b) of Fig. 14). There is a significant difference in the evolution of the Lyapunov function. For the 1 PN case, there is a double kink which corresponds to the close encounters before and after the ejection. Nevertheless, the 2 PN evolution presents a single long kink (see Fig. 14-(c)). Fig. 15 shows a comparison between the Newtonian, 1 PN and 2 PN evolutions for 11 nearby configurations. The initial parameters are ι = 0 and e 3 = 0 with r 3 ∈ {279.95, 279.96, . . . , 280.05}. Similarly to the case presented in Fig. 13, the relative change in the initial separation r 3 between the configurations is ∼ 0.004%. The set of configurations includes the solution presented in Fig. 14 r 3 = 280. The evolutions exhibit sensitive dependence on initial conditions. In the three cases, Newtonian, 1 PN and 2 PN, relatively nearby initial parameters 2 In the three-body problem literature, an ejection refers to a temporal large separation of one of the bodies which finally returns [21]. produce a significant change in the orbits. Moreover, the inclusion of PN corrections produced a noticeable change in the trajectories. For the parameters that we have studied, the inclusion of spin correction has a small influence in the final outcome. However, like in the 2 PN case, for some of the initial parameters, the evolution is completely different. For the initial parameters ι = 0, r 3 = 320 and e 3 = 0.39, the 1 PN orbit ends with the strong interaction between particles 1 and 3 (see Fig. 16-(a)). In contrast, the spinning leading order 1 PN for the spin configuration-a (28) results with the escape of particle 3 (see Fig. 16-(b)). For the escaping orbit, the evolution produces a Lyapunov function which after t = 10 5 decreases linearly. The spinning case shows initially the same trend. However, the strong interaction produces a short kink at the end of the evolution (see Fig. 16-(c)).
We explored a different set of spin parameters. Here we present one of the results of configuration-b given by Fig. 17 shows the result for the initial parameters ι = π/4, r 3 = 242, e 3 = 0.2 for the 1 PN and spinning leading order 1 PN cases. In this evolution, the non-spinning case results in the escape of particle 3. For the spinning particles, the result is the strong interaction between particles   is similar to (a) but for the leading order spinning 1 PN configuration-a (initial spin given by (28)). The lower panel (c) shows the evolution of the Lyapunov indicator for both cases. 1 and 3. The Lyapunov function shows the characteristic exponential growth of the difference vector during the close interaction. Table IV shows the evolution of the eccentricity of the inner binary e b . The eccentricity is computed by calculating numerically the value of the apo-apsis and peri-apsis of the inner binary i.e., computing a sequence of local maxima and minima of r 12 . The eccentricity is given by Notice that the initial value of e b in every case is relatively small but not zero due to the perturbation of the third body. However, the final value of the eccentricity depends on the final state of the evolution. For regular orbits, the change is relatively small (see the standard deviation of rows 2 and 3 in Table IV). On the other hand, the final eccentricity of the escaping orbits is relatively large during the whole evolution (compare the final eccentricity and the arithmetic mean of the corresponding rows). The orbit with resonance exhibits the largest increment in eccentricity. The final value before the evolution stops is close to 0.5. For PN evolutions, we did not find evolutions with this kind of resonance. However, a different set of parameters may lead to similar behavior. For a strong interaction between two of the bodies, the resulting final eccentricity is not meaningful since there are two possible outcomes in the full evolution scenario. One possibility is an ejection or escape of the lighter body leading to a increment in the eccentricity of the inner binary similar to the cases presented. The other possibility is a successive merger of the three bodies.

IV. DISCUSSION
We performed a numerical study of stability and chaos of a hierarchical three compact object configuration. The configuration is composed by a binary system which is perturbed by a third, lighter body positioned initially far away from the binary. For a given inner binary, we explored the influence of the third body taking as free parameters the angle of the osculating orbital planes, the initial apo-apsis and eccentricity of the third body. Using the basin boundary method, a total of 2,960,000 simulations were analyzed. From the total, 47.3% corresponded to Newtonian evolutions, 40.54% to (0+1) PN simulations, 4.05% to (0+1+2.5) PN evolutions and 8.11% are simulations of spinning particles including leading order in spin and 1 PN terms in the orbital part.
A total of twenty-five basin boundary maps were produced. The basin boundaries exhibit the characteristics of fractal sets: some degree of self-similarity and an increasing complexity under magnification. Fractal basin boundaries are an indicator of chaos in dynamical systems. By measuring the fractal dimension of the basin boundaries it is possible to determine the uncertainty exponent α. The property of the exponent is that for α ≈ 0 the dynamical system is chaotic and for α = 1 the system is regular. The values of α for the fractal basin boundaries considered here are between 0.02 and 0.26. The 1 PN set of data produces an exponent slightly larger than the Newtonian case (the maximum difference is 1.6%). On the other hand, the osculating angle ι has a strong influence in the uncertainty exponent. The exponent decreases as a linear oscillatory function in the range [0, π/2]. The difference between the planar case (ι = 0) and the case where the third body starts from a direction perpendicular to the orbital plane of the inner binary (ι = π/2) is 7.1% for the Newtonian simulations and 8.1% for the 1 PN case.
In addition to the uncertainty exponent, we quantified the percentage of stable orbits, escapes and strong interactions. The configurations selected contain between 0% and 9.3% of stable orbits. The distribution of escapes are between 40.3% and 77.6%. The strong interactions oscillate between 20.3% and 82%. A remarkable case was found at ι = π/2 where the Newtonian simulation presents a resonance which drops the number of stable orbits to 0.1%. Opposite to the Newtonian case, a dynamic including 1 PN or higher corrections eliminates the resonance. The number of stable orbits increases to 9.3%. The presence of resonances in the three-body problem is fundamental for understanding the chaotic properties of the system, see e.g. [18,21,80] and references therein.
By looking at specific orbits, it is possible to notice some of the different couplings between the inner and external binaries. For most of the orbits, the 1 PN terms have the strongest effect on the dynamics. The inclusion of the spinning particles or 2 PN corrections produces small differences in the orbits when the lighter body has a quick encounter with the inner binary which is followed by an escape or a strong interaction. In that cases the orbits are essentially identical. However, the coupling of the orbits for some of the configurations can produce a significant change in the final outcome. We presented five representative examples.
The spin can produce a significant change in the orbits producing precession of the orbital planes. However, the change on the dynamics is not strong enough to modify the final outcome. On the other hand, gravitational radiation, in general, produces a small change in the energy to be significant in the short time scale of the evolution. Similarly to the spinning case, gravitational radiation has a small effect in the asymptotic behavior.
Considering the scenario of a binary system in a galactic core or a region with high density of compact objects, we expect to have the following results. Between a 40% and 70% probability that the lighter body escapes the system depending on the relative orbital planes. The remanent binary may result in eccentricities between 0.1 and 0.2. Between 23% and 58% of the bodies may lead to a merger with one of the components of the inner binary; contributing to the growth of the compact objects. The probability of finding a body in a stable orbit is between 1% and 10%, producing a perturbation which leads to small eccentricities of the inner binary of around 0.03. Nevertheless, the combined effect of many bodies of comparable size can significantly change the results.
More general statements about the chaotic behavior of three compact bodies are possible but require an extensive parameter study. Other configurations include, for example, a characterization of the stability and chaos based on the mass ratios, inner binary separation, magnitude and direction of the spin or initial eccentricity of the inner binary. We consider this a topic for future study.