Issue 
A&A
Volume 545, September 2012



Article Number  A88  
Number of page(s)  13  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201219183  
Published online  12 September 2012 
Dynamical evolution of the Gliese 436 planetary system
Kozai migration as a potential source for Gliese 436b’s eccentricity
^{1}
UJFGrenoble 1 / CNRSINSU, Institut de Planétologie et d’Astrophysique de
Grenoble (IPAG) UMR 5274,
38041
Grenoble,
France
email: Herve.Beust@obs.ujfgrenoble.fr
^{2}
European Organization for Astronomical Research in the Southern
Hemisphere (ESO), Casilla 19001, Santiago 19, Chile
Received:
7
March
2012
Accepted:
1
August
2012
Context. The closein planet orbiting GJ 436 presents a puzzling orbital eccentricity (e ≃ 0.14) considering its very short orbital period. Given the age of the system, this planet should have been tidally circularized a long time ago. Many attempts to explain this were proposed in recent years, either involving abnormally weak tides, or the perturbing action of a distant companion.
Aims. In this paper, we address the latter issue based on Kozai migration. We propose that GJ 436b was formerly located further away from the star and that it underwent a migration induced by a massive, inclined perturber via Kozai mechanism. In this context, the perturbations by the companion trigger highamplitude variations to GJ 436b that cause tides to act at periastron. Then the orbit tidally shrinks to reach its present day location.
Methods. We numerically integrate the 3body system including tides and general relativity correction. We use a modified symplectic integrator as well as a fully averaged integrator. The former is slower but accurate to any order in semimajor axis ratio, while the latter is first truncated to some order (4th) in semimajor axis ratio before averaging.
Results. We first show that starting from the presentday location of GJ 436b inevitably leads to damping the Kozai oscillations and to rapidly circularizing the planet. Conversely, starting from 5–10 times further away allows the onset of Kozai cycles. The tides act in peak eccentricity phases and reduce the semimajor axis of the planet. The net result is a twofold evolution, characterized by two phases: a first one with Kozai cycles and a slowly shrinking semimajor axis, and a second one once the planet gets out of the Kozai resonance characterized by a more rapid decrease. The timescale of this process appears in most cases much longer than the standard circularization time of the planet by a factor of 50 or above.
Conclusions. This model can provide a solution to the eccentricity paradox of GJ 436b. Depending on the various orbital configurations (mass and location o the perturber, mutual inclination, etc.), it can take several Gyr to GJ 436b to achieve a full orbital decrease and circularization. According to this scenario, we could be witnessing today the second phase of the scenario, where the semimajor, axis is already reduced while the eccentricity is still significant. We then explore the parameter space and derive in which conditions this model can be realistic given the age of the system. This yields constraints on the characteristics of the putative companion.
Key words: planetary systems / planets and satellites: dynamical evolution and stability / celestial mechanics / stars: individual: Gliese 436 / methods: numerical / planetstar interactions
© ESO, 2012
1. Introduction
The Mdwarf GJ 436 has been the subject of growing interest in recent years. This star is known to host a closein Neptunemass planet (Gl 436b Butler et al. 2004) that was furthermore observed to undergo transit (Gillon et al. 2007).
The monitoring of the transits of GJ 436b helped constraining its orbital solution. Noticeably, it appears to have a significant nonzero eccentricity: e = 0.14 ± 0.1 (Maness 2007), furthermore refined to e = 0.14 ± 0.01 (Demory et al. 2007). This eccentricity is abnormally high for a small period planet (P ≃ 2.64 days Ballard et al. 2010). With such a small orbital period, tidal forces are expected to circularize the orbit within much less than the present age of the system (1–10 Gyr Torres 2007). Tidal forces seem indeed to be at work in this system. The detection of the secondary transit of the planet enabled Demory et al. (2007) to derive a brightness temperature of T = 712 ± 36 K, with flux reradiated across the dayside only. Conversely, a stellar irradiation/thermal reradiation balance leads for to an equilibrium temperature for the planet T_{eq} = 642 K, assuming T_{eff} = 3350 K for GJ 436 and zero albedo for the planet. According to Deming et al. (2007), the temperature difference indicates tidal heading of GJ 436b, but this could alternatively be due to the 8 μm sampling the atmosphere in a hot band pass, if the planetary atmosphere does not radiate like a blackbody.
Many theories were proposed in recent years to account for the residual eccentricity of GJ 436b. The most straightforward one is that the tides are sufficiently weak and/or the age of the system is small enough to allow a regular tidal circularization not to be achieved yet. This idea was proposed by Mardling (2008) for GJ 436b, after Trilling (2000) had suggested that this accounts for any closein exoplanet observed with significant eccentricity. It is in fact related to our poor knowledge of the planet’s tidal dissipation Q_{p}. Q_{p} is a dimensionless parameter related to the rate of energy dissipated per orbital period by tidal forced oscillations (Barker & Ogilvie 2009). The smaller Q_{p}, the more efficient the tidal dissipation is. To lowest order, the circularization timescale t_{circ} of an exoplanet reads (Goldreich & Soter 1966; Jackson et al. 2008) (1)where M_{∗} is the stellar mass, a is the planet’s orbital semimajor axis, m_{p} its mass and R_{p} its radius. Assuming Q_{p} = 10^{5}, R = 27 600 km, a = 0.0287 AU, m = 23.2 M_{⊕} and M = 0.452 M_{⊙} (Mardling 2008), this formula gives t_{circ} = 4.7 × 10^{7} yr, which is obviously less than the age of the system. But Q_{p} is very badly constrained. For Neptune, it is estimated between ~10^{4} and 3.3 × 10^{5} (Banfield & Murray 1992; Zhang & Hamilton 2008). We can thus consider 10^{5} as a standard likely value for GJ 436b, but Mardling (2008) argues that Q_{p} could be as high as a few 10^{6}. In this context, assuming the lower bound (1 Gyr) for the age of the system, the circularization could not be achieved yet. This depends however on the starting eccentricity at time zero. Wisdon (2008) performed analytical calculations of energy dissipation rates for synchronized bodies at arbitrary eccentricities and obliquities. His concluded that for a large enough starting eccentricity (≳ 0.4, as will be the case below), the circularization timescale should be significantly reduced with respect to Eq. (1). In this context, the conclusions by Mardling (2008) may no longer hold. This shows also that conclusions drawn on asymptotic formulas like Eq. (1) must be treated with caution. Numerical work is required.
Alternatively, many authors proposed that the eccentricity of GJ 436b is sustained by perturbations by an outer massive planet, despite a standard Q_{p} value. Such a long period companion was initially suspected as Maness (2007)’s radial velocity data revealed a linear drift. More data did not confirm the trend, which is now believed to be spurious. Montagnier et al. (in prep.) give instead observational uppers limits for the putative companion that rule out any companion in the Jupiter mass range up to a few AUs.
Perturbers can be either resonant or nonresonant. Looking for a planet locked in meanmotion resonance with GJ 436b appears a promising idea, as resonant configuration usually trigger larger eccentricity modulations. This was suggested by Ribas et al. (2008), who fitted in the residuals of the radial velocity data a ~5 M_{⊕} planet locked in 2:1 mean motion resonance with GJ 436b. However, this additional planet was invalidated furthermore. Dynamical calculations by Mardling (2008) showed that this planet cannot sustain the eccentricity of GJ 436b unless Q_{p} ≳ 10^{6}. Moreover, Alonso et al. (2008) showed that a ~5 M_{⊕} planet in 2:1 resonance with GJ 436b should trigger transit time variations that should have been detected yet. Such variations are unseen today with a high level of confidence (Pont et al. 2009). Actually the constraints deduced from radial velocity monitoring, transit time and geometry monitoring (Ballard et al. 2010) almost rule out an additional large planet locked in a loworder meanmotion resonance with GJ 436b (P < 8.5 days), and to have a Jovian like planet up to a few AUs.
Conversely, Tong & Zhou (2009) investigated the effects of secular planetary perturbations in eccentricity pumping, considering both nonresonant (secular) and resonant configurations. They found some perturbers configurations that could account for the present eccentricity of GJ 436b while still being compatible with the observational constraints, but with no tides. Incorporating the tides leads inevitably to damp all the eccentricity modulation and to circularize the orbit. Meanwhile, Batygin et al. (2009) reanalyzed carefully this issue. They found that as expected, in most cases the eccentricity of GJ 436b quickly falls to zero thanks to tidal friction, but that this effect can be considerably delayed if the 2 planet system lies initially in a specific configuration where the eccentricities of the two planets are locked at stationary points in the secular dynamics diagram. In this case, they show that the circularization time can be as high as ~8 Gyr despite a standard Q_{p} value.
To date, the model by Batygin et al. (2009) is the only one able to explain the high present day eccentricity of GJ 436b with a standard Q_{p} value. This model could appear as not generic, as it requires a specific initial configuration, but Batygin & Morbidelli (2011) showed that in the framework of Hamiltonian planetary dynamics with additional dissipative forces (which is the case here), these points tend to behave as attractors. Hence various initial configurations can lead to reach such a point with a further dynamical evolution like described by Batygin et al. (2009). As is shown below, in the model presented here, we exactly encounter such a configuration where many different routes can lead to such a stationary point (the transition between Phase 1 and Phase 2, see Sect. 4).
The purpose of this paper is to present an alternate model based on Kozai mechanism, assuming again a distant perturber. This model can be viewed as complementary to the Batygin et al. (2009) study. Kozai mechanism is a major dynamical effect in noncoplanar systems that can trigger eccentricity modulations up to very high values. We make a review on this effect, without and with tidal effects taken into account and describe our approach in Sects. 2 and 3 respectively. In Sect. 4, we present an application to the case of GJ 436. We first present calculations starting from the present day orbit of GJ 436b, with the result that even Kozai mechanism cannot overcome the damping effect of tidal friction. Then we present other simulations where GJ 436b is initially put much further away from the star than today. We show that Kozai mechanism pumps GJ 436b’s eccentric regularly up to high values where tides become active at periastron. The result is a decay of the orbit that drives it to its present day location but on a timescale that can be considerably longer than the standard circularization time. Hence even after several Gyr, GJ 436b’s eccentricity can still be significant. In Sect. 5, we describe a parametric study of this scenario and derive in which conditions it is compatible with the present day situation of GJ 436b. Our conclusions are presented in Sect. 6.
2. Kozai mechanism
2.1. General features
Kozai mechanism was first described by Kozai (1962), initially for comets on inclined orbit (≳50°) with respect to the ecliptic. Under the effect of secular planetary perturbations (mainly arising from Jupiter), their orbit is subject to a periodic evolution that drives it to lower inclination but very high eccentricity, while the semimajor axis remains constant. As pointed out by Bailey et al. (1992), this mechanism is responsible for the origin of most sungrazing comets in our solar system. In the circular restricted 3body problem, the Kozai Hamiltonian that describes this secular evolution is obtained by a double averaging of the original Hamiltonian over the orbital motions of the planet and of the particle (Kinoshita & Nakai 1999).
Assuming zero eccentricity for the primaries, the secular motion of the particle is characterized, to quadrupolar approximation (expansion of the secular Hamiltonian to 2nd order in powers of the semimajor axis ratio) by the following constants of motions (Kinoshita & Nakai 1999; Krymolowski & Mazeh 1999; Ford et al. 2000; Beust & Dutrey 2006) where e stands for orbital the eccentricity of the particle, i for its inclination with respect to the orbital plane of the primaries and ω for the particle’s argument of periastron relative to this plane. The first constant (F) is the reduced Hamiltonian itself and H is the reduced vertical component of the angular momentum, which is a secular constant thanks to the axial symmetry of the averaged problem. Combining both we can eliminate the inclination (4)Contour plots of this expression in (ω,e) space for a given H value are called Kozai diagrams. The Kozai mechanism regime is characterized by small values of H (high initial inclinations). Figure 1 shows such a diagram corresponding to H = 0.05. We note the two islands of libration around ω = 90° and ω = 270° (or −90°) characteristic for Kozai resonance. In a pure 3body Newtonian dynamics, the particle just secularly moves along one of these curves. It is then obvious to see that its evolution is characterized by strong periodic eccentricity increases. In the same time e increases, i decreases to keep H constant. If for a given particle ω librates around ± 90°, it is said to be trapped in Kozai resonance. But even if ω circulates, its evolution is characterized by the same eccentricity/inclination modulation described above, the maximum e values being reached when ω = ± 90°.
As quoted above, this regime is valid for small H values. If the particle has initial zero eccentricity, it can be shown analytically (Kinoshita & Nakai 1999) that the minimum inclination required for Kozai mechanism to start is arccos ≃ 39.23°. Alternatively, it can operate at smaller inclination if the initial eccentricity is larger. But a some point in the secular evolution, it will reach an inclination above this threshold.
This dynamical mechanism is also active in hierarchical triple systems (Harrington 1968; Söderhjelm 1982; Eggleton et al. 1998; Krymolowski & Mazeh 1999; Ford et al. 2000), the inner orbit being subject to a similar evolution under secular perturbations by the outer body. This is particularly relevant for triple systems because they are often noncoplanar. Similarly, it is expected to be active on planetary orbits in binary star systems, as soon as the orbit of the planets are significantly non coplanar with that of the binary (i ≳ 40°, see above). This can affect the stability of the planetary orbits, especially if several planets are present (Innanen et al. 1997; Malmberg et al. 2007). It was actually invoked to explain the high eccentricity of some extrasolar planets in binary systems (Holman et al. 1997; Libert & Tsiganis 2009).
2.2. Kozai migration
This mechanism could also apply to the GJ 436 case, as a way to sustain the eccentricity of GJ 436b. But tides inevitably lead to a circularization within the standard timescale (see simulations below and also application to the TyCra system in Beust et al. 1997).
To render Kozai mechanism active despite the presence of tides, one must lower the strength of the tides. This can be done assuming a higher Q_{p} as suggested by Mardling (2008) (see Introduction), but this is even not enough. Even with an arbitrary large Q_{p}, there is still forced apsidal precession due to the planetary and stellar quadrupole moments and to general relativity (see Sect. 3). With the present day semimajor axis of GJ 436b, this is enough to override the Kozai mechanism.
Then only way to damp all these effects to allow the onset of Kozai cycles is to assume a larger semimajor axis for the planet during a first evolutionary phase. The purpose of the present paper is to investigate this idea. Of course, the present orbital semimajor axis of GJ 436b is well constrained from radialvelocity monitoring, but our idea is to suggest that it used to be larger in the past, typically 5–10 times larger than today. In this context, all perturbing effects listed are too weak to prevent the onset of Kozai cycles. Then when the eccentricity increases, the periastron drops, because as long as tides are not active, the semimajor axis remains a secular invariant. Hence tidal friction starts at periastron. This leads to a gradual decrease of the orbit (i.e., the semimajor axis), which causes it to shrink and to finally reach its present location. This process is called Kozai migration. It has been described by Eggleton & KiselevaEggleton (2001); Wu & Murray (2003) and more recently by Fabrycky & Tremaine (2007). It is viewed today as one of the various migration mechanism that can act in planetary systems. Our aim is to investigate this scenario in the context of GJ 436b. We show that although it cannot prevent the final circularization of the planet, under suitable assumptions it can considerably slow down the circularization process, thus explaining how the eccentricity of GJ 436b could still be significant a few Gyr after the formation of the system in spite of strong tidal friction.
3. Tidal forces and Nbody integration
3.1. Basic tidal forces
Our goal is to test the scenario of Kozai migration to the GJ 436b case. This needs to be done via numerical integration. The various tidal forces are described in Mardling & Lin (2002). We recall briefly this theory here. Let us consider i) a star with mass m_{s}, radius R_{s}, an apsidal constant k_{s} and a dissipation parameter Q_{s}; ii) a planet with mass m_{p}, radius R_{p}, an apsidal constant k_{p} and a dissipation parameter Q_{p}. The tidal forces acting on the planet are

