Chaos in two black holes with next-to-leading order spin-spin interactions

We take into account the dynamics of a complete third post-Newtonian conservative Hamiltonian of two spinning black holes, where the orbital part arrives at the third post-Newtonian precision level and the spin-spin part with the spin-orbit part includes the leading-order and next-to-leading-order contributions. It is shown through numerical simulations that the next-to-leading order spin-spin couplings play an important role in chaos. A dynamical sensitivity to the variation of single parameter is also investigated. In particular, there are a number of \textit{observable} orbits whose initial radii are large enough and which become chaotic before coalescence.


Introduction
Massive binary black-hole systems are likely the most promising sources for future gravitational wave detectors. The successful detection of the waveforms means using matchedfiltering techniques to best separate a faint signal from the noise and requires a very precise modelling of the expected waveforms. Post-Newtonian (PN) approximations can satisfy this requirement. Up to now, high-precision PN templates have already been known for the non-spin part up to 3.5PN order (i.e. the order 1/c 7 in the formal expansion in powers of 1/c 2 with c being the speed of light) [1,2], the spin-orbit part up to 3.5PN order including the leading-order (LO, 1.5PN), next-to-leading-order (NLO, 2.5PN) and nextto-next-to-leading-order (NNLO, 3.5PN) interactions [3][4][5], and the spin-spin part up to 4PN order consisting of the LO (2PN), NLO (3PN) and NNLO (4PN) couplings [6][7][8][9]. a e-mail: xwu@ncu.edu.cn However, an extremely sensitive dependence on initial conditions as the basic feature of chaotic systems would pose a challenge to the implementation of such matched filters, since the number of filters required to detect these waveforms is exponentially large with increasing detection sensitivity. This has led to some authors focusing on research of chaos in the orbits of two spinning black holes. Chaos was firstly found and confirmed in the 2PN Lagrangian approximation of comparable mass binaries with the LO spin-orbit and LO spin-spin effects [10]. Moreover, it was reported in [11] that the presence of chaos should be ruled out in these systems because no positive Lyapunov exponents could be found. As an answer to this claim, Refs. [12,13] obtained some positive Lyapunov exponents and pointed out these zero Lyapunov exponents of [11] due to the less rigorous calculation of the Lyapunov exponents of two nearby orbits with unapt rescaling. In fact, the conflicting results on Lyapunov exponents are because the two papers [11,12] used different methods to compute their Lyapunov exponents, as was mentioned in [14]. Ref. [11] computed the stabilizing limit values of Lyapunov exponents, and Ref. [12] worked out the slopes of fit lines. This is the so-called doubt regarding different chaos indicators causing two distinct claims on the chaotic behavior. Besides this, there was second doubt on different dynamical approximations making the same physical system have distinct dynamical behaviors. The 2PN harmonic-coordinates Lagrangian formulation of the two-black hole system with the LO spin-orbit couplings of one spinning body allows chaos [15], but the 2PN ADM (Arnowitt-Deser-Misner) -coordinates Hamiltonian does not [16,17]. Levin [18] thought that there is no formal conflict between them since the two approaches are not exactly but approximately equal, and different dynamical behaviors between the two approximately related systems are permitted according to dynamical systems theory. Seen from the canonical, conjugate spin coordinates [19], the former non-integrability and the latter integrability are clearer. As extensions, both any PN conservative Hamiltonian binary system with one spinning body and a conservative Hamiltonian of two spinning bodies without the constraint of equal mass or with the spin-orbit couplings not restricted to the leading order are still integrable. Recently, [20,21] argued the integrability of the 2PN Hamiltonian without the spin-spin couplings and with the NLO and/or NNLO spin-orbit contributions included. On the contrary, the corresponding Lagrangian counterpart with spin effects limited to the spin-orbit interactions up to the NLO terms exhibits the stronger chaoticity [22]. Third doubt relates to different dependence of chaos on single dynamical parameter or initial condition. The description of the chaotic regions and chaotic parameter spaces in [15] are inconsistent with that in [23]. The different claims are regarded to be correct according to the statement of [24] that chaos does not depend only on single physical parameter or initial condition but a complicated combination of all parameters and initial conditions.
It is worth emphasizing that the spin-spin effects are the most important source for causing chaos in spinning compact binaries, but they were only restricted to the LO term in the published papers on research of the chaotic behavior. It should be significant to discuss the NLO spin-spin couplings included to a contribution of chaos. For the sake of this, we shall consider a complete 3PN conservative Hamiltonian of two spinning black holes, where the orbital part is up to the 3PN order and the spin-spin part as well as the spin-orbit part includes the LO and NLO interactions. In this way, we want to know whether the inclusion of the NLO spin-spin couplings have an effect on chaos, and whether there is chaos before coalescence of the binaries.

