Addressing the statistical mechanics of planet orbits in the solar system

TL;DR

This paper applies statistical mechanics to model the chaotic orbital dynamics of planets in the solar system, successfully reproducing probability distributions of orbital elements for major planets.

Contribution

It introduces a statistical mechanics framework to describe planetary orbit distributions, accurately modeling giant planets and inner planets without free parameters.

Findings

Microcanonical ensemble reproduces giant planet orbit statistics.

Equiprobability ansatz models Earth and Venus PDFs accurately.

Reasonable agreement with Mars and Mercury inclination PDFs.

Abstract

The chaotic nature of planet dynamics in the solar system suggests the relevance of a statistical approach to planetary orbits. In such a statistical description, the time-dependent position and velocity of the planets are replaced by the probability density function (PDF) of their orbital elements. It is quite natural to set up this kind of approach in the framework of statistical mechanics. In the present paper I focus on the collisionless excitation of eccentricities and inclinations by gravitational interactions in a planetary system, the prototype of such a dynamics being the future planet trajectories in the solar system. I thus address the statistical mechanics of the planetary orbits in the solar system and try to reproduce the PDFs numerically constructed by Laskar (2008). I show that the microcanonical ensemble of the Laplace-Lagrange theory accurately reproduce the statistics of the giant planet orbits. To model the inner planets I then investigate the ansatz of equiprobability in the phase space constrained by the secular integrals of motion. The eccentricity and inclination PDFs of Earth and Venus are reproduced with no free parameters. Within the limitations of a stationary model, the predictions also show a reasonable agreement with Mars PDFs and that of Mercury inclination. The eccentricity of Mercury demands in contrast a deeper analysis. I finally revisit Laskar's random walk approach to the time dependence of the inner planet PDFs. Such a statistical theory could be combined with direct numerical simulations of planet trajectories in the context of planet formation, which is likely to be a chaotic process.

Equations

$$\mathcal{H}=\sum _{k=1}^N\left(\frac{\text{@vec}{p}_k^2}{2{\mu }_k}-G\frac{m_0{m}_k}{r_k}\right)+\sum _{k=1}^N\sum _{l=k+1}^N\left(\frac{\text{@vec}{p}_k\cdot \text{@vec}{p}_l}{m_0}-G\frac{m_k{m}_l}{|\text{@vec}{r}_k-\text{@vec}{r}_l|}\right)$$

$$\begin{array}{rl}{\Lambda }_k& ={\mu }_k\sqrt{G(m_0+m_k)a_k}\\ P_k& ={\Lambda }_k(1-\sqrt{1-e_k^2})\\ Q_k& =2{\Lambda }_k\sqrt{1-e_k^2}{\sin }^2(i_k/2)\end{array}$$

$$⟨\mathcal{H}⟩=\frac{1}{(2\pi {)}^N}{\int }_0^{2\pi }d{M}_1\cdots {\int }_0^{2\pi }d{M}_N\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_{\text{sec}}$$

$${\mathcal{H}}_{\text{sec}}=-\sum _{k=1}^N\sum _{l=k+1}^N\frac{G{m}_k{m}_l}{(2\pi {)}^2}{\int }_0^{2\pi }{\int }_0^{2\pi }\frac{d{M}_{k}d{M}_l}{|\text{@vec}{r}_k-\text{@vec}{r}_l|}$$

$${\mathcal{H}}_{\text{sec}}=-\sum _{k=1}^N\sum _{l=k+1}^N\frac{G{m}_k{m}_l}{(2\pi {)}^2}{\int }_0^{2\pi }{\int }_0^{2\pi }\frac{d{M}_{k}d{M}_l}{|\text{@vec}{r}_k-\text{@vec}{r}_l|}-\sum _{k=1}^N\frac{3G^2{m}_0^2{m}_k}{c^2{a}_k^2\sqrt{1-e_k^2}}$$

$${\mathcal{H}}_{\mathrm{LL}}=-\frac{\text{@vec}{x}^{\mathrm{T}}\text{@tens}A\text{@vec}x+\text{@vec}{y}^{\mathrm{T}}\text{@tens}A\text{@vec}y+\text{@vec}{v}^{\mathrm{T}}\text{@tens}B\text{@vec}v+\text{@vec}{z}^{\mathrm{T}}\text{@tens}B\text{@vec}z}{2}$$

$$\begin{array}{rl}{x}_k& =\sqrt{{\Lambda }_k}{e}_k\sin p_k\\ y_k& =\sqrt{{\Lambda }_k}{e}_k\cos p_k\end{array}$$

$$\text{@tens}A=\text{@tens}{O}_A{D}_A{O}_A^{\text{/}\mathrm{T}},\text{@tens}B=\text{@tens}{O}_B{D}_B{O}_B^{\text{/}\mathrm{T}}$$

$$\begin{array}{rl}\text{@vec}x& =\text{@tens}{O}_A\text{@vec}{x}^{\prime }\\ \text{@vec}y& =\text{@tens}{O}_A\text{@vec}{y}^{\prime }\end{array}$$

$$\begin{array}{rl}{x}_k^{\prime }& =\sqrt{2P_k^{\prime }}\sin p_k^{\prime }\\ y_k^{\prime }& =\sqrt{2P_k^{\prime }}\cos p_k^{\prime }\end{array}$$

$${\mathcal{H}}_{\mathrm{LL}}=-\sum _{k=1}^N\left(g_k{P}_k^{\prime }+s_k{Q}_k^{\prime }\right)$$

$$L_x$$

$$L_y$$

$$C$$

$$C=\sum _{k=1}^N{P}_k+Q_k=\sum _{k=1}^N{P}_k^{\prime }+Q_k^{\prime }$$

$$L_x$$

$$L_y$$

$$v_5^{\prime }=\sum _{k=1}^N\text{@tens}{O_B^{\text{/}\mathrm{T}}}_{5k}{v}_k=-\frac{L_x}{2\sqrt{{\sum }_{k=1}^N{\Lambda }_k}}$$

$$z_5^{\prime }=\sum _{k=1}^N\text{@tens}{O_B^{\text{/}\mathrm{T}}}_{5k}{z}_k=-\frac{L_y}{2\sqrt{{\sum }_{k=1}^N{\Lambda }_k}}$$

$$\rho (\text{@vec}{P}^{\prime },\text{@vec}{Q}^{\prime },\text{@vec}{p}^{\prime },\text{@vec}{q}^{\prime })=\frac{\delta (q_5^{\prime }-\overline{q_5^{\prime }})}{(2\pi {)}^{2N-1}}\prod _{k=1}^N\delta (P_k^{\prime }-\overline{P_k^{\prime }})\delta (Q_k^{\prime }-\overline{Q_k^{\prime }})$$

$$\rho (e_k)=e_k{\int }_0^{+\infty }duu{J}_0(e_{k}u)\prod _{n=1}^N{J}_0({\gamma }_{kn}u)$$

$$\rho (\text{@vec}{P}^{\prime },\text{@vec}{Q}^{\prime },\text{@vec}{p}^{\prime },\text{@vec}{q}^{\prime })\propto \delta (C-\overline{C})f$$

$$f(\text{@vec}{P}^{\prime },\text{@vec}{Q}^{\prime },\text{@vec}{p}^{\prime },\text{@vec}{q}^{\prime })=\delta (q_5^{\prime }-\overline{q_5^{\prime }})\prod _{k=5}^8\delta (P_k^{\prime }-\overline{P_k^{\prime }})\delta (Q_k^{\prime }-\overline{Q_k^{\prime }})$$

