DYNAMICAL EVOLUTION AND SPIN–ORBIT RESONANCES OF POTENTIALLY HABITABLE EXOPLANETS: THE CASE OF GJ 581d

The habitability of an emerging class of super-Earths—exoplanets with masses greater than that of the Earth but lower than that of Uranus—depends on a combination of physical parameters. Among these, the intensity of irradiation from the host star is crucial. A favorable rate of irradiation permits water to be available in liquid form. The irradiation intensity depends on the luminosity of the star, the size of the planet's orbit, and, to a lesser degree, on the orbital eccentricity. The chemical composition of the planet's atmosphere influences the average temperature and the temperature variations on the surface. For example, estimates show that a certain minimal amount of CO₂ in the atmosphere of GJ 581d is required to keep water above the freezing point on the surface, while GJ 581c is likely to have experienced a runaway greenhouse event (von Paris et al. 2010; Hu & Ding 2011). Planet GJ 581d, which is the main target of our study, may be located on the outer edge of the habitable zone, according to von Braun et al. (2011).

We begin our study by addressing the still-controversial problem of the composition of this planetary system. In Section 2, we confirm that only the four originally detected planets (b through e) are real and detectable at the 0.99 confidence level by our fully automated, fast grid search algorithm. The two additional planets, f and g, proposed by Vogt et al. (2010) are not found with any combination of the published radial velocity (RV) data.

The orbit estimation technique employed in our paper is designed, in particular, to produce accurate and robust results for the eccentricity of each detected planet—an important asset for the subsequent dynamical simulations. The dynamical stability and the presence of chaos in the orbital motion of the four planets is investigated in Section 3 by long-term integrations, assuming zero inclination and coplanar orbits. The system is found to be stable for the long term but strongly chaotic, with the eccentricity and semimajor axis of planet d varying little over gigayears. This allows us to investigate, in Section 4, the spin–orbit dynamics of GJ 581d with its orbital parameters fixed. The planet is assumed to have a terrestrial rheology and a near-zero obliquity.

Although this question has attracted a lot of attention in the literature, some doubts seem to be lingering, and a few papers have been published on subtle features of the planets that may well be nonexistent. The first four planets were discovered with the same instrument, HARPS at La Silla Observatory: GJ 581b (Bonfils et al. 2005), GJ 581c and d (Udry et al. 2007), and GJ 581e (Mayor et al. 2009). The combined series of observations with HARPS, published in the latter paper, spanned 1570 days, while a typical single measurement precision was about 1 m s⁻¹.

At the first stage of our investigation, we processed the HARPS data with a simple benchmark planet detection algorithm. This algorithm is an iterative process of optimization, based on a simple grid search in the space of free parameters. The grid search is applied only to the fitting of the nonlinear parameters, i.e., of the orbital period, eccentricity, and periastron time. The other parameters are determined by a direct least-squares solution. At the first stage of fitting, the most significant sinusoidal variation in the RV data is determined by a corrected Lomb–Scargle periodogram method (Lomb 1976; Scargle 1981). Our upgrades to the method mainly concern the way the common offset of the RV points is treated. Instead of the traditional subtraction of the mean RV from the data prior to the periodogram computation, we fit the entire set of model functions [1, cos (2πt_j/p_i), sin (2πt_j/p_i)] for each trial period p_i, where t_j are the times of observations. This seemingly trivial modification allows us to update the amplitudes of already detected planets once a new planet is detected, and to keep the original RV measurements unchanged throughout the data reduction cycle. The mean RV is inevitably biased when a non-integer number of waves is present in the data. Using this method, subtraction of the mean RV from the data is only legitimate when no periodic signals are present in the data, which defeats the purpose of periodogram analysis. It can be proven that the periodogram value in the classical Lomb's algorithm is equivalent to the absolute change of the χ² before and after the least-squares fit of the cos - and sin -terms. However, fitting these terms to the data corrupted by the subtraction of the mean RV gives rise to side lobes and biases in the periodogram. The downside of this modification is that a direct computation of the post-fit χ² by Scargle's formula is no longer possible. Instead, the complete least-squares solution of the matrix equation should be computed for each trial period. This complication, however, becomes negligible, given present-day computing power.

Once the period of the most prominent sinusoidal variation is determined, a separate grid search is implemented to determine the best-fitting eccentricity. This technique is based on the theory of harmonic decomposition of orbital motion (e.g., Brouwer & Clemence 1961; Konacki et al. 2002). The observed RV variation due to the orbital motion of a single planet can be written as

$$\begin{eqnarray} \dot{z}&=&{-}C\frac{2\pi }{P}\sum _{l=1}^\infty l\,\alpha _l(e)\cos (l{\cal {M}})\nonumber\\ && +H\frac{2\pi }{P}\sqrt{1-e^2}\sum _{l=1}^\infty l\,\phi _l(e)\sin (l{\cal {M}}), \end{eqnarray} \tag{ 1 }$$

where e is the eccentricity, C and H are the Thielle–Innes constants, P is the orbital period determined from the periodogram, ${\cal {M}}$ is the mean anomaly, and α_l(e) and ϕ_l(e) are coefficients that only depend on e and define the relative magnitudes of harmonics. Although the coefficients can be computed via the Bessel functions of the first kind, we find it practical to compute them numerically using Kepler's equation at the required resolution in e between zero and one, and to store the resulting values of the coefficients in a table. Further processing is then reduced to solving the overdetermined equations

$$\begin{equation} \left[ 1 \sum _lp_l\sin (l{\cal {M}}_j) \sum _lq_l\cos (l{\cal {M}}_j) \right]{\bf s}=\dot{z}_j, \end{equation} \tag{ 2 }$$

where p_l = lα_l(e) and $q_l=l\phi _l(e)\sqrt{1-e^2}$ are tabulated harmonic coefficients, and the integer j = 1, 2, ..., N serves to number the available data points. A least-squares solution of the problem is performed for a grid of e, and results in the selection of the value minimizing the post-fit χ². The solution vector $\bf s$ includes three elements: s(1) is the RV offset, while s(2) and s(3) are the fitting coefficients from which the Thielle–Innes constants C and H can be derived. Parallel to this, a separate grid search is carried out for the phase of orbital motion, i.e., for the mean anomaly at the first observation.

Finally, we have to estimate the confidence of the orbital solution. We are using an adaptation of the significance F-test, which was carefully tested and verified by Monte Carlo simulations on a large number of real and artificial planetary systems. In this case, the null hypothesis for the F-test is that the data set contains a constant offset and a random, uncorrelated noise. Fitting of a harmonic [cos (2πt_j/p_i), sin (2πt_j/p_i)] to this set by the least-squares method will generate a χ² change corresponding to a random realization of the F-statistic, F_i, which is distributed as F_PDF. The confidence level for a detection at p_i is simply the probability of F ⩽ F_i, which is defined by the known cumulative distribution function F_CDF. It is customary to accept a detection if the confidence is greater than 0.99, i.e., if the observed reduction in χ² can happen within the null hypothesis in less than 1 trial out of 100. The confidence of detection is computed as

$$\begin{equation} P_{\rm conf}=1-\left[1-F_{\rm CDF}\left(\frac{\chi _1^2-\chi _2^2}{d_2} \frac{d_1}{\chi _2^2}\right)\right]f_{\rm pow}, \end{equation} \tag{ 3 }$$