1.
the acceleration due to the quadrupole moment of the planetdue to the distorsion produced by the presence of the star:(5)where Ω_{p} is the spin vector of the planet, r is the distance between the planet and the star, and u = r/r, a unit vector parallel to the radius vector of the orbit;

2.
the acceleration produced by the tidal damping of the planet (6)where a is the semimajor axis and n the mean motion of the orbit.
Of course we have symmetric accelerations γ_{QD,s} and γ_{TF,s} for the star. We must add to these forces the relativistic postNewtonian acceleration (7)where η = m_{p}m_{s}/(m_{p} + m_{s})^{2} and c is the speed of light. Finally, the equation of motion describing the evolution of the orbit of the planet is (8)We must also add an equation governing the evolution of the spin vectors. For the planet we have (Mardling & Lin 2002) (9)and a symmetric equation for . Note that we typically have γ_{QD,s} ≪ γ_{QD,p} and γ_{TF,s} ≪ γ_{TF,p} by a factor m_{p}/m_{s}, showing that the star itself is very little affected. To lowest order for a synchronized planet, we have γ_{QD,p}/γ_{TF,p} ≃ 11Q_{p}/6. With Q_{p} ≳ 10^{5}. We see that γ_{QD,p} dominates the orbital evolution. However, γ_{TF,p} (the damping term) is responsible for the circularization and must be retained. Also, we have to lowest order (10)with the numerical values for the parameters of GJ 436b quoted in Mardling (2008), this ratio is ≃1.6. This shows that relativistic effects are expected to be comparable to tides and must not be neglected. However, as pointed out by Mardling & Lin (2002), the only secular effect of the relativistic potential is to produce an apsidal advance. This also holds for γ_{QD,p} as both effects are conservative. They do not affect the circularization process due to γ_{TF,p}.
3.2. Inclusion into Nbody integration and averaging
As we want to investigate the combined effect of tides and planetary perturbations, these forces must be included into a Nbody integration scheme. Our basic scenario is an inner planet (GJ 436b) subject to tidal interaction from the star, and perturbed by an outer, inclined and more massive planet (with mass m′ and semimajor axis a′) that is itself not subject to tides.
A major difficulty is that we have to deal with phenomena with very different timescales. In the context of GJ 436, the period of Kozai cycles is typically ~10^{5} yr, while the tidal circularization acts on timescales ≳ 10^{8} yr. Hence an efficient Nbody integrator is required.
A standard way to overcome this difficulty is to average all interaction terms over the orbital motions of both bodies. This can concern the interaction 3body Hamiltonian between the planet and the perturber, as well as the tidal terms. This way only the secular part of the perturbations is retained.
The averaging of the tides and of the relativistic terms is done only over the orbital period of the inner planet. To do this, it is convenient to derive the effect of the tides on the specific angular momentum vector and the RungeLenz vector (Mardling & Lin 2002; Fabrycky & Tremaine 2007; Eggleton & KiselevaEggleton 2001) (11)Each additional acceleration γ will induce a rate of change of e and h given by (12)These equations are then averaged over the orbital period of the inner planet. Using Hansen coefficients, this can be done in closed form. Explicit formulas are given by Mardling & Lin (2002).
The averaging of the interaction Hamiltonian with the outer perturber over the orbital period of both bodies is more complex, as this cannot be done in closed form. One has to expand the corresponding terms in ascending powers of the semimajor axis ratio a/a′ between the planet and the perturber (using Legendre polynomials). Assuming a/a′ ≪ 1, it is truncated to some given order and then averaged over both orbits. Wu & Murray (2003) and Fabrycky & Tremaine (2007) truncate the Hamiltonian to second order (quadrupolar), while Ford et al. (2000) recommend third order (octupolar). Mardling & Lin (2002) give explicit formulas up to third order (for coplanar orbits though). As explained below, the level of accuracy required by our study forced us retain terms up to 4th order in a/a′, i.e. one level beyond octupolar. We used Maple to derive analytic formulas up to this level.
Another way to handle the 3body interaction is to use symplectic integration. The main difficulty is to combine it with tidal forces. The theory of symplectic integration is based on the Hamilton equations of motion that apply to any Nbody problem. Its background is for instance described in Saha & Tremaine (1992) and Chambers (1999). The key idea is to split the Hamiltonian H into pieces, say H = H_{A} + H_{B}, each of them being exactly integrable, and then to integrate each part separately alternatively. Most common symplectic integrators (Chambers 1999; Levison & Duncan 1994) are based on the following second order scheme: at each timestep τ, we must:

