© ESO, 2017

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. Laskar 1997; Malhotra 1999; 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 of the sheared-sheet approximation and the restriction to the planar case, this ansatz permits the analytical computation of 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 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 mergings. However, one could ask for a more fundamental approach, exploiting for instance 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 actually be set up by focusing on how eccentricities and inclinations are excited by gravitational interactions in an early, collisionless final phase of terrestrial planet formation and by 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 box to assure a bounded phase space. 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 one 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 ≫ 1 regime is generally inappropriate and so is the employment of a mean field approach¹. 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 (Gauss 1818; Hill 1882; Touma et al. 2009; Kondratyev 2012). In such a system, the rings interact with each other via secular gravitational interactions. Their eccentricities and inclinations relax via exchange of angular momentum, while their semimajor axes are constant and hence are their Keplerian energies (Rauch & Tremaine 1996). Even if there is no theorem assuring that this secular dynamics approaches the actual planet motion, this approach can nevertheless represent the source of important results (Arnold 1978). Indeed, the first clear indication 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 by emphasizing the relevance of a statistical theory of planetary orbits in the context of planet formation.

2. Secular planet dynamics

Using standard notation (e.g. Morbidelli 2002), the Hamiltonian of the N = 8 solar system planets can be written as $\mathrm{\mathscr{H}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\left(\frac{{{p}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathit{k}}}\mathrm{-}\mathit{G}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{k}}}{{\mathit{r}}_{\mathit{k}}}\right)\mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathit{k}\mathrm{+}\mathrm{1}}^{\mathit{N}}\left(\frac{{{p}}_{\mathit{k}}\mathrm{·}{{p}}_{\mathit{l}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{-}\mathit{G}\frac{{\mathit{m}}_{\mathit{k}}{\mathit{m}}_{\mathit{l}}}{\mathrm{|}{{r}}_{\mathit{k}}\mathrm{-}{{r}}_{\mathit{l}}\mathrm{|}}\right)\mathit{,}$(1)where heliocentric canonical variables r_k,p_k are employed, m₀ is the Sun mass, m_k the planet masses, μ_k = m₀m_k/ (m₀ + m_k) the reduced masses, and G the gravitational constant. As usual, the index k lists the planets by increasing semimajor axis, from Mercury to Neptune. The first term on the right side of Eq. (1)is the Hamiltonian of the Keplerian motion, while the second term 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), as follows: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{\Lambda }}_{\mathit{k}}\mathrm{=}{\mathit{\mu }}_{\mathit{k}}\sqrt{\mathit{G}\mathrm{(}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathit{k}}\mathrm{)}{\mathit{a}}_{\mathit{k}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{P}}_{\mathit{k}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{k}}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\mathrm{)}\\ \begin{array}{c}\end{array}{\mathit{Q}}_{\mathit{k}}\mathrm{=}\mathrm{2}{\mathrm{\Lambda }}_{\mathit{k}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}{\sin }^{\mathrm{2}}\mathrm{(}{\mathit{i}}_{\mathit{k}}\mathit{/}\mathrm{2}\mathrm{)}\end{array}& & \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\lambda }}_{\mathit{k}}\mathrm{=}{\mathit{M}}_{\mathit{k}}\mathrm{+}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathrm{\Omega }k\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{p}}_{\mathit{k}}\mathrm{=}\mathrm{-}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathrm{\Omega }k\\ \begin{array}{c}\end{array}{\mathit{q}}_{\mathit{k}}\mathrm{=}\mathrm{-}{\mathrm{\Omega }}_{\mathit{k}}\mathit{.}\end{array}\end{array}$(2)Keplerian orbital elements are used in these definitions: a_k is the planet orbit semimajor axis, e_k the eccentricity, i_k the inclination, M_k the mean anomaly, ω_k the argument of perihelion, and Ω_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, $\mathrm{⟨}\mathrm{\mathscr{H}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{N}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}{\mathit{M}}_{\mathrm{1}}\hspace{0.17em}\mathrm{·}\mathrm{·}\mathrm{·}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}{\mathit{M}}_{\mathit{N}}\hspace{0.17em}\hspace{0.17em}\mathrm{\mathscr{H}}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{+}{\mathrm{\mathscr{H}}}_{\sec }\mathit{,}$(3)where ${\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{-}\mathit{G}\mathrm{(}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}{\mathit{\mu }}_{\mathit{k}}^{\mathrm{3}}\mathit{/}\mathrm{(}\mathrm{2}{\mathrm{\Lambda }}_{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{-}\mathit{G}{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{k}}\mathit{/}\mathrm{(}\mathrm{2}{\mathit{a}}_{\mathit{k}}\mathrm{)}$ is the total Keplerian energy and ${\mathrm{\mathscr{H}}}_{\sec }\mathrm{=}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathit{k}\mathrm{+}\mathrm{1}}^{\mathit{N}}\frac{\mathit{G}{\mathit{m}}_{\mathit{k}}{\mathit{m}}_{\mathit{l}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}{\mathit{M}}_{\mathit{k}}\mathrm{d}{\mathit{M}}_{\mathit{l}}}{\mathrm{|}{{r}}_{\mathit{k}}\mathrm{-}{{r}}_{\mathit{l}}\mathrm{|}}$(4)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 (e.g. Morbidelli 2002). The action variables Λ_k are integrals of motion of the secular dynamics (Arnold 1978), and so are the semimajor axes a_k. As a consequence, ℋ₀ is an additive constant that 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_kdq_kdP_kdQ_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 ${\mathrm{\mathscr{H}}}_{\sec }\mathrm{=}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathit{k}\mathrm{+}\mathrm{1}}^{\mathit{N}}\frac{\mathit{G}{\mathit{m}}_{\mathit{k}}{\mathit{m}}_{\mathit{l}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}{\mathit{M}}_{\mathit{k}}\mathrm{d}{\mathit{M}}_{\mathit{l}}}{\mathrm{|}{{r}}_{\mathit{k}}\mathrm{-}{{r}}_{\mathit{l}}\mathrm{|}}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{3}{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{m}}_{\mathit{k}}}{{\mathit{c}}^{\mathrm{2}}{\mathit{a}}_{\mathit{k}}^{\mathrm{2}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}}\mathit{,}$(5)where c is the speed of light. A short derivation of the general relativistic contribution is presented in Appendix A for future reference.