Third post-Newtonian order Hamiltonian approach
It is too difficult to strictly describe the dynamics of a system of two mass comparable spinning black holes in general relativity. Instead, the PN approximation method is often used. Suppose that the two bodies have masses m 1 and m 2 with m 1 ≤ m 2 . Other mass parameters are the total mass M = m 1 + m 2 , the reduced mass µ = m 1 m 2 /M, the mass ratio β = m 1 /m 2 and the mass parameter η = µ/M = β /(1 + β ) 2 . As to other specified notations, a 3-dimensional vector r represents the relative position of body 1 to body 2, its unit radial vector is n = r/r with the radius r = |r|, and p stands for the momenta of body 1 relative to the centre. The momenta, distances and time t are respectively measured in terms of µ, M and M. Additionally, geometric units c = G = 1 are adopted. The two spin vectors are S i = S iŜi (i = 1, 2) with unit vec-torsŜ i and the spin magnitudes S i = χ i m 2 i /M 2 (0 ≤ χ i ≤ 1). In ADM coordinates, the system can be expressed as the di-mensionless conservative 3PN Hamiltonian H(r, p, S 1 , S 2 ) = H o (r, p) + H so (r, p, S 1 , S 2 ) +H ss (r, p, S 1 , S 2 ). (1) In the following, we write its detailed expressions although they are too long. For the conservative case, the orbital part H o does not include the dissipative 2.5PN term (which is the leading order radiation damping level) but the Newtonian term H N and the PN contributions H 1PN , H 2PN and H 3PN , that is, As given in [25], they are The spin-orbit part H so is linear functions of the two spins. It is the sum of the LO spin-orbit term H LO so and the NLO spin-orbit term H NLO so , i.e.
Ref. [5] gave their expressions where the related notations are and the Newtonian-looking orbital angular momentum vector is The constant terms in g and g * correspond to the LO part, and the others, the NLO part. Similarly, the spin-spin Hamiltonian H ss also consists of the LO spin-spin coupling term H LO ss and the NLO spin-spin coupling term H NLO ss , namely, The first sub-Hamiltonian reads [25] with S 0 = S + S * . The second sub-Hamiltonian is made of three parts, They are written as [7,8] Here, p 1 = −p 2 = p. In a word, the conservative Hamiltonian (1) up to the 3PN order is not completely given until Eq. (15) appears. Clearly, Hamiltonian (1) does not depend on any mass but the mass ratio. The evolutions of position r and momentum p satisfy the canonical equations of the Hamiltonian (1): The spin variables vary with time according to the following relations Besides the two spin magnitudes, there are four conserved quantities in the Hamiltonian (1), including the total energy E = H and three components of the total angular momentum vector J = L + S. A fifth constant of motion is absent, so the Hamiltonian (1) is non-integrable. 1 Its high nonlinearity seems to imply that it is a richer source for chaos. Next, we shall search for chaos, and particularly investigate the effect of the NLO spin-spin interactions on the dynamics of the system.

Detection of chaos before coalescence
With numerical simulations, we use some chaos indicators to describe dynamical differences between the NLO spinspin couplings excluded and included. The appropriate ones of the indicators are selected to study dependence of chaos on single parameter when the NLO spin-spin couplings are included. Finally, we expect to find chaos before coalescence by estimating the Lyapunov and inspiral decay times.

