© ESO, 2017

1. Introduction

The numerical integration of the secular equations of motion for the eight planets of the solar system including the moon (Laskar 1989) has shown that the solar system is chaotic and its Lyapunov time τ_i has been estimated around 10 Myr. The chaos is significant enough such that a single integration of the equations of motion is not at all representative of the state of the solar system after 10 Myr. One property of chaotic motion is that it increases exponentially any difference in the initial positions of the planets. Therefore, an error of few meters in the initial positions of the planets leads, on a tens of Myr timescale, to a complete indetermination on the actual position of the planets. This has been confirmed by more and more accurate numerical integrations of the solar system dynamics without averaging and including perturbations induced by the largest asteroids (Laskar et al. 2011a; Fienga et al. 2009).

It is a natural question to address the effect of external sources of chaoticity on the dynamics of the planets. Those external sources could be for example, the comets, the asteroids, or even the radiative pressure from the sun. From all those external sources, a computation of orders of magnitude shows that the asteroids should be the main ones. The asteroid belt is a very good candidate for white noise on the planetary dynamics, because asteroids, especially the larger ones Ceres, Pallas, and Vesta, are large enough to exert a non negligible gravitational interaction on the planets. Moreover numerical studies by Laskar et al. (2011b) have shown that their dynamics is chaotic with a Lyapunov time τ_ϕ ≃ 10⁴ yr much smaller than the Lyapunov time of 10 Myr of the solar system. Numerical simulations of the dynamics of the solar system including the asteroids, compared to others without asteroids (Laskar et al. 2011b), show that they can indeed change the secular behavior of planetary orbits. It has, however been observed that simulations with asteroids do not affect significantly the Lyapunov time of the solar system (Laskar et al. 2011b).

Nevertheless, a systematic investigation of the influence of asteroids on planets has not been done yet. We have numerical evidence that asteroids perturbs planetary motion on tens of Myr timescale, but we don’t know which physical parameters do control this timescale. The aim of this paper is to answer this question, and to give a theoretical support to numerical simulations. We find that the dominant effect of asteroids is on planetary mean longitudes. The first result is that, if one would neglect secular evolutions of the orbital parameters, the error on the longitude of Mars would be described by a superdiffusion, with standard deviation growing in the following way: ${{\mathrm{(}}^{}\mathit{t}\mathit{/}{\mathit{\tau }}_{\mathrm{diff}}{\mathrm{)}}^{}}^{\mathrm{3}\mathit{/}\mathrm{2}}$, with ${\mathit{\tau }}_{\mathrm{diff}}\mathrm{\propto }{\left({\left(\frac{{\mathit{m}}_{\mathrm{a}}}{{\mathit{M}}_{\mathrm{S}}}\right)}^{\mathrm{2}}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}}{{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(1)where is Newton’s constant of gravitation, m_a,M_S are the masses of the asteroid and the sun respectively, a_p is the semi-major axis of Mars, and τ_ϕ is the Lyapunov time of the asteroid. Altogether, the magnitude of τ_diff is of order of 20 Myr. The second result is that the asteroids influence on the secular dynamics leads to a superdiffusion of the orbital parameters eccentricities and inclinations, with standard deviation growing like ${{\mathrm{(}}^{}\mathit{t}\mathit{/}{\mathit{\tau }}_{\sec }{\mathrm{)}}^{}}^{\mathrm{3}\mathit{/}\mathrm{2}}$, with ${\mathit{\tau }}_{\sec }\mathrm{\propto }{\left({\left(\frac{{\mathit{m}}_{\mathrm{a}}}{{\mathit{M}}_{\mathrm{S}}}\right)}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\sec }^{\mathrm{2}}}{{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(2)where ν_sec is the typical order of magnitude of the secular frequencies of the Laplace–Lagrange system, about a few arcsec/yr. We can then estimate that τ_sec is of order of 10 Gyr.

At a conceptual level, it may be useful to discuss the meaning of stochasticity for deterministic dynamical systems. First, for times much longer than the inverse of the largest Lyapunov exponent of a dynamical system (what we call the Lyapunov time), chaotic dynamics is hardly distinguishable from stochastic dynamics. One can consider that two deterministic trajectories in the dynamical system differing in their initial conditions are like two independent realizations of a stochastic process. Precise mathematical results can be obtained in several frameworks. For instance let us consider the case of two coupled deterministic systems evolving on distinct fast and slow timescales, assuming that the fast deterministic system is chaotic. One can prove mathematically that in the limit of a large timescale separation, the effect of the fast deterministic system is equivalent to the effect of a white noise which properties are related to the statistics of the fast chaotic dynamics. Unveiling the required hypothesis for such results to be valid is a difficult and fascinating part of ergodic theory (see for instance (Kifer 2004, 2009) and reference therein). The key message, which has been commonly accepted in statistical mechanics for several decades, is that under generic conditions, there is no fundamental differences between deterministic complex chaotic dynamics and stochastic ones, beyond the very difficult and subtle mathematical problems (see for instance Gallavoti 2008; Gallavotti & Cohen 1995; Ruelle 1980). As an illustration of this idea, Pavliotis & Stuart (2008) have shown that a pendulum coupled to the chaotic dynamics of the Lorentz system is equivalent to a stochastic process. Another classical example is the relation between Brownian motion and the deterministic dynamics of atoms: a large particle coupled to a large number of smaller particles has a stochastic dynamics which is very accurately described by the Langevin equation for Brownian motion. In this latter example, the small particles act as a white noise force on the large particle.

The case of the coupling between fast stochastic dynamics and slow dynamics is less technical from a mathematical point of view than the deterministic one, while equivalent at a formal level. The related tools, known as stochastic averaging are described in classical textbooks (Pavliotis & Stuart 2008; Gardiner 1985; Freidlin & Wentzell 1984). When chaotic dynamics is what we call “mixing”, which means that the system looses the memory of its initial condition fast enough, it is statistically equivalent to a white noise when it is considered on a timescale much longer than the mixing time. The properties of the white noise are given by a Green–Kubo formula, as we explain in Sect. 4. This result will be the main theoretical starting point of our analysis.

The importance of stochastic approaches for the solar system dynamics has been understood quite recently. For instance Laskar (2008) numerically integrated 1000 trajectories of planets of the solar system with small differences in the initial conditions in order to compute numerically the evolution of the probability distribution functions of the eccentricities and inclinations of the main planets. More recent work include (Tremaine 2015; Mogavero 2017). Mogavero tried to reproduce the probability distributions of planets numerically computed by Laskar (2008) using simply the hypothesis of equiprobability in phase space and taking into account the integrals of motions. The work of Batygin et al. (2015) about the chaotic motion of Mercury is, to our knowledge, the first one in the context of celestial mechanics that models a chaotic Hamiltonian dynamics with a stochastic process on a long timescale, where the order of magnitude of the stochastic process properties are discussed precisely. The case of the stochastic effect of the asteroids on the planets has motivated the theoretical work of Lam & Kurchan (2014). They assumed that the effect of asteroids on planetary Keplerian motions is equivalent to that of a white noise acting on an integrable system. In the present paper, we show that the white noise assumption is an exact asymptotic result when the timescale considered is much larger than the Lyapunov time τ_ϕ of asteroids, and we give the theoretical expression and the order of magnitude of the white noise process. If the white noise has an amplitude δ, Thu Lam and Kurchan have shown that the Lyapunov time of the perturbed system was scaling like ${{}^{\mathrm{\right(}}\frac{\mathrm{1}}{\mathit{\delta }}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{3}}$. The first technical part of our paper is very similar to that of Lam & Kurchan (2014). Rather than studying the Lyapunov exponent, we study the diffusion of planetary mean longitudes, which is more relevant if one is interested in the evolution of the probability distribution functions. The relevant timescales are however similar: the Lyapunov time (the inverse of the Lyapunov exponent) and the diffusion time τ_diff are both of the order of 20 Myr (see Table 2). The second technical part of our paper connects correctly for the first time the deterministic dynamics of the asteroids with the model with white noises, using stochastic averaging. This allows us to compute correctly for the first time the order of magnitude of δ, and to actually compute the relevant timescales for the solar system. We have applied these results to both the diffusion of the planetary mean longitudes and the orbital elements for the Laplace–Lagrange dynamics.

In Sect. 2 we give the precise mathematical formulation of the question we are interested in. Section 2 may thus appear a little abstract from the point of view of celestial mechanics, but it is essential to understanding the scope of the present work. As we explain in this introduction, the scope is larger than simply answering the question of long-term influence of asteroids and could be applied to many dynamical systems. In Sect. 3 we introduce a simplified Hamiltonian model to describe the dynamics of Mars under the influence of asteroids, and we justify why this model is relevant to computing orders of magnitude. The stochastic method is explained in Sect. 4. This section is a little technical, though most of the difficulties are detailed in Appendix A. Section 5 contains the main results of the paper, the computation of the timescale τ_diff. We then investigate the influence of asteroids on planetary secular motion in Sect. 6 and we conclude our discussion in Sect. 7. In order to help the reader, we have gathered the main notations appearing in our article in Appendix C.

2. The theoretical framework to compare intrinsic and extrinsic chaos

In this section, we formulate the problem in a more precise mathematical set-up. The present discussion goes far beyond the particular problem of the influence of asteroids on planetary dynamics. The influence of asteroids is one simple application of our results, but the mathematical framework should find many other applications in celestial mechanics.

We have considered a dynamical system of the form $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{\eta }{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{ϵg}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{.}$(3)The system has two small parameters η ≪ 1 and ϵ ≪ 1, x is a multidimensional vector. The field f₀ + ηf₁ depends only on x, it can be considered as the intrinsic dynamics of x. The field g depends both on x and the time, it acts as an external perturbation on the system.

We assumed that the zeroth order dynamics defined by $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(4)is an integrable Hamiltonian dynamics. When the small perturbation ηf₁ is added, we assumed that the intrinsic dynamics $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{\eta }{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}$(5)is chaotic with a Lyapunov time τ_i (the Lyapunov time is defined as the inverse of the largest Lyapunov exponent). In the present paper, we consider two cases for the dynamics (5). Section 3 will present the case in which the integrable dynamics f₀ is the Keplerian dynamics of planets, ηf₁ being the perturbative function coming from the gravitational interactions between the planets. In Sect. 6, f₀ corresponds to the secular system of Laplace–Lagrange and ηf₁ gathers all non-linear terms of order two and more in eccentricities and inclinations. From the results of the work of Laskar (1989), we know that the full non-linear secular dynamics is chaotic with an intrinsic Lyapunov time τ_i of the order of 10 Myr.

External sources may add a small perturbation on the intrinsic dynamics (5) described by the term ϵg(x,t). In the present paper for example, we have considered the asteroids as an external source of noise for planets of the solar system. The intrinsic dynamics is already chaotic, and therefore the complete system (3) is also chaotic. On the other hand, the integrable system perturbed by the external source of noise $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{ϵg}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}$(6)is chaotic with a Lyapunov time τ_e. In the case (6), the chaos is due to the external perturbation.

As we explained in the introduction, if the perturbation g(x,t) satisfies a mixing condition, which means that its time correlation function for any fixed value of x decreases fast enough, the deterministic system of Eqs. (6) is equivalent to a Langevin process $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{\xi }\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}$(7)where ξ(x,t) is a white noise, which amplitude can be expressed with the properties of g(x,t) with a Green-Kubo formula. This point will be discussed in Sect. 4. Once the system (6) has been written as a stochastic process (7), it is quite easy to give the order of magnitude of τ_e. Equations (5) and (6) define two different regimes depending on the Lyapunov times τ_i and τ_e:

• 1.

  The regime τ_i ≪ τ_e defines a regime of intrinsic chaos. On a time of order of the internal Lyapunov time τ_i, the effect of the external perturbation ϵg is small. Thus, the probability distributions of the variable x are essentially the same, to leading order in ϵ, for the full system (3) and for the intrinsic dynamics (5).

• 2.

  If on the contrary τ_e ≪ τ_i, then the external perturbation creates chaos in the integrable system, in the sense that the system looses the memory of its initial condition before the intrinsic chaos can develop. For intermediate times between τ_e and τ_i, the complete dynamics (3) can thus be described by a stochastic process, and the probability distributions may be strongly influenced by the external perturbation.

Using the present framework, our question formulates in a very simple way. Let the complete set of Eqs. (3) be the equations for planetary motion of the solar system perturbed by the asteroid belt, are we in the regime of dominant intrinsic chaos with τ_i ≪ τ_e or in the regime of external source of chaos with τ_e ≪ τ_i? What is then the order of magnitude of τ_e describing the interactions between asteroids and planets?

3. A simplified model for planet-asteroid interactions

For times smaller than the Lyapunov time of the solar system, the secular motion of planetary orbits is very accurately described by the periodic solution of the Laplace–Lagrange equations. Except for the smallest planet Mercury, the planetary orbital elements eccentricities and inclinations remain very small (less than 0.15 for the eccentricity, and less than ten degrees for the inclination) in the Myr timescale (Laskar 2008). The computations to be performed in the following could be done without fundamental difficulties for elliptic and inclined trajectories. However solving these equations for the elliptic motion is technically much more tedious than for restricted planar and circular motions, and it would not change the orders of magnitude to leading order in eccentricities and inclinations. To study the order of magnitude of the perturbation induced by the asteroids on the planets, we have thus introduced a simplified model where the orbits of celestial bodies are circular and coplanar. To describe the motion of the planet, we thus kept only the two orbital elements mean longitude λ and semi-major axis a.

The main perturbation will be on planet Mars whose orbit is the closest to the asteroid belt, but our method can be straightforwardly applied to other planets, using the correct orbital elements. Most of the mass of the asteroids comes only from the contribution of the largest ones Ceres, Vesta, and Pallas (about 58% of the total mass of the asteroids). The physical properties of those three asteroids is summarized in Table 1. In our simplified model, we will thus retain only the planet Mars and one asteroid. Asteroids have chaotic motions because of gravitational perturbation by the planets, but also because of interactions between each other as shown by Laskar et al. (2011b). They are thus not independent. Yet, their motions are decorrelated after the Lyapunov time, and the perturbations induced by each of the asteroids can thus be considered as independent for times larger than the Lyapunov time. We will add the individual contributions of each asteroid in the final result to obtain the right order of magnitude. The simplified model of the planet Mars perturbed by one asteroid is described by the Hamiltonian $\mathit{H}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}{\mathit{m}}_{\mathrm{p}}}{\mathrm{2}{\mathit{a}}_{\mathrm{p}}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}\mathit{m}}{\mathrm{2}\mathit{a}}\mathrm{-}\frac{\mathrm{𝒢}\mathit{m}{\mathit{m}}_{\mathrm{p}}}{\left|{\mathit{a}}_{\mathrm{p}}{\mathrm{e}}^{\mathrm{i}\mathit{\lambda }}\mathrm{-}\mathit{a}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }}\right|}\mathit{,}$(8)where M_S,m_p,m stand for the masses of the Sun, Mars, and the asteroid respectively, λ,ϕ are the mean longitudes of Mars and the asteroid, and a_p and a their respective semi-major axis. The Sun is considered as fixed. This means that we only retain the direct term of the perturbative function and we take the real mass m instead of the reduced mass $\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{=}\frac{\mathrm{1}}{\mathit{m}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{S}}}$ as should be done rigorously.