When eccentricities and inclinations are sufficiently small, e_k,sin(i_k/ 2) ≪ 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, ${\mathrm{\mathscr{H}}}_{\mathrm{LL}}\mathrm{=}\mathrm{-}\frac{{{x}}^{\mathrm{T}}\mathrm{A}{x}\mathrm{+}{{y}}^{\mathrm{T}}\mathrm{A}{y}\mathrm{+}{{v}}^{\mathrm{T}}\mathrm{B}{v}\mathrm{+}{{z}}^{\mathrm{T}}\mathrm{B}{z}}{\mathrm{2}}\mathit{,}$(6)where I have introduced the Poincaré canonical variables at the lowest order in eccentricity², $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{x}}_{\mathit{k}}\mathrm{=}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}{\mathit{e}}_{\mathit{k}}\sin \mathit{p}k\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{y}}_{\mathit{k}}\mathrm{=}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}{\mathit{e}}_{\mathit{k}}\cos \mathit{p}k\end{array}& & \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{v}}_{\mathit{k}}\mathrm{=}\mathrm{2}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}\sin \mathrm{(}{\mathit{i}}_{\mathit{k}}\mathit{/}\mathrm{2}\mathrm{)}\sin \mathit{q}k\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{z}}_{\mathit{k}}\mathrm{=}\mathrm{2}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}\sin \mathrm{(}{\mathit{i}}_{\mathit{k}}\mathit{/}\mathrm{2}\mathrm{)}\cos \mathit{q}k\end{array}\end{array}$(7)and for instance x^T = (x₁,...,x_N). The elements of the matrices A and B are given in Appendix B. These matrices are real and symmetric and can therefore be diagonalized through orthogonal matrices O_A and O_B, $\mathrm{A}\mathrm{=}{\mathrm{O}}_{\mathrm{A}}{\mathrm{D}}_{\mathrm{A}}{\mathrm{O}}_{\mathrm{A}}^{\mathrm{T}}\mathit{,} \mathrm{B}\mathrm{=}{\mathrm{O}}_{\mathrm{B}}{\mathrm{D}}_{\mathrm{B}}{\mathrm{O}}_{\mathrm{B}}^{\mathrm{T}}\mathit{.}$(8)The elements of the diagonal matrices D_A and D_B are the eigenvalues of A and B, respectively, and the columns of O_A and O_B are their normalized eigenvectors. New action-angle coordinates P^′,Q^′,p^′,q^′ can be introduced by combining the following canonical transformations: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{x}\mathrm{=}{\mathrm{O}}_{\mathrm{A}}{{x}}^{{\prime }}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{y}\mathrm{=}{\mathrm{O}}_{\mathrm{A}}{{y}}^{{\prime }}\end{array}& & \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{v}\mathrm{=}{\mathrm{O}}_{\mathrm{B}}{{v}}^{{\prime }}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{z}\mathrm{=}{\mathrm{O}}_{\mathrm{B}}{{z}}^{{\prime }}\end{array}\end{array}$(9)$\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{2}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\sin {\mathit{p}}_{\mathit{k}}^{\mathrm{\prime }}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{2}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\cos {\mathit{p}}_{\mathit{k}}^{\mathrm{\prime }}\end{array}& & \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{v}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{2}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}\sin {\mathit{q}}_{\mathit{k}}^{\mathrm{\prime }}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{z}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{2}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}\cos {\mathit{q}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{.}\end{array}\end{array}$(10)Employing the new variables, the LL Hamiltonian becomes${\mathrm{\mathscr{H}}}_{\mathrm{LL}}\mathrm{=}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{g}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}{\mathit{s}}_{\mathit{k}}{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{,}$(11)where g_k and s_k are the eigenvalues of the matrices A and B, respectively. According to the corresponding Hamilton equations, the actions ${\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}$ and ${\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}$ are integrals of motion, while the angles ${\mathit{p}}_{\mathit{k}}^{\mathrm{\prime }}$ and ${\mathit{q}}_{\mathit{k}}^{\mathrm{\prime }}$ 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).

Generally speaking, the dynamics of a planetary system obeys the conservation of mechanical energy, momentum, and angular momentum. The secular dynamics verifies the conservation of momentum in the barycentric reference system in a trivial way: 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 LL dynamics (e.g. Kocsis & Tremaine 2011). As formula (11)shows, the LL Hamiltonian is a function of the action variables only. Because these are integrals of motion, the conservation of energy follows. The angular momentum content of a planetary system can be described by the following quantities: $\begin{array}{ccc}{\mathit{L}}_{\mathit{x}}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Lambda }}_{\mathit{k}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\sin {\mathit{i}}_{\mathit{k}}\sin {\mathrm{\Omega }}_{\mathit{k}}\\ {\mathit{L}}_{\mathit{y}}& \mathrm{=}& \mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Lambda }}_{\mathit{k}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\sin {\mathit{i}}_{\mathit{k}}\cos {\mathrm{\Omega }}_{\mathit{k}}\\ \mathit{C}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Lambda }}_{\mathit{k}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}\right)\mathit{,}\end{array}$(12)where x,y,z are the Cartesian coordinates related to the definition of the Keplerian elements i_k,ω_k, and Ω_k (Morbidelli 2002)³. The variables L_x and L_y are the components of the angular moment lying on the xy-reference plane, while C is the angular momentum deficit (AMD) (Laskar 1997, 2000), i.e. 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 $\mathit{C}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{Q}}_{\mathit{k}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{,}$(13)where I have used the definitions (2)and the fact that the matrices O_A and O_B, appearing in the transformation (9), are orthogonal. As the action variables P^′ and Q^′ are conserved in LL dynamics, the AMD is also a conserved quantity⁴. Conservation of L_x and L_y derives from the properties of the matrix B. As shown in Appendix B, one of the eigenvalues s_k is zero and ${}^{\mathrm{(}}\sqrt{{\mathrm{\Lambda }}_{\mathrm{1}}}\mathit{,}\mathit{...}\mathit{,}{\sqrt{{\mathrm{\Lambda }}_{\mathit{N}}}}^{\mathrm{)}}{\hspace{0.17em}}^{\mathrm{􏼮}}\sqrt{{\sum }_{\mathit{k}}{\mathrm{\Lambda }}_{\mathit{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₅. Moreover, neglecting terms of order three in eccentricities and inclinations, one has $\begin{array}{ccc}{\mathit{L}}_{\mathit{x}}& \mathrm{=}& \mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{2}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}{\mathit{v}}_{\mathit{k}}\\ {\mathit{L}}_{\mathit{y}}& \mathrm{=}& \mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{2}\sqrt{{\mathrm{\Lambda }}_{\mathit{k}}}{\mathit{z}}_{\mathit{k}}\mathit{.}\end{array}$(14)All this implies that in LL dynamics the following identities are verified: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{5}}^{\mathrm{\prime }}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{{\mathrm{O}}_{\mathrm{B}}^{\mathrm{T}}}_{\mathrm{5}\mathit{k}}{\mathit{v}}_{\mathit{k}}\mathrm{=}\mathrm{-}\frac{{\mathit{L}}_{\mathit{x}}}{\mathrm{2}\sqrt{{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Lambda }}_{\mathit{k}}}}& & \\ {\mathit{z}}_{\mathrm{5}}^{\mathrm{\prime }}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{{\mathrm{O}}_{\mathrm{B}}^{\mathrm{T}}}_{\mathrm{5}\mathit{k}}{\mathit{z}}_{\mathit{k}}\mathrm{=}\mathrm{-}\frac{{\mathit{L}}_{\mathit{y}}}{\mathrm{2}\sqrt{{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Lambda }}_{\mathit{k}}}}\mathit{,}& & \end{array}$(15)where I have used the fact that O_B is orthogonal, ${\mathrm{O}}_{\mathrm{B}}^{-1}\mathrm{=}{\mathrm{O}}_{\mathrm{B}}^{\mathrm{T}}$, and its columns are the normalized eigenvectors of B. As the frequency s₅ is zero, ${\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}$ does not change with time. Therefore, ${\mathit{v}}_{\mathrm{5}}^{\mathrm{\prime }}$ and ${\mathit{z}}_{\mathrm{5}}^{\mathrm{\prime }}$ are integrals of motion and hence are L_x and L_y.