1.
advance H_{B} by a half step τ/2;

2.
advance H_{A} for the full time step τ;

3.
readvance H_{B} by τ/2.
In planetary and hierarchical Nbody problems, H_{A} corresponds typically to a collection of independent Keplerian Hamiltonians and H_{B} to the disturbing function. We will name this scheme B(τ/2)A(τ)B(τ/2). Our goal here is to include tides to it, i.e. dissipative forces. Of course with dissipation the integration, as well as the real problem, is no longer symplectic. This does not prevent the integration method from being efficient, at least for the Nbody forces. Typically, with the second order scheme described above, a timestep of ~1/20 of the smallest orbital period is enough in standard problems to ensure a conservation of the energy to ~10^{6} relative level (Beust 2003). Including dissipative forces leads to a non preservation of the energy, but the same timestep can be kept. There have been actually many successful attempts in recent years to include dissipative forces into symplectic integrators (Malhotra 1994; Kehoe et al. 2003; Ćuk & Burns 2004; Marzari & Weidenschilling 2006).
Fig. 2 Secular evolution of the eccentricity of GJ 436b starting from the presentday orbit and assuming a standard perturber (m′ = 0.1 M_{Jup}, T′ = 10 yr) in different conditions (see text): a) perturbing planet, but no tides nor relativity; b) tides and relativity but no perturbing planet; c) perturbing planet, relativity but no tides; d) perturbing planet, relativity and tides. 
The most straightforward way to proceed is to include the integration of the various forces given above into the B(τ/2) substeps (the disturbing part of the integration). However this turns out not to be efficient in our case. The reason is that the tidal forces are highly dependant on the distance r between the planet and the star (∝ r^{4} and even steeper). Therefore, when considering significantly eccentric orbits, keeping the standard symplectic timestep leads to an insufficient mapping of the tidal forces close to periastron.
To overcome this difficulty, we used averaged formulas only for the tidal and relativistic terms. Then we add to our symplectic scheme two substeps T(τ/2) of integration of these averaged equations for e and h, in the order B(τ/2)T(τ/2)A(τ)T(τ/2)B(τ/2):

1.
advance H_{B} by a half step τ/2;

2.
integrate the averaged tidal equations by τ/2;

3.
advance H_{A} for the full time step τ;

4.
reintegrate the averaged tidal equations by τ/2;

5.
readvance H_{B} by τ/2.
At each T(τ/2) substep, we integrate simultaneously the evolution of the spin vectors of the bodies via Eq. (9). The advantage of this approach is that it avoids the use of truncated expansions, while still being efficient. We modified this way our symplectic integrator HJS dedicated to hierarchical systems (Beust 2003). All computations presented in Sect. 4 have been done using this integrator. Its use is nevertheless more computing time consuming than a fully averaged integrator. This is why we also built a fully averaged integrator for our parametric study (Sect. 5).
To check the validity of our integration method, we also made some comparisons with integrations performed with a conventional Bulirsh & Stoer integrator. Of course the latter integrations were much more time consuming so that they were done over a much less extended timespan. In all cases the agreement appeared excellent. An example is shown in Fig. 3 (eccentricity plot) where we have overplotted in red the result of the conventional integration. Both curves are almost identical, so that the red curve barely appears under the black one.
4. Application to Gliese 436
4.1. Integration starting from the present state
The first element we want to test is how planetary perturbations may affect the present orbit of GJ 436b. As mentioned above, given the present semimajor axis a = 0.0287 AU, tides are expected to dominate the dynamical evolution. We therefore performed four runs starting from the presentday orbital configuration with different assumptions:

run (a): with a distant perturber, but not taking into account thetides nor the general relativistic postNewtonian correction(hereafter GR), i.e. using the standard HJS code;

run (b): taking into account tides and GR, but with no perturbing planet (run (b) in Fig. 2);

run (c): with the same perturbing planet as in run (a), but taking now into account the relativistic postNewtonian correction, but no tides;