In physical units, it is difficult to see where is the small parameter in the problem. Therefore, we rescaled all physical variables. The mass of the sun was taken as the unit mass, M_S = 1, and the reduced mass $\mathit{ϵ}\mathrm{:}\mathrm{=}\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{S}}}$ of the asteroid is thus very small (Table 1). Then we changed the units for time and length, the astronomic unit became the new length scale, and we chose the units of time such that the actual Keplerian period of Mars was 2π. From the relation $\frac{{\mathit{T}}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}}\mathrm{=}\frac{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}}$ we obtained the new time unit $\frac{{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}}{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}}\mathrm{=}\mathrm{1.}$ Denoting ${\mathit{n}}_{\mathrm{p}}\mathrm{:}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{T}}_{\mathrm{p}}}$ the Keplerian pulsation of Mars, this means that n_p = 1. Finally, we let $\mathrm{\Lambda :}\mathrm{=}{\mathit{ϵ}}_{\mathrm{p}}\sqrt{{\mathit{a}}_{\mathrm{p}}}$ be the canonical momentum associated to λ, we perform the canonical change of variables $\mathit{H}\mathrm{\leftarrow }\frac{\mathit{H}}{{\mathit{ϵ}}_{\mathrm{p}}}\mathit{,}\mathrm{\Lambda }\mathrm{\leftarrow }\frac{\mathrm{\Lambda }}{{\mathit{ϵ}}_{\mathrm{p}}}$.

In the work of Laskar et al. (2011b), the orbits of asteroids have been shown to be chaotic and the Lyapunov time τ_ϕ on their longitudes has been computed numerically. The Lyapunov times for the three main asteroids are given in Table 1. In the present model in which the planet Mars is perturbed by one asteroid, it is the chaos of the asteroid’s motion that breaks the periodic Keplerian regular motion of Mars. Of course the chaotic motion of the asteroid does not come from the influence of Mars alone, but rather from interaction and close encounters with other asteroids and on the influence of the giant planets. The retroactive influence of Mars on the asteroid can be considered to be negligible compared to the influence of Jupiter or Saturn, and the interaction with Mars cannot substantially change the characteristics of the orbit of the asteroid, in particular its Lyapunov time. This is why we have not solved the equations of motion for the asteroid. To capture the physical phenomena coming from the gravitational interaction between Mars and the asteroid, the trajectory of the asteroid has to be considered as an input function in the model which properties come from more precise numerical studies. We thus took the semi-major axis a as a constant and the phase of the asteroid ϕ(t) as a function of time with correlation time of the order of the Lyapunov time τ_ϕ. The trajectory of Mars in our model is thus a functional of the trajectory ϕ(t) of the asteroid. The reader should always bear in mind that those strong hypothesis are used with the aim of giving orders of magnitude and not precise quantitative results.