where χ²₁ and χ²₂ are the values of χ² before and after the orbital fit, respectively; meanwhile, d₁ and d₂ are the corresponding numbers of degrees of freedom. The f_pow multiplier is required to take into account the fact that we are testing multiple periods p_i and are selecting the one that delivers the greatest reduction in χ². If the probability of a single realization of F to be less than a certain limit F_lim is F_CDF(F_lim), then the probability of the largest F among n independent trials to be less than F_lim is F_CDF(F_lim)ⁿ. Thus, the f_pow multiplier in Equation (3) is the number of independent frequencies for a given periodogram search window. For irregularly spaced observation times and search intervals that often exceed the Nyquist boundaries, evaluation of this number is nontrivial. We employ the following estimate, which is analogous to the Nyquist limit: f_pow = ceil (ΔT/2/s_mean), with s_mean being the mean separation between the RV data points.

Handling of multiple planet detections is another important feature of the algorithm. Once a planet with period P_m is detected, the corresponding terms [∑_lp_lcos (2π(t_j − t_0m)/P_m), ∑_lq_lsin (2π(t_j − t_0m)/P_m)] are permanently added to the model. Then, a search for another planet is performed outside the already detected frequencies, on the original RV data (which are never altered in this algorithm). With m planets already detected, the size of the model is 2m + 1. Also, the number of degrees of freedom d₁ is reduced by 5 with each detected planet. The coefficients C and H and the RV offset are re-adjusted each time a new planet is detected, which provides for a better decoupling of planets in complex systems, especially of planets on commensurable or resonant orbits. Figure 1 shows the initial χ² periodogram of the HARPS data prior to any planet detection. The dips corresponding to the planets detected from these data are marked with arrows and planet designations. The strongest signal is associated with planet b, which shows up first. The planets are detected in their historical order: b, c, d, and e. The next strongest signal in the cleaned periodogram corresponds to a period of 391 d. Formally, though, it should be rejected because the confidence level is only 0.975.

Zoom In Zoom Out Reset image size

The orbital and physical parameters estimated from the fits of the four detected planets are given in Table 1. Our periods are in close agreement with the estimates by Mayor et al. (2009), but e and ω (the longitude of the periastron) are markedly different. The likely reason is that these parameters for the planets b and e were "fixed" at 0 in Mayor et al. (2009), which is not advisable for multiple systems. In particular, ignoring the detectable eccentricity in a least-squares adjustment retains the periodogram side lobes and harmonics associated with the previously detected planets and can distort the results for other planets. Our estimate of e for planet d is 0.27, which is lower than the 0.38 ± 0.09 obtained in Mayor et al. (2009). Our results confirm that the planets c and e may be in the 4:1 orbital resonance, and the closeness of our estimated ω lends more credence to this possibility.

Table 1. Orbital Parameters of the Four-planet GJ 581 System (b–e)

Planet	P	Mass	a	e	ω	${\cal {M}}_0$
	(days)	(M☉)	(AU)		(°)	(°)
... b	5.37	4.8 × 10−5	0.041	0.01	288.2	52
... c	12.91	1.8 × 10−5	0.073	0.09	349.7	160
... d	66.98	1.8 × 10−5	0.218	0.27	180.9	92
... e	3.15	5.5 × 10−6	0.028	0.13	327.2	356

Download table as:  ASCIITypeset image

Our results generally agree with the work by Tuomi (2011), who re-analyzed the HARPS and HIRES data together using a Bayesian approach and found evidence for the existence of only four planets. However, their estimated eccentricities are consistent with zero for all four planets. We too applied our benchmark algorithm to the combination of HARPS and HIRES data, introducing separate RV offsets for the two data sets. This resulted in the detection of only two planets, b and c, but with significantly weaker signals. The planet d was rejected due to insufficient confidence. Finally, the algorithm completely failed to reproduce the results obtained by Vogt et al. (2010) from their data. In their most recent update on the GJ 581 system, Vogt et al. (2012) used the yet unpublished, much expanded set of RV data from the HARPS (Forveille et al. 2011). They found that a fifth planet (f) can be detected if the eccentricity of the 66.7 day planet (d) is set to zero. Using the same data set, but not rejecting any data points from it, we were able to confirm this statement with our detection algorithm by forcing the eccentricity of the four planets b–e to zero. A 32.1 day planet with a projected mass of 1.99M_E emerges at a confidence level above 99%. It is likely therefore that we are dealing with the second harmonic in the signal of an eccentric planet, which can be confused with a separate, nearly commensurate planet. A potentially successful method for resolving this controversy is to look for the third harmonic of the 67 day signal, which should be generated by an eccentric planet with that period. Indeed, we find a dip in the χ² periodogram at p_j ≈ 22.5 days after fitting out the four planets at zero eccentricity, but the presence of this dip cannot be taken as conclusive evidence. At e = 0.27, the amplitudes of the first three harmonics scale as 1:0.256:0.073, hence, the expected amplitude of the third harmonic is only 2 × 0.073 ≃ 0.15 m s⁻¹, and so it should drown in the observational noise. Vogt et al. (2012) make a strong argument that a system of four eccentric planets can be dynamically unstable. We therefore set out to establish that the original four-planet system as determined from the HARPS measurements is dynamically viable.

The preceding studies of the dynamical status and evolution of the planetary system GJ 581 have been mainly concerned with verification of its physical stability. Beust et al. (2008) integrated the system of three planets (b, c, and d) for 10⁸ yr at different inclinations to the line of sight (which changes the estimated planetary masses by a factor of sin ⁻¹i) and for different orbital eccentricities. Mayor et al. (2009) performed similar integrations for the four-planet system (b–e) with an updated period of planet d. In all of these simulations, the orbits were assumed to be coplanar. The systems inferred from the RV data proved stable for inclinations down to i ≃ 30°. At smaller inclinations, the masses of the super-Earth planets would become so large that the inner planet GJ 581e would have been ejected quickly. Beust et al. (2008) also computed the maximum Lyapunov exponents and concluded that planet d is less chaotic than the inner planets. Within a wide range of initial parameters, the system is chaotic but long-term stable, which indicates that the regions of bounded chaos (Laskar 1990) are extensive.

Using the values listed in Table 1 as initial parameters, and assuming the orbits to be coplanar, we performed multiple simulations of the dynamical evolution of the five-body GJ 581 system. The symplectic integrator HNBody, version 1.0.7, (Rauch & Hamilton 2002) was utilized with the symplectic option. We chose a time step of 5 × 10⁻⁵ yr (0.018 days) and integrated the trajectories of all five bodies over 100 million years. As previous studies have found, the system appears to be quite stable. The eccentricities and semimajor axes show small periodic variations over the integration time. For GJ 581d, these variations over the first 10,000 yr of integration are presented in Figure 2. The eccentricity of this planet oscillates, seemingly forever, within a narrow range around 0.27. The semimajor axis, too, varies very little.

Zoom In Zoom Out Reset image size

