Eccentricity generation in hierarchical triple systems with non-coplanar and initially circular orbits
aa r X i v : . [ a s t r o - ph . E P ] A ug Eccentricity generation in hierarchical triple systemswith non-coplanar and initially circular orbits
Nikolaos Georgakarakos ( [email protected] ) School of Mathematics, Edinburgh University, Mayfield Road, Edinburgh EH9 3JZ,UK
Abstract.
In a previous paper, we developed a technique for estimating the innereccentricity in coplanar hierarchical triple systems on initially circular orbits, withcomparable masses and with well separated components, based on an expansion ofthe rate of change of the Runge-Lenz vector. Now, the same technique is extended tonon-coplanar orbits. However, it can only be applied to systems with I < . ◦ or I > . ◦ , where I is the inclination of the two orbits, because of complicationsarising from the so-called ’Kozai effect’. The theoretical model is tested againstresults from numerical integrations of the full equations of motion. Keywords:
Celestial mechanics, stellar dynamics, binaries:general
1. Introduction
In a previous paper (Georgakarakos, 2002), we derived a formula forestimating the eccentricity of the inner binary of a hierarchical triplesystem with well separated components and initially circular orbits.However, the derivation was based on the assumption that the orbitsof the three bodies were on the same plane. In the current paper, weextend the derivation to the more general case of non-coplanar orbits.Generally, the method is the same as in the coplanar regime, althoughmore terms are included in the derivation of the short period equa-tions to improve the accuracy of the model. There are also some minorchanges in the calculation of the secular contribution to the eccentricity.Other recent work on the dynamics of hierarchical triple systemincludes the work done by Ford et al. (2000), Kiseleva et al. (1998) andKrymolowski and Mazeh (1999).
2. Theory
The theoretical model is constructed in the same way as in the copla-nar case: we derive expressions for the short period terms by usingthe definition of the Runge-Lenz vector, while the secular evolution isstudied by means of canonical perturbation theory. The combinationof the short period and secular part of the eccentricity is achieved by c (cid:13) Kluwer Academic Publishers. Printed in the Netherlands. pap2cmdarev.tex; 3/07/2018; 0:04; p.1
Nikolaos Georgakarakos considering the eccentricity (inner or outer) to consist of a short periodand a long period (secular) component, i.e. e = e short + e sec (one canpicture this by recalling the expansion of the disturbing function in solarsystem dynamics, where the perturbing potential is given as a sum ofan infinite number of cosines of various frequencies). Thus, consideringthe eccentricity to be initially zero leads to e short = − e sec (initially),which implies that, although the eccentricity is initially zero, the shortperiod and secular eccentricity may not be.2.1. Calculation of the short-period contribution to theeccentricity
The equation of motion of the inner binary, using the Jacobi notation( r for the relative position vector of the inner binary and R the vectorfrom the centre of mass of m and m to the outer mass m ), is:¨ r = − G ( m + m ) r r + F , (1)where F , the perturbation to the inner binary motion, is F = Gm ( R − µ r | R − µ r | − R + µ r | R + µ r | ) = Gm ∂∂ r ( 1 µ | R − µ r | ++ 1 µ | R + µ r | ) (2)with µ i = m i m + m , i = 1 , . Now, since we are dealing with hierarchical triple systems with wellseparated components, implying that r/R is small, the inverse distancesin equation (2) can be expressed as:1 | R − µ r | = 1 R ∞ X n = o (cid:16) µ rR (cid:17) n P n (cos θ )and 1 | R + µ r | = 1 R ∞ X n = o (cid:16) − µ rR (cid:17) n P n (cos θ ) , where P n are the Legendre polynomials and θ is the angle between thevectors r and R . Expanding to third order, the perturbation becomes F = Gm ∂∂ r (cid:18)
32 ( r · R ) R − r R − µ − µ )2 ( r · R ) R ++ 3( µ − µ )2 r ( r · R ) R (cid:19) . (3) pap2cmdarev.tex; 3/07/2018; 0:04; p.2 ccentricity generation in HTS P ), while the other two come from the octupole term ( P ).Using now the definition of the Runge-Lenz vector, we can obtainan expression for the inner eccentricity. The inner eccentric vector e is given by e = − r r + 1 µ ( ˙ r × h ) , (4)where h = r × ˙ r and µ = G ( m + m ) . Assuming that r · ˙ r = 0, i.e. theinner binary remains nearly circular (bear in mind that all bodies are oninitially circular orbits), differentiating equation (4) and substitutingfor F we obtain:˙ e = Gm µR (cid:20)(cid:18) r · R )( ˙ r · R ) R − µ − µ ) ( r · R ) ( ˙ r · R ) R ++ 3( µ − µ ) r ( ˙ r · R ) R (cid:19) r + (cid:18) r − r · R ) R + 152 ( µ −− µ ) ( r · R ) R −
92 ( µ − µ ) r ( r · R ) R (cid:19) ˙ r (cid:21) . (5)Now, choosing a frame of reference such that the initial plane of theinner orbit is our reference plane and the line of nodes is initially onthe x-axis, with the positive direction of the x-axis pointing at theascending node of the outer orbit, the Jacobi vectors can be representedapproximately in polar form as r = a (cos n t, sin n t,
0) and R = a (cos ( n t + φ ) , sin ( n t + φ ) cos I, sin ( n t + φ ) sin I ) , where a and a are the semi-major axes of the inner and outer orbit respectively and φ is the initial relative phase of the two binaries. After integrating, thecomponents x and y of the eccentric vector become (expanding thecoefficients in powers of X and retaining the four leading terms, with X being the period ratio of the two orbits): x = m M X (cid:20) P x21 ( t ) + 1 X P x22 ( t ) + m ∗ (cid:18) X P x31 ( t ) + 1 X P x32 ( t ) (cid:19)(cid:21) ++ C x (6) y = m M X (cid:20) P y21 ( t ) + 1 X P y22 ( t ) + m ∗ (cid:18) X P y31 ( t ) + 1 X P y32 ( t ) (cid:19)(cid:21) ++ C y (7)where P i ( t ) are given in the appendix and m ∗ = m − m ( m + m ) M . (8) pap2cmdarev.tex; 3/07/2018; 0:04; p.3 Nikolaos Georgakarakos M is the total mass of the system and C x and C y are constants ofintegration. The semi-major axes and mean motions were treated asconstants in the above calculation.2.2. Calculation of the secular contribution to theeccentricity
Secular terms cannot be obtained by the method of Sect. 2.1, because,for an eccentric outer binary, those terms appear as a linear functionof time in the expansion of the eccentric vector and therefore, they arevalid for limited time. Thus, the secular contribution to the eccentricityis studied by means of the Von Zeipel method.The doubly averaged over short period terms for both orbits Hamil-tonian is (Marchal, 1990, with changes of notation): H = − Gm m a S − G ( m + m ) m a T + Q + Q , (9)where Q = 18 Gm m m a ( m + m ) a (1 − e ) [ − − e + 3 sin I (1 − e ++5 e sin g S )] , (10) Q = 15 Gm m m ( m − m ) a e S e T m + m ) a (1 − e ) [(sin g S sin g T cos I ++ cos g S cos g T )(4 + 3 e − I (1 − e + 7 e sin g S )) −− − e ) sin I cos I sin g S sin g T ] . (11)The subscripts S and T denote the inner and outer longer period orbitrespectively. The first term in the Hamiltonian is the Keplerian energyof the inner binary, the second term is the Keplerian energy of the outerbinary, while the other two terms represent the interaction between thetwo binaries. The Q term comes from the P Legendre polynomial andthe Q term comes from the P Legendre polynomial. There are alsoterms which arise from the canonical transformation, but they are ofsmaller order than the P term.By using Hamilton’s equations (see Marchal, 1990 for the derivationof equations of motion involving the P term), we can now derive theaveraged equations of motion of the system (see appendix), which areconsistent with the ones in Ford et al. (2000), except that the P termin the Hamiltonian in Ford et al. (2000) has the wrong sign. The samesign error appears in Krymolowski and Mazeh (1999). After some ex-ploratory numerical integrations of the five equations of motion of the pap2cmdarev.tex; 3/07/2018; 0:04; p.4 ccentricity generation in HTS x S , y S and e T abovethe first order) and keeping the dominant term in the equation for ˙ g T produced the following simpler system of differential equations:d g T d τ = A d x S d τ = − By S + Ce T sin g T (12)d y S d τ = Dx S − Ee T cos g T where x S = e S cos g S , y S = e S sin g S ,A = cos I + 12 β (4 − I ) , B = 2 − I + β cos I,C = 516 α cos I (4 −
15 sin I ) , D = 2 + β cos I, E = 516 α (4 − I ) ,α = m − m m + m a S a T , β = m m M m ( m + m ) ( a S a T ) , d τ = 34 G m a S a ( m + m ) d t. The solution to system (12) is: g T ( τ ) = Aτ + g T0 (13) x S ( τ ) = K cos √ BDτ + K sin √ BDτ ++ AC + BEBD − A e T cos ( Aτ + g T0 ) (14) y S ( τ ) = K r DB sin √ BDτ − K r DB cos √ BDτ ++ AE + CDBD − A e T sin ( Aτ + g T0 ) , (15)where K , K are constants of integration and g T0 is the initial valueof g T . It should be pointed out here that although the initial outereccentricity is zero, the secular one is not and consequently an initialouter secular argument of pericentre can be defined. However, the aboveapproximation to secular motion is not valid if 39 . ◦ < I < . ◦ ,since the inner eccentricity is expected to become significant due to pap2cmdarev.tex; 3/07/2018; 0:04; p.5 Nikolaos Georgakarakos the Kozai effect (Kozai, 1962) and therefore the assumption that theeccentricity remains small is invalid in this case.2.2.1.
Calculation of the initial outer secular eccentricity
As was stated in the previous section, the outer secular eccentricityremains almost constant. Thus, the only thing that remains now is toobtain an estimate for the initial outer secular eccentricity. This canbe done in the following way: First, we find an expression for the shortperiod outer eccentricity, by following the same procedure as we did inSect. (2.1), but this time we do it for the outer orbit. The equation ofmotion of the outer binary is¨ R = − GM (cid:18) µ R + µ r | R + µ r | + µ R − µ r | R − µ r | (cid:19) (16)and eventually we obtain, to leading order, for the components of theouter short-period eccentric vector: x = M ∗ X [cos I ( 1516 cos ( n t + φ ) −
716 cos (3 n t + 3 φ )) ++ 716 cos (3 n t + 3 φ ) −
316 cos ( n t + φ )] + C x (17) y = M ∗ X [cos I ( 2116 sin ( n t + φ ) −
716 sin (3 n t + 3 φ )) ++ cos I ( 716 sin (3 n t + 3 φ ) −
916 sin ( n t + φ ))] + C y (18) z = M ∗ X [cos I sin I ( 2116 sin ( n t + φ ) −
716 sin (3 n t + 3 φ )) ++ sin I ( 716 sin (3 n t + 3 φ ) −
916 sin ( n t + φ ))] + C z , (19)where M ∗ = m m ( m + m ) M . (20)Suppose now that the outer secular eccentric vector is e T = ( x T , y T , z T ).Then, the constants C x , C y and C z in equations (17), (18) and(19) can be replaced by x T , y T and z T , since δx T δx << y T − y and z T − z . Considering that the outer binary is initiallycircular, i.e. e out = 0, the secular outer eccentric vector will initially be: e T = ( x T , y T , z T )where x T = − M ∗ X [cos I ( 1516 cos φ −
716 cos 3 φ ) + 716 cos 3 φ − pap2cmdarev.tex; 3/07/2018; 0:04; p.6 ccentricity generation in HTS −
316 cos φ ] (21) y T = − M ∗ X [cos I ( 2116 sin φ −
716 sin 3 φ ) + cos I ( 716 sin 3 φ −
916 sin φ )] (22) z T = − M ∗ X [cos I sin I ( 2116 sin φ −
716 sin 3 φ ) + sin I ( 716 sin 3 φ −−
916 sin φ )] . (23)(24)Thus, the initial outer secular eccentricity is: e T = q x + y + z . (25)2.3. A formula for the inner eccentricity
The expressions that were derived in paragraphs (2.1) and (2.2) for theshort period and secular contribution to the inner eccentric vector, canbe combined to give an expression for the total eccentricity in the sameway we obtained the estimate for the initial outer secular eccentricity,i.e. by replacing the constants in equations (6) and (7) by equations (14)and (15), since the latter evolve on a much larger timescale. This yields(the minus sign in front of x S and y S has to do with the orientation ofthe coordinate system used in Marchal, 1990): x in = x − C x − x S (26) y in = y − C y − y S (27)The constants K and K in equations (14) and (15) are determinedby the fact that the inner eccentricity is initially zero.Averaging over time and over the initial relative phase φ , the aver-aged square inner eccentricity will be given by: e = < x + y > = m M X (cid:20) − I + 36964 cos I ++ DB ( 1916 −
52 cos I + 3 cos I ) + BD cos I + 1 X [ 1011144 ++ 47124 cos I + 543144 cos I + 499 DB cos I + BD ( 12136 ++ 119 cos I + 19 cos I )] + 1 X [ 453 cos I + 936 cos I + pap2cmdarev.tex; 3/07/2018; 0:04; p.7 Nikolaos Georgakarakos + DB ( 76 cos I + 143 cos I ) + BD ( 113 cos I + 23 cos I )] ++ m ∗ X [ 12758192 + 365258192 cos I − I + 768758192 cos I ++ DB ( 139258192 cos I − I + 256258192 cos I ) + BD ( 4258192 −− I + 31258192 cos I )] + m ∗ X [ 138519131072 + 62289131072 cos I ++ 121185131072 cos I + 102375131072 cos I + DB ( 54333131072 − I ++ 119025131072 cos I ) + BD ( 94113131072 cos I − I ++ 10125131072 cos I )] + m ∗ X [ 124958192 cos I − I ++ 243758192 cos I + DB ( 2554516384 cos I − I + 4312516384 cos I ) ++ BD ( − I − I + 562516384 cos I )] (cid:21) − m m ∗ M ∗ M X ×× [( 3351024 cos I − I + 17751024 cos I )(1 + DB ) AC + BEBD − A ++( 1551024 cos I − I + 8751024 cos I )(1 + BD ) AE + CDBD − A ] −− m m ∗ M ∗ M X [( 2194096 − I + 37954096 cos I )(1 + DB ) ×× AC + BEBD − A + ( 6874096 cos I − I + 15754096 cos I ) ×× (1 + BD ) AE + CDBD − A ] + M ∗ X [ ( AC + BE ) ( BD − A ) (1 + DB )( 29512 −− I + 137512 cos I ) + ( AE + CD ) ( BD − A ) (1 + BD ) ×× ( 65512 cos I − I + 245512 cos I )] + 12 M ∗ X ×× (cid:20) ( AC + BE ) ( BD − A ) + ( AE + CD ) ( BD − A ) (cid:21) ( 29256 − I −− I + 245256 cos I ) . (28)The above formula is expected to be rather inaccurate in situationswhere the system parameters yield very small values for the quantity BD − A , i.e. when we are near to a secular resonance. Although thequantities √ BD and A are approximations to the secular frequencies of pap2cmdarev.tex; 3/07/2018; 0:04; p.8 ccentricity generation in HTS BD − A = 0, which reduces to3(1 − I )(1 − β ) − β sin I = 0 , (29)could roughly identify the location of the secular resonance. However,in order to get a more accurate solution to the problem, the inclusion ofmore terms is necessary in system (12) or in the averaged Hamiltonian.
3. Comparison with numerical results
In order to test the range of applicability of the theory developed in theprevious sections, we integrated the full equations of motion numeri-cally, using a symplectic integrator with time transformation (Mikkola,1997).The code calculates the relative position and velocity vectors of thetwo binaries at every time step. Then, by using standard two bodyformulae, we computed the orbital elements of the two binaries. In-tegrations were arranged such that the writing index
Iwr was 1, theaverage number of steps per inner binary period
N S was 60, the methodcoefficients a a icor was 1. In all simulations, we confined ourselves to systemswith mass ratios within the range 10 : 1 since, “among stellar triples,mass ratios are rare outside a range of approximately 10 : 1, althoughsuch systems would be inherently difficult to recognise” (Eggleton andKiseleva, 1995); and initial period ratio X ≥
10. We also used unitssuch that G = 1 and m + m = 1 and we always started the integra-tions with a = 1. In that system of units, the initial conditions for thenumerical integrations were as follows: r = 1 , r = 0 , r = 0 R = a cos φ, R = a sin φ cos I, R = a sin φ sin I ˙ r = 0 , ˙ r = 1 , ˙ r = 0˙ R = − r Ma sin φ, ˙ R = r Ma cos φ cos I, ˙ R = r Ma cos φ sin I, where r and R are the relative position vectors of the inner and outerorbit respectively. pap2cmdarev.tex; 3/07/2018; 0:04; p.9 Nikolaos Georgakarakos
Short period effects
First we tested the validity of equations (6) and (7). The integrationsand comparison with the analytical results were done for φ = 90 ◦ , i.e.the outer binary was ahead of the inner one at right angles. However,this does not affect the qualitative understanding of the problem at all.These results are presented in Table 1, which gives the percent-age error between the averaged, over time, numerical and theoretical e in (the theoretical eccentricity was obtained by evaluating equations(6) and (7) everytime we had an output from the symplectic integra-tor; both averaged numerical and theoretical eccentricities were calcu-lated by using the trapezium rule). The integrations were performedover one outer orbital period time span (in our system of units, theinitial outer orbital period is T out = 2 πX , where X is the initialperiod ratio). For each pair ( m , X ) in Table 1, there are five en-tries, corresponding, from top to bottom, to the following inner bina-ries: m = 0 . − m = 0 . m = 0 . − m = 0 . m = 0 . − m = 0 . m = 0 . − m = 0 . m = 0 . − m = 0 .
5. A dash in Table 1 de-notes that the analogy among the masses was outside the range 10 : 1.The inclination of the two orbits is I = 20 ◦ .Most of the results show a rather significant error for systems withstrong perturbation to the inner binary (small X -large m ). However,the error drops considerably as we move to larger values of X (the errorbecomes less than 10% for all systems with X ≥ m = 0 . m = 4, X = 20, I = 20 ◦ , φ = 90 ◦ and the integration time span is oneouter orbital period ( T out = 125 . Short and long period effects
Next, we tested equation (28), which accounts for the short period andsecular effects to the inner eccentricity. The formula was comparedwith results obtained from integrating the full equations of motionnumerically. These results are presented in Table 2, which gives thepercentage error between the averaged, over time and initial phase φ ,numerical e and equation (28). In the case of a system with noticeablesecular evolution, the error is accompanied by the period of the longoscillation of the eccentricity, which is the same as the integration time pap2cmdarev.tex; 3/07/2018; 0:04; p.10 ccentricity generation in HTS Table I. Percentage error between the av-eraged numerical and averaged theoreti-cal e in . The theoretical model is basedon equations (6) and (7). For all sys-tems, I = 20 ◦ and φ = 90 ◦ . Each line cor-responds to a different inner binary pair(0 . − . , . − . , . − . , . − . . − . m \ X
10 15 20 25 30 500.05 - - - - - -- - - - - -- - - - - -- - - - - -5.5 2.6 1.5 1 0.7 0.30.09 5.1 2.5 1.6 1.1 0.9 0.45.6 2.7 1.7 1.2 0.9 0.45.9 2.8 1.7 1.2 0.9 0.46.1 2.9 1.8 1.2 0.9 0.46.2 3 1.8 1.2 0.9 0.40.5 10 5.4 3.6 2.7 2.1 1.210.4 5.5 3.6 2.6 2.1 1.110.7 5.6 3.6 2.6 2 111 5.7 3.7 2.7 2.1 111.2 5.9 3.8 2.8 2.2 1.11 13.9 7.5 5 3.7 3 1.714.2 7.6 5 3.7 2.9 1.614.5 7.7 5 3.7 2.8 1.514.7 7.9 5.1 3.7 2.9 1.515 8.1 5.3 3.9 3 1.61.5 - - - - - -16.7 8.9 5.9 4.3 3.4 1.816.9 9 5.9 4.3 3.4 1.817.3 9.3 6 4.4 3.4 1.817.7 9.5 6.3 4.6 3.6 1.92 - - - - - -18.5 9.9 6.5 4.8 3.7 218.7 10 6.5 4.8 3.7 1.919 10.2 6.7 4.9 3.8 219.4 10.4 6.9 5.1 4 2.12.6 - - - - - -- - - - - -20.2 10.8 7 5.2 4 2.120.6 11 7.2 5.3 4.1 2.121 11.3 7.5 5.5 4.3 2.23 - - - - - -- - - - - -21 11.2 7.3 5.4 4.2 2.221.4 11.5 7.5 5.4 4.3 2.221.8 11.7 7.7 5.7 4.4 2.33.4 - - - - - -- - - - - -- - - - - -22 11.8 7.7 5.6 4.4 2.322.4 12.1 8 5.8 4.6 2.44 - - - - - -- - - - - -- - - - - -22.9 12.2 8 5.8 4.6 2.323.2 12.5 8.3 6 4.7 2.54.5 - - - - - -- - - - - -- - - - - -- - - - - -23.8 12.8 8.4 6.2 4.8 2.55 - - - - - -- - - - - -- - - - - -- - - - - -24.2 13 8.6 6.3 4.9 2.6 pap2cmdarev.tex; 3/07/2018; 0:04; p.11 Nikolaos Georgakarakos
Figure 1.
Inner eccentricity against time for a system with m = 0 . m = 4, X = 20, I = 20 ◦ and φ = 90 ◦ . The integration time span is one outer orbital period( T out = 125 . T in = 2 π . span, while the rest of the systems were integrated over one outer orbitalperiod (denoted by a dash), since there was not any noticeable secularevolution when those systems were integrated over longer time spans.For each pair ( m , X ) in Table 2, there are three pairs of entries,corresponding, from top to bottom, to 10 ◦ , ◦ and 30 ◦ inclination. Anegative entry means that formula (28) is an overestimate to the numer-ical result. Each system was numerically integrated for φ = 0 ◦ − ◦ with a step of 45 ◦ . After the end of each simulation, e was averagedover time using the trapezium rule and after the integrations for all φ were done, we averaged over the relative initial phase by using therectangle rule. The integrations were also done for smaller steps in φ (10 ◦ , 1 ◦ and 0 . ◦ ), but there was not any difference in the outcome.All the integrations presented in Table 2 were done for m = 0 . m = 0 .
8, but similar results are expected for the other inner binarymass ratios.For systems with m = 0 .
09, there is a significant discrepancy be-tween the numerical results and the theoretical model. This is clearlydemonstrated in Figs. 2 and 3. It is easily noted in Fig. 3, which isa plot based on equations (26) and (27), that the long period andamplitude of the oscillation are larger than the ones obtained from thenumerical integrations (Fig. 2). This is due to the fact that the system pap2cmdarev.tex; 3/07/2018; 0:04; p.12 ccentricity generation in HTS √ BD − A = 0. However, for systems with m = m there is no such concern, since the long term evolution of thesystem is independent of the outer argument of pericentre g T , as canbe seen in Sect. (2.2).The rest of the numerical results are generally in good agreementwith our theory. The error is just above or drops under 10% for allmasses and inclinations with X ≥
15. Similar results are expected forother inner mass ratios.
4. Conclusion
We have extended the method of getting an estimate for the innereccentricity in hierarchical triple systems on initially circular orbitsto situations where the stars are in non-coplanar orbits. Again, theequations developed throughout this paper, give reasonable results forthe parameter ranges discussed. The addition of two more terms inthe equations for short term evolution has made our theoretical modelmore accurate, especially for systems where the perturbation to theinner binary is rather strong (the equations in the coplanar case onlyinclude the dominant P and P terms). Our future aim is to completethat type of calculation by deriving a formula for systems with eccentricouter binaries. AppendixA.
Short period components of the eccentric vector: P x21 ( t ) = 118 cos n t + 18 cos 3 n t + 316 cos (( n − n ) t − φ ) ++ 316 cos (( n + 2 n ) t + 2 φ ) + 116 cos ((3 n − n ) t − φ ) ++ 116 cos ((3 n + 2 n ) t + 2 φ ) + cos I [ 98 cos (( n − n ) t − φ ) −−
98 cos (( n + 2 n ) t + 2 φ ) + 18 cos ((3 n − n ) t − φ ) − pap2cmdarev.tex; 3/07/2018; 0:04; p.13 Nikolaos Georgakarakos
Table II. Percentage error between the averaged numerical e andequation (28) for systems with m = 0 . m = 0 .
8. The error isaccompanied by the integration time span. A dash denotes that thesystem was integrated for one outer orbital period. For each m − X pair we have three entries corresponding, from top to bottom, to aninclination of 10 ◦ , 20 ◦ and 30 ◦ respectively. m \ X
10 15 20 25 30 500.09 -219.8 -51.3 -20.6 -12 -6.6 -0.724000 62000 105000 160000 215000 510000-245.1 -50.3 -21.5 -9.3 -5.9 0.333500 81000 142000 200000 280000 650000-1383.8 -396.3 -75.5 -33.9 -15.8 -2.680000 210000 350000 550000 700000 16400000.5 14.6 7.1 3.8 2.8 1.8 0.1- - - 3500 5200 1500012.5 6.3 3.7 2.3 1.5 -0.1550 1400 2400 4100 6000 170008.2 2.8 1.5 0.7 0.1 -0.4650 1800 3200 5200 7800 225001 19.4 10.9 7.3 5.3 4.1 2.4- - - - - 1000014.8 8.2 5.5 4.2 3.2 1.6- 900 1800 2700 4200 122007.4 3.2 1.6 0.9 0.3 0.1450 1100 2100 3400 5200 150001.5 22 12.5 8.6 6.5 5.1 3- - - - - 900016.8 9 6.4 4.9 3.9 2- 800 1400 2200 3400 99008.7 2.3 1 0.8 -0.1 -0.4300 900 1800 2800 4400 130002 23.7 13.4 9.2 7 5.7 3.1- - - - - -18.1 9.8 6.8 5.1 4.2 2.2- 600 1250 2000 3000 90008.1 2.6 0.8 0.7 0.1 0.1300 800 1600 2500 3800 11000 pap2cmdarev.tex; 3/07/2018; 0:04; p.14 ccentricity generation in HTS Figure 2.
Secular resonance for a system with m = 0 . m = 0 . X = 10, I = 30 ◦ and φ = 90 ◦ . The outer orbital period is T out = 62 .
8. The graphs comefrom the numerical integration of the full equations of motion. The top graph is amagnification of the first peak of the bottom graph. −
18 cos ((3 n + 2 n ) t + 2 φ )] + cos I [ −
158 cos n t −
18 cos 3 n t ++ 1516 cos (( n − n ) t − φ ) + 1516 cos (( n + 2 n ) t + 2 φ ) ++ 116 cos ((3 n − n ) t − φ ) + 116 cos ((3 n + 2 n ) t + 2 φ )] (30) pap2cmdarev.tex; 3/07/2018; 0:04; p.15 Nikolaos Georgakarakos
Figure 3.
Secular resonance for a system with m = 0 . m = 0 . X = 10, I = 30 ◦ and φ = 90 ◦ based on equations (26) and (27). The outer orbital periodis T out = 62 .
8. Note the long period and large amplitude of the oscillation. P x22 ( t ) = 38 cos (( n − n ) t − φ ) −
38 cos (( n + 2 n ) t + 2 φ ) ++ 124 cos ((3 n − n ) t − φ ) −
124 cos ((3 n + 2 n ) t + 2 φ ) ++ cos I [ 94 cos (( n − n ) t − φ ) + 94 cos (( n + 2 n ) t + 2 φ ) ++ 112 cos ((3 n − n ) t − φ ) + 112 cos ((3 n + 2 n ) t + 2 φ )] ++ cos I [ 158 cos (( n − n ) t − φ ) −
158 cos (( n + 2 n ) t + 2 φ ) ++ 124 cos ((3 n − n ) t − φ ) −
124 cos ((3 n + 2 n ) t + 2 φ )] (31) P x31 ( t ) = cos I [ 2564 cos (3 n t + 3 φ ) − n t + φ )] ++ cos I [ 22564 cos ( n t + φ ) − n + 3 φ )] (32) P x32 ( t ) = − n − n ) t − φ ) − n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ ) −− n − n ) t − φ ) − n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ ) + pap2cmdarev.tex; 3/07/2018; 0:04; p.16 ccentricity generation in HTS
17+ cos I [ − n − n ) t − φ ) + 2132 cos ((2 n + n ) t + φ ) −− n − n ) t − φ ) + 1532 cos ((2 n + 3 n ) t + 3 φ ) −− n − n ) t − φ ) + 45512 cos ((4 n + n ) t + φ ) −− n − n ) t − φ ) + 45512 cos ((4 n + 3 n ) t + 3 φ )] ++ cos I [ 75128 cos ((2 n − n ) t − φ ) + 75128 cos ((2 n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ ) ++ 45512 cos ((4 n − n ) t − φ ) + 45512 cos ((4 n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ )] ++ cos I [ 4564 cos ((2 n − n t − φ ) − n + n ) t + φ ) −− n − n ) t − φ ) + 1564 cos ((2 n + 3 n ) t + 3 φ ) ++ 45512 cos ((4 n − n ) t − φ ) − n + n ) t + φ ) −− n − n ) t − φ ) + 15512 cos ((4 n + 3 n ) t + 3 φ )] (33) P y21 ( t ) = −
78 sin n t + 18 sin 3 n t − n − n ) t − φ ) −− n + 2 n ) t + 2 φ ) + 116 sin ((3 n − n ) t − φ ) ++ 116 sin ((3 n + 2 n ) t + 2 φ ) + cos I [ −
98 sin (( n − n ) t − φ ) ++ 98 sin (( n + 2 n ) t + 2 φ ) + 18 sin ((3 n − n ) t − φ ) −−
18 sin ((3 n + 2 n ) t + 2 φ )] + cos I [ 38 sin n t −
18 sin 3 n t −−
316 sin (( n − n ) t − φ ) −
316 sin (( n + 2 n ) t + 2 φ ) ++ 116 sin ((3 n − n ) t − φ ) + 116 sin ((3 n + 2 n ) t + 2 φ )] (34) P y22 ( t ) = −
158 sin (( n − n ) t − φ ) + 158 sin (( n + 2 n ) t + 2 φ ) ++ 124 sin ((3 n − n ) t − φ ) −
124 sin ((3 n + 2 n ) t + 2 φ ) + pap2cmdarev.tex; 3/07/2018; 0:04; p.17 Nikolaos Georgakarakos + cos I [ −
94 sin (( n − n ) t − φ ) −
94 sin (( n + 2 n ) t + 2 φ ) ++ 112 sin ((3 n − n ) t − φ ) + 112 sin ((3 n + 2 n ) t + 2 φ )] ++ cos I [ −
38 sin (( n − n ) t − φ ) + 38 sin (( n + 2 n ) t + 2 φ ) ++ 124 sin ((3 n − n ) t − φ ) −
124 sin ((3 n + 2 n ) t + 2 φ )] (35) P y31 ( t ) = 2564 sin (3 n t + 3 φ ) − n t + φ ) ++ cos I [ 7564 sin ( n t + φ ) − n + 3 φ )] (36) P y32 ( t ) = 3364 sin ((2 n − n ) t − φ ) + 3364 sin ((2 n + n ) t + φ ) ++ 1564 sin ((2 n − n ) t − φ ) + 1564 sin ((2 n + 3 n ) t + 3 φ ) −− n − n ) t − φ ) − n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ ) ++ cos I [ 51128 sin ((2 n − n ) t − φ ) − n + n ) t + φ ) ++ 75128 sin ((2 n − n ) t − φ ) − n + 3 n ) t + 3 φ ) −− n − n ) t − φ ) + 45512 sin ((4 n + n ) t + φ ) −− n − n ) t − φ ) + 45512 sin ((4 n + 3 n ) t + 3 φ )] ++ cos I [ − n − n ) t − φ ) − n + n ) t + φ ) ++ 1532 sin ((2 n − n ) t − φ ) + 1532 sin ((2 n + 3 n ) t + 3 φ ) ++ 45512 sin ((4 n − n ) t − φ ) + 45512 sin ((4 n + n ) t + φ ) −− n − n ) t − φ ) − n + 3 n ) t + 3 φ )] ++ cos I [ − n − n ) t − φ ) + 45128 sin ((2 n + n ) t + φ ) ++ 15128 sin ((2 n − n ) t − φ ) − n + 3 n ) t + 3 φ ) ++ 45512 sin ((4 n − n ) t − φ ) − n + n ) t + φ ) − pap2cmdarev.tex; 3/07/2018; 0:04; p.18 ccentricity generation in HTS − n − n ) t − φ ) + 15512 sin ((4 n + 3 n ) t + 3 φ )] (37)The complete secular equations of motion: dx S dτ = 5(1 − e ) sin I − x (1 − x − y ) y S + 516 α e T (1 − e ) (1 − x −− y ) [sin g T cos I (4 + 3( x + y ) − I (1 − x + 6 y )) −− − x − y ) sin I cos I sin g T + 2(3 + 5 sin I )( y sin g T ×× cos I + x S y S cos g T ) + 20 sin I cos Iy sin g T −
70 sin I ( y ×× sin g T cos I + x S y S cos g T )] + 516 αβ e T (1 − e ) [ y sin g T (4 ++3( x + y ) − I (1 − x + 6 y )) + 10 cos I ( y sin g T ×× cos I + x S y S cos g T )(1 − x + 6 y ) + 10(1 − x − y ) y ×× sin g T (2 cos I − sin I )] + 516 α e T (1 − e ) − x − y ) ×× [ y sin g T cos I (4 + 3( x + y ) − I (1 − x + 6 y )) ++10( y sin g T cos I + x S y S cos g T ) cos I (1 − x + 6 y ) ++10(1 − x − y ) y sin g T (2 cos I − sin I ) cos I ] −− β (1 − e ) y S cos I (1 − x + 4 y ) − − e ) y S ×× − x + 3 y (1 − x − y ) (38) dy S dτ = − − e ) sin I x S y (1 − x − y ) y S + 516 α e T (1 − e ) (1 −− x − y ) [ − cos g T (4 + 3( x + y ) − I (1 − x ++6 y )) − I )( x S y S sin g T cos I + x cos g T ) −−
20 sin I cos Ix S y S sin g T ] + 516 αβ e T (1 − e ) [ − x S y S sin g T ×× (4 + 3( x + y ) − I (1 − x + 6 y )) −
10 cos I ( x S ×× y S sin g T cos I + x cos g T )(1 − x + 6 y ) − − x −− y ) x S y S sin g T (2 cos I − sin I )] + 516 α e T (1 − e ) ×× − x − y ) [ − x S y S sin g T cos I (4 + 3( x + y ) − I ×× (1 − x + 6 y )) − x S y S sin g T cos I + x cos g T ) cos I (1 −− x + 6 y ) − − x − y ) x S y S sin g T (2 cos I − sin I ) × pap2cmdarev.tex; 3/07/2018; 0:04; p.19 Nikolaos Georgakarakos × cos I ] + β (1 − e ) x S cos I (1 − x + 4 y ) + 1(1 − e ) x S ×× − x + 3 y (1 − x − y ) (39) dg T dτ = 12 β (1 − e ) [4 + x + 11 y − I (1 − x + 4 y )] ++ 1(1 − x − y ) (1 − e ) cos I (1 − x + 4 y ) − αβ ×× e (1 − e ) e T [( y S sin g T cos I + x S cos g T )(4 + 3( x + y ) −− I (1 − x + 6 y )) − − x − y ) sin I cos I ×× y S sin g T ] + ( 516 α e T (1 − e ) − x − y ) + 516 αβ ×× e T (1 − e ) cos I )[ − y S sin g T (4 + 3( x + y ) − I (1 −− x + 6 y )) − y S sin g T cos I + x S cos g T ) cos I (1 − x ++6 y ) − − x − y ) y S sin g T (2 cos I − sin I )] (40) de T dτ = 516 αβ (1 − e ) [( y S cos g T cos I − x S sin g T )(4 + 3( x + y ) −− I (1 − x + 6 y )) − − x − y ) y S sin I ×× cos I cos g T ] (41) dIdτ = − x S ˙ x S + y S ˙ y S (1 − x − y ) ( 1tan I (1 − x − y ) + β sin I (1 − e ) ) −− e T ˙ e T (1 − e ) ( 1 β sin I (1 − x − y ) + 1tan I (1 − e ) ) (42) Acknowledgements
The author is grateful to Prof. Douglas Heggie for all the useful dis-cussions on the context of this paper. The author also thanks SeppoMikkola, who kindly provided the code for integrating hierarchicaltriple systems. pap2cmdarev.tex; 3/07/2018; 0:04; p.20 ccentricity generation in HTS References
Eggleton, P. and Kiseleva, L.: 1995, An empirical condition for stability ofhierarchical triple systems, Astroph. J. , 640-645.Ford, E.B., Kozinsky, B. and Rasio, F.A.: 2000, Secular evolution of hierarchicaltriple systems, Astroph. J. , 385-401.Georgakarakos, N.: 2002, Eccentricity generation in hierarchical triple systems withcoplanar and initially circular orbits, MNRAS , 559-566.Kiseleva, L., Eggleton, P. and Mikkola, S.: 1998, Tidal friction in triple stars,MNRAS , 292-302.Kozai, Y.: 1962, Secular perturbations of asteroids with high inclination andeccentricity, Astron. J. , 591-598.Krymolowski, Y. and Mazeh, T.: 1999, Studies of multiple stellar systems -II.Second-order averaged Hamiltonian to follow long-term orbital modulations ofhierarchical triple systems, MNRAS , 720-732.Marchal, C.: 1990, The Three-Body Problem, Elsevier Science Publishers, theNetherlands.Mikkola S.: 1997, Practical symplectic methods with time transformation for thefew-body problem, CeMDA , 145-165.Press, W., Teukolsky S., Vetterling W. and Flannery B.: 1996, Numerical RecipesIn Fortran 77, 2nd ed., Cambridge Univ. Press, New York., 145-165.Press, W., Teukolsky S., Vetterling W. and Flannery B.: 1996, Numerical RecipesIn Fortran 77, 2nd ed., Cambridge Univ. Press, New York.