The Hamiltonian of the simplified model we will study is $\mathit{H}\mathrm{\left(}\mathrm{\Lambda }\mathit{,\lambda ,t}\mathrm{\right)}\mathrm{:}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathrm{\Lambda }}^{\mathrm{2}}}\mathrm{-}\frac{\mathit{ϵ}}{\left|{\mathrm{\Lambda }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{i}\mathit{\lambda }}\mathrm{-}\mathit{a}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right|}\mathrm{·}$(9)The Hamiltonian (9) does not conserve the total energy and angular momentum of the asteroid and Mars. But to first order in ϵ, energy and angular momentum are the ones given by the Keplerian orbit and depend thus only on the canonical momentum Λ. Our model is consistent if the change in energy and angular momentum occurs on a timescale much larger than the time we are considering for the perturbation of Mars. This point will be checked a posteriori in Sect. 5 where we obtain the time τ_diff over which the perturbation on Mars becomes large.

In the final paragraph of this section, we discuss more precisely in which sense the position of Mars should be considered as a stochastic variable associated to a probability distribution. What is usually done in numerical simulations of chaotic planetary motion is to choose a large number of initial conditions differing only by a small shift in the initial position. Because of chaos, the different trajectories do not stay close together but separate exponentially fast on a timescale given by the Lyapunov exponent. After a sufficiently long time compared to the Lyapunov time, the positions of all trajectories give a distribution. This distribution is an estimation of all possible positions that could be reached from a uniform distribution on a very small set of initial conditions. In this sense, it is a probability distribution. In the simplified model (9), we have not fixed the initial position of the asteroid. We will thus study the motion of Mars for different possible realizations of the function ϕ(t) and we take for ϕ(0) a uniform distribution over the range [0,2π]. We made this choice because the incertitude on the longitude of the asteroid becomes complete after a time large enough compared to the Lyapunov time τ_ϕ. The position of Mars has a probability distribution because it is conditioned by the realizations of the stochastic function ϕ(t). We emphasize that the trajectories of Mars will not separate exponentially with time, as would have been the case if we had just considered a set of very close initial conditions for the position of Mars and one single function ϕ(t). In our model, ϕ(t) is a random function with a probability distribution P[ϕ]. Given this probability distribution, we want to obtain the probability distribution of the canonical variables Λ,λ. Instead of an exponential separation, we will obtain a diffusive behavior as shown in the next section.

4. The stochastic longitude evolution as a diffusion

From the Hamiltonian (9), we obtain the set of Hamilton equations for λ,Λ as $\begin{array}{ccc}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{t}}& \mathrm{=}& \\ \frac{\mathrm{d\Lambda }}{\mathrm{d}\mathit{t}}& \mathrm{=}& \mathrm{-}\mathit{ϵ}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{(}\mathrm{\Lambda }\mathit{,\lambda }\mathrm{-}\mathit{\varphi }\mathrm{(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$(10)where we have introduced the Keplerian pulsation ${\mathit{n}}_{\mathrm{p}}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{\Lambda }}^{\mathrm{3}}}$ and the gravitational interaction $\mathit{G}\mathrm{(}\mathrm{\Lambda }\mathit{,\lambda }\mathrm{-}\mathit{\varphi }\mathrm{)}\mathrm{:}\mathrm{=}\frac{-1}{\left|{\mathrm{\Lambda }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathit{\varphi }\mathrm{)}}\mathrm{-}\mathit{a}\right|}$. For technical reasons, it is convenient to consider the variable p: = n_p(Λ) instead of Λ and write the Hamilton Eqs. (10) as $\begin{array}{ccc}\mathit{\lambda ̇}& \mathrm{=}& \\ \mathit{ṗ}& \mathrm{=}& \mathrm{-}\mathit{ϵ}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{(}\mathrm{\Lambda }\mathit{,\lambda }\mathrm{-}\mathit{\varphi }\mathrm{(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathrm{·}\end{array}$(11)It is physically clear that to zeroth order in ϵ, the motion of (p,λ) is simply a linear flow, with p = p(0) and λ(t) = λ(0) + p(0)t. Equations are then invariant with the change of variables p ← p−p(0), λ ← λ−p(0)t, provided the function ϕ(t) is also changed as ϕ(t) ← ϕ(t)−p(0)t. This change of variables means that we integrate out the Keplerian motion of Mars and the new function ϕ(t) represents the difference between the mean longitude of the asteroid and the Keplerian mean longitude of Mars.

Next is the important technical step in the calculation: from Eqs. (11), it is not obvious on which timescale the integrable motion of Mars will be perturbed. It should scale with ϵ but we still do not know the precise scaling at this step of the calculation. Following the method proposed in Lam & Kurchan (2014), we thus introduced a priori an exponent α> 0 and rescale the time according to t′ = ϵ^αt. The variables λ and p were also rescaled according to ${\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{\lambda }{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}$ and ${\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathit{\alpha }}}\mathit{p}{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{a}}}^{\mathrm{\left)}}$ and ${\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{\Lambda }{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}$. The equations for the rescaled variables are $\begin{array}{ccc}\mathit{\lambda ̇}\mathrm{\prime }& \mathrm{=}& \\ \mathit{ṗ}\mathrm{\prime }& \mathrm{=}& \mathrm{-}{\mathit{ϵ}}^{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\alpha }}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\left({\mathrm{\Lambda }}^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }\left(\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}\right)\right)\mathrm{·}\end{array}$(12)To go one step further, we had to find the asymptotic behavior of the two functions $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}$ and $\frac{\mathit{\partial G}}{\mathit{\partial }\mathrm{\Lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}$ in the limit ϵ → 0. This limit is not at all trivial. The averaging principle, which is classically used in celestial mechanics to obtain the secular equations, states that the oscillating function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}$ should be averaged over the distribution of the fast angle ϕ. However, it happens that the function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ is periodic in ϕ, and its average over a uniform distribution of ϕ is zero. The averaging principle only tells us that the semi-major axis is invariant to main order in ϵ. The semi-major axis is an adiabatic invariant because it is a conserved quantity of the Keplerian fast motion, this result is known in celestial mechanics as the theorem of Laplace–Lagrange. To obtain a non trivial variation of the semi-major axis, we have to go to next order, beyond the averaging principle.

For general fast oscillating functions $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$, there is no asymptotic expansion beyond the averaging principle. However, if the function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ satisfies the two following hypotheses:

• The function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ has mixing properties, which means that its correlation time is finite.

• The function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ has zero mean.

then the limit of $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}$ exists and is a Gaussian stochastic process. This is an exact result, not an approximation. The principle used to find the limit of $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ when ϵ goes to zero is called the homogenization process (Gardiner 1985; Pavliotis & Stuart 2008). The reader should be aware that the hypothesis of zero mean is crucial. If it is not satisfied (but this is not the situation we are dealing with), the limit process is not a Gaussian process, and the white noise is only an approximation (Freidlin 1978; Veretennikov 2000; Bouchet et al. 2016).

In the case of Eqs. (12), the first hypothesis is satisfied because of the presence of the asteroid mean longitude ϕ(t). The hypothesis on ϕ is that it has a mixing chaotic dynamics and its correlation time is assimilated to the Lyapunov time τ_ϕ. The second hypothesis is satisfied because of the theorem of Laplace–Lagrange. If $\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathit{\alpha }}}\mathrm{\gg }{\mathit{\tau }}_{\mathit{\varphi }}$, the theory of stochastic averaging (Gardiner 1985, Chap. 8) shows that the function $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}$ is equivalent (it is said to be equivalent in law) to a stochastic process $\begin{array}{ccc}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\left({\mathrm{\Lambda }}^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }\left(\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}\right)\right)& \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\mathrm{~}}& \end{array}$(13)where ξ(t) is the normal Gaussian white noise, ⟨ξ(t)ξ(t′)⟩ = δ(t−t′), and the coefficients B and D are given in terms of the correlation function of $\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}$ (see Gardiner 1985, p. 189) $\begin{array}{ccc}\mathit{B}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \\ \mathit{D}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{t}⟨\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\left({\mathrm{\Lambda }}^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }\mathrm{(}{\mathit{t}}^{\mathrm{)}}\right)\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\left({\mathrm{\Lambda }}^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }\mathrm{(}{\mathrm{0}}^{\mathrm{)}}\right)⟩\mathit{.}\end{array}$(14)It should be noticed that the function G depends on the variable λ−ϕ, and as a result, the coefficients expressed in Eqs. (14) do not depend on λ′. The function $\frac{\mathit{\partial G}}{\mathit{\partial }\mathrm{\Lambda }}$ on the contrary, is not periodic in ϕ and a simple averaging principle is enough to give the equivalent, $\frac{\mathit{\partial G}}{\mathit{\partial }\mathrm{\Lambda }}{}^{\mathrm{(}}\mathrm{\Lambda }^{\mathrm{\prime }}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{-}\mathit{\varphi }{{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{ϵ}}^{\mathit{\alpha }}}^{\mathrm{\left)}}}^{\mathrm{)}}\underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\mathrm{~}}{}^{\mathrm{⟨}}\frac{\mathit{\partial G}}{\mathit{\partial }\mathrm{\Lambda }}^{\mathrm{⟩}}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}$, where the average is done over ϕ.