In order to quantify the chaos in GJ 581, we employed the sibling simulation technique developed by Hayes (2007, 2008) to investigate the behavior of the outer solar system with slightly different initial conditions. Two sibling trajectories are generated by perturbing the initial semimajor axis of the planet GJ 581b by a factor of 10⁻¹⁴. The distance between the unperturbed planets and their siblings is then computed as a function of time. Exponential divergence between the trajectories indicates that the system is chaotic. We find that chaos sets in within 5000 yr, with a characteristic Lyapunov time of only ∼30 yr. The distance between siblings is shown for planet GJ 581d in Figure 3. In this case, the epoch of exponential divergence is close to 4000 yr, with a characteristic Lyapunov time of only 38 yr. The dominating planet b, on the other hand, displays two epochs of exponential divergence, with Lyapunov times of 18 and 277 yr.

Zoom In Zoom Out Reset image size

Long-term dynamical stability coupled with strong but bounded chaos is not uncommon among the currently known systems of multiple exoplanets. For example, a similar chaotic behavior was found for the dynamically robust system 55 Cnc, where a hidden commensurability may exist for the dominating 14-day planet. This makes the solar system, whose chaos is much slower, to stand out. In the context of our study, the important conclusion is that the parameters of GJ 581d listed in Table 2 can be safely taken as being constant over extended intervals of time.

Table 2. Default Parameters of GJ 581 and Planet d

Parameter	Value
ξ	... $\frac{2}{5}$
R	... 1.7 REarth
Mplanet	... 7.1 MEarth
Mstar	... 0.31 M☉
a	... 0.218 AU
e	... 0.27
(B − A)/C	... 5 × 10−5
Porb	... 67 days
τM	... 50 yr
μ	... 0.8 × 1011 kg m−1 s−2
α	... 0.2

Download table as:  ASCIITypeset image

A planet is subject to the torque $\stackrel{\rightarrow }{{\cal {T}}}^{\, \rm {_{(TRI)}}}$ due to the planet's permanent triaxiality, and to the torque $\stackrel{\rightarrow }{{\cal {T}}}^{\, \rm {_{(TIDE)}}}$ generated by tidal deformation. In our model, we shall assume that both these torques are exerted upon the planet by its host star only. We shall thus neglect the torques exerted upon the planet by its moons or by other planets.

4.1. The Equation of Motion

Consider a planet with a mean radius R and mass M_planet and treat it as the primary. The tide-raising perturber (the star of mass M_star) will be regarded as the secondary effectively orbiting the primary.

The principal moments of inertia of the planet will be denoted as A, B, and C, assuming that A < B < C. The maximal moment of inertia will read as C = ξM_planetR², with the numerical coefficient ξ assuming the value of 2/5 in the homogeneous-sphere limit. These and other notations used in our study are listed in Table 3.

Table 3. Symbol Key

Notation	Description
ξ	... Moment of inertia coefficient of the planet
R	... Radius of the planet
${\cal {T}}^{\rm {{(TRI)}}}$	... Triaxiality-caused torque acting on the planet
${\cal {T}}^{\rm {{(TIDE)}}}$	... Tidal torque acting on the planet
Mplanet	... Mass of the planet
Mstar	... Mass of the star
a	... Semimajor axis of the planet
r	... Instantaneous distance between the planet and the star
f	... True anomaly of the planet
e	... Orbital eccentricity
${\cal {M}}$	... Mean anomaly of the planet
C	... The maximal moment of inertia of the planet
B	... The middle moment of inertia of the planet
A	... The minimal moment of inertia of the planet
n	... Mean motion, i.e., $\sqrt{G(M_{{\rm planet}}+M_{{\rm star}})/a^{3}}$
G	... Gravitational constant, =66468 m3 kg−1 yr−2
τM	... Viscoelastic characteristic time (Maxwell time)
τA	... Anelastic characteristic time
μ	... Unrelaxed rigidity modulus
J	... Unrelaxed compliance
α	... The Andrade parameter

Download table as:  ASCIITypeset image

Suppose that the planet is rotating about its major-inertia axis z, the one corresponding to the maximal moment of inertia C. Let θ be the sidereal angle of this rotation, reckoned from an arbitrary line fixed in inertial space (say, the line of apsides) to the axis x of the largest elongation, i.e., the axis corresponding to the minimal moment of inertia A of the planet.

Rotation will then be described by the equation

$$\begin{equation} \ddot{\theta }= \frac{{\cal {T}}^{\rm {{(TRI)}}}_z+{\cal {T}}^{\rm {{(TIDE)}}}_z}{C}=\frac{{\cal {T}}^{\rm {{(TRI)}}}_z+ {\cal {T}}^{\rm {{(TIDE)}}}_z}{\xi M_{{\rm planet}}R^2}, \end{equation} \tag{ 4 }$$

with the subscript z serving to denote the polar components of the torques. High-accuracy integration of Equation (4) requires a step much smaller than the orbital period.

4.2. The Triaxiality-caused Torque

We approximate the polar component of the torque with its quadrupole part given by¹

$$\begin{eqnarray} {\cal {T}}^{\rm {{(TRI)}}}_z&=&\frac{3}{2} (B-A) \frac{{\it GM}_{{\rm star}}}{r^3} \sin 2\psi \end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray} &&\nonumber \\ &\approx &{-} \frac{3}{2} (B-A) n^2 \frac{a^3}{r^3} \sin 2(\theta -f). \end{eqnarray} \tag{ 5b }$$

Here, G stands for the Newton gravitational constant, r denotes the distance between the centers of mass of the two bodies, and M_star stands for the mass of the star (which, in our setting, is playing the role of perturbing secondary effectively orbiting the tidally deformed primary). The mean motion is given by $n\equiv \sqrt{G(M_{{\rm planet}}+M_{{\rm star}})/a^{3}}\approx \sqrt{GM_{{\rm star}}/a^3}$, with the understanding that the star is much more massive than the planet.

Our notations for angles are given in Figure 4.

Zoom In Zoom Out Reset image size

The angle ψ displays the separation between the planetocentric direction toward the star and the principal axis x of the maximal elongation (the minimal-inertia axis). It is equal to the difference between the angles f and θ made by these two directions onto an arbitrary line fixed in inertial space. If this fiducial line is chosen to be parallel to the line of apsides connecting the star with the empty focus, then f will be the planet's true anomaly, as seen from the star, while θ will be the sidereal angle (which we agreed to reckon from the line of apsides.)

Employing (5b) in practical computations, it is advantageous to approximate the functions (a/r)³sin 2f and (a/r)³cos 2f with truncated series of the Hansen's coefficients (e.g., Branham 1990), which in turn can be approximated with series of the Bessel functions of the first kind.

4.3. The Tidal Torque

Due to a several-orders-of-magnitude difference in the intensities of internal friction between a terrestrial planet and a star, we shall take into account only the tides exerted by the star on the planet, and shall neglect the tides exerted on the star.

The technique of calculation of tidal torques had until recently remained in a somewhat embryonic state, requiring extensive corrections. On the one hand, it had long been a common habit to combine the expression for the tidal torque with unphysical rheological models. As a result, much of the hitherto obtained results on the spin history of terrestrial planets (e.g., Celletti et al. 2007; Heller et al. 2011) have to be reexamined, because they were based on very ad hoc rheologies incompatible with the behavior of realistic solids. On the other hand, in many publications, the description of the tidal torque was marred by a mathematical mistake which proliferated through many papers and ended up in textbooks. Explanation and correction of that long-standing oversight is presented in Efroimsky & Makarov (2012).