run (d): with the same planet taking into account both GR and tides.
In all cases, the perturbing planet was chosen with mass m′ = 0.1 M_{Jup} and semimajor axis a′ = 3.583 AU (i.e., with orbital period T′ = 10 yr) as a typical convenient set of parameters (see Sect. 5), compatible with the constraints of Montagnier et al. (in prep.) based on radial velocity residuals and photometry. In addition we fix the initial tilt angle i_{0} between the two orbits to 75° to be able to generate a strong Kozai resonance. We assume the stellar, planetary and tidal parameters given by Mardling (2008) and listed in the Introduction.
In Fig. 2 we show the secular evolution of the eccentricity of GJ 436b over 2 × 10^{7} yr in the four cases. In case (a), we see as expected a large amplitude modulation characteristic for the Kozai resonance created by the perturbing planet. In case (b), we note a gradual circularization with a timescale ≃ 2 × 10^{8} yr. In case (c) (same as (a) but with GR taken into account), there is no Kozai resonance anymore. The relativistic precession actually smoothes the Kozai mechanism. The eccentricity is basically constant with high frequency, small amplitude oscillations. In run (d) (with tides and relativity), the high frequency oscillation is superimposed to a tidal decrease at the same rate as in case (b). The main conclusion we can derive from these runs is that the perturbing planet could in principle induce a strong Kozai resonance, but that tides and GR override this effect. Even GR alone (run (c)) is enough to cancel out the Kozai resonance. The reason is that the secular precession induced by GR is larger than that induced by the Kozai resonance. As shown by Fabrycky & Tremaine (2007), GR adds an extra term to the secular precession of the argument of periastron ω, and subsequently, the maximum eccentricity of the Kozai cycles is significantly reduced. This can be understood considering Eq. (4). With an extra precession for ω, the cos2ω term in the averaged Hamiltonian F can be considered as rapidly varying, so that it secularly averages to zero. Only the first term can be kept in F. This makes the Hamiltonian independent of all angles, making all related actions constant, and hence e and i are constant. This is what we get here.
If we consider the dissipative case with tides, the net result in our case is that the tidal circularization process remains almost unchanged (run (d)) with respect to the case with no perturbing planet (run (b)).
The Kozai resonance appears thus unable to sustain the high present eccentricity of GJ 436b over the age of the system. The only ways to solve this paradigm are either to increase the strength of the Kozai resonance either to lower the tides and the relativistic effects. The former way would imply to increase the mass of the perturbing planet by at least two orders of magnitude; this would render it incompatible with the constraints of Montagnier et al. (in prep.). The latter way can be achieved assuming a larger semimajor axis for GJ 436b. Of course the present semimajor axis is well constrained, but we may assume that it used to be larger in the past. This is why we come to the idea of Kozai migration.
4.2. The Kozai migration scenario
4.2.1. General description
Fig. 3 Secular evolution over 2 × 10^{7}yr of GJ 436b assuming an initial value a = 0.287 AU and a perturber with m′ = 0.1 M_{Jup} and a′ = 3.583 AU (T′ = 10 yr): eccentricity (top left), semimajor axis (top right), rotation speed with respect to the orbital meanmotion (Ω/n, bottom left), and obliquity angle (bottom right). On the eccentricity plot we overplot in red the result of the same integration, but using a conventional Bulirsh & Stoer integrator. 
Fig. 4 Same evolution as in Fig. 3, but over 2 × 10^{10}yr: eccentricity (top left), semimajor axis (top right), rotation speed with respect to the orbital meanmotion (Ω/n, bottom left), mutual inclination between the two orbits (bottom right). On all plots, the black curve corresponds to our symplectic integration. The color curves on the eccentricity plot are calculated with a fully averaged integrator, with interaction Hamiltonian truncated up to various orders: 2nd order (blue), 3rd order (green), and 4th order (red). 
We describe now into details a typical run where we assume an initial semimajor axis for GJ 436b significantly larger than the present value, and where both GR and tides are taken into account, with unchanged tidal parameters with respect to Sect. 4.1. We assume for the perturber the same parameters as in Sect. 4.1 (a′ = 3.583 AU, m′ = 0.1 M_{Jup}, i_{0} = 75°), but we take as initial semimajor axis for GJ 436b a_{0} = 0.287 AU (i.e, 10 times its presentday value).
Figure 3 shows the secular evolution of GJ 436b in this context over 2 × 10^{7} yr. We see that the eccentricity is subject to a large amplitude modulation characteristic for Kozai resonance. With a_{0} = 0.287 AU, apsidal precession arising from GR and tides is no longer able to override the Kozai resonance. But at peak eccentricity (e ≃ 0.8) in the Kozai cycles, the periastron drops to 0.05–0.06 AU, i.e., only twice the presentday semimajor axis. As a result, tides are at work in peak eccentricity phases. This can be seen in Fig. 3. We note a rapid reduction of the obliquity of GJ 436b relative to its orbital plane (initially fixed to 30°), so that after a few 10^{7} yr, the spin axis of GJ 436b remains perpendicular to its orbital plane. We also note a spinning up of the planet’s rotation. In Fig. 3 (bottomleft), we plot as a function of time the ratio of the rotation velocity GJ 436b to its orbital meanmotion (Ω/n), so that 1 means synchronism. This ratio was initially fixed to 2, and we note a gradual increase up to 10 at t = 2 × 10^{7} yr. This reveals a spinning up of the planet, as the orbital meanmotion (i.e. the semimajor axis) is not subject to so drastic changes in the same time. Finally, we note a very small decrease of the semimajor axis over the timespan considered. These effects are clearly due to the tides at work in high eccentricity phases. This shows up obviously in the staircaselike shape of the temporal evolution of the parameters displayed. Actually they are subject to some evolution only in the peak eccentricity phases, and remain unchanged inbetween. The spinning up of the planet is also a consequence of the tides in high eccentricity phases. The rotation of GJ 436b tries indeed to tidally synchronize with the orbital angular velocity at periastron, which is very different from the meanmotion due to the high eccentricity. With a peak eccentricity of ~0.8, synchronism with angular velocity at periastron means actually Ω/n ≃ 16. Finally, the semimajor axis decrease reveals a loss of orbital energy due to tidal friction.
We come now to describing the same simulation, but over a much longer time span, i.e. 20 Gyr (Fig. 4). This timescale is intentionally taken larger than the age of the system and even of the Universe to show the possible future evolution. We display as in Fig. 3 the eccentricity, the semimajor axis, the Ω/n ratio, and the mutual inclination between the two orbits. On the eccentricity plot, the color curves are calculated using various averaged integrators (see below). The one corresponding to our symplectic integration in the black one. We clearly see a twofold evolution, first before 10 Gyr and then after. In the first phase, Kozai resonance is still active as can be seen from the high amplitude eccentricity oscillations. In the same time, the semimajor axis is subject to a slight decrease. The eccentricity oscillations do not remain identical though. Actually the peak eccentricity slightly increases but a more striking fact is that the bottom eccentricity of the Kozai cycles gradually increases to finally reach the peak eccentricity after ~10 Gyr. At this point (beginning of the second phase) the eccentricity oscillations stop, which means that the GJ 436b gets out of the Kozai resonance. The semimajor axis drops then much more quickly and the eccentricity decreases to zero within a few extra Gyr. During almost all the evolution (apart from the beginning), the rotation of GJ 436b remains synchronized with the angular velocity at periastron in high eccentricity phases. The increase of Ω/n during the first phase is related to the gradual increase of the peak eccentricity. In the second phase Ω/n gradually decreases towards 1. At the end of the simulation, GJ 436b is tidally synchronized and circularized, but the final semimajor axis is several times smaller than the initial one. The Kozai migration is ended.
This twofold behavior in Kozai migration process was already noted by Wu & Murray (2003) and Fabrycky & Tremaine (2007). Actually our results of Fig. 4 are very similar to the corresponding ones of Fabrycky & Tremaine (2007).
4.2.2. The first phase
During the first phase, the planet still undergoes Kozai oscillations, but the bottom eccentricity of the cycles gradually increases to finally reach the peak one. This can be explained by a shift in Kozai diagram.
In a pure 3body Newtonian dynamics, the planet just secularly moves along one of the level curves in a Kozai diagram similar to that of Fig. 1. In the ωcirculating case, the curves are explored from left to right with dω/dt > 0, and in the ωlibrating case (the genuine “resonant” orbits), the curves are explored clockwize around the stable point. Now, both GR and tides add an extra, positive component to dω/dt (Fabrycky & Tremaine 2007) at peak eccentricity in the cycles. This results in a slight shift towards right in the Kozai diagram each time the planet passes through a peak eccentricity phase. The planet jumps to another nearby curve located to the right. It is then obvious to see from Fig. 1 that this change induces a decrease of the bottom eccentricity in the ωlibrating case and an increase in the ωcirculating case. Now it is also obvious that thanks to these shifts in diagram at each eccentricity peak, a planet initially in the ωlibrating case will be pushed to move to ωcirculating after a while. Consequently the bottom eccentricity of the planet is expected to eventually first decrease during the ωlibrating phase and then to further increase in the subsequent ωcirculating regime. This is exactly what is reported by Fabrycky & Tremaine (2007) in their Fig. 1. After a short decrease, the bottom eccentricity of the cycles increases gradually to 1. In our Fig. 4, we see that the bottom eccentricity (black curve) of the cycles does not immediately start to increase, but we do not see an initial decrease like in Fabrycky & Tremaine (2007). In fact, the more realistic integration is actually ours. The reason is that the integration of Fabrycky & Tremaine (2007) was made on the basis of a quadrupolar expansion of the Hamiltonian, as well as the Kozai diagram in Fig. 1. It is thus not surprising to see them in perfect agreement. But our integration, made using a symplectic integrator, is not truncated and is accurate to any order in a/a′. Ford et al. (2000) showed that retaining higher order terms, in particular up to the octupolar level (3rd order in expansion in a/a′) considerably improves the accuracy of the integration, especially in the vicinity of the separatrix between the two regimes, which is the case discussed here.
In any case, when the planet enters the ωcirculating regime, the bottom eccentricity of the cycles is forced to increase.
4.2.3. The second phase and the link with the Batygin et al. (2009) model
When the bottom eccentricity reaches the peak one, then the planet moves along a more or less straight line in the Kozai diagram. Tides and GR are now active permanently as the eccentricity is permanently high. This causes the precession of ω to increase. This forces the planet to start a more sharp tidal decrease that makes it get out of the Kozai resonance. In this second phase, the planet gets gradually tidally circularized.
An important outcome of this scenario is that it can provide a solution to the eccentricity paradox of GJ 436b. We note in Fig. 4 that the first phase of the secular evolution lasts several Gyr. Hence GJ 436b may have first spent a long time in the first phase before starting a sharp tidal decrease. Moreover, we note that in the second phase, even if the semimajor axis shrinks quickly, it still takes a few more Gyr to the eccentricity to decrease downto zero. This could appear surprising, as we are now away from Kozai resonance. In Fig. 4, the situation of the planet at ~13 Gyr is similar to the one we studied in Fig. 2. In the latter case, the eccentricity of GJ 436b decreases downto zero in ~2 × 10^{8} yr, but in Fig. 4 it still takes a much longer time from 13 Gyr to reach zero. How could we explain this discrepancy?
The difference comes from the dynamical configuration of the 3body system. It is in fact closely related to the quasistationary solution described by Batygin et al. (2009). Batygin et al. (2009) studied the coplanar 3body system. They showed that there exist two apsidal fixed points in this system, characterized by de/dt = 0, and ϖ − ϖ′ = 0 or ϖ − ϖ′ = π, where ϖ and ϖ′ are the longitude of periastron of GJ 436b and the perturber respectively. There is always one stable point among these two configurations around which librations can occur. In that case, ϖ − ϖ′ just librates around its equilibrium value, and de/dt = 0 periodically vanishes. As a result, the decrease of the eccentricity is much slower. Batygin et al. (2009) showed that adding tides to this peculiar system considerably increases the tidal circularization time, because the angular momentum transfer associated with tidal dissipation is now shared among the two interacting orbits. This is of course a peculiar situation, but as mentioned above, these points tend to behave as attractors.
The situation we describe here differs from that of Batygin et al. (2009) because we do not have a coplanar configuration. But the starting point of the second phase is also characterized by a stationary point with de/dt = 0. At the end of the first phase (Kozai cycles), when the bottom eccentricity of the cycles equals the peak one, GJ 436b moves along a straight line in the Kozai diagram. Hence we have de/dt = 0. Contrary to Batygin et al. (2009), here ϖ − ϖ′ circulates during the second phase. The reason is that in the coplanar problem, there is no quadruplar component to de/dt. The expression of de/dt given by Batygin et al. (2009) arises indeed from the octupolar expansion. But at high inclination, the quadrupolar component to de/dt dominates, so that the constraints are different. But the common feature is that the starting point of both evolutions is characterized by a stationary point in eccentricity in the pure secular 3body evolution. The result is a huge increase of the tidal circularization time by a factor 50 at least, like in the evolution described by Batygin et al. (2009). Besides, the apsidal alignment criterion could be used as a diagnostic. If an apsidally aligned companion is discovered in the future in the radial velocity residuals of the star, it would be a strong indication in favor of a coplanar configuration like described by Batygin et al. (2009); conversely if this is not the case, the companion could be considered as highly inclined with good probability.
5. Exploration of the parameter space
5.1. Description of the results
We claim that the Kozai migration process we describe could fit the present day configuration of the GJ 436 system. GJ 436b could fairly well be today in the middle of the second phase, with an already reduced semimajor axis and a slowly decreasing eccentricity. In Fig. 4, the situation of the planet at ~13 Gyr can be compared to the present day orbital configuration of GJ 436b.
Of course 13 Gyr is too long compared to the age of the GJ 436 system. But the timescale of the twofold evolution can vary considerably depending on the initial parameters of the model, so that adjusting them to fit to any given age less than the age of the Universe is possible.
We come now to exploring the parameter space of this scenario to try to derive in which case it can be compatible with the observation. We tested various orbital configurations, exploring the parameter space. The important parameters are