Altogether, we can give a stochastic equivalent of the system of Eqs. (12) $\begin{array}{ccc}\mathit{\lambda ̇}\mathrm{\prime }& \mathrm{=}& \\ \mathit{ṗ}\mathrm{\prime }& \mathrm{=}& \mathrm{-}{\mathit{ϵ}}^{\mathrm{1}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{\alpha }}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}\left[\mathit{B}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{+}\sqrt{\mathit{D}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{\right)}}\mathit{.}\mathit{\xi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\right]\mathit{.}\end{array}$(15)In the first of Eqs. (15), the perturbation appears at order ϵ^1−α, whereas in the second equation it is ${\mathit{ϵ}}^{\mathrm{1}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{\alpha }}$. As α> 0, the term in the equation for p′ is dominant. We chose the value $\mathit{\alpha }\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}$, because it was the only choice of α that gives a non trivial finite limit to the set of Eqs. (15). We were then able to drop the term ${\mathit{ϵ}}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}{}^{\mathrm{⟨}}\frac{\mathit{\partial G}}{\mathit{\partial }\mathrm{\Lambda }}^{\mathrm{⟩}}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}$ in the first equation. Forgetting the primes for the variables λ′ and p′ we finally write the stochastic equations for the dynamics of Mars perturbed by an asteroid $\begin{array}{ccc}\mathit{\lambda ̇}& \mathrm{=}& \\ \mathit{ṗ}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathit{B}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\sqrt{\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}}\mathit{\xi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}\end{array}$(16)The last set of equations will describe to main order the diffusion of λ over a timescale $\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{2}\mathit{/}\mathrm{3}}}$. This result could seem counterintuitive: a glance at the initial Eqs. (10) could suggest that the perturbation should grow as $\frac{\mathrm{1}}{\mathit{ϵ}}$, or at least as an entire power of $\frac{\mathrm{1}}{\mathit{ϵ}}$. The strange exponent $\frac{\mathrm{2}}{\mathrm{3}}$ comes from the fact that the perturbation over the semi-major axis Λ² is amplified on the mean longitude λ by Keplerian motion. The system (16) has an exact solution because Λ should be considered as a constant variable in the equation. In fact, the initial canonical variable Λ(t) evolves on a much longer timescale of order $\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{2}}}$, and on the timescale $\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{2}\mathit{/}\mathrm{3}}}$, it has thus only very small variations. From the integration of Eq. (16), we deduce that the probability distribution of λ is a Gaussian law. We were then able to compute its variance with respect to the realizations of the noise ξ. The term $\mathrm{-}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathit{B}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}$ will give a deterministic contribution on λ scaling like t². The interesting part comes from the white noise, because it makes the probability distribution of λ spread over time. We have $\mathit{p}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathit{B}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathit{t}\mathrm{+}\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\sqrt{\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}}\mathit{W}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ where W is the standard Brownian motion (also called the Wiener process), and $\begin{array}{ccc}{\mathit{V}}_{\mathit{\lambda }}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{:}\mathrm{=}& {E}\mathrm{\left[}{\mathit{\lambda }}^{\mathrm{2}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathrm{-}{E}\mathrm{\left[}\mathit{\lambda }\mathrm{\right(}\mathit{t}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\\ & \mathrm{=}& {\left(\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\right)}^{\mathrm{2}}\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{\int }\mathrm{\int }\mathrm{d}\mathit{s}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}{E}\mathrm{\left[}\mathit{W}\mathrm{\right(}\mathit{s}\mathrm{\left)}\mathit{W}\mathrm{\right(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\\ & \mathrm{=}& {\left(\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\right)}^{\mathrm{2}}\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{\int }\mathrm{\int }\mathrm{d}\mathit{s}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\begin{array}{c}\inf \end{array}\mathrm{\left(}\mathit{s,}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}}{\left(\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\right)}^{\mathrm{2}}\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}{\mathit{t}}^{\mathrm{3}}\mathit{.}\end{array}$This proves that λ is “superdiffusive”, because the variance is proportional to t³ instead of t for standard Brownian motion. Coming back to the original time, we were then able to define a diffusive timescale and write the mean variation Δλ as $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{diff}}& \mathrm{=}& \\ \mathrm{\Delta }\mathit{\lambda }\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \end{array}$A final remark will conclude this subsection: the asymptotic result (13) shows that the interaction with chaotic asteroids is equivalent to a white noise force acting on the planet, on a timescale much larger than the Lyapunov time τ_ϕ of the asteroid. Equations (14) give the properties of the noise, and is the starting point to study the order of magnitude of the noise amplitude, which is done in the next section. Our computation thus gives a theoretical ground to the model studied in Lam & Kurchan (2014) of an integrable dynamics perturbed by a white noise, and the correct order of magnitude for B and D in Eq. (14).

Comparison with Taylor-Aris dispersion.

Although the results of this section could seem very technical at first sight, the underlying physical mechanism is very simple. The superdiffusion and the scaling ${\mathit{t}}^{\frac{\mathrm{3}}{\mathrm{2}}}$ of the dispersion has been known for a long time in the hydrodynamics community where it is called “Taylor–Aris” dispersion. Let us describe this phenomenon to illustrate the result (17)–(18). Consider a flow in a channel with linear velocity profile v(y) = ny as displayed in Fig. 1. At t = 0, particles are released in the middle of the channel at y = 0. Particles will diffuse along the x-direction and the y-direction. Along the y-direction, we have a simple diffusion and ${}^{\mathrm{⟨}}\mathit{y}^{\mathrm{2}}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{⟩}}\mathrm{\propto }\mathit{t}$. In the x-direction, the flow will amplify the diffusion: particles going up at y feel a velocity ny and are carried forward in the x-direction, whereas those going done at -y feel the velocity -ny and are carried backward in the x-direction. Altogether, if at time t particles have diffused over a range Δy, the dispersion along x will scale like Δx = nΔyt, and since $\mathrm{\Delta }\mathit{y}\mathrm{\propto }\sqrt{\mathit{t}}$ we find that $\mathrm{\Delta }\mathit{x}\mathrm{\propto }\mathit{n}{\mathit{t}}^{\frac{\mathrm{3}}{\mathrm{2}}}$. A diffusion scaling with ${\mathit{t}}^{\frac{\mathrm{3}}{\mathrm{2}}}$ is an example of what we call a superdiffusion, it is illustrated in Fig. 1. The main result we prove in this paper is that the same mechanism happens for planets. The gravitational interaction with the chaotic asteroid causes a diffusion of the semi-major axis a_p. a_p is the equivalent of the y coordinate in Taylor-Aris dispersion. The mean longitude λ circulates at angular velocity – or Keplerian pulsation-n_p. The Keplerian motion amplifies the diffusion on λ exactly as the flow on Fig. 1 amplifies the dispersion on the x-axis.

Fig. 1 Dispersion of particles in a channel at t = 5 (top) and t = 10 (bottom). The flow is indicated with arrows.

5. Orders of magnitude for the diffusion coefficient, the diffusion and Lyapunov timescales

Scaling of the diffusion coefficient.

We aim to evaluate the order of magnitude of the typical diffusion time for the mean longitude of Mars, given by the formula (17) together with the theoretical expression of the diffusion coefficient D(Λ) in Eq. (14). The difficult task comes of course from the diffusion coefficient, because it involves the correlation function of the derivative of G. The full computation is reported in Appendix A. The computation depends only on an averaging over the phase of the asteroid. As the motion of the asteroid is a chaotic function for which we do not have an analytic expression, we must do an hypothesis on how the phase ϕ of the asteroid differs from a simple Keplerian motion. We have to take into account that the Lyapunov time of the phase is given by τ_ϕ. To perform the computation, we assumed that the perturbation of the phase of the asteroid is similar to a Brownian motion $\mathit{W}\left(\frac{\mathit{t}}{{\mathit{\tau }}_{\mathit{\varphi }}}\right)$. However we strongly emphasize that this particular ansatz for the perturbation of the asteroid is chosen to perform analytic calculations, but while the exact result will depend on the particular expression of ϕ, the order of magnitude of the result will not. Our result will thus give the correct order of magnitude of the diffusion coefficient as a function of the correlation time τ_ϕ of ϕ.

The important result of the calculation in Appendix A shows how the diffusion coefficient D(Λ) scales with τ_ϕ. We explained in Sect. 4 that diffusion only occurs if the motion of the asteroid decorrelates fast enough. It is therefore natural to expect that the diffusion coefficient is larger when τ_ϕ is smaller. On the contrary, if the motion of the asteroid is regular, there is no diffusion at all, the diffusion coefficient should be zero for infinite τ_ϕ. A rough order of magnitude for the diffusion coefficient D(Λ) from the formula (A.7) of Appendix A gives, in non dimensional variables$\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{\propto }\frac{{\left|\mathit{G}\right|}^{\mathrm{2}}}{\mathrm{(}{\mathit{n}}_{\mathrm{p}}\mathrm{-}\mathit{\nu }{\mathrm{)}}^{\mathrm{2}}{\mathit{\tau }}_{\mathit{\varphi }}}\mathit{,}$(19)where |G| is the typical order of magnitude of the function G, n_p is the Keplerian pulsation of Mars, ν is the Keplerian pulsation of the asteroid, and τ_ϕ is the Lyapunov time of the asteroid, all expressed in unit of time period of Mars. We can summarize this result by saying that the effect of the chaotic asteroid on the planet Mars is similar to a white noise process acting on the semi-major axis a_p of Mars. In the equation for a_p, the white noise term due to the asteroid gives $\frac{\mathrm{d}{\mathit{a}}_{\mathrm{p}}}{\mathrm{d}\mathit{t}}\mathrm{=}\sqrt{\mathit{D}}\mathit{\xi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(20)where D can be written in dimensional units $\mathit{D}\mathrm{\propto }\mathrm{4}{\left(\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{S}}}\right)}^{\mathrm{2}}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}{\mathit{a}}_{\mathrm{p}}}{\mathrm{|}{\mathit{a}}_{\mathrm{p}}\mathrm{-}\mathit{a}{\mathrm{|}}^{\mathrm{2}}\mathrm{(}{\mathit{n}}_{\mathrm{p}}\mathrm{-}\mathit{\nu }{\mathrm{)}}^{\mathrm{2}}{\mathit{\tau }}_{\mathit{\varphi }}}\mathrm{·}$(21)Combining the expressions (17) and (19) gives the non-dimensional expression for the typical time τ_diff after which Mars looses the memory of its initial longitude ${\mathit{\tau }}_{\mathrm{diff}}\mathrm{=}{\left(\frac{\mathrm{1}}{\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}}{\left(\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}\right)}^{\mathrm{2}}\frac{{\left|\mathit{G}\right|}^{\mathrm{2}}}{\mathrm{(}{\mathit{n}}_{\mathrm{p}}\mathrm{-}\mathit{\nu }{\mathrm{)}}^{\mathrm{2}}{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(22)In order to derive this result, we have chosen the scaling of time, length and mass such that all parameters are of the order of one, except the non-dimensional mass $\mathit{ϵ}\mathrm{=}\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{S}}}$ of the asteroid, and the non-dimensional Lyapunov time of the asteroid τ_ϕ ≃ 10⁴ expressed in the units of time period of Mars. To give the order of magnitude in non dimensional variables, Expression (22) can thus be simplified in ${\mathit{\tau }}_{\mathrm{diff}}\mathrm{=}{\left(\frac{{\mathit{ϵ}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(23)from which we deduce the dimensional expression of τ_diff in seconds ${\mathit{\tau }}_{\mathrm{diff}}\mathrm{\propto }{\left({\left(\frac{{\mathit{m}}_{\mathrm{a}}}{{\mathit{M}}_{\mathrm{S}}}\right)}^{\mathrm{2}}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}}{{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(24)

Assumption of timescale separation.

Our model relies on two hypotheses. First, in order for the white noise limit in Eq. (13) to be valid, we have to assume that a timescale separation exists between the correlation time of the asteroid and the diffusion time. Second, the semi-major axis should be constant to first order in ϵ during the diffusion over the mean longitude. We can summarize the timescale separation as the following inequalities ${\mathit{\tau }}_{\mathit{\varphi }}\mathrm{\ll }{\mathit{\tau }}_{\mathrm{diff}}\mathrm{\ll }{\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{p}}}\mathit{,}$(25)where we have introduced τ_ap a typical time of variation of the semi-major axis.

To check the first inequality, we evaluated the ratio $\frac{{\mathit{\tau }}_{\mathit{\varphi }}}{{\mathit{\tau }}_{\mathrm{diff}}}$ with Eq. (1). We find that $\frac{{\mathit{\tau }}_{\mathit{\varphi }}}{{\mathit{\tau }}_{\mathrm{diff}}}\mathrm{\propto }{\left(\mathit{ϵ}\frac{{\mathit{\tau }}_{\mathit{\varphi }}}{{\mathit{T}}_{\mathrm{p}}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}$ where T_p is the Keplerian period of Mars. With the values of Table 1$\mathit{ϵ}\mathrm{=}\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{S}}}\mathrm{\approx }{\mathrm{10}}^{-9}$ and $\frac{{\mathit{\tau }}_{\mathit{\varphi }}}{{\mathit{T}}_{\mathrm{p}}}\mathrm{\approx }{\mathrm{10}}^{\mathrm{4}}$, then the ratio is smaller than 10^-3 and the first hypothesis is fulfilled. Equation (20) together with (21) also show that a typical diffusion time for the semi-major axis a_p should be ${\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{p}}}\mathrm{\propto }\frac{{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}}{\mathit{D}}$. This gives the order of magnitude τ_ap ≈ 10²² yr and the ratio $\frac{{\mathit{\tau }}_{\mathrm{diff}}}{{\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{p}}}}$ is of the order of $\frac{{\mathit{\tau }}_{\mathrm{diff}}}{{\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{p}}}}\mathrm{\approx }{\mathrm{10}}^{-15}$. As a consequence, we check that our simplified model (9) is consistent.

The model does not conserve energy and angular momentum. The conservation laws could be taken into account in principle by studying the back reaction of Mars on the asteroid at next order in the coupling between Mars and the asteroid. This would add a drift term in the stochastic Eqs. (16) to balance the long term variations of energy and angular momentum. However the mean values of energy and angular momentum are expressed in terms of the semi-major axis a_p which evolves on a timescale τ_ap much longer than the timescale of interest (τ_diff), that’s why the integrals of motion can be considered as constant. As a consequence the drift term is irrelevant on this timescale. We thus conclude that the timescale separation (25) ensures that our model has no bias because of diffusion of energy and angular momentum, on the range of timescales of interest.

Action of several asteroids.

With the results of Appendix A, we are able to give the order of magnitude for the diffusion time of the mean longitude of Mars. On timescales longer that their respective Lyapunov time, the motion of the asteroids can be considered as statistically independent. The perturbations of the asteroids on Mars can thus be considered separately. This allows us to give an estimation of the total perturbation.

Table 2 is the first important quantitative result of this work. It gives an estimation of the time we have to wait before seeing a noticeable influence of the three largest asteroids on the mean longitude of Mars. The diffusion time τ_diff depends strongly both on the mass of the asteroid and on the Lyapunov time τ_ϕ of the asteroid. In particular we see that Ceres and Vesta seem to have similar effects on Mars because Ceres is much larger, but is also much less chaotic than Vesta.

Separation of the trajectories of Mars: super diffusion versus exponential separation.

We also report in the last column of Table 2 the extrinsic Lyapunov time τ_e of planets perturbed by the asteroid belt. It is defined as the inverse of the largest Lyapunov exponent of the dynamical system (10). Physically, it corresponds to the typical time of exponential separation of two close initial conditions in Eq. (10) for the mean longitude of Mars. We do not report details about the computation of τ_e because the method is very similar to the computation of τ_diff and is already explained in the work of Lam & Kurchan (2014). We explain in the following the physical meaning of the two times τ_e and τ_diff.

Consider a number of numerical simulations of the dynamics of Mars and the asteroid belt. If one takes different but very close initial conditions for the planet Mars and the same initial conditions for the asteroids, the trajectories of Mars separate exponentially fast with Lyapunov time τ_e. This means that for times smaller than τ_e we have $\mathrm{\Delta }\mathit{\lambda }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }\mathrm{\Delta }{\mathit{\lambda }}_{\mathrm{0}}\exp \left(\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{e}}}\right)\mathit{,}$(26)where Δλ is the difference in mean longitude between two trajectories. On the contrary, if one simulates the system with different initial conditions for the asteroids, but the same initial condition for the planet Mars, the trajectories of Mars separate as $\mathrm{\Delta }\mathit{\lambda }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }{\left(\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{·}$(27)because of the superdiffusion mechanism. Now, the simulation can start with both different initial conditions for Mars and the asteroids. Table 2 shows that τ_e and τ_diff are of the same order of magnitude. Then the trajectories of Mars will also separate following the power law $\mathrm{\Delta }\mathit{\lambda }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }\mathrm{\Delta }{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}{\left(\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(28)because the superdiffusion mechanism overcomes the exponential divergence for small times t<τ_e. This is illustrated in Fig. 2. This explains why the superdiffusion mechanism is more relevant than exponential separation for the computation of the probability distribution of planetary mean longitudes.

Fig. 2 Separation between the trajectories of Mars with close initial conditions. The superdiffusion mechanism represented by the blue curve gives a separation scaling as a power-law ${{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}^{\mathrm{\left)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}$. The chaos in the system separates two trajectories exponentially fast as displayed by the red curve. Here we have chosen an initial separation Δλ0 ≈ 0.01.

The first conclusion we can give at this stage is that the characteristic timescale over which asteroids affect the planet longitudes is of order of 10 Myr. In order to find this order of magnitude, we have made a computation with circular and coplanar orbits (see Sect. 3). However we stress that taking into account of the eccentricity or the inclination of the planetary orbit would not change this order of magnitude. Moreover, taking into account the secular motion of the eccentricities and inclinations would neither change this order of magnitude. Indeed the mechanism that causes longitude diffusion is the perturbation by the asteroids of the planetary semi-major axis, while on secular timescales the semi-major axis is an adiabatic invariant and does not evolve due to the secular evolution of the planets.

This superdiffusion timescale for the planetary longitude, of order of 10 Myr, is of the same order of magnitude as the Lyapunov time τ_i for the orbital parameters of the inner solar system, as given by Laskar (1989). It is thus a natural question to study the stochastic effects of the asteroids on the secular orbital parameters themselves. This is the subject of next section.

6. Asteroid influence on secular motions

In Sect. 5, we explained the effect of the gravitational influence of asteroids on planetary mean longitudes is a superdiffusion, which means that the standard deviation of the mean longitude grows with time as ${{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}^{\mathrm{\left)}}}^{\frac{\mathrm{3}}{\mathrm{2}}}$. However an incertitude on the exact position of the planet on its orbit does not change the secular equations, and the evolution of the orbital parameters inclinations i and eccentricities e of planetary orbits. The change of the orbital parameters are the most important properties that affect many applications, for instance climate. We are thus naturally led to consider the influence of asteroids on the orbital parameters.

The secular equations with the contribution of asteroids can be described very formally by a Hamiltonian composed of three terms

${\mathit{H}}_{\sec }\mathrm{=}{\mathit{H}}_{\mathrm{LL}}\mathrm{+}{\mathit{H}}_{\mathrm{NL}}\mathrm{+}{\mathit{H}}_{\mathrm{a}}\mathit{,}$where H_LL is the Hamiltonian of Laplace–Lagrange with all quadratic terms in e and i, H_NL gathers the terms of order larger than three in e and i and H_a describes the interaction with the asteroids. We are exactly in the framework described in Sect. 2: the secular Hamiltonian H_LL + H_NL without asteroids is chaotic with an intrinsic Lyapunov time τ_i. The Hamiltonian H_LL + H_a describes the integrable motion of Laplace–Lagrange perturbed by the gravitational interaction with asteroids. The perturbation by asteroids H_a is an external source of chaos in the integrable system, and causes exponential separation of trajectories with close initial conditions over a time τ_e. We want to know which time is larger between τ_i and τ_e. This amounts to determining if the external source of chaos dominates compared to the intrinsic chaos of the secular system. The other related question is to know wether a superdiffusion mechanism occurs on the secular orbital parameters, similar to the one that happens for planetary mean longitudes, and on which timescale τ_sec the superdiffusion mechanism perturbs the secular motion.

As we have done in the last sections, here we will focus only on orders of magnitude. The present section is again somehow technical, and is mostly an extension of the techniques already developed in Sect. 4. The reader is welcome to skip it and go directly to the conclusions presented in the next section.

Let α be the set of all canonical variables eexp(jϖ), $\sin {}^{\mathrm{\right(}}\frac{\mathit{i}}{\mathrm{2}}^{\mathrm{\left)}}\exp \mathrm{\left(}\mathit{j}\mathrm{\Omega }\mathrm{\right)}$ for all planets of the solar system, and their complex conjugates, where $\mathit{j}\mathrm{=}\sqrt{-1}$, ϖ is the longitude of the perihelion, and Ω is the longitude of the ascending node (see e.g., Murray & Dermott 1999, p. 48). α is a vector of dimension 32, taking into account the eight planets. The Laplace–Lagrange system gives $\mathit{\alpha ̇}\mathrm{=}\mathit{A}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathit{\alpha }$(29)where Λ is now the eight-dimensional vector of the canonical momentums associated to planetary mean longitudes. In the secular equations, Λ is a constant because the mean longitudes λ do not appear any more in the Hamiltonian after the averaging procedure. The gravitational interaction with asteroids is described by the Hamiltonian ${\mathit{h}}_{\mathrm{a}}\mathrm{:}\mathrm{=}\sum _{\mathit{p}}\frac{\mathrm{-}\mathrm{𝒢}{\mathit{m}}_{\mathrm{p}}\mathit{m}}{{}_{\mathrm{|}}{r}_{\mathrm{p}}\mathrm{-}{{r}}_{\mathrm{|}}}$, and H_a is the secular contribution of this Hamiltonian. The perturbation induced by asteroids perturbs the dynamics (29) by two different mechanisms:

• 1.

  The asteroid motion gives a white noise term in the equationon the semi-major axis (see Eq. (20)), and thus on${\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{:}\mathrm{=}{\mathit{m}}_{\mathrm{p}}\sqrt{\mathrm{𝒢}{\mathit{M}}_{\mathrm{S}}{\mathit{a}}_{\mathrm{p}}}$ as was explained in Sect. 4. The matrix A of the secular system (29) should not be considered as a constant in the secular equations, because we have to take into account the diffusion over Λ. We show in the following that the vector α of orbital parameters is stochastic and follows a multidimensional superdiffusion, through a mechanism very similar to the superdiffusion of the mean longitudes.

• 2.

  The secular Hamiltonian H_a creates a direct additional term in the equation. It creates terms involving the derivatives of H_a in the secular equations. In principle, those terms could be analyzed by the method developed is Sect. 4. If we write H_a in action-angle variables, we can isolate the average contribution (coming from an averaging procedure over the secular angles), and the diffusive contribution, beyond the averaging principle. The computation of the diffusion coefficient involves the correlation time of the orbital parameters of asteroid’s orbits. We do not know precisely the order of magnitude of this time, that’s why it seems difficult to give the order of magnitude of this diffusion coefficient. We have assumed that this direct diffusion mechanism is smaller than the superdiffusion mechanism coming from diffusion of the planetary semi-major axes.

Let A(t) = A(Λ(t)) and u(t) = e^{− A₀t}α(t), where A₀: = A(0). We thus have u(0) = α(0). With this change of variable, we integrated out the fast motion coming from A₀. u(t) satisfies the differential equation

$\mathit{u̇}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}\mathit{t}}\mathit{\delta A}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}\mathit{t}}\mathit{u}\mathrm{\left(}\mathit{t}\mathrm{\right)}$where we have set δA(t) = A(t)−A₀. We expected δA(t) to be small because it comes from the variation of Λ. We remember that ϵ = m/M_s quantify the asteroid effect on the planets. When ϵ goes to zero, we thus obtain δA → 0. This allowed us to perform a computation to order one in ϵ and we have $\mathit{u}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }\mathit{\alpha }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}\mathit{s}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}\mathit{s}}\mathit{\delta A}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}\mathit{s}}\mathit{\alpha }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(30)The difference δu: = u(t)−α(0) comes from the diffusive behavior of the set of variables α(t). δA is small, and the first order is given by a linear expansion of the matrix A w.r.t Λ. We have $\mathit{\delta A}\mathrm{=}\frac{\mathit{\partial A}}{\mathit{\partial }\mathrm{\Lambda }}\mathit{\delta }\mathrm{\Lambda }$. According to Sect. 4, δΛ(t) is a N-dimensional Brownian motion – N being the number of planets – with magnitude $\mathit{ϵ}\sqrt{\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}}$ where D is now a diagonal matrix gathering the coefficients D_p(Λ_p) for each planet p. We took for δΛ the expression $\mathit{\delta }\mathrm{\Lambda }\mathrm{=}\mathit{ϵ}\sqrt{\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}}\mathit{W}\mathrm{\left(}\mathit{t}\mathrm{\right)}$, which is valid only on the Myr timescale. Finally, we estimated the diffusion time of the set α by computing explicitly the quantity E[δu^†δu] († stands for the complex conjugated and transposed of a matrix or a vector). The computation is shown Appendix B. We have shown that E[δu^†δu] has a component growing like t³ which is the signature of a superdiffusive behavior.