Fig. 1 Eccentricity and inclination PDFs of the solar system giant planets. The PDFs of L08 are plotted in solid lines, while those predicted by the microcanonical distribution function (16)are shown in dashed lines. The inclinations are computed with respect to the J2000 ecliptic.

3. Probability density functions 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 remaining lifetime of the Sun 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 time step. 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 those of the inner planets are shown. Unfortunately, the data behind the L08 PDFs are not currently available (Laskar, priv. comm.). 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 planets, 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 dynamics 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 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 quasi-periodic approximation of their dynamics can reproduce the numerical PDFs very accurately. 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 a non-zero mean and a variance that increases linearly with time, as in a diffusion process.

Fig. 2 Eccentricity and inclination PDFs of the solar system inner planets. The PDFs of L08 are plotted in solid lines at three different times, t = 500 Myr, 2.5 Gyr, and 5 Gyr, to illustrate chaotic diffusion (the PDF broadens and its peak decreases as time increases; the inclination PDFs are not the original PDFs, but the Rice distributions fitted by Laskar, see Sect. 3 for explanations). The PDFs predicted by the microcanonical ensemble of the Laplace-Lagrange theory, Eq. (16), are shown in dashed line. The inclinations are computed with respect to the J2000 ecliptic.

4. Giant planets