the mass m′ and semimajor axis a′ (or equivalently the orbital period T′) of the perturber;

the initial semimajor axis a_{0} of GJ 436b;

the initial mutual inclination i_{0} between the two orbits.
Note that instead of i_{0}, the actual critical parameter is the reduced vertical component of the specific angular momentum , which is a secular constant of motion (see Sect. 2.1). But in all our simulations we took an initial eccentricity e_{0} = 0.01 for GJ 436b, so that h is unambiguously related to i_{0}.
Other parameters such as the eccentricity of the perturber, the initial rotation speed and rotation axis orientation of GJ 436b turned out to have only a minor influence on the global behavior of the system, so that exploring the space of the major parameters is enough for our purpose. In particular, in all the following, the initial eccentricity e′ of the perturber was set to 0.03 as a typical standard value.
The key outcome to monitor is the transition time between the two phases in the scenario (hereafter t_{tr}), which occurs close to 10 Gyr in Fig. 4. This timescale appeared to be extremely sensitive to the parameters, ranging from less than 10^{5} yr in some extreme cases to hundreds of Gyr in other configurations. In all cases, most of the orbital decrease and circularization that occurs after t_{tr} is achieved with ~1.5t_{tr}.
In our exploration, we consider as likely any t_{tr} value between 1 Gyr and 10 Gyr. Shorter values are not compatible with a not yet circularized GJ 436b today, and with longer values, the star is not old enough to have entered the second phase yet. Our goal was to derive which sets of parameters yield a convenient t_{r}. To do this, we first used our symplectic code in a few cases, and then our averaged integrator to save computing time. As mentioned above, the averaging process of the Hamiltonian due to the perturber cannot be done in closed form. One has to first expand it in ascending powers of a/a′, to truncate it to some given order and then average over both orbits. Fabrycky & Tremaine (2007) truncate to second order (quadrupolar), while Ford et al. (2000) recommends third order (octupolar). We checked that being in correct agreement with the symplectic integration (in particular in the monitoring of t_{tr}, which is a very sensitive parameter) requires to retain terms up to 4th order in a/a′, i.e. one level beyond octupolar. The use of a truncated integrator does not change the global twofold behavior reported above, but it affects its timescale. This is illustrated in Fig. 4a (eccentricity plot) where we have superimposed to our symplectic integration (black curve) the same eccentricity evolution, but calculated with averaged integrators truncated at various orders: second (blue), third (green) and 4th (red). The use of the quadrupolar truncation generally leads to an error on t_{tr} of ~50%. Most of the time it overestimates it. This can be understood easily: in a pure Newtonian dynamics (no tides, no relativity, such as in Fig. 2a), all Kozai cycles are strictly identical in the quadrupolar approximation. The motion is strictly periodic, as can be seen from Fig. 1. This is due to the absence of the perturber’s eccentricity and argument of periastron in the interaction Hamiltonian (Eq. (2)). This is no longer true at higher order, as shown by Krymolowski & Mazeh (1999); Ford et al. (2000). Due to secular changes of the perturber’s orbit, the Kozai cycles are no longer identical (see Fig. 3a for instance). The total number of Kozai cycles achieved within a given (long) timespan can thus vary depending on the theory assumed. As the efficiency of the Kozai migration process depends on the integrated number of peak eccentricity phases reached, this explains the variability of t_{tr}. We see from Fig. 4a that the use of a third order (octupolar) theory (green curve) considerably enhances the agreement with the symplectic integration. A ~5% agreement is reached if we go to 4th order (red curve). Or course better agreement could be reached if we used even higher order formulas. But in that case, the averaged formulas are so complex that the gain in computing time with respect to the symplectic integration is small. As a 5% accuracy on t_{tr} is enough for our analysis, we decided to use 4th order theory.
Fig. 5 Result of the exploration of the parameter space for four (a_{0},i_{0}) combinations. a) a_{0} = 0.287 AU, i_{0} = 75°; b) a_{0} = 0.1435 AU, i_{0} = 75°; c) a_{0} = 0.287 AU, i_{0} = 85°; d) a_{0} = 0.574 AU, i_{0} = 85°. In each plot, the upper observational limits (radial velocities and imaging) are represented as a solid line; the oblique dotted line corresponds to the limit t_{GR} = t_{Kozai} as given by Matsumura et al. (2009). In the upper left plot, the dashed lined are m′ ∝ T′ − 1/3 and m′ ∝ T′^{2} power laws (see text). 
Figure 5 shows a summary of the result of this exploration of the parameter space. Each plot corresponds to the use of a fixed set of values (a_{0},i_{0}) where we let the other parameters (T′,m′) vary. According to the simulation results, we divide the (T′,m′) plane into different regions depending on the resulting value of t_{tr}. Regions marked “Fast Kozai” correspond to cases where t_{tr} < 1 Gyr. In these regions, the Kozai migration process should have ended a long time ago, and we would expect GJ 436b to be synchronized and circularized today. This situation logically applies to massive and/or closein perturbers (topleft parts of the panels), when the disturbing action is stronger and the period of Kozai cycles smaller. Conversely, towards the bottom parts of the panels (small and/or distant perturbers), we have regions marked as “Slow Kozai”, thus corresponding to t_{tr} > 10 Gyr. In such cases, t_{tr} is just too long to allow GJ 436b to lie in the second phase of the Kozai migration process today. This cannot match the observation, as during the first phase of the evolution (Kozai cycles), the semimajor axis only undergoes a very small decrease an cannot reach the present day value. We also note in some cases regions marked “No Kozai” corresponding to very low mass and distant companions. In such cases, the perturbing action of the companion is so weak that GR (even at the considered distance) is strong enough to prevent the onset of Kozai cycles. Thus GJ 436b’s eccentricity never grows and it can never reach its presentday location. To illustrate this fact, we add to all plots in Fig. 5 an oblique line, which corresponds in each case to t_{GR} = t_{Kozai}, where t_{GR} and t_{Kozai} are characteristic times for GR and for Kozai oscillations respectively. These timescales read (Matsumura et al. 2009) Here we have computed these timescales assuming e = 0. Note that the real period of the Kozai cycles equals t_{Kozai} within a factor of order unity (Ford et al. 2000; Beust & Dutrey 2006). It is straightforward to see that t_{GR} = t_{Kozai} implies m′ ∝ T^{′2} (with m′ ≪ M), hence the oblique dotted lines in Fig. 5. No Kozai oscillations are to be expected below this limit (t_{GR} < t_{Kozai}), which is confirmed by our numerical exploration. The “No Kozai” regions extend indeed somewhat above the dotted lines the our plots. This is due to the higher order effects that are not taken into account in the above simple formulas.
In the middle of the fast and slow Kozai regions, a convenient strip appears where the value of t_{tr} allows GJ 436b to lie presently in the middle of the second phase of the Kozai migration process. Note that for a given set of parameters, the probability of witnessing this phase today is not small, as the second phase lasts typically ~0.5t_{tr}.
To all plots in Fig. 5 we superimpose a solid line that corresponds to the observational upper limits of Montagnier et al. (in prep.) for any additional planet to GJ 436. The left oblique part corresponds to the radial velocity limit and the horizontal upper part to the imaging (adaptive optics) limit. Consequently, any set of (T′,m′) values located above this line must be rejected, and for an initial set of parameters (a_{0},i_{0}) to be possible, the convenient strip must cross the region below the observational limits. In all plots, the lower T′ value (in xaxis) was set according to the stability limit of the system. Systems with smaller T′ values are unstable. The stability was tested for these orbits with the unaveraged symplectic integrator. The upper limit in T′ was set at the point where the convenient strip definitely gets above the observational imaging limit, so that any convenient companion located further away would have been already detected. So, for any given set of initial conditions (a_{0},i_{0}), the convenient values for T′ and m′ must appear on the plot.
For all displayed cases, the convenient strip appears twofold. Towards lower T′ values (closein perturbers), the convenient m′ values are a slowly decreasing function of T′, while for more distant perturbers, they increase rapidly. This can be understood with a careful analysis of the Kozai mechanism. The efficiency of the Kozai migration scenario depends actually on how close GJ 436b can get to the star during the first phase (hence the maximum eccentricity it can reach) and the time it spends at high eccentricity to allow tides to be at work. For closein perturbers, the process is limited by the maximum eccentricity it can reach. In fact in Fig. 5, closein perturbers in the convenient strip have masses comparable to that of GJ 436b (≃ 0.074 M_{Jup}), and have semimajor axes only a few times more than GJ 436b initially. Hence the initial angular momentum ratio between the perturber and GJ 436b is of order unity or a bit more. Now, due to total angular momentum conservation (we neglect exchanges with spins), we necessarily have (15)where l′ and l stand for the angular momenta of both bodies at peak eccentricity, and i for the mutual inclination in the same phase. Assume now (the perturber is not affected) and (the semimajor axis is only little affected in the first phase). This equation simplifies to (16)Solving this equation for e shows that for moderate values, the eccentricity cannot go beyond a maximum value that is significantly below 1, while for larger value, any value as close possible to 1 is allowed. In other words, the angular momentum of GJ 436b cannot reach zero because that of the perturber cannot grow to account for the whole angular momentum of the system. Let us assume for instance i_{0} = 75°, i = 30° (typical values) and . We derive e ≃ 0.74, which is high, but significantly below 1 to avoid GJ 436b to be efficiently tidally braked in peak eccentricity phases. The efficiency of the Kozai migration process is therefore limited by the maximum eccentricity value than can be reached, hence by the value of . Similar t_{tr} correspond to similar values. Assuming now and l_{0} ∝ mT^{1/3} (low initial eccentricities), a constant t_{tr} leads m′ ∝ T^{′−1/3}, hence the decrease of the convenient strip towards closein perturbers.
For more distant perturbers, the ratio is large enough not to limit the peak eccentricity of the Kozai cycles. The efficiency of the process is now limited by the amount of time spent in high eccentricity phases, as tides are only active there. We thus expect t_{tr} ∝ t_{Kozai}, as the longer t_{tr}, the less time is spent within a given time span in high eccentricity phases. As t_{Kozai} ∝ (M/m′) × (T^{′2}/T), a fixed value for t_{tr} will correspond to a constant ratio T^{′2}/m′. This explains the increase of the convenient strip towards distant perturbers. A more distant perturber just needs to be more massive to generate the same Kozai migration as a closer one.
To illustrate this simple analysis, we have drawn in Fig. 5a two dashed lines corresponding to arbitrary m′ ∝ T^{′ − 1/3} and m′ ∝ T^{′2} power laws. We note that these power laws described above are well followed in the simulations for both regimes, as the convenient strip quite faithfully remains parallel to the dashed lines sketched. To summarize, the convenient strip roughly follows a dual power law regime, with m′ ∝ T^{′ − 1/3} for small T′ and m′ ∝ T^{′2} for large T′. Its exact location in (T′,m′) space depends on the initial set of parameters (a_{0},i_{0}), as can be seen from Fig. 5.
We also note that in all cases, the convenient strip appears quite narrow, although this effect is enhanced in Fig. 5 by the logarithmic scale. In fact, for a given period T′, Kozai migration is very sensitive to the mass m′ of the perturber, so that m′ must assume a well defined value within a factor ~1.5 to allow t_{tr} to fall in the right range. However, for a wide range of initial values sets (a_{0},i_{0}) (see discussion below), there is a convenient strip in (T′,m′), which is compatible with the observational constraints, so that the probability of the scenario is not small.
5.2. Discussion
We come now to compare the location of the convenient strip between various configurations. If we compare Figs. 5a and 5c (same a_{0} but different i_{0}), we note that for i_{0} = 85° the convenient strip is located lower than for i_{0} = 75°. Here again this may be easily understood. With a larger i_{0} the Kozai mechanism is stronger, so that GJ 436b more easily gets a high peak eccentricity and the tides act more efficiently. Hence a less massive perturber is required to achieve the same result.
If we compare configurations with same i_{0} but different a_{0} (Figs. 5a and b on the one hand, and c and d on the other hand) we can see how the convenient strip moves in the (T′,m′) diagram if we vary a_{0}. We see that increasing a_{0} brings the convenient strip up in the decreasing part of the diagram (m′ ∝ T^{′ − 1/3}), and down in the increasing part of the diagram (m′ ∝ T^{′2}). If for instance we superimpose Figs. 5a and b, the convenient strip for a_{0} = 0.1435 AU falls below that for a_{0} = 0.287 AU for small T′, then crosses it and get above for large T′. This can be explained. In the decreasing part, the efficiency of the Kozai mechanism is controlled by the angular momentum ratio . for a larger a_{0}, l_{0} is larger and the angular momentum ratio is smaller, so the Kozai mechanism is less efficient. We thus need a larger (and hence a larger m′) to get the same results. This is why the convenient strip is brought up for a larger a_{0}. In the increasing part of the diagram, the Kozai mechanism is controlled by t_{Kozai}. As we know that t_{Kozai} ∝ T^{1}, a larger a_{0} leads naturally to a smaller t_{Kozai} and to a more efficient Kozai migration process. Hence a less massive companion (larger m′) for the same T′ is required to achieve the same t_{Kozai}. This is why the convenient strip is brought down if we increase a_{0}.
This regime tends nevertheless not to be valid for all values of a_{0}. For larger initial orbits, the whole convenient strip is shifted upwards if we increase a_{0}. For instance, for a_{0} = 0.574 AU and i_{0} = 75° (configuration not displayed in Fig. 5), the convenient strip is located much higher than with a_{0} = 0.287 AU, so that is it barely compatible with the observational constraints. Similarly, with i_{0} = 85°, the convenient strip is shifted upwards for a_{0} ≳ 0.6 AU. This behavior is due to the fact that the efficiency of the tides are related to the minimum periastron reached rather than the maximum eccentricity. As in the first phase, the semimajor axis of GJ 436b is only little affected, for a larger a_{0} a higher eccentricity is necessary to reach the required periastron range.
Now, is it possible to set absolute limits to the unknown parameters? Setting a lower limit to i_{0} is easy. As described above, if we crease i_{0} the convenient strip gets higher in (T′,m′) space. At some point it crosses the observational limits so that the considered (a_{0},i_{0}) must be rejected. We tested several other runs with i_{0} = 65° and i_{0} = 70°. In all cases the convenient strip (if any) appears above the observational limits. It is then possible to stress that i_{0} must be at least 70–75°. Note also that the Kozai mechanism is symmetrical with respect to i_{0} = 90°. We checked that runs assuming i_{0} = 95° and i_{0} = 105° behave exactly the same way as corresponding runs with i_{0} = 85° and i_{0} = 75° respectively. It is thus possible to put an upper limit to i_{0} to 105–110°. The initial configuration of the system must have been close to perpendicular within less than 20°. This is of course valid for small initial eccentricity e_{0}. If e_{0} is high, then i_{0} may take different values. In the general case, a more relevant constraint must be derived for . 70° ≲ i_{0} ≲ 90° for small e_{0} means h ≲ 0.34, which is valid irrespective of e_{0}.
It is also possible to set a lower limit to a_{0}. Looking at Fig. 5, we see that for lower a_{0} the convenient strip gets shifted towards left, so that at some point it should cross the radial velocity observational limit. This occurs for a_{0} ≃ 0.1 AU. Hence we can stress that the initial value a_{0} must have been larger.
Giving an upper limit to a_{0} is less easy. We first can notice than choosing a larger a_{0} implies a much more restricted compatible range for i_{0}. As noted above, for i_{0} = 75°, a_{0} = 0.287 AU is possible, but a_{0} = 0.574 AU is barely excluded (too high convenient strip). For i_{0} = 85° nevertheless, we tested that even fairly large a_{0} values (up ~10 AU) we can find a convenient strip in (T′,m′) space that falls still below the observational limits. We may thus set an absolute upper limit to a_{0} ≃ 10 AU. However, we must point out that for such values, we are in an extreme Kozai regime with a very high peak eccentricity. Even if such a situation cannot be rejected, whether the onset of such a regime with nearly orthogonal orbits is likely is questionable. Hence even if the absolute compatible range for a_{0} is 0.1–10 AU, we may stress that the range 0.1–0.5 AU is more probable.
As pointed out above, in all cases, the convenient strip in (T′,m′) space to allow t_{tr} to fall in the suitable time range is quite narrow. This does not tell that the situation we are describing is improbable, as suitable configurations are possible for a wide range of initial configurations. This tells that for any given set of values (a_{0},i_{0},T′) the mass of the companion m′ must fall around a well defined value within a factor ~1.5. We nevertheless note that taking larger a_{0} values leads to generate narrower convenient strips. This shows up already in Fig. 5d (a_{0} = 0.574 AU), and it is even more striking for larger a_{0}. This is another reason why we think that a configuration with a moderate a_{0} is more likely.
Our parametric study also shows that the evolution we describe is not the only one possible. First the initial inclination i_{0} must be high to have a Kozai migration. Second, if the perturber’s mass m′ is too small, there is no Kozai migration, and if it is too high, Kozai migration is extremely rapid. Kozai migration was also suspected to account for the misaligned systems that have been measured from RossiterMcLaughlin effect. But among all the systems for which an obliquity with respect to the stellar axis has been measured, only ~30% show an obliquity larger than 45° (Hébrard et al. 2011). This shows that Kozai migration is not a universal process. The limited size of the convenient strip in our parametric study supports this idea. Systems like GJ 436 with eccentric closein planets, although not unrealistic, are not expected to be numerous.
6. Conclusion
The Kozai migration scenario is a generic process that causes inward migration of planets in non coplanar planetary systems. This evolution is characterized by two phases, a first one with perturbed Kozai cycles, and a second one with a more drastic decay of the semimajor axis. In any case, the characteristic time for this evolution can be considerably longer than the standard tidal circularization time, providing this way a solution to the issue of the high present eccentricity of GJ 436b. If we combine our study and that of Montagnier et al. (in prep.), various suitable initial configurations can be found. They all require a high initial mutual inclination (≳ 75°), an initial orbit for GJ 436b several times wider than today, perturber masses ranging between less than 0.1 M_{Jup} and 50 M_{Jup}, with orbital periods ranging between a few years up to several hundreds. This model is also generic enough to be invoked to explain some other puzzling cases of high eccentricity closein exoplanets. But for a given set of initial conditions, the orbital period and the mass of the perturber must assume a fairly well constrained relationship, otherwize the suspected process is either too rapid or nonexistent. This could then explain why system like GJ 436 are not numerous.
Note that the constraint derived for i_{0} is valid for small initial eccentricity e_{0}. The real condition is h ≲ 0.3, which can be achieved either with low e_{0} and high i_{0}, but also with low i_{0} and high e_{0}. For small i_{0}, h ≲ 0.3 implies e_{0} ≳ 0.95. Although this is actually just a matter of fixing the starting point on the Kozai cycles, this has a different meaning for what concerns the formation of such a planetary system. Planetary systems are usually thought to form in circumstellar disk with initial coplanar configurations (see e.g. theoretical studies by Pierens & Nelson 2008; Crida et al. 2009). The only way to imagine an initially inclined configuration is to assume that both planets formed independently and were assembled further to build the GJ 436 system. This could happen for instance if GJ 436 was initially member of a multiple stellar system. Such systems are indeed rarely coplanar. Alternatively, the system could have formed coplanar but the inner planet could have been driven to high eccentricity. There are many ways to achieve this. This can be done by planetplanet scattering events or chaotic diffusion. Laskar (1996) indeed showed (with application to Mercury, see also Wu & Lithwick 2011) that in any planetary system, the inner planet may be subject to drastic eccentricity increases due to secular chaos. Similarly, recent work on the socalled Nice model for the formation of the solar system (Tsiganis et al. 2005; Levison et al. 2011) have demonstrated the role of planetplanet scattering among the giant planets. Meanmotion resonances can also drive inner bodies to high eccentricities. This mechanism was invoked to explain the suspected falling evaporated bodies (FEB) scenario in the βPictoris system (Beust & Morbidelli 1996, 2000).
In conclusion, an initially nearly coplanar but eccentric configuration of the system appears more realistic than a highly inclined one. Note that this does not change anything to the further evolution described in this paper, as Kozai cycles cause the system to rapidly oscillate between these two extreme configurations. Even if the system is initially nearly coplanar, it evolves to highly inclined within half a Kozai cycle. The relevant issue is actually the definition of the starting point. According to Batygin & Morbidelli (2011), in all systems potentially subject to Kozai resonance, Kozai cycles are first frozen as long as self gravity and strong planetplanet interactions induce a large enough apsidal precession. But at some point when the system relaxes, Kozai cycles can start, which constitutes the “starting point” in our scenario.
The meanmotion resonance model that applies to the FEB scenario can also provide a valuable starting point. The only requirement is that both planets are initially in a meanmotion resonance like 4:1 ou 3:1, with a small initial mutual inclination (a few degrees) and a small eccentricity (≳ 0.05) for the outer, more massive planet. If this is fulfilled, the inner body is subject to an eccentricity increase virtually up to ~1. It was shown by Beust & Valiron (2007) that whenever it reaches the high eccentricity state (the FEB state in the case of βPictoris), it undergoes high amplitude inclination oscillations that can be viewed as Kozai cycles inside the meanmotion resonance. This is a kind of resonant version of our model that could constitute a full coherent scenario driving the system from a coplanar and circular configuration to strong Kozai resonance.
The reality of this scenario will nevertheless be difficult to test observationally. Even if the two planets are initially in meanmotion resonance, as soon the semimajor axis of GJ 436b starts to evolve thanks to tidal dissipation, they leave the resonance. The present day orbital configuration keeps no memory of the initial resonance. Assuming that we are today in the middle of the second phase with a still significant eccentricity and an already reduced semimajor axis (by a typical factor 10), we are now unable to distinguish between initially resonant and nonresonant configurations. We also note from Fig. 4 that the present day mutual inclination between the two planets should be in any case ~20°, which makes it difficult also to distinguish observationally between coplanar (Batygin et al. 2009) and Kozai models (this work). However, as explained above, the exact characterization of the apsidal behavior of the planets can lead to constraints on their mutual inclination. Basically, if ϖ − ϖ′ is close to 0 or π, this would constitute a strong clue in favor of the coplanar model of Batygin et al. (2009). Otherwize, this would indicate an inclined configuration.
Before that, a crucial issue is to continue the radial velocity and photometric follow up of GJ 436. If we confirm the presence of a more distant companion, and if we are able to significantly constrain its orbit, then we will be able to test the proposed scenarii into more details. This would not only resolve the eccentricity paradox of GJ 436b, but it would also provide clues towards the system’s past evolutionary history.
Acknowledgments
All computations presented in this paper were performed at the Service Commun de Calcul Intensif de l’Observatoire de Grenoble (SCCI) on the supercomputer funded by Agence Nationale pour la Recherche under contracts ANR07BLAN0221, ANR2010JCJC050401 and ANR2010JCJC050101. We gratefully acknowledge financial support from the French Programme National de Planétologie (PNP) of CNRS/INSU. We thank our referee, Konstantin Batygin, for valuable suggestions to improve the discussion of the paper.
References
 Alonso, R., Barbieri, M., Rabus, M., et al. 2008, A&A, 487, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bailey, M. E., Chambers, J. E., & Hahn, G. 1992, A&A, 257, 315 [NASA ADS] [Google Scholar]
 Ballard, S., Christiansen, J. L., Charbonneau, D., et al. 2010, ApJ, 716, 1047 [NASA ADS] [CrossRef] [Google Scholar]
 Banfield, D., & Murray, N. 1992, Icarus, 99, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Barker, A. J., & Ogilvie, G. I. 2009, MNRAS, 395, 2268 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Laughlin, G., Meschiari, S., et al. 2009, ApJ, 699, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2011, Celest. Mech. 111, 219 [Google Scholar]
 Beust, H., Corporon, P., Siess, L., Forestini, M., & Lagrange, A.M. 1997, A&A, 320, 478 [Google Scholar]
 Beust, H. 2003, A&A, 400, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beust, H., & Dutrey, A. 2006, A&A, 446, 137 [Google Scholar]
 Beust, H., & Morbidelli, A. 1996, Icarus, 120, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Beust, H., & Morbidelli, A. 2000, Icarus, 143, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Beust, H., & Valiron, P. 2007, A&A, 466, 201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Butler, R. P., Vogt, S. S., & Marcy, G. W. 2004, ApJ, 617, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ćuk, M., & Burns, J. A. 2004, Icarus, 167, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B. O., Gillon, M., Barman, T. et al., 2007, A&A, 475, 1125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deming, D., Harrington, J., Laughlin, G., et al. 2007, ApJ, 667, L199 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P., Kiseleva, L. G., & Hut, P. W. 1998, ApJ, 499, 853 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P., & KiselevaEggleton, L. 2001, ApJ, 562, 1012 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Pont, F., Demory, B.O., et al. 2007, A&A, 472, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Harrington, R. S. 1968, AJ, 73, 190 [NASA ADS] [CrossRef] [Google Scholar]
 Hébrard, G., Ehrenreich, D., Bouchy, F., et al. 2011, A&A, 527, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Holman, M., Touma, J., & Tremaine, S. 1997, Nature, 386, 254 [NASA ADS] [CrossRef] [Google Scholar]
 Innanen, K. A., Zheng, J. Q., Mikkola, S., & Valtonen, M. J. 1997, AJ, 113, 1915 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396 [NASA ADS] [CrossRef] [Google Scholar]
 Kehoe, T. J. J., Murray, C. D., & Porco, C. C. 2003, AJ, 126, 3108 [NASA ADS] [CrossRef] [Google Scholar]
 Kinoshita, H., & Nakai, H. 1999, Celest. Mech., 75, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Laskar, J. 1996, Celest. Mech., 64, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorný, D., & Gomes, R. 2011, AJ, 142, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2009, A&A, 493, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malhotra, R. 1994, Celest. Mech., 60, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Malmberg, D., Davies, M. B., & Chambers, J. E. 2007, MNRAS, 377, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Mardling, R. A. 2008 [arXiv:0805.1928] [Google Scholar]
 Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., & Weidenschilling, S. J. 2006, BAAS, 38, 505 [NASA ADS] [Google Scholar]
 Maness, H. L., Marcy, G. W., Ford, E. B., et al. 2007, PASP, 119, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., Takeda, G., & Rasio, F. A. 2009, in Transiting Planets, eds. F. Pont, D. Sasselov, & M. Holman, Proc. IAU Symp. 253, 189 [Google Scholar]
 Pierens, A., & Nelson, R. P. 2008, A&A, 482, 333 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pont, F., Gilliland, R. L., Knutson, H., Holman, M., & Charbonneau, D. 2009, MNRAS, 393, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Ribas, I., FontRibera, A., & Beaulieu, J.P. 2008, ApJ, 677, L59 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Saha, P., & Tremaine, S. 1992, AJ, 104, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Söderhjelm, S. 1982, A&A, 107, 54 [NASA ADS] [Google Scholar]
 Tong, Xiao, Zhou, JiLin 2009, Sci. Chin. Ser. G, 52, 4 [CrossRef] [Google Scholar]
 Torres, G. 2007, ApJ, 671, L65 [NASA ADS] [CrossRef] [Google Scholar]
 Trilling, D. E. 2000, ApJ, 537, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wisdom, J. 2008, Icarus, 193, 637 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, K., & Hamilton, D. P. 2008, Icarus, 193, 267 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Contour plots of the Kozai Hamiltonian (4) in a (ω,e) space for H = 0.05. 

