Stability and chaos of hierarchical three black hole configurations
SStability and chaos of hierarchical three black hole configurations
Pablo Galaviz
1, 2, ∗ School of Mathematical Science, Monash University, Melbourne, VIC 3800, Australia Theoretical Physics Institute, University of Jena, 07743 Jena, Germany (Dated: April 16, 2019)We study the stability and chaos of three compact objects using post-Newtonian (PN) equationsof motion derived from the Arnowitt-Deser-Misner-Hamiltonian formulation. We include termsup to 2.5 PN order in the orbital part and the leading order in spin corrections. We performednumerical simulations of a hierarchical configuration of three compact bodies in which a binarysystem is perturbed by a third, lighter body initially positioned far away from the binary. Therelative importance of the different PN orders is examined. The basin boundary method and thecomputation of Lyapunov exponent were employed to analyze the stability and chaotic propertiesof the system. The 1 PN terms produced a small but noticeable change in the stability regions ofthe parameters considered. The inclusion of spin or gravitational radiation does not produced asignificant change with respect to the inclusion of the 1 PN terms.
PACS numbers: 04.25.Nx, 04.70.Bw, 05.45.Pq, 05.45.Jn
I. INTRODUCTION
The dynamics of compact objects play an importantrole in the evolution of galaxies and other stellar sys-tems. Post-Newtonian (PN) techniques are a useful toolfor modeling the dynamics of multiple compact objects.In astrophysical models, the gravitational radiation is in-cluded via an effective force which considers 1 PN and2.5 PN corrections, and in some cases stellar dynamicalfriction. Hierarchical three black hole configurations in-teracting in a galactic core have been the focus of recentstudies by several authors. For example, intermediate-mass black holes with different mass ratios are consid-ered in [1–3]; simulations of dynamical evolution of tripleequal-mass supermassive black holes in a galactic nucleiwere performed in [4, 5]. Other astrophysical applicationsof multiple black hole simulations include three-bodykicks [6, 7] and binary-binary encounters (see e.g. [8–12]).In the context of the final parsec problem [13, 14], three-body interactions are considered as a mechanism that candrive a binary black holes system to a separation belowone parsec.Since the 1990s there has been an increasing effort todetect and study compact objects. The most likely sourceof gravitational waves are binary compact objects. Re-cently, it was shown that the probability of more than twoblack holes to interact in the strongly relativistic regimeis, not surprisingly, very small [15]. For practical pur-poses, the creation of gravitational waveform templatesfor gravitational wave detectors is naturally focused onbinary systems. Binary systems can produced complicatewaveforms when taking into account spinning black holesand eccentric orbits, e.g. [16]. On the other hand, triplesystems consisting of a binary black hole and a star (e.g.,a white dwarf) are potential sources of electromagneticcounterparts associated with the mergers [17–20]. ∗ [email protected]
The chaotic behavior of triple systems is well known inthe Newtonian case (see [21] and references therein). Forbinaries, it is known that chaos appears when using cer-tain post-Newtonian approximations for spinning binarysystems (see [22–28]). In the present work, we studiedthree-body systems with PN methods, where the maintechnical novelty is the inclusion of the 2.5 PN terms inthe orbital dynamics and leading order of spin-orbit andspin-spin terms [29–31].Recently, similar PN techniques were applied to thegeneral relativistic three-body problem. Periodic solu-tions were studied using the 1 PN and 2 PN approxima-tions in [32–34]. Examples of three compact bodies ina collinear configuration were considered in [35, 36], andLagrange’s equilateral triangular solution was studied in-cluding 1 PN effects in [37]. In [38], the stability of theLagrangian points in a black hole binary system was stud-ied in the test particle limit, where the gravitational radi-ation effects were modeled by a drag force. The waveformcharacterization of hierarchical non-spinning three-bodyconfigurations using up to 2.5 PN terms was presentedin [39].Close interaction and merger of black holes require nu-merical relativistic simulations. The first complete simu-lations, using general-relativistic numerical evolutions ofthree black holes, were presented in [40, 41]. These recentsimulations showed that three compact object dynamicsdisplay a qualitatively different behavior than Newtoniandynamics. In [42], the sensitivity of fully relativistic evo-lutions of three and four black holes to changes in theinitial data was examined, where the examples for threeblack holes are some of the simpler cases already dis-cussed in [40, 41]. The apparent horizon and the eventhorizon of multiple black holes have been studied in [43–45]. Although fully general-relativistic simulations areavailable, they are constrained by the number of orbitsand separations between black holes.The paper is organized as follows: in Sec. II, we sum-marize the equation of motion up to 2.5 post-Newtonian a r X i v : . [ g r- q c ] F e b approximation for three spinning bodies. This is followedby a discussion of the chaos indicators that we used tocharacterize the triple system. In Sec. III A, we describethe numerical techniques used to solve the equation ofmotion, and we present some results for test cases. Theperturbation of a binary system by a third object is pre-sented in Sec. III B, where we performed numerical ex-periments in order to study the stability and chaos of thetriple system. The conclusions are presented in Sec. IV. A. Notation and units
We employed the following notation: (cid:126)x = ( x i ) denotesa point in the three-dimensional Euclidean space R andletters a, b, . . . are the particle labels. We defined (cid:126)r a := (cid:126)x − (cid:126)x a , r a := | (cid:126)r a | , ˆ n a := (cid:126)r a /r a ; for a (cid:54) = b , (cid:126)r ab := (cid:126)x a − (cid:126)x b , r ab := | (cid:126)r ab | and ˆ n ab := (cid:126)r ab /r ab ; here, | · | denotesthe length of a vector. The mass parameter of the a -thparticle is denoted by m a with M = (cid:80) a m a . Summationruns from 1 to 3. The linear momentum vector is denotedby (cid:126)p a . A dot over a symbol, ˙ (cid:126)x , means the total timederivative, and the partial differentiation with respect to x i is denoted by ∂ i .In order to simplify the calculations, it is useful to de-fine dimensionless variables (see e.g. [46]). We used as ba-sis quantities for the Newtonian and post-Newtonian cal-culation the gravitational constant G , the speed of light c and the total mass of the system M . Using derivedconstants for time τ = M G/c , length l = M G/c , lin-ear momentum P = M c , spin S = M G/c and energy E = M c , we construct dimensionless variables. Thephysical variables are related to the dimensionless vari-ables by means of scaling. Denoting with capital lettersthe physical variables with the standard dimensions andwith lowercase the dimensionless variables. We definedfor a particle a its position (cid:126)x a := (cid:126)X a / l , linear momentum (cid:126)p a := (cid:126)P a / P and mass m a = M a /M (notice that m a < ∀ a ). II. EVOLUTION METHODA. Equations of motion
In the ADM post-Newtonian approach, it is possible tosplit the orbital and spin contribution to the Hamiltonian(see e.g. [47, 48]) H ( (cid:126)x a , (cid:126)p a , (cid:126)s a ) = H ( (cid:126)x a , (cid:126)p a ) Orb + H ( (cid:126)x a , (cid:126)p a , (cid:126)s a ) Spin , (1)where each component forms a series with coefficientswhich are inverse powers of the speed of light. We in-cluded terms up to 2.5 PN contributions where the or-bital part for three bodies is given by [34, 39, 49] H Orb = H (0)N + H (1)PN + H (2)PN + H (2 . . (2) Here each term of the Hamiltonian is labeled by a super-script n that denotes the PN order (powers of c − n ) anda subscript which distinguished between the Newtonianterms and the PN components. The spin contributionforms a power series where the leading order is given by[50–53] (see [31] for the next-to-leading order terms) H Spin = H LOSO + H LOSS + H LOS . (3)In this case, the superscripts denote the leading orderterms (LO) while the subscript distinguishes the kindof interaction: spin-orbit (SO), spin(a)-spin(b) (SS) orspin(a)-spin(a) (S ). We specified the initial spin of theparticles by using a dimensionless parameter χ a ∈ [0 , θ a and φ a . The initial spinof the particles is given by (cid:126)s a (0) = χ a m a (cos φ a sin θ a ˆ x + sin φ a sin θ a ˆ y + cos θ a ˆ z ) , (4)where ˆ x , ˆ y and ˆ z are the unitary basis vectors in Carte-sian coordinates. Using the Hamiltonian (1), the equa-tions of motion are ˙ x ia = ∂H∂p ia , (5) − ˙ p ia = ∂H∂x ia , (6)˙ s ia = ∂H∂s ja s ka (cid:15) ijk . (7)The first term in (2) is the Hamiltonian for n -particlesinteracting under Newtonian gravity H (0) N = 12 n (cid:88) a (cid:126)p a m a − n (cid:88) a,b (cid:54) = a m a m b r ab . (8)The explicit form of the PN terms in (2) and the equationof 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 tothe equations of motion derived from H (0)N , H (0)N + H (1)PN , H (0)N + H (1)PN + H (2)PN and H (0)N + H (1)PN + H (2)PN + H (2 . , respec-tively. We denoted as radiative Newtonian and radiative1 PN the equations of motion derived from H (0)N + H (2 . and H (0)N + H (1)PN + H (2 . , respectively. B. Chaos indicators
There are a number of methods for diagnosing chaosin a dynamical system: the Poincar´e surface (for systemwith less than 3 degrees of freedom), the Kolmogorov-Sinai entropy, the characteristic Lyapunov exponent andthe fractal basin boundary method, among others (seee.g. [54–56]). For our analysis, we employed the frac-tal basin boundary method and the Lyapunov exponent.In this section, we describe the implementation of bothmethods.The basin boundary method and the Lyapunov expo-nents are able to characterize different types of chaos.The basin boundary characterizes the sensitive depen-dence on initial conditions for nearby orbits with differ-ent attractors on the phase-space . However, it does notprovide information about the orbit itself. Alternatively,the Lyapunov indicator characterizes the regularity orchaos of an orbit in a way that is independent from theattractor.
1. Fractal basin boundary
A dynamical system usually exhibits a transient be-havior followed by an asymptotic regime. The time-asymptotic behavior defines attracting sets in the phase-space which are called attractors. The attractor can beperiodic, quasi-periodic or chaotic. The closure of theset of initial conditions, which has a particular attrac-tor, is called basin of attraction. Basin boundaries for atypical dynamical system can be either smooth or fractal[56, 57]. A system with fractal basin boundaries is sensi-tive to initial uncertainty. Given (cid:15) >
0, we call the initialcondition point P ∗ unsafe if there is another initial con-dition point P inside a neighborhood of size (cid:15) centered at P ∗ which converges to a different attractor. The fraction f of unsafe points is related to (cid:15) by f ( (cid:15) ) ∼ (cid:15) D − d f , (9)where D is the dimension of the phase-space and d f isthe dimension of the basin boundary (see [56, 58] for adetailed discussion). It is useful to define an uncertaintyexponent α = D − d f , which satisfies the relation 0 <α ≤
1. The value of α quantifies the presence of chaosin a basin boundary, α ≈ α ≈ box-dimension whichis defined as d box := lim δ → ln( N δ ( A ))ln(1 /δ ) , (10)where A is a nonempty bounded set in the d E -dimensional Euclidean space and N δ ( A ) is the smallestnumber of sets in a δ -cover of A . A δ -cover of A isa collection of sets of diameter δ whose union contains A . For a physical or numerical experiment, the fractalbasin boundary has finite resolution. A graphical rep-resentation of the basin boundary is given by an imageformed by pixels of size l res . In order to estimate the di-mension of the basin boundary, we compute the quotient An attractor is a set towards which a dynamical system asymp-totically 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 typicalbehavior of the resulting data is well approximated by astraight line. A linear regression method is used to fita linear function to the resulting data (where δ n = nl res and n ∈ { , , . . . , } ). Therefore, the value of the fittedfunction for δ = 0 gives an estimate value of the dimen-sion. We tested the method by computing the dimensionof non-fractal one-dimensional sets (i.e., a regular basinboundary) and Sierpinski’s gasket and carpet.
2. Lyapunov exponent
Lyapunov exponents are another way for characteriz-ing chaos in a dynamical system (see e.g. [54]). For anautonomous n -dimensional system˙ (cid:126)X = (cid:126)F ( (cid:126)X ) , (11)a given reference solution (cid:126)X ∗ and a nearby solution (cid:126)X = (cid:126)X ∗ + δ (cid:126)X . While the evolution of the difference δ (cid:126)X isgiven by the linearized equation δ ˙ (cid:126)X = J ( (cid:126)X ∗ ) · δ (cid:126)X, (12)where J ( (cid:126)X ∗ ) i,j := ∂ j F i ( (cid:126)X ∗ ) is the Jacobian matrix of (cid:126)F evaluated at the reference solution (cid:126)X ∗ . The solution of(12) is given by [60] δ (cid:126)X ( t ) = e J t · δ (cid:126)X (0) . (13)The Lyapunov exponents are defined by the eigenvaluesΛ i ( t ) of the distortion matrix Λ := e ( J + J T ) t λ i := lim t →∞ t ln Λ i ( t ) . (14)Sensitive dependence on initial conditions exists when λ i >
0. Therefore, in order to distinguish between reg-ular and chaotic orbits, it is sufficient to compute theprincipal Lyapunov exponent λ p := max( λ i ). Calculat-ing λ p directly is computationally expensive. However,by using the evolution of the difference, δ (cid:126)X , it is possibleto obtain an estimation of λ p λ ∗ p := lim t →∞ lim δX (0) → λ ∗ p ( t, δX (0)) , (15)where, λ ∗ p ( t, δX (0)) := 1 t ln (cid:18) δX ( t ) δX (0) (cid:19) , (16)and δX ( t ) = | δ (cid:126)X | . We will refer as a Lyapunov in-dicator to λ ∗ p in order to distinguish it from the prin-cipal Lyapunov exponent and the Lyapunov function, λ ∗ p ( t, δX (0)). The norm, | δ (cid:126)X | , plays an important rolein the computation of the Lyapunov indicator. Partic-ularly, in the case of general relativity, additional gaugeeffects can affect the estimation of the Lyapunov indica-tor [23, 55, 61]. Results from literature were applied togeneral relativistic dynamic of test particles and to bi-nary systems in the center of mass frame. Additionally,the estimation of λ p is often computed using the two-particle method instead of the evolution of the difference(see e.g. [62]). However, we did not noticed significantdifferences in the computation of λ ∗ p given the gauge ef-fects. We present the results in terms of the Cartesianlength of δ (cid:126)X .In practice, the Lyapunov indicator is computed ap-proximating the limits, using a sufficient small δ (cid:126)X (0)and sufficient large integration time. In our simulations,we used | δ (cid:126)X (0) | = 10 − . The integration time dependson the system (see Sec. III). The evolution is performedby solving numerically the system of equations (5)-(7),therefore (cid:126)X = ( (cid:126)x a , (cid:126)p a , (cid:126)s a ) T is a vector with 27 compo-nents. The right hand side of (11) is the vector (cid:126)F = (cid:18) ∂H∂p ia , − ∂H∂x ia , ∂H∂s ja s ka (cid:15) ijk (cid:19) , (17)and the equation of motion for the difference δ (cid:126)X =( δ(cid:126)x a , δ(cid:126)p a , δ(cid:126)s a ) is δ ˙ x ia = δ (cid:126)X · ∇ ∂H∂p ia , (18) − δ ˙ p ia = δ (cid:126)X · ∇ ∂H∂x ia , (19) δ ˙ s ia = δ (cid:126)X · ∇ ∂H∂s ja s ka (cid:15) ijk , (20)where we define the auxiliary operator δ (cid:126)X · ∇ := (cid:88) b =1 3 (cid:88) j =1 (cid:32) δx jb ∂∂x jb + δp jb ∂∂p jb + δs jb ∂∂s jb (cid:33) . (21) III. SIMULATIONS AND RESULTSA. Numerical methods
We solved the equations of motion numerically usingthe GNU Scientific Library (GSL) [63] and the
Ollip-tic code infrastructure [42]. We generated the righthand side (RHS) of the equations of motion (5)-(7) with
Mathematica (cid:126)x a , (cid:126)p a , (cid:126)s a ) was set to 10 − andto 10 − for the variation variables ( δ(cid:126)x a , δ(cid:126)p a , δ(cid:126)s a ). Weused as an indicator of accuracy the estimate of the localerror provided by the GSL routines and the conserva-tion of the Hamiltonian for the simulations which do notincluded the radiation term. For our stability analysis, we solved the system of equa-tions (5)-(7) and (17)-(19) simultaneously. In order toperform the simulations in an efficient way, we adaptedthe number of variables depending on the problem. Theparameters to consider are the number of bodies nb , themagnitude of the spin (cid:126)s a and whether we are solving aplanar or non-planar problem. As we mentioned before,the RHS of Eqs. (5)-(7) were generated by a Mathe-matica script. For the Newtonian case and 1 PN cases,it is convenient to compute the analytical expression forthe RHS of Eqs. (18)-(20). The analytical expression forthe 2 PN and the spinning cases produces large compu-tational source files. The source file for the 2 PN is of25 megabytes while the inclusion of the spin generatesfiles over a hundred megabytes. The compilation andoptimization of large files is not practical, since compil-ers run out of memory even in the 16 GB of RAM serverthat we tried. Therefore, for the 2 PN order and the spin-ning case we computed the RHS of Eqs. (18)-(20) numer-ically. In order to estimate the errors in the numericalcomputation of the Jacobian, we compared the resultsobtained for the Newtonian and 1 PN cases employingboth methods (see Sec. III A 1). The system (18)-(20) isthe Jacobian matrix of (17). The problem reduces to thecomputation of J i,j = ∂F i ( (cid:126)X ∗ ) ∂X j . (22)A simple 2nd, 4th and even 6th order finite differencemethod with a fix step size h does not produce accuratesolutions. The main issue is the nature of the compo-nents, for example, the position variables require a dif-ferent optimal step-size than the momentum or the spinvariables. We implemented an adaptive method basedon the GSL differentiation routine. The method pro-duces a solution with an estimated error below 10 − (seeSec. III A 1).An important issue in the numerical integration of athree-body system arises when two of the bodies are veryclose to each other. In the case of adaptive step sizemethods, it is necessary to reduce the step size in or-der to resolve properly the orbits in the close interactionphase. A usual approach to deal with this problem is toperform a regularization of the equations of motion, seee.g., [65–69] and references therein. However, in our sim-ulations, we included a different criteria. Regardless ofthe equation of motion employed, we monitored the ab-solute value of each conservative part of the Hamiltonian(1) relative to the sum of the absolute values H % i := 100 (cid:32) | H i || H (0)N | + | H (1)PN | + | H (2)PN | + L ( H Spin ) (cid:33) . (23)The simulations stop when the contribution of the firstpost-Newtonian correction is larger than 10%. We willrefer to H %1 >
10% as strong interaction. Empirical re-sults showed that, for a similar configuration, the resultedwaveform exhibits a growth characteristic of the mergerphase [39]. The dynamics after a strong interaction maylead to the merger of two of the bodies or to the escapeof the lighter body. In order to determine safely whetherthere is a merger or not, it is necessary to employ fullnumerical relativistic simulation.Additionally, we considered the Newtonian escape en-ergy of each body in respect of the other two components E a esc := 12 µ a (cid:18) (cid:126)p a m a − (cid:126)p ∗ a m ∗ a (cid:19) − µ a | (cid:126)x a − (cid:126)x ∗ a | , (24)where µ a := m a m ∗ a , m ∗ a := (cid:80) b (cid:54) = a m b , (cid:126)x ∗ a :=( m ∗ a ) − (cid:80) b (cid:54) = a m b (cid:126)x b and (cid:126)p ∗ a := (cid:80) b (cid:54) = a (cid:126)p b . The simulationstops when E aesc > S := (cid:80) a | (cid:126)x a | is 100 times the initial size. Dur-ing a close encounter of two of the bodies it is common tohave 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 er-rors. Here we summarize the results of these tests.
1. Numerical tests
Our first test was done with Lagrange’s triangle solu-tion using Newtonian equations of motion. In Lagrange’ssolution each body is sitting in one corner of an equilat-eral triangle (see e.g. [70]). We set the sides of such trian-gle to r = r = r = 100, the mass ratio to 1:10:100,and the eccentricity to zero. Then each body follows acircular orbit (with different radii) around the center ofmass. This configuration is unstable and after 14 periodsone of the bodies escapes from the system. The unsta-ble property is reflected in the Lyapunov indicator λ ∗ p .Fig. 1 (a) shows log ( λ ∗ p ) computed using two methodsfor the prescription of the RHS of Eqs. (17)-(19). Thefirst method is the analytic expression (labeled J a ) andthe second, the numerical solution (labeled J n ). The twomethods produce solutions that cannot be distinguish-able within the plot. The maximum relative difference is4 × − . The Lyapunov indicator exhibits the charac-teristic behavior of a chaotic system. After t = 10 , λ ∗ p ( t )oscillates around a constant value. Fig. 1 (b) shows therelative variation of the Hamiltonian∆ H := H (0) − H ( t ) H (0) , (25)for the same solution. The conservation of the Hamilto-nian has a similar behavior either using J a or J n . TheHamiltonian exhibits a variation of around 10 − beforethe evolution stops.Fig. 1 (c) illustrates the L norm of the estimatederrors for Lagrange’s triangle solution. By looking atthe errors, the difference between the solution producedby J n and J a becomes clear. The numerical solutiongenerates errors which are two orders of magnitude largerthan the analytical Jacobian. However, even when using J a , the accuracy is good, with errors below 1 . × − . -4-3.8-3.6-3.4-3.2-3 2 2.5 3 3.5 4 4.5 5 l og ( λ ∗ p ) (a)-1.2-0.9-0.6-0.30 2 2.5 3 3.5 4 4.5 5 ∆ H (b) [ × − ] L ( E i ) (c) log ( t ) [ × − ] J a J n FIG. 1. Newtonian Lagrange’s triangle solution. The up-per panel (a) shows the evolution of the Lyapunov indicator λ ∗ p ( t, δX (0)), the middle panel (b) shows the conservation ofthe Hamiltonian ∆ H and the lower panel (c) shows the L norm of the estimated errors. In every panel, the solid linedenotes the solution obtained using the analytical Jacobianand the dashed line, the result provided by the numericalcomputation of the Jacobian (see Sec. III A). A second test was done for a stable system. Using the1 PN equations of motion, we computed the Lyapunovindicator λ ∗ p , the relative variation of the Hamiltonian∆ H and the L norm of the estimated errors for theH´enon’s criss-cross solution [32, 71, 72]. We used initialparameters given by (cid:126)x (0) = 1 . λ ˆ x, (cid:126)p (0) = 3 − / · . λ − ˆ y,(cid:126)x (0) = − . λ ˆ x, (cid:126)p (0) = − − / · . λ − ˆ y,(cid:126)x (0) = − . λ ˆ x, (cid:126)p (0) = 3 − / · . λ − ˆ y, where ˆ x , ˆ y and ˆ z are the unitary basis vectors in Carte-sian coordinates, and λ is a scaling factor (for our sim-ulation λ = 10). Notice that for this test we used theparameters given in [73] with the scaling factor λ , anddoing a change of variables from initial velocity to ini-tial momentum. Therefore, we are not including post-Newtonian corrections to the initial parameters. As inthe previous case, λ ∗ p does not show significant differenceswhen computed using J a or J n (see Fig. 2). The valueof λ ∗ 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 is8 × − . The Hamiltonian is conserved with a maximumvariation of 5 × − . The L norm of the estimated er-rors shows noisy errors using J n while J a exhibits a flattrend. -6-5.5-5-4.5-4-3.5-3 2 3 4 5 6 7 l og ( λ ∗ p ) (a)-4.8-3.6-2.4-1.201.2 2 3 4 5 6 7 ∆ H (b) [ × − ] L ( E i ) (c) log ( t ) [ × − ] J a J n FIG. 2. 1 PN H´enon Criss-cross solution. The upper panel (a) shows the evolution of the Lyapunov indicator λ ∗ p ( t, δX (0)),the middle panel (b) shows the conservation of the Hamilto-nian ∆ H and the lower panel (c) shows the L norm of theestimated errors. In every panel, the solid line denotes the so-lution obtained using the analytical Jacobian and the dashedline, the result provided by the numerical computation of theJacobian. A third test was done with a known chaotic binarysystem. We used a chaotic eccentric configuration withmass ratio 3:2 and initial parameters: (cid:126)r (0) = 50ˆ x(cid:126)p (0) = − (cid:126)p (0) = 0 . y + 0 . z,(cid:126)s (0) = − . x − . y + 0 . z,(cid:126)s (0) = − . x + 0 . y − . z. (26)A similar system, with a different unit convention, wasstudied in [74] (Sec. IV-C1). Following [74], we solvedthe equation of motion using up to 2 PN correction forthe orbital part and leading order in the spin. As wasemphasized in [74], in this configuration, the two com-pact bodies are rather close to each other. In order toobtain an accurate physical description, it is necessaryto include, for the orbital part, corrections higher than2 PN and, for the spinning part, higher than the leadingorder. The 1 PN contribution to the Hamiltonian (22)during a close encounter of the bodies can reach 50%. Atthis point it is important to recall one of the conclusionsof [74]. When using up to 2 PN Hamiltonian with leadingorder spin contributions, there is no evidence of chaoticsystems for configurations which are physically valid tothe PN order considered. Therefore, for comparison, thistest was performed without using the stopping criteria H %1 > -6-5-4-3 2 3 4 5 6 7 l og ( λ ∗ p ) (a)-0.6-0.300.30.60.9 2 3 4 5 6 7 ∆ H (b) [ × − ] L ( E i ) (c) log ( t ) [ × − ] | ~ s | , | ~ s | = FIG. 3. Spinning binaries with an eccentric orbit. Lyapunovfunction λ ∗ p ( t, δX (0)) (a) , conservation of the Hamiltonian∆ H (b) and the L norm of the estimated errors, for thesystem with initial parameters given by (26) (solid line) andthe same configuration except for (cid:126)s (0) = (cid:126) system was evolved using the parameters given in (26),and for reference, a non-chaotic orbit with the same ini-tial condition but with (cid:126)s (0) = (cid:126) (a) shows the Lyapunov indicator for the chaotic configura-tion (solid line) and the arithmetic mean computed for t > . The estimated value of the Lyapunov indicatoris (cid:104) λ ∗ p (cid:105) = 10 − (dotted line) with a standard deviation of1 . × − , corresponding to a relative variation of 13%.In contrast, the regular solution shows that λ ∗ p decreasesalmost monotonically after t = 10 . In both simulations,the Hamiltonian is numerically conserved below 10 − (seeFig. 3- (b) ), and the estimated errors are of order 10 − for the chaotic solution and 10 − for the regular solu-tion. The main source for the errors, in the case of thechaotic solution, is the numerical calculation of the Ja-cobian matrix.We performed, additional tests already presented inSec. 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 binarycompact object system due to a third smaller compactobject. As a basic configuration, we studied a Jacobiansystem with mass ratio 10:20:1. The inner binary sys-tem had initial apo-apsis ( i.e. maximum separation of a m m m i Π ext Π in FIG. 4. Hierarchical system (see [39]). Initial configurationof the inner and external binaries. The initial momentum ofthe third body is given by considering the external binary asa Newtonian binary. Shown are the osculating orbital planesΠ in and Π ext for inner and external binary orbits. The twoplanes are inclined by an angle ι . Keplerian orbit) r b (0) = 130 and eccentricity e b (0) = 0.Although there are methods in post-Newtonian dynam-ics to specify the initial parameters of a binary systemwith a given eccentricity (see e.g. [74–79]), in this studythe 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 consideringonly the Newtonian dynamics of a non-perturbed binary,where the eccentricity refers to the Newtonian case. Inthis approach, we view the third compact body and thecenter of mass of the inner binary as a new binary (wewill refer to it as the external binary). The bodies startfrom a configuration where the initial radial vector (cid:126)r isperpendicular to initial vector position (cid:126)x of the exter-nal body (see Fig. 4). We denote the inclination anglebetween the osculating orbital planes Π in and Π ext by ι (see Fig. 4). To be more specific, the initial position andmomentum of each body is given by (cid:126)x = m m b r b ˆ x − m r (ˆ y cos i + ˆ z sin i ) ,(cid:126)x = m m b r b ˆ x − m r (ˆ y cos i + ˆ z sin i ) ,(cid:126)x = m b r (ˆ y cos i + ˆ z sin i ) ,(cid:126)p = − p b ˆ y − m m b p ˆ x,(cid:126)p = p b ˆ y − m m b p ˆ x,(cid:126)p = p ˆ x, (27)where m b := m + m , p b = m m (cid:112) (1 − e b ) / ( m b r b ) and p = m b m (cid:112) (1 − e ) /r . Notice that we set the mo-mentum in terms of the eccentricity and the apo-apsis.Therefore, a higher eccentricity given to the external bi-nary implies a stronger perturbation of the inner binary.The main goal of this study was to characterize thestability and chaos of a hierarchical system as function ofthe initial apo-apsis r , the eccentricity e for the external binary and the inclination angle ι between the osculatingorbital planes. We explored the influence of the post-Newtonian corrections.
1. Final state survey
An analysis of the asymptotic behavior of the systemas function of three parameters r , e and ι was per-formed. For a given osculating angle ι , we produceda map which is a subset of the space of configurations( r , e ). We define three possible outcomes which char-acterize the asymptotic behavior of the system: the es-cape of one of the bodies, a strong interaction H %1 > t f = 2 × .We assign a color to each outcome in the map of ini-tial configurations ( r , e ): dark gray (color online) foran escape, white for a strong interaction and black for astable configuration. Strictly speaking, the result is nota basin boundary map since the outcomes are not attrac-tors in the phase space and the parameters ( r , e ) arenot coordinates of the phase space. However, Eqs. (27)gives the link between ( r , e ) and the coordinates of thephase space. On the other hand, a system where one ofthe bodies escapes will be attracted to a hyper-plane inthe phase space with one of the spatial coordinates at in-finity; a strong interaction may become an escape or themerger of two of the bodies. A merger can be representedby a hyper-plane where two of the bodies share the samevalues of position and momentum. For a stable system, itis hard to define a single attractor because it can exhibita different asymptotic behavior for t > t f . For simplicity,we refer to the boundary between the various color in themap as basin boundary since it can capture some of thefeatures of a basin boundary.The choice of t f = 2 × is motivated by the dy-namics of a similar Jacobian system presented in [39].Including up to 2.5 PN terms, the inner binary arrivesto the merger phase after t mgr ≈ . × . The assump-tion is that if a conservative system remains stable untilthat time, then the corresponding radiative system willremain stable until the merger of the inner binary. Weselected a shorter time (i.e., t f < t mgr ) due to compu-tational costs. However, a set of numerical experimentsfor different values of t f suggest that a non-stable orbitis manifest before t f .We explored the parameter space r , e using 4 regionsdefined by R : [200 , × [0 , , ∆ r = 1 , ∆ e = 0 . R : [200 , × [0 , . , ∆ r = 0 . , ∆ e = 0 . R : [275 , × [0 . , . , ∆ r = 0 . , ∆ e = 0 . R : [295 , × [0 . , . , ∆ r = 0 . , ∆ e = 0 . . The domain of each region is the Cartesian product andits resolution is denoted by ∆ r , ∆ e . Each map wasdone computing 300 ×
400 orbits except for region R that contains 200 ×
400 orbits. Fig. 5 shows the resultfor a Newtonian simulation and osculating angle i = 0. FIG. 5. Fractal basin boundary (Newtonian). Result of theasymptotic behavior of initial conditions for the parameters r , e and i = 0 for region R (a) , R (b) , R (c) and R (d) . The black points denote stable orbits, dark gray points(blue color online) are escapes of the lighter body and whitepoints denote a strong interaction of two of the bodies.TABLE I. Newtonian orbits. Here d box is the box dimen-sion of the basin boundaries presented in Fig. 5, the columns‘Regular’, ‘Escape’ and ‘ H > box Regular Escape H > R . ± .
015 8.2% 32.4% 59.4% R . ± .
013 1.0% 41.3% 57.7% R . ± .
015 0.0% 24.7% 75.3% R . ± .
007 0.0% 27.2% 72.8%
Notice that R ⊂ R ⊂ R ⊂ R , i.e. each R i is a magni-fication of the region R i − (in the figure, the sub-regionsare indicated by rectangles). Fig. 5 shows the character-istic behavior of a fractal set; every magnification revealsa more complex structure from which it is possible todistinguish some degree of self-similarity. Table I sum-marizes the quantitative results for the basin boundary.Notice that the box-dimension reflects the fact that thesets are not perfectly self-similar. We choose R and R to perform additional simulations.A comparison using Newtonian and 1 PN orbits for re-gions R and R is presented in Fig. 6. There is a cleardifference in the set of stable orbits. In the Newtoniancase, there are two disjoint sets where the 1 PN systemexhibits a single set (compare the black point in panels (a) and (b) of Fig. 6). The magnification R shows dif-ferences in the escape and strong interaction zones. The FIG. 6. Fractal basin boundary (Newtonian and 1 PN). Com-parison of regions R and R for Newtonian and 1 PN orbits.Panels (a) and (b) show R for Newtonian and 1 PN respec-tively. Similarly, panels (c) and (d) show R for Newtonianand 1 PN, respectively. (c) and (d) of Fig. 6). Our comparison includes aset of simulations for leading order spinning 1 PN parti-cles. We employed maximally spinning particles χ a = 1and θ a = 0 , φ a = 0 except by θ = π/ , θ = 2 π/ configuration-a ) is given by (cid:126)s = m ˆ x,(cid:126)s = − m ˆ x,(cid:126)s = m ˆ z, (28)Opposite to the non-spinning case, the orbits for our con-figuration of spinning particles are not coplanar. How-ever, the change in the final states is negligible, the re-sulting map is almost indistinguishable to the 1 PN case(for brevity, we do not display the basin boundaries whichare similar to panels (b) and (d) of Fig. 6). The quan-titative analysis of the six basin boundaries under com-parison is displayed in Table II. The difference on thesubstructures between the Newtonian case and the 1 PNcase is reflected by the box-dimension (smaller value forthe 1 PN case). In the region R , there is a small changein the distribution of ‘Regular’, ‘Escape’ and strong in-teractions between Newtonian and 1 PN. Nevertheless, inregion R there are noticeable differences in the percent-ages of ‘Escape’ and strong interactions. The similaritybetween the spinning and the non-spinning 1 PN case ap-pears for the dimension and the distribution of the threeoutcomes (with minor differences in the percentages). TABLE II. Comparison of Newtonian, 1 PN, radiative 1 PNand leading order spinning 1 PN orbits. Here column PNtakes value 0 for Newtonian, 0+1 for 1 PN and 0+1+2.5 forradiative 1 PN. The column Sp takes value 0 for non-spinningparticles and 1 for the spinning ones. The meaning of theremaining columns is similar to Table I. See the text for detailsabout the parameters of the spinning case (Fig. 6 and 7 showsthe non-spinning cases).PN Sp d box
Regular Escape H > R . ± .
013 1.0% 41.3% 57.7%0+1 0 1 . ± .
014 1.4% 40.3% 58.3%0+1 1 1 . ± .
014 1.2% 40.7% 58.1%0+1+2.5 0 1 . ± .
014 1.7% 40.0% 58.3% R . ± .
015 0.0% 24.7% 75.3%0+1 0 1 . ± .
015 0.0% 18.0% 82.0%0+1 1 1 . ± .
015 0.0% 18.1% 81.9%FIG. 7. Basin boundary of region R . The dynamic includes(0+1+2.5) PN terms. In order to investigate the influence of the gravitationalradiation, we performed an evolution which contains the2.5 PN terms. The inclusion of PN terms is computa-tionally expensive. Therefore, we used a radiative 1 PNHamiltonian (i.e., we include the Newtonian, 1 PN and2.5 PN terms in (2)). For similar configurations, the dif-ference in the dynamics between the full 2.5 PN systemand the radiative 1 PN is relatively small [39]. Fig. 7shows the fractal basin boundary and Table II the corre-sponding quantitative measures. The radiative term doesnot seem to produce a significant change in respect of theinclusion of 1 PN terms (compare Fig. 7 with Fig. 6- (b) ).Particularly, the box-dimension of the boundary is thesame, and there is only a small difference in the distribu-tion of the percentages of escape and strong interactions.The osculating angle ι has a strong influence in the re-gion of the initial parameters R . For the Newtoniancase, we produced a set of eight additional maps forregion R where the osculating angle takes the values FIG. 8. Fractal basin boundary of region R as function ofthe osculating angle ι (Newtonian). The meaning of the coloris the same as in previous plots (see e.g. Fig. 6). From thetop to the bottom and from the left to the right ι takes thevalues 0 , π/ , π/ , π/ , π/ , π/ , π/ , π/ , π/ R as function ofthe osculating angle ι (1 PN). The meaning of the color is thesame as in previous plots (see e.g. Fig. 6). From the top tothe bottom and from the left to the right ι takes the values0 , π/ , π/ , π/ , π/ , π/ , π/ , π/ , π/ α π/ π/ π/ π/ π/
16 3 π/ π/ π/ (b) ι α Newtonian α Nfit α PNfit
FIG. 10. Uncertainty exponent α as function of the osculat-ing angle ι for the region R . The upper panel (a) shows theresult for the Newtonian case, the solid line is the fitting func-tion α Nfit ( ι ) = a N ι + b N + c N sin ϕ N ι . The lower panel (b) showsthe corresponding result for the 1 PN simulations, where thesolid line is α PNfit ( ι ) = a PN ι + b PN + c PN sin ϕ PN ι . The fittingparameters are given in Eq. (29). ι = nπ/ n ∈ { , , . . . , } . The results, including thecase ι = 0, are presented in Fig. 8. An analogous set ofsimulations was produced using the 1 PN equation of mo-tion. Fig. 9 shows the corresponding basin boundaries.Table III shows the quantitative results for both sets.Notice that for i > π/ α . A relatively simple func-tion α fit ( ι ) = a i ι + b i + c i sin ϕ i ι fits the data. The fitparameters are a N = − . ± . , a = − . ± . ,b N = 0 . ± . , b = 0 . ± . ,c N = 0 . ± . , c = 0 . ± . ,ϕ N = 7 . ± . , ϕ = 6 . ± . . (29)The slope of the linear part of the 1 PN is 1.6 timeslarger than the Newtonian one. However, the oscillatorypart is nearly 0.7 times the 1 PN value for the Newtoniancase. This result shows that the chaotic properties of thishierarchical configuration increases in both cases, reach-ing the maximum at ι = π/
2. From the basin bound-ary figures, we can notice an increasing number of unsafepoints. There is an evident difference for ι = π/ l og ( a i n , e x t ) (a)Newtonian1.822.22.42.6 l og ( a i n , e x t ) (b)1 PN-6-5-4-3 2 3 4 5 6 7 l og ( λ ∗ p ) (c) log ( t ) inner externalNewtonian 1 PN FIG. 11. Semi-major axis and Lyapunov function for initialparameters ι = π/ , r = 340 , e = 0. The upper panel (a) shows the semi-major axis a of the inner and external binariesas function of time for the Newtonian case. The middle panel (b) is similar to (a) but for the 1 PN case. The lower panel (c) shows the evolution of the Lyapunov indicator for bothcases. Notice that the quantities are given in log scales.
2. Analysis of orbits
We used the basin boundaries as a guide for a detailedstudy of specific orbits. Given the region R , an osculat-ing angle ι , two different basin boundaries Ω( r , e ) andΩ ∗ ( r , e ), it is possible to produce a new map D thatshows the differences between Ω( r , e ) and Ω ∗ ( r , e ).We constructed the map as follows; D ( r , e ) takes value1 if Ω( r , e ) = Ω ∗ ( r ∗ , e ∗ ) in a neighborhood of ( r , e )of size at least 6d r × e , i.e., for r ∗ = r + i d r , e ∗ = e + i d e and i ∈ {− , − , . . . , } . If the previ-ous condition is not satisfied, then the D ( r , e ) is set to0. D ( r , e ) provides a way for identifying zones wherethe two basin boundaries differed in a relatively largeregion. From the analysis of this map and the basinboundaries, we explored several initial parameters wherewe compared the orbits for different configurations of thePN equations of motion. In the following, we describe 5representative comparisons.For the region R and ι = π/ e < ( r − /
320 ends when H > TABLE III. Measures of the basin boundaries of Figs. 8 and 9. The osculating angle is denoted by ι . The meaning of theremaining columns is similar to Table II.Newtonian 1 PN ι d box Regular Escape H >
10% d box
Regular Escape H > . ± .
013 1.0% 41.3% 57.7% 1 . ± .
014 1.4% 40.3% 58.3% π/
16 1 . ± .
013 1.0% 43.0% 56.0% 1 . ± .
014 0.7% 41.6% 57.7% π/ . ± .
013 1.1% 49.2% 49.7% 1 . ± .
013 0.9% 48.2% 50.9%3 π/
16 1 . ± .
012 1.9% 62.3% 35.8% 1 . ± .
012 1.5% 58.0% 40.5% π/ . ± .
012 2.7% 67.9% 29.4% 1 . ± .
012 2.2% 64.2% 33.6%5 π/
16 1 . ± .
012 2.7% 74.3% 23.0% 1 . ± .
012 3.2% 67.9% 28.9%3 π/ . ± .
013 2.2% 77.5% 20.3% 1 . ± .
012 4.3% 70.8% 24.9%7 π/
16 1 . ± .
012 1.5% 73.9% 24.6% 1 . ± .
012 6.1% 69.3% 24.6% π/ . ± .
011 0.1% 57.1% 42.8% 1 . ± .
011 9.3% 67.4% 23.3% growth on the eccentricity of the inner binary induced bythe external body (see Fig. 11- (a) ). This effect is pro-duced by a resonance between the orbital dynamic of theinner and external binary (see [80] for a detailed descrip-tion). For the 1 PN case the orbits remain stable duringthe simulation (see Fig. 11- (b) ). In both cases, the evolu-tion of the Lyapunov indicator does not showed evidenceof exponential growth behavior (see Fig. 11- (c) ). Noticethat in the Newtonian case, the strong interaction is be-tween the components of the inner binary. The behaviorof the Lyapunov function suggests that the strong inter-action between the heavy components of the system doesnot follow chaotic trajectories.For the region R and ι = 0, there is a significantdifference between the Newtonian and the 1 PN basinboundaries for stable orbits (i.e., the black dots in panels (a) and (b) of Fig. 6). The Newtonian configuration re-mains stable without noticeable instability. On the otherhand, the coupling between the orbits in the 1 PN caseresults in the escape of the third body (see panels (a) and (b) of Fig. 12). The Lyapunov function shows for1 PN dynamics the characteristic behavior of escape or-bits. During a close encounter between the lighter bodyand one of the heavy components of the inner binary,the difference vector has an exponential growth whichis followed by a regular behavior. The regular behavioris 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 anindication of sensitivity to initial conditions. In general,the resulting escape orbit is chaotic; a small change inthe initial condition produces a significant change in theescape direction.A comparison between 1 PN and 2 PN dynamicsshowed that almost identical orbits are produced by anevolution where the lighter body has a quick encounterwith the inner binary which is followed by an escape ora strong interaction. Fig. 13 displays a comparison be-tween the Newtonian, 1 PN and 2 PN evolutions for 11nearby configurations. The initial parameters are ι = 0 l og ( r ab ) (a)Newtonian1.82.433.64.24.8 l og ( r ab ) (b)1 PN-6.3-5.6-4.9-4.2-3.5-2.8 2 3 4 5 6 7 l og ( λ ∗ p ) (c) log ( t ) a=1,b=2 a=1,b=3 a=2,b=3Newtonian 1 PN FIG. 12. Relative separation and Lyapunov function for ini-tial parameters ι = 0 , r = 320 , e = 0. The upper panel (a) shows the relative separation r ab between particles as func-tion of time for the Newtonian case. The middle panel (b) is similar to (a) but for the 1 PN case. The lower panel (c) shows the evolution of the Lyapunov indicator for both cases. and e = 0 . r ∈ { . , . , . . . , . } . Therelative change in the initial separation r between theconfigurations is ∼ . (a) ), exhibits a chaotic behavior. Every simula-tion ends with the escape of the lighter body. There aresignificant differences in the trajectories. On the otherhand, the 1 PN and 2 PN evolutions (panels (b) and (c) of Fig. 13) produce a quick escape which are similarfor each initial condition. The maximum final differencebetween the orbits is 1 .
3% for the 1 PN case and 1 . l og ( | ~ x | ) (a)Newtonian r (0) = + . · i l og ( | ~ x | ) (b)1 PN1.82.433.64.24.8 2 3 4 5 6 7 l og ( | ~ x | ) (c) log ( t ) i = − i = − i = − i = − i = − i = i = i = i = i = i = FIG. 13. Evolution of the coordinate position | (cid:126)x | for 11nearby configurations. Initial parameters ι = 0 , e = 0 . r ∈ { . , . , . . . , . } . The upper panel (a) showsNewtonian case, the middle panel (b) shows the 1 PN evolu-tion and the lower panel (c) the 2 PN case. between the 1 PN and 2 PN reference evolution r = 270is 6%.However, for successive encounters the evolution is sig-nificantly different. Fig. 14 shows the resulted orbits forinitial parameters ι = 0 , r = 240 , e = 0. For theseinitial parameters, the resulting outcome is the escapeof the lighter body. Nevertheless, in the 1 PN evolu-tion there is an ejection of the lighter body before thefinal escape. In contrast, the ejection is not present inthe 2 PN orbit (compare panels (a) and (b) of Fig. 14).There is a significant difference in the evolution of theLyapunov function. For the 1 PN case, there is a dou-ble kink which corresponds to the close encounters beforeand after the ejection. Nevertheless, the 2 PN evolutionpresents 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 = 0 with r ∈ { . , . , . . . , . } . Similarly to the casepresented in Fig. 13, the relative change in the initialseparation r between the configurations is ∼ . r = 280. The evolutions exhibit sensitive de-pendence on initial conditions. In the three cases, Newto-nian, 1 PN and 2 PN, relatively nearby initial parameters In the three-body problem literature, an ejection refers to a tem-poral large separation of one of the bodies which finally returns[21]. l og ( r ab ) (a)1 PN1.82.433.64.24.8 l og ( r ab ) (b)2 PN-5-4.5-4-3.5-3 2 3 4 5 6 7 l og ( λ ∗ p ) (c) log ( t ) a=1,b=2 a=1,b=3 a=2,b=31 PN 2 PN FIG. 14. Relative separation and Lyapunov function for ini-tial parameters ι = 0 , r = 280 , e = 0. The upper panel (a) shows the relative separation r ab between particles as functionof time for the 1 PN case. The middle panel (b) is similarto (a) but for the 2 PN case. The lower panel (c) shows theevolution of the Lyapunov indicator for both cases. produce a significant change in the orbits. Moreover, theinclusion of PN corrections produced a noticeable changein the trajectories.For the parameters that we have studied, the inclusionof spin correction has a small influence in the final out-come. However, like in the 2 PN case, for some of theinitial parameters, the evolution is completely different.For the initial parameters ι = 0 , r = 320 and e = 0 . (a) ). In contrast, the spin-ning 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 Lya-punov function which after t = 10 decreases linearly.The spinning case shows initially the same trend. How-ever, the strong interaction produces a short kink at theend of the evolution (see Fig. 16- (c) ).We explored a different set of spin parameters. Herewe present one of the results of configuration-b given by (cid:126)s = m ˆ x,(cid:126)s = m ˆ x,(cid:126)s = m ˆ y. (30)Fig. 17 shows the result for the initial parameters ι = π/ , r = 242 , e = 0 . l og ( | ~ x | ) (a)Newtonian r (0) = + . · i l og ( | ~ x | ) (b)1 PN1.82.433.64.24.8 2 3 4 5 6 7 l og ( | ~ x | ) (c) log ( t ) i = − i = − i = − i = − i = − i = i = i = i = i = i = FIG. 15. Evolution of the coordinate position | (cid:126)x | for 11nearby configurations. Initial parameters ι = 0 , e = 0 and r ∈ { . , . , . . . , . } . The upper panel (a) showsNewtonian case, the middle panel (b) shows the 1 PN evolu-tion and the lower panel (c) the 2 PN case. The bold-blackline is the reference solution presented in Fig. 14. l og ( r ab ) (a)1 PN1234 l og ( r ab ) (b)1 PN, Slo-a-5-4.5-4-3.5-3 2 3 4 5 6 7 l og ( λ ∗ p ) (c) log ( t ) a=1,b=2 a=1,b=3 a=2,b=31 PN 1 PN, Slo-a FIG. 16. Relative separation and Lyapunov function for ini-tial parameters ι = 0 , r = 320 , e = 0 .
39. The upper panel (a) shows the relative separation r ab between particles asfunction of time for the 1 PN case. The middle panel (b) 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 bothcases. l og ( r ab ) (a)1 PN11.522.53 l og ( r ab ) (b)1 PN, Slo-b-5-4.5-4-3.5-3 2 3 4 5 6 7 l og ( λ ∗ p ) (c) log ( t ) a=1,b=2 a=1,b=3 a=2,b=31 PN 1 PN, Slo-b FIG. 17. Relative separation and Lyapunov function for ini-tial parameters ι = π/ , r = 242 , e = 0 .
2. The upperpanel (a) shows the relative separation r ab between parti-cles as function of time for the 1 PN case. The middle panel (b) is similar to (a) but for the spinning leading order 1 PN configuration-b (initial spin given by (30)). The lower panel (c) shows the evolution of the Lyapunov indicator for bothcases.TABLE IV. Evolution of the eccentricity of the inner binary e b . Here we refer to the orbits in terms of the figure givenin the column ‘Fig.’, e b (0) is the initial eccentricity of theinner binary and e b ( t f ), the final eccentricity computed beforethe corresponding simulation finish. (cid:104) e b (cid:105) and σ ( e b ) are thearithmetic mean and standard deviation of e b respectively.The column ‘FS’ refers to the final state of the evolution. e b (0) e b ( t f ) (cid:104) e b (cid:105) σ ( e b ) FS1 11- (a) H > (b) (a) (b) (a) (b) (a) H > (b) (a) (b) H > e b . The eccentricity is computed by calculat-ing numerically the value of the apo-apsis and peri-apsis4of the inner binary i.e., computing a sequence of localmaxima and minima of r . The eccentricity is given by e b = r ap − r per r ap + r per . (31)Notice that the initial value of e b in every case is relativelysmall but not zero due to the perturbation of the thirdbody. However, the final value of the eccentricity dependson the final state of the evolution. For regular orbits,the change is relatively small (see the standard deviationof rows 2 and 3 in Table IV). On the other hand, thefinal eccentricity of the escaping orbits is relatively largeduring the whole evolution (compare the final eccentricityand the arithmetic mean of the corresponding rows). Theorbit with resonance exhibits the largest increment ineccentricity. The final value before the evolution stops isclose to 0.5. For PN evolutions, we did not find evolutionswith this kind of resonance. However, a different set ofparameters may lead to similar behavior. For a stronginteraction between two of the bodies, the resulting finaleccentricity is not meaningful since there are two possibleoutcomes in the full evolution scenario. One possibilityis an ejection or escape of the lighter body leading to aincrement in the eccentricity of the inner binary similar tothe cases presented. The other possibility is a successivemerger of the three bodies. IV. DISCUSSION
We performed a numerical study of stability and chaosof a hierarchical three compact object configuration. Theconfiguration is composed by a binary system which isperturbed by a third, lighter body positioned initiallyfar away from the binary. For a given inner binary, weexplored the influence of the third body taking as free pa-rameters the angle of the osculating orbital planes, theinitial apo-apsis and eccentricity of the third body. Usingthe basin boundary method, a total of 2,960,000 simula-tions were analyzed. From the total, 47.3% correspondedto Newtonian evolutions, 40.54% to (0+1) PN simula-tions, 4.05% to (0+1+2.5) PN evolutions and 8.11% aresimulations of spinning particles including leading orderin spin and 1 PN terms in the orbital part.A total of twenty-five basin boundary maps were pro-duced. The basin boundaries exhibit the characteristicsof fractal sets: some degree of self-similarity and an in-creasing complexity under magnification. Fractal basinboundaries are an indicator of chaos in dynamical sys-tems. By measuring the fractal dimension of the basinboundaries it is possible to determine the uncertaintyexponent α . The property of the exponent is that for α ≈ α = 1the system is regular. The values of α for the fractalbasin boundaries considered here are between 0.02 and0.26. The 1 PN set of data produces an exponent slightlylarger than the Newtonian case (the maximum differenceis 1.6%). On the other hand, the osculating angle ι has a strong influence in the uncertainty exponent. The ex-ponent decreases as a linear oscillatory function in therange [0 , π/ ι = 0) and the case where the third body starts from adirection perpendicular to the orbital plane of the innerbinary ( ι = π/
2) is 7.1% for the Newtonian simulationsand 8.1% for the 1 PN case.In addition to the uncertainty exponent, we quantifiedthe percentage of stable orbits, escapes and strong inter-actions. The configurations selected contain between 0%and 9.3% of stable orbits. The distribution of escapesare between 40.3% and 77.6%. The strong interactionsoscillate between 20.3% and 82%. A remarkable casewas found at ι = π/ ACKNOWLEDGMENTS
It is a pleasure to thank Todd Oliynyk, MarkFisher, Bernd Br¨ugmann, Cynnamon Dobbs and JenniferFern´andez for valuable discussions and comments on themanuscript. This work was supported in part by ARCgrant DP1094582, DFG grant SFB/Transregio 7 and byDLR grant LISA Germany. [1] K. Gultekin, M. C. Miller, and D. P. Hamilton, AIP Conf.Proc. , 135 (2003), arXiv:astro-ph/0306204[2] K. Gultekin, M. C. Miller, and D. P. Hamilton, TheAstrophysical Journal , 221 (2004), arXiv:astro-ph/0402532[3] K. Gultekin, M. C. Miller, and D. P. Hamilton, TheAstrophysical Journal , 156 (2006), arXiv:astro-ph/0509885[4] M. Iwasawa, Y. Funato, and J. Makino, The Astrophys-ical Journal , 1059 (2006), arXiv:astro-ph/0511391[5] L. Hoffman and A. Loeb, Mon. Not. R. Astron. Soc. ,957 (2007), arXiv:astro-ph/0612517[6] L. Hoffman and A. Loeb, The Astrophysical Journal ,L75 (2006), arXiv:astro-ph/0511242[7] A. Gualandris, S. Portegies Zwart, and M. S. Sipior, Mon.Not. R. Astron. Soc. , 223 (Oct. 2005), arXiv:astro-ph/0507365[8] S. Mikkola, Mon. Not. R. Astron. Soc. , 1107 (Jun.1983), http://adsabs.harvard.edu/abs/1983MNRAS.203.1107M [9] S. Mikkola, Mon. Not. R. Astron. Soc. , 115 (Mar.1984), http://adsabs.harvard.edu/abs/1984MNRAS.207..115M [10] M. C. Miller and D. P. Hamilton, The AstrophysicalJournal , 894 (2002), arXiv:astro-ph/0202298[11] P. Hein¨am¨aki, A&A , 795 (Jun. 2001)[12] M. Valtonen and S. Mikkola, Annu. Rev. Astron. Astro-phys. , 9 (1991)[13] M. Milosavljevi´c and D. Merritt, AIP Conference Pro-ceedings , 201 (2003), http://link.aip.org/link/?APC/686/201/1 [14] M. J. Valtonen, S. Mikkola, P. Heinamaki, and H. Valto-nen, Astrophys. J. Suppl. , 69 (Nov. 1994)[15] P. Amaro-Seoane and M. Dewi Freitag, ArXiv e-prints(2010), arXiv:1009.1870 [astro-ph.CO][16] J. Levin and H. Contreras, ArXiv e-prints(Sep. 2010),arXiv:1009.2533 [gr-qc][17] N. Seto and T. Muto, Phys. Rev. D , 103004 (May2010)[18] N. Seto and T. Muto, ArXiv e-prints(2011),arXiv:1105.1845 [astro-ph.CO][19] N. Stone and A. Loeb, ArXiv e-prints(2010),arXiv:1004.4833 [astro-ph.CO]%%CITATION =1004.4833;%%[20] X. Chen, A. Sesana, P. Madau, and F. K. Liu, Astrophys.J , 13 (2011), http://stacks.iop.org/0004-637X/729/i=1/a=13 [21] M. J. Valtonen and H. Karttunen, The three-body problem (Cambridge University Press, New York, 2006) ISBN 0-521-85224-2 (hardcover) [22] J. Levin, Phys. Rev. Lett. , 3515 (Apr 2000), http://link.aps.org/doi/10.1103/PhysRevLett.84.3515 [23] N. J. Cornish and J. Levin, Phys. Rev. D , 024004 (Jul.2003)[24] J. D. Barrow and J. Levin(Mar. 2003),arXiv:nlin/0303070 [nlin.CD][25] J. Levin, Phys. Rev. D , 124027 (Dec. 2006), arXiv:gr-qc/0612003[26] A. Gopakumar and C. K¨onigsd¨orffer, Phys. Rev. D ,121501 (Dec. 2005), arXiv:gr-qc/0511009[27] W. Han, Gen. Rel. Grav. , 1831 (Sep. 2008),arXiv:1006.2229 [gr-qc][28] J. D. Schnittman and F. A. Rasio, Phys. Rev. Lett ,121101 (Sep. 2001), arXiv:gr-qc/0107082[29] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 064032 (2008), arXiv:0711.1048 [gr-qc][30] J. Steinhoff, S. Hergt, and G. Sch¨afer, Phys. Rev. D ,081501 (Apr 2008), arXiv:0712.1716 [gr-qc][31] J. Hartung and J. Steinhoff, Phys. Rev. D , 044008(Feb 2011), arXiv:1011.1179 [gr-qc][32] C. Moore, Phys. Rev. Lett. , 3675 (Jun 1993)[33] T. Imai, T. Chiba, and H. Asada, Phys. Rev. Lett. ,201102 (May 2007)[34] C. O. Lousto and H. Nakano, Class. Quantum Grav. ,195019 (2008), arXiv:0710.5542 [gr-qc][35] K. Yamada and H. Asada, Phys. Rev. D , 104019 (Nov2010), arXiv:1010.2284 [gr-qc][36] K. Yamada and H. Asada, ArXiv e-prints(2010),arXiv:1011.2007 [gr-qc][37] T. Ichita, K. Yamada, and H. Asada, ArXiv e-prints(2010), arXiv:1011.3886 [gr-qc][38] J. D. Schnittman, ArXiv e-prints(2010), arXiv:1006.0182[astro-ph.HE][39] P. Galaviz and B. Br¨ugmann, Phys. Rev. D , 084013(Apr 2011)[40] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D77 , 101501(R) (2008), arXiv:0710.0879 [gr-qc][41] C. O. Lousto and Y. Zlochower, Phys. Rev.
D77 , 024034(2008), arXiv:0711.1165 [gr-qc][42] P. Galaviz, B. Br¨ugmann, and Z. Cao, Phys. Rev. D ,024005 (Jul 2010), arXiv:1004.1353 [gr-qc][43] P. Diener, Class. Quantum Grav. , 4901 (2003),arXiv:gr-qc/0305039[44] G. Jaramillo and C. O. Lousto, ArXiv e-prints(2010),arXiv:1008.2001 [gr-qc][45] M. Ponce, C. Lousto, and Y. Zlochower, ArXiv e-prints(2010), arXiv:1008.2761 [gr-qc][46] T. Szirtes, Applied Dimensional Analysis and Modeling ,2nd ed. (Butterworth-Heinemann, 2006)[47] L. Blanchet, Living Reviews in Relativity (2006), [48] M. Maggiore, Gravitational Waves. Vol. 1: Theory andExperiments (Oxford University Press, Oxford, 2007)ISBN ISBN-13: 978-0-19-857074-5[49] G. Sch¨afer, Physics Letters A , 336 (1987)[50] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 044024 (2000), gr-qc/9912092[51] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 084011 (2000), gr-qc/0005034[52] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev. D62 , 021501 (2000), erratum-ibid. 63, 029903, (2000),gr-qc/0003051[53] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Lett.B , 147 (2001), gr-qc/0105038[54] J. P. Eckmann and D. Ruelle, Rev. Mod. Phys. , 617(Jul 1985)[55] X. Wu and Y. Xie, Phys. Rev. D , 103012 (May 2008)[56] S. W. McDonald, C. Grebogi, E. Ott, and J. A.Yorke, Physica D: Nonlinear Phenomena , 125(1985), ISSN 0167-2789, [57] J. Farmer, E. Ott, and J. A. Yorke, Physica D: Non-linear Phenomena , 153 (1983), ISSN 0167-2789, [58] H.-O. Peitgen, H. J¨urgens, and D. Saupe, Chaos andFractals (Springer, 1993)[59] K. Falconer,
Fractal Geometry: Mathematical Founda-tions and Applications , 1st ed. (John Wiley & Sons, 1990)[60] V. I. Arnold,
Ordinary Differential Equations (The MITPress, Cambridge, Ma., 1978) ISBN 0262510189[61] X. Wu and T.-Y. Huang, Physics Letters A , 77(2003), ISSN 0375-9601, [62] X. Wu, T.-Y. Huang, and H. Zhang, Phys. Rev. D ,083001 (Oct 2006)[63] M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman,P. Alken, M. Booth, and F. Rossi, GNU Scientific LibraryReference Manual , 3rd ed. (Network Theory Ltd., 2009)ISBN 0954612078[64] Wolfram Research, Inc.,
Mathematica , version 7.0 ed.(Wolfram Research, Inc., 2008) [65] C. A. Burdet, Zeitschrift Angewandte Mathematik undPhysik , 434 (May 1967)[66] D. C. Heggie, Celestial Mechanics and Dynamical As-tronomy , 217 (Oct. 1974)[67] S. Mikkola and S. J. Aarseth, Celestial Mechanics andDynamical Astronomy , 375 (Dec. 1989)[68] S. Mikkola and S. J. Aarseth, Celestial Mechanics andDynamical Astronomy , 439 (Nov. 1993)[69] S. Mikkola and S. J. Aarseth, Celestial Mechanics andDynamical Astronomy , 197 (Sep. 1996)[70] H. Goldstein, C. P. Poole, and J. Safko, Classical Me-chanics (Addison Wesley, 2001) ISBN 0-201-65702-3[71] M. H´enon, Celestial Mechanics and Dynamical Astron-omy , 267 (1976)[72] M. Nauenberg, Physics Letters A , 93 (Dec. 2001),arXiv:nlin/0112003v2[73] C. Moore and M. Nauenberg, ArXiv e-prints(2008),arXiv:math/0511219v2[74] M. D. Hartl and A. Buonanno, Phys. Rev. D , 024027(Jan 2005)[75] B. Walther, B. Br¨ugmann, and D. M¨uller, Phys. Rev. D79 , 124040 (2009), arXiv:arXiv: 0901.0993 [gr-qc] [gr-qc][76] S. Husa, M. Hannam, J. A. Gonz´alez, U. Sperhake,and B. Br¨ugmann, Phys. Rev.
D77 , 044037 (2008),arXiv:arXiv:0706.0904v1 [gr-qc] [gr-qc][77] W. Tichy, B. Br¨ugmann, M. Campanelli, and P. Diener,Phys. Rev. D , 064008 (2003), gr-qc/0207011[78] E. Berti, S. Iyer, and C. M. Will, Phys. Rev. D , 061503(2006), gr-qc/0607047[79] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. Scheel, Class. Quant. Grav. , S59(2007), arXiv:gr-qc/0702106[80] R. Mardling, in The Cambridge N-Body Lectures , Lec-ture Notes in Physics, Vol. 760, edited by S. Aarseth,C. Tout, and R. Mardling (Springer Berlin / Heidelberg,2008) pp. 59–96, 10.1007/978-1-4020-8431-7 3,, Lec-ture Notes in Physics, Vol. 760, edited by S. Aarseth,C. Tout, and R. Mardling (Springer Berlin / Heidelberg,2008) pp. 59–96, 10.1007/978-1-4020-8431-7 3,