We shall limit our treatment to the simpler case, where the planet is not too close to the star (R/a ≪ 1), the obliquity of the planet is small (i ≃ 0), and its eccentricity e is not very large. As demonstrated in Appendix A, under these assumptions the polar component of the secular part of the tidal torque can be approximated by

$$\begin{eqnarray} \big\langle {\cal {T}}_z^{\rm {{(TIDE)}}}\big\rangle _{ {{ { {{l=2}}}}}} &=& \frac{3}{2} GM_{{\rm star}}^{2}R^5a^{-6}\sum _{q=-1}^{7} G^{2}_{ {{ {{20\mbox{\it {q}}}}}}}(e) k_2(\omega _{ {{ {{220\mbox{\it {q}}}}}}})\nonumber\\ &&\times \sin |\epsilon _2(\omega _{ {{ {{220\mbox{\it {q}}}}}}})| \,{\rm Sgn}\,(\omega _{220q}) +O(e^8\epsilon)+O({\it i}^2\epsilon), \nonumber\\ \end{eqnarray} \tag{ 6 }$$

where R denotes the radius of the planet, $\dot{\theta }$ stands for its rotation rate, and a, e, and i are the semimajor axis, the eccentricity, and the obliquity of the perturber (the star) as seen from the planet. The notation G_lpq(e) stands for the so-called eccentricity polynomials that coincide with the Hansen polynomials X^{(− l − 1), (l − 2p)}_{(l − 2p + q)}(e).

This expression first appeared, without proof, in the work by Goldreich & Peale (1966), who summed over all integer values of q. For the reasons explained in Appendix A below, we limit the summation to the range from q = −1 through q = 7. This truncation effectively limits the Taylor expansions of eccentricity-dependent functions to terms of the order of up to e⁷, inclusive. The resulting relative precision of our calculations is approximately 0.1%, because the largest and the smallest terms of the so-truncated sum contain G₂₀₀(0.27) = 0.82202 and G₂₀₇(0.27) = 0.01392, respectively; meanwhile, all of the eccentricity polynomials G_20q(e) outside the range of summation (i.e., for q outside of the range −1, ..., 7) assume much smaller values.

In Equation (6), both the dynamical Love number k₂ and the phase lag ₂ are functions of the tidal Fourier mode, which in this simplified case reads as

$$\begin{eqnarray} &&\omega _{ {{220\mbox{\it {q}}}}}=(2+q)n-2\dot{\theta }, \end{eqnarray} \tag{ 7 }$$

with a more general expression for ω_lmpq given in Appendix A. As usual, n denotes the mean motion.

The dynamical Love numbers k₂(ω_220 q) are positive definite, while the sign of each lag ₂(ω_220 q) coincides with that of the tidal mode ω_220 q, as explained in Appendix A. Therefore, each factor k₂(ω_220 q)sin ₂(ω_220 q) can be written as k₂(ω_220 q)sin |₂(ω_220 q)| Sgn (ω_220 q).

Through Equation (7), each of this type of factor becomes a function of the planetary spin rate $\dot{\theta }$:

$$\begin{eqnarray} &&\nonumber k_2(\omega _{ {{220\mbox{\it {q}}}}})\sin |\epsilon _2(\omega _{ {{220\mbox{\it {q}}}}})| \,{\rm Sgn}\,(\omega _{ {{220\mbox{\it {q}}}}})= k_2(2(n-\dot{\theta })+qn)\nonumber\\ &&\qquad\times \sin |\epsilon _2(2(n-\dot{\theta })+qn)| \,{\rm Sgn}\,(2(n-\dot{\theta }) +qn). \end{eqnarray} \tag{ 8 }$$

Accordingly, the entire sum (6) can be regarded as a function of $\dot{\theta }$. The mean motion and eccentricity act as parameters, with a much slower evolution than that of $\dot{\theta }$.

In each term of series (6), the factor k₂sin ₂, expressed as a function of $\dot{\theta }$, has the shape of a kink. In Figure 5, the dotted line² illustrates the behavior of the factor

$$\begin{eqnarray} &&k_2(\omega _{ {{ {{2202}}}}}) \sin \epsilon _2(\omega _{ {{ {{2202}}}}})= k_2(4n-2\dot{\theta }) \sin \epsilon _2(4n-2\dot{\theta })\nonumber\\ &&\quad= k_2(4n-2\dot{\theta }) \sin |\epsilon _2(4n-2\dot{\theta })| \,{\rm Sgn}\,(4n-2\dot{\theta }). \end{eqnarray} \tag{ 9 }$$

The origin of the kink shape is explained in Appendix B, with references provided. From the physical viewpoint, the emergence of this shape is natural, in that each term should transcend zero and change its sign smoothly when the spin rate goes through the appropriate spin–orbit resonance.

Zoom In Zoom Out Reset image size

However, in Equation (6), the kink-shaped factors are accompanied by different multipliers G²_20q(e). So the sum (6), as a function of $\dot{\theta }$, will be a superposition of many kinks with different magnitudes, centered at different resonances (9 kinks, if the sum over q goes from −1 to 7). The resulting curve will cross the horizontal axis in points close to the resonances $\dot{\theta }=n\left(1+{ q}/{ 2}\right)$, but not exactly in these resonances—see the solid line in Figure 5.

Goldreich & Peale (1968) were the first to point out that in the vicinity of a particular resonance q', i.e., when $\dot{\theta }$ is close to n(1 + q'/2), the right-hand side of Equation (6) can be conveniently broken down into two parts. One part is the q = q' term, an odd function vanishing at $\dot{ \theta }=n(1+{ q^{\prime }}/{ 2})$. As was demonstrated later by Efroimsky (2012a, 2012b), for realistic rheologies, this function has the shape of a kink—see the dotted line in Figure 5. Another part, called the bias, is comprised of the rest of the sum. So the bias is the contribution of all the q ≠ q' modes to the values assumed by the torque in the vicinity of the q = q' resonance. For eccentricities that are not too large, the bias is usually negative in value. The bias is a very slowly changing function, which can, to a good approximation, be treated as constant.

The q = q' term by itself is an odd function, and it goes through nil at exactly $\dot{\theta }=n(1+{ q^{\prime }}/{ 2})$. However, the bias displaces the location of zeroes. In Figure 5, the torque (depicted with a solid line) is mostly defined by the term with q = 2 (rendered by the dotted line). However, the curve is shifted down due to the bias which is defined mainly by the right slope of the q = 1 kink located to the left. As the right slope of the q = 1 kink is negative, the q = 2 kink in Figure 5 is shifted slightly down, and the zero is located slightly to the left of the point $\dot{\theta }=2n$. The crossing point in this case is at $\dot{\theta }/n=1.999976481$. Such a miniscule shift of the equilibrium away from the resonance frequency does not have any practical consequences, because the net nonzero tidal torque is compensated for by a tiny secular triaxial torque, as discussed in more detail in Makarov & Efroimsky (2012).

The shifts of the tidal-equilibrium frequencies at resonances are larger for bodies with a lower Maxwell time, e.g., for bodies whose mantles contain a large fraction of partial melt. Still, no matter how shifted the tidal equilibrium happens to be, the mean rotation is exactly resonant due to the presence of the compensating triaxiality-caused secular torque.