$$\rho (\text{@vec}{P}^{\prime },\text{@vec}{Q}^{\prime },\text{@vec}{p}^{\prime },\text{@vec}{q}^{\prime })\propto \delta (C-\overline{C})\delta ({\mathcal{H}}_{\mathrm{LL}}-\overline{{\mathcal{H}}_{\mathrm{LL}}})f$$

$$\rho (\text{@vec}{P}^{\prime },\text{@vec}{Q}^{\prime },\text{@vec}{p}^{\prime },\text{@vec}{q}^{\prime })\propto \delta (C-\overline{C})\delta ({\mathcal{H}}_{\mathrm{LL}}-\overline{{\mathcal{H}}_{\mathrm{LL}}})\delta (P_1^{\prime }-\overline{P_1^{\prime }})f$$

$$\rho (h_k)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }du{e}^{-i{h}_{k}u}\prod _{n=1}^N{J}_0\left({\gamma }_{kn}u\right)$$

$$\text{@vec}r=\text{@vec}{r}_{\sec }+\text{@vec}{r}_{\mathrm{chaos}}$$

$${\rho }_{\sec }(x_k^{\prime },y_k^{\prime })=\frac{\delta \left[r-{\left(2\overline{P_k^{\prime }}\right)}^{1/2}\right]}{2\pi {\left(2\overline{P_k^{\prime }}\right)}^{1/2}}$$

$$\text{@vec}{\dot{r}}_{\mathrm{chaos}}=\sqrt{2D_k}\text{@vec}\eta (t)$$

$$⟨\text{@vec}\eta (t)⟩$$

$$⟨{\eta }_x(t){\eta }_x(t^{\prime })⟩$$

$$⟨{\eta }_x(t){\eta }_y(t^{\prime })⟩$$

$$\text{@vec}{r}_{\mathrm{chaos}}(t=0)=\text{@vec}0$$

$$\rho (\text{@vec}r,t)={\int }_{{\mathbb{R}}^2}d\text{@vec}z{\rho }_{\sec }(\text{@vec}z){\rho }_{\mathrm{chaos}}(\text{@vec}r-\text{@vec}z,t)$$

Full text

1 Introduction

The chaotic dynamics of solar system planets (Laskar, 1989; Sussman & Wisdom, 1992) raises the question of a statistical approach to planetary orbits. When chaos is significant, a single integration of the equations of motion is not representative of the entire possible dynamics. A description based on the probability density function (PDF) of the planet orbital elements, instead of their time-dependent position and velocity, is essentially more suitable. Such a statistical approach could be combined with usual direct numerical integrations, especially in the context of planet formation, which is likely to be a chaotic process (e.g., Hoffmann et al., 2017).

The role of statistical mechanics in assessing the PDF of the planet orbital elements naturally emerges. Recently, the statistical mechanics of terrestrial planet formation has been addressed by Tremaine (2015). In the spirit of a packed planetary systems hypothesis (Fang & Margot, 2013), he advances the equiprobability of phase-space configurations verifying a certain criterion for the long-term stability of planetary systems. Within the limitations due to the sheared-sheet approximation and the restriction to the planar case, this ansatz allows to analytically compute any property derivable from the complete N-planet distribution function, such as the distributions of planet eccentricities and semimajor axis differences. The predictions show an encouraging agreement with N-body simulations and data from the Kepler catalogue, while significant discordances arise in the case of radial-velocity data and solar system planets. The simplicity of Tremaine’s ansatz is fundamental to describe at once the two main steps of the final giant-impact phase of terrestrial planet formation: the excitation of embryo eccentricities and inclinations through mutual gravitational interactions and the subsequent collisions and merging. However, one could ask for a more fundamental approach, such as one exploiting the integrals of motion (e.g., Kocsis & Tremaine, 2011; Touma & Tremaine, 2014; Petrovich & Tremaine, 2016), in line with the spirit of equilibrium statistical mechanics (Landau & Lifshitz, 1969b). Such an approach can be actually set up by focusing on how eccentricities and inclinations are excited by gravitational interactions in a early, collisionless final phase of terrestrial planet formation, and postponing the analysis of collisions to subsequent treatments. As there is no fundamental physical difference between this problem and that of future planetary trajectories in the solar system (Laskar, 1996), the latter can provide a benchmark to any statistical theory of planetary orbits to be applied to exoplanet systems. In the present study I address the problem of setting up the statistical mechanics of planetary trajectories in the solar system.

The statistical mechanics of gravitating systems is notoriously challenging (Padmanabhan, 1990). Firstly, the system has to be confined in a spherical box to assure a bounded phase space, similarly to ideal gases. Then, the non-extensive nature of gravitational energy, due to the long-range character of Newtonian potential, breaks the equivalence between the microcanonical and canonical ensembles, which is typical of systems like neutral gases and plasmas. The long-range nature of gravity also prevents the existence of a thermodynamic limit. In addition to this, the short-range singularity of the gravitational potential requires to take into account the specific nature of small-scale interactions to guarantee the existence of microcanonical equilibrium states. Finally, the number $N$ of planets in a typical planetary system, even during the giant-impact phase of its formation process, is limited. Therefore, the $N\gg 1$ regime is generally inappropriate and so is the employment of a mean field approach111The lack of the $N\rightarrow+\infty$ limit could seem to question the usefulness of employing statistical mechanics. Actually, the foundation for such an approach relies on the chaotic nature of planet dynamics, independently of their number $N$.. As pointed out by Touma & Tremaine (2014) mainly in the context of stellar systems, one effective way to construct the statistical mechanics of orbits in a planetary system is to average the planet motion over its fastest timescale, that is the orbital period. As a result of this averaging procedure, each planet is replaced by a massive ring following its Keplerian orbit, whose linear mass density is inversely proportional to the planet orbital velocity, a so-called Gaussian ring (Kondratyev, 2012). In such a system, the rings interact with each other via secular gravitational interactions. Their eccentricities and inclinations relax by exchange of angular momentum, while their semi-major axes are constant, and so are their Keplerian energies (Rauch & Tremaine, 1996). Even if there is no theorem assuring this secular dynamics approaches the actual planet motion, it can nevertheless represent the source of important results (Arnold, 1978). Indeed, the first clear indications of a chaotic planetary motion in the solar system came from averaged equations of motion (Laskar, 1989). Moreover, by using the same secular dynamics, Laskar (2008, from now on L08) numerically computed the PDFs of eccentricity and inclination of the solar system planets. In the present paper I address the problem of reproducing these PDFs through a statistical mechanics approach.

The paper is structured as follows. In Sect. 2 I briefly recall how secular planet dynamics can be introduced in its simplest form. Then, in Sect. 3 I describe the PDFs of eccentricity and inclination of the solar system planets calculated numerically in L08. I introduce in Sect. 4 the microcanonical ensemble of the Laplace-Lagrange theory to reproduce the statistics of the giant planet orbits. The ansatz of equiprobability in the phase space constrained by the secular integrals of motion is presented in Sect. 5 for the inner planets. Finally, I revisit in Sect. 6 the random walk ansatz of L08 to account for the time dependence of the inner planet PDFs. I conclude emphasizing the relevance of a statistical theory of planetary orbits in the context of planet formation.