The timescale for the superdiffusion of orbital elements e,i,ϖ,Ω writes in non-dimensional units ${\mathit{\tau }}_{\sec }\mathrm{\simeq }{\left(\frac{\mathrm{1}}{\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}}\begin{array}{c}\sum \end{array}_{\mathit{l,}\mathrm{p}}{\mathit{D}}_{\mathrm{p}}{\left|{\mathit{\alpha }}_{\mathit{l}}^{†}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathit{\alpha }}_{\mathit{l}}\right|}^{\mathrm{2}}\right)}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}}\mathit{,}$(31)where the α_l are the eigenvectors of the matrix of secular frequencies A₀. As one can see, the mechanism of superdiffusion of orbital elements is similar to the superdiffusion of the mean longitudes.

Comparing Expression (17) for τ_diff and Expression (31) for τ_sec, one sees that the nonisochronic parameter $\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}}{\mathit{\partial }\mathrm{\Lambda }}$ is replaced here by the elements of the derivative $\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }\mathrm{\Lambda }}$. The latter matrix is essentially the matrix of the derivatives of secular frequencies. Therefore, the mechanism of superdiffusion for the secular orbital parameters is similar to the one for the mean longitudes, because in both cases it comes from a diffusion of the eigenfrequencies of the integrable motion. τ_diff comes from the diffusion of the Keplerian pulsation n_p, and τ_sec comes from the diffusion of the eigenfrequencies g and f.