According to secular dynamics, the orbital motion of the giant planets is well approximated by a quasi-periodic 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 to 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 quasi-periodic dynamics is that resulting from the LL Hamiltonian, Eq. (11). Such a motion is ergodic and the asymptotic stationary probability density is given by the microcanonical ensemble (Arnold 1978). It is straightforward to write down the corresponding PDF employing the integrals of motion (Landau & Lifshitz 1969b), $\mathit{\rho }\mathrm{\left(}{{P}}^{{\prime }}\mathit{,}{{Q}}^{{\prime }}\mathit{,}{{p}}^{{\prime }}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta }\mathrm{(}{\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}}\mathrm{)}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}}\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\mathit{\delta }\mathrm{(}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{)}\hspace{0.17em}\mathit{\delta }\mathrm{(}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{)}\mathit{,}$(16)where δ stands for the Dirac delta and the bar indicates the initial value of the corresponding variable⁵. This PDF reflects the conservation of the action variables P^′,Q^′ in the LL dynamics, while the angle variables are uniformly distributed over the (2N−1)-dimensional torus T^2N−1; the angle ${\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}$ is conserved because the angular momentum is conserved (see Sect. 2). The following integral formula for the PDFs of eccentricity and inclination can be obtained (Mardia 1972): $\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathit{k}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{u}\hspace{0.17em}\hspace{0.17em}\mathit{u}\hspace{0.17em}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathit{u}\mathrm{\right)}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}{\mathit{\gamma }}_{\mathit{kn}}\mathit{u}\mathrm{\right)}\mathit{,}$(17)where J₀ is the Bessel function of the first kind and order zero and ${\mathit{\gamma }}_{\mathit{kn}}\mathrm{=}{{\mathrm{O}}_{\mathrm{A}}}_{\mathit{kn}}\mathrm{(}\mathrm{2}\overline{{\mathit{P}}_{\mathit{n}}^{\mathrm{\prime }}}\mathit{/}{\mathrm{\Lambda }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{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 (dozens of kyr to few Myr for the solar system planets). Since the PDFs of Laskar 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 the PDFs found by Laskar. The agreement between the two curves is indeed very satisfactory for the inclination PDFs, with a very good match 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 small. However, the terms of the fourth order in e and i in the Hamiltonian (5)yield non-linearities in the corresponding Hamilton equations and produce higher order harmonics in the eccentricity and inclination time dependence⁶. The amplitude of these additional harmonics is generally much smaller than that of the LL modes, but some of the higher order harmonics can nevertheless contribute to a significant level. Bretagnon (1974, Tables 12 and 13) reports the amplitudes of these higher order harmonics. From the Bretagnon tables it is clear that their contribution is bigger for eccentricities than inclinations and for Uranus and Neptune than for 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 the PDFs found by Laskar 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 carried out by means of a perturbation approach to the Hamiltonian (1)based on successive quasi-identity canonical transformations and normal forms (Morbidelli 2002)⁷.

Fig. 3 Eccentricity and inclination PDFs of the solar system inner planets. The same as Fig. 2, except that the dashed lines represent here the predictions of the ansatz (18).

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 the orbits of Mercury and Mars 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 the PDFs found by Laskar in case of inner planets is to consider what the distribution function (16)predicts for their eccentricities and inclinations. 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 general agreement between the two curves. Nevertheless, it is interesting to note that, contrary to what could emerge from the analysis by Laskar in L08, the difference in shape between the giant and inner planet PDFs is not fundamental. Figure 2 shows that Eq. (16)based on 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, which justifies the strong differences between the corresponding PDFs found by Laskar and Eq. (16)in Fig. 2.

Since the effect of chaos is relevant for the long-term dynamics of the inner planets, I searched for a more pertinent statistical description than one based on quasi-periodic motion. The chaotic nature of planet dynamics allows for sharing of angular momentum between the 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 to assume that, by consequence, the motion of the inner planets uniformly visits the subset of the phase space defined by the conservation of angular momentum L and energy ℋ_sec. Clearly, the fact that chaotic diffusion finally leads to such a uniform phase-space measure cannot be strictly true over 5 Gyr, as the time dependence of the inner planet PDFs is still evident at that time. More probably, a 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). The hypothesis of equiprobability in the phase space is the best estimate that one can make without a deeper analysis of the dynamics (Jaynes 1957). It is also 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 seems nevertheless appropriate to the 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, it turns out that the total angular momentum of the inner planets is fairly independent of that of the giant planets. More precisely, Laskar (1997) has shown that the flux of AMD from the large reservoir of the outer planets to the inner planets 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 ${\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}$ and ${\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}$ as the AMD content in the corresponding LL mode. These considerations suggest fixing the action variables corresponding to the modes k ∈ { 5,6,7,8 }, as these are the only modes that are relevant for the dynamics of the outer planets⁸. Indeed, with this choice the conservation of the total AMD reduces to that of the quantity ${\sum }_{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{4}}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}$, 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 ∈ { 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 that 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, contrary to 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 is added later. The ansatz corresponding to the above considerations is written $\begin{array}{ccc}& & \mathit{\rho }\mathrm{\left(}{{P}}^{{\prime }}\mathit{,}{{Q}}^{{\prime }}\mathit{,}{{p}}^{{\prime }}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{\propto }\mathit{\delta }\mathrm{(}\mathit{C}\mathrm{-}\overline{\mathit{C}}\mathrm{)}\hspace{0.17em}\mathit{f}\\ & & \mathit{f}\mathrm{\left(}{{P}}^{{\prime }}\mathit{,}{{Q}}^{{\prime }}\mathit{,}{{p}}^{{\prime }}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{=}\mathit{\delta }\mathrm{(}{\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}}\mathrm{)}\underset{\mathit{k}\mathrm{=}\mathrm{5}}{\overset{\mathrm{8}}{\mathrm{􏽙}}}\mathit{\delta }\mathrm{(}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{)}\hspace{0.17em}\mathit{\delta }\mathrm{(}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{)}\mathit{.}\end{array}$(18)Conservation of L_x and L_y to second order in eccentricities and inclinations is contained in this distribution function as the variables ${\mathit{Q}}_{\mathrm{5}}^{\mathrm{\prime }}$ and ${\mathit{q}}_{\mathrm{5}}^{\mathrm{\prime }}$ 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 Appendix D. Moreover, I present an analytical characterization of these PDFs in Appendix C. The predictions of the ansatz (18)are compared in Fig. 3 with the PDFs found by Laskar. The probability density (18)reproduces Venus and Earth PDFs very satisfactorily. 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 parameter⁹. For what concerns the PDFs of Mars and that of the inclination of Mercury, 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 the PDFs found by Laskar. 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 PDFs. In particular, the mean of the curves found by Laskar 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.

Fig. 4 Eccentricity PDF of the Earth. The same as Fig. 3 with the addition of the prediction of Eq. (19)in dash-dotted line.

One can add to the ansatz (18) the conservation of secular energy at quadratic order in eccentricities and inclinations (see Eq. (11)), $\mathit{\rho }\mathrm{\left(}{{P}}^{{\prime }}\mathit{,}{{Q}}^{{\prime }}\mathit{,}{{p}}^{{\prime }}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{\propto }\mathit{\delta }\mathrm{(}\mathit{C}\mathrm{-}\overline{\mathit{C}}\mathrm{)}\hspace{0.17em}\mathit{\delta }\mathrm{(}{\mathrm{\mathscr{H}}}_{\mathrm{LL}}\mathrm{-}\overline{{\mathrm{\mathscr{H}}}_{\mathrm{LL}}}\mathrm{)}\hspace{0.17em}\mathit{f}\mathit{.}$(19)In 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 ≤ 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 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.

Fig. 5 Inclination PDF of the Earth. The same as Fig. 3 with the addition of the prediction of Eq. (19)in dash-dotted line.

Fig. 6 Eccentricity PDF of Mercury. The same as Fig. 3 with the addition of the predictions of Eq. (19)in dash-dotted line and Eq. (20)in dotted line.

The major disagreement between the predictions of equiprobability in the phase space and the curves of Laskar occurs for the eccentricity of Mercury, which is the less massive planet and has the highest typical eccentricity and inclination¹⁰. 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 ${\mathit{P}}_{\mathrm{1}}^{\mathrm{\prime }}$ to its initial value, i.e. $\mathit{\rho }\mathrm{\left(}{{P}}^{{\prime }}\mathit{,}{{Q}}^{{\prime }}\mathit{,}{{p}}^{{\prime }}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{\propto }\mathit{\delta }\mathrm{(}\mathit{C}\mathrm{-}\overline{\mathit{C}}\mathrm{)}\hspace{0.17em}\mathit{\delta }\mathrm{(}{\mathrm{\mathscr{H}}}_{\mathrm{LL}}\mathrm{-}\overline{{\mathrm{\mathscr{H}}}_{\mathrm{LL}}}\mathrm{)}\hspace{0.17em}\mathit{\delta }\mathrm{(}{\mathit{P}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{-}\overline{{\mathit{P}}_{\mathrm{1}}^{\mathrm{\prime }}}\mathrm{)}\hspace{0.17em}\mathit{f,}$(20)the predicted PDF of Mercury eccentricity starts to reproduce the characteristic mean value of the corresponding curves of Laskar, as shown in Fig. 6. This is related to the fact that the chaotic dynamics of Mercury eccentricity keeps a clear memory of the dominant LL mode ${\mathit{P}}_{\mathrm{1}}^{\mathrm{\prime }}$¹¹. However, it is clear that more involved analyses, taking into account higher order terms in the Hamiltonian (5), are needed to reproduce the numerical PDFs closely.

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 non-linear 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. Laskar justifies the choice of this particular distribution by stating that, because of chaotic diffusion, Poincaré variables (7)could behave similar to 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 inclination degrees of freedom. According to the microcanonical distribution function (16), the PDF of the rectangular variable h_k = e_ksin(p_k) (i.e. marginalized over the remaining Poincaré variables) can be shown to be¹²$\mathit{\rho }\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{u}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{h}}_{\mathit{k}}\mathit{u}}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathit{J}}_{\mathrm{0}}\left({\mathit{\gamma }}_{\mathit{kn}}\mathit{u}\right)\mathit{,}$(21)where $\mathit{i}\mathrm{=}\sqrt{-1}$ is the imaginary unit, J₀ the Bessel function of the first kind and order zero, and ${\mathit{\gamma }}_{\mathit{kn}}\mathrm{=}{{\mathrm{O}}_{\mathrm{A}}}_{\mathit{kn}}\mathrm{(}\mathrm{2}\overline{{\mathit{P}}_{\mathit{n}}^{\mathrm{\prime }}}\mathit{/}{\mathrm{\Lambda }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. Depending on the coefficients γ_kn, the PDF (21)can be double-peaked and strongly different from a Gaussian distribution. Moreover, since ρ(−h_k) = ρ(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_sec of few Myr, the motion of the inner planets is well reproduced by a quasi-periodic time function. On a Lyapunov timescale, t_Lyap ≃ 5 Myr, chaos starts to decorrelate the actual motion from the initial quasi-periodic approximation. On a much more longer timescale, t ≫ t_Lyap, the precise nature of the chaotic terms in the Hamiltonian (5)is possibly irrelevant in assessing the general form of the inner planet PDFs. Considering such a case, I propose an ansatz in which the stochastic dynamics of the Poincaré variables ${\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}$ and ${\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}$ is given by the following superposition of a secular term and a chaotic term: ${r}\mathrm{=}{{r}}_{\sec }\mathrm{+}{{r}}_{\mathrm{chaos}}\mathit{,}$(22)where ${r}\mathrm{=}\mathrm{\left(}{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{,}{\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{\right)}$¹³. The stationary random variable r_sec describes the statistics of the LL dynamics and corresponds to the microcanonical ensemble (16), ${\mathit{\rho }}_{\sec }\mathrm{\left(}{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{,}{\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta }\left[{\mathit{r}}_{\mathit{k}}\mathrm{-}{\mathrm{(}\mathrm{2}{\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]}{\mathrm{2}\mathit{\pi }{\mathrm{(}\mathrm{2}{\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,}$(23)where ${\mathit{r}}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}{{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{2}}\mathrm{+}{{\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{2}}$. The time-dependent chaotic component is modelled by the following Langevin equation: ${{ṙ}}_{\mathrm{chaos}}\mathrm{=}\sqrt{\mathrm{2}{\mathit{D}}_{\mathit{k}}}\hspace{0.17em}{\eta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(24)where η is a white noise, i.e. $\begin{array}{ccc}& & \mathrm{⟨}{\eta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{0}\\ & & \mathrm{⟨}{\mathit{\eta }}_{\mathit{x}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{\eta }}_{\mathit{x}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)} \mathrm{⟨}{\mathit{\eta }}_{\mathit{y}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{\eta }}_{\mathit{y}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)}\\ & & \mathrm{⟨}{\mathit{\eta }}_{\mathit{x}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{\eta }}_{\mathit{y}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{0}\end{array}$(25)for all times t,t′ ≥ 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, ${{r}}_{\mathrm{chaos}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{0.}$(26)Since the random variables r_sec and r_chaos are independent, the PDF of r is given by the convolution of the PDF of r_sec, ρ_sec, with that of r_chaos, ρ_chaos. Therefore one has $\mathit{\rho }\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{R^{\mathrm{2}}}\mathrm{d}{z}\hspace{0.17em}{\mathit{\rho }}_{\sec }\mathrm{\left(}{z}\mathrm{\right)}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{chaos}}\mathrm{(}{r}\mathrm{-}{z}\mathit{,t}\mathrm{)}\mathit{.}$(27)The PDF of the random variable satisfying Eqs. (24)–(26)is the Green function of the Fokker-Planck equation $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{chaos}}\mathrm{\left(}{\xi }\mathit{,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}{\mathit{D}}_{\mathit{k}}\hspace{0.17em}{{\nabla }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{chaos}}\mathrm{\left(}{\xi }\mathit{,t}\mathrm{\right)}\mathit{.}$(28)This Green function is ${\mathit{\rho }}_{\mathrm{chaos}}\mathrm{\left(}{\xi }\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{\mathit{D}}_{\mathit{k}}\mathit{t}}\exp \left(\mathrm{-}\frac{{{\xi }}^{\mathrm{2}}}{\mathrm{4}{\mathit{D}}_{\mathit{k}}\mathit{t}}\right)\mathrm{·}$(29)Substituting Eqs. (23)and (29)in Eq. (27)one gets $\mathit{\rho }\mathrm{\left(}{\mathit{r}}_{\mathit{k}}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{r}}_{\mathit{k}}}{{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}}\exp \left(\mathrm{-}\frac{{\mathit{r}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{m}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}}\right){\mathit{I}}_{\mathrm{0}}\left(\frac{{\mathit{r}}_{\mathit{k}}\hspace{0.17em}{\mathit{m}}_{\mathit{k}}}{{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}}\right)\mathit{,}$(30)where I₀ indicates the modified Bessel function of the first kind and order zero, ${\mathit{m}}_{\mathit{k}}\mathrm{=}\mathrm{(}\mathrm{2}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}\mathrm{2}{\mathit{D}}_{\mathit{k}}\mathit{t}$. Therefore, one obtains a Rice distribution for the variable ${\mathit{r}}_{\mathit{k}}\mathrm{=}\mathrm{(}{{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{2}}\mathrm{+}{{\mathit{y}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. In the case where the variables x_k and y_k are dominated by one particular LL mode, i.e. | γ_kl | ≫ | γ_kn | for a certain index l and for all n ≠ l, the PDF of the eccentricity e_k turns out to be a Rice distribution with parameters ${\mathit{m}}_{\mathit{k}}\mathrm{=}\mathrm{\left|}{{\mathrm{O}}_{\mathrm{A}}}_{\mathit{kl}}\mathrm{\right|}\hspace{0.17em}\mathrm{(}\mathrm{2}\overline{{\mathit{P}}_{\mathit{l}}^{\mathrm{\prime }}}\mathit{/}{\mathrm{\Lambda }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}\mathrm{2}{\mathit{D}}_{\mathit{k}}\mathrm{(}{{\mathrm{O}}_{\mathrm{A}}}_{\mathit{kl}}{\mathrm{)}}^{\mathrm{2}}\mathit{t}\mathit{/}{\mathrm{\Lambda }}_{\mathit{k}}$. In this limiting case one thus recovers Laskar’s ansatz. This is because 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, m_k 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 m_k given 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_k clearly emerges in the present framework.

As shown in Table 2, the PDF (30)casts some light on the fits of Laskar 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 quasi-periodic 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 $\mathrm{⟨}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{r}}_{\mathit{k}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathrm{2}\mathrm{=}{\mathit{m}}_{\mathit{k}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{+}{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{+}\mathrm{2}{\mathit{D}}_{\mathit{k}}\mathit{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 either the conservation of AMD or 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 values. 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 r_sec and r_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 one 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 one 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, following Eq. (5). The problem would then be how to sample the corresponding microcanonical ensemble in an efficient way. 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 small to allow for further closeencounters. 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), the Laskar model is not unique because it requires an ad hoc prescription for the random evolution of the orbital elements between collisions. Statistical descriptions, such as those in Sects. 4 and 5 of this paper, can provide such a prescription on the basis of 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. Evidently, collisions have to be taken into account when addressing planet formation, while they are negligible for the future 5 Gyr of the solar system as long as one is concerned with the bulk of the planet PDFs. Finally, such a statistical theory would be even more pertinent to asteroids because of their large number. Moreover, the eccentricity distribution in the main asteroid belt (Malhotra & Wang 2017) seems to share similarities with the inner planet PDFs treated in this study. This should constitute the subject of future studies.

Acknowledgments

I am grateful to the two anonymous referees for their careful reading of the paper and the very pertinent points. I also acknowledge the fruitful discussions with A. Morbidelli, J. Touma, J. Laskar, J. Teyssandier, N. Cornuault, and L. Pittau. I finally thank J.-P. Beaulieu for his support and help in reviewing the manuscript.

References

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)² the Lagrangian of Sun and planets is given by (Landau & Lifshitz 1975) $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{L}\mathrm{=}{\mathit{L}}_{\mathrm{0}}\mathrm{+}\mathrm{3}\sum _{\mathit{a}}{\sum _{\mathit{b}}}^{\mathrm{\prime }}\frac{{\mathit{m}}_{\mathit{a}}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{G}{\mathit{m}}_{\mathit{b}}}{{\mathit{c}}^{\mathrm{2}}{\mathit{r}}_{\mathit{ab}}}\mathrm{+}\sum _{\mathit{a}}\frac{{\mathit{m}}_{\mathit{a}}{\mathit{v}}_{\mathit{a}}^{\mathrm{4}}}{\mathrm{8}{\mathit{c}}^{\mathrm{2}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{-}\sum _{\mathit{a}}{\sum _{\mathit{b}}}^{\mathrm{\prime }}\frac{\mathit{G}{\mathit{m}}_{\mathit{a}}{\mathit{m}}_{\mathit{b}}}{\mathrm{4}{\mathit{c}}^{\mathrm{2}}{\mathit{r}}_{\mathit{ab}}}\left[\mathrm{7}\hspace{0.17em}{{v}}_{\mathit{a}}\mathrm{·}{{v}}_{\mathit{b}}\mathrm{+}\mathrm{(}{{v}}_{\mathit{a}}\mathrm{·}{{n}}_{\mathit{ab}}\mathrm{)}\mathrm{(}{{v}}_{\mathit{b}}\mathrm{·}{{n}}_{\mathit{ab}}\mathrm{)}\right]\\ \begin{array}{c}\end{array} \mathrm{-}\sum _{\mathit{a}}{\sum _{\mathit{b}}}^{\mathrm{\prime }}{\sum _{\mathit{c}}}^{\mathrm{\prime }}\frac{{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathit{a}}{\mathit{m}}_{\mathit{b}}{\mathit{m}}_{\mathit{c}}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}{\mathit{r}}_{\mathit{ab}}{\mathit{r}}_{\mathit{ac}}}\mathit{,}\end{array}& & \end{array}$(A.1)where ${\mathit{L}}_{\mathrm{0}}\mathrm{=}{\sum }_{\mathit{a}}{\mathit{m}}_{\mathit{a}}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{+}{\sum }_{\mathit{a}}{\sum }_{\mathit{b}}^{\mathrm{\prime }}\mathit{G}{\mathit{m}}_{\mathit{a}}{\mathit{m}}_{\mathit{b}}\mathit{/}\mathrm{2}{\mathit{r}}_{\mathit{ab}}$ is the Newtonian Lagrangian. In the equation above a,b,c ∈ { 0,1,...,8 }, r_ab = | r_a−r_b |, n_ab = (r_a−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₀ ≪ 1 for i ∈ { 1,2,...,8 }. I then simplify Eq. (A.1)by keeping among the relativistic terms only those of order (v/c)², meaning that I neglect terms of order (m_i/m₀)(v/c)² and smaller. One obtains L = L₀ + δL, with $\mathit{\delta L}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathit{m}}_{\mathit{i}}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\left(\frac{{\mathit{v}}_{\mathit{i}}^{\mathrm{4}}}{\mathrm{4}}\mathrm{+}\mathrm{3}\frac{\mathit{G}{\mathit{m}}_{\mathrm{0}}{\mathit{v}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{r}}_{\mathit{i}}}\mathrm{-}\frac{{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}}\right)\mathit{,}$(A.2)where r_i is the Heliocentric position of planet i. The general relativistic correction to the planetary Hamiltonian is thus given by δℋ = −δL (Landau & Lifshitz 1969a). To obtain the secular correction one has to average δℋ over the Keplerian orbits of the non-interacting planets. Neglecting corrections of the order m_i/m₀, along a Keplerian orbit one has ${\mathit{v}}_{\mathit{i}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{G}{\mathit{m}}_{\mathrm{0}}}{{\mathit{a}}_{\mathit{i}}}\left(\frac{\mathrm{2}{\mathit{a}}_{\mathit{i}}}{{\mathit{r}}_{\mathit{i}}}\mathrm{-}\mathrm{1}\right)\mathit{,}$(A.3)where a_i is the semimajor axis of planet i. Moreover, averaging over an orbital period one obtains $\mathrm{⟨}{\mathit{r}}_{\mathit{i}}^{-1}\mathrm{⟩}\mathrm{=}{\mathit{a}}_{\mathit{i}}^{-1}$ and $\mathrm{⟨}{\mathit{r}}_{\mathit{i}}^{-2}\mathrm{⟩}\mathrm{=}{\mathit{a}}_{\mathit{i}}^{-2}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$, with e_i the eccentricity of planet i. Substituting Eq. (A.3)in δℋ and taking the average, one finally finds $\mathrm{⟨}\mathit{\delta }\mathrm{\mathscr{H}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}}{{\mathit{c}}^{\mathrm{2}}{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}\left[\frac{\mathrm{15}}{\mathrm{8}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right]\mathrm{·}$(A.4)Discarding the first term as it is constant in the secular dynamics, one obtains that the dominant secular contribution of general relativity is ${\mathrm{\mathscr{H}}}_{\sec \mathit{,}\hspace{0.17em}\mathrm{GR}}\mathrm{=}\mathrm{-}{\sum }_{\mathit{i}}\mathrm{3}{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}\mathit{/}{\mathit{c}}^{\mathrm{2}}{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$.

Appendix B: Laplace-Lagrange Hamiltonian

The matrix A appearing in the quadratic Hamiltonian (6)can be expressed as ${\mathrm{A}}_{\mathit{kl}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{k}\ne \mathit{l}\\ & \mathrm{if}\mathit{k}\mathrm{=}\mathit{l,}\end{array}$(B.1)where I have defined $\begin{array}{ccc}{\mathit{I}}_{\mathit{nkl}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}\mathrm{d}\mathit{\theta }\frac{\cos \mathrm{\left(}\mathrm{2}\mathit{n\theta }\mathrm{\right)}}{{\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{a}}_{\mathit{k}}{\mathit{a}}_{\mathit{l}}{\sin }^{\mathrm{2}}\mathit{\theta }\mathit{/}\mathrm{(}{\mathit{a}}_{\mathit{k}}\mathrm{-}{\mathit{a}}_{\mathit{l}}{\mathrm{)}}^{\mathrm{2}}\mathrm{]}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathrm{·}& & \end{array}$Similarly, the matrix B is given by ${\mathrm{B}}_{\mathit{kl}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{k}\ne \mathit{l}\\ & \mathrm{if}\mathit{k}\mathrm{=}\mathit{l}\mathit{.}\end{array}$(B.2)From formula (B.2)it is easy to check that matrix B verifies the following relation: $\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{B}}_{\mathit{kl}}\sqrt{{\mathrm{\Lambda }}_{\mathit{l}}}\mathrm{=}\mathrm{0.}$(B.3)This implies that one of the s_k eigenvalues is zero and ${}^{\mathrm{(}}\sqrt{{\mathrm{\Lambda }}_{\mathrm{1}}}\mathit{,}\mathit{...}\mathit{,}{\sqrt{{\mathrm{\Lambda }}_{\mathit{N}}}}^{\mathrm{)}}{\hspace{0.17em}}^{\mathrm{􏼮}}\sqrt{{\sum }_{\mathit{k}}{\mathrm{\Lambda }}_{\mathit{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 Delaunay 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 ω is then written$\frac{\mathit{\omega }\mathrm{\left(}{\mathit{C}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}\mathit{N}}}\mathrm{=}\left[\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathrm{\Lambda }}_{\mathit{k}}}\mathrm{d}{\mathit{P}}_{\mathit{k}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{k}}\mathrm{-}{\mathit{P}}_{\mathit{k}}\mathrm{)}}\mathrm{d}{\mathit{Q}}_{\mathit{k}}\right]\mathit{\delta }\left({\mathit{C}}_{\mathrm{0}}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{Q}}_{\mathit{k}}\right)\mathit{,}$(C.1)where N = 4 in the present analysis, C₀ is the initial total AMD of the inner planets, and the factor (2π)^2N comes from the integration over the angle variables p_k and q_k. Since C₀< Λ_k for all k (Laskar 1997, Table 1), one can restrict the above integration to the hypercube [ 0,C₀] ^2N. To calculate the density of states one can therefore write ω(C₀) = dΣ(C₀)/dC₀, with $\frac{\mathrm{\Sigma }\mathrm{\left(}{\mathit{C}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{\left(}\mathrm{2}\mathit{N}\mathrm{\right)}}}\mathrm{=}\left[\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{C}}_{\mathrm{0}}}\mathrm{d}{\mathit{P}}_{\mathit{k}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{C}}_{\mathrm{0}}}\mathrm{d}{\mathit{Q}}_{\mathit{k}}\right]\hspace{0.17em}\hspace{0.17em}\mathit{H}\left({\mathit{C}}_{\mathrm{0}}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{Q}}_{\mathit{k}}\right)\mathit{,}$(C.2)where H is the Heaviside step function. Performing the integrations by iteration, one obtains Σ(C₀) = (2π)^2NC₀^2N/ (2N) ! and then $\mathit{\omega }\mathrm{\left(}{\mathit{C}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}\mathit{N}}\frac{{{\mathit{C}}_{\mathrm{0}}}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}}{\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}\mathrm{·}$(C.3)The joint microcanonical PDF of P_k and Q_k (i.e. marginalized over the remaining modified Delaunay variables) is then given by $\mathit{\rho }\mathrm{\left(}{\mathit{P}}_{\mathit{k}}\mathit{,}{\mathit{Q}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}\mathit{N}}}{\mathit{\omega }\mathrm{\left(}{\mathit{C}}_{\mathrm{0}}\mathrm{\right)}}\frac{\mathrm{d}}{\mathrm{d}{\mathit{C}}_{\mathrm{0}}}\left[\frac{\mathrm{(}{\mathit{C}}_{\mathrm{0}}\mathrm{-}{\mathit{P}}_{\mathit{k}}\mathrm{-}{\mathit{Q}}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}\right]$(C.4)for P_k + Q_k ≤ C₀, otherwise it is zero. Thus one obtains $\mathit{\rho }\mathrm{\left(}{\mathit{P}}_{\mathit{k}}\mathit{,}{\mathit{Q}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}\mathrm{)}}{{{\mathit{C}}_{\mathrm{0}}}^{\mathrm{2}}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{Q}}_{\mathit{k}}}{{\mathit{C}}_{\mathrm{0}}}\right)}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{3}}$(C.5)for P_k + Q_k ≤ C₀. Switching to the variables e_k and i_k, one finds $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathit{,}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}\mathrm{)}}{{\mathit{\eta }}_{\mathit{k}}^{\mathrm{2}}}{\mathit{e}}_{\mathit{k}}\sin {\mathit{i}}_{\mathit{k}}\hspace{0.17em}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }{\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}}{{\mathit{\eta }}_{\mathit{k}}}\right)}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{3}}\end{array}& & \end{array}$(C.6)for ${\mathit{\eta }}_{\mathit{k}}^{-1}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}\right)\mathrm{\le }\mathrm{1}$, with η_k = C₀/ Λ_k. The PDF of e_k is then given by $\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{i}}^{\mathit{\star }}}\mathit{d}{\mathit{i}}_{\mathit{k}}\hspace{0.17em}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathit{,}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}\mathit{,} {\mathit{e}}_{\mathit{k}}\mathrm{\le }\sqrt{\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}$(C.7)where $\cos {\mathit{i}}^{\mathit{\star }}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. Performing the integral, one obtains $\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}{{\mathit{\eta }}_{\mathit{k}}}\frac{{\mathit{e}}_{\mathit{k}}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}}{\left[\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}^{-1}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\right)\right]}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}}$(C.8)for ${\mathit{e}}_{\mathit{k}}\mathrm{\le }\sqrt{\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}$. Similarly, to obtain the PDF of i_k one has to calculate $\mathit{\rho }\mathrm{\left(}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{e}}^{\mathit{\star }}}\mathrm{d}{\mathit{e}}_{\mathit{k}}\hspace{0.17em}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathit{,}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}\mathit{,} \cos {\mathit{i}}_{\mathit{k}}\mathrm{\ge }\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}\mathit{,}$(C.9)where e^⋆ = [ 1−(1−η_k)²/ cos²i_k] ^1/2. Performing the integration, one finds $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{\rho }\mathrm{\left(}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\tan {\mathit{i}}_{\mathit{k}}\hspace{0.17em}{\mathrm{[}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}^{-1}\mathrm{(}\mathrm{1}\mathrm{-}\cos {\mathit{i}}_{\mathit{k}}{\mathrm{)}}^{\mathrm{]}}}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}}\hspace{0.17em}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }\left[\frac{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}{{\mathit{\eta }}_{\mathit{k}}}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}^{-1}\mathrm{(}\mathrm{1}\mathrm{-}\cos {\mathit{i}}_{\mathit{k}}\mathrm{)}}{\cos {\mathit{i}}_{\mathit{k}}}\right]\end{array}& & \end{array}$(C.10)for cosi_k ≥ 1−η_k. As in the present study one has η_k ≪ 1 (Laskar 1997, Table 1), one can simplify these PDFs obtaining $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}{\mathit{e}}_{\mathit{k}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}{{\mathit{\eta }}_{\mathit{k}}}{\mathit{e}}_{\mathit{k}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\eta }}_{\mathit{k}}}\right)}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}}\mathit{,} {\mathit{e}}_{\mathit{k}}\mathrm{\le }\sqrt{\mathrm{2}{\mathit{\eta }}_{\mathit{k}}}\\ \mathit{\rho }\mathrm{\left(}{\mathit{i}}_{\mathit{k}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{1}}{{\mathit{\eta }}_{\mathit{k}}}{\mathit{i}}_{\mathit{k}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{i}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\eta }}_{\mathit{k}}}\right)}^{\mathrm{2}\mathit{N}\mathrm{-}\mathrm{2}}\mathit{,} {\mathit{i}}_{\mathit{k}}\mathrm{\le }\sqrt{\mathrm{2}{\mathit{\eta }}_{\mathit{k}}\mathit{.}}\end{array}$In this limit one has that $\begin{array}{ccc}\mathrm{⟨}{\mathit{e}}_{\mathit{k}}\mathrm{⟩}& \mathrm{=}& \\ \mathrm{⟨}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}\mathrm{⟩}& \mathrm{=}& \end{array}$where Γ is the Gamma function, and the average AMD of planet k is thus given by $\mathrm{⟨}{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{Q}}_{\mathit{k}}\mathrm{⟩}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{k}}\frac{\mathrm{⟨}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}_{\mathit{k}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{2}}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{0}}}{\mathit{N}}\mathrm{·}$(C.15)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 (C.11)and (C.12)present the same overall structure as those predicted by the ansatz (18). In particular, as Eqs. (C.13)and (C.14)show, the mean and standard deviation of eccentricity and inclination are proportional to ${\mathit{\eta }}_{\mathit{k}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\propto }{\mathit{m}}_{\mathit{k}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$.

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 $\mathrm{(}{\mathit{L}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{L}}_{\mathit{y}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{/}{\mathit{L}}_{\mathit{z}}$.

Appendix D: Ansatz of equiprobability: sampling

Conservation of AMD and energy in Eqs. (18)and (19)is equivalent to the two following conditions: $\begin{array}{ccc}{\mathit{C}}_{\mathrm{0}}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{4}}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}\\ {\mathit{H}}_{\mathrm{0}}& \mathrm{=}& \mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{4}}\mathrm{(}{\mathit{g}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}{\mathit{s}}_{\mathit{k}}{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{,}\end{array}$with ${\mathit{C}}_{\mathrm{0}}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{4}}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{+}\overline{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}$ and ${\mathit{H}}_{\mathrm{0}}\mathrm{=}\mathrm{-}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{4}}{}^{\mathrm{(}}\mathit{g}_{\mathit{k}}\overline{{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}}\mathrm{+}{\mathit{s}}_{\mathit{k}}{\overline{{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}}}^{\mathrm{)}}$, the bar indicating the initial value of the corresponding variable (see Sect. 4 and footnote 5). These equations can be rewritten as $\begin{array}{ccc}\mathrm{1}& \mathrm{=}& \\ \mathit{h}& \mathrm{=}& \end{array}$where ${\mathit{a}}_{\mathit{k}}\mathrm{=}{\mathit{P}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{/}{\mathit{C}}_{\mathrm{0}}$, ${\mathit{b}}_{\mathit{k}}\mathrm{=}{\mathit{Q}}_{\mathit{k}}^{\mathrm{\prime }}\mathit{/}{\mathit{C}}_{\mathrm{0}}$ and h = −H₀/C₀. Moreover, I recall that a_k,b_k ≥ 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. (D.3). To evaluate the PDF (19)a_k and b_k have to verify both Eqs. (D.3)and (D.4).

Uniform sampling of N positive real numbers that add up to unity, as in Eq. (D.3), 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₁,...,u_N−1 }, with u_i−1 ≤ u_i for i ∈ { 2,...,N−1 }. Then, the numbers { u₁,u₂−u₁,...,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. (D.3)and (D.4)there is probably no direct sampling method. However, I present here an efficient algorithm with rejection. Let us assume that h ≥ 0, as in the present study. A similar algorithm can be constructed in case h is negative. One can multiply Eq. (D.3)by g_max = max_k { a_k }, subtract Eq. (D.4)and finally divide by g_max−h (one can assume this to be different from zero, otherwise the sampling is trivial). One obtains $\mathrm{1}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{4}}\frac{{\mathit{g}}_{\max }\mathrm{-}{\mathit{g}}_{\mathit{k}}}{{\mathit{g}}_{\max }\mathrm{-}\mathit{h}}{\mathit{a}}_{\mathit{k}}\mathrm{+}\frac{{\mathit{g}}_{\max }\mathrm{-}{\mathit{s}}_{\mathit{k}}}{{\mathit{g}}_{\max }\mathrm{-}\mathit{h}}{\mathit{b}}_{\mathit{k}}\mathit{.}$(D.5)The coefficients of a_k and b_k in Eq. (D.5)are positive and the term a_k^⋆, such that g_k^⋆ = g_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^⋆ = 1−S to obtain a set of numbers following the target distribution. Otherwise, the numbers are rejected and Eq. (D.5)sampled again.