2 Secular planet dynamics

Using standard notation (Morbidelli, 2002), the Hamiltonian of the $N=8$ solar system planets can be written as

[TABLE]

where heliocentric canonical variables $\@vec{r}_{k},\@vec{p}_{k}$ are employed, $m_{0}$ is the Sun mass, $m_{k}$ the planet masses, $\mu_{k}=m_{0}m_{k}/(m_{0}+m_{k})$ the reduced masses and $G$ the gravitational constant. As usual, the index $k$ lists the planets by increasing semi-major axis, from Mercury to Neptune. The first term on the right side of equation (1) is the Hamiltonian of the Keplerian motion, while the second one contains the gravitational interactions among planets. Therefore, it is worthwhile to introduce a set of action-angle variables for the Keplerian motion, e.g. the modified Delaunay variables (Morbidelli, 2002):

[TABLE]

Keplerian orbital elements are used in these definitions: $a_{k}$ is the planet orbit semi-major axis, $e_{k}$ the eccentricity, $i_{k}$ the inclination, $M_{k}$ the mean anomaly, $\omega_{k}$ the argument of perihelion and $\Omega_{k}$ the longitude of node. As long as planets do not experience close encounters and do not lie near mean motion resonances, secular dynamics can be introduced by averaging the Hamiltonian (1) over the fastest motion timescales, i.e. the planet orbital periods. This corresponds to average the Hamiltonian over the mean anomalies $M_{k}$:

[TABLE]

where $\mathcal{H}_{0}=\sum_{k=1}^{N}-G(m_{0}+m_{k})^{2}\mu_{k}^{3}/(2\Lambda_{k}^{2})=\sum_{k=1}^{N}-Gm_{0}m_{k}/(2a_{k})$ is the total Keplerian energy and

[TABLE]

are the orbit-averaged gravitational interactions between planets (the kinetic contribution in the second term on the right side of Eq. (1) averages out to zero over a Keplerian orbit). This averaging procedure corresponds to constructing a secular normal form to first order in planetary masses (Morbidelli, 2002). The action variables $\Lambda_{k}$ are integrals of motion of the secular dynamics (Arnold, 1978), and so are the semi-major axes $a_{k}$. As a consequence, $\mathcal{H}_{0}$ is an additive constant which may be dropped from the Hamiltonian. The single-planet secular phase space is four-dimensional and compact, with the canonical volume element given by $dp_{k}dq_{k}dP_{k}dQ_{k}$. The secular contribution of general relativity is of primary importance for the long-term dynamics of the solar system inner planets (e.g., L08). The secular Hamiltonian accounting for the leading relativistic correction is given by

[TABLE]

where $c$ is the speed of light. A short derivation of the general relativistic contribution is presented in the Appendix A for future reference.

When eccentricities and inclinations are sufficiently small, $e_{k},\sin(i_{k}/2)\ll 1$, it is valuable to develop the secular Hamiltonian in a power series of these variables. Neglecting terms of the fourth order and smaller, one obtains the Laplace-Lagrange (LL) linear secular theory:

[TABLE]

where I have introduced the Poincaré canonical variables,

[TABLE]

and, for instance, $\@vec{x}^{\rm T}=(x_{1},\dots,x_{N})$. The elements of the matrices $\@tens{A}$ and $\@tens{B}$ are given in the Appendix B. These matrices are real and symmetric and can therefore be diagonalized through orthogonal matrices $\@tens{O_{A}}$ and $\@tens{O_{B}}$,

[TABLE]

The elements of the diagonal matrices $\@tens{D_{A}}$ and $\@tens{D_{B}}$ are the eigenvalues of $\@tens{A}$ and $\@tens{B}$, respectively, and the columns of $\@tens{O_{A}}$ and $\@tens{O_{B}}$ are their normalized eigenvectors. New action-angle coordinates $\@vec{P^{\prime}},\@vec{Q^{\prime}},\@vec{p^{\prime}},\@vec{q^{\prime}}$ can be introduced by combining the following canonical transformations:

[TABLE]

[TABLE]

Employing the new variables, the LL Hamiltonian becomes:

[TABLE]

where $g_{k}$ and $s_{k}$ are the eigenvalues of the matrices $\@tens{A}$ and $\@tens{B}$, respectively. According to the corresponding Hamilton equations, the actions $P^{\prime}_{k}$ and $Q^{\prime}_{k}$ are integrals of motion, while the angles $p^{\prime}_{k}$ and $q^{\prime}_{k}$ change linearly with time at constant frequencies $-g_{k}$ and $-s_{k}$, respectively. Therefore, the time dependence of the Poincaré variables $x_{k},y_{k},v_{k},z_{k}$ turns out to be the superposition of $N$ harmonics with different frequencies and amplitudes. In the present study I refer to each of these harmonics as a LL mode and I follow the conventional ordering of the frequencies $g_{k}$ and $s_{k}$ (Morbidelli, 2002, Table 7.1).

2.1 Integrals of motion

Generally speaking, the dynamics of a planetary system obeys to the conservation of mechanical energy, momentum and angular momentum. The secular dynamics introduced above automatically verifies the conservation of momentum: as Keplerian orbits are closed, the average of planet momenta over mean anomalies identically vanishes. Therefore, the dynamics resulting from the secular Hamiltonian (5) must conserve energy and angular momentum. This can be easily shown to be the case in the LL dynamics (e.g., Kocsis & Tremaine, 2011). As formula (11) shows, the LL Hamiltonian is a function of the action variables only. These being integrals of motion, the conservation of energy follows. The angular momentum content of a planetary system can be described by the following quantities:

[TABLE]

where $x,y,z$ are the Cartesian coordinates related to the definition of the Keplerian elements $i_{k},\omega_{k}$ and $\Omega_{k}$ (Morbidelli, 2002)222Typically, the $x\,y$-reference plane is chosen to be perpendicular to the total angular momentum vector. However, in this study I choose the ecliptic plane as the reference one, since this is the choice made in L08.. $L_{x}$ and $L_{y}$ are the components of the angular moment lying on the $x\,y$-reference plane, while $C$ is the angular momentum deficit (AMD) (Laskar, 1997, 2000), the difference between the angular momentum that planets would have on coplanar, circular orbits and $L_{z}$, the component of the total angular momentum perpendicular to the reference plane. The AMD can be expressed as

[TABLE]

where I have used the definitions (2) and the fact that the matrices $\@tens{O_{A}}$ and $\@tens{O_{B}}$, appearing in the transformation (9), are orthogonal. As the action variables $\@vec{P^{\prime}}$ and $\@vec{Q^{\prime}}$ are conserved in the LL dynamics, the AMD is also a conserved quantity. Conservation of $L_{x}$ and $L_{y}$ derives from the properties of the matrix $\@tens{B}$. As shown in the Appendix B, one of the eigenvalues $s_{k}$ is zero and $\left.\left(\!\!\sqrt{\Lambda_{1}},\dots,\!\sqrt{\Lambda_{N}}\right)\,\right/\!\!\sqrt{\sum_{k}\Lambda_{k}}$ is, up to a sign, the corresponding normalized eigenvector. By convention (Morbidelli, 2002, Table 7.1), the null eigenvalue is chosen to be $s_{5}$. Moreover, neglecting terms of order three in eccentricities and inclinations, one has

[TABLE]

All this implies that in the LL dynamics the following identities are verified:

[TABLE]

where I have used the fact that $\@tens{O_{B}}$ is orthogonal, $\@tens{O_{B}^{\/-1}}=\@tens{O_{B}^{\/\rm T}}$, and its columns are the normalized eigenvectors of $\@tens{B}$. As the frequency $s_{5}$ is zero, $q^{\prime}_{5}$ does not change with time. Therefore, $v^{\prime}_{5}$ and $z^{\prime}_{5}$ are integrals of motion and so are $L_{x}$ and $L_{y}$.

3 The PDFs of Laskar (2008)

In L08 Laskar performs a statistical analysis of the future planet orbits in the solar system. He employs 1001 integrations of secular equations over the 5 Gyr of the Sun’s remaining lifetime before its red giant phase. These integrations differ for the initial conditions, which are obtained with small variations of the initial Poincaré variables of the VSOP82 solution (Bretagnon, 1982). The total integration time is divided in 250 Myr intervals and statistics are performed over each interval recording the state of the 1001 solutions with a 1 kyr timestep. Normalized PDFs of planet eccentricities and inclinations are then estimated. The PDFs of the giant planets are plotted in Fig. 1, while in Fig. 2 are shown those of the inner ones333Unfortunately, the data behind the L08 PDFs are not currently available (Laskar, private communication). Therefore, I have extracted these curves from the original plots through WebPlotDigitizer (Rohatgi, 2017). This works fine for the giant planets PDFs and the eccentricity PDFs of the inner ones, as Laskar plotted them separately at three different times. Because of chaotic diffusion, such an extraction is infeasible for the inclination PDFs of the inner planets, which Laskar plotted at several different times on the top of each other. In this case I use as a reference Laskar’s fits instead of the original numerical PDFs, even if some differences exists between these curves, especially in the tails.. The dynamics considered by Laskar is more accurate than the one I introduced in Sect. 2, as Laskar’s secular Hamiltonian contains terms of order two in planetary masses and six in eccentricities and inclinations. This is why he conjectures that, even though his PDFs are obtained through secular equation, they should be nevertheless close to those arising from the full, non-averaged, dynamics. While analysing these PDFs, what strikes Laskar is the very different shape between the giant planets curves and those of the inner planets. The PDFs of the outer planets are characterized by two peaks and restricted to a certain range of eccentricities and inclinations. In contrast, those of the inner planets have zero value at $e=0$ and $i=0$, a single peak and a continuous decaying at large eccentricities and inclinations. On the basis of these differences, Laskar takes two different approaches in explaining his numerical PDFs. By means of a frequency analysis applied to the outer planet motion, he shows that the PDF of a quasiperiodic approximation of their dynamics can reproduce very accurately the numerical PDFs. On the other hand, to take into account the significant chaotic diffusion existing in the statistics of the inner planets, Laskar fits their PDFs through a Rice distribution (Rice, 1945). Indeed, he assumes that, because of chaotic diffusion, the Poincaré variables (7) become practically independently Gaussian distributed over long timescales, with non-zero mean and a variance that increases linearly with time, like in a diffusion process.

4 Giant planets

According to secular dynamics, the orbital motion of the giant planets is well approximated by a quasiperiodic time function (L08). Indeed, chaotic diffusion is very limited for the outer planets, their PDFs virtually do not change with time. A certain chaotic behaviour arises from the full equations of motion (Sussman & Wisdom, 1992), but according L08 this is irrelevant to the overall structure of the PDFs. Taking advantage of this regularity of the giant planets, it is straightforward to predict the PDF of their Keplerian orbital elements. Neglecting chaotic diffusion, the most basic quasiperiodic dynamics is that resulting from the LL Hamiltonian, Eq. (11). Such a motion is ergodic and the asymptotic stationary probability density is the microcanonical one (Arnold, 1978). It is straightforward to write down the corresponding PDF employing the integrals of motion (Landau & Lifshitz, 1969b):

[TABLE]

where $\delta$ stands for the Dirac delta and the bar indicates the initial value of the corresponding variable444In the present study these initial values are chosen to match the initial conditions of the VSOP82 solution (Bretagnon, 1982, Tables 4 and 5), since they are the ones employed in L08.. This PDF reflects the conservation of the action variables $\@vec{P^{\prime}},\@vec{Q^{\prime}}$ in the LL dynamics, while the angle variables are uniformly distributed over the (2N$$-$$1)-dimensional torus $\mathbf{T}^{2N-1}$ (the angle $q^{\prime}_{5}$ is conserved because the angular momentum is, see Sect. 2.1). An integral formula for the PDFs of eccentricity and inclination can be obtained (Mardia, 1972),

[TABLE]

where $J_{0}$ is the Bessel function of the first kind and order zero and $\gamma_{kn}=\@tens{O_{A}}_{kn}(2\overline{P^{\prime}_{n}}/\Lambda_{k})^{1/2}$. A similar formula holds for the inclination $i_{k}$. Otherwise, these PDFs can be very rapidly estimated by direct sampling of Eq. (16) and using transformations (10), (9) and (7). The distribution function (16) is time-independent and describes the statistics of planet orbits over timescales much longer than the timescale of the LL dynamics (dozen of kyr to few Myr for the solar system planets). Since Laskar’s PDFs are constructed over 250-Myr intervals, one can expect the statistics of the giant planets resulting from Eq. (16) to agree with that of L08. Therefore, I plot in Fig. 1 the PDFs predicted by the microcanonical distribution (16), along with Laskar’s PDFs. The agreement between the two curves is indeed very satisfactory for the inclination PDFs, with a very good matching for Jupiter and Saturn and minor differences in the case of Uranus and Neptune. The agreement is also good for the eccentricity PDFs of Jupiter and Saturn, while some major discordances arise in the eccentricity PDFs of Uranus and Neptune. Even though the shape of the PDFs is correctly reproduced by the microcanonical density (16) in both cases, the endpoints of the eccentricity interval over which the PDFs vary are somewhat different from L08. However, the reason for such a discrepancy is clear. The LL Hamiltonian (11) employed to derive the microcanonical distribution (16) is valid up to the third order in eccentricities and inclinations. This regime is meaningful for the giant planets of the solar system, whose eccentricities and inclinations are quite moderate. However, terms of the fourth order in $e$ and $i$ in the Hamiltonian (5) yield non-linearities in the corresponding Hamilton equations and produce high-order harmonics in the eccentricity and inclination time dependence555These higher order terms also affect amplitude and frequency of the LL modes. It is worthwhile to remember that the secular Hamiltonian (5) is itself valid to the first order in planetary masses. Contributions from higher order terms in $m_{k}$ slightly adjust amplitude and frequency of the harmonics.. The amplitude of these additional harmonics is generally much smaller than that of the LL modes, but some of them can nevertheless contribute to a significant level. Bretagnon (1974, Tables 12 and 13) reports the amplitudes of these high-order harmonics. From Bretagnon’s tables it is clear that their contribution is bigger for the eccentricities than the inclinations and for Uranus and Neptune than Jupiter and Saturn. This is in agreement with Fig. 1. Moreover, the amplitudes of the higher order harmonics are roughly in accord with the size of the discrepancies in Fig. 1. Generally speaking, the agreement between the predictions of the microcanonical distribution (16) and Laskar’s PDFs is very satisfactory when one realizes the simplicity of the LL approximation. Moreover, a better agreement can be obtained extending the analysis of Sect. 2 to higher order terms in eccentricities, inclinations and masses. In principle, this could be done by means of a perturbation approach to the Hamiltonian (1) based on successive quasi-identity canonical transformations and normal forms (Morbidelli, 2002).