To evaluate an order of magnitude of τ_sec, we did the rough assumption that the coefficients ${\mathit{\alpha }}_{\mathit{l}}^{†}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathit{\alpha }}_{\mathit{l}}$ are essentially of the order of magnitude of the secular frequencies, that is, of a few arcseconds/yr. This is the major difference with the superdiffusion of mean longitudes. In Expression (21), the Keplerian frequency is of order of the arc/yr. Because of the difference between the frequencies of Keplerian motion and the secular eigenfrequencies g and f for the solar system, the timescale τ_sec for a superdiffusion of the orbital parameters is expected to be much larger than the diffusion time of the mean longitudes. With dimensional variables, the order of magnitude of the time τ_sec is expressed as ${\mathit{\tau }}_{\sec }\mathrm{\propto }{\left({\left(\frac{{\mathit{m}}_{\mathrm{a}}}{{\mathit{M}}_{\mathrm{S}}}\right)}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\sec }^{\mathrm{2}}}{{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(32)and we find that τ_sec should be typically larger than 10 Gyr. This time is also the order of magnitude of the Lyapunov time τ_e associated with the chaos created by asteroids on the secular system of Laplace–Lagrange.

We thus conclude that τ_e is much larger than τ_i, the Lyapunov time for the eccentricities and inclinations of the internal planets. The asteroids provide thus a very small perturbation to the intrinsically chaotic secular dynamics.

7. Conclusions

In this paper, we have investigated the influence of asteroids on the long-term dynamics of the solar system. With the technique of stochastic averaging, we have shown that the chaotic behavior of the main asteroids in the asteroid belt: Ceres, Vesta, and Pallas, perturbs the long term dynamics of planets, and give a stochastic contribution to the motion of planets. We have explained that the main effect is a superdiffusion of planet mean longitudes, that is, a variance of the mean longitudes growing as ${{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}^{\mathrm{\left)}}}^{\mathrm{3}}$ instead of $\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}$ for normal diffusion. An order of magnitude of τ_diff is given by Eq. (1). We have found that for the planet Mars perturbed by the asteroid belt, the superdiffusion becomes significative after 10 Myr. This order of magnitude should be the same for all planets in the inner solar system. This result thus confirms on a theoretical ground that no accurate planetary ephemerides can be elaborated for times larger that 10 Myr, as was already observed in the numerical simulations of Laskar et al. (2011b).