When the planet's spin rate approaches a resonance $2\dot{\theta }=(2+q)n$, with an integer q, the planet can be captured, or can traverse the resonance, depending on the specific trajectory in the phase space. A method for estimating capture probabilities was developed by Goldreich & Peale (1966), for two simplified models of the tidal torque. A clearer explanation of this theory was presented in a later work by Goldreich & Peale (1968), to which we shall refer.

In one model considered in these papers, the torque was assumed to be linear in the tidal frequency. To be more precise, in the Fourier expansion of the torque over tidal modes, each torque component was set to be proportional to the appropriate tidal mode³—see formula (19) in Goldreich & Peale (1968).

Another model addressed in the referenced above two works model was the constant-torque. Specifically, in the Fourier expansion of the torque over tidal modes, each torque component was set to be a constant multiplied by the sign of the corresponding mode—see Equation (29) in Goldreich & Peale (1968).

In their treatment of both models, Goldreich & Peale (1968) only took into account the secular, orbit-averaged, component of the tidal torque, and ignored the existence of an oscillating component. Below, we shall test the validity of this approximation.

The pivotal ideas and formulae of the capture theory are explained in short in Appendix C below. Here, we shall employ some of those formulae, though with an important difference. Instead of the toy models introduced in Goldreich & Peale (1968), we shall rely on a realistic rheology of solids. As we shall see, capture probabilities are sensitive to changes in (at least some of the) rheological parameters. For example, the probabilities turn out to be more sensitive to the mantle's Maxwell time than to the planet's triaxiality.

The principal result of the analysis carried out by Goldreich & Peale (1968) is their estimate of the probability of capture into an arbitrary resonance q = q', i.e., into a steady rotation at the rate of $\dot{\theta }=n(1+{ q^{\prime }}/ { 2})$. The probability is given by

$$\begin{equation} P_{\rm capt}=\frac{ 2}{ {1+{ 2 \pi V}/{ \int _{-\pi }^{\pi }W(\dot{\gamma })d\gamma }}}, \end{equation} \tag{ 10 }$$

where the new, auxiliary variable γ is defined through

$$\begin{eqnarray} &&\gamma \equiv 2 \theta - (2+q^{\prime }){\cal {M}}, \end{eqnarray} \tag{ 11 }$$

its time derivative thus being the negative double of the tidal mode $\omega _{ {{220\mbox{\it {q}}^{\prime }}}}$:

$$\begin{eqnarray} &&\dot{\gamma }=- (2+q^{\prime })n+2\dot{\theta }=-\omega _{220\mbox{\it {q}}^{\prime }}. \end{eqnarray} \tag{ 12 }$$

This mode vanishes in the q = q' resonance, so we may say that this resonance corresponds to $\dot{\gamma }=0$. The odd function $W(\dot{\gamma })$ included in expression (10) is called kink and is given by

$$\begin{eqnarray} &&W(\dot{\gamma }) = - K G^{2}_{ {{220\mbox{\it {q}}^{\prime }}}} k_2(\dot{\gamma }) \sin |\epsilon _2(\dot{\gamma })| \, {\rm Sgn}\, (\dot{\gamma }), \end{eqnarray} \tag{ 13 }$$

where K is a positive constant. Remember that our definition of γ is in agreement with that used by Makarov (2012) and is twice that introduced in Goldreich & Peale (1966, 1968).

Also included in formula (10) is a function⁴

$$\begin{eqnarray} \nonumber V&=&K\sum _{q\ne q^{\prime }}G^{2}_{220\mbox{\it {q}}} k_2(\omega _{ {{220\mbox{\it {q}}}}}) \sin | \epsilon _2(\omega _{ {{220\mbox{\it {q}}}}})| \, {\rm Sgn}\, (\omega _{ {{220\mbox{\it {q}}}}})\\ &&\nonumber \\ &=&K\sum _{q\ne q^{\prime }}G^{2}_{220\mbox{\it {q}}} k_2((q-q^{\prime })n) \sin | \epsilon _2((q-q^{\prime })n)| \, {\rm Sgn}\, (q-q^{\prime })\nonumber\\ \end{eqnarray} \tag{ 14 }$$

called bias. The physical meaning of V and $W(\dot{\gamma })$ is explained in Appendix C. The Appendix also explains how in the above integral, the dependence of $\dot{\gamma }$ upon γ should be set so that the integral could be evaluated.

It is important that the above expressions for V and $W(\dot{\gamma })$ contain the mode-dependent factor k₂(ω_220 q)sin ₂(ω_220 q). The functional form of its mode-dependence is defined by the self-gravitation and rheology of the planet (Efroimsky 2012a). This observation helps us to mark the point at which our analysis will diverge from that carried out by Goldreich & Peale (1968): our choice of a realistic rheology will yield for k₂sin ₂ a mode-dependence different from both cases addressed in Goldreich & Peale (1968).

Recently, Makarov (2012) generalized Goldreich & Peale's (1968) treatment for a more physical tidal torque, i.e., for a torque whose terms bear a mode-dependence stemming from the properties of realistic solids. The rheology of silicates and ices is well described by the Andrade model, insofar as the tidal frequency is sufficiently high and the inner friction is dominated by anelasticity. At lower frequencies, the friction is dominated by viscosity, and the behavior of the material becomes closer to that of the Maxwell body (Efroimsky 2012a, 2012b). The resulting mode-dependencies of the terms of the torque are presented below in Appendix B. These mode-dependencies are more complicated than the two models considered by Goldreich & Peale (1968).

In our computations, we ignore the oscillatory components of the tidal torque, and consider only the secular part which can be obtained by averaging over one orbital period—a convenient simplification is to be justified below. Thus, we begin in the spirit of Goldreich & Peale (1968), i.e., we use formula (10), with the functions W and V provided by formulae (13) and (14), correspondingly. Following Makarov (2012), we then insert into those formulae the realistic mode-dependence of k₂sin ₂, written down and explained in Appendix B. Finally, we perform a brute force numerical check of whether the neglect of the oscillating part is legitimate.

The resulting capture probabilities as functions of eccentricity are presented in Figure 6 for two values of the Maxwell time: τ_M = 50 yr (the left plot) and τ_M = 500 yr (the right plot). While 500 yr is the current Maxwell time of Earth's mantle, the value of 50 yr is deemed to represent a slightly⁵ warmer planet, one with a lower viscosity.

Zoom In Zoom Out Reset image size

From the plots, we see that probabilities of capture are sensitive to the Maxwell time. Warmer planets with lower τ_M have more chances of being captured at spin–orbit resonances, as they loose their spin angular momentum. In particular, for τ_M = 50 yr and (B − A)/C = 5 × 10⁻⁵, the capture probabilities are 1 in resonance 3:2, 0.897 in 2:1, 0.383 in 5:2, 0.133 in 3:1, and 0.040 in 7:2. Therefore, if GJ 581d has had enough time to de-spin into a state of spin–orbit equilibrium at the current value of eccentricity, then the probabilities of the planet now being in resonance are the following: 0.05 in resonance 3:2, 0.46 in 2:1, 0.32 in 5:2, 0.13 in 3:1, and 0.04 in 7:2.