5 Inner planets

The analysis of the giant planet orbits is based on the assumption that the effect of chaos on their motion is largely negligible. This is not the case for the inner planets. Laskar (1994) already showed that Mercury and Mars orbits are affected by a significant chaotic diffusion over 5 Gyr, while this diffusion is moderate for Venus and the Earth. A first interesting step in analysing the structure of Laskar’s PDFs of the inner planets is to consider what the distribution function (16) predicts for their eccentricity and inclination. This is illustrated in Fig. 2, where I also plot the corresponding L08 PDFs at three different times, $t=$ 500 Myr, 2.5 Gyr and 5 Gyr, to illustrate chaotic diffusion. As one might expect, there is no agreement between the two curves. Nevertheless, it is interesting to note that, contrary to what could emerge from Laskar’s analysis in L08, the difference in shape between the giant and inner planet PDFs is not fundamental. Fig. 2 shows that Eq. (16) based on the LL dynamics is a priori able to produce the general shape of the inner planet PDFs, i.e. the zero value and positive linear slope at $e=0$ and $i=0$, and the continuous decaying at large eccentricities and inclinations. This is particularly evident in the Venus and Earth eccentricity curves and is due to the strong linear coupling existing between several LL modes in the Earth and Venus dynamics. Analogue strong couplings are not present in the LL dynamics of the giant planets. This explains the different shape of their PDFs, with two peaks and a zero probability outside a certain range (Fig. 1). This is also the case for Mercury, what justifies the strong differences between the corresponding Laskar’s PDFs and Eq. (16) in Fig. 2.

5.1 The ansatz of equiprobability of phase-space states

Since the effect of chaos is relevant for the long-term dynamics of the inner planets, I search for a more pertinent statistical description than one based on quasiperiodic motion. The chaotic nature of planet dynamics allows for sharing of angular momentum between the LL action variables, which would otherwise be constant. This chaotic diffusion of angular momentum is already mentioned in Laskar (1996) and has been recently illustrated in Wu & Lithwick (2011). The simplest statistical approach is assuming that, by consequence, the motion of the inner planets uniformly visits the subset of the phase space defined by the conservation of angular momentum $\@vec{L}$ and energy $\mathcal{H}_{\mathrm{sec}}$. Clearly, the fact that chaotic diffusion finally lead to such an uniform phase-space measure cannot be strictly true, as the inner planet PDFs depend on time. More probably, an uniform exploration could indeed occur adiabatically on a smaller subset of the phase space than that allowed by the conservation of angular momentum and energy only. Such a domain could be in principle determined by studying the higher order terms in the Hamiltonian (5). However, the hypothesis of equiprobability in the phase-space is still the best estimate that one can make without a deeper analysis of the dynamics (Jaynes, 1957). It so simple that it is worthwhile to investigate its predictions before undertaking more involved studies. Indeed, the idea of equipartition between LL modes has already been considered from a dynamical perspective by Wu & Lithwick (2011). Moreover, even though it leads to a stationary distribution function, such an ansatz is nevertheless quite appropriate to Earth and Venus, since their PDFs diffuse very slowly over time, and can be still useful to understand the main factors contributing to the overall structure of Mercury and Mars PDFs.

Since the giant planets are well described by the microcanonical distribution (16), in which all the LL action variables are fixed, it is clear that an equiprobability simply based on conservation of the total angular momentum and energy of all the planets would not produce reasonable results. To achieve a global description of giant and inner planets at once, one needs to statistically disconnect them. To this end, I note that the total angular momentum of the inner planets is quite independent of that of the giant planets. More precisely, Laskar (1997) has shown that the flux of angular momentum deficit (AMD) from the large reservoir of the outer planets to the inner ones is moderate over 5 Gyr. As a first approximation, one can therefore assume that the total AMD of the inner planets is constant. Moreover, from Eq. (13) one can interpret the action variables $P^{\prime}_{k}$ and $Q^{\prime}_{k}$ as the AMD content in the corresponding LL mode. These considerations suggest to fix the action variables corresponding to the modes $k\in\{5,6,7,8\}$, as they are the only ones which are relevant for the dynamics of the outer planets666As stated in Sect. 2, I adopt the conventional ordering of the LL modes (Morbidelli, 2002, Table 7.1)..

Indeed, with this choice the conservation of the total AMD reduces to that of the quantity $\sum_{k=0}^{4}P^{\prime}_{k}+Q^{\prime}_{k}$, whose initial value is very close to the initial total AMD of the inner planets, as shown in Table 1. The remaining action variables, corresponding to the modes $k\in\{1,2,3,4\}$, can then be allowed to vary stochastically according to equiprobability. The statistics of the giant planet orbits would therefore be virtually identical to the one shown in Fig. 1.

Among the integrals of motion, the AMD is of particular interest. For instance, Laskar (1997, 2000) has highlighted its role in the orbital configuration of planetary systems. An important reason for this relevance is that, differently from Eqs. (11) and (14), the simple Eq. (13) is valid to all orders in eccentricities and inclinations. I then start to investigate the predictions of equiprobability based on AMD conservation only. Conservation of energy will be added later. The ansatz corresponding to the above considerations reads:

[TABLE]

Conservation of $L_{x}$ and $L_{y}$ to second order in eccentricities and inclinations is contained in this distribution function as the variables $Q^{\prime}_{5}$ and $q^{\prime}_{5}$ are set to their initial values (see Eq. (15)). I do not have a simple analytical formula for the PDFs of eccentricity and inclination. However, they can be calculated numerically very rapidly by direct sampling, as illustrated in the Appendix D. Moreover, I present an analytical characterization of these PDFs in the Appendix C. The predictions of the ansatz (18) are compared in Fig. 3 with Laskar’s PDFs. They reproduce very satisfactorily Venus and Earth PDFs. A detailed comparison in the case of Earth is shown in Figs. 4 and 5. The agreement is remarkable if one considers the simplicity of the ansatz and that it does not contain any free parameter777To fit his PDFs Laskar employs three free parameters for each planet and each orbital element (eccentricity or inclination).. For what concerns the PDFs of Mars and that of Mercury inclination, the predictions are still reasonable as the differences with the numerical PDFs are of the same order of magnitude of the PDF diffusion over time. Indeed, in these cases one could a priori interpret the predicted curves as the long-term, stationary limit of Laskar’s PDFs. However, this interpretation breaks down when one considers the eccentricity PDF of Mercury, as shown in detail in Fig. 6. In this case, the predicted PDF is substantially different from the numerical one. In particular, the mean of Laskar’s curve is not correctly reproduced by the ansatz. It is clear that equiprobability, as applied in Eq. (18), is far to approximate the actual Mercury dynamics.

One can add to the ansatz (18) the conservation of secular energy at quadratic order in eccentricities and inclinations (see Eq. (11)),

[TABLE]

In the Appendix D I suggest an algorithm to efficiently sampling this distribution function. The addition of energy conservation does not change in a substantial way the predictions shown in Fig. 3. The principal consequence is that the average AMD in a single LL mode depends now, for $k\leq 4$, on the mode considered via the specific values of the frequencies $g_{k}$ and $s_{k}$. In contrast, Eq. (18) predicts an average AMD which is independent of the particular mode (see Appendix C). This slightly shifts the predicted PDFs, without changing the validity of the above considerations. This is shown in the case of Earth in Figs. 4 and 5, and in Fig. 6 for Mercury.