Within the Myr timescale, the motion over the orbit is averaged out and the dynamics is described by the secular equations. We have studied the effect of asteroids on secular motions and we have found that the mechanism of superdiffusion also happens on the planet eccentricities and inclinations. The time τ_sec that characterizes this superdiffusion, given by Eq. (2), has been estimated to 10 Gyr. The intrinsic chaotic dynamics of the secular equations occurs within a time τ_i of order of 10 Myr. As a consequence the evolution of the planetary eccentricities and inclinations probability distributions, over one Gyr timescale, is completely dominated by the intrinsic secular chaos. At this stage, the conclusion is that asteroids have an influence on secular motion so small that it should change the statistical distributions of eccentricities and inclinations computed with the secular equations by Laskar (2008), only through a tiny perturbation.

Acknowledgments

We thank J. Laskar for his helpful advice during the preliminary stage of this work. The research leading to these results received funding from the European Research Council under the European Union’s seventh Framework Program (FP7/2007-2013 Grant Agreement No. 616811).

References

Appendix A: Computation of the diffusion coefficient

We aim to give an order of magnitude of the Lyapunov time (17) of the longitude λ of Mars feeling the perturbation of a large asteroid of the asteroid belt. The perturbative function gives $\mathit{G}\mathrm{(}\mathrm{\Lambda }\mathit{,\lambda }\mathrm{-}\mathit{\varphi }\mathrm{)}\mathrm{=}\frac{-1}{\left|{\mathrm{\Lambda }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathit{\varphi }\mathrm{)}}\mathrm{-}\mathit{a}\right|}$(A.1)where Λ is the conjugated momentum of the mean longitude λ. In the following, we always assume that Λ²<a. The Fourier expansion of G has been given in the general case by Laskar & Boué (2010), and writes in our simplified model $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \\ {\mathit{f}}_{\mathit{n,q}}& \mathrm{=}& \frac{\mathrm{\left(}\mathrm{2}\mathit{q}\mathrm{\right)}\mathrm{!}\mathrm{(}\mathrm{2}\mathit{n}\mathrm{-}\mathrm{2}\mathit{q}\mathrm{)}\mathrm{!}}{{\mathrm{2}}^{\mathrm{2}\mathit{n}}\mathrm{(}\mathit{q}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}\mathrm{[}\mathrm{(}\mathit{n}\mathrm{-}\mathit{q}\mathrm{)}\mathrm{!}{\mathrm{]}}^{\mathrm{2}}}\end{array}$(A.2)where N is a cut off to stop the Fourier expansion. In our calculations, we took N = 20.