Comparisons
Numerical methods are convenient to study nonlinear dynamics of the Hamiltonian (1). Symplectic integrators are efficient numerical tools since they have good geometric and physical properties, such as the symplectic structure conserved and energy errors without secular changes. However, they cannot provide high enough accuracies, and the computations are expensive when the mixed symplectic integration algorithms [21,26] with a composite of the second-order explicit leapfrog symplectic integrator and the second-order implicit midpoint rule are chosen. In this sense, we would prefer to adopt an 8(9) order Runge-Kutta-Fehlberg algorithm of variable time steps. In fact, it gives such high accuracy to the energy error in the magnitude of about order 10 −13 ∼ 10 −12 when integration time reaches 10 6 , as shown in Fig. 1. Here, orbit 1 we consider has initial conditions (p(0); r(0)) = (0, 0.39, 0; 8.55, 0, 0), which correspond to the initial eccentricity e 0 = 0.30 and the initial semi-major axis a 0 = 12.2. Other parameters and initial spin angles are respectively β = 0.79, χ 1 = χ 2 = χ = 1.0, θ i = 78.46 • and φ i = 60 • , where polar angles θ i and azimuthal angles φ i satisfy the relationsŜ i = (cos φ i sin θ i , sin φ i sin θ i , cos θ i ), as commonly used in physics. The NLO spin-spin couplings are not included in Fig. 1(a), but in Fig. 1(b). It can be seen clearly that the inclusion of the NLO spin-spin couplings with a rather long expression decreases only slightly the numerical accuracy. Therefore, our numerical results are shown to be reliable although the energy errors have secular changes.
We apply several chaos indicators to compare dynamical behaviors of orbit 1 according to the two cases without and with the NLO spin-spin couplings. The method of Poincaré surface of section can provide a clear description of the structure of phase space to a conservative system whose phase space is 4 dimensions. As a point to note, it is not suitable for such a higher dimensional system (1). Fortunately, power spectra, Lyapunov exponents and fast Lyapunov in-dicators would work well in finding chaos regardless of the dimensionality of phase space.

Power spectrum analysis
Power spectrum analysis reveals a distribution of various frequencies ω of a signal x(t). It is the Fourier transformation where i is the imaginary unit. In general, the power spectra X(ω) are discrete for periodic and quasi periodic orbits but continuous for chaotic orbits. That is to say, the classification of orbits can be distinguished in terms of different features of the spectra. On the basis of this, we know through Fig. 2 that the orbit seems to be regular when the NLO spinspin couplings are not included, but chaotic when the NLO spin-spin couplings are included. Notice that the method of power spectra is only a rough estimation of the regularity and chaoticity of orbits. More reliable chaos indicators are strongly desired.

Lyapunov exponents
The maximum Lyapunov exponent is used to measure the average separation rate of two neighboring orbits in the phase space and gives quantitative analysis to the strength of chaos. Its calculations are usually based on the variational method and the two-particle method [27]. The former needs solving the variational equations as well as the equations of motion, and the latter needs solving the equations of motion only. Considering the difficulty in deriving the variational equations of a complicated dynamical system, we pay attention to the application of the latter method. In the configuration space, it is defined as [28] λ = lim where |∆ r(0)| and |∆ r(t)| are the separations between the two neighboring orbits at times 0 and t, respectively. The initial distance cannot be too big or too small, and 10 −8 is regarded as to its suitable choice in the double precision [27]. For the sake of the overflow avoided, renormalizations from time to time are vital in the tangent space. A bounded orbit is chaotic if its Lyapunov exponent is positive, but regular when its Lyapunov exponent tends to zero. In this way, we can know from Fig. 3 that orbit 1 is regular for the case without the NLO spin-spin couplings, but chaotic for the case with the NLO spin-spin couplings. Of course, it takes much computational cost to distinguish between the ordered and chaotic cases.

Fast Lyapunov indicators
A quicker method to find chaos than the method of Lyapunov exponents is a fast Lyapunov indicator (FLI). This indicator that was originally considered to measure the expansion rate of a tangential vector [29] does not need any renormalization, while its modified version dealing with the use of the two-particle method [30] does. The modified version is of the form Its computation is based on the following expression: where k denotes the sequential number of renormalization. The FLI of Fig. 4(a) corresponding to Fig. 3(a) increases algebraically with logarithmic time log 10 t, and that of Fig.  4(b) corresponding to Fig. 3(b) does exponentially with logarithmic time. The former indicates the character of order, but the latter, the feature of chaos. Only when the integration time adds up to 1 × 10 5 , can the ordered and chaotic behaviors be identified clearly for the use of FLI unlike the application of Lyapunov exponent. There is a threshold value of the FLIs between order and chaos, 5. Orbits whose FLI are larger than 5 are chaotic, whereas those whose FLIs are less than 5 are regular. The above numerical comparisons seem to tell us that chaos becomes easier when the NLO spin-spin terms are included. This sounds reasonable. As claimed in [20,21], the system (1) is integrable and not at all chaotic when the spinspin couplings are turned off. The occurrence of chaos is completely due to the spin-spin couplings, which include particularly the NLO spin-spin contributions leading to a sharp increase in the strength of nonlinearity. In fact, we employ FLIs to find that there are other orbits (such as orbits 2-5 in Table 1), which are not chaotic for the absence of the NLO spin-spin couplings but for the presence of the NLO spin-spin couplings. In addition, the strength of the chaoticity of orbits 6-8 increases. As a point to illustrate, the other initial conditions beyond Table 1 are those of orbit 1; the starting spin unit vectors of orbit 2 are those of orbit 1, and those of orbits 3-8 are θ 1 = 84.26 • , φ 1 = 60 • , θ 2 = 84.26 • and φ 2 = 45 • . Hereafter, only the dynamics of the complete Hamiltonian (1) with the NLO spin-spin effects included is focused on.