The major disagreement between the predictions of equiprobability in the phase-space and Laskar’s curves occurs for the eccentricity of Mercury, which is the less massive planet and has the highest typical eccentricity and inclination888This is a crucial point, as this study takes into account conservation of $L_{x}$, $L_{y}$ and $\mathcal{H}_{\mathrm{sec}}$ only at quadratic order in eccentricities and inclinations (see Eqs. (11) and (15)).. As suggested above, a kind of equiprobability ansatz could still be valuable if applied to a more restricted phase space domain than that defined by the conservation of the total inner planet AMD and energy. If, as a matter of speculation, one restricts the equiprobability domain by setting the action variable $P^{\prime}_{1}$ to its initial value, i.e.

[TABLE]

the predicted PDF of Mercury eccentricity starts to reproduce the characteristic mean value of the corresponding Laskar’s curve, as shown in Fig. 6. This is related to the fact that the chaotic dynamics of Mercury eccentricity keeps a strong memory of the dominant LL mode $P^{\prime}_{1}$999As Boué et al. (2012) and Batygin et al. (2015) show, the dynamics of Mercury treated as a test particle in the gravitational field of the other planets, whose orbits are fixed, can be described by a time-independent Hamiltonian, which is therefore an integral of motion and keeps memory of the initial conditions.. However, it is clear that a much more involved analysis, taking into account high order terms in the Hamiltonian (5), is needed to closely reproduce the numerical PDF.

6 Random walk approach to chaotic diffusion

Reproducing the PDFs of the inner planets is particularly difficult because of the time dependence of these curves. The complexity of describing diffusion of the integrals of motion in weakly nonlinear systems, starting from a given partition between linear modes, is perhaps best illustrated by the famous Fermi–Pasta–Ulam–Tsingou problem (Fermi et al., 1955). As stated in Sect. 3, in L08 the inner planet PDFs are fitted by means of a Rice distribution (Rice, 1945), whose parameters depend on time to account for the curve diffusion. The choice of this particular distribution is justified by Laskar’s suggestion that, because of chaotic diffusion, Poincaré variables (7) could behave like independent Gaussian random variables over a long timescale. In what follows I try, in the simplest way, to frame this ansatz in the context of the present paper, showing the difficulties of such a random walk approach.

In the following discussion I focus on the degrees of freedom related to eccentricity, but a similar analysis applies to the inclinations ones. According to the microcanonical distribution function (16), the PDF of the rectangular variable $h_{k}=e_{k}\sin(p_{k})$ (i.e. marginalized over the remaining Poincaré variables) can easily be shown to be

[TABLE]

where $i=\sqrt{-1}$ is the imaginary unit, $J_{0}$ the Bessel function of the first kind and order zero, and $\gamma_{kn}=\@tens{O_{A}}_{kn}(2\overline{P^{\prime}_{n}}/\Lambda_{k})^{1/2}$. Depending on the coefficients $\gamma_{kn}$, the PDF (21) can be double-peaked and strongly different from a Gaussian distribution. Moreover, since $\rho(-h_{k})=\rho(h_{k})$, the mean value of $h_{k}$ is zero. It seems doubtful whether chaotic diffusion, acting on the LL modes, can drive Poincaré variables to the Gaussianity suggested by Laskar. Therefore, I set up the random walk approach on a slight different basis. On a secular timescale $t_{\mathrm{sec}}$ of few Myr, the motion of the inner planets is well reproduced by a quasiperiodic time function. On a Lyapunov timescale, $t_{\mathrm{Lyap}}\simeq 5$ Myr, chaos starts to decorrelate the actual motion from the initial quasiperiodic approximation. On a much more longer timescale, $t\gg t_{\mathrm{Lyap}}$, the precise nature of the chaotic terms in the Hamiltonian (5) is maybe irrelevant in assessing the general form of the inner planet PDFs. Considering such a case, I propose the following ansatz, in which the stochastic dynamics of Poincaré variables $x^{\prime}_{k}$ and $y^{\prime}_{k}$ is given by the superposition of a secular term and a chaotic one modelled by means of a white noise:

[TABLE]

where $\@vec{r}=(x^{\prime}_{k},y^{\prime}_{k})$101010Since the action variables $\@vec{P^{\prime}}$ are conserved in the LL dynamics, it is reasonable to consider them as a starting point. However, as the boundary condition $P^{\prime}_{k}\geq 0$ applies, the most simple and straightforward approach is to use Poincaré variables $x^{\prime}_{k}$ and $y^{\prime}_{k}$ which are defined on the entire real line.. The stationary random variable $\@vec{r}_{\mathrm{sec}}$ describes the statistics of the LL dynamics and corresponds to the microcanonical ensemble (16),

[TABLE]

where $r^{2}={x^{\prime}_{k}}^{2}+{y^{\prime}_{k}}^{2}$. The time-dependent chaotic component is modelled by the following Langevin equation:

[TABLE]

where $\@vec{\eta}$ is a white noise, i.e.

[TABLE]

for all times $t,t^{\prime}\geq 0$. For simplicity I have assumed an isotropic diffusion in the phase space of the LL mode $k$. The diffusion coefficient $D_{k}$ is a free parameter in this model. Moreover, I choose the initial condition on the chaotic term to match the LL dynamics,

[TABLE]

Since the random variables $\@vec{r}_{\mathrm{sec}}$ and $\@vec{r}_{\mathrm{chaos}}$ are independent, the PDF of $\@vec{r}$ is given by the convolution of the PDF of $\@vec{r}_{\mathrm{sec}}$, $\rho_{\mathrm{sec}}$, and that of $\@vec{r}_{\mathrm{chaos}}$, $\rho_{\mathrm{chaos}}$. Therefore one has

[TABLE]

The PDF of the random variable satisfying Eqs. (24), (25) and (26) is the Green function of the Fokker-Planck equation

[TABLE]

This Green function is

[TABLE]

Substituting Eqs. (23) and (29) in Eq. (27) one gets

[TABLE]

where $I_{0}$ indicates the modified Bessel function of the first kind and order zero, $m=(2\overline{P^{\prime}_{k}})^{1/2}$ and $\sigma^{2}=2D_{k}t$. Therefore, one obtains a Rice distribution for the variable $r$. In the case where the variables $x_{k}$ and $y_{k}$ are dominated by one particular LL mode, i.e. $|\gamma_{kl}|\gg|\gamma_{kn}|$ for a certain index $l$ and for all $n\neq l$, the PDF of the eccentricity $e_{k}$ turns out to be a Rice distribution too, with parameters $m=|\@tens{O_{A}}_{kl}|\,(2\overline{P^{\prime}_{l}}/\Lambda_{k})^{1/2}$ and $\sigma^{2}=2D_{k}(\@tens{O_{A}}_{kl})^{2}t/\Lambda_{k}$. In this limiting case one thus recovers Laskar’s ansatz. This is due to the fact that the Rice distribution describes the envelope of a sine wave plus Gaussian noise (Rice, 1945), as already mentioned in L08. However, within the present framework, the parameter $m$ is not a free parameter as in L08, but it is set up by the initial LL dynamics.