In the text 
Fig. 2 Secular evolution of the eccentricity of GJ 436b starting from the presentday orbit and assuming a standard perturber (m′ = 0.1 M_{Jup}, T′ = 10 yr) in different conditions (see text): a) perturbing planet, but no tides nor relativity; b) tides and relativity but no perturbing planet; c) perturbing planet, relativity but no tides; d) perturbing planet, relativity and tides. 

In the text 
Fig. 3 Secular evolution over 2 × 10^{7}yr of GJ 436b assuming an initial value a = 0.287 AU and a perturber with m′ = 0.1 M_{Jup} and a′ = 3.583 AU (T′ = 10 yr): eccentricity (top left), semimajor axis (top right), rotation speed with respect to the orbital meanmotion (Ω/n, bottom left), and obliquity angle (bottom right). On the eccentricity plot we overplot in red the result of the same integration, but using a conventional Bulirsh & Stoer integrator. 

In the text 
Fig. 4 Same evolution as in Fig. 3, but over 2 × 10^{10}yr: eccentricity (top left), semimajor axis (top right), rotation speed with respect to the orbital meanmotion (Ω/n, bottom left), mutual inclination between the two orbits (bottom right). On all plots, the black curve corresponds to our symplectic integration. The color curves on the eccentricity plot are calculated with a fully averaged integrator, with interaction Hamiltonian truncated up to various orders: 2nd order (blue), 3rd order (green), and 4th order (red). 

In the text 
Fig. 5 Result of the exploration of the parameter space for four (a_{0},i_{0}) combinations. a) a_{0} = 0.287 AU, i_{0} = 75°; b) a_{0} = 0.1435 AU, i_{0} = 75°; c) a_{0} = 0.287 AU, i_{0} = 85°; d) a_{0} = 0.574 AU, i_{0} = 85°. In each plot, the upper observational limits (radial velocities and imaging) are represented as a solid line; the oblique dotted line corresponds to the limit t_{GR} = t_{Kozai} as given by Matsumura et al. (2009). In the upper left plot, the dashed lined are m′ ∝ T′ − 1/3 and m′ ∝ T′^{2} power laws (see text). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.