Mercury’s Chaotic Secular Evolution as a Subdiffusive Process

Mercury’s orbit can destabilize, generally resulting in a collision with either Venus or the Sun. Chaotic evolution can cause g 1 to decrease to the approximately constant value of g 5 and create a resonance. Previous work has approximated the variation in g 1 as stochastic diffusion, which leads to a phenomological model that can reproduce the Mercury instability statistics of secular and N-body models on timescales longer than 10 Gyr. Here we show that the diffusive model significantly underpredicts the Mercury instability probability on timescales less than 5 Gyr, the remaining lifespan of the solar system. This is because g 1 exhibits larger variations on short timescales than the diffusive model would suggest. To better model the variations on short timescales, we build a new subdiffusive phenomological model for g 1. Subdiffusion is similar to diffusion but exhibits larger displacements on short timescales and smaller displacements on long timescales. We choose model parameters based on the behavior of the g 1 trajectories in the N-body simulations, leading to a tuned model that can reproduce Mercury instability statistics from 1–40 Gyr. This work motivates fundamental questions in solar system dynamics: why does subdiffusion better approximate the variation in g 1 than standard diffusion? Why is there an upper bound on g 1, but not a lower bound that would prevent it from reaching g 5?

The inherent unpredictability of chaotic dynamical systems like the Solar System makes it necessary to describe the long-term evolution statistically.Statistical descriptions of how the phase space distribution of an ensemble of systems evolves over time are relatively welldeveloped for simple area-preserving planar maps and 2 degree-of-freedom systems (e.g., Mackay et al. 1984;Meiss 1992;Zaslavsky 2002), though even in this simplest case there remain important unresolved questions (e.g., Meiss 2015).Given the lack of theoretical understanding of chaotic transport in systems with a larger number of degrees of freedom, phenomenological models such as those discussed in this paper can serve as useful tools for describing complex real-world systems and can potentially provide clues for better understanding the underlying chaotic dynamics.
Several diffusive phenomenological models have been used to approximate the Solar System dynamics and predict the probability of Mercury instability events as a function of time.Woillez & Bouchet (2020) analyzed the simplified secular Hamiltonian of Batygin et al. (2015), which considers Mercury to be massless and approximates the other planets as quasiperiodic, and they identified the slowly varying component of this Hamiltonian as driving Mercury's dynamics.They approximated the dynamics as diffusive with constant diffusivity, a reflecting upper boundary, and an absorbing lower boundary that signifies Mercury instability events.
Later, Mogavero & Laskar (2021) speculated that the diffusive model might apply to the long-term variation of g 1 itself (see section 2.1 for the definition of g i as the Solar System's secular eigenfrequencies).They applied the diffusive model using g 5 , which is effectively constant (Hoang et al. 2021), as the absorbing lower boundary at which Mercury instability occurs (Fig. 1).They tuned the upper boundary and diffusivity to produce a reasonable approximation to the Mercury instability probability statistics of a secular model on timescales longer than 10 Gyr, when at least 4% of the simulations have gone unstable.The 10 Gyr timescale is longer than the future lifespan of the Sun, so their model can only be interpreted as an abstract investigation of the Solar System dynamics.
Recently, Brown & Rein (2023) compared the g 1 diffusive model (Mogavero & Laskar 2021) to the 5 Gyr N -body simulations performed by Laskar & Gastineau (2009).They found what appeared to be reasonable correspondence; however, in their Figures 2-4, they plotted one minus the probability of a Mercury instability event, which obscures the difference between small probabilities spanning orders of magnitude.Brown & Rein (2023) also compared the g 1 diffusive model to their own N -body simulations with general relativity artificially either fully or partially disabled.Similar to the N -body simulations in this paper, they approximated general relativity as a simple potential.Their plots suggest that the diffusive model provides a qualitatively reasonable approximation of the evolution of the Mercury instability probability when it is above ≈5%, but the plots obscure smaller probabilities.In summary, the g 1 diffusive model (Fig. 1) can be tuned to approximate the Mercury instability probability produced by more complex secular or N -body models, as long as the Mercury instability probability exceeds ≈5%.In this paper we will show that the g 1 diffusive model significantly underpredicts the Mercury instability probability over the next 5 Gyr period in which the Sun will remain on the main sequence.We find that the discrepancy primarily results from the diffusive model producing variations of g 1 that are too small on timescales less than ∼0.3 Gyr.To better model the short-time variations of the g 1 trajectories, we propose a g 1 subdiffusive model (see Henry et al. 2010, for an introduction to subdiffusion).The subdiffusive model is consistent with the work of Hoang et al. (2021), who found that g 1 subdiffuses by fitting a power law to the standard deviation of g 1 in a large ensemble of secular models as a function of time.We fit the parameters of our subdiffusive model to the mean square displacement of g 1 , the probability density function (pdf) of g 1 from N -body simulations, and the fraction of N -body simulations that reach the instability threshold.This produces a fiveparameter stochastic subdiffusive model that closely approximates N -body Mercury instability probabilities as small as ∼ 10 −4 (the lowest value we are able to estimate with N -body results due to limited sample size) and as high as ∼ 0.5 (the highest value we are able to estimate with N -body results due to limited run lengths).The g 1 subdiffusive model is successful at reproducing a variety of statistics from the N -body code, suggesting that despite its simplicity, the model captures important aspects of the relevant dynamics.We calculate g 1 statistics from the 2750-member Fix dt ensemble of simulations produced by Abbot et al. (2023).The simulations contain all 8 Solar System planets and no moons, asteroids, or comets.The simulations use the WHFAST integration scheme (Rein & Tamayo 2015) from the REBOUND N -body code (Rein & Liu 2012), which is a Wisdom-Holman scheme (WH) (Wisdom & Holman 1991).The only parameterized physics scheme is an approximation of general relativity with a modified position-dependent potential (Nobili & Roxburgh 1986), which is implemented as the gr potential scheme in REBOUNDx (Tamayo et al. 2020).We initialized the simulations with Solar System conditions on February 10, 2018 taken from the NASA Horizons database and then added a uniform grid of perturbations to Mercury's x-position, each separated by 10 cm.We used a fixed time step of √ 10 ≈ 3.16 days, which we demonstrated was sufficiently small to produce converged Mercury instability statistics (Abbot et al. 2023), and ran the simulations for 5 Gyr.The trajectories of eccentricity and g 1 are nearly identical among the simulations over the first 138 Myr of the simulation and then begin to noticeably separate from each other on the order of ≈1%.This is longer than the typical quoted value of ∼50 Myr for Solar System orbital calculations to diverge (e.g., Laskar et al. 2011;Zeebe 2017), possibly because Mercury's eccentricity and g 1 take longer to diverge than other variables.
The g i frequencies are defined through the following eigenfunction expansion (Murray & Dermott 1999, Ch. 7): where the index i ranges over the planets i = 1, 2, . . ., 8, e i denotes the eccentricity of each planet, and ϖ i denotes the longitude of perihelion.The values (M i,k ) 1≤i,k≤8 are the coefficients in the eigenfunction expansion, and the terms α k = g k t + β k describe the angles of the oscillations.
The first-order approximation with constant parameter values (1)-( 2) is accurate over the shortest period of oscillation min i 2π/g i , but degrades thereafter.Indeed, when we fit the approximation to the output of a more complex N -body simulation, the values M ik , g k , and β k are all slowly changing as a function of time.Despite the imprecision, the first-order approximation leads to qualitative insights, as it correctly indicates the possibility for destabilization when two of the g i frequencies approach resonance.
To calculate g 1 , we use the frequency modified fourier transform routine ( Šidlichovskỳ & Nesvornỳ 1997) from the celmech package (Hadden & Tamayo 2022), which requires a number of output times that is a power of 2. We find that 2 10 output times are sufficient to obtain stable g i estimates from data, so we divide the Fix dt ensemble data into blocks of 10.24 Myr, with each block containing 2 10 output times.Our calculated values of g 1 are therefore the average value over each 10.24 Myr segment.We apply frequency modified fourier transform to e 1 exp(iϖ 1 ) to calculate the four Fourier modes with the highest amplitude for a given 10.24 Myr segment and chose g 1 to be the g i closest to g 1 from the previous segment.

Model used to calculate 40 Gyr instability statistics
As a new contribution of this paper, we perform 1000 extensions of the Fix dt simulations and make them publicly available at https://archive.org/download/LongRun.These extensions have exactly the same parameters as the Fix dt simulations described above, but we run them for a total of 40 Gyr.As in Abbot et al. (2021Abbot et al. ( , 2023)), we define a Mercury instability event as occurring when Mercury passes within 0.01 AU of Venus and stop the simulations at that point.In Abbot et al. (2023), we showed that after 10 12 time steps (10 Gyrs), roundoff relative error is of order 10 −5 for the semimajor axis and order 10 −9 for the energy.Roundoff error is growing as ∼ t 0.5 at 10 Gyr, so it should remain small for the 40 Gyr simulations we performed.

Ensemble of ensembles used to compute 5 Gyr instability statistics
To compute Mercury instability statistics on a timescale of less than 5 Gyr, we use the ensemble of N -body ensembles constructed by Abbot et al. (2023), which includes the 2501-member Laskar & Gastineau (2009) ensemble, the 1600-member Zeebe (2015a) ensemble, as well as both the 2750-member Var dt and 2750-member Fix dt ensembles of Abbot et al. (2023), for a total of 9601 members.

Diffusive and subdiffusive models
Both the diffusive and subdiffusive models are defined using fractional Brownian motion.Fractional Brownian motion is a mean-zero Gaussian process which we denote by W (t) at each time t ≥ 0. It is defined as the unique mean-zero Gaussian process which starts from W (0) = 0 and has increments satisfying (3) at all times t, s.Here, α ∈ (0, 1) is the Hurst parameter, which is α = 1/2 for a standard diffusion, whereas α ∈ (0, 1/2) for a subdiffusion and α ∈ (1/2, 1) for a superdiffusion (Henry et al. 2010, Sec. 1.2).As a result of the scaling relation (3), a subdiffusion exhibits larger variations over short timescales than a standard diffusion.
The simple diffusion model of Mogavero & Laskar (2021, Sec. 8.2) states that the g 1 value at time t (units of Gyrs) satisfies: where W (t) is a standard Brownian motion, r > 0 is a scaling factor, g 1 (0) = g 0 1 is the initial condition, g ↑ 1 > g 0 1 is a reflecting upper boundary condition, and is the first time the process hits the upper boundary.The process advances forward until hitting the lower boundary g ↓ 1 at a random time Then, instability occurs in the model (Fig. 1).Mogavero & Laskar (2021, Sec. 8.2) motivate their use of diffusion and the reflecting upper boundary by observing that the pdf of g 1 has Gaussian tail behavior at low g 1 values, but drops off sharply at high g 1 values (Hoang et al. 2021).Woillez & Bouchet (2017) provide additional motivation, by referencing the theory of slow-fast dynamical systems (Gardiner 2009, Ch. 6), in which the evolution of a slow variable can sometimes be modeled as a diffusive SDE.However, Hoang et al. (2021) argue that the time evolution of g 1 under the secular equations matches a subdiffusion more closely than a standard diffusion.
Next, we propose a more general model in which we allow W (t) to be a fractional Brownian motion with a Hurst parameter α not necessarily equal to 1/2.In the general model, the g 1 value at time t is given by the same equations ( 4)-( 6) but the standard Brownian motion is replaced by a fractional Brownian motion.As before, the process advances forward until hitting the lower boundary g ↓ 1 .Then, instability occurs.To simulate from the diffusive and subdiffusive models, we do not invoke equation ( 4) directly, because this leads to a discretization error of size O(δ α ) for a time step δ (McGlaughlin & Chronopoulou 2016).Therefore, we pursue an alternative, more efficient discretization strategy (Wada & Vojta 2018).First, we generate a random vector containing the values of the fractional Brownian motion at discrete output times W = W (δ), W (2δ), . . ., W (N δ) , using the algorithm of Dietrich & Newsam (1997) as implemented in the stochastic package for python (Flynn 2022).Next, we set g 1 (t = 0) = g 0 1 and apply the recursive update formula This formula directly invokes a fractional Brownian motion that is reflected off an upper boundary.The formula yields the exact correct distribution without any discretization error for α = 1/2 (Doob 1953, pg. 393) and closely approximates the distribution for all α ∈ (0, 1).We apply equation ( 7) with a time step δ = 0.01 Gyr in our simulations, and our code is available at https://knowledge.uchicago.edu/record/10351.

Problems with the g 1 diffusive model
In this subsection, we point out several issues with the g 1 diffusive model.As a starting point, the g 1 diffusive model is not effective at predicting Mercury instability probabilities over physically realistic timescales.On timescales longer than ∼10 Gyr, the model matches with the instability probabilities from secular model simulations (Mogavero & Laskar 2021) and from our 40 Gyr  N -body simulations (Fig. 2).However, on the 5 Gyr timescale of the future of the Solar System, the diffusive model underpredicts Mercury instability events by a factor of 3-1000 (Fig. 2).One possible explanation for the fact that the g 1 diffusion produces too few Mercury instability events on shorter timescales (<10 Gyr) is that subdiffusion (Henry et al. 2010) better approximates the chaotic evolution of g 1 .The main characteristic of subdiffusion is that it exhibits larger displacements on short timescales and smaller displacements on long timescales than diffusion.To investigate this idea quantitatively, we calculate the mean square displacement ⟨|∆g 1 | 2 ⟩ across a range of time offsets ∆t and consider a scaling relationship where α ∈ (0, 1) is the Hurst parameter that can be chosen to match the data.As explained in section 2.4, α = 1/2 corresponds to standard diffusion, whereas α ∈ (0, 1/2) corresponds to subdiffusion and α ∈ (1/2, 1) corresponds to superdiffusion.The top panel of Fig. 3 shows that α for the N -body model is much less than 1/2 (we will show below that α ≈ 0.26), and the diffusive model produces a mean square displacement for g 1 that is too small on timescales less than ∼0.3 Gyr.
Next, let us compare the pdf of g 1 between the N -body model and the diffusive model (Fig. 3 bottom).The diffusive model as tuned by Mogavero & Laskar (2021) overpredicts by an order of magnitude the probability that g 1 has a value less than 5 ′′ yr −1 .The unrealistic low values of g 1 occur in the diffusion model because the lower boundary on g 1 is set to be g 5 = 4.257 ′′ yr −1 .However, while the main physical mechanism for a Mercury instability event is a g 1 -g 5 resonance (Batygin et al. 2015), the g 5 resonance might cause non-diffusive behavior as g 5 is approached.For example, the g 1 trajectories that lead to Mercury instability events in Figure 5 of Mogavero & Laskar (2021) often show large, erratic jumps from a value of ∼ 5 ′′ yr −1 to g 5 as the instability event occurs.When we allow the lower boundary on g 1 to be a free parameter in the best-fit diffusive model, the fit for small g 1 improves by about a factor of two, but a large discrepancy remains in order for the model to match the instability probability.To conclude this subsection, the g 1 diffusive model has the following defects: 1.It underpredicts the Mercury instability probability on timescales less than 10 Gyr.
2. It produces too small variations in g 1 on timescales less than ∼0.3 Gyr.
3. It leads to a scaling of the mean square displacement ⟨|∆g 1 | 2 ⟩ with the time offset ∆t that does not fit the N -body simulations.
4. As formulated by Mogavero & Laskar (2021), it assumes that g 1 must diffuse all the way to g 5 to produce an instability event, which is not the case.

An improved g 1 subdiffusive model
We now apply our new subdiffusive model that addresses the limitations of the g 1 diffusive model.To begin, we fix the initial value g 0 1 = 5.60 ′′ yr −1 , because this is the g 1 value at which the N -body simulations start to diverge.Then, we make the modeling assumption that the mean square displacement scales as a power law over short timescales ∆t.This modeling assumption fits the mean square displacement data for the N -body simulations with α = 0.26 and r = 0.16 ′′ yr −1 (see Fig. Next, we impose an absorbing lower boundary g ↓ 1 and a reflecting upper boundary g ↑ 1 .These boundary conditions fit the probability density function (pdf) of g 1 from the N -body simulations for g ↓ 1 = 4.85 ′′ yr −1 , g ↑ 1 = 5.79 ′′ yr −1 .See the bottom panel of Fig. 4. We also tried a soft upper boundary similar to the potential energy of a coiled spring; although it improved the pdf somewhat, it did not significantly improve the overall cost, so we do not include it here to reduce the number of free parameters.
As the main result of our modeling efforts, the g 1 subdiffusive model yields an accurate approximation of the Mercury instability probability statistics for the N -body model both from 1-5 Gyr and from 5-40 Gyr (Fig. 5).The ability of the g 1 subdiffusive model to reproduce the Mercury instability statistics on timescales less than 10 Gyr is primarily due to larger variations in g 1 on timescales less than 0.3 Gyr and represents a significant advance beyond the previous g 1 diffusive model.

Tuning Algorithm
To fine-tune the values of α, r, g ↓ 1 , and g ↑ 1 in the subdiffusive model presented in section 3.2, we introduce the cost function where the components of the cost function are defined as follows: • C prob is the sum of differences in the instability probabilities p for the subdiffusive model and the N -body model.The sum is performed over all times t = 0, 0.01, . . ., 40 Gyr • C log is the sum of differences in the weighted log instability probabilities log p/t for the subdiffusive model and the N -body model.The sum is performed over all times t = 0, 0.01, . . ., 40 Gyr for which p > 0 in the N -body model.
• C disp is the sum of differences in the weighted log mean square displacements log⟨|∆g 1 | 2 ⟩/∆t for the subdiffusive model and the N -body model.The sum is performed over all time lags ∆t = 0.01, . . ., 5.0 Gyr.
In summary, the cost function penalizes disagreements in the long-time instability probabilities (C prob component), the short-time instability probabilities (C log component), the mean square displacements (C disp component), and the pdfs of g 1 position (C hist component).These four components are scaled by the appropriate orders of magnitude so that they contribute similarly to the overall cost.
To minimize the overall cost C, we perform a grid search over possible parameter values and calculate probabilities based on 10 6 realizations of the stochastic model, leading to the results listed in  Laskar (2021).The "diffusive" and "subdiffusive" models fix the starting parameter g 0 1 = 5.60 ′′ yr −1 and adapt the parameters r, g ↓ 1 , and g ↑ 1 to minimize the cost.The "diffusive" model fixes α = 0.50, while the "subdiffusive" model adapts α to minimize the cost.
subdiffusive Hurst parameter α = 0.26 dramatically reduces the cost.Compared to the original Mogavero & Laskar (2021) model, the best subdiffusive model has a 3× smaller cost.The subdiffusive model leads to the most improvement in C log , which is weighted toward earlier times, and C disp , as expected since it was motivated by the poor fit of the diffusive model to the mean square displacement.See Table 1 for a list of the five converged parameter values for the subdiffusive model (four of which are free parameters).

DISCUSSION
The success of the g 1 subdiffusive model at approximating so many characteristics of the N -body model sets new challenges for the planetary dynamics community: First, why does g 1 subdiffuse rather than diffuse?Second, what is the physical cause of the restoring upper boundary on g 1 ?Given that the upper boundary prevents g 1 from resonating with g 2 ≈ 7.45 ′′ yr −1 , why isn't there a restoring lower boundary on g 1 that prevents it from reaching g 5 and thereby prevents Mercury instability events?
Clues to the answers to these questions may lie in the properties of a conservative, Hamiltonian system.The phase space of such systems are generically comprised of a mixture of regular and chaotic trajectories.The dynamics of two degree-of-freedom Hamiltonian systems are equivalent to those of area-preserving maps via the construction of Poincarè return maps.The mixed phase space of two-dimensional area-preserving maps consists of elliptic periodic orbits surrounded by KAM curves that constitute "islands" embedded in a chaotic "sea" (e.g., Lichtenberg & Lieberman 1992).The KAM curves of these regular islands form strict barriers for trajectories in the chaotic regime.The "stickiness" of the borders of these regular islands could lead to behavior that can effectively be described as subdiffusion, similar to how diffusion on fractal materials can lead to subdiffusion (Henry et al. 2010).In higher dimensions, the surviving KAM tori of regular trajectories would no longer impose strict topological constraints on the phase space accessible to chaotic orbits, but may still limit the range of excursions in a way that can be described by a soft, spring-like boundary.Of course this discussion is highly speculative, and more detailed research is needed to satisfactorily explain the dynamical properties of the Solar System we have identified in this paper.
The g 1 subdiffusive model is an improvement over the g 1 diffusive model, but it does not produce identical behavior to the N -body model.It is useful for improving understanding of the N -body model, not for replacing it.For example, the N -body trajectories show intermittency, transitioning from sustained quiescent periods to sustained active periods (Fig. 6).Periods of relative quiescence could be associated with proximity to islands of regularity.Finally, it is important to remember that an N -body code is not a perfect representation of reality either.The subdiffusive model fits Mercury instability statistics from the N -body model much better than the diffusive model on timescales less than 5 Gyr, but we do not have enough N -body Mercury instability events on timescales less than 2 Gyr to thoroughly test the subdiffusive model on the shortest timescales.One option would be to generate large enough ensembles (possibly with 10 6 members) using a high-order secular code (e.g., Mogavero & Laskar 2021) to estimate the probability of a Mercury instability event on these short timescales.Alternatively, more N -body Mercury instability examples on shorter timescales could be obtained using Diffusion Monte Carlo (Ragone et al. 2018;Webber et al. 2019;Ragone & Bouchet 2021;Abbot et al. 2021) or action minimization (E et al. 2004;Plotkin et al. 2019;Woillez & Bouchet 2020;Schorlepp et al. 2023) rare event schemes, aided by machine learning predictor functions (Ma & Dinner 2005;Chattopadhyay et al. 2020;Finkel et al. 2021;Miloshevich et al. 2022;Finkel et al. 2023).

CONCLUSIONS
The main conclusions of this paper are: 1.The g 1 diffusive model significantly underpredicts the Mercury instability probability relative to an N -body model on timescales less than 5 Gyr, which is the physically relevant timescale for the future of the Solar System.The underprediction results from the fact that the g 1 diffusive model produces too small variations of g 1 on timescales less than ∼0.3 Gyr.
2. We are able to fit N -body Mercury instability statistics on timescales of less than 5 Gyr as well as longer timescales using the g 1 subdiffusive model.We tune the model using the mean square displacement of g 1 , probability density function (pdf) of g 1 from N -body simulations, and the fraction of N -body simulations that reach the instability threshold.

Figure 1 .
Figure 1.Schematic diagram of the g1 diffusion model for Mercury instability events.g1 is initialized at g 0 1 and diffuses.g ↑ 1 is a hard, reflecting upper boundary and g ↓ 1 is an absorbing lower boundary.If g1 reaches g ↓ 1 , a Mercury instability event occurs.
Model used to calculate g 1 statistics

Figure 3 .
Figure 3. Top: Mean square displacement of g1 as function of time offset ∆t for the N -body model (black), the Mogavero & Laskar (2021) diffusive model (magenta), and the best fit diffusive model (yellow).Bottom: Probability density function (pdf) of g1 from the N -body model (black) and the two diffusive models (magenta and yellow).The pdfs include all g1 data over 5 Gyr from the 2750-member N -body ensemble and an equivalent generated diffusive ensemble.The values of g5 and g 0 1 for the Mogavero & Laskar (2021) diffusive model are indicated on the plot.

Figure 4 .
Figure 4. Top: Mean square displacement of g1 as function of time offset (∆t) for the N -body model (black) and the subdiffusive model described in Secs.3.2-3.3(purple).Bottom: Probability density function of g1 from the N -body model (black) and the subdiffusive model (purple).The pdfs include all g1 data over 5 Gyr from the 2750-member N -body ensemble and an equivalent generated subdiffusive ensemble.The values of g5 and g 0 1 are indicated on the plot.

Figure 5 .
Figure 5. Probability that Mercury's orbit becomes unstable as a function of time for the N -body model (black) and the subdiffusive model described in Secs.3.2-3.3(purple).This figure should be compared with Fig. (2), which gives the equivalent information for the Mogavero & Laskar (2021) diffusive model.

Figure 6 .
Figure 6.Left: A sample trajectory of g1 from the N -body model that demonstrates intermittency: sustained quiescent periods alternate with sustained active periods.Right: A sample trajectory of g1 from the subdiffusive model for comparison.

Table 1 .
Parameters and their values for the diffusive and subdiffusive models.The parameter values for the ML21 diffusive model are taken fromMogavero & Laskar (2021), but rewritten using our variables.

Table 2 .
The table shows how the diffusive model ofMogavero &  Laskar (2021)can be slightly improved by optimizing the parameters r, g ↓ 1 , and g ↑ 1 .In contrast, adopting the

Table 2 .
Cost function C for different diffusive and subdiffusive models.The ML21 model is taken from Mogavero &