I compare in Table 2 the value of the parameter $m$ in L08 (Table 6) for the eccentricity PDFs to the prediction of Eq. (30) when one considers for each planet only the LL mode with the biggest amplitude. The physical meaning of the parameter $m$, which is free in L08, clearly emerges in the present framework.

As shown in Table 2, the PDF (30) casts some light on Laskar’s fits in L08. Nevertheless, it cannot reproduce the numerical PDFs accurately, even if the diffusion coefficients $D_{k}$ are free parameters. There are a couple of reasons for that. Firstly, the initial distribution (23) is based on the LL dynamics. However, terms of higher order in eccentricities and inclinations in the Hamiltonian (5) are of primary importance for the inner planets. To have a faithful representation of their initial quasiperiodic dynamics one has to go beyond the LL approximation. The other principal limitation of the present framework is that PDF (30) cannot reproduce the slightly decreasing peak positions of Mars PDFs. Indeed, one has $\langle P^{\prime}_{k}\rangle=\langle r^{2}\rangle/2=m^{2}/2+\sigma^{2}=\overline{P^{\prime}_{k}}+2D_{k}t$, which means that the average content of AMD in the LL mode $k$ is an increasing function of time. Therefore, the PDFs of eccentricity and inclination can only diffuse towards increasing mean values of the corresponding random variable. The reason for this behaviour is that the random walk considered above does not take into account neither the conservation of AMD nor that of energy. Indeed, if one assumes that the AMD of the inner planets is conserved over time, it is understandable that Mercury, Earth and Venus PDFs diffuse towards increasing mean values, while Mars PDFs diffuse towards decreasing ones. Therefore, it would be a considerable improvement to set up a random walk exploiting the secular integrals of motion. Such a modification would end in a feedback dynamics between $\@vec{r}_{\mathrm{sec}}$ and $\@vec{r}_{\mathrm{chaos}}$ in Eq. (22).

7 Conclusions

Motivated by the chaotic dynamics of solar system planets and the stochastic nature of planet formation, I address in the present paper the construction of a statistical description of planetary orbits. I suggest that such an approach can be based on the statistical mechanics of secular dynamics. Considering the solar system as a benchmark, I try to reproduce the PDFs of eccentricity and inclination calculated numerically by Laskar (2008). I show that the statistics of the giant planet orbits is well reproduced by the microcanonical ensemble of the Laplace-Lagrange linear secular dynamics. This is particularly relevant as such a theory can be directly applied to a generic planetary system. Minor discrepancies in this description are connected to higher-than-quadratic terms in the planet Hamiltonian. Their contribution could be taken into account by improving the perturbation approach via successive quasi-identity canonical transformation and normal forms. I then try to reproduce the statistics of the inner planet orbits through the ansatz of equiprobability in the phase-space constrained by the secular integrals of motion, namely angular momentum and energy. Within the limitations of a stationary model, such an ansatz allows to easily and accurately reproduce the structure of Venus and Earth PDFs without any free parameter. However, major discrepancies rise in the case of Mercury eccentricity. I finally show the difficulties of constructing a random walk to model the chaotic diffusion of the inner planet PDFs, following the original ansatz of Laskar (2008). Within certain limitations, the PDF I obtain allows to illustrate the physical meaning of one of the free parameters in Laskar’s ansatz.

It is clear that the above description of the inner planet statistics has to be improved. One has to taken into account higher-than-quadratic terms in the planet Hamiltonian. Generally speaking, one could think of using the full expression of the secular energy, Eq. (5). The problem would then be how efficiently sample the corresponding microcanonical ensemble. Another viable approach could be to start constructing an ad hoc model for Mercury, since it shows the most relevant discrepancies with respect to the proposed ansatz. Such an approach could consist in considering Mercury as a test particle in the gravitational field of the other solar system planets, whose orbits are fixed (Boué et al., 2012; Batygin et al., 2015). With a selection of the most important higher order terms relevant for the non-linear dynamics of Mercury, important improvements could be achieved in reproducing the PDFs of its eccentricity and inclination. With respect to the random walk approach presented in Sect. 6, a considerable improvement would be taking into account the conservation of angular momentum deficit and energy. Even if such a model is not completely predictive, as the diffusion in the phase space is described by free parameters, it could be nevertheless really instructive in judging how the conservation of the integrals of motion constrains the relaxation of the Laplace-Lagrange action variables.

The relevance of a statistical description of planet orbits in a generic planetary system is particularly manifest if one considers the model of planet spacing by Laskar (2000) (see also Laskar & Petit, 2017). In this approach to planet formation, Laskar describes the collisionless dynamics of planetary embryos by a stochastic variation in their orbital elements that conserves the total angular momentum deficit. Conservation of mass and linear momentum is taken into account to model the result of collisions, during which the total angular momentum deficit decreases. Collisions stop when the total angular momentum deficit is too smalls to allow for further close encounters. With these simple assumptions, Laskar is able to show how the structure of a generic planetary system derives from the initial mass distribution of planetary embryos. As already pointed out by Tremaine (2015), Laskar’s model is not unique because it requires an ad hoc prescription for the random evolution of the orbital elements between collisions. Statistical descriptions like those in Sects. 4 and 5 can provide such a prescription based on solid physical arguments. A statistical approach to planetary orbits could then be useful in studying planet formation, in combination with standard N-body simulation of planetary trajectories.

Acknowledgements.

I am grateful to A. Morbidelli, J. Touma, J. Laskar, J. Teyssandier, N. Cornuault and L. Pittau for the fruitful discussions. I would also like to thank J.-P. Beaulieu for his support and help in reviewing the manuscript.

Appendix A General relativistic correction

I present a short derivation of the leading general relativistic contribution to the secular planetary Hamiltonian. At the order $(v/c)^{2}$ the Lagrangian of Sun and planets is given by (Landau & Lifshitz, 1975)

[TABLE]

where $L_{0}=\sum_{a}m_{a}v^{2}_{a}/2+\sum_{a}{\sum}^{\prime}_{b}Gm_{a}m_{b}/2r_{ab}$ is the Newtonian Lagrangian. In the equation above $a,b,c\in\{0,1,\dots,8\}$, $r_{ab}=|\@vec{r}_{a}-\@vec{r}_{b}|$, $n_{ab}=(\@vec{r}_{a}-\@vec{r}_{b})/r_{ab}$ and the prime symbol means that the terms $b=a$ and $c=a$ must be excluded from the summation. I also recall that $m_{i}/m_{0}\ll 1$ for $i\in\{1,2,\dots,8\}$. I then simplify Eq. (31) by keeping among the relativistic terms only those of order $(v/c)^{2}$. I thus neglect terms of order $(m_{i}/m_{0})(v/c)^{2}$ and smaller. One obtains $L=L_{0}+\delta L$, with

[TABLE]

where $r_{i}$ is the Heliocentric position of planet $i$. The general relativistic correction to the planetary Hamiltonian is thus given by $\delta\mathcal{H}=-\delta L$ (Landau & Lifshitz, 1969a). To obtain the secular correction one has to average $\delta\mathcal{H}$ over the Keplerian orbits of the non-interacting planets. Neglecting corrections of the order $m_{i}/m_{0}$, along a Keplerian orbit one has

[TABLE]

