The non-resonant, relativistic dynamics of circumbinary planets
MMon. Not. R. Astron. Soc. , 1–18 (2008) Printed 1 November 2018 (MN L A TEX style file v2.2)
The non-resonant, relativistic dynamics of circumbinary planets
Cezary Migaszewski (cid:63) and Krzysztof Go´zdziewski (cid:63) † Toru´n Centre for Astronomy, Nicolaus Copernicus University, Gagarin Str. 11, 87-100 Toru´n, Poland
Accepted 2010 September 13. Received 2010 September 10; in original form 2010 June 30
ABSTRACT
We investigate the non-resonant, 3-D (spatial) model of the hierarchical system composed ofpoint-mass stellar (or sub-stellar) binary and a low-mass companion (a circumbinary planetor a brown dwarf). We take into account the leading relativistic corrections to the Newtoniangravity. The secular model of the system relies on the expansion of the perturbing Hamiltonianin terms of the ratio of semi-major axes α , averaged over the mean anomalies. We foundthat the low-mass object in a distant orbit may excite large eccentricity of the inner binarywhen the mutual inclination of the orbits is larger than about of 60 deg. This is related tostrong instability caused by a phenomenon which acts similarly to the Lidov-Kozai resonance(LKR). The secular system may be strongly chaotic and its dynamics unpredictable over thelong-term time scale. Our study shows that in the Jupiter– or brown dwarf– mass regimeof the low-massive companion, the restricted model does not properly describe the long-termdynamics in the vicinity of the LKR. The relativistic correction is significant for the parametricstructure of a few families of stationary solutions in this problem, in particular, for the directorbits configurations (with the mutual inclination less than 90 degrees). We found that thedynamics of hierarchical systems with small α ∼ .
01 may be qualitatively different in therealm of the Newtonian (classic) and relativistic models. This holds true even for relativelylarge masses of the secondaries.
Key words: celestial mechanics – secular dynamics – analytical methods – stationary solu-tions – extrasolar planetary systems
The extrasolar planets are discovered routinely . Recently, aboutof 500 low-mass companions to stars of different spectral types areknown. Most of them are bounded to single stars. Moreover, there isalso a growing number of planetary candidates in binaries as wellas multi-stellar systems (see Eggenberger 2010, for the statisticalproperties of planets in binaries). Generally, following nomencla-ture in (Rabl & Dvorak 1988), we can consider two classes of suchmultiple systems. In the satellite case or the S-type configuration,a planet revolves around one of the primaries in the binary, and thesecond primary is much more distant. In the cometary or circumbi-nary configuration (C-type from hereafter), the planet has a wideorbit around the inner, massive binary.Current theories of planet formation in multiple stellar sys-tems (e.g., Takeda et al. 2009, and references therein) show that theinclination of a planetary orbit to the orbital plane of the binary maybe non-zero in the C-type and the S-type systems. It is well knownthat in the S-type configurations, when the mutual inclination of cir-cular orbits is larger than the critical value i crit ∼ ◦ , the inner orbit (cid:63) E-mail: [email protected]† E-mail: [email protected] undergoes large–amplitude oscillations of the eccentricity, whichis in anti-phase with the mutual inclination. This dynamical phe-nomenon is well known as the Kozai (or Lidov–Kozai) resonance(Lidov 1962; Kozai 1962). We will call it the LKR from hereafter.Keeping in mind the two types of possible orbital configurations,this instance of the LKR may be also understood as the inner LKR (see, e.g., Krasinskii 1972, 1974). Actually, many authors explainlarge eccentricities of some planetary candidates due to forcing bya distant star or massive companion (see, e.g., Tamuz, O., et al.,2008; Fabrycky & Tremaine 2007). In fact, the amplification of ec-centricity and inclination may also appear in the C-type systems.The critical inclination is then ∼ ◦ , and it may be attributed to the outer LKR (see, e.g., Krasinskii 1972, 1974; Migaszewski &Go´zdziewski 2010), and this will be addressed further in this work.In the literature, the problem is most often considered as therestricted one, which means that the planet is a mass-less particlenot perturbing the motion of the binary. It has been studied in manypapers, with different analytical and numerical techniques. The re-stricted model help us to simplify the analysis, nevertheless, theassumption of negligible influence of the planet on the motion ofprimaries may be not valid if the planet is large. In fact, low-massobjects in a few Jupiter-mass range are quite common. If the mutualinteractions are significant, as we will show further in this work, thebinary orbit may be strongly perturbed by a distant, relatively mas- c (cid:13) a r X i v : . [ a s t r o - ph . E P ] O c t C. Migaszewski and K. Go´zdziewski sive planet or a brown dwarf moving in inclined, wide orbit, even ifthe mass of inner companion is ten times larger than the perturber.In this work, we focus on the unrestricted C-type problem bymeans of the secular model in terms of the semi-major axes ratio, α . We assume that α is small ( < . ∼
100 times longer. Hence the short-term mean motionresonances (MMRs) are not present in the system. The secular evo-lution of the mean orbits, depending on the mutual interactions, iseven much longer, and spans Myr time-scale. To follow the orbitalevolution in all these time-scales in details, one could integrate theequations of motion numerically. Unfortunately, this direct, brute-force approach requires too large CPU overhead.Because we are primarily interested in the long-term evolutionof the hierarchical systems, the problem may be simplified withthe help of the averaging principle. The Hamiltonian of the hier-archical system may be splitted onto integrable Keplerian part ofthe inner binary, and for the perturbation part of the mutual inter-action with the low-mass companion. Because the system is non-resonant, the perturbing Hamiltonian may be averaged out over themean anomalies or the mean longitudes, which play the role of thefast angles. After the averaging, we obtain the secular Hamiltoniandescribing the long-term evolution of the mean, slowly varying or-bits. To obtain the secular model, we extend a simple averagingscheme (the so called mixed anomalies method ) in our earlier pa-per (Migaszewski & Go´zdziewski 2008) to the non-coplanar case.The perturbing Hamiltonian is expanded in power series with re-spect to α , and then these series may be averaged out term by term.We derived such averaged expansion of the secular Hamiltonian tothe 10-th order in α . It is a more general version of the third-order(octupole-level) theory, studied in the planetary context in (Fordet al. 2000; Lee & Peale 2003) and of the integrable, second order(quadrupole-level) approximation in many papers (e.g., Harring-ton 1968; Krasinskii 1972; Lidov & Ziglin 1976; Ferrer & Os´acar1994; Farago & Laskar 2010, and references therein).Our work is closely related to remarkable study byMichtchenko et al. (2006), to which we will refer many times, aswell as to our earlier paper (Migaszewski & Go´zdziewski 2009a)devoted to the analysis of equilibria in the three-dimensional prob-lem of two planets. Moreover, this work extends these papers intwo important aspects. The averaging of the secular problem isdone analytically, which simplifies and optimizes the computa-tions, helping us to avoid numerical artifacts. Here, we also con-sider a generalization of the Newtonian model, by accounting forthe leading non-Newtonian point-mass corrections to the perturb-ing Hamiltonian, i.e., the relativistic, post-Newtonian (PN) correc-tion. It can be also easily averaged out over the mean anomalies.The PN corrections are particularly important if the frequenciesof the slow angles, which they induce, are comparable with thefrequencies causes by the Newtonian interactions (e.g., Adams &Laughlin 2006; Fabrycky & Tremaine 2007). We shown already(Migaszewski & Go´zdziewski 2009b) that accounting for the rela-tivistic corrections in the co-planar case of two-planet system maylead to qualitatively different global view of the phase space in bothmodels. Hence, the PN correction might be also important in the 3Dproblem. Indeed, as we will show below, this apparently subtle ef-fect induces global, qualitative changes of the structure of the phasespace.It should be emphasised here, that a large number of physi-cal and orbital parameters fully characterising planetary configura- tions contradicts our desire to study the problem in possibly qual-itative way, with the help of particular geometric tools. Hence, werestrict the work to specific ranges of these parameters, focusingon a “typical” binary with relatively large mass ratio of the pri-maries, as well as the circumbinary object in Jupiter/brown-dwarfmass range. Moreover, considering corrections to the Newtonianpoint-mass gravity, we only consider the relativistic effects, which,in turn, limit the orbital parameters of the binary. The conserva-tive and dissipative tidal distortions are neglected here, though theymight dominate in compact binaries, or in configurations with veryhot-Jupiter or super-Earth planets (Fabrycky & Tremaine 2007;Mardling 2007; Batygin et al. 2009; Ragozzine & Wolf 2009;Mardling 2010). In the range of semi-major axes ∼ .
025 au, theplanetary tidal bulge raises apsidal rotation of the inner orbit whichmay may reach a few degrees per year, exceeding the effects of gen-eral relativity and the rotational stellar quadrupole by more than anorder of magnitude (Ragozzine & Wolf 2009). However, in gen-eral, as we explain below, the rotational distortions introduce extra-degrees of freedom to the model (assuming that the stellar and plan-etary spins may be arbitrary), that cannot be treated in terms of ge-ometric tools natural to investigate two-degrees of freedom Hamil-tonian dynamics. Still, although the tidal effects could be basicallytreated in this formalism too, it would introduce new parameters(bodies’ radii, Love numbers), hence we postpone investigations ofthis more general and complex model in future papers. Overall, aswe show below, in the parameter ranges investigated here (semi-major axis of the binary ∼ . . N -BODIES We consider the general, spatial model of the N -planet sys-tem around a central star. It may be described in terms of theHamiltonian written with respect to canonical Poincar´e variables(Michtchenko et al. 2006) in the form of H = H kepl + H pert , where H kepl = N ∑ i = (cid:18) p i β ∗ i − µ ∗ i β ∗ i r i (cid:19) , (1)stands for the integrable part comprising of the direct sum of the rel-ative, Keplerian motions of point-mass secondaries m i , i = ,..., N , c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets with respect to the primary mass m . We also define the mass pa-rameters µ ∗ i = k ( m + m i ) , where k is the Gauss gravitational con-stant, and β ∗ i = ( / m i + / m ) − are the so called reduced masses.The term H pert stands for the perturbing function of the Keplerianmotions. We assume that the perturbation is a sum of two terms: H pert = H NG + H GR , (2)where H NG is related to a small Newtonian mutual interactions be-tween m i and m j , and we assume that H pert / H kepl (cid:28)
1. That maybe accomplished either by keeping m i small (then we have the planetary regime ) or permitting that secondary masses are rela-tively large (even comparable with the central object) and simul-taneously requiring large separations between particular orbits (the stellar regime ). The term of H GR is for the leading general relativity(PN) corrections to the potential of the central star and the inner-most companion. Here we focus on the non-resonant systems, withwell separated orbits, hence we may neglect the relativistic post-Newtonian perturbations of the outer bodies. If the semi-major axesratios α i , j = a i / a j < . H NG per-turbation may be written as follows: H NG = N − ∑ i = N ∑ j > i (cid:18) − k m i m j ∆ i , j (cid:124) (cid:123)(cid:122) (cid:125) direct part + p i · p j m (cid:124) (cid:123)(cid:122) (cid:125) indirect part (cid:19) , (3)where r i , i = ,..., N are for the position vectors of m i relative tothe central body, p i are for their conjugate momenta relative to the barycenter of the whole ( N + ) -body system, ∆ i , j = (cid:107) r i − r j (cid:107) de-note the relative distance between bodies i and j .After (Richardson & Kelly 1988), or developing the PNHamiltonian from the general Lagrangian in (Brumberg 2007), thepost-Newtonian potential in the PN formulation, H GR ≡ β ∗ H (cid:48) GR ,where H (cid:48) GR has the following form: H (cid:48) GR = γ P + γ P r + γ ( r · P ) r + γ r , (4)with coefficients γ , γ , γ , γ : γ = − ( − ν ) c , γ = − µ ∗ ( + ν ) c , γ = ( µ ∗ ) c , γ = − µ ∗ ν c , where c stands for the speed of light in a vacuum, µ ∗ = k ( m + m ) , ν ≡ m m / ( m + m ) , P is the astrocentric momentum of theinnermost secondary (normalized through the reduced mass): P = v + c (cid:20) γ ( v · v ) v + γ r v + γ r ( r · v ) r (cid:21) , (5)and v ≡ ˙ r stands for the astrocentric velocity of the innermost ob-ject (still, assuming that the relativistic corrections from the otherbodies in the system are neglected). Hence, P = v with the accuracyof O ( c − ) and then the relativistic Hamiltonian is conserved up tothe order of O ( c − ) .It is well known that the equations of motion of the N -bodysystem with N (cid:62) mean orbital elements.To perform the averaging, the perturbation must be expressedin terms of the canonical action–angle variables ( I , φ ) : H ( I , φ ) = H kepl ( I ) + H pert ( I , φ ) , (6) and we assume that H pert ( I , φ ) ∼ ε H kepl ( I ) , where ε (cid:28) l i ≡ M i , L i = β ∗ i (cid:112) µ ∗ i a i , g i ≡ ω i , G i = L i (cid:113) − e i , (7) h i ≡ Ω i , H i = G i cos I i , where M i are the mean anomalies, a i stand for canonical semi-major axes, e i are the eccentricities, I i denote inclinations, ω i arethe arguments of pericenter, and Ω i denote the longitudes of as-cending node. The Hamiltonian of the N -planet system written interms of the these Delaunay variables (Eq. 7) has the form of: H = − N ∑ i = ( µ ∗ i ) ( β ∗ i ) L i + H pert ( L i , l i , G i , g i , H i , h i ) (cid:124) (cid:123)(cid:122) (cid:125) i = ,..., N . (8)In this formulation, l i are the fast angles. They may be eliminatedthrough the averaging that is accomplished with: H sec = ( π ) N (cid:90) π ... (cid:90) π (cid:124) (cid:123)(cid:122) (cid:125) i = ,..., N H pert d M ... d M N . (9)We should remember here that H sec is valid only if (1) H pert ∼ ε H kepl (where ε (cid:28) unperturbed Keplerian orbits is equivalent to performing the firststep of the perturbation calculus (Ferraz-Mello 2007), (2) there isno mixed resonances between the inner binary and the outer com-panion (e.g., between slow frequencies of the inner orbit and themean motion of the outer orbit). We checked that planetary systemsstudied in this paper obey these assumptions within respective pa-rameter ranges. These calculations rely on the averaged model, andwill be given below (see the end of Sect. 3). Because the secularHamiltonian H sec does not depend on mean anomalies M i , the con-jugate momenta L i are integrals of the secular problem. Obviously,the mean semi-major axes are also constant, hence we get rid of N -degrees of freedom (DOF). N -BODIES In (Migaszewski & Go´zdziewski 2008), we describe a simplescheme of the averaging the perturbing function H pert (Eq.2) incoplanar case, which makes use of the very basic properties ofthe Keplerian motion, the mixed anomalies algorithm. This methodmay be easily applied to the 3D problem. At first, we consider the direct part of the mutual interaction between the planets, H NG (Eq.3). The indirect part averages out to a constant and does not con-tribute to the secular dynamics (Brouwer & Clemence 1961).The secular Hamiltonian may be written as a sum of termsrepresenting mutual interactions between all pairs of bodies i < j ,where a i < a j : (cid:10) H NG (cid:11) = N − ∑ i = N ∑ j > i (cid:68) H ( i , j ) NG (cid:69) . (10)For a particular pair of planets, we calculate the following integral: (cid:68) H ( i , j ) NG (cid:69) = ( π ) (cid:90) π (cid:90) π − k m i m j ∆ i , j d M i d M j . (11)Hence, the problem may be reduced to averaging the inverse of c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski the distance ∆ i , j between two particular planets over their meananomalies: ∆ i , j = (cid:113) r i + r j − r i · r j , (12)where r i and r j must be expressed in a common reference frame F . The same vectors, written in the orbital reference frames F i ofeach planet, are r i (cid:12)(cid:12) F i = [ x i , y i , ] T , and expressed in the commonreference frame, they have the form of: r i = A i r i (cid:12)(cid:12) F i . Here, the rota-tion matrix A i ≡ A i ( ω i , Ω i , I i ) is the matrix product of elementaryEulerian rotations (Murray & Dermott 2000): A i ( ω i , Ω i , I i ) = P ( − Ω i ) P ( − I i ) P ( − ω i ) . Formulae 12 may be rewritten as follows: ∆ i , j = r j (cid:115) − r j r i · r j r j + (cid:18) r i r j (cid:19) . (13)Following (Migaszewski & Go´zdziewski 2008), we express the ra-dius vector r i of the inner body in each planetary pair with respectto the eccentric anomaly. The radius vector of the outer body in thepair is parametrised by the true anomaly. This choice implies that ∆ − i , j expanded in Taylor series with respect to α is a trigonometricpolynomial. To compute the integral in Eq. 11, we also change theintegration variables: d M i ≡ I i dE i , and d M j ≡ J j d f j , where auxiliary functions appear: I i = − e i cos E i , J j = ( − e j ) / ( + e j cos f j ) . (14)Finally, the averaged mutual perturbation has the same generalform as in the coplanar case (Migaszewski & Go´zdziewski 2008): (cid:68) H ( i , j ) NG (cid:69) = − k m i m j a j × (15) × (cid:34) + (cid:113) − e j ∞ ∑ l = X li , j R ( i , j ) l ( e i , e j , ω i , ω j , Ω i , Ω j , I i , I j ) (cid:35) , although explicit expressions for R ( i , j ) l are obviously different inthe 3D model. The zeroth-order term in Eq. 15 is reduced to aconstant and does not influence the secular equations of motion.Also the first order term vanishes identically. The remaining terms R ( i , j ) l have rather complex form. In the simplest three-body sys-tem ( i ≡ , j ≡ ∆Ω = ± π and G sin I = G sin I (see e.g.,Michtchenko et al. 2006). It is also natural to introduce the mutualinclination, i mut ≡ I + I . Then the R ( i , j ) l -terms of the order of 2and 3 may be identified with the quadrupole and octupole terms,respectively (see, e.g., Ford et al. 2000; Lee & Peale 2003; Farago& Laskar 2010). The quadrupole-level term is the following: R ( , ) = D (cid:16) + e (cid:17) − e D cos2 ω , (16)where C I ≡ cos i mut , and D = ( C I − ) / , D = C I − . The third order (octupoloe-level) term reads as follows: R ( , ) = − D e e cos ∆ϖ (cid:16) e + (cid:17) − D D e e cos ( ω − ω ) − D D e e cos ( ω + ω ) (17) + D e e cos ( ω + ω )+ D e e cos ( ω + ω ) , where coefficients D j are: D = ( + C I ) / , D = − C I , D = − C I + C I + C I − , D = ( C I + C I − C I − ) / . Equations 16 and 17 are written similarly to terms appearing in thecoplanar model [see equations (22) and (23) in (Migaszewski &Go´zdziewski 2008)]. Clearly, if i mut = D = D = D = D = D = D =
0, and formulae R ( i , j ) , R ( i , j ) coincide with thoseones of in the coplanar problem.An explicit expansion of H sec shows that the quadrupole-orderterm in α introduces the evolution of eccentricity e , and in this ap-proximation, the outer eccentric e is constant. The variation of theouter eccentricity may be only introduced through the third order(octupole) and higher terms. Indeed, up to the quadrupole approxi-mation, the secular Hamiltonian does not depend on ω (the cyclicangle), and the eccentricity of the outer body becomes an additionalintegral of motion. In this case, the problem can be reduced to oneDOF and is integrable (Lidov & Ziglin 1976). This feature has beenaccounted for in many recent papers, moreover, the apparently sub-tle third-order perturbation to the Keplerian model, or the first orderperturbation to the integrable quadrupole-order approximation mayintroduce qualitative changes of the dynamics.We calculated the secular expansion (Eq. 15) up to the 10-thorder . One should be aware that by increasing the order of thisexpansion, we do not necessarily improve the approximation of thesecular model of the real system, because this model is still lim-ited by the first order perturbation theory. In Section 3.2 we willexamine the accuracy of the secular expansion in more details.Finally, the averaged relativistic correction possesses the sameform as in coplanar case (Migaszewski & Go´zdziewski 2008).Moreover, we include this perturbation only to the mutual inter-action of masses m and m : (cid:10) H GR (cid:11) = − ( µ ∗ ) ( β ∗ ) c L G + const , (18)as it was explained above.Having the averaged model in hand, we may calculate the sec-ular frequencies of the inner companion, and compare them withthe mean motion of the outer object ( n ). For the relativistic ad-vance of the inner periastron we have: f , rel n = µ ∗ c a (cid:115) µ ∗ µ ∗ α − / − e , This expansion is available on request in the form of a raw MATHE-MATICA input file; it will be also available on-line, after publishing thismanuscript. c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets and for the apsidal motion forced by the mutual interaction of theinner and outer companion (in the quadrupole approximation): f , mut n = (cid:115) µ ∗ µ ∗ m m α / (cid:113) − e Λ + m m α ( − e ) Λ , where Λ , are the following functions of the geometric elements: Λ = ( − e ) (cid:104) ( C I − ) − ( C I − ) cos2 ω (cid:105) + C I (cid:104) ( + e ) − e cos2 ω (cid:105) , Λ = C I (cid:104) ( + e ) − e cos2 ω (cid:105) . Assuming now that µ ∗ , ∼ k m (the central mass dominates), wemay obtain the following estimates of the secular frequencies interms of the characteristic units in our model: the relativistic fre-quency relative to n is f , rel n = × − (cid:18) m m (cid:12) (cid:19)(cid:18) . a (cid:19)(cid:18) . α (cid:19) / − e , while the mutually forced frequency relative to n : f , mut n = × − (cid:18) m m J (cid:19)(cid:18) m (cid:12) m (cid:19)(cid:16) α . (cid:17) / × ( − e ) / (cid:113) − e Λ + × − (cid:18) m m J (cid:19)(cid:18) m (cid:12) m (cid:19)(cid:16) α . (cid:17) ( − e ) Λ . These frequencies may be compared with the tidal apsidal fre-quency induced by the primary and the inner body bulge (Mi-gaszewski & Go´zdziewski 2010, in preparation): f , tid n = (cid:26) × − (cid:18) m m J (cid:19)(cid:18) m (cid:12) m (cid:19)(cid:18) R R (cid:12) (cid:19) (cid:18) k L , . (cid:19) + × − (cid:18) m m (cid:12) (cid:19)(cid:18) m J m (cid:19)(cid:18) R R J (cid:19) (cid:18) k L , . (cid:19)(cid:27) × (cid:18) . a (cid:19) (cid:18) . α (cid:19) / + e / + e / ( − e ) , where R is the stellar radius, R is the radius of the inner sec-ondary and k L , , k L , are tidal Love numbers of these bodies. Letus choose a reference model through setting characteristic parame-ters of a ∼ . α ∼ . m ∼ m (cid:12) , m ∼ m J , R ∼ R (cid:12) , R ∼ R J . Assuming that the bodies are modeled by polytropes withindices of 3 and 2, respectively, we compute their Love numbers, ∼ .
03 for the primary, and ∼ .
15 for the inner secondary. Thensetting e ∼
0, we obtain that the relativistic frequency is compara-ble with the mutual, Newtonian frequency, while the mean motionof the outer secondary is orders of magnitude larger than both ofthem (hence no mixed resonances are possible). Simultaneously,the assumptions of the averaging principle are well fulfilled. Thisguarantees that the evolution of the mean (secular) system closelyfollows the real configuration over the time scale of order ∼ / ε ,where ε is the small parameter of the perturbation.Under the same assumptions, the tidal frequency is ordersof magnitude smaller than the leading frequencies of the mutual(Newtonian) and relativistic corrections. This means that the tidaleffect is negligible, indeed, as far as the model parameters do notstrongly deviate from the characteristic values, as defined above. Because the general planetary N -body problem is very complex, werestrict the further analysis to its simplest, non-trivial case of threebodies (“non-trivial” in the sense of its non-integrability). We shallconsider configurations of the host star and two planets or C-typesystems comprising of a binary and a more distant body (a planet).We recall that the secular Hamiltonian H sec of the three bodyproblem does not depend on M , M , therefore the conjugate ac-tions ( L , L ) are constant. The Hamiltonian H written in theLaplace reference frame depends on ∆Ω = Ω − Ω ≡ ± π only,not on Ω and Ω separately. Hence, the following canonical trans-formation (e.g., Michtchenko et al. 2006): ( ω , G )( ω , G )( Ω , H )( Ω , H ) ⇒ ( ω , G )( ω , G ) (cid:0) θ = ( Ω + Ω ) , J = H + H (cid:1)(cid:0) θ = ( Ω − Ω ) , J = H − H (cid:1) (19)removes Ω , Ω from the secular Hamiltonian. After this transfor-mation it does not depend on θ , therefore J ≡ | C | = C = const,where C is the total angular momentum of the system. Moreover, θ = ± π / = const (after the Jacobi reduction of nodes) and J may be expressed as a function of G , G and J in the followingform: J = ( G − G ) / J . Therefore, for constant values of the angular momentum J ≡ C andthe secular energy H sec , the secular dynamics are reduced to twoDOF Hamiltonian system. Instead of the total angular momentum C , we shall use the so called Angular Momentum Deficit ( AMD )(Laskar 2000):
AMD = L + L − C , or its normalized value of A ≡ AMD / ( L + L ) , A ∈ [ , ] (Mi-gaszewski & Go´zdziewski 2009a). Because L , L and C are inte-grals of the secular system, the relativistic correction, Eq. 18 doesnot change A , and the DOF number does not change. Because (cid:10) H GR (cid:11) depends on G only, thus it affects only the temporal evo-lution of ω . We note here that the perturbation induced by thequadrupole moment of the star, which was discussed (Migaszewski& Go´zdziewski 2009b) in the coplanar case, also depends on H inthe 3D problem, i.e. on the orbital inclination to the equatorial planeof the star. This perturbation introduces an additional frequency to Ω and then ∆Ω is not constant anymore. It means that we couldnot perform the reduction of nodes and the secular problem wouldhave three DOF. This also means that the Laplace reference frame,defined in terms of the total orbital angular momentum, does notpossess any constant orientation in space. Being aware of this prob-lem, we do not consider the dynamical flattening of the star and/orof the innermost planet. The two DOF model is then less general butthe dynamics are better tractable, thanks to the geometrical tools,which can be applied to study this basic, low-dimensional problem.If we fix the secular Hamiltonian in the form of H sec ≡ H sec ( G , G , ω , ω ) , then A may be considered as a free param-eter of the secular model. Moreover, the phase space is four-dimensional , and to represent the phase space of the system glob-ally in terms of two-dimensional sections, which are easy to visual-ize, we follow a concept of the representative plane of initial condi-tions (Michtchenko & Malhotra 2004), the Σ -plane from hereafter.The Σ -plane may be chosen in different ways, although all repre-sentations may be fixed and defined through the following condi- c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski tions: ∂ H sec ∂ω = , ∂ H sec ∂ω = . (20)These conditions imply that all phase trajectories of the secular sys-tem cross the Σ -plane (Michtchenko et al. 2006; Libert & Henrard2007), see also our explanation in (Migaszewski & Go´zdziewski2009a). In accord with the symmetries in the secular 3D model, thesolutions to these equations are the following four pairs of angles: ( ω , ω ) = { ( , ) ; ( , ± π ) ; ( ± π / , ± π / ) ; ( ± π / , ∓ π / ) } , (21)that also define four distinct quarters of the Σ -plane, numbered withRoman numbers II, I, IV, and III, respectively, see (Michtchenkoet al. 2006) for details. Let us note that no other combinations of theangles are permitted by Eqs. 20. This feature of the secular systemflows from the symmetry of the secular Hamiltonian with respectto the apsidal lines of the mean orbits, and may be also justifiedby the explicit form of the equations of motion derived from theexpansion of the perturbing Hamiltonian, see (Michtchenko et al.2006; Migaszewski & Go´zdziewski 2009a) for details.In this sense, the Σ -plane may be thought as an analogue ofthe Poincar´e cross section. The conditions fixing the characteristicplane may be also rewritten as follows:cos ω = cos ω = ∪ sin ω = sin ω = . Further, we shall use three, basically equivalent geometric repre-sentations of the Σ -plane, which cover certain combinations of thequarters (the solution pairs of the pericenter arguments): • the P S -plane is defined with cos ω = cos ω =
0, and P S = { x = e sin ω , y = e sin ω , e ∈ [ , ) , e ∈ [ , ) } , (22) • the P C -plane is defined with sin ω = sin ω =
0, and P C = { x = e cos ω , y = e cos ω , e ∈ [ , ) , e ∈ [ , ) } , (23) • and, finally, two incarnations of the S -plane: S = { x = e cos ∆ϖ , y = e cos2 ω , e ∈ [ , ) , e ∈ [ , ) } , (24) S ’ = { x = e cos2 ω , y = e cos ∆ϖ , e ∈ [ , ) , e ∈ [ , ) } , (25)that was defined originally in (Michtchenko et al. 2006). It can beshown, that the P S - and P C -planes carry out the same informationas the S -plane. However, the S -plane has a discontinuity along the x -axis, and the former representations are sometimes more conve-nient to the analysis of solutions of the secular system. We left a test of the accuracy of the secular expansion to this end,because the introduced Σ -planes are useful to illustrate the resultsof this test in a global manner. We select initial conditions in the S -plane, and the secular energy computed with the help of the an-alytic expansion is compared with the results of numerical averag-ing developed in (Gronchi & Milani 1998; Michtchenko & Mal-hotra 2004), which are exact up to the numerical quadrature error.We consider the non-relativistic case only, because the secular rel-ativistic correction is exact (with the first non-zero post-Newtonianterm included), and it does not influence the precision of analyticformulae. Figure 1 shows the levels of H sec , marked with solidcurves in the S -plane. Each panel is for a different order of theanalytic approximation. The relative difference between values of the mean Hamiltonians, derived through the analytic (“A”) and nu-merical (“N”) algorithms may be defined as follows: ∆ l ≡ (cid:13)(cid:13)(cid:13)(cid:13) H A sec ( l ) − H N sec H N sec (cid:13)(cid:13)(cid:13)(cid:13) , (26)where l is the order of the analytic expansion in α . The resultsof this comparison are illustrated in Fig. 1 that shows the lev-els of ∆ l computed for a hierarchical system with α = .
04 and µ ≡ m / m =
5. The quadrupole-level model reproduces the sec-ular Hamiltonian as the even function with respect to both x ≡ e { sin , cos } ω and y ≡ e { sin , cos } ω . The higher order approx-imations of H sec broke this symmetry. We have shown in (Mi-gaszewski & Go´zdziewski 2009a) that the shape of H sec signifi-cantly depends on α . This is more important in the spatial prob-lem, because even for relatively small α , the quadrupole model dis-torts the structure of the phase-space (see Sect. 4.1). An inspectionof Fig. 1 reveals that the octupole model reconstructs the secularHamiltonian and its shape in the S -plane very well, because thelargest deviations ∆ ∼ − appear only for e ∼
1. In other partsof the representative plane, ∆ ∼ − . The high-order expansionsare obviously even more precise. Following estimates of the secu-lar frequencies in the relevant parameter ranges (see Sect. 3), theaveraging principle assures us that the orbital evolution of the sec-ular system follows closely the “real” orbits. Hence, we may quitesafely skip a comparison of the results from both approaches bydirect numerical integrations.The octupole model introduces the first order perturbation tothe integrable quadrupole model, hence, because it is very precisein the range of small α , we further focus on this most simple, non-trivial case. In this work, we focus on the simplest class of solutions that arethe equilibria (or stationary solutions). These solutions imply thebasic structure of the phase space. By determining their stability,we may derive the local structure of neighboring phase trajecto-ries through relatively simple analysis. The equilibria of the non-resonant system in terms of the quadrupole approximations werestudied in the past, (e.g., Krasinskii 1974; Lidov & Ziglin 1976).In our earlier work (Migaszewski & Go´zdziewski 2009a), we clas-sified families of equilibria emerging in the two-planet, non-planarproblem, in terms of the quasi-analytic averaging, and basically ex-act secular model. We studied a few families of equilibria knownin the literature, e.g., the zero-eccentricity solutions and the Lidov-Kozai resonance. We also found new families of these equilibria,in particular, the so called chained orbits solution that could behardly derived with the perturbative approach, although the resultsin Gronchi & Milani (1998) might be applied here. This work fo-cus on the planetary regime of parameters µ , α , the mass ratio µ was restricted to the range of [ . , ] and α ∈ [ . , . ] . More-over, we explored the whole permitted range of A ∈ [ , ] . Yet welearned that the semi-analytic approach has serious disadvantages.All calculations, including the integrations of the equations of mo-tions, and the stability analysis must be performed with the help ofnumerical algorithms. This may introduce large errors and hindersthe analytic, qualitative analysis of the problem.In this section, we extend the study of stationary solutions in(Migaszewski & Go´zdziewski 2009a) to a wider range of the massratio, µ >
2. Simultaneously, we consider smaller α < .
1. Thatmakes it possible to apply the analytic model described and testedin Sect. 3. Because planets emerge most likely from remnants of a c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets Figure 1.
The precision of the analytic theory in terms of ∆ l , Eq. 26, which is color-coded in the S -plane. Parameters of the planetary system are as follows: α = . , µ = , A = .
32. Panels from (a) to (d) are for the expansions in α of the second, third, sixth, and tenth order, respectively. Figure 2.
Levels of H sec at the P S -plane calculated for α = .
01 and µ =
20. Panels from a) to f) are for different values of A = { . , . , . , . , . , . } , respectively. These values of A imply the mutual inclination at the origin i = { ◦ , ◦ , ◦ , ◦ , ◦ , ◦ } , respec-tively. The intensity of shaded ares encodes the mutual inclination in prescribed ranges, darker shade is for larger i mut . The inclinations ranges are i mut = ( ◦ , ◦ ) , ( ◦ , ◦ ) , ( ◦ , ◦ ) , ( ◦ , ◦ ) , ( ◦ , ◦ ) , ( ◦ , ◦ ) in subsequent panels. thin protoplanetary disk, we also restrict our attention to mutual in-clinations up to i mut ∼ ◦ (direct orbits). In that range, we may findfamilies of equilibria related to the zero-eccentricity orbits, and theLK resonance classified in (Migaszewski & Go´zdziewski 2009a)as solutions of family IVa, accompanied by families IIIa, IIIb, andIVb+. Our analysis is also restricted to the initial conditions in quar-ters III and IV of the S -plane, in which these solutions may only“reside”. This region of the parameter space has been studied (e.g.,Krasinskii 1974) with the quadrupole-level model which is some-how the next to trivial, non-interacting Keplerian approximation ofthe three-body orbits. However, as we show below, this approxima-tion may introduce artifacts due to generally non-realistic symme-try of the secular Hamiltonian. Obviously, to avoid this problem, higher order expansions are required. We also attempt to extend theresults of Ford et al. (2000), who applied the octupole-level theoryto the analysis of hierarchical planetary systems with very small α ∼ . H sec at the P S -plane for α = . µ =
20 and a few different values of A . Panel 2a was derivedfor A = .
08, and it reveals the zero-eccentricity equilibrium. Forlarger A = .
12 (Fig. 2b) the figure-eight structure of Lidov–Kozaiequilibrium appears, and is labeled with IVa. We recall, that themutual inclination of circular orbits that corresponds to the LKbifurcation will be called the critical inclination , i crit although, in c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
Figure 3.
Families of stationary solutions (IVa, IIIa, IIIb, IVb+). Semi-major axes ratio α = .
01, mass ratio µ =
20. Dark dots are for stable equi-libria, grey dots are for unstable equilibria. Crossed and dotted circles markbifurcations. Small arrows shows the direction of particular stationary solu-tion with increasing A . Green dots are for the positions of equilibria for par-ticular values of A = . , . , . , . , .
25 (the energy levels are shownin Fig. 2). Solutions are for classic model, and were obtained with the helpof the octupole theory. general, any such value of the mutual inclination that leads to abifurcation of equilibria has the sense of being “critical” (Krasin-skii 1972, 1974). For larger value of A = .
16 (Fig. 2c) the LKR“moves” towards larger e (simultaneously, e ∼ A = .
19 (Fig. 2d), these structures still expand, andfor A = . P S -plane, nevertheless it does not correspond to abifurcation of this equilibrium, see (Migaszewski & Go´zdziewski2009a). Solution IVa moves towards larger e , see Fig. 2f).The parametric paths of these equilibria in terms of A , maybe depicted in the P S -plane (Fig. 3). Black and grey curves are forthe stable and unstable equilibria, respectively. The relevant fami-lies of equilibria are labeled with Roman numbers and Latin letters.The direction of “motion” of particular solutions along the A -axisis marked with arrows. For a reference, positions of the equilibriafor a few discrete values of A = . , . , . , . , .
25 (corre-sponding to subsequent panels in Fig. 2), are marked with greendots and labeled. Following a particular evolution path of equilib-rium IVa, we see that it appears for A ∼ . A increases, this solution moves along e ∼ e . For A ∼ .
17, it reaches the maximal e ∼ . e with simultaneous increase of e . When A increases even more, this solution reaches the borderof convergence of the analytic expansion. We call that border theanti-collision line [see (Migaszewski & Go´zdziewski 2008)].The parametric evolution of equilibria in quadrant III is morecomplex. The first non-trivial stationary solution appears for A ∼ .
15, which is unstable, saddle point IIIa. It evolves along e ∼ e increases to large values. For A ∼ . ( e ∼ . , e ∼ . ) . When A increases, the solution IVb+ movestowards e → e →
0. Simultaneously, equilibrium IIIb isshifted towards increasing e and decreasing e . Then also IIIa in- creases e and leaves off the e = A ∼ . ( e ∼ . , e ∼ . ) . Following (Michtchenko et al. 2006), we shown in (Migaszewski& Go´zdziewski 2009a) that the averaged 3D model may involvestrongly chaotic motions, in the secular time-scale. To study in de-tails, how the model parameters influence the structure of the S -plane, and how it relates to the long-term chaotic phenomena, weapply the fast indicator approach. Among many numerical toolsof this kind, we choose the so called coefficient of the diffusion offundamental frequencies σ introduced by (Laskar 1990). To checkwhether a phase space trajectory of a quasi-integrable Hamiltoniansystem is regular (quasi-periodic) or irregular (chaotic), one inte-grates the equations of motion over two subsequent intervals oftime, e.g., [ , T ] and [ T , T ] . Next, we resolve the frequencies inthe discrete orbital signal with the help of refined FFT analysis(Laskar 1990), obtaining two estimates of a given frequency, say f T and f T , over these two intervals of time. The coefficient of thediffusion of the fundamental frequency is then defined through: σ = (cid:13)(cid:13)(cid:13)(cid:13) f T f T − (cid:13)(cid:13)(cid:13)(cid:13) . Clearly, if the signal does not change over time, σ ∼ σ issignificantly different from 0, the trajectory is chaotic and regardedunstable. In our calculations, we used a variant of the frequencyanalysis developed by (ˇSidlichovsk´y & Nesvorn´y 1996), which iscalled the Frequency Modified Fourier Transform (FMFT). We alsoused publicly available code of the FMFT algorithm, kindly pro-vided by David Nesvorn´y on his personal web-page .Because the secular evolution is associated with ω i angles,we compute the σ coefficient on the basis of complex time-series { G i ( t ) expi ω i ( t ) } ,... i = , , where i is the imaginary unit. In thissignal, the osculating eccentricity and pericenter argument for eachplanet represent rescaled canonical action-angle variables. Hence,resolving its Fourier components, we may determine the leadingamplitudes (proper eccentricities) and the fundamental frequenciesof pericenter angles.Next, we did massive integrations of the secular equations ofmotion. The initial conditions were selected at the grid of 200 × P S -plane. At each point of the dynamical map, weintegrated the secular trajectory over ∼ secular periods withrespect to the smaller frequency (typically, one of the fundamentalfrequencies is much larger than the other one). Having the com-puted phase trajectory, we then find an estimate σ , as well as themaximal eccentricities attained by both orbits during the integra-tion time (the so called max e indicator) as well as the amplitudeof variation of the mutual inclination ∆ i ≡ ( max i mut − min i mut ) at-tained during the integration time. These geometrical characteris-tics are very useful to understand the source of instabilities indi-cated and detected by the diffusion coefficient σ .Figures 4–6 illustrate the results derived for the same systemthat energy levels are shown in Fig. 2. Subsequent figures are for A = . , . , .
2, respectively (see Figs. 2b,c,e for the respec-tive levels of H sec ). (We note, that for A = .
08 the view of thephase space is basically very simple and the motions are regular ∼ davidn/ c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets Figure 4.
Dynamical maps for the classic (point-mass Newtonian) model, shown in the P S -plane. Panels from a) to d) are for the coefficient of the diffusionof fundamental frequencies ( σ ), maximal eccentricity of the inner and outer planet (max e , max e ), respectively, and the amplitude of variation of the mutualinclination ( ∆ i ). Dots mark areas of librations of ω (panel b), ω (panel c), and ∆ϖ (panel d). Stationary solutions are marked with circles: dotted circles arefor stable equilibria (corresponding to the maximum of the secular Hamiltonian), crossed circles are for unstable equilibria, empty circles are for the linearlystable equilibria. These solutions are labeled, in accord with (Migaszewski & Go´zdziewski 2009a). Parameters of the system: α = . , µ = , A = .
12, seealso Fig. 2b.
Figure 5.
Dynamical maps for the Newtonian, point mass model in the P S -plane, A = .
16. See the caption to Fig. 4 for more details, and also Fig. 2c. everywhere). The right-hand panel of each figure is for σ . The dy-namical character of the phase-space trajectories is color-coded:black color means quasi-regular evolution of the secular system( σ ∼ − –10 − ), and yellow colour is for strongly chaotic mo-tions ( σ (cid:62) − ). The σ -maps reveals that almost the whole phasespace is filled up with regular motions, and chaos appears only insome small regions in the P S -plane. For a reference, the coordinates of equilibria IVa, IIIa, IIIb, IVb+ and 0 (the equilibrium at the ori-gin) are marked with circles. Dotted circle is for Lyapunov stablesolution IVa, crossed circles are for unstable solutions 0, IIIa andIIIb, respectively. Equilibrium IVb+ is linearly stable and is markedwith open circle. Clearly, trajectories close to unstable equilibriaIIIa and 0 are chaotic. Let us note that unstable equilibrium IIIa liesin a very narrow, chaotic zone too. c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
Figure 6.
Dynamical maps for the Newtonian, point-mass model, in the P S -plane, A = .
20. See the caption to Fig. 4 for more details, also Fig. 2e.
Besides dynamical maps for σ , Figs. 4–6 illustrate geometriccharacteristics or orbits in the P S -plane, in terms of max e , indica-tor for the inner and the outer orbit, respectively, and for the ampli-tude of variation of the mutual inclination, ∆ i . In all these dynam-ical maps, the green dashed contours mark the mutual inclinationequal to the critical inclination of the LKR bifurcation (Krasinskii1972, 1974). Also regions, in which the secular angles ω , ω , and ∆ϖ oscillate, are shown: the gray/white/black dots surrounded byappropriate curves bound initial conditions that lead to librations of ω around ± π / e -maps, librations of ω around ± π / e -maps, and librations ∆ϖ around 0 , π in the ∆ i -maps.We note that ω librates in almost all regular trajectories with theinitial i mut > i crit . If this condition of the Lidov-Kozai mechanismis fulfilled then ω librates around ± π /
2. In chaotic zones of the S -plane, these angles do not oscillate, even if condition i mut > i crit holds true. Angle ω librates around π / e ∼
0, for A = . , .
16. This region extends significantly for larger A = .
2. In such a case, there are three zones in which both angles ω and ω librate: a region connected to equilibrium IVa, the neighborhoodof linearly stable solution IVb+, and a region of small e (thoughthere is no stationary solution associated with these librations). Thelater area is surrounded by strongly chaotic motions shown in the σ -plane.Clearly, the equilibria permitted by relatively large A affectthe secular dynamics, which become very complex. This may bebetter seen in the max e -maps. Subsequent panels in Figs. 4–6 showthat the inner eccentricity may be strongly excited if i mut > i crit . Forrelatively small A = .
12, the inner orbit, which is initially quasi-circular, becomes moderately eccentric, with max e ∼ .
3. Notethat in this specific case, the inner body is 20 times more mas-sive then the outer companion. Simultaneously, as Fig. 4 shows, theouter eccentricity e is not amplified. We may conclude that in theregime of the outer LKR, when the inner body is much more mas-sive than the outer body, the secular evolution may lead to strongperturbations of the inner orbit. Hence, even the mass hierarchy isreversed, a strong excitation of the inner eccentricity is still possi- ble, similarly to the LKR in the restricted problem. This dynamicalphenomenon may explain a large eccentricity of the binary, whena distant, low-mass and dark companion cannot be detected due tovery long orbital period.The dynamical maps in Figs. 4–6 reveal some zones, in whichthe outer eccentricity is strongly amplified. This eccentricity is ob-viously constant in terms of the quadrupole model. This feature ofthe octupole model shows that the small perturbation may causeextended, geometric changes of the mean orbits. The amplifica-tion of e , with simultaneously almost constant relative inclination, ∆ i ∼ const (see Figs. 6c,d) seems appear due to the bifurcation ofequilibria IVb+ and IIIb for A ∼ . The origin of chaotic secular dynamics may be explained by thepresence of separatrices, which encompass different types of mo-tions, librations and rotations of angles ω , ω and ∆ϖ . A classifi-cation of these libration modes in the secular problem of two plan-ets modes has been introduced in (Michtchenko et al. 2006). Theseparatrices appear due to equilibria and periodic solutions in theintegrable, or close to integrable secular models, which might beunderstood as the quadrupole approximation and/or the co-planarconfiguration.Here, we found a simple explanation of the mechanism gen-erating chaos, which is in fact the same as in the perturbed pen-dulum. To show this, let us recall that the expansion of H sec tothe second order in α is the celebrated integrable quadrupole-levelmodel. The energy levels of this model are illustrated in Fig. 7a.Because e is the integral of motion, the representative plane maybe constructed similarly to the integrable co-planar problem, S (cid:48)(cid:48) ≡{ e cos2 ω × e } . In the region marked in green colour, any con-stant level of e crosses the energy curve in two turning points limit-ing the range of variation of e . If the e level is tangent to the givenlevel of the energy, the dynamics must be then confined to fixed e ,hence we obtain an equilibrium point (stable or unstable). A set of c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets these equilibria for increasing, fixed values of e is marked by twothick curves that meet around of ( e = , e ∼ . ) . This is the bi-furcation point, at which two families of equilibria emerge — thestable branch of the LKR and the unstable origin ( e = , e = ) .For a given value of constant e , we may then find the max-imum range of e , which corresponds to librations of the canoni-cal angle ω , and this value determines the position of separatrix,see Fig. 7b for the separatrix shown in the phase diagram in the { e cos2 ω × e sin2 ω } -plane, corresponding to the energy levelmarked by the blue curve on the levels diagram (Fig. 7a). Collectingsuch points along increasing e , we construct the red dashed curvein Fig. 7a, showing the separatrix region globally in the whole e -range (see the shaded area in Fig. 7a). The dashed line that marksthe unstable origin, might be also understood as the other branch ofthe global separatrix, which, for any value of e , corresponds to theequilibrium point.Now, considering the perturbed quadrupole-model in terms ofthe octupole-level secular expansion, we may expect that chaoticmotions will appear close to the separatrix, due to the perturbation.Indeed, this is illustrated in two panels of Fig. 8. The left panel isfor the frequency diffusion σ computed for initial conditions in theS (cid:48)(cid:48) ≡ { e cos2 ω × e } -plane. This plane shows the energy levelsof the octupole model, with stable (thick, solid curve) and unstable(thick, dashed line) equilibria in the quadrupole model over-plotted.The red curve is for the global separatrix of the LK resonance in thequadrupole model, and the green region is for the librations of ω in the integrable case. Clearly, chaotic orbits are found along thebranch of unstable equilibria, as well as close to the red separatrixcurve. This may be better seen in Poincar´e cross section computedfor a fixed energy level marked with blue curve that is shown inthe right-hand panel of Fig. 8b. This panel in fact comprises two sections ∆ϖ = π in the { e cos2 ω × e sin2 ω } -plane, for ˙ ∆ϖ > ∆ϖ < ∆ϖ around π and simultaneous circulationof angle ω , which is classified as mode in (Michtchenko et al.2006, see their Fig. 3) . Further, the plots of fundamental frequen-cies computed along the x -axis of the S ’-plane in Fig. 8a, which areshown in Fig. 9, reveal that indeed, in this region the fundamentalfrequencies decrease to 0, indicating the separatrix region.This example shows clearly that apparently very small per-turbation (recall that α = .
01) to the integrable model introducesextended chaotic behavior and qualitative change of the secular dy-namics of the system.This interpretation is helpful to understand the source ofchaotic motions shown in all dynamical maps (e.g., Figs. 4–6) –usually, they appear close to the separatrices associated with unsta-ble equilibria or resonances of the secular angles (unstable periodicorbits).Figure 10 shows the temporal evolution of the mean orbitswith initial conditions close to the origin in Fig. 4 (see its captionfor the details). The eccentricity of the inner orbit is significantlyperturbed while the outer orbit is almost unchanged. The mutualinclination evolves in anti-phase with e , and ω librates around π /
2. These are typical features of the Lidov-Kozai resonance, stillwe recall that in this instance, the smaller outer body forces the LKcycles on the orbits of much more (20 times) massive than the innercomponent of the binary.
Figure 10.
The mean orbits with the following initial condition: m = (cid:12) , m =
200 m J , m =
10 m J , a = . a =
30 au, e = . e = . ω = π / ω = π / i mut = ◦ . A = .
12. Panels from thetop to the bottom present the time-evolution of e , e , i mut , ω , ω , respec-tively. The dynamical maps and their analysis show that the equilibriaconstrain the secular dynamics of the perturbed model. In termsof the Newtonian, point-mass formulation, the stationary solutionsdepend on parameters α and µ . Limiting our survey to i mut < π / µ , covering a transition from the planetary regime(small µ ) to the circumbinary case (large µ ) than in (Migaszewski& Go´zdziewski 2009a). Yet the assumption of small α makes itpossible to use the analytic formulation of the secular Hamiltonian,which is very precise in terms of the octupole-level approximation,instead of the numerical approach in (Migaszewski & Go´zdziewski2009a).The parameter dependence of the equilibria is illustrated inFig. 11 that shows the P S -plane (note that due to the symmetry, onlythe upper half-plane is illustrated, see also Fig. 3). We set α = . µ ∈ [ . , , , , , , ] , and then positions of the equilibria maybe traced along increasing A , which might be understood as thenatural curve parameter.In fact, instead of µ , we choose a new parameter β (see Krasin-skii 1974; Migaszewski & Go´zdziewski 2010) β ( µ , α ) ≡ L / L ∼ µ √ α , that better reflects the dependence of the dynamics on both µ and α than on one of these parameters itself (we did a few numerical testsfor different α which confirm this scaling very well). Actions L i c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
Figure 7.
The left panel is for the levels of the secular energy (thin curves) depicted in the representative plane of the quadrupole problem, S ” ≡ { e cos2 ω × e } -plane, α = . µ = A = .
12. Recall that this problem is integrable and e ≡ constant is a parameter. The black, thick curve marks the LKR equilibriumfor different energies and e integral. The red, dotted curve marks the separatrix between librations and rotations of ω for fixed values of integrals. The shaded(green) zone indicates librations of ω around π / The right hand panel is for the phase diagram in the { e cos2 ω × e sin2 ω } -plane and a fixed energylevel marked by blue curve in panel a. The separatrix and the region of librations of ω are also marked with green colour. Figure 8.
The left panel is for the secular energy levels (thin curves) shown in the S (cid:48) -plane of the octupole model, α = . µ = A = .
12. Shaded areasare for initial conditions leading to chaotic evolution of the secular model. The red curve is for the separatrix of the LKR resonance in the quadrupole model,the black thick curve is for the equilibria in this model.
The right hand panel is for the Poincar´e cross-section ∆ϖ = π in the { e cos2 ω × e sin2 ω } -plane.It corresponds to a fixed energy level marked with the blue curve in panel a. The separatrix of the quadrupole problem and a region of librating ω are alsomarked. Figure 9.
The spectrum of fundamental frequencies computed along the x -axis of Fig. 8, for the same set of parameters. Abscissas, in which the frequenciesapproach zero, indicate the separatrices of secular resonances of the respective angles, ω (panel a), and ∆ϖ (panel b). Parameters are: m = m (cid:12) , m = J , m = J , a = . a =
10 au. c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets reflect the angular momentum partition between both components.If β ∼ ( µ ∼ ) , both secondaries are dynamically equivalent, if β > µ >
5) then the inner body “dominates” dynamically in thesystem, and if β < µ and β parameters.In the right-hand half-plane (quarter IV, ω = ω = π / i mut < π /
2, only one solution appears that is the LK resonance. Forsmall µ , it evolves with increasing A along the axis of e ∼ e = e ∼
1. After the second bifurcation of the LKR (Fig. 2e),two solutions emerge. One of them is the LKR for i mut > π / e ∼ e .For µ (cid:62)
2, the LKR does not reach e = e and increasing e . For larger mass ratio, themaximal e is smaller. The families of LKR for µ = . , β <
1, hence the dynamical hierarchy isreversed, and the eccentricity of the inner, more massive body isforced by the outer companion. The position of the LKR familymoves to the region of smaller e , and for large enough µ , this solu-tion is confined to e ∼ e . The direction ofparametric evolution of this equilibrium, corresponding to increas-ing A , is marked with an arrow. We see that for all µ this directionis the same.The view of the left-hand half-plane is more complex. As weshown in (Migaszewski & Go´zdziewski 2009a), the first solutionappearing in quarter III of the P S -plane is unstable equilibrium IIIaemerging due to the second bifurcation of the origin ( e = e = ) (marked with I − in our paper). Then, with increasing A , this so-lution moves along e ∼ e . For some value of A (e.g., ∼ .
195 for α = . µ =
20; see Fig. 3) solutions IIIb andIVb+ appear in the same point of the phase space (this is illustratedwith crossed circle in the left-hand half-plane of Fig. 3). SolutionIVb+ then moves along e → e →
0. Simultaneously, solu-tion IIIb evolves towards a point marked with dotted circle. Solu-tion IIIa also moves to that point and for some critical A both solu-tions merge and vanish. Figure 11 reveals that the parametric pathsof equilibria depend on parameter µ ( α = .
04 was fixed in thistest). The path of the LKR becomes closer to e ∼ µ (cid:62)
25. Wemay also notice that the bifurcation points (crossed and dotted cir-cles, respectively) tend to each other with increasing mass ratio,which is consistent with the transition between the planetary andthe circumbinary regime. For µ (cid:62)
25, these two points are merged.In this particular case, the structure of equilibria in quarter III ofthe P S -plane is more simple than in the general case because onlyone solution persists, a stable equilibrium corresponding to nearlycircular outer orbit.For different α < .
1, the general, global view of the fami-lies of equilibria is very similar to the results presented in Fig. 11.In fact, as we noticed previously, the parametric evolution of theequilibria is reflected by parameter β ( µ , α ) , and basically does notdepend on the individual values of α and µ . The results presented in the previous section illustrate the alreadywell known feature of the classic model (see, e.g., Michtchenko& Malhotra 2004; Michtchenko et al. 2006). In the approxima-tion of small masses, the secular dynamics of the Newtonian 2-planet, hierarchical model depend on the semi-major axes ratioand planetary mass ratio, and not on individual system parame-ters, ( a , a , m , m ) . In our works (Migaszewski & Go´zdziewski2009b, 2010), we shown that this feature is not preserved in themore general model, including relativistic, rotational and tidal cor-rections to the Newtonian point-mass interactions. In these settings,the secular dynamics depend on the individual semi-major axes, aswell as individual planetary masses. Because the overall structureof the phase space is characterized by the equilibria, we now at-tempt to show that deviations between these equilibria in the classicand relativistic models become more important when the system di-mension scales down ( a , a decrease when α is const), and masses m , m are smaller, when their ratio µ is kept constant.The differences between the two coplanar models manifestitself through the shapes and localization of stationary solutionsand depend on the individual masses and semi-major axes (Mi-gaszewski & Go´zdziewski 2009b). Now we can observe the samefeature in the spatial planetary system, corrected for the relativisticperturbation. Figure 12 presents families of equilibria in the samemanner as Fig. 11 (due to the symmetry, only the y -positive half-plane of P S -plane is presented). Families of stationary solutions inthe classic model are drawn with blue and violet curves for sta-ble and unstable solutions, respectively. Solutions in the relativisticmodel are plotted with black (stable equilibria) and grey (unsta-ble equilibria) curves. In this experiment, we varied the individualmasses of the secondary bodies, still keeping their ratio µ =
10 as aconstant. The masses are changed between ( m =
100 m J , m =
10 m J ) to ( m = J , m = . J ) , and the primary mass is m = (cid:12) .The results illustrated in Fig. 12 confirm that even for largesecondary masses, the parametric curves of stationary solutions inthe realm of the classic model depends weakly on the individualmasses. Yet similarly to the co-planar case, curves of equilibriacalculated with the relativistic corrections differ significantly fromthe results obtained for the classic model. The deviations betweenboth models are more significant for smaller masses of the sec-ondary bodies. For instance, the family IVa moves to the range ofmuch smaller e . Solutions of this type exist even for very small m and m , which is just not possible in the Newtonian model (Mi-gaszewski & Go´zdziewski 2009a). When the mutual perturbationsbetween secondaries are small enough, the critical inclination in therelativistic model, that leads to the LK bifurcation, becomes largerthan π / i crit increases with decreasing masses.In quadrant III of the P S -plane ( ω = − ω = ± π / m ∼ .
07 m J ) the parametric paths divide into two parts.The top part, characterized by larger e , comprises of two unstableequilibria emerging from one bifurcation point, which then meet inanother bifurcation point. The bottom part of the equilibria curverepresents a saddle point, that changes its stability, from unstableIIIa-type solution to the stable solution IVb+. This branch is verysimilar to equilibria in the classic system with large µ (or large β >
5) observed already in Fig. 11, however this takes place forquite a different value of β = P S -plane, Fig. 13. c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
Figure 11.
Families of stationary solutions in the P S -plane. The semi-major axes ratio α = . µ = { . , , , , , , } (each curve is labeled accordinglywith the value of µ ). Dark dots are for stable equilibria, grey dots are for unstable equilibria. Crossed and dotted circles mark the positions on the equilibriacurves where solutions bifurcate (see the text for details). Small arrows show the directions of the evolution of particular equilibrium with increasing A .Parameter β ≡ L / L (see the text). Stationary curves corresponding to β = Figure 12.
Families of stationary solutions presented in the P S -plane, calculated for α = . , µ = , a = . , a = . , m = (cid:12) and varied m = , , . , . , . , . , . , . J . Equilibria in the classic model are compared with the equilibria in the relativistic model. The mass of the outerbody m labels each particular curve. Black dots are for stable equilibria, grey dots are for unstable equilibria of the relativistic model. Equilibria of the classicmodel are plotted with blue and violet dots for stable and unstable solutions, respectively. Positions of equilibria were calculated with the help of the octupoletheory. Each panel in this figure is calculated for the same parameters µ , α , A (but in this case, particular values of masses and semi-majoraxes are varied). Let us recall, that Fig. 13a shows the phase spacecalculated in classic model, while subsequent panels 13b,c,d arefor relativistic model, and different masses and semi-major axes,labeled in subsequent panels.If the masses are large (Fig. 13b), we can see four elliptic points separated by four saddles (let us recall that the P S -plane isredundant, and the energy levels are reflected with respect to theorigin, thus in fact we have only two unique elliptic points and onlytwo saddles). The elliptic points may be identified with solutionsIVa ( ω = ω = ± π /
2) and IIIb ( ω = − ω = ± π / e ∼
0) and IVb+ ( e ∼ µ is kept constant), the structure surround- c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets ing solution IIIb becomes smaller and moves towards solution IIIa.Simultaneously, the saddle point IVb+ tends to the origin. For themasses small enough, (Fig. 12), solutions IIIa and IIIb merge andvanish, while IVb+ “falls” into the origin.Figures 14–16 shows the dynamical maps for the relativis-tic model, and are constructed in the same manner as maps inFigs. 4–6. Subsequent figures correspond to masses and A usedto plot the energy levels in Figs. 13a,c,d respectively: Fig. 14 isfor classic model, Figs. 15 and 16 are for the relativistic modelwith m =
10 m J and m = .
85 m J , respectively. The mass ra-tio is µ =
10 in both instances. The order of panels and symbol-coding of equilibria, as well as coding libration zones of the angles ω , ω , and ∆ϖ are the same as in Figs. 4–6; in particular, dot-ted, crossed and empty circles mark the Lyapunov stable, unstableand linearly stable equilibria, respectively. Clearly, the overall viewof the phase space is different in all cases. The regions of chaoticmotions (yellow areas in the σ -maps) obtained for the classic andrelativistic model are significantly different. Also the dependenceon the masses of secondaries in the relativistic models is evident.The structure of chaotic/regular secular evolution is reflected inmax e , -maps and through librational regions of ω and ω . Wenote, that in the case of the regular solutions, ω librates around ± π / i mut > i crit . In some part of this region also ω libratesaround ± π / e = .
45 and e = . e is strongly amplified, com-pared to the variations in the relativistic model. Still, this is not arule. Inspecting the bottom row of Fig. 16, we can find regions inthe P S -plane, in which, for the same initial condition, max e be-comes larger in the relativistic model than in the Newtonian model.Also the secularly chaotic configurations appear in quite differentzones of the phase space in both models. Moreover, the relativis-tic corrections may transform the regular evolution in the classicmodel into the chaotic evolution in the relativistic systems, and viceversa . Remarkably, the configuration illustrated in Fig. 14 has verylarge masses m =
100 m J and m =
10 m J . In this work, we attempt to show that the global features of the sec-ular dynamics of the 3-D, non-resonant planetary system dependqualitatively on the apparently subtle relativistic corrections to theNewtonian gravity. A lesson, which we learned studying the co-planar case (Migaszewski & Go´zdziewski 2009b) is that the non-Newtonian point-mass interactions might be very important for theglobal dynamics, because, in contrary to intuition one may have,the corrections to the Newtonian interactions might be not small ,as compared to the mutual point-to-point gravity. The numericalanalysis of this multi-parameter problem is complex, hence the un-derlining idea of this paper lies in a construction of possibly pre-cise analytical model. Essentially simple averaging of the perturb-ing Hamiltonian in (Migaszewski & Go´zdziewski 2008), makes itpossible to derive the analytic secular theory up to the order of 10in the semi-major axes ratio α . The accuracy of this model may becompared with the results of the numerical averaging of the pertur- Figure 17.
Evolution of the mean orbits in the following system: m = (cid:12) , m =
100 m J , m =
10 m J , a = . a = e = . e = . ω = π / ω = π / i mut = ◦ . A = . e , e , i mut , ω , ω , respectively.The grey curves are for the classic secular model, the black curves are forthe model including the relativistic corrections to the potential of the innerbinary. bation (Michtchenko & Malhotra 2004). Moreover, for a class of hi-erarchical systems considered in our paper, already the third-ordermodel is precise enough to find and investigate the qualitative fea-tures of the system. In this work, we focus on three-body configura-tions, e.g., a star and two massive planets or a binary stars and oneplanet. “Fortunately”, the second-order model is integrable, hencethe octupole-level approximation might be considered as the firstorder perturbation to this analytically soluble case. This makes itpossible to understand the sources of instabilities appearing in thefull (non-averaged) model.The averaging over mean anomalies reduces the dynamics toa system having two DOF, which then may be investigated withthe help of rich geometrical tools, like the Poincar´e cross sectionand the representative planes of initial conditions introduced inMichtchenko et al. (2006). Unfortunately, all additional correctionsthat increase the DOF number must be neglected here, as for in-stance the rotational quadrupole moment of the star and/or of theplanets. This is the price that must be paid for the possibly globalmodel of the dynamics. In the assumed range of a and masses,the relativistic “corrections” in fact compare with the Newtonianpoint-mass mutual interactions, and are much larger than other per-turbations, like the tidal and rotational distortions of the bodies.The analysis of the secular frequencies introduced by various cor-rections justify that the model only takes into account the PN per-turbation to the potential of star and the inner secondary.Having the two DOF model, we investigate the simplest class c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
Figure 13.
Energy levels presented in the P S -plane, and calculated for the following parameters α = . , µ = , A = . , i = ◦ , a = . , a = . , m = m (cid:12) . A sequence of panels demonstrates the view of the phase space when masses of the secondary bodies decrease (their ratio is constant).From panel b) to panel d) ( m , m )[ m J ] are ( , ) , ( . , . ) , ( , ) . Shaded areas mark ranges of i mut larger than 60 ◦ , ◦ , (the light– and dark– grey,respectively). Figure 14.
Dynamical maps for the classic, Newtonian model in the P S -plane. The model parameters are α = . , µ = , A = . , i = ◦ , a = . , a = . , m = m (cid:12) , m =
100 m J , m =
10 m J . See the text and caption to Fig. 4 for the details. c (cid:13) , 1–18 he relativistic dynamics of circumbinary planets Figure 15.
The dynamical maps for the octupole model with relativistic corrections. Masses of the secondary bodies are m = J , m = J . See thetext and caption to Fig. 4 for more details. Figure 16.
The dynamical maps for the octupole model with relativistic corrections. Masses of the secondary bodies are m = . J , m = . J . See thetext and caption to Fig. 4 for more details. of solutions that are the equilibria. We focus on the Lidov-Kozairesonance (LKR) in systems characterized by large range of themass ratio µ . This part extends the results derived for small µ in(Migaszewski & Go´zdziewski 2009a). We found that even muchsmaller outer body (planetary mass m with respect to sub-stellarvalues of m ) moving in a wide, highly inclined ( i mut > ∼ ◦ ) or-bit may significantly perturb the inner orbit. In turn, the restrictedmodel of the circumbinary planet, when we assume that the planet does not influence the binary, is not generally valid, even if the in-ner mass m is 100 times larger than the outer body m .We also studied the parametric structure of families of partic-ular equilibria classified as IVa, IIIa, IIIb, IVb+ in (Migaszewski& Go´zdziewski 2009a) for small α . Thanks to this assumption, theanalytic model makes it possible to investigate the transition be-tween the planetary regime (small µ ) and the circumbinary regime( µ ∼ , c (cid:13) , 1–18 C. Migaszewski and K. Go´zdziewski
A particularly interesting feature of the ocupole model is theappearance of the secular chaos . We found that if i mut exceeds thecritical inclination i crit ( µ , α ) , the long term evolution of the systemmay be strongly chaotic, leading to large amplification of the eccen-tricities. In the regular regions of the phase space, the mean angle ω librates around ± π /
2. In some parts of these regions also thesecond secular angle ω librates around ± π /
2. The initial condi-tions satisfying i mut > i crit lead to strong amplification of the innereccentricity e . Simultaneously, for the same values of the angu-lar momentum of the system, we may observe strong amplifica-tion of e in some region of the phase space, with almost constant relative inclination of the orbits. This behaviour may be attributedto unstable equilibrium IIIb emerging in the secular system. Theamplification of e happens not only for librations of ω around ± π / [ , π ] . Thedynamical maps reveal that the primarily source of the chaotic mo-tions are the unstable equilibria and unstable periodic orbits in thefull system, following the appearance of separatrices in the inte-grable, quadrupole-level model.Thanks to the simple analytic model, the influence of relativis-tic correction H GR on the global secular dynamics of the problemmight be clearly demonstrated. A simple proof of this influence isprovided by the analysis of the equilibria in the perturbed model.The differences between the Newtonian and relativistic models arelarger when the mutual interactions between the secondaries areweaker, e.g., when companion masses are smaller. Yet the dynam-ics are basically very simple up to the limit of the critical inclina-tion, when the first bifurcation of the origin ( e = e =
0) occurs,and this feature of the dynamics is preserved in both models.We stress that although the analysis is done for specific, dis-crete mass ratios, the results are valid as far the assumptions ofthe averaging theorem are fulfilled and the corrections besides gen-eral relativity are negligible. We demonstrated that similarly to theco-planar problem, the global 3-D dynamics of the classic, Newto-nian model essentially depend only on the ratios of semi-major axesand masses of the secondaries. Hence, although we consider mostlythe circumbinary configurations, the results may be are easily ex-trapolated to the “typical” planetary regime, investigated already(Michtchenko et al. 2006). Moreover, when the PN corrections areadded to the model, the dynamics are much more complex. Stillthe global picture of the phase space is determined by the ratio ofthese corrections and the Newtonian mutual interactions. Hence, ifthe system scales down, while this ratio is roughly preserved, thestructure of the phase space, determined by stationary solutions ofthe secular system should not essentially change.Our approach may be generalized for other perturbations, likethe rotational and conservative tidal distortion of the bodies in thesystem. Unfortunately, in the most general case, the dimension ofthe hierarchical system cannot be generally reduced to two DOF.Moreover, these perturbations lead to even more interesting andintriguing dynamics, which we investigated in the co-planar andspatial case (e.g., Fabrycky & Tremaine 2007; Mardling 2007;Ragozzine & Wolf 2009; Migaszewski & Go´zdziewski 2009b). Wework on a global approach, suitable for the 3D systems, aiming topublish these results in future papers.
ACKNOWLEDGEMENTS
We thank the anonymous referee for comments that improved themanuscript. This work is supported by the Polish Ministry of Sci- ence and Higher Education, through grants N/N203/402739 and92/N-ASTROSIM/2008/0. CM is a recipient of the stipend of theFoundation for Polish Science (programme START, edition 2010).
REFERENCES
Adams F. C., Laughlin G., 2006, International Journal of ModernPhysics D, 15, 2133Arnold V. I., Kozlov V. V., Neishtadt A. I., 1993, Dynamical sys-tems III. Mathematical aspects of classical and celestial mechan-ics. Encyclopaedia of mathematical sciences, Springer VerlagBatygin K., Bodenheimer P., Laughlin G., 2009, ApJL, 704, L49Brouwer D., Clemence G. M., 1961, Methods of celestial mechan-ics. New York: Academic Press, 1961Brumberg V., 2007, Celestial Mechanics and Dynamical Astron-omy, 99, 245Eggenberger A., 2010, in K. Go´zdziewski, A. Niedzielski, &J. Schneider ed., EAS Publications Series Vol. 42. pp 19–37Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298Farago F., Laskar J., 2010, MNRAS, 401, 1189Ferraz-Mello S., ed. 2007, Canonical Perturbation Theories - De-generate Systems and Resonance Vol. 345 of Astrophysics andSpace Science LibraryFerrer S., Os´acar C., 1994, Celestial Mechanics and DynamicalAstronomy, 60, 187Ford E. B., Kozinsky B., Rasio F. A., 2000, ApJ, 535, 385Gronchi G. F., Milani A., 1998, Celestial Mechanics and Dynam-ical Astronomy, 71, 109Harrington R. S., 1968, AJ, 73, 190Kozai Y., 1962, AJ, 67, 579Krasinskii G. A., 1972, Celestial Mechanics, 6, 60Krasinskii G. A., 1974 Vol. 62 of IAU Symposium. pp 95–116Laskar J., 1990, Icarus, 88, 266Laskar J., 2000, Physical Review Letters, 84, 3240Lee M. H., Peale S. J., 2003, ApJ, 592, 1201Libert A.-S., Henrard J., 2007, Icarus, 191, 469Lidov M. L., 1962, Planetary and Space Science, 9, 719Lidov M. L., Ziglin S. L., 1976, Celestial Mechanics, 13, 471Mardling R. A., 2007, MNRAS, 382, 1768Mardling R. A., 2010, MNRAS, 407, 1048Michtchenko T. A., Ferraz-Mello S., Beaug´e C., 2006, Icarus,181, 555Michtchenko T. A., Malhotra R., 2004, Icarus, 168, 237Migaszewski C., Go´zdziewski K., 2008, MNRAS, 388, 789Migaszewski C., Go´zdziewski K., 2009a, MNRAS, 395, 1777Migaszewski C., Go´zdziewski K., 2009b, MNRAS, 392, 2Migaszewski C., Go´zdziewski K., 2010, in K. Go˙zdziewski,A. Niedzielski, & J. Schneider ed., EAS Publications SeriesVol. 42, . pp 385–391Murray C. D., Dermott S. F., 2000, Solar System Dynamics. Cam-bridge Univ. PressRabl G., Dvorak R., 1988, A&A, 191, 385Ragozzine D., Wolf A. S., 2009, ApJ, 698, 1778Richardson D. L., Kelly T. J., 1988, Celestial Mechanics, 43, 193Takeda G., Kita R., Rasio F. A., 2009, in IAU Symposium Vol. 253of IAU Symposium, . pp 181–187Tamuz, O., et al., 2008, A&A, 480, L33ˇSidlichovsk´y M., Nesvorn´y D., 1996, Celestial Mechanics andDynamical Astronomy, 65, 137 c (cid:13)000