For τ_M = 500 yr and (B − A)/C = 5 × 10⁻⁵, the capture probabilities are 1 at resonance 3:2, 0.568 in 2:1, 0.188 in 5:2, 0.054 in 3:1, and 0.014 in 7:2. Then, with a probability of 0.35, the planet is currently entrapped in resonance 3:2, 0.43 in 2:1, 0.18 in 5:2, 0.05 in 3:1, and 0.01 in 7:2. The 2:1 resonance is the most likely state of the planet in both cases. However, for τ_M = 50 yr, the second likeliest state is 5:2, while for τ_M = 500 yr it is 3:2.

The probabilities of capture depend on the degree of triaxiality through the parameter (B − A)/C in the equation of separatrix (C6). This equation tells us that larger values of (B − A)/C entail greater libration amplitudes of $\dot{\gamma }$. The kink $W(\dot{\gamma })$, which is the antisymmetric part of the secular torque, becomes smaller with $\dot{\gamma }$ growing outside the resonance, see Figure 5. Therefore, the integral $\int _{-\pi }^{\,\pi } W(\dot{\gamma })d\gamma$ in Equation (C11) is expected to become smaller for larger (B − A)/C. Accordingly, the capture probability is to become smaller for larger (B − A)/C. This is confirmed by computations of capture probabilities for τ_M = 50 yr, (B − A)/C = 2 × 10⁻⁴, and the other parameters as given in Table 2. The capture probabilities are 1 in the resonance 3:2, 0.776 in 2:1, 0.301 in 5:2, 0.097 in 3:1, and 0.028 in 7:2. Each of these probabilities is slightly smaller than its counterpart taken for the smaller⁶ triaxiality (B − A)/C = 5 × 10⁻⁵. At the same time, the relatively small differences between the capture probabilities into the same resonance, for different values of (B − A)/C, imply that the spin–orbit resonances are more sensitive to τ_M than to (B − A)/C.

As a way of spot-checking these theoretical results, we conducted brute-force simulations of the GJ 581d system. The equation of motion, incorporating both the tidal and triaxial torques acting on the planet, was numerically integrated 40 times for τ_M = 50 yr, (B − A)/C = 5 × 10⁻⁵, e = 0.27, the initial spin rate $\dot{\theta }(0)=2.51n$, the initial mean anomaly ${\cal {M}}(0)=0$, and the initial sidereal angle θ(0) = πi/40, with i = 0, 1, ..., 39. This method of estimation tacitly assumes that the initial sidereal angle at a fixed rate of rotation can assume any values with equal probability (which is less than obvious). With this assumption accepted, simulations spanning 7000 yr, with a step of 1.5 × 10⁻³ yr, resulted in 14 captures into the 5:2 resonance and 26 passages. The estimated probability of capture is thus 14/40 = 0.35, which is surprisingly close to the theoretical estimate 0.383. To make the numerical and theoretical estimations consistent, the former included only the secular part of the torque (6). This points to one of the potential weaknesses of the semi-analytical derivation of probabilities. The oscillatory terms of the tidal torque are ignored altogether.⁷ Furthermore, even though the deformity of the planet (B − A)/C enters the computation of capture probabilities, the cyclic variations of the spin rate caused by the triaxial torque are not involved in any way. In reality, however, the smooth sinusoidal separatrix trajectories defined by Equation (C6) are superimposed with a jitter caused by the harmonics of the time-dependent terms of the torque.

Figure 7 illustrates the role of the oscillatory torques that are ignored in the theoretical calculation of capture probability. The graph renders the behavior of the quantity $\dot{\gamma }/2$, which is proportional to the energy of rotation, as a function of γ for two full libration cycles, one going forward toward the point of resonance ($\dot{\gamma }>0$) and the other going back toward $\dot{\gamma }=0$ beyond the point of resonance ($\dot{\gamma }<0$). Due to the dissipation of energy by the tides, the latter curve is systematically lower than the former one, but the difference is so small that it cannot be seen on the graph. Besides, the oscillatory terms of the net torque make the curves jittery. Although the amplitude of the jitter diminishes in the vicinity of the resonance, it may appear to be significant for the outcome of this process. According to Goldreich & Peale (1968), capture occurs when the first post-resonance minimum of $\dot{\gamma }^2/2$ reaches 0. If the jitter is superimposed with the smooth separatrix trajectory, the chances to bump into 0 seem greater than without jitter. The inset in Figure 7 shows in much greater detail this important segment of the post-resonance libration, in the axes $\lg (\dot{\gamma }/2)$ versus γ. The rotational energy comes very close to 0 at the lower extent of oscillations, never quite reaching it. And indeed, in this simulation, the planet traversed the 2:1 resonance.

Zoom In Zoom Out Reset image size

To better understand the influence of oscillating terms of tidal torque on the chances of the planet to be captured, we repeated the 40 high-accuracy integrations described in the previous paragraph, with the same initial parameters, though this time including the entire set of periodic terms. We found that the results for individual simulations often changed, i.e., what resulted in a capture for a given set of parameters became a passage, and vice versa. Surprisingly, we recorded 14 captures out of 40 trials, yielding the same probability of 0.35. From this small-scale experiment, we see that while the outcome of a particular integration may change because of the jitter of the tidal torque, the overall probability of capture is unlikely to depend on it.

The growing number of detected systems of multiple exoplanets and the impressive quality of observational data collected for them allow astronomers to perform analysis of the probable dynamical states and evolution of these remote worlds at a level of detail that was unthinkable just several years ago. Still, this analysis is riddled with difficulties and uncertainties. The planets of GJ 581 present a challenge for both observational practice and interpretational theory. The story of two fictitious planets tells a lesson about the hazards of combining, without proper caution, the data from two instruments with their own sets of systematic errors. It also calls for a certain unification of the planet detection techniques or, at least, for an easily accessible and well-tested standard detection algorithm available as a web application. It is fine to apply a variety of different detection methods to the same data, but the standard algorithm should always be checked, and discrepancies, if any, should be investigated and reported. N-body integration of detected systems should also be standard, reducing the probability of error. The system of GJ 581 proves to be remarkably stable, with the four bona fide planets remaining on their orbits despite the strong evidence of chaos. The characteristic Lyapunov times are very short compared to the dynamical lifetime of the system.

We have explored the rotation history of the planet GJ 581d that is assumed to have composition similar to that of the terrestrial planets of the Solar system.

Contrary to the previous publications on the subject, which were based on ad hoc, simplistic tidal theories, we find that this planet cannot be captured in a synchronous or pseudo-synchronous rotation if it began its evolution from faster, prograde initial spin rates. Instead, for a plausible range of parameters, the most likely state of the planet is the 2:1 spin–orbit resonance. In this state, the day on GJ 581d should be 67 Earth's days long, which bodes for an inhospitable environment, though the potential habitability of this planet cannot be ruled out just on climatic considerations. The case of 2:1 spin–orbit resonance was considered in the simulations of a hypothetical atmosphere on GJ 581d by Wordsworth et al. (2011) whose modeling confirmed the possibility of liquid water being present on the surface, under some favorable conditions.

