Chaotic Excitation and Tidal Damping in the GJ 876 System

In the absence of external gravitational perturbations, the tidal evolution of eccentricity occurs on a characteristic timescale (Yoder & Peale 1981):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{1}{e}\displaystyle \frac{{de}}{{dt}}\right)}^{-1}=\left(\displaystyle \frac{1}{n}\right)\left(\displaystyle \frac{2}{21}\right)\left(\displaystyle \frac{m}{{M}_{\star }}\right){\left(\displaystyle \frac{a}{R}\right)}^{5}\left(\displaystyle \frac{Q}{{k}_{2}}\right)={\tau }_{e},\end{eqnarray} \tag{ 1 }$$

where n is the mean motion, a is the semimajor axis, R is the radius of the orbiting body, Q is the tidal dissipation quality factor, and k₂ is the Love number. We note that in addition to eccentricity damping, tidal dissipation also facilitates semimajor axis decay. However, this process unfolds on a much longer characteristic timescale: ${\tau }_{a}={\tau }_{e}/2\,{e}^{2}\sim 30\,{\tau }_{e}$ and can be neglected for the moment (our numerical experiments will treat tidal dissipation self-consistently).

Expression (1) indicates that in the presence of only tidal forcing, a planet's eccentricity will exhibit nearly exponential decay. Calculating the true circularization timescale of an exoplanet directly using this expression, however, is abortive because first-principles based estimates of the tidal quality factor remain elusive, with lack of constraints on the planetary composition further aggravating the uncertainty. Instead, the oft-cited estimates tidal Q are based on observations of bodies within our own solar system (Goldreich & Soter 1966), with no faithful point of comparison for the tidal quality factor of 876 d, a super-Earth.

To motivate further calculations, let us first consider a scenario where the only relevant forcing on eccentricity comes from tidal dissipation. The dynamics would then be described by Equation (1), and the currently observed high eccentricity would be due to a tidal dissipation timescale on the order of the lifetime of the system, about 10 Gyr. Using the mass and orbital parameters from Table (1), we may draw upon the exoplanet mass–radius relationship obtained by Weiss et al. (2013) to estimate a planetary radius of R ∼ 3 R_⊕. Then, using a Love number of k₂ ∼ 0.4 (Dermott et al. 1988; Banfield & Murray 1992), we obtain a tidal quality factor on the order of Q ≳ 3 × 10⁶. This is a value considerably larger than that of solar system gas giants, so either 876 d is more than an order of magnitude less dissipative than these planets, or other factors must be important to the evolution of 876 d's eccentricity. Let us consider this alternative.

Table 1.  Adopted Orbital Fit of the GJ 876 System.

	M (M⊙)	a (au)	e	${ \mathcal M }$ (deg)	ϖ (deg)
⋆	0.37
d	2.25 × 10−5	0.02	0.12	301.87	219.10
c	7.99 × 10−4	0.13	0.25	139.32	115.96
b	2.54 × 10−3	0.21	0.03	74.62	110.86
e	4.98 × 10−5	0.35	0.03	184.11	129.56

Download table as:  ASCIITypeset image

Our first numerical experiment aims to examine the evolution of GJ 876 d's eccentricity in consideration of tidal forcing from its host star and the periodic gravitational perturbations of the giant planets in the system. We do this by considering the system in a nonchaotic regime by removing 876 e from our calculations. This renders the evolution of planets b and c strictly periodic, because the stochastic evolution of 876 e defines the chaotic nature of the system (Batygin et al. 2015a). At the same time, planet e's low mass and significant distance from 876 d result in negligible direct influence on the evolution of 876 d's orbit. Thus, the results of these simulations elucidate the behavior of 876 d's eccentricity in the absence of chaotic effects.

In our numerical experiments, we used shorter circularization times than those implied by conventional theory, so that many dynamical timescales can be simulated with less computation time. We note that this adjustment is valid as long as the timescales associated with angular momentum exchange among the planets remain much shorter than the circularization timescales we have adopted.

The results of a representative simulation are shown in Figure 1. Clearly, this calculation yields no explanation for the high eccentricity of planet d. That is, in presence of tidal forcing, no significant excitations in eccentricity exist. Instead, the behavior of 876 d's eccentricity is governed by the exponential decay from tides with small periodic perturbations from the outer planets. Thus, strictly periodic gravitational perturbations from the outer giant planets alone cannot explain the high eccentricity that is observed.

Zoom In Zoom Out Reset image size

We now consider the full system and its chaotic behavior. Unlike the simulation shown in Figure 1, here we model the full planetary system without considering tidal forces. Our results, shown in Figure 2, demonstrate that chaotic forcing alone is enough to drive eccentricity of 876 d to the observed value on a timescale of a few million years, much less than the lifetime of the system.

Zoom In Zoom Out Reset image size