We have the Fourier decomposition of G in the form $\mathit{G}\mathrm{=}\sum _{\mathit{k}}\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{k}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathit{\varphi }\mathrm{)}}$. Our aim is to evaluate the quantity $E{}^{\mathrm{[}}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{(}\mathrm{0}{\mathrm{\left)}}^{\mathrm{\right]}}$ where for simplicity we used the shortcoming G(t) = G(Λ,λ(t)−ϕ(t)). The reader should bear in mind that λ(t) is given by the unperturbed dynamics, because to compute the fast motions we should “freeze” the slow variables. If we freeze Λ in the dynamics of λ, we simply get the Keplerian motion. Therefore λ(t) = λ₀ + nt. The chaotic dynamics of the asteroid is modeled by ϕ(t) = ϕ₀ + ν_at + δϕ(t) where δϕ should account for the chaotic diffusion on a timescale τ_ϕ. We will give its expression later. We thus have $\begin{array}{ccc}E\left[\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\right]& \mathrm{=}& \begin{array}{c}\sum \end{array}_{\mathit{k,}{\mathit{k}}^{\mathrm{\prime }}}\mathit{k}{\mathit{k}}^{\mathrm{\prime }}\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{(}{\mathit{k}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{\ast }}\\ & & \mathrm{\times }E\mathrm{[}{\mathrm{e}}^{\mathrm{i}\mathrm{\left(}\mathit{k\lambda }\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}\mathit{\lambda }\mathrm{(}\mathrm{0}\mathrm{\left)}\mathrm{\right)}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{\left(}\mathit{k\varphi }\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}\mathit{\varphi }\mathrm{(}\mathrm{0}\mathrm{\left)}\mathrm{\right)}}\mathrm{]}\\ & \mathrm{=}& \begin{array}{c}\sum \end{array}_{\mathit{k,}{\mathit{k}}^{\mathrm{\prime }}}\mathit{k}{\mathit{k}}^{\mathrm{\prime }}\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{(}{\mathit{k}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{\ast }}{\mathrm{e}}^{\mathrm{i}\mathit{kt}\mathrm{(}\mathit{n}\mathrm{-}{\mathit{\nu }}_{\mathrm{a}}\mathrm{)}}\mathit{...}\\ & & \mathit{...}E\mathrm{[}{\mathrm{e}}^{\mathrm{i}\mathrm{(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{-}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{)}\mathrm{(}\mathit{k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}\mathrm{)}}\mathrm{]}E\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\delta \varphi }\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{]}\mathit{.}\end{array}$(A.3)We have two averages to perform. The invariant measure for the initial conditions λ₀ and ϕ₀ is the uniform measure over the range [0,2π]. Thus the term $E{}^{\mathrm{[}}\mathrm{e}^{\mathrm{i}\mathrm{(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{-}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{)}\mathrm{(}\mathit{k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}\mathrm{)}}{}^{\mathrm{]}}\mathrm{=}{\mathit{\delta }}_{\mathit{k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}}$. What is then the expected value of e^{− ikδϕ}? For a general chaotic trajectory, it should be a complicated function, decreasing with the timescale τ_ϕ. We chose for δϕ the Brownian motion $\mathit{W}\left(\frac{\mathit{t}}{{\mathit{\tau }}_{\mathit{\varphi }}}\right)$ which has Gaussian statistics. It becomes $\begin{array}{ccc}E\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\delta \varphi }\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{]}& \mathrm{=}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi t}\mathit{/}{\mathit{\tau }}_{\mathit{\varphi }}}}\mathrm{\int }\mathrm{d}\mathit{x}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kx}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}\mathit{t}\mathit{/}{\mathit{\tau }}_{\mathit{\varphi }}}}\\ & \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}}{\mathrm{2}{\mathit{\tau }}_{\mathit{\varphi }}}{\mathit{k}}^{\mathrm{2}}}\mathit{.}\end{array}$Finally, the expression of our correlation function is $E\left[\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\frac{\mathit{\partial G}}{\mathit{\partial \lambda }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\right]\mathrm{=}\begin{array}{c}\sum \end{array}_{\mathit{k}}{\mathit{k}}^{\mathrm{2}}{\left|\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right|}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}}{\mathrm{2}{\mathit{\tau }}_{\mathit{\varphi }}}{\mathit{k}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{i}\mathit{k}\mathrm{(}\mathit{n}\mathrm{-}{\mathit{\nu }}_{\mathrm{a}}\mathrm{)}\mathit{t}}\mathit{.}$(A.6)Finally we integrate this expression over time according to Eq. (14), and we obtain $\mathit{D}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{=}\begin{array}{c}\sum \end{array}_{\mathit{k}}{\mathit{k}}^{\mathrm{2}}{\left|\begin{array}{c}\mathrm{􏽢}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right|}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{(}{\mathit{\nu }}_{\mathrm{a}}\mathrm{-}\mathit{n}{\mathrm{)}}^{\mathrm{2}}{\mathit{\tau }}_{\mathit{\varphi }}}\frac{\mathrm{1}}{{\left(\frac{\mathit{k}}{\mathrm{2}\mathrm{(}{\mathit{\nu }}_{\mathrm{a}}\mathrm{-}\mathit{n}\mathrm{)}{\mathit{\tau }}_{\mathit{\varphi }}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathrm{·}$(A.7)Expressions (A.2) and (A.7) allow us to compute the magnitude of the diffusion coefficient and then to estimate the diffusion timescale of the mean longitude of Mars.

Appendix B: Diffusion of secular variables

We consider again Eq. (30) which gives $\mathit{\delta u}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }\mathit{ϵ}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}\mathit{s}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}\mathit{s}}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\sqrt{\mathit{D}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{\right)}}{\mathit{W}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}\mathit{s}}{\mathit{u}}_{\mathrm{0}}\mathit{,}$(B.1)with a summation on the index p. We want to compute $\begin{array}{ccc}E\left[\mathit{\delta }{\mathit{u}}^{†}\mathit{\delta u}\right]\mathrm{=}{\mathit{ϵ}}^{\mathrm{2}}{”}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}\mathit{s}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}E\mathrm{[}{\mathit{W}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{W}}_{{\mathrm{p}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{s}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{]}}\sqrt{\mathit{D}\mathrm{(}{\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{)}\mathit{D}\mathrm{(}{\mathrm{\Lambda }}_{{\mathrm{p}}^{\mathrm{\prime }}}\mathrm{)}}\mathit{...}& & \\ \mathit{...}{\mathit{u}}_{\mathrm{0}}^{†}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}^{†}{\mathit{s}}^{\mathrm{\prime }}}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}^{†}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}^{†}{\mathit{s}}^{\mathrm{\prime }}}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}\mathit{s}}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{{\mathrm{p}}^{\mathrm{\prime }}}}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}\mathit{s}}{\mathit{u}}_{\mathrm{0}}\mathit{.}& & \end{array}$(B.2)Then we used the fact that A₀ is an anti-hermician operator, and therefore ${\mathit{A}}_{\mathrm{0}}^{†}\mathrm{=}\mathrm{-}{\mathit{A}}_{\mathrm{0}}$. The Brownian motions W_p are independent, thus we have $E{}^{\mathrm{[}}\mathit{W}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{W}}_{{\mathrm{p}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{s}}^{\mathrm{\prime }}{\mathrm{\left)}}^{\mathrm{\right]}}\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{p}\mathrm{-}{\mathit{p}}^{\mathrm{\prime }}\mathrm{)}\inf \mathrm{\left(}\mathit{s,}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}$. We get $\begin{array}{ccc}E\left[\mathit{\delta }{\mathit{u}}^{†}\mathit{\delta u}\right]& \mathrm{=}& {\mathit{ϵ}}^{\mathrm{2}}{”}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}\mathit{s}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\begin{array}{c}\inf \end{array}\mathrm{\left(}\mathit{s,}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{D}\mathrm{(}{\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{)}{\mathit{u}}_{\mathrm{0}}^{†}\\ & & \mathrm{\times }\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}{\mathit{s}}^{\mathrm{\prime }}}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}^{†}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathrm{e}}^{\mathrm{-}{\mathit{A}}_{\mathrm{0}}\mathrm{(}\mathit{s}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathrm{e}}^{{\mathit{A}}_{\mathrm{0}}\mathit{s}}{\mathit{u}}_{\mathrm{0}}\mathit{.}\end{array}$(B.3)We caution the reader that in the preceding expression, A₀ does not commute with its derivative $\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }\mathrm{\Lambda }}$. Let α_i be a set of eigenvectors of A₀ with eigenfrequencies jν_i, such that e^A₀sα_i = e^jν_isα_i. We decomposed u₀ on this set of eigenvectors, u₀ = ∑ c_iα_i. Finally, we introduced the closure relation $\mathrm{Id}\mathrm{=}\sum {\mathit{\alpha }}_{\mathrm{i}}{\mathit{\alpha }}_{\mathrm{i}}^{†}$. We have $\begin{array}{ccc}E\left[\mathit{\delta }{\mathit{u}}^{†}\mathit{\delta u}\right]\mathrm{=}{\mathit{ϵ}}^{\mathrm{2}}\begin{array}{c}\sum \end{array}_{\mathit{i,k,l}}{\mathit{c}}_{\mathit{k}}^{\mathrm{\ast }}{\mathit{c}}_{\mathit{l}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}\mathit{s}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\begin{array}{c}\inf \end{array}\mathrm{\left(}\mathit{s,}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{D}\mathrm{(}{\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{)}{\mathit{\alpha }}_{\mathit{k}}^{†}{\mathrm{e}}^{\mathrm{-}\mathit{j}{\mathit{\nu }}_{\mathit{k}}{\mathit{s}}^{\mathrm{\prime }}}\mathit{...}& & \\ \mathit{...}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}^{†}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathit{\alpha }}_{\mathrm{i}}{\mathrm{e}}^{\mathrm{-}\mathit{j}{\mathit{\nu }}_{\mathrm{i}}\mathrm{(}\mathit{s}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}}{\mathit{\alpha }}_{\mathrm{i}}^{†}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathrm{e}}^{\mathit{j}{\mathit{\nu }}_{\mathit{l}}\mathit{s}}{\mathit{\alpha }}_{\mathit{l}}\mathit{.}& & \end{array}$We have to evaluate a double integral of the form ^!dsds′inf(s,s′)e^{i(ν_l−ν_i)s}e^{− i(ν_k−ν_i)s′}, which does not present any difficulty. The result is that this double integral scales with t unless i = k = l and in that case the growth scales with t³. This means that $E\left[\mathit{\delta }{\mathit{u}}^{†}\mathit{\delta u}\right]\underset{\mathit{\nu t}\mathrm{\gg }\mathrm{1}}{\mathrm{~}}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}}\begin{array}{c}\sum \end{array}_{\mathit{i}}{\left|{\mathit{c}}_{\mathrm{i}}\right|}^{\mathrm{2}}\mathit{D}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathrm{p}}\mathrm{\right)}{\left|{\mathit{\alpha }}_{\mathrm{i}}^{†}\frac{\mathit{\partial }{\mathit{A}}_{\mathrm{0}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathrm{p}}}{\mathit{\alpha }}_{\mathrm{i}}\right|}^{\mathrm{2}}{\mathit{t}}^{\mathrm{3}}\mathit{,}$(B.6)this asymptotic behavior is valid for times t large in front of $\frac{\mathrm{1}}{{\mathit{\nu }}_{\sec }}$ where ν_sec is the typical order of magnitude of the frequencies of the Laplace–Lagrange system, and t is small compared to $\frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{2}}\parallel \mathit{D}\left(\mathrm{\Lambda }\right)\parallel }$. The latter assumption is easily satisfied because ϵ< 10^-9 and the former assumption is satisfied for times larger than a Myr.

Appendix C: Notations

All Tables

All Figures

	Fig. 1 Dispersion of particles in a channel at t = 5 (top) and t = 10 (bottom). The flow is indicated with arrows.
In the text

Fig. 2 Separation between the trajectories of Mars with close initial conditions. The superdiffusion mechanism represented by the blue curve gives a separation scaling as a power-law ${{}^{\mathrm{\right(}}\frac{\mathit{t}}{{\mathit{\tau }}_{\mathrm{diff}}}^{\mathrm{\left)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}$. The chaos in the system separates two trajectories exponentially fast as displayed by the red curve. Here we have chosen an initial separation Δλ0 ≈ 0.01.

In the text