The next likeliest equilibrium states are the 3:2 or 5:2 resonances, depending on the temperature and viscosity of the mantle (much less on the planet's triaxiality).

At the same time, in the event that the initial rotation of the planet was retrograde, the most probable final state is synchronous rotation.

It is our pleasure to thank Jérémy Leconte for a fruitful discussion on the topic of the paper, and for referring us to the work by Wordsworth et al. (2011).

A laborious calculation (presented in Efroimsky 2012a, 2012b) demonstrates that the factors k_l(ω_lmpq)sin _l(ω_lmpq) can be expressed via the real and imaginary parts of the complex compliance of the mantle material, and the mass and radius of the planet. In this way, the mode-dependence of these factors is defined by both the rheology and self-gravitation of the planet. The factors come out to be odd functions, which is very natural, since each such factor (and the term of the torque, which contains this factor) must change its sign when the appropriate commensurability is transcended.

Being odd, the factors can then be written down as k_l(ω_lmpq)sin |_l(ω_lmpq)| Sgn (ω_lmpq), with the product k_l(ω_lmpq)sin |_l(ω_lmpq)| as an even function of the tidal mode. In other words, this product may be regarded as a function not of the tidal mode ω_lmpq but of its absolute value χ_lmpq = |ω_lmpq|, which is the physical forcing frequency of tidal oscillations in the mantle. All in all, we have

$$\begin{eqnarray} k_l(\omega _{ {{lmpq}}}) \sin \epsilon _l(\omega _{ {{lmpq}}})&=&k_l(\omega _{ {{lmpq}}}) \sin |\epsilon _l (\omega _{ {{lmpq}}})| \, {\rm Sgn}\, (\omega _{ {{lmpq}}})\nonumber \\ &=&k_l(\chi _{ {{lmpq}}}) \sin |\epsilon _l (\chi _{ {{lmpq}}})| \, {\rm Sgn}\, (\omega _{ {{lmpq}}}). \nonumber\\ \end{eqnarray} \tag{ B1 }$$

The calculation in Efroimsky (2012a, 2012b) provides the following frequency dependence for these factors:

$$\begin{eqnarray} &&k_l(\omega _{ {{lmpq}}})\sin \epsilon _l(\omega _{ {{lmpq}}})\nonumber\\ &&\quad=\frac{3}{2({\it l}-1)}\frac{-A_lJ{\cal {I}}{\it {m}}[\bar{J}(\chi)]}{({\cal {R}}{\it {e}}[\bar{J}(\chi)]+A_lJ)^2+({\cal {I}}{\it {m}}[\bar{J}(\chi)])^2} \, {\rm Sgn}\, (\omega _{ {{lmpq}}}), \nonumber\\ \end{eqnarray} \tag{ B2 }$$

where χ is a shortened notation for the frequency χ_lmpq, and the coefficients A_l are given by

$$\begin{eqnarray} &&A_{\it l} \equiv \frac{ {(2{\it {l}}^{2}+4{\it {l}}+3){\mu }}}{ {{\it {l}}\mbox{g} \rho R}}=\frac{ {3(2{\it {l}}^{2}+4{\it {l}}+3){\mu }}}{ {4{\it {l}} \pi G\rho ^2R^2}}, \end{eqnarray} \tag{ B3 }$$

with R, ρ, μ, and g being the radius, mean density, unrelaxed rigidity, and surface gravity of the planet, and G being the Newton gravitational constant.

The functions ${\cal R}{\it e} [\bar{J}(\chi)]$ and ${\cal I}{\it m} [\bar{J}(\chi)]$ are the real and imaginary parts of the complex compliance $\bar{J}(\chi)$ of the mantle. These are rendered by the formulae

$$\begin{eqnarray} &&{\cal R}{\it e} [ \bar{J}(\chi)]=J+J(\chi \tau _{A})^{-\alpha }\cos \left(\frac{\alpha \pi }{2}\right) \Gamma (\alpha +1) \end{eqnarray} \tag{ B4 }$$

and

$$\begin{eqnarray} &&{\cal I}{\it m} [ \bar{J}(\chi)]=-J (\chi \tau _{M})^{-1}-J(\chi \tau _{A})^{-\alpha }\sin \left(\frac{\alpha \pi }{2}\right)\Gamma (\alpha +1), \nonumber\\ \end{eqnarray} \tag{ B5 }$$

where J = 1/μ is the unrelaxed compliance of the mantle, and α is a numerical parameter assuming values of about 0.3 for solid silicates and about 0.14–0.2 for partial melts. In our computations, we used α = 0.2.

Among the rheological parameters entering (B4–B5) is the viscoelastic (Maxwell) time τ_M, which is the ratio of the mantle's viscosity η and rigidity μ. In the present geological epoch, the Maxwell time of the terrestrial mantle is about 500 yr. For warmer mantles, it may be much shorter, given the exponential temperature-dependence of the viscosity.

Another characteristic time entering the above expressions is the anelastic (Andrade) time τ_A. Referring the reader to Efroimsky (2012a, 2012b) for details, we would mention that below some threshold, frequency anelasticity ceases to play a major role in the internal friction, giving way to viscosity. Thus, the mantle's behavior becomes closer to that of a Maxwell body. Mathematically, this means that below the threshold the parameter τ_A increases rapidly as the frequency goes down. So only the first term in (B4) and the first term in (B5) survive, and we are thus left with the complex compliance of a Maxwell material.

In our computations, we treated τ_A in the same way as in Makarov (2012): we kept τ_A = τ_M at the frequencies above the threshold (which was set to be 1 yr⁻¹, just like in the solid Earth case). For frequencies lower than that, we set τ_A to grow exponentially with the decrease in the frequency, so the rheological model approached the Maxwell one in the low-frequency limit. Numerical simulation has demonstrated that the resulting capture probabilities are not very sensitive to how quickly the switch to the Maxwell model takes place. For more details, see Makarov (2012).

In a computer code, it is easier to divide both the numerator and denominator of (B2) by J²:

$$\begin{eqnarray} &&k_l(\omega _{ {{lmpq}}})\sin \epsilon _l(\omega _{ {{lmpq}}})=\frac{3}{2({\it l}-1)}\frac{-A_l{\cal {I}}}{\left({\cal {R}} +A_l\right)^2+{\cal {I}}^{ {{2}}}} \, {\rm Sgn}\, (\omega _{ {{lmpq}}}),\nonumber\\ \end{eqnarray} \tag{ B6 }$$

where ${\cal R}$ and ${\cal I}$ are the dimensionless real and imaginary parts of the complex compliance:

$$\begin{eqnarray} &&{\cal R}=1+(\chi \tau _{A})^{-\alpha }\cos \left(\frac{\alpha \pi }{2}\right) \Gamma (\alpha +1), \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} &&{\cal I}=-(\chi \tau _{M})^{-1}-(\chi \tau _{A})^{-\alpha }\sin \left(\frac{\alpha \pi }{2}\right)\Gamma (\alpha +1). \end{eqnarray} \tag{ B8 }$$

These dependencies were also used in Makarov (2012) to explore the spin–orbit dynamics of a Mercury-like planet.

This section offers a brief explanation of the capture theory developed by Goldreich & Peale (1966, 1968).

As was mentioned in Section 4.3, and explained at length in Appendix A above, a good approximation for the polar component of the tidal torque is provided by expression (6). The tidal-mode-dependent factors entering that expression can be expressed as functions of the spin rate $\dot{\theta }$, see formula (8). Following Goldreich & Peale (1968), in the vicinity of each particular resonance q', i.e., when $\dot{\theta }$ is close to n(1 + q'/2), it is instrumental to break down the right-hand side of Equation (6) into two parts, one being the q = q' kink, another part being the bias. Constituted by the inputs from all the q ≠ q' modes, the bias is smooth and virtually constant in the vicinity of the q' resonance.