With this insight into the importance of chaotic forcing, we perform our tidal dissipation numerical experiment again, this time accounting for the presence of 876 e. The behavior of planet d's eccentricity in this case is shown in Figure 3. While the dominant exponential decay of the eccentricity is similar, the periodic perturbations that manifested in the first experiment are replaced by chaotic excitations in the second experiment. Correspondingly, at sufficiently low eccentricities, these chaotic excitations dominate the behavior over the exponential decay of the tidal forcing.

Zoom In Zoom Out Reset image size

To quantify this relationship between the tidal and chaotic forcing, we ran a series of simulations sequentially increasing τ_e. A range of 20 different tidal circularization timescales were taken, and for each one, the simulation was performed over 10 circularization timescales. The chaotic nature of the system is acutely sensitive to initial conditions, making it difficult to extract precise information from the system at different circularization timescales—that is, each simulation is only representative of the true evolution of the system in a statistical sense. To account for these effects, 10 simulations were performed at every timescale, with a randomly generated initial mean anomaly for 876 d's orbit.² The results are shown in Figure 4, where rms eccentricity of planet d (measured over the last seven circularization timescales) is depicted as a function of τ_e. Clearly, as circularization time is increased, the mean eccentricity attained by planet d also grows. Quantifying this relationship is the focus of the next section.

Zoom In Zoom Out Reset image size

Numerical experiments presented in the preceding section suggest that the average eccentricity acquired by planet d after a period of chaotic equilibration is a monotonically increasing function of τ_e. To understand this relationship from analytic grounds, in this section we construct a simplified, perturbative model for the inner-most planet's motion. Given the large separation between the semimajor axes of planets d and c (as well as b) and the lack of a low-order commensurability between the orbital periods, it is sensible to assume that the interactions between these objects are secular in nature. Accordingly, in order to obtain a handle on the dynamical evolution of planet d, we employ Lagrange–Laplace perturbation theory (Murray & Dermott 1999).

To second order in the orbital eccentricities, the secular Hamiltonian that governs GJ 876 d's evolution reads:

$$\begin{eqnarray}{ \mathcal H } & = & \displaystyle \frac{-1}{\sqrt{{ \mathcal G }{Ma}}}\displaystyle \frac{{ \mathcal G }m^{\prime} }{4\,a^{\prime} }{\left(\displaystyle \frac{a}{a^{\prime} }\right)}^{2}[{b}_{3/2}^{(1)}\,{\rm{\Gamma }}\\ & & -{b}_{3/2}^{(2)}e^{\prime} \,\sqrt{2\,{\rm{\Gamma }}}\cos (\xi \,t-\gamma )]\end{eqnarray} \tag{ 2 }$$

where b's are Laplace coefficients, ${\rm{\Gamma }}=1-\sqrt{1-{e}^{2}}\simeq {e}^{2}/2$ and γ = − ϖ are scaled Poincare´ action-angle variables corresponding to planet d, while the primed quantities correspond to the perturbing body (i.e., planet c or b). For definitiveness, we assume that e' is nearly constant and the rate of the perturbing planet's pericenter regression, ξ, is time-independent. We note that both of these assumptions are justified by numerical simulations and simultaneously apply to planets b and c, as these objects have co-linear apsidal lines (Correia et al. 2010; Batygin et al. 2015a).

To leading order in the semimajor axis ratio, Hamilton's equation for planet d's longitude of perihelion is:

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma }{{dt}}=\displaystyle \frac{\partial { \mathcal H }}{\partial {\rm{\Gamma }}}\simeq -\displaystyle \frac{3}{4}n\,\left(\displaystyle \frac{m^{\prime} }{M}\right){\left(\displaystyle \frac{a}{a^{\prime} }\right)}^{3}+{ \mathcal O }{\left(\displaystyle \frac{a}{a^{\prime} }\right)}^{4},\end{eqnarray} \tag{ 3 }$$

where n is the mean motion of planet d, and we have expanded the Laplace coefficient as a hypergeometric series to leading order. Using Equation (3), it is trivial to check that $| d\gamma /{dt}| \ll \xi$. This disparity in precession frequencies explains why no coherent angular momentum exchange can be seen in Figure 1. That is, resonantly induced rapid regression of planet b and c's pericenter essentially renders the secular harmonic of Hamiltonian (2) a short-periodic effect that can be averaged over.

Accounting for the tidal decay of the orbital eccentricity (e.g., Hut 1981), the equation of motion for planet d's angular momentum deficit reads:

$$\begin{eqnarray}\displaystyle \frac{d{\rm{\Gamma }}}{{dt}} & = & -\displaystyle \frac{\partial {H}}{\partial \gamma }+\displaystyle \frac{\partial {\rm{\Gamma }}}{\partial e}{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{\mathrm{tide}}\\ & \simeq & \displaystyle \frac{15}{8\sqrt{2}}n{\left(\displaystyle \frac{a}{a^{\prime} }\right)}^{4}\,e^{\prime} \,\sqrt{{\rm{\Gamma }}}\sin (\xi \,t)-2\displaystyle \frac{{\rm{\Gamma }}}{{\tau }_{e}},\end{eqnarray} \tag{ 4 }$$

where we have neglected the slow evolution of γ and adopted ξ as the harmonic circulation frequency. Averaged over timescales much longer than 2π/ξ, chaotic variation in e' sin(ξ t), that arises from perturbations due to planet e, can be crudely modeled as a drift-free Weiner process, ${ \mathcal W }$ (Batygin et al. 2015b). Under this prescription, Equation (4) reduces to a readily integrable stochastic differential equation:

$$\begin{eqnarray}&&d\langle {\rm{\Gamma }}\rangle =2\sqrt{{ \mathcal D }\langle {\rm{\Gamma }}\rangle }d{ \mathcal W }-2\displaystyle \frac{\langle {\rm{\Gamma }}\rangle }{{\tau }_{e}}{dt},\end{eqnarray} \tag{ 5 }$$

where $\langle {\rm{\Gamma }}\rangle$ denotes the time-averaged angular momentum deficit, and ${ \mathcal D }$ is an effective diffusion coefficient that encompasses the chaotic nature of the perturbing orbits.

The standard deviation of the distribution function arising from Equation (5) is given by

$$\begin{eqnarray}&&{\sigma }_{\langle {\rm{\Gamma }}\rangle }=\sqrt{2\,\langle {{\rm{\Gamma }}}_{0}\rangle \,{ \mathcal D }\,{\tau }_{e}(\exp (-2t/{\tau }_{e})-\exp (-4t/{\tau }_{e}))},\end{eqnarray} \tag{ 6 }$$

where $\langle {{\rm{\Gamma }}}_{0}\rangle \ne 0$ is an (arbitrary) initial condition. The function (7) has a global maximum at t = log(2)τ_e/2, and asymptotically approaches zero as $t\to \infty$. This long-term behavior is unphysical and stems from the fact that $\langle {\rm{\Gamma }}\rangle =0$ is an absorbing boundary in Equation (5), while in reality it is reflective (i.e., a path that encounters $\langle {\rm{\Gamma }}\rangle =0$ within the framework of Equation (5) will remain at $\langle {\rm{\Gamma }}\rangle =0$ for all subsequent time, whereas real circular orbits can become eccentric³ ). Here, we circumvent this artificial limitation simply by evaluating ${\sigma }_{\langle {\rm{\Gamma }}\rangle }$ at its global maximum. Then, the standard deviation of the stochastic process (5) becomes

$$\begin{eqnarray}&&{({\sigma }_{\langle {\rm{\Gamma }}\rangle })}_{\max }=\sqrt{\langle {{\rm{\Gamma }}}_{0}\rangle \,{ \mathcal D }\,{\tau }_{e}/2}.\end{eqnarray} \tag{ 7 }$$

The key result of the above analysis is that the average eccentricity of planet d scales as the quarter-root of the circularization timescale:

$$\begin{eqnarray}&&\langle e\rangle \mathop{\propto }\limits_{\sim }{\tau }_{e}^{1/4}.\end{eqnarray} \tag{ 8 }$$

Note that this scaling was derived exclusively under the assumption of secular perturbations arising from a chaotic companion and should therefore be applicable even if planet d experiences a limited degree of inward tidal migration over the lifetime of the system. The relationship (8) is shown as a blue dashed line in Figure 4, and clearly provides an adequate match to our suite of numerical simulations, allowing us to extrapolate the existing results to more realistic values of ${\tau }_{e}$.

The preceding theory can be combined with the results of the numerical experiment described in Section 2.2 to estimate the lower bound of GJ 876 d's tidal quality factor. Specifically, within the context of our N-body simulations, the maximum value of the eccentricity attained by planet d over the modeled 10τ_e integration period is well fit by a quarter-root function:

$$\begin{eqnarray}&&\langle e\rangle =9.1\times {10}^{-4}{\left(\displaystyle \frac{{\tau }_{e}}{\mathrm{years}}\right)}^{1/4}.\end{eqnarray} \tag{ 9 }$$

This relationship, shown in Figure 5, can readily be extrapolated to the observed eccentricity to determine an approximate lower bound on the tidal circularization timescale of GJ 876 d. Using the parameters described in Table (1), we estimate that the maximum circularization timescale is of order ${\tau }_{e}\sim 3\times {10}^{8}$ years. Applying this value to Equation (1) and assuming a planetary radius of 3 R_⊕ (Weiss et al. 2013) and a Love number of k₂ = 0.4 as before (Dermott et al. 1988; Banfield & Murray 1992), we obtain a lower bound on the tidal quality factor of order Q ∼ 10⁵.

Zoom In Zoom Out Reset image size