Lyapunov and inspiral decay times
Taking β = 0.5, the initial conditions and the initial unit spin vectors of orbit 1 as reference, we start with the spin parameter χ at the value 0.2 that is increased in increments of 0.01 up to a final value of 1 and draw dependence of FLI on χ in Fig. 5(a). This makes it clear that chaos occurs when χ ≥ 0.7. Precisely speaking, the larger the spin magnitudes get, the stronger the chaos gets. Note that this dependence of chaos on χ relies typically on the choice of the initial conditions, the initial unit spin vectors and the other parameters. As claimed in [24], there is different dependence of chaos on χ if the choice changes. On the other hand, taking the initial spin angles of orbit 3, fixing the spin parameter χ = 0.90 and the initial conditions (p(0); r(0)) = (0, 0.39, 0; 8.4, 0, 0), which correspond to the initial eccentricity e 0 = 0.28 and the initial semi-major axis a 0 = 11.6, we study the range of the mass ratio β beginning at 0.5 and ending at 1 in increments of 0.01. At once, dependence of FLI on β can be described in Fig. 5(b). There is chaos when β ≤ 0.86 and chaos seems easier for a smaller mass ratio.
As in the panel (a), this result is given only under the present initial conditions, initial unit spin vectors and other parameters.
Do the above-mentioned chaotic behaviors occur before the merger of the binaries? To answer it, we have to compare the Lyapunov time T λ = 1/λ (i.e. the inverse of the Lyapunov exponent) with the inspiral decay time T d , estimated by [31] and γ = 64m 1 m 2 M/5. When T λ is less than T d (or λ T d > 1), chaos would be observed. Because T λ = 3.0 × 10 3 and T d = 1.3 × 10 3 for orbit 1, the chaoticity can not be seen before the merger. Values of λ T d for orbits 2-8 are listed in Table  1. Clearly, only chaotic orbit 8 is what we expect. Besides these, we plot two panels (a) and (b) of Fig. 6 regarding dependence of Lyapunov exponent on single parameter, which correspond respectively to Figs 5(a) and 5(b). Here are two facts. First, the results in Fig. 6 are the same as those in Fig.  5. Second, lots of chaotic orbits whose Lyapunov times are many times greater than the inspiral times should be ruled out, and there are only a small quantity of desired chaotic orbits left. In order to make the accuracy of the PN approach better, we should choose orbits whose initial radii are larger enough than roughly 10M. All chaotic orbits in Table 2 are expected. Notice that the other initial conditions of these orbits beyond this table are y = z = p x = p z = 0, and the starting spin angles are still the same as those of orbit 3. Although an orbit has a large initial radius, it may still be chaotic when its initial eccentricity is high enough. This supports the result of [23] that a high eccentric orbit can easily yield chaos.

Conclusions
This paper is devoted to studying the dynamics of the complete 3PN conservative Hamiltonian of spinning compact binaries in which the orbital part is accurate to the 3PN order and the spin-spin part as well as the spin-orbit part includes the LO and NLO contributions. Because of the high nonlinearity, the NLO spin-spin couplings included give rise to the occurrence of strong chaos in contrast with those excluded. By scanning single parameter with the FLIs, we obtained dependence of chaos on the parameter. It was shown sufficiently that chaos appears easier for larger spins or smaller mass ratios under the present considered initial conditions, starting unit spin vectors and other parameters. So does for a smaller initial radius. In spite of this, an orbit with a large initial radius is still possibly chaotic if its initial eccentricity is high enough. Above all, there are some observable chaotic orbits whose initial radii are suitably large and whose Lyapunov times are less than the corresponding inspiral times.