The q = q' term is an odd function of $\omega _{ {{220\mbox{\it {q}}^{\prime }}}}$. While Goldreich & Peale (1968) assumed this term to be either a constant multiplied by Sgn$(\omega _{ {{220\mbox{\it {q}}^{\prime }}}})$ (the constant-torque model) or a constant multiplied by the mode $\omega _{ {{220\mbox{\it {q}}^{\prime }}}}$ itself (the linear in frequency model), Makarov (2012) endowed this term (in fact, all terms) with a realistic mode-dependence originating from the properties of actual rocks.

It is convenient to introduce an auxiliary variable

$$\begin{eqnarray} &&\gamma \equiv 2 \theta - (2+q^{\prime }){\cal {M}} \end{eqnarray} \tag{ C1 }$$

that vanishes when the long axis of the planet points toward the star at the perigee. Then, the q = q' resonance will correspond to $\dot{\gamma }=0$ because

$$\begin{eqnarray} &&\dot{\gamma }=- (2+q^{\prime })n+2\dot{\theta }=-\omega _{ {{220\mbox{\it {q}}^{\prime }}}}. \end{eqnarray} \tag{ C2 }$$

Remember that our γ coincides with that employed by Makarov (2012), and is twice the quantity γ introduced in Goldreich & Peale (1968).

In the absence of tidal friction, the equation of motion near the said resonance looks like

$$\begin{eqnarray} &&C \ddot{\gamma } + 3 (B-A) \frac{M_{{\rm star}}}{M_{{\rm star}}+M_{{\rm planet}}} n^2 G_{20\mbox{{q}}^{\prime }}(e) \sin \gamma = 0. \qquad \end{eqnarray} \tag{ C3 }$$

Multiplying this by $\dot{\gamma }$, with a subsequent integration over time t, provides the first integral of motion,

$$\begin{eqnarray} &&C \frac{{\dot{\gamma }}^{2}}{2} - 3 (B-A) \frac{M_{{\rm star}}}{M_{{\rm star}}+M_{{\rm planet}}} n^2 G_{20\mbox{{q}}^{\prime }}(e) \cos \gamma = E^{\prime },\qquad \end{eqnarray} \tag{ C4 }$$

whose value depends on the initial conditions. The latter equation can also as be rewritten

$$\begin{eqnarray} &&C \frac{{\dot{\gamma }}^{2}}{2} - 3 (B-A) \frac{M_{{\rm star}}}{M_{{\rm star}}+M_{{\rm planet}}} n^2 G_{20\mbox{{q}}^{\prime }}(e) 2 \cos ^2\frac{\gamma }{2} = E,\nonumber\\ \end{eqnarray} \tag{ C5 }$$

with E differing from E' by a γ-independent constant. As demonstrated in Goldreich & Peale (1968), the vanishing of the integral of motion E corresponds to the separatrix dividing rotations from librations. So the equation of this separatrix is¹¹

$$\begin{equation} \dot{\gamma }= 2n\left[\frac{3 (B-A)}{C} \frac{M_{{\rm star}}}{M_{{\rm star}}+M_{{\rm planet}}} G_{20\mbox{\it {q}}^{\prime }}(e)\right]^{{1}/{2}}\cos \frac{\gamma }{2}. \end{equation} \tag{ C6 }$$

In terms of $\dot{\gamma }$, the q = q' term of the tidal torque reads as¹²

$$\begin{eqnarray} &&W(\dot{\gamma }) = - K G^{2}_{ {{220\mbox{\it {q}}^{\prime }}}} k_2(\dot{\gamma }) \sin |\epsilon _2(\dot{\gamma })| \, {\rm Sgn}\, (\dot{\gamma }), \end{eqnarray} \tag{ C7 }$$

where K is a positive constant. To record the bias in the vicinity of this resonance, it is convenient to express an arbitrary tidal mode via the resonant one:

$$\begin{eqnarray} \nonumber \omega _{ {{220\mbox{\it {q}}}}}&=&(2+q)n-2\dot{\theta }=(2+q^{\prime })n-2\dot{\theta }+(q-q^{\prime })n \nonumber \\ &=&{-} 2\dot{\gamma }+(q-q^{\prime })n\approx (q-q^{\prime })n, \end{eqnarray} \tag{ C8 }$$

where we recall that $\dot{\gamma }$ vanishes at the point of resonance. Taking (C8) into account, we can write the bias as

$$\begin{eqnarray} \nonumber V&=&K\sum _{q\ne q^{\prime }}G^{2}_{220\mbox{\it {q}}} k_2(\omega _{ {{220\mbox{\it {q}}}}}) \sin | \epsilon _2(\omega _{ {{220\mbox{\it {q}}}}})| \, {\rm Sgn}\, (\omega _{ {{220\mbox{\it {q}}}}}) \nonumber \\ &=&K\sum _{q\ne q^{\prime }}G^{2}_{220\mbox{\it {q}}} k_2((q-q^{\prime })n) \sin | \epsilon _2((q-q^{\prime })n)| \, {\rm Sgn}\, (q-q^{\prime }), \nonumber\\ \end{eqnarray} \tag{ C9 }$$

For a planet that is slowing down, an estimate for the capture probability is derived from the consideration of two librations around the point of resonance $\dot{\gamma }=-\omega _{220\mbox{\it {q}}^{\prime }}=0$, i.e., of the last libration with a positive $\dot{\gamma }$ and the first libration with a negative $\dot{\gamma }$. Goldreich & Peale (1968) assumed that the energy offset from zero at the beginning of the last libration above the resonance is uniformly distributed between 0 and $\Delta E=\int \langle T\rangle \dot{\gamma }dt$. This assumption provided them with the following estimate for the probability of the capture:

$$\begin{equation} P_{\rm capt}=\frac{\delta E}{\Delta E}, \end{equation} \tag{ C10 }$$

where δE is the total change of the kinetic energy at the end of the libration below the resonance.

Thus, $\langle T\rangle \dot{\gamma }$ should be integrated over one cycle of libration, to obtain ΔE, and should be integrated over two librations symmetrically around the resonance $\omega _{ {{220\mbox{\it {q}}^{\prime }}}}=0$ to obtain δE. As a result, the odd part of the tidal torque at q = q' doubles in the integration for δE, whereas the bias vanishes. Both of these components are involved in the computation of ΔE.

This makes capture probability look like

$$\begin{eqnarray} P_{\rm capt}&=&\frac{2 \int _{-\pi }^{\pi }W(\dot{\gamma })d\gamma }{\int _{-\pi }^{\pi }\left[W(\dot{\gamma })+V\right]d\gamma } \approx \frac{2 \int _{-\pi }^{\pi }W(\dot{\gamma })d\gamma }{\int _{-\pi }^{\pi }W(\dot{\gamma })d\gamma +2\pi V}\nonumber\\ &=& \frac{ 2}{ {1+{ 2 \pi V}/{ \int _{-\pi }^{\pi }W(\dot{\gamma })d\gamma }}}. \end{eqnarray} \tag{ C11 }$$

The integral in this equation can be evaluated if we further assume, following Goldreich & Peale (1968), that in the vicinity of resonance the trajectory follows the singular separatrix solution (C6) that corresponds to the disappearance of the integral of motion E.