All Tables

All Figures

	Fig. 1 Eccentricity and inclination PDFs of the solar system giant planets. The PDFs of L08 are plotted in solid lines, while those predicted by the microcanonical distribution function (16)are shown in dashed lines. The inclinations are computed with respect to the J2000 ecliptic.
In the text

Fig. 2 Eccentricity and inclination PDFs of the solar system inner planets. The PDFs of L08 are plotted in solid lines at three different times, t = 500 Myr, 2.5 Gyr, and 5 Gyr, to illustrate chaotic diffusion (the PDF broadens and its peak decreases as time increases; the inclination PDFs are not the original PDFs, but the Rice distributions fitted by Laskar, see Sect. 3 for explanations). The PDFs predicted by the microcanonical ensemble of the Laplace-Lagrange theory, Eq. (16), are shown in dashed line. The inclinations are computed with respect to the J2000 ecliptic.

In the text

	Fig. 3 Eccentricity and inclination PDFs of the solar system inner planets. The same as Fig. 2, except that the dashed lines represent here the predictions of the ansatz (18).
In the text

	Fig. 4 Eccentricity PDF of the Earth. The same as Fig. 3 with the addition of the prediction of Eq. (19)in dash-dotted line.
In the text

	Fig. 5 Inclination PDF of the Earth. The same as Fig. 3 with the addition of the prediction of Eq. (19)in dash-dotted line.
In the text

	Fig. 6 Eccentricity PDF of Mercury. The same as Fig. 3 with the addition of the predictions of Eq. (19)in dash-dotted line and Eq. (20)in dotted line.
In the text