where $a_{i}$ is the semi-major axis of planet $i$. Moreover, averaging over an orbital period one obtains $\langle r_{i}^{-1}\rangle=a_{i}^{-1}$ and $\langle r_{i}^{-2}\rangle=a_{i}^{-2}(1-e^{2}_{i})^{-1/2}$, with $e_{i}$ the eccentricity of planet $i$. Substituting Eq. (33) in $\delta\mathcal{H}$ and taking the average, one finally finds

[TABLE]

Discarding the first term as it is constant in the secular dynamics, one obtains that the dominant secular contribution of general relativity is $\mathcal{H}_{\mathrm{sec,\,GR}}=-3G^{2}m^{2}_{0}m_{i}/c^{2}a^{2}_{i}(1-e^{2}_{i})^{1/2}$.

Appendix B Laplace-Lagrange Hamiltonian

The matrix $\@tens{A}$ appearing in the quadratic Hamiltonian (6) is given by the following expression:

[TABLE]

where I have defined

[TABLE]

Similarly, the matrix $\@tens{B}$ is given by

[TABLE]

From formula (36) it is easy to check that matrix $\@tens{B}$ verifies the following relation:

[TABLE]

This implies that one of the $s_{k}$ eigenvalues is zero and $\left.\left(\!\!\sqrt{\Lambda_{1}},\dots,\!\!\sqrt{\Lambda_{N}}\right)\,\right/\!\!\sqrt{\sum_{k}\Lambda_{k}}$ is, up to a sign, the corresponding normalized eigenvector.

Appendix C Ansatz of equiprobability: analytical description

An analytical characterization of the PDFs of eccentricity and inclination predicted by Eq. (18) can be obtained by writing the conservation of the total AMD of the inner planets, as motivated in Sect. 5, through the modified Delauney variables (2). Compared to that of Sect. 5, such an approach does not include the statistics of giant planet orbits. The microcanonical density of states $\omega$ then reads

[TABLE]

where $N=4$ in the present analysis, $C_{0}$ is the initial total AMD of the inner planets and the factor $(2\pi)^{2N}$ comes from the integration over the angle variables $p_{k}$ and $q_{k}$. Since $\Lambda_{k}\geq 3C_{0}/2$ for all $k$ (Laskar, 1997, Table 1), one can restrict the above integration to the hypercube $[0,C_{0}]^{2N}$. To calculate the density of states one can therefore write $\omega(C_{0})=d\Sigma(C_{0})/dC_{0}$, with

[TABLE]

where $H$ is the Heaviside step function. Performing the integrations by iteration, one obtains $\Sigma(C_{0})=(2\pi)^{2N}{C_{0}}^{2N}/2N!$ and then

[TABLE]

The joint microcanonical PDF of $P_{k}$ and $Q_{k}$ (i.e. marginalized over the remaining modified Delauney variables) is then given by

[TABLE]

for $P_{k}+Q_{k}\leq C_{0}$, otherwise it is zero. Thus one obtains

[TABLE]

for $P_{k}+Q_{k}\leq C_{0}$. Switching to the variables $e_{k}$ and $i_{k}$, one finds

[TABLE]

for $\eta^{-1}\left(1-\sqrt{1-e_{k}^{2}}\cos i_{k}\right)\leq 1$, with $\eta=C_{0}/\Lambda_{k}$. The PDF of $e_{k}$ is then given by

[TABLE]

where $\cos i^{\star}=(1-\eta)/(1-e_{k}^{2})^{1/2}$. Performing the integral, one obtains

[TABLE]

for $e_{k}\leq\sqrt{1-(1-\eta)^{2}}$. Similarly, to obtain the PDF of $i_{k}$ one has to calculate

[TABLE]

where $e^{\star}=[1-(1-\eta)^{2}/\cos^{2}i_{k}]^{1/2}$. Performing the integration, one finds

[TABLE]

for $\cos i_{k}\geq 1-\eta$. As in the present study one has $\eta\ll 1$ (Laskar, 1997, Table 1), one can simplify these PDFs obtaining

[TABLE]

In this limit one has that

[TABLE]

where $\Gamma$ is the Gamma function, and the average AMD of planet $k$ is thus given by

[TABLE]

Therefore one obtains the equipartition of AMD between planets and between the eccentricity and inclination degrees of freedom. Even though they do not coincide, the PDFs (48) and (49) present the same overall structure as those predicted by the ansatz (18).

One can also include in the above analytical treatment the conservation of the angular momentum components $L_{x}$ and $L_{y}$ at quadratic order in eccentricities and inclinations (see Eq. (14)). While the shape of the eccentricity PDF is unaffected, that of the inclination PDF turns out to be somewhat influenced by the parameter $(L_{x}^{2}+L_{y}^{2})^{1/2}/L_{z}$.

Appendix D Ansatz of equiprobability: sampling

Conservation of AMD and energy in Eqs. (18) and (19) is equivalent to the two following conditions:

[TABLE]

with $C_{0}=\sum_{k=1}^{4}\overline{P^{\prime}_{k}}+\overline{Q^{\prime}_{k}}$ and $H_{0}=-\sum_{k=1}^{4}\left(g_{k}\overline{P^{\prime}_{k}}+s_{k}\overline{Q^{\prime}_{k}}\right)$. These equations can be rewritten as

[TABLE]

where $a_{k}=P^{\prime}_{k}/C_{0}$, $b_{k}=Q^{\prime}_{k}/C_{0}$ and $h=-H_{0}/C_{0}$. Moreover, I recall that $a_{k},b_{k}\geq 0$, $g_{k}>0$ and $s_{k}<0$ . The PDF (18) can be therefore evaluated by an uniform sampling of $a_{k}$ and $b_{k}$ verifying Eq. (55). To evaluate the PDF (19) $a_{k}$ and $b_{k}$ have to verify both Eq. (55) and Eq. (56).

Uniform sampling of $N$ positive real numbers that add up to unity, as in Eq. (55), can be performed via a direct sampling algorithm. One starts by sampling $N-1$ random numbers uniformly from the interval $[0,1]$. After sorting them, one obtains the set $\{u_{1},\dots,u_{N-1}\}$, with $u_{i-1}\leq u_{i}$ for $i\in\{2,\dots,N-1\}$. Then, the numbers $\{u_{1},u_{2}-u_{1},\dots,u_{N-1}-u_{N-2},1-u_{N-1}\}$ follow the target distribution.

To uniformly sampling positive real numbers $a_{k}$ and $b_{k}$ from Eqs. (55) and (56) there is probably no direct sampling method. However, I present here an efficient algorithm with rejection. Let’s assume that $h\geq 0$, as in the present study. A similar algorithm can be constructed in case $h$ is negative. One can multiply Eq. (55) by $g_{\mathrm{max}}=\max_{k}\{a_{k}\}$, subtract Eq. (56) and finally divide by $g_{\mathrm{max}}-h$ (one can assume this to be different from zero, otherwise the sampling is trivial). One obtains

[TABLE]

The coefficients of $a_{k}$ and $b_{k}$ in Eq. (57) are positive and the term $a_{k^{\star}}$, such that $g_{k^{\star}}=g_{\mathrm{max}}$, does not appear in it. Therefore, the other $N-1$ numbers can be sampled from this equation using the direct sampling algorithm presented above. If their sum $S$ turns out to be smaller or equal to unity, then one can take $a_{k^{\star}}=1-S$ to obtain a set of numbers following the target distribution. Otherwise, the numbers are rejected and Eq. (57) sampled again.