11institutetext: Astronomie et Systèmes Dynamiques, Institut de Mécanique
Céleste et de Calcul des Éphémérides
CNRS UMR 8028, Observatoire de Paris, Université PSL, Sorbonne Université, 77
Avenue Denfert-Rochereau, 75014 Paris, France
11email<EMAIL_ADDRESS>
# Chaotic diffusion of the fundamental frequencies
in the Solar System

Nam H. Hoang  Federico Mogavero  Jacques Laskar

The long-term variations of the orbit of the Earth govern the insolation on
its surface and hence its climate. The use of the astronomical signal, whose
imprint has been recovered in the geological records, has revolutionized the
determination of the geological time scales (e.g. Gradstein & Ogg 2020).
However, the orbital variations beyond 60 Myr cannot be reliably predicted
because of the chaotic dynamics of the planetary orbits in the Solar System
(Laskar 1989). Taking into account this dynamical uncertainty is necessary for
a complete astronomical calibration of geological records. Our work addresses
this problem with a statistical analysis of 120 000 orbital solutions of the
secular model of the Solar System ranging from 500 Myr to 5 Gyr. We obtain the
marginal probability density functions of the fundamental secular frequencies
using kernel density estimation. The uncertainty of the density estimation is
also obtained here in the form of confidence intervals determined by the
moving block bootstrap method. The results of the secular model are shown to
be in good agreement with those of the direct integrations of a comprehensive
model of the Solar System. Application of our work is illustrated on two
geological data: the Newark-Hartford records and the Libsack core.

###### Key Words.:

Chaos – Diffusion - Celestial mechanics – Methods: statistical – Solar System
– Milankovitch cycles

##  1 Introduction

Milankovitch (1941) hypothesized that some of the past large changes of
climate on the Earth originate from the long-term variations of its orbital
and rotational elements. These variations are imprinted along the
stratigraphic sequences of the sediments. Using their correlations with an
orbital solution (Laskar et al. 2004, 2011a; Laskar 2020), some of the
geological records can be dated with precision. This method, named
astrochronology, has becomes a standard practice in the stratigraphic
community, and proves to be a powerful tool to reconstruct the geological
timescale (e.g Gradstein et al. 2004, 2012; Gradstein & Ogg 2020).

The climate rhythms found in the geological records are directly related to
the Earth’s precession constant and to the fundamental secular frequencies of
the Solar System: the precession frequencies $(g_{i})_{i=1,8}$ of the planet
perihelia and the precession frequencies $(s_{i})_{i=1,8}$ of their ascending
nodes. The evolution of these fundamental frequencies is accurately determined
up to 60 Myr (Laskar et al. 2004, 2011a; Laskar 2020). Beyond this limit, even
with the current highest precision ephemerides, it is hopeless to obtain a
precise past history of the Solar System simply by numerical integration. This
limit does not lie in the precision of determination of the initial conditions
but originates in the chaotic nature of the Solar System (Laskar 1989, 1990;
Laskar et al. 2011b). However, because the astronomical signal is recorded in
the geological data, it appears to be possible to trace back the orbital
variations in the geological record and thus to constrain the astronomical
solutions beyond their predictability horizon (Olsen & Kent 1999; Ma et al.
2017; Olsen et al. 2019). Nevertheless, a deterministic view of the Solar
System is no longer reliable beyond 60 Myr, and a statistical approach should
be adopted. Geological constraints likewise should be retrieved in a
statistical setting. In this spirit, a recent Bayesian Markov Chain Monte
Carlo approach has been proposed to provide geological constraints on the
fundamental frequencies (Meyers & Malinverno 2018). For such a Bayesian
approach to give any meaningful constraint, proper prior distributions of the
fundamental secular frequencies are required and therefore a statistical study
of the orbital motion of the Solar System planets is needed. This constitutes
the motivation for the present study.

Laskar (2008) performed the first statistical analysis of the long-term
chaotic behavior of the planetary eccentricities and inclinations in the Solar
System. Mogavero (2017) reconsidered the problem from the perspective of
statistical mechanics. Our study is a follow-up of (Laskar 2008). We study
fundamental frequencies instead of orbital elements because they are more
robust and are closer to the proxies that can be traced in the geological
records. This study is based on the numerical integrations of 120 000
different solutions of the averaged equations of the Solar System over 500
Myr, among which 40 000 solutions were integrated up to 5 Gyr. The initial
conditions of the solutions are sampled closely around a reference value, that
is compatible with our present knowledge of the planetary motion.

##  2 Dynamical model

###  2.1 Secular equations

We use the secular equations of motions of (Laskar 1985, 1990, 2008, and
references therein). They are obtained by series expansion in planetary
masses, eccentricities and inclinations, and through analytical averaging of
second order over the rapidly-changing mean longitudes of the planets. The
expansion is truncated at second order with respect to the masses and degree 5
in eccentricities and inclinations. The equations include corrections from
general relativity and Earth-Moon gravitational interaction. This leads to the
following system of ordinary differential equations:

$\frac{d\omega}{dt}=\sqrt{-1}\\{\Gamma+\Phi_{3}(\omega,\bar{\omega})+\Phi_{5}(\omega,\bar{\omega})\\},$
(1)

where $\omega=(z_{1},\dots,z_{8},\zeta_{1},\dots,\zeta_{8})$ with
$z_{k}=e_{k}\exp(\varpi_{k})$, $\zeta_{k}=\sin(i_{k}/2)\exp(\Omega_{k})$. The
variable $\varpi_{k}$ is the longitude of the perihelion, $\Omega_{k}$ is the
longitude of the ascending node, $e_{k}$ is eccentricity and $i_{k}$ is
inclination. $\Phi_{3}(\omega,\bar{\omega})$ and
$\Phi_{5}(\omega,\bar{\omega})$ are the terms of degree 3 and 5. The $16\times
16$ matrix $\Gamma$ is the linear Laplace-Lagrange system, which is slightly
modified to make up for the higher order terms in the outer Solar System. With
adapted initial condition, the secular solution is very close to the solution
of direct integration over 35 Myr (Laskar et al. 2004). The major advantage of
the secular system over direct integration is speed. Numerically integrating
averaged equations is 2000 times faster than non-averaged ones, due to the
much larger step size: 250 years instead of 1.8265 days. It is thus desirable
to employ the secular equations to study the statistics of the Solar System.
Yet, we will also compare their predictions to those of a non-averaged
comprehensive dynamical model in Sect. 7.

###  2.2 Frequency Analysis

We employed the frequency analysis (FA) technique proposed by Laskar (1988,
1993) to extract the fundamental secular frequencies from the integrated
solutions. The method find a quasi-periodic approximation
$f^{\prime}(t)=\sum_{k=1}^{N}a_{k}e^{i\nu_{k}t}$ of a function $f(t)$ over a
time span interval $[0,T]$. It first finds the strongest mode, which
corresponding to the maximum of the function:

$\phi(\sigma)=\langle f(t)|e^{i\sigma
t}\rangle=\frac{1}{T}\int^{T}_{0}\chi(t)f(t)e^{-i\sigma t}dt,$ (2)

where $\chi(t)$ is a weight function which improves the precision of the
maximum determination; it is chosen to be the Hanning window filter, that is
$\chi(t)=1+\cos(\pi t/T)$. The next step is the Gram-Schmidt
orthogonalization. The complex amplitude $a_{1}$ of the first frequency
$\nu_{1}$ is calculated by orthogonal projection of function $f(t)$ on
$e^{i\nu_{1}t}$. This mode is then subtracted from the function $f(t)$ to get
a new function $f_{1}(t)=f(t)-a_{1}e^{\nu t}$. The process is then repeated
with this newly-obtained function until desired N strongest modes are
obtained. This technique works very well for weakly-chaotic systems like the
Solar System when variables can be decomposed into quasi-periodic modes over a
sufficiently short period of time. It has been proven that this algorithm
converges towards the true frequencies much faster than the classical Fast
Fourier Transform (Laskar 2005). Therefore, it is a good tool to study the
chaotic diffusion of the fundamental frequencies. In this work, we used a
routine naftab written in the publicly-available computer algebra software
TRIP (Gastineau & Laskar 2011), developed at IMCCE, to apply directly the
frequency analysis.

To extract the fundamental secular frequencies of the Solar System, we apply
the FA to the proper variables $(z_{i}^{\bullet},\zeta_{i}^{\bullet})_{i=1,8}$
of the secular equations (Laskar 1990). Each fundamental frequency is obtained
as the frequency of the strongest Fourier component of the corresponding
proper variable. To track the time evolution of the frequencies, the FA is
applied with a moving interval whose sliding step is 1 Myr. The interval sizes
are 10 Myr, 20 Myr and 50 Myr.

##  3 Estimation of probability density functions

The samples of this study consist of the secular frequencies of the
astronomical solutions which are obtained by integrating the secular equations
(Eq. 1) from very close initial conditions. Due to this initial proximity, the
correlation of the solutions in the samples lasts for a long period of time,
but will eventually diminish. Our objective is to obtain a robust estimation
of the marginal probability density functions (PDFs) from these correlated
samples. In fact, this correlation is the main motivation for our use of the
estimation methods in this section. Details of our samples are described in
the first part of Sect. 4.

We use Kernel Density Estimation (KDE) to estimate the time-evolving marginal
PDFs of the fundamental frequencies of the Solar System. In addition, the
statistical uncertainty of our density estimations (i.e. PDF estimations) is
measured by moving block bootstrap (MBB). To our knowledge, this application
of MBB for KDE of a time-evolving sample, whose correlation changes overtime,
is new. Therefore, in order to ensure the validity of the method, we have
carried out several tests in Sect. 3.3.

###  3.1 Kernel Density Estimation

We chose KDE, also known as the Parzen–Rosenblatt window method, as our
preferred non-parametric estimator of the PDFs of the fundamental frequencies
of the Solar System because of its high convergence rate and its smoothing
property (Rosenblatt 1956; Parzen 1962).

We briefly present the method here. Let
$\mathbf{X}=\\{X_{1},X_{2},\dots,X_{n}\\}$ be a univariate independent and
identically distributed (i.i.d.) sample drawn from an unknown distribution P
with density function $p(x)$. Then the KDE of the sample is defined as:

$\hat{p}_{h}(x|\mathbf{X})=\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{x-X_{i}}{h}\right),$
(3)

where $K$ is a non-negative kernel function, and $h$ is bandwith of the KDE.
The choice of bandwidth $h$ is much more important than the choice of kernel
$K$, which is chosen to be the standard normal distribution in this paper. In
this work we shall consider bandwidths of the following form:

$BW_{\beta}=0.9\min\left(\hat{\sigma},\frac{IQR}{1.34}\right)n^{-\beta},$ (4)

where $\hat{\sigma}$ is the standard deviation of the sample, IQR is its
interquartile range and $\beta$ is a constant of choice. The bandwidth with
$\beta=1/5$ corresponds to the Silverman’s rule of thumb (Silverman 1986).
$BW_{\beta=1/5}$ is a slightly modified version of the optimal choice of
bandwidth for Gaussian distributed data for better adaption to non-Gaussian
data. Bias error and variance error of KDE with this bandwidth will be in the
same order of magnitude. Undersmoothing, i.e choosing a smaller bandwidth,
shrink the bias so that the total error is dominated by the variance error,
which can then be estimated by bootstrap method (Hall et al. 1995); the common
value of $\beta$ for undersmoothing is $1/3$.

When the sample is identically distributed but correlated, kernel density
estimation is still valid under some mixing conditions (Robinson 1983; Hart
1996). Indeed, in case of observations that are not too highly correlated, the
dependence among the data that fall in the support of the kernel function $K$
can be actually much weaker than it is among the entire sample. This principle
is known as _whitening by windowing_ (Hart 1996). Therefore, the correlation
in the sample does not invalidate the employ of kernel density estimation, and
only impacts the variability of the estimation (see Sect. 3.2). With regard to
the choice of bandwidth, Hall et al. (1995) suggested that using the
asymptotically optimal bandwidth for the independent data is still a good
choice, even for some strongly dependent data sequences. The samples generated
by a chaotic measure-preserving dynamical system (as in the case of the
numerical integration of the Solar System) resembles those of mixing
stochastic processes, hence the theory of KDE should also be applicable for
dynamical systems, although the formulation might be different (Bosq & Guégan
1995; Maume-Deschamps 2006; Hang et al. 2018).

###  3.2 Moving Block Bootstrap

Since the seminal paper by Efron (1979), bootstrap has become a standard
resampling technique to evaluate the uncertainty of a statistical estimation.
Bias and variance errors of KDE with the choice of bandwidth $BW_{1/5}$ (Eq.
4) are in the same order of magnitude, hence one should either undersmooth,
i.e. choose a smaller bandwidth, to minimize the bias (Hall 1995) or use an
explicit bias correction with appropriate studentizing (Cheng et al. 2019;
Calonico et al. 2018). However, applying naively the i.i.d. bootstrap
procedure on dependent data could underestimate the true uncertainty, because
all the dependence structure would be lost just by “scrambling” data together
(Hart 1996). To remedy the problem, Kunsch (1989) and Liu et al. (1992)
independently devised the moving block bootstrap (MBB), which later becomes
standard practice to evaluate the uncertainty in dependent data (See Kreiss &
Lahiri (2012) for a review). Although the MBB for smooth functional has been
intensively studied, the literature on MBB for KDE of dependent data is very
limited. Recently, Kuffner et al. (2019) has formulated an optimality theory
of the block bootstrap for KDE under weak dependence assumption. They proposed
both undersmoothing and explicit bias correction scheme to obtain the sampling
distribution of the KDE. However, good tuning parameters, which are generally
difficult to find if the data are from an unknown distribution, is required to
provide a decent result. In this paper, we propose to overcome this problem
with an inductive approach: the optimal tuning parameters obtained in a known
model are tested on different models, and then extrapolated to our study - the
Solar System.

#### $\mathbf{Procedure\ of\ MBB}$

We briefly describe the undersmooth MBB for the kernel density estimation
method (the detailed version can be found in (Kuffner et al. 2019)). Suppose
$\mathbf{X}=\\{X_{1},X_{2},\dots,X_{n}\\}$ is a dependent sample of a mixing
process with underlying density function $p(x)$. The KDE of the sample is
$\hat{p}_{h}=\hat{p}_{h}(x|\mathbf{X})$; the hat overhead a given quantity
denotes its estimated value. We use MBB to estimate the distribution of
$\delta(x)=\hat{p}_{h}(x)-p(x)$, where $x\in\boldsymbol{\Omega}$, and
$\boldsymbol{\Omega}$ is the domain of interest. Let $l$ be an integer
satisfying $1\leq l\leq n$. Then $B_{i,l}=\\{X_{i},X_{i+1},\dots,X_{i+l-1}\\}$
with $i\in\\{1,\dots,n-l+1\\}$ denotes all the possible overlapping block of
size $l$ of the sample. Suppose $l$ divides $n$ for simplicity, then $b=n/l$.
MBB samples are obtained by selecting $b$ blocks randomly with replacement
from $\\{B_{1,l},\dots,B_{n-l+1,l}\\}$. Serial concatenation of $b$ blocks
will give $n$ bootstrap observations:
$\mathbf{X}^{*}_{l}=\\{B^{*}_{1,l},\dots,B^{*}_{b,l}\\}$. By choosing $l$
large enough (preferably larger than the correlation length), MBB sample can
retain the structure of the sample dependence. For $k>0$, the KDE of the
bootstrap sample is $\hat{p}^{*}_{k,l}=\hat{p}_{k}(x|\mathbf{X}^{*}_{l})$ and
its expectation
$\mathbb{E}[\hat{p}^{*}_{k,l}]=\hat{p}_{k}(x|B_{1,l},\dots,B_{n-l+1,l})$.
Let’s define

$\delta^{*}_{k,l}(x)=\sqrt{\frac{k}{h}}(\hat{p}^{*}_{k,l}-\mathbb{E}[\hat{p}^{*}_{k,l}]),$
(5)

so that if $h$ is chosen properly to reduce the bias to be asymptotically
negligible with respect to the stochastic variation, then the MBB distribution
$P(\delta^{*}_{k,l}(x)|\mathbf{X})$ is a consistent estimator of the error
distribution $P(\delta(x))$ when $h\rightarrow 0$, $nh\rightarrow\infty$,
$k\rightarrow 0$, $lk\rightarrow\infty$ and $n/l\rightarrow\infty$. Note that
if $l=1$, then $k=h$ and MBB reverts to the undersmoothing procedure for i.i.d
sample studied by Hall et al. (1995). The efficiency of this estimator depends
sensitively on two tuning parameters $l$ and $k$. We are interested in the
uncertainty of the KDE, which is characterised by the confidence interval
$CI_{1-\alpha}(x)$ and the confidence band $CB_{1-\alpha}$, which are defined
as:

$P(|\delta(x)|<CI_{1-\alpha}(x))=1-\alpha,\\\ $ (6)
$P(|\delta(x)|<CB_{1-\alpha}\forall x\in\boldsymbol{\Omega})=1-\alpha,$ (7)

where $\alpha$ denote the level of uncertainty; for example, $\alpha=0.05$
denotes $95\%$ CI.

They can all be estimated by the MBB distribution
$P(\delta^{*}_{k,l}(x)|\mathbf{X})$ as:

$\displaystyle P(|\delta^{*}_{k,l}(x)|$
$\displaystyle<\widehat{CI}_{1-\alpha}(x))=1-\alpha,$ (8) $\displaystyle
P(|\delta^{*}_{k,l}(x)|$ $\displaystyle<\widehat{CB}_{1-\alpha}\forall
x\in\boldsymbol{\Omega})=1-\alpha,$ (9)

In practice, we will also use CB and CI without the hat overhead to denote the
estimated value for short.

Our choice of the parameters of the MBB procedure $l$ and $k$ is based on the
effective sample size $n_{\text{eff}}$, defined as (Kass et al. 1998):

$n_{\text{eff}}=\frac{n}{1+2\sum_{k=1}^{\infty}\rho(k)},$ (10)

where $\rho(k)$ is the sample autocorrelation of lag $k$. The block length $l$
is chosen by the sample correlation size $l_{\text{corr}}$ as

$l=l_{\text{corr}}:=\frac{n}{n_{\text{eff}}},$ (11)

and the bootstrap bandwidth is parametrized as

$k=h(c_{0}+(1-c_{0})l_{\text{corr}}^{-\gamma}),$ (12)

where $\gamma$ and $c_{0}$ are two optimizing constants. The reason for this
choice of parametrization is two folds. First, when
$l_{\text{corr}}\rightarrow 1$, the sample become independent, then
$k\rightarrow h$. Secondly, the rate of change of $k$ with respect to $l$
should be greater when $l_{\text{corr}}$ is small than when it is large.
Therefore when $l_{\text{corr}}\gg 1$, the optimal value of $k$ should be
quite stable. We also observe experimentally that the optimal value of $k$ is
indeed relatively stable at around $2h$ as long as $l=l_{\text{corr}}\gg 1$.
So we simply choose $\gamma=1$ and $c_{0}=2$. This choice of parameters turns
out to be quite robust as it will be demonstrated by the two numerical
experiments in Sect. 3.3.

The literature on KDE and MBB focuses on stationary, weakly dependent
sequences. The data in our case, however, are different: they are not strictly
stationary; the formulation of mixing condition might be different (Hang et
al. 2018); the correlation in the sample is not constant but evolving with
time; finally and most importantly, the data structure is different. Our
sample units, which are the orbital solutions, are ordered based on their
initial distances in phase space. The solutions evolve over time but their
order remains unchanged. Therefore statistical notions such as correlation and
stationarity should be considered within this framework, for the sample at
fixed values of time. Because of the differences presented, an optimality
theory, which is not currently available, might be needed for this case.
However, we assume it would not differ significantly from the orthodox
analysis, and that a decent working MBB procedure could be obtained with some
good choice of parameters. This will be tested in the section below.

###  3.3 Numerical experiments

We perform the KDE/MBB procedure on two numerical experiments to ensure its
validity. The double Gaussian mixture model (DGMM) with different degrees of
correlation is used to calibrate the algorithm of KDE/MBB method, that is to
calibrate and test the tuning parameters (Eqs. 11-12). Then the resulting
algorithm is applied on the Fermi Map as an example of a real dynamical system
for assessment.

Figure 1:  Kernel density estimation (red line) with bandwidth $h=BW_{1/3}$
and its 90% pointwise confidence interval (red band) from a MCMC sample of
$n=10^{4}$ units with $l_{\text{corr}}=140$ of double Gaussian mixture model
(Eq. 13 \- black line).

#### $\mathbf{Double\ Gaussian\ Mixture\ Model}$

The calibration of KDE/MBB was done on Markov chain Monte Carlo (MCMC)
sequences that sample a double Gaussian distribution. Our particular choice of
the DGMM, inspired by Cheng et al. (2019), is

$f_{\text{DG}}(x)=0.6\phi(x)+0.4\phi(x-4),$ (13)

where $\phi(x)$ is the standard normal distribution. The MCMC sequence
$X_{1},\dots,X_{n}$ is obtained by the Metropolis-Hasting algorithm, with a
Gaussian proposal distribution whose standard deviation $\sigma_{p}$
characterises the correlation length of the sequence (Metropolis et al. 1953;
Hastings 1970). We can then either vary $\sigma_{p}$ or perform thinning to
obtain a sequence of a desired correlation length. The size of each MCMC
sequence is $n=10^{4}$. The initial state of the sequence is directly sampled
from the distribution $f_{\text{DG}}$ itself, so that burn-in is not necessary
and the whole sequence is usable.

On each MCMC sequence, we apply the KDE/MBB procedure with 500 MBB samples and
the parameters specified above (Eqs. 11-12) with the undersmoothing bandwidth
$h=BW_{1/3}$, to get the uncertainty estimation (standard error, confidence
interval and confidence band). Figure 1 shows an example of the density
estimation of a MCMC sequence of $10^{4}$ units, with correlation length
$l_{\text{corr}}=140$; the sequence is generated by choosing
$\sigma_{p}=0.25$. It is clear in the figure that the true PDF (in black) lies
inside the range of the 90$\%$ CI estimated by MBB.

To assess the precision of the MBB uncertainty estimation, at each value of
$l_{\text{corr}}$, we sample 1000 different MCMC sequences, and apply the
KDE/MBB process on each of them to obtain their estimation of the confidence
interval. We thus obtain a distribution of confidence interval estimated by
MBB, which are depicted by the bands of mean $\pm$ standard deviation (Fig.
2). With the knowledge of the true PDF (Eq. 13), the true values of CI and CB
are obtained from this collection of KDEs at each correlation length (Eqs.
6-7). Figure 2 compares the true values of the pointwise 90% CI with the
distribution of its estimated values by MBB at various correlation lengths. We
see that the $90\%$ confidence interval of KDE estimated by MBB is accurate
across a very wide range of $l_{\text{corr}}$, that is the true values almost
always lie inside the estimation bands.

Figure 2: Estimation of the confidence interval (90% CI) of KDEs from
different samples with different correlation lengths. Solid lines are the true
values and the bands denote the distribution of MBB estimations (mean $\pm$
standard deviation), both are computed from $10^{3}$ different samples, each
composes of $n=10^{4}$ units generated from the double Gaussian mixture model
by MCMC algorithm with the same $l_{\text{corr}}$.

#### $\mathbf{Fermi\ Map}$

The second test for the MBB scheme is a 2-dimensional chaotic Hamiltonian
system: The Fermi map (Fermi 1949; Murray et al. 1985). This simple toy model,
originally designed to model the movement of cosmic particles, describes the
motion of a ball bouncing between two walls, one of which is fixed and the
other oscillating sinusoidally. The equations of the Fermi map reads:

$\displaystyle u_{t+1}$ $\displaystyle=u_{t}+\sin\psi,$ (14)
$\displaystyle\psi_{t+1}$ $\displaystyle=\psi_{t}+\frac{2\pi
M}{u_{t+1}}\pmod{2\pi},$

Figure 3: Density estimation of Fermi Map. The black line (“true” density) is
the KDE with bandwidth $h=BW_{1/5}$ from a sample of the velocity $u$ of
$n=10^{6}$ solutions of Fermi Map at t = 30. Red line denotes the KDE with
bandwidth $h=BW_{1/3}$ with its 95% confidence band estimated by MBB (red
band) from a sample of $n=10^{3}$ solutions at the same time. Blue-dashed line
represents the KDE with bandwidth $h=BW_{1/5}$ from the same sample.

where $u_{t}$ and $\psi_{t}$ are the normalized velocity of the ball and the
phase of the moving wall right before the $t^{\text{th}}$ collision,
respectively; the stochastic parameter of the system $M$ is chosen here to be
$10^{4}$. The system is studied in the region of large-scale chaos. Our sample
is obtained by evolving Eq. 14 from $n=10^{3}$ initial conditions. The initial
conditions $u_{0},\psi_{0}$ are drawn from uniform distributions on
$[90-10^{-5},90+10^{-5}]$ and $[0,2\pi]$, respectively; then they sorted by
ascending value of $\psi_{0}$. The sorting step is imperative to quantify the
sample dependence as the initially-ordered neighboring solutions still tend to
be correlated afterward and the autocorrelation function can be calculated
from ordered samples. Large initial phase variation will guarantee that
chaotic diffusion will be immediately perceptible and the PDF of the velocity
will center around its initial distribution. At each time $t$, we compute the
KDE of the sample with the bandwidth $BW_{1/3}$. The KDE uncertainty is
estimated by applying the MBB 400 times with the same parameters specified
above (Eqs. 11-12). The whole process is then repeated 300 times with
different sets of initial conditions. In the Fermi Map experiment, the true
analytical form of the density is not available, a numerical “true” density,
which is obtained by calculating the KDE with bandwidth $BW_{1/5}$ from a
$n=10^{6}$ sample, is used instead to assess the validity of the MBB
uncertainty estimation. From this “true” density and the KDEs from 300
samples, we are able to determine the true confidence band (Eq. 6-7).

Figure 3 shows an example of KDE and its confidence band estimated by MBB at
$t=30$. Although the estimated confidence band is valid, as the true curve
lies completely in the band, the KDE itself looks quite jagged. The jagged KDE
is the result of the undersmooth bandwidth $h=BW_{1/3}$ so the bias is
dominated by the variance error. In fact, with the rule-of-thumb bandwidth
$h=BW_{1/5}$, we can have a smoother KDE (cf. blue-dashed line of Fig. 3).
This choice is valid because its uncertainty is always smaller than that
estimated with $h=BW_{1/3}$. Therefore we choose to use the uncertainty
estimation computed with $h=BW_{1/3}$ as an upper bound of the uncertainty of
a KDE with the rule of thumb bandwidth. For example, a $95\%$ CB of a KDE with
$h=BW_{1/3}$ will represent an upper bound of a $95\%$ CB of a KDE with
$h=BW_{1/5}$ (Fig. 4). The same also applies for the confidence interval.

Also from Fig. 4, we can see that the estimation of the confidence band
follows well the true value. This is quite remarkable because, first we can
extend the use of MBB to measure the uncertainty of density estimation from
solutions of chaotic dynamical system, where the correlation of the sample
defined by initial distances in phase space changes over time, and secondly,
our simple choice of MBB parameters appears to work well across very different
models.

Figure 4: Comparison of the confidence band estimation with its true value.
The blue region represents the mean $\pm$ 3 standard deviations of the
confidence bands estimated by MBB of KDEs with bandwidth $h=BW_{1/3}$ from 300
different samples, each consisting in $n=10^{3}$ solutions of the Fermi Map;
the true values of confidence band are also calculated from the KDEs with
bandwidth $h=BW_{1/5}$ (red dotted line) and KDEs with bandwidth $h=BW_{1/3}$
(blue dotted line) from the same 300 different samples.

Figure 5: Density estimation with 95% pointwise confidence interval of the
fundamental frequencies of the four sets of the batch $\\{Z_{i}\\}_{i=1,4}$.
Frequency values are obtained by FA over an interval of 20 Myr centering at
150 Myr in the future.

###  3.4 Combining samples

After having obtained a well-tested uncertainty estimator of the KDE for
correlated data, we proceed with a practical question: given several samples
of different correlation length $l_{\text{corr}}$, how can we efficiently
combine them? Let’s assume we have $m$ samples of the same size n:
$\mathbf{X}_{1}=\\{X^{1}_{1},X^{1}_{2},\dots,X^{1}_{n}\\},\dots,\mathbf{X}_{m}=\\{X^{m}_{1},X^{m}_{2},\dots,X^{m}_{n}\\}$,
where the correlation within each sample is different. The KDE of these $m$
samples are $\hat{p}_{1},\dots,\hat{p}_{m}$, and their pointwise standard
error can be estimated by MBB as $\hat{\sigma}_{1},\dots,\hat{\sigma}_{m}$.

If all the samples conform to the same probability density $p$, we can simply
use the inverse variance weighting to get our combined KDE:

$\hat{p}_{\text{wa}}=\frac{\sum_{i}\hat{\sigma}_{i}^{-2}\hat{p}_{i}}{\sum_{i}\hat{\sigma}_{i}^{-2}},$
(15)

where $\hat{p}_{i}$ and $\hat{\sigma}_{i}^{2}$ are the individual KDE and its
pointwise estimated variance, respectively. The variance of the weigted
average will be
$\hat{\sigma}_{\text{wa}}^{2}=(\sum_{i}\hat{\sigma}_{i}^{-2})^{-1}$. With this
choice, the variance of the weighted average $\hat{p}_{\text{wa}}$ will be
minimized and we can thereby have the most accurate estimation of the density
$p$ (Hartung et al. 2011).

This inverse variance weighting is only applicable if we assume that all our
samples follow the same underlying distribution, as for example in the double
Gaussian mixture model. In the case of the Solar System, different samples
come from different sets of initial conditions. The densities evolving from
different initial conditions will be different. Although they might converge
towards a common density function after a large period of time, the assumption
that different samples follow the same distribution is not generally true.
However, if the differences between the density estimations from different
samples are small compared to their in-sample density uncertainty, then we can
assume a common true density that all the sample share, and the inverse
variance weighted mean and its minimized variance will in this case capture
this assumed density function.

A more conservative approach consists in not trying to estimate the common
true density, but rather capture the variability stemming from different
samples and combine it with their in-sample uncertainty. In this case, a
combined KDE may be introduced as a pointwise random variable, whose
distribution function is given by:

$\hat{p}_{\text{c}}\sim
m^{-1}\sum_{i=1}^{m}P(\hat{p}_{i}(x)-\delta^{*}_{i}(x)|\mathbf{X}_{i}),$ (16)

where $\delta^{*}_{i}(x)$ is defined in Eq. (5) and
$P(\hat{p}_{i}(x)-\delta^{*}_{i}(x)|\mathbf{X}_{i})$ is the MBB estimation of
the distribution of the pointwise KDE from sample $\mathbf{X}_{i}$. If the m
samples are good representatives of any other sample taken from a certain
population, then $\hat{p}_{\text{c}}$ is a reasonable choice for the KDE from
an arbitrary sample generated from that same population. We sample
$\hat{p}_{\text{c}}$ in a pointwise manner. For each value of $x$, the same
number of realizations are drawn from each of the m MBB distributions, which
are assumed here to be Gaussian for practical purpose. It should be noticed
that a realization of $\hat{p}_{\text{c}}(x)$ is not a continuous curve. When
needed, we can choose the pointwise median of $\hat{p}_{\text{c}}$ as the
nominal continuous combined KDE.

##  4 Application to the Solar System

The goal of this paper is to have a consistent statistical description of the
propagation of the dynamical uncertainty on the fundamental secular
frequencies of the Solar System induced by its chaotic behavior; that is,
simply put, to obtain their time-evolving marginal PDF. We first sample the
initial orbital elements of the Solar System planets closely around the
reference value, then numerically integrate the secular equations (Eqs. 1)
from those initial conditions to obtain a sample of orbital solutions. KDE is
then used to estimate the marginal PDF of the frequencies of the sample at
fixed value of time and then the MBB method is applied to estimate the
uncertainty of the density estimation.

The evolution of our sample can be divided into two stages: Lyapunov
divergence and chaotic diffusion. In the first stage, because the initial
density is extremely localized around the reference value, all solutions
essentially follow the reference trajectory, the difference between the
solutions is very small but diverges exponentially with a characteristic
Lyapunov exponent of $\sim 1/5\ \text{Myr}^{-1}$ (Laskar 1989). The solutions
in the first stage are almost indistinguishable and the correlation between
them is so great that regardless of how many solutions in the sample we
integrated, the effective size of the sample is close to one. The second stage
begins when the differences between the solutions are large enough to become
macroscopically visible. The Lyapunov divergence saturates and gives place to
chaotic diffusion. The correlation between the solutions starts to decrease,
the distribution of the sample settles and the memory of the initial
conditions fades. It can take several hundreds of millions of years for the
sample to forget its initial configuration. Contrary to the exponential growth
in the first stage, the dispersion of the samples expands slowly with a power
law in time (see Fig. 10). The time boundary between the two stages depends on
the dispersion of the initial conditions, the wider they are, the faster the
second stages come and vice versa. If they are chosen to represent the
uncertainty of the current ephemeris, the second stage should take place
around 60 Myr in the complete model of the Solar System (Laskar et al. 2011a,
b).

In this section, we focus on the statistical description of the fundamental
frequencies of the Solar System in the second stage. The aim is to have a
valid estimation of time-evolving PDF of the frequencies beyond 60 Myrs.
However, the PDF evolution generally depends on the choice of initial
conditions. Moreover, the simplification of the secular equations compared to
the complete model of the Solar System could, a priori, provide results which
are not sufficiently accurate.

Figure 6: Density estimation with 95% pointwise confidence interval of the
fundamental frequencies of the 12 sets coming from the 3 batches $\\{X_{i}\\}$
(yellow colors), $\\{Y_{i}\\}$ (blue colors) an $\\{Z_{i}\\}$ (green colors).
Estimations of sets from the same batch have the same color. Frequencies are
obtained by FA over an interval of 20 Myr centering at 100 Myr in the past
(first and third row) and 200 Myr in the past (second and fourth row). The
time reference is the time of $\\{Z_{i}\\}$, while the solutions of
$\\{X_{i}\\}$, $\\{Y_{i}\\}$ are shifted ahead 30 Myr and 20 Myr respectively.
Figure 7: Density estimation with 95% pointwise confidence interval of the
fundamental frequencies of the 12 sets coming from the 3 batches $\\{X_{i}\\}$
(yellow colors), $\\{Y_{i}\\}$ (blue colors) an $\\{Z_{i}\\}$ (green colors)
and from the 2 500 solutions of complete model of Laskar & Gastineau (2009)
(red colors), which is denoted as La09. Estimations of sets from the same
batch have the same color. Frequencies are obtained by FA over an interval of
20 Myr centering at 150 Myr (first and third rows) and 350 Myr (second and
fourth rows) in the future. The time reference is the time of $\\{Z_{i}\\}$,
while the solutions of $\\{X_{i}\\}$, $\\{Y_{i}\\}$, La09 are shifted 30 Myr,
20 Myr and -70 Myr respectively.

###  4.1 Choice of initial conditions

Solutions | Offsets | $\epsilon$
---|---|---
$X_{i}$ |  $-5000\,\epsilon$ to $5000\,\epsilon$ | $10^{-10}$
$Y_{i}$ |  $-5000\,\epsilon$ to $5000\,\epsilon$ | $10^{-8}$
$Z_{i}$ |  $-5000\,\epsilon$ to $5000\,\epsilon$ | $10^{-11}$
Table 1: Offsets of the initial eccentricities of the 4 planets: $\\{$Mercury,
Venus, Earth, Mars$\\}$, which corresponds to $i=\\{1,2,3,4\\}$. Different
integrations correspond to offsets of $N\epsilon$ in eccentricity of a single
planet for $N=-5\times 10^{3},\dots,+5\times 10^{3}$, while other variables
are kept to their nominal values.

For a complete model of the Solar System, the initial conditions should be
sampled in such a way that they are representative of the current planetary
ephemeris uncertainty. Nevertheless, for a simplified secular model (Eq. 1),
the difference with the complete model is greater than the ephemeris
uncertainty. An optimized secular solution follows the complete solution
initially but will depart from it long before 60 Myr (Laskar et al. 2004). A
direct adaptation of the current planetary ephemeris uncertainty into the
initial conditions of the secular model can thus be misleading. Therefore, we
adopt a more cautious approach, that is to study first the effect of sampling
the initial conditions on the PDF estimation.

Initial conditions can be sampled in many ways, especially in a high-
dimensional system like the Solar System. Our choice of initial conditions is
quite particular, but they encompass different possible ways of sampling
initial conditions (Table 1). There are 3 batches: $\\{X_{i}\\}$,
$\\{Y_{i}\\}$ and $\\{Z_{i}\\}$, which corresponds to 3 different variation
sizes $\epsilon$ of initial conditions. Each batch composes of 4 different
sets of samples. Each set contains 10 000 initial conditions, where
eccentricity of the associated planet is linearly-spaced from the reference
value, with the spacing $\epsilon$ corresponding to the batch it belongs to.
We then integrated the secular equations (Eqs. 1) from these initial
conditions 500 Myr to the past and 500 Myr to the future. For the batch
$\\{Z_{i}\\}$, the integration time is 5 billion years in both directions. The
frequencies are then extracted using frequency analysis (section 2.2). It
should be noted that we do not aim to obtain the joint probability
distribution of all the fundamental frequencies, but rather their individual
marginal PDFs, i.e. the PDF of one frequency at a time. The marginal PDFs of
the frequencies are estimated by the KDE with rule-of-thumb bandwidth
($BW_{1/5}$), upper bounds of their 95% confidence intervals (denoted for
short as 95% CI) are measured by the MBB method with the bandwidth
$h=BW_{1/3}$ and the optimized parameters (Eqs. 11-12) from 1 000 MBB samples.

We first compare the evolution of the density of the four sets in each batch
in section 4.2. In section 4.3, the second test is performed to compare the
statistics between the batches. Robustness of secular statistics will be
assessed by these two tests, which will additionally shed light on the
initial-condition-dependence aspect of the statistics. All the density
estimations from these sets will then be compared with that of 2 500 complete
solutions obtained in the previous work of (Laskar & Gastineau 2009), to test
the accuracy of the secular statistics. It should be reminded that this
numerical experiment needed 8 millions hours of CPU time, the output of which
was saved and could thus be used in the present study.

When comparing two sets of different sizes, because the rates of divergence in
the first stage are similar, the wider set reaches the chaotic diffusion phase
faster than the more compact set, hence it is essentially diffusing ahead for
a certain time in the second stage. Therefore, in order to have a relevant
comparison, a proper time shift is introduced to compensate for this effect.
We shift $\\{X_{i}\\}$ and $\\{Y_{i}\\}$ ahead 30 Myr and 20 Myr respectively,
while keeping the time of $\\{Z_{i}\\}$ as reference. This choice is motivated
by the fact that the transition to the chaotic diffusion of $\\{Z_{i}\\}$ is
around 50 - 60 Myr, which is indicated by direct integration of the Solar
System (Laskar et al. 2011a, b).

###  4.2 First test: different samples of the same variation size $\epsilon$

Comparing the density estimation of different sets in the same batch is the
first test of prediction robustness from our secular model. The evolution of
the density estimations will be sensitive to how initial conditions are
sampled and therefore this initial-condition sensitivity must be quantified
for a valid prediction. In the first test, we compare the time-evolving PDFs
whose initial conditions are sampled with the same variation size $\epsilon$
but in different variables.

The result is quite clear. The different sets of the same batch slowly lose
the memory of their initial differences due to chaos and then converge towards
the same distribution. This convergence is illustrated by Fig. 5, which shows
that the density estimation of $(g)_{i=1,4}$ of the 4 sets of the batch
$\\{Z_{i}\\}$ nearly overlap each other at 150 Myr in the future. The rates of
convergence of different batches are different. So that although $\\{X_{i}\\}$
and $\\{Y_{i}\\}$ exhibit the same behavior as $\\{Z_{i}\\}$, they converge
differently with disparate rates: At around 100-150 Myr in the future, the
density estimations of the frequencies of $\\{X_{i}\\}$ nearly overlap each
other, this occurs at 150-200 Myr for $\\{Y_{i}\\}$, depending on the
frequency. Interestingly, for the samples that are integrated in the past, the
rates of convergence are higher, and the overlap generally happens at around
$-100$ Myr (see Fig. 6).

###  4.3 Second test: different samples of different variation sizes

Comparing the density estimation of the 3 batches $\\{X_{i}\\}$, $\\{Y_{i}\\}$
and $\\{Z_{i}\\}$ is our second test of robustness. Although the initial
conditions of the 3 batches are varied around the same reference values, the
ways they are sampled are different from each other, since the variation sizes
are different. Differences in the initial variation sizes mean that the
batches will enter the diffusion stage at different times and also at
different points in the phase space, so that the convergence between batches,
if exists, would take longer. The result of our test is summarised by the
density estimation of the frequencies at two times in the past, $-100$ Myr and
$-200$ Myr (Fig. 6). At $-100$ Myr, the density estimations from different
sets of each batch cluster around each other as described in the previous
test. Each batch forms a cluster of density estimations and the differences
between 3 clusters are noticeable. Moreover, the estimation uncertainty,
depicted by the colored band, is quite large. Fast forward 100 Myrs of chaotic
mixing, at $-200$ Myr, density estimations of the frequencies spread out,
estimation uncertainty shrinks, and more importantly, differences between the
3 batches are much smaller, and they will continue to diminish even further as
time goes. In the opposite time direction, the same phenomenon is observed but
the rate of convergence between the batches is slower (Fig. 7). At 150 Myr in
the future, the density estimation of frequencies of $\\{X_{i}\\}$ and
$\\{Z_{i}\\}$ has practically converged when those of $\\{Y_{i}\\}$ is still
trying to, yet differences between 3 batches are noticeable. However, the
estimations of all 12 sets from 3 batches nearly overlap each other at 350
Myr, which demonstrates that the effect of different initial samplings
vanishes by chaotic mixing. It should be noted that when looking at some
specific properties of the PDF such as means and variances, the differences
between the sets are small. For example, the differences between the means of
the PDF estimation of the 12 sets are generally smaller than $0.1\
^{\prime\prime}/\text{yr}$ for most of the fundamental frequencies at $-100$
Myr; at $-200$ Myr, the differences diminish to be twice as small.

###  4.4 Final test: comparison with the complete model

Figure 8: Density estimation with $95\%$ pointwise confidence interval of the
fundamental frequencies of the set $Z_{1}$, which are obtained by FA over an
interval of 10 Myr (blue colors), 20 Myr (yellow colors) and 50 Myr (green
colors) centering at $-150$ Myr (top row) and $-350$ Myr (bottom row).  Figure
9: Density estimation with 95% pointwise confidence interval of the
fundamental frequencies of the solar System from the $\\{Z_{i}\\}$ solutions
in the past at 500 Myr (red band), 1 Gyr (green band) and 4 Gyr (blue band);
the dashed curves with the corresponding colors denote their fitting
distribution (Eqs. 19 and Table 2). The frequencies are obtained by FA over an
interval of 20 Myr. Figure 10: Evolution of the parameters of the skew
Gaussian distribution in the fitting model (Eqs. 19) for density estimation of
the fundamental frequencies of the Solar System of $\\{Z_{i}\\}$ solutions in
the past (blue line) and in the future (green line). The bands denote $\pm 2$
standard error. The dashed orange lines denote the power-law fits for
$\sigma_{0}^{2}$ and the linear fits for $\mu_{0}$ and $\alpha_{0}$.

The convergence of the density estimations from different sets of initial
conditions (seen in the first two tests) does not guarantee convergence
between the secular model and the complete model of the Solar System. As a
result, whether the secular statistics will resemble that of the complete
model is still not clear. In this final test, we respond to this question by
comparing the density estimation from our secular solutions with those
obtained from the 2 500 solutions of the complete model integrated in the
future (Laskar & Gastineau 2009). The initial conditions of the complete
solutions are sampled in a way similar to the present work, that is one
variable (the semi-major axis of Mercury) is linearly-spaced with a spacing of
0.38 mm and a range of about one meter. We shift the complete solution
backward 70 Myrs to adjust for the difference in the initial variation sizes
and also for the difference between the initial divergence rates.

The density estimations of the frequencies of the 2 500 complete solutions are
shown along with the estimations from the secular model in Fig. 7. The
estimations of the complete model are either close to or overlap those of the
secular model even at 150 Myr. Both models even predict the same minor
features in the frequency density, for example, the second peaks of $g_{3}$ or
the tails of $g_{4}$. The origins of these features are related to the
resonances associating with $g_{3}$ and $g_{4}$, so having these same features
indicates that the secular model could capture the resonance dynamics of the
complete model. At 350 Myr, for most of the frequencies, the differences
between the results of the two models are very small, especially for some
frequencies like $g_{1}$ where it is difficult to distinguish the two models.
If we look at some specific properties of the PDF like its mean for example,
the differences between the secular model and the complete model are at the
same order with the variability of results from the same secular model. These
differences of means of the PDFs are generally smaller than $0.05\
\arcsec/\text{yr}$ at 350 Myr. Although some minor differences between the two
models are still visible, especially in $s_{1}$ for example, these differences
diminish as time goes. This convergence between the two models strongly
suggests the compatibility of the secular system with the direct integration
of the realistic model of the Solar System used in (Laskar & Gastineau 2009).

###  4.5 A complementary test on frequency analysis

The fundamental frequencies of the Solar System are central in our work, and
the method to obtain them is thus essential. In this section, we briefly
examine the frequency analysis method: FA (Sect. 2.2). FA searches a quasi-
periodic approximation of the solution of the Solar System with constant
frequencies over a time window $\Delta t$. So unique frequencies at time $t$
are extracted from an oscillating sequence in the time interval $[t-\Delta
t/2,t+\Delta t/2]$. For a quasi-periodic solution, the longer we choose the
time interval $\Delta t$, the more accurate the extracted frequencies are.
Nevertheless, the fundamental frequencies of the Solar System are expected to
vary over a few Lyapunov times, that is 5 Myr. Therefore we have a trade-off
between the extraction accuracy and time variation of the frequencies when
choosing $\Delta t$: when $\Delta t$ is too large, the obtained frequency will
tend to be the average of its variation over the same period. We chose $\Delta
t=20\ \text{Myr}$ as the standard FA interval. In some circumstances that
requires detection of rapid changes in frequencies such as the resonance
transition, a smaller $\Delta t$ is more favorable (see section 6.3).

The extracted frequencies are sensitive to the choice of $\Delta t$, yet,
their density estimation is relatively robust. Figure 8 compares the density
estimation of the eccentricity frequencies $(g_{i})_{1,4}$, which are
extracted by FA with 3 different $\Delta t$ at two different times. The
differences are generally small but still notable for $g_{2}$ and $g_{4}$ at
150 Myr, and they will diminish as time goes.

| $\mu_{0}\ [\arcsec/\text{yr}]$ | $a\ [\arcsec/\text{yr}]^{2}$ | $b$ | $\alpha_{0}$ | ${}\mu_{1}\ [\arcsec/\text{yr}]{}$ | ${}\sigma_{1}^{2}\ [\arcsec/\text{yr}]^{2}{}$ | ${}A_{1}{}$
---|---|---|---|---|---|---|---
${}\rule{0.0pt}{10.00002pt}g_{1}{}$ | ${}5.759+0.006\,T{}$ | ${}3.37\cdot 10^{-2\phantom{0}}{}$ | ${}0.52{}$ | ${}-2.25-0.50\,T{}$ |  |  |
${}g_{2}{}$ | ${}7.448-0.004\,T{}$ | ${}4.17\cdot 10^{-4\phantom{0}}{}$ | ${}0.70{}$ | ${}1.38+0.21\,T{}$ |  |  |
${}g_{3}{}$ | ${}17.269+0.002\,T{}$ | ${}6.63\cdot 10^{-3\phantom{0}}{}$ | ${}0.43{}$ |  |  |  |
${}g_{4}{}$ | ${}17.896+0.005\,T{}$ | ${}6.88\cdot 10^{-3\phantom{0}}{}$ | ${}0.41{}$ |  | ${}17.6755{}$ | ${}0.0034{}$ | ${}0.110-0.012\,T{}$
${}s_{1}{}$ | ${}-5.652-0.032\,T{}$ | ${}2.68\cdot 10^{-2\phantom{0}}{}$ | ${}0.83{}$ | ${}1.12+0.16\,T{}$ |  |  |
${}s_{2}{}$ | ${}-6.709+0.030\,T{}$ | ${}1.20\cdot 10^{-1\phantom{0}}{}$ | ${}0.76{}$ | ${}-2.94-1.23\,T{}$ |  |  |
${}s_{3}{}$ | ${}-18.773+0.009\,T{}$ | ${}2.86\cdot 10^{-2\phantom{0}}{}$ | ${}0.56{}$ | ${}-3.40-0.08\,T{}$ | ${}-18.5256{}$ | ${}0.0028{}$ | ${}0.023{}$
${}s_{4}{}$ | ${}-17.707+0.013\,T{}$ | ${}1.19\cdot 10^{-2\phantom{0}}{}$ | ${}0.68{}$ | ${}-1.73-0.28\,T{}$ |  |  |
${}\rule{0.0pt}{10.00002pt}g_{5}{}$ | ${}4.257454-2.1\cdot 10^{-6}\ T{}$ | ${}4.63\cdot 10^{-10}{}$ | ${}0.88{}$ |  |  |  |
${}g_{6}{}$ | ${}2.824523-1.4\cdot 10^{-4}\ T{}$ | ${}1.40\cdot 10^{-6\phantom{0}}{}$ | ${}0.84{}$ |  |  |  |
${}g_{7}{}$ | ${}3.087957-1.2\cdot 10^{-6}\ T{}$ | ${}4.80\cdot 10^{-10}{}$ | ${}1.11{}$ |  |  |  |
${}g_{8}{}$ | ${}6.730237\phantom{\,\ -1.0\cdot 10^{-2}\ T}{}$ | ${}9.89\cdot 10^{-11}{}$ | ${}1.49{}$ |  |  |  |
${}s_{6}{}$ | ${}-2.634787+1.5\cdot 10^{-5}\ T{}$ | ${}1.21\cdot 10^{-8\phantom{0}}{}$ | ${}0.85{}$ |  |  |  |
${}s_{7}{}$ | ${}-2.992527\phantom{\,\ -1.0\cdot 10^{-2}\ T}{}$ | ${}8.31\cdot 10^{-10}{}$ | ${}1.39{}$ |  |  |  |
${}s_{8}{}$ | ${}-6.917366\phantom{\,\ -1.0\cdot 10^{-2}\ T}{}$ | ${}2.93\cdot 10^{-11}{}$ | ${}1.47{}$ |  |  |  |
Table 2: Linear and power law fits for the time evolution of the parameters
(Fig. 10) of the skew Gaussian mixture model (Eq. 19) for the fundamental
frequencies of the Solar System. Col. 1 contains the considered secular
frequencies. In col. 2, we show the linear fits of $\mu_{0}$ which represent
the center of the distribution. The power law fit of $\sigma^{2}_{0}$ has the
form $\sigma^{2}_{0}(t)=aT^{b}$ (Eq. 20), where $T=t/(1\,\text{Gyr})$, $a$ and
$b$ are given in col. 3 and col. 4 respectively. Linear fits of the skewness
parameter $\alpha_{0}$ are given in col. 5. The last three columns show linear
fits of the secondary mode of $g_{4}$ and $s_{3}$ (Eq. 19). The parameter
$\sigma_{0}^{2}$ is fitted from 200 Myr to 5 Gyr in the past, while all the
other ones are fitted from 500 Myr.

##  5 Parametric fitting

From Sect. 4, we see that the long-term PDFs of the secular frequencies
possess a distinct Gaussian-like shape that flattens as time passes. It is
interesting to approximate these densities by a simple parametric model, so
that its parameters can characterise the shape of the density and summarize
its evolution. The model with its fitting parameters can also be used as an
approximation of the numerical densities for later application. For this
purpose, we will use the density estimation of the secular frequencies of the
Solar System from the batch $\\{Z_{i}\\}$, which composes of 40 000 different
orbits over 5 Gyr. The inner fundamental frequencies (i.e. the frequencies of
the inner planets $(g_{i},s_{i})_{i=1,4}$) are obtained by FA over an interval
of 20 Myr. For the outer fundamental frequencies (i.e. the frequencies of the
outer planets $(g_{i})_{i=5,8}$, $(s_{i})_{i=6,8}$), the FA interval is 50
Myr. The frequency $s_{5}$ is zero due the conservation of the total angular
momentum of the Solar System.

###  5.1 Skew Gaussian mixture model

Laskar (2008) found that the 250 Myr-averaged marginal PDF of eccentricity and
inclination of the inner Solar System planets are described quite accurately
by the Rice distribution, which is essentially the distribution of the length
of a 2D vector when its individual components follow independent Gaussian
distributions. In in our case, the density estimations of the fundamental
frequencies of the Solar System resemble Gaussian distributions, but many of
them get skewed as time passes, this is especially true for the inner
frequencies. To account for this skewness, we propose the skew normal
distribution as the fitting distribution to the density estimation of the
frequencies:

$f_{\mu_{0},\sigma_{0},\alpha_{0}}(x)=\frac{2}{\sigma_{0}}\phi\left(\frac{x-\mu_{0}}{\sigma_{0}}\right)\Phi\left(\alpha\left(\frac{x-\mu_{0}}{\sigma_{0}}\right)\right),$
(17)

where $\alpha$ is the parameter characterising the skewness; $\phi(x)$ denotes
the standard normal probability density distribution with mean $\mu_{0}$ and
standard deviation $\sigma_{0}$, and $\Phi(x)$ is its cumulative distribution
function given by

$\Phi(x)=\int_{-\infty}^{x}\phi(t)\
dt=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)\right],$
(18)

where erf denotes the error function.

Some of the frequencies, interestingly, has several secondary modes in their
density estimation apart from their primary one. Most of the secondary modes,
if exists, are quite small compared to the primary one. They are also often
short-lived, most of them emerge at the beginning of the diffusion stage and
disappear quickly after. Therefore they are not included in this parametric
model, whose aim is to fit the long-term PDF of the fundamental frequencies.
However, some persist for a long time and have a small but non-negligible
amplitude. To account for these secondary modes, we simply add Gaussian
functions to the skew Gaussian distribution and adjust for their amplitudes so
that the fitting distribution is ultimately the skew Gaussian mixture model:

$f(x)=A_{0}f_{\mu_{0},\sigma_{0},\alpha_{0}}(x)+\sum_{i=1}^{m}A_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2}),$
(19)

where $\sum_{i=0}^{m}A_{i}=1$, and m is the number of secondary modes. The
secondary modes are much smaller than the primary mode:
$A_{0}\gg\sum_{i=1}^{m}A_{i}$. In our case here, $m=1$ for $g_{4}$ and
$s_{3}$, while the asymptotic secondary modes for the other frequencies can be
considered negligible.

For the outer fundamental frequencies, we do not observe significant skewness
or secondary modes in their density estimations. Therefore, the density
estimations of the outer fundamental frequencies can be approximated by a
simple Gaussian distribution.

The density estimations of the frequencies are shown at 3 well separated times
in Fig. 9. The diffusion is clearly visible as the density estimations get
more and more disperse over time. Moreover, the density estimations get more
skew as time goes. For $g_{3}$ and $g_{4}$, the skewness of the density
estimations is small, so we assume $\alpha_{0}=0$ (Eq. 17) for these
frequencies for simplicity. The fundamental frequencies of the outer planets
are very stable. Their variations are much smaller then those of the inner
frequencies. Taking as an example the most unstable outer frequencies $g_{6}$,
its standard deviation at $-4$ Gyr is only about $2\cdot 10^{-3}\
\arcsec/\text{yr}$, while that of the most stable inner frequency $g_{2}$ is
15 times larger. Therefore, when considering a combination involving an inner
fundamental frequency and an outer one, the latter can be effectively regarded
as a constant.

The result of our fitting model is also plotted in Fig. 9. It is remarkable
that the density estimations of the frequencies are well approximated by the
fitting curve over 3 different epochs. It should be noticed that the base of
our fitting model - the skew Gaussian distribution only has 3 parameters. Only
for some frequencies, $g_{4},s_{3}$ for examples, more parameters are needed.
Nevertheless, the additional parameters only account for the minor features
and 3 parameters are sufficient to represent the bulk of the density
estimations over a long time scale.

Figure 11: Density estimation with 95% pointwise confidence interval of the
period of $g_{2}-g_{5}$ (top row) and $s_{3}-s_{6}$ (bottom row) from the
combining solutions from $\\{X_{i}\\}$, $\\{Y_{i}\\}$ an $\\{Z_{i}\\}$.
Frequencies are obtained by FA over an interval of 20 Myr centering at $-100$
Myr (left column), $-210$ Myr (middle column) and $-400$ Myr (right column).
Figure 12: Density estimation with 95% pointwise confidence interval (blue
colors) of the secular frequencies (right column) and their combinations (left
column) of the combined solutions from $\\{X_{i}\\}$, $\\{Y_{i}\\}$ and
$\\{Z_{i}\\}$, along with the corresponding values extracted from from the
La2010d solution (yellow lines) and Newark-Hartford data (red lines) (Olsen et
al. 2019) and the current values (green lines). The frequencies are obtained
by FA over an interval of 20 Myr centering at 210 Myr in the past. The time
reference is the time of $\\{Z_{i}\\}$. Figure 13: Density estimation with 95%
pointwise confidence interval (red colors) of $g_{4}-g_{3}$ (left column) and
$g_{4}$ (right column) at $-90$ Myr (first row) and $-300$ Myr (second row) of
$\\{Z_{i}\\}$, overlaid by their histograms (gray blocks) for better
visualisation. Frequencies values are obtained by frequency analysis over an
interval of 10 Myr.

Figure 14: Percentage of the solutions (black line) whose $g_{4}-g_{3}>0.9$ at
each time (left) and over a 10 Myr interval (right) and their 90% CI
estimation uncertainty (larger gray band) and their standard errors (smaller
gray band). The time interval of Libsack record, which is 80 Myr $-90$ Myr, is
denoted by the red band in the left plot and by a red line in the right plot.

###  5.2 Evolution of the parameters

The parameters of our fitting models are extracted by the method of least
squares, implemented by the routine curve_fit in the scipy package in Python.
To retrieve the statistical distribution (mean and standard deviation) of the
parameters of a given model, we implement a bootstrap approach based on Eq.
(16), with the assumption that pointwise standard errors of the KDE estimated
by MBB are independent. We remark that, with such an assumption, the variance
of the fitting parameters tends to be underestimated.

The time evolution of the mean of the parameters of the fitting models is
shown in Fig. 10, along with their $\pm 2$ standard error, for both past and
future. It turns out that the evolution of $\sigma_{0}^{2}$ is robustly fitted
from 200 Myr to 5 Gyr in the past by a power law function

$\sigma_{0}^{2}(t)=a\ T^{b},$ (20)

where $T=t/(1\,\text{Gyr})$, as shown in Fig. 10. For all the other
parameters, we perform a linear fit. All these fits are summarised in Table 2.

The differences between past and future evolutions are small and generally
tend to decrease with time. Therefore, the fit for the parameters in the past
given in Table 2 should also be representative of their future evolution. In
general, the parameters follow relatively smooth curves with distinct
tendencies. The skewness parameters $\alpha_{0}$, increasing in absolute
value, show that the PDFs of the inner fundamental frequencies get more and
more skew over time. The center of the distributions, indicated by $\mu_{0}$,
does not change significantly compared to $\alpha_{0}$ and $\sigma_{0}^{2}$
(Fig. 10). The secondary modes of $g_{4}$ and $s_{3}$ are also quite stable.

The diffusion of the frequencies is quantified by the increasing
$\sigma_{0}^{2}$, which is closely linked to their distribution variance. As
the exponents $b$ of the power laws in Table 2 are different from unity, the
chaotic diffusion of the fundamental frequencies turns out to be an anomalous
diffusion process. Interestingly, all the inner frequencies clearly undergo
subdiffusion, that is the exponent of the power law $b$ is smaller than 1
($b=1$ corresponding to Brownian diffusion). Therefore, an extrapolation of
the variance of the inner frequencies based on the assumption of linear
diffusion over a short time interval would generally lead to its
overestimation over larger times. On the contrary, the exponents $b$ are
either smaller or larger than unity for the outer frequencies. It should be
noted that, because the variations of the outer frequencies are very small,
the value of the corresponding exponents $b$ might be overestimated due to the
finite precision of frequency analysis.

##  6 Geological application

The aim of this project is to have a reliable statistical picture of the
secular frequencies of the Solar System beyond 60 Myr. It is interesting to
put recent geological results in this astronomical framework. First, this can
be used as a geological test of our study, and secondly the application
provides a glimpse of how astronomical data could be used in a
cyclostratigraphy study. We first show how the uncertainty of a widely-used
dating tool in astrochronology can be quantified, then we apply our work on
two recent geological results, which are from the Newark-Hartford data (Olsen
et al. 2019) and Libsack core (Ma et al. 2017).

###  6.1 Astronomical metronomes

Although it is not possible to recover the precise planetary orbital motion
beyond 60 Myr, some astronomical forcing components are stable and prominent
enough so that they can be used to calibrate geological records in the
Mesozoic Era or beyond (see (Laskar 2020) for a review). The most widely-used
is the 405-kyr eccentricity cycle $g_{2}-g_{5}$, which is the strongest
component of the evolution of Earth’s eccentricity. Recently the inclination
cycle $s_{3}-s_{6}$ has also been suggested for the time calibration of
stratigraphic sequences (Boulila et al. 2018; Charbonnier et al. 2018).
Although $s_{3}-s_{6}$ is not the strongest among the obliquity cycles, it is
quite isolated from other cycles and thus easy to be identified (Laskar 2020).
The main reason that $g_{2}-g_{5}$ or possibly $s_{3}-s_{6}$ can be used as
astronomical metronomes is their stability. Indeed, their uncertainty has to
be small to be reliably used for practical application.

In previous work, the uncertainty of the frequency combination was derived
from the analysis of a few solutions (Laskar et al. 2004, 2011a; Laskar 2020).
Here we use a much larger number of solutions and the KDE/MBB method to derive
both the PDF of the frequencies and their statistical errors. Starting from
the results of all the sets of orbital solutions (Sect. 4), we produce the
compound density estimation of the fundamental frequencies and their relevant
combinations following the conservative approach in Eq. (16). These data are
archived with the paper.

With this result, we can reliably estimate the uncertainty of the astronomical
metronomes. The uncertainty of the cycles $g_{2}-g_{5}$ and $s_{3}-s_{6}$ is
in fact the density width of the period of these frequency combinations, which
is shown in Fig. 11. As time goes, the two metronomes get more uncertain as
their density spreads due to the chaotic diffusion, but they are still
reliable. At 400 Myr in the past, the relative standard deviations (ratio of
standard deviation to mean) of $g_{2}-g_{5}$ and $s_{3}-s_{6}$ metronomes are
approximately $0.4\%$ and $1.26\%$, respectively.

###  6.2 Newark-Hartford data

In (Olsen et al. 2019), the astronomical frequencies were retrieved from a
long and well-reserved lacustrine deposit. The frequency signals of the
Neward-Hartford data are very similar to that of the astronomical solution
La2010d, which was taken among 13 available astronomical solutions (Laskar et
al. 2011a). Having 120 000 astronomical solutions at hand, we can derive a
more precise statistical analysis of this result. The fundamental frequencies
from the geological record are obtained in the following way. The data were
dated by the relatively stable 405-kyr cycle ($g_{2}-g_{5}$) and additionally
verified by a zircon U–Pb-based age model (Olsen et al. 2019). A frequency
analysis (Sect. 2.2) was performed on the geological data as well as on the
eccentricity of the astronomical solution La2010d to retrieve their strongest
frequencies. The frequency analysis of geological data are cross-checked with
that of astronomical solution to identify its orbital counterpart. For
example, the third strongest frequency $12.989\ ^{\prime\prime}/\text{yr}$
from Newark data is identified with the second strongest frequency $12.978\
^{\prime\prime}/\text{yr}$ from La2010d, which is $g_{3}-g_{5}$. The frequency
$g_{3}$ is then derived by summing the obtained combination with $g_{5}$,
which is nearly constant. The same goes for $g_{1},g_{2}$ and $g_{4}$.
Definite values of astronomical fundamental frequencies are thus determined.
Unfortunately, the uncertainty estimation of the geological frequencies is yet
not available.

Our density estimation of the frequencies is obtained by combining all our
samples in a conservative approach specified in section 3.4. Figure 12
compares the secular frequencies (right) and their combinations (left)
extracted from NH data with our density estimation around the same time in the
past. For $(g)_{i=1,3}$, both the geological and La2010d lie relatively close
to the main peak of their corresponding density while their values for $g_{4}$
are near the tail. The frequency combinations tell the same story, as the
geological values involving $g_{4}$ are off from the main peak. Yet, they are
all consistent with our density estimation as there is a non-negligible
possibility of finding a secular solution that agrees with the geological
data. It should be noted that certain frequencies are significantly
correlated, $g_{3}$ and $g_{4}$ for example. We cannot assume that they are
independent, therefore we have to calculate the density estimation of their
combination directly.

Given the unavailability of the uncertainty in geological frequencies, the
probability of finding the geological frequencies in a numerical orbital
solution cannot be directly obtained. However, we can use La2010d, which is
the solution from the complete model of the Solar System that matches the best
with NH data, as the benchmark for our secular statistics. There are several
criteria to determine how good a solution is, i.e. how well it could match
with the geological data. A simple and rather straightforward criterion that
we will use is $\delta=\sqrt{\frac{1}{4}\sum_{i=1}^{4}(g_{i}-g_{i}^{*})^{2}}$,
where $g_{i},g_{i}^{*}$ are the frequencies from astronomical solution and
geological data respectively. A better-suited solution will have smaller
$\delta$ and vice versa. We found that in the range from $-200$ Myr to $-220$
Myr, out of the 120 000 solutions, there are around 5000 solutions (accounting
for roughly $4.2\%$ of the total number) having smaller $s$ than that of
La2010d at -210 Myr, which is the value originally used to compared with the
geological data. It should be noted that La2010d is one of 13 available
complete solutions. The 95% confidence interval of the probability of
obtaining such a good matching solutions from the complete model of the Solar
System is thus $(1.37\%,33.31\%)$ (Wilson 1927). Therefore, with the criterium
$\delta$, our result is statistically compatible with that of (Olsen et al.
2019).

###  6.3 Libsack core

Laskar (1990, 1992) presented several secular resonances to explain the origin
of chaos in the Solar System. In particular, the argument of the resonance
$(s_{4}-s_{3})-2(g_{4}-g_{3})$ is currently in a librational state, that is

$(s_{4}-s_{3})-2(g_{4}-g_{3})=0$ (21)

and moves out to the rotational state around $-50$ Myr. The dynamics can even
switch to the librational state of a new resonance:

$(s_{4}-s_{3})-(g_{4}-g_{3})=0.$ (22)

This transition corresponds to a change from the 2:1 resonance to the 1:1
resonance of two secular terms $g_{4}-g_{3}$ and $s_{4}-s_{3}$. Because of
this resonance, the

Ma et al. (2017) found a sudden change of the period of a long cycle from 2.4
Myr to 1.2 Myr in the Libsack core of the Cretaceous basin from around $-90$
Myr to $-83$ Myr. This change was also visible in La2004 astronomical
solution, and the long cycle was attributed to frequency combination
$g_{4}-g_{3}$, visible from the spectrum of the eccentricity of the Earth.
Although the exact value before and especially after the transition is not
clear, the change of the period is visible from the band power of the core
(Fig. 1 of Ma et al. (2017)).

This change of $g_{4}-g_{3}$ observed in the Libsack core corresponds to a
transition from the resonance $(s_{4}-s_{3})-2(g_{4}-g_{3})$, which is the
resonance that the Solar System is currently at, to the resonance
$(s_{4}-s_{3})-(g_{4}-g_{3})$. With a large number of astronomical solutions,
we could understand better this phenomenon. The transition is usually very
fast: the frequency changes quickly to another value and then reverts back
also as quickly. Therefore, to study the transition we use a smaller window
for the frequency analysis, 10 Myr instead of 20 Myr. Figure 13 shows the
density estimation of $g_{4}-g_{3}$ at 90 Myr as well as 300 Myr in the past.
Both has a principle population in the range $[0.4,0.8]\
^{\prime\prime}/\text{yr}$, and a small but non-insignificant one centered
around $1.0\ ^{\prime\prime}/\text{yr}$, which corresponds to the small chunk
of $g_{4}$ centered at $18.2\ ^{\prime\prime}/\text{yr}$. The transition
observed is a sudden jump of frequency from the main population to the
secondary one of $g_{4}$, and therefore $g_{4}-g_{3}$ as well.

The size of the secondary population is defined as the proportion of the
solutions whose $g_{4}-g_{3}>0.9\ \arcsec/\text{yr}$ and denoted as
$\mathbb{P}(g_{4}-g_{3}>0.9\ ^{\prime\prime}/\text{yr})$. The rate of
transition is defined as the proportion of the solutions whose
$g_{4}-g_{3}>0.9\ ^{\prime\prime}/\text{yr}$ over 10 Myr and denoted as
$\mathbb{P}_{10Myr}(g_{4}-g_{3}>0.9\ \arcsec/\text{yr})$. Both are shown in
Fig 14. During the predictable period that is from now to $-50$ Myr, no
transition is observed. After $-50$ Myr, a transition can occur, its rate
rises until 100 Myr, when the percentage of secondary population stay
relatively stable at around $2.1\%\pm 0.5\%$ at a time, and the rate of the
transition could also be determined to be $8\%\pm 1\%$ every 10 Myr during
this period. At 80-90 Myr, when the transition was detected in the Libsack
core, with our numerical solutions, that rate of transition during this period
is found to be $7.7\%\pm 4.5\%$.

##  7 Conclusion

In this work, we give a statistical description of the evolution of the
fundamental frequencies of the Solar System beyond 60 Myr, that is beyond the
predictability horizon of the planetary motion, with the aim to quantify the
uncertainty induced by its chaotic behaviour. The base of our analysis is 120
000 orbital solutions of the secular model of the Solar System. The PDF of the
frequencies is estimated by the kernel density estimator, its uncertainty
evaluated by the moving block bootstrap; both methods are tested by numerical
experiments.

We have benchmarked the secular model by sampling the initial conditions in
different ways, then compare the density estimation of their solutions between
themselves and finally with the complete model. The results are two folds.
First, regardless of how initial conditions are sampled, their density
estimation will converge towards a single PDF; after this overlap, a robust
estimation is guaranteed. Secondly, the density estimation of the secular
model is compatible with that of the complete model of the Solar System. This
agreement means that results of the secular model, with superior computational
simplicity, can be used for application on geological data.

We observe that the density estimations of the fundamental frequencies can be
well fitted by skew Gaussian mixture models. The time evolution of the
parameters $\sigma_{0}^{2}$, related to the frequency variances, follows
power-law functions. Interestingly for the inner fundamental frequencies, the
exponents of such power laws are all smaller than 1, which indicates that they
undergo subdiffusion processes.

We show several examples of how this result can be used for geological
applications. First, the uncertainty of any astronomical frequency signal is
fully quantified, so that, for example, a proper quantitative response can be
now made for the question of how stable are the astronomical metronomes. With
this statistical framework, previous results from geological records beyond 60
Myr can also be interpreted with a more comprehensive approach. A more
quantitative answer, not only about the possibility but also about the
probability of occurrence of an astronomical signal in geological data can be
made. Apart from these direct applications, a more systematic approach could
fully make use of the density estimation of frequencies. The method TimeOpt by
Meyers & Malinverno (2018), for example, shows that it is possible to combine
both the uncertainty from astronomical signals and geological records to
derive an effective constraint for both. In fact, any similar Bayesian method
could use the density estimation of frequencies as proper priors.

###### Acknowledgements.

NH is supported by the PhD scholarship of CFM Foundation for Research. This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Advanced
Grant AstroGeo-885250) and from the French Agence Nationale de la Recherche
(ANR) (Grant AstroMeso ANR-19-CE31-0002-01).

## References

  * Bosq & Guégan (1995) Bosq, D. & Guégan, D. 1995, Statistics & probability letters, 25, 201
  * Boulila et al. (2018) Boulila, S., Vahlenkamp, M., De Vleeschouwer, D., et al. 2018, Earth and Planetary Science Letters, 486, 94
  * Calonico et al. (2018) Calonico, S., Cattaneo, M. D., & Farrell, M. H. 2018, Journal of the American Statistical Association, 113, 767
  * Charbonnier et al. (2018) Charbonnier, G., Boulila, S., Spangenberg, J. E., et al. 2018, Earth and Planetary Science Letters, 499, 266
  * Cheng et al. (2019) Cheng, G., Chen, Y.-C., et al. 2019, Electronic Journal of Statistics, 13, 2194
  * Efron (1979) Efron, B. 1979, Ann. Statist., 7, 1
  * Fermi (1949) Fermi, E. 1949, Physical review, 75, 1169
  * Gastineau & Laskar (2011) Gastineau, M. & Laskar, J. 2011, ACM Communications in Computer Algebra, 44, 194
  * Gradstein & Ogg (2020) Gradstein, F. & Ogg, J. 2020, in Geologic Time Scale 2020 (Elsevier), 21–32
  * Gradstein et al. (2012) Gradstein, F. M., Ogg, J. G., Schmitz, M., & Ogg, G. 2012, The geologic time scale 2012 (elsevier)
  * Gradstein et al. (2004) Gradstein, F. M., Ogg, J. G., Smith, A. G., Bleeker, W., & Lourens, L. J. 2004, Episodes, 27, 83
  * Hall et al. (1995) Hall, P., Horowitz, J. L., & Jing, B.-Y. 1995, Biometrika, 82, 561
  * Hang et al. (2018) Hang, H., Steinwart, I., Feng, Y., & Suykens, J. A. 2018, The Journal of Machine Learning Research, 19, 1260
  * Hart (1996) Hart, J. D. 1996, Journal of nonparametric statistics, 6, 115
  * Hartung et al. (2011) Hartung, J., Knapp, G., & Sinha, B. K. 2011, Statistical meta-analysis with applications, Vol. 738 (John Wiley & Sons)
  * Hastings (1970) Hastings, W. K. 1970, Biometrika, 57, 97
  * Kass et al. (1998) Kass, R. E., Carlin, B. P., Gelman, A., & Neal, R. M. 1998, The American Statistician, 52, 93
  * Kreiss & Lahiri (2012) Kreiss, J.-P. & Lahiri, S. N. 2012, in Handbook of statistics, Vol. 30 (Elsevier), 3–26
  * Kuffner et al. (2019) Kuffner, T. A., Lee, S. M.-S., & Young, G. A. 2019, arXiv preprint arXiv:1909.02662
  * Kunsch (1989) Kunsch, H. R. 1989, The annals of Statistics, 1217
  * Laskar (1985) Laskar, J. 1985, Astronomy and Astrophysics, 144, 133
  * Laskar (1988) Laskar, J. 1988, Astronomy and Astrophysics, 198, 341
  * Laskar (1989) Laskar, J. 1989, Nature, 338, 237
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Laskar (1992) Laskar, J. 1992, in Symposium-International Astronomical Union, Vol. 152, Cambridge University Press, 1–16
  * Laskar (1993) Laskar, J. 1993, Physica D: Nonlinear Phenomena, 67, 257
  * Laskar (2005) Laskar, J. 2005, in Hamiltonian Systems and Fourier Analysis: New Prospects For Gravitational Dynamics, ed. D. Benest, C. Froeschle, & E. Lega (Cambridge Scientific Publishers Ltd), –, arXiv: math/0305364
  * Laskar (2008) Laskar, J. 2008, Icarus, 196, 1
  * Laskar (2020) Laskar, J. 2020, in Geologic Time Scale 2020 (Elsevier), 139–158
  * Laskar et al. (2011a) Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011a, Astronomy & Astrophysics, 532, A89
  * Laskar & Gastineau (2009) Laskar, J. & Gastineau, M. 2009, Nature, 459, 817
  * Laskar et al. (2011b) Laskar, J., Gastineau, M., Delisle, J.-B., Farrés, A., & Fienga, A. 2011b, Astronomy & Astrophysics, 532, L4
  * Laskar et al. (2004) Laskar, J., Robutel, P., Joutel, F., et al. 2004, Astronomy & Astrophysics, 428, 261
  * Liu et al. (1992) Liu, R. Y., Singh, K., et al. 1992, Exploring the limits of bootstrap, 225, 248
  * Ma et al. (2017) Ma, C., Meyers, S. R., & Sageman, B. B. 2017, Nature, 542, 468
  * Maume-Deschamps (2006) Maume-Deschamps, V. 2006, Stochastics and Dynamics, 6, 535
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, The journal of chemical physics, 21, 1087
  * Meyers & Malinverno (2018) Meyers, S. R. & Malinverno, A. 2018, Proceedings of the National Academy of Sciences, 115, 6363
  * Milankovitch (1941) Milankovitch, M. 1941, Koniglich Serbische Akademice Beograd Special Publication, 132
  * Mogavero (2017) Mogavero, F. 2017, Astronomy & Astrophysics, 606, A79
  * Murray et al. (1985) Murray, N. W., Lieberman, M. A., & Lichtenberg, A. J. 1985, Physical Review A, 32, 2413
  * Olsen & Kent (1999) Olsen, P. E. & Kent, D. V. 1999, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357, 1761
  * Olsen et al. (2019) Olsen, P. E., Laskar, J., Kent, D. V., et al. 2019, Proceedings of the National Academy of Sciences, 201813901
  * Parzen (1962) Parzen, E. 1962, The annals of mathematical statistics, 33, 1065
  * Robinson (1983) Robinson, P. M. 1983, Journal of Time Series Analysis, 4, 185
  * Rosenblatt (1956) Rosenblatt, M. 1956, Proceedings of the National Academy of Sciences of the United States of America, 42, 43
  * Silverman (1986) Silverman, B. W. 1986, Density estimation for statistics and data analysis, Vol. 26 (CRC press)
  * Wilson (1927) Wilson, E. B. 1927, Journal of the American Statistical Association, 22, 209


# Resistive MHD Modelling of Quasi-Single Helicity State in the KTX Regimes

###### Abstract

The potential formation of the quasi-single-helicity (QSH) state in the Keda
Torus eXperiment (KTX) is investigated in resistive MHD simulations using the
NIMROD code. We focus on the effects of finite resistivity on the mode
structure and characteristics of the dominant linear and nonlinear resistive
tearing-mode in a finite $\beta$, cylindrical configuration of reversed field
pinch model for KTX. In the typical resistive regimes of KTX where Lundquist
number $S=5\times 10^{4}$, the plasma transitions to a steady QSH state after
evolving through an initial transient phase with multiple helicities. The
dominant mode of the QSH state develops from the dominant linear tearing mode
instability. In lower $\beta$ regime, the QSH state are intermittent and short
in duration; in higher $\beta$ regime, the QSH state persists for a longer
time and should be more observable in experiment.

###### pacs:

52.30.Cv, 52.55.HC, 52.55.Tn, 52.65.Kj

Bing Luoa,b, Ping Zhua,b,c, Hong Lib, Wandong Liub, and KTX teamb

  * June 2017

##  1 Introduction

The reversed-field pinch (RFP) is a toroidal magnetic confinement device. Its
main difference from tokamak is its toroidal field $B_{\phi}$, which is of the
same order of magnitude as the poloidal field $B_{\theta}$ and becomes
reversed near plasma edge [1]. Since the safety factor $q(r)$ in RFP is always
less than one, the RFP plasma can easily become unstable to resistive-kink
modes.

During the last two decades, the QSH states have been discovered in four
different RFP devices [2-5]: RFX in Padova [6], MST at the University of
Wisconsin-Madison [7], EXTRAP T2R in Stockholm [8] and RELAX in Kyoto [9].
Significant efforts within the RFP community have been devoted to the study of
quasi-single helicity (QSH) states ever since [10-12]. The main characteristic
of the QSH state is the presence of an inner resonant dominant mode with
poloidal periodicity $m=1$ and toroidal periodicity $n_{D}\geq{3R_{0}/2a}$.
Its amplitude is typically several times larger than those of secondary modes.
As a consequence, the QSH state plasma core has a helical deformation that
generates a 3D structure like the stellarator [11]. QSH states with either
Double Axis (DAx) or Single Helical Axis (SHAx) have been found in experiments
[13,14]. The appearance of the QSH state reduces magnetic fluctuations and
improves the RFP confinement, especially for the SHAx state, which has an
inner electron temperature transport barrier [14,15]. The RFX-mod experimental
statistics shows that the QSH state properties correlate with Lundquist number
S, and higher plasma current leads to longer persistence of the QSH state
[15-17]. Lately, the helical magnetic perturbation from the external active
coils has been used to excite QSH with specific dominant modes in RFX-mod
[18].

The proposition that RFP plasma could exist in a pure single helicity SH state
was put forward in 1983, based on two-dimensional numerical simulations
[19,20], The 2D simulations of pure (SH) states show the existence of two
topologies for the corresponding magnetic surfaces: with or without a magnetic
separatrix [12]. Later, numerical simulations of the RFP using 3D resistive
magnetohydrodynamics codes SPECYl [21,22] also find MH and QSH states. Those
simulations indicate the transition from MH to QSH state is controlled by the
Hartmann number $H\propto\frac{1}{\sqrt{{\eta}{\nu}}}$ [23], where $\eta$ is
resistivity, and $\nu$ is viscosity. If H is large, the system is in the MH
state; when $H<3000$, the system displays temporal intermittency with a
laminar phase of the QSH state [24,25]. The scheme of using the helical
boundary condition generated by active coils to excite and control the
specific dominate modes of QSH state in RFP has been studied using nonlinear
3D MHD simulations [18,26,27].

In this work, we use the 3D full MHD code NIMROD [28] to evaluate the
possibilities of achieve QSH state in the newly constructed RFP device KTX,
which is a middle-sized torus, with major radius $R_{0}=1.4$ m, minor radius
$a=0.4$ m, and a designed maximum plasma current of 1 MA [29]. With 96 saddle
coils installed outside of its copper shell, KTX is mainly designed to explore
novel magnetic confinement regimes and advanced feedback control schemes for
MHD instabilities, along with 3D physics. For the first time, our NIMROD
simulations have demonstrated the potential formation of QSH state in the RFP
configuration of KTX regimes. In particular, we find the connection between
the QSH state and the dominant linear tearing mode on KTX. We further
investigate the condition of QSH state formation and find a new parameter
dependence of QSH, which is different from previous studies.

The rest of the paper is organized as follows. In Section 2, the single-fluid
MHD model and a typical KTX equilibrium are described; in Section 3, the
linear MHD instability of KTX is analysed including the growth rate and mode
structure; in Section 4, the 3D QSH state is demonstrated in a nonlinear
simulation, where a saturated magnetic island is observed to form and persist.
Finally, Section 5 is devoted to a summary and discussion.

##  2 MHD model and equilibrium

The numerical simulations in this work are performed with the single-fluid MHD
model implemented in the NIMROD code. The single-fluid MHD equations can be
written as follows:

$\displaystyle~{}\frac{\partial\rho}{\partial{t}}+\nabla\cdot(\rho\textbf{u})=0$
(1)
$\displaystyle~{}\rho(\frac{\partial}{\partial{t}}+\textbf{v}\cdot\nabla)\textbf{v}=\textbf{J}\times\textbf{B}-\nabla
p+\rho\nu\nabla^{2}\textbf{v}$ (2)
$\displaystyle~{}\frac{{N}}{(\gamma-1)}(\frac{\partial}{\partial{t}}+\textbf{v}\cdot\nabla)T=-p\nabla\cdot\textbf{v}-\nabla\cdot\textbf{q}$
(3)
$\displaystyle~{}\frac{\partial\textbf{B}}{\partial{t}}=\nabla\times(\textbf{v}\times\textbf{B}-\eta\textbf{J})$
(4) $\displaystyle~{}\nabla\times\textbf{B}=\mu_{0}\textbf{J}$ (5)
$\displaystyle\textbf{q}=-N(\chi_{\parallel}\nabla_{\parallel}T+\chi_{\perp}\nabla_{\perp}T)$
(6)

where $\rho$, $N$, v, $p$, j, B, q, $\gamma$, $\eta$, $\nu$, and
$\chi_{\parallel}$ ($\chi_{\perp}$) are the plasma mass density, number
density, velocity, pressure, current density, magnetic field, heat flux,
specify heat ratio, resistivity, viscosity, and parallel (perpendicular)
thermal conductivity, respectively. The Lundquist number is defined as
$S$=$\tau_{R}$/$\tau_{A}$ with Alfven time $\tau_{A}$=$a$/$V_{A}$
($V_{A}$=$B_{0}/{\sqrt{\mu_{0}\rho)}}$, resistivity diffusion time
$\tau_{R}$=$\mu_{0}a^{2}$/$\eta$, and $B_{0}$ is the initial magnetic field
strength on the magnetic axis.

In this work, we consider the low-$\beta$ regime of KTX plasmas, in contrast
to the zero $\beta$ regime considered in the SPECY1 simulations [22,25]. Its
initial configuration can be approximated as a cylindrical force-free RFP
equilibrium, which can be specified using the following q profile (Fig. 1)
[30]:

$q(r)=q_{0}[1-1.8748(r)^{2}+0.83232(r)^{4}]\\\ $ (7)

The above RFP equilibrium has been commonly adopted in previous studies
[21,31,32]. For the KTX device, the equilibrium parameters at the magnetic
axis are $q_{0}=2a/3R_{0}=0.18$, $\Theta_{0}=\mu_{0}aJ_{\phi 0}/2B_{0}=1.46$,
$B_{0}=1.0$ T, $N=1.2\times 10^{19}$ $m^{-3}$, $J_{\phi_{0}}=6.15\times
10^{3}$ $MA/m^{2}$, resistivity $\eta$ and viscosity $\nu$ share the same
normalized radial profile
$\frac{{\eta(r)}}{\eta_{0}}=\frac{{\nu(r)}}{\nu_{0}}$=${[1+(r/a)^{20}]}^{2}$.
In our simulations, $3\times 10^{-7}<\eta_{0}<3\times 10^{-4}$, $1\times
10^{4}<S<1\times 10^{7}$, $\nu_{0}=1\times 10^{-3}$, and $10^{3}<H<7\times
10^{4}$. For an arbitrary initial perturbation, the MHD equations (1)-(6) are
numerically solved in the cylindrical geometry, along with the periodic
boundary conditions in the axial direction and the ideal wall boundary
conditions in the radial direction.

##  3 Linear MHD instability analysis

For benchmark purpose, two kinds of grids for cylindrical geometry, namely
“circular” and “rectangular”, are been set up for the simulations (Fig. 2).
Simulations results from the “circular” and the “rectangular” grids have been
compared to verify their correctness. For the same equilibrium physical
parameters, the time evolutions of the magnetic energy of the m=1, n=9 mode
yield the same growth rate from those two grids. The mode structures from the
two grids are also consistent with each other (Fig. 3). From now on, all
results reported in this paper are obtained from simulations on the “circular”
grid.

The relationship between the growth rate $\gamma$ and the toroidal mode number
$n$ shown in Fig.4a for different Lundquist number $S$, clearly indicates that
the instability is located in plasma core, and the n=8 mode is the most
unstable for this equilibrium. The normalized growth rate $\gamma\tau_{A}$ as
a function of $S$ for the n=8 mode shown in Fig.4b, suggests that the
instability is likely the resistive tearing mode in nature, since the
asymptotic scaling $\gamma\tau_{A}$ $\propto$ $S^{-0.6}$ obtained from
calculation agrees with the traditional theory [33,34].

The characteristic perturbations and the structure of instability are further
analysed. Fig. 5 shows the radial profiles of the normalized perturbation
$v_{r}$ and $b_{r}$ along the mid-plane in the cross section of the n=8 mode
for $S=1\times 10^{6}$. Here $v_{r}$ and $b_{r}$ are the radial components of
the perturbed velocity and magnetic filed, respectively. The resonant surface
locates at r=0.163 m, which is the location of the 1/8 resonant surface based
on the q profile. Across that resonant surface, $v_{r}$ is anti-symmetric, and
$b_{r}$ is symmetric, which is same as the parity of tearing-mode [35]. The 2D
mode structure in the poloidal plane can be more clearly seen in Fig. 6, which
indicates that the instability is dominated by the m=1 tearing mode.

##  4 QSH state from nonlinear simulations

For nonlinear simulations, we use the parameter ${\eta}_{0}=6.0\times 10^{-5}$
$\Omega m$, ${\nu}_{0}=1\times 10^{-3}$ $m^{2}/s$, density $N=1.2\times
10^{19}$ $m^{-3}$, $B_{0}=1.0$ T, $H=1.5\times 10^{3}$, $\beta=0.05$, and
magnetic Lundquist number $S=5.0\times 10^{4}$ on the equilibrium magnetic
axis. The nonlinear evolution of magnetic energy is dominated by $n=7-12$
modes (Fig.7a). At t=0.35 ms, the instability mode evolves towards a saturated
phase. During $t=0.35-1.1$ ms, the modes n=7, 8, and 9 compete strongly, but
at t=1.1 ms, the energy of the perturbation begins to reside in the n=8 mode.
After t=2.0 ms, with the decrease in the magnetic energy of the n=7 and 9
modes, the plasma relaxes to the QSH state, which is dominated by the n=8
mode. After t=5.6 ms, the plasma transitions from the QSH to the MH state. It
is customary to describe the width of the toroidal spectrum of m=1 mode by
using the $N_{s}$ parameter, which is defined as
$N_{s}=[\sum_{n=1}(W_{1,n}/\sum_{n^{\prime}=1}{W_{1,n^{\prime}}})^{2}]^{-1}$,
where $W_{1,n}$ is the magnetic energy of the (m=1, n) mode. A pure SH
spectrum corresponds to $N_{s}=1$; for pure QSH state, $N_{s}$ may be less
than 2. The time evolution of $N_{s}$ verifies that during those periods of
QSH as indicated from magnetic energy evolution of each in modes(Fig. 7a), the
$N_{s}$ number is indeed less than 2 (Fig. 7b).

The emergence of QSH can be further demonstrated in the time history of the
Poincare plot in the poloidal plane(Fig.8). Here the plot at t=0 ms
corresponds to the equilibrium state, t=1.5 ms the MH state, and the plot for
$t=2.3-3.9$ ms correspond to a phase of the QSH state. A clearly shaped
magnetic island appears and persists during the QSH state. Around t=5.6-6.2
ms, the island disappears and the plasma transitions from the QSH to the MH
state. These Poincare plots show that the magnetic stochastic region is
located in the plasma interior from $r=0.1$ to $0.2$ m, corresponding to the
location of the $n=7-12$ inner modes. The size of the magnetic island is
determined by the ratio of the secondary modes within a certain range. The
amplitudes of the m=1 modes for different toroidal numbers are compared
between the MH states at t=0.71 ms, 1.5 ms, 6.2 ms and the QSH states at t=3.1
ms, t=3.9 ms and t=5.16 ms (Fig.9), which again shows that the QSH state is
dominated by the n=8 mode. Since the n=8 mode has the largest linear growth
rate, the nonlinear simulations results suggest that the QSH state may have
developed mainly from the dominant linear tearing mode.

The formation mechanism for QSH has long been a subject of theory and
experimental studies. SPECY1 code simulation indicates that Hartmann number
$H$ and helical magnetic perturbation can determine the QSH state [24-26].
Statistical results from RFX-mod experiment indicate that plasma current and
Lundquist number $S$ correlate to the presence of the QSH state. In contrast,
our simulation suggests that $\beta$ is the governing parameter for the
formation of the QSH state in the KTX configuration and parameter regimes. The
time dependence of the spectral spread $N_{s}$ for different $\beta$ shows
clearly the correlation (Fig.10). For lower $\beta$, the QSH state appears
intermittently in a shorter duration of time; for higher $\beta$, the QSH
state persists in a longer duration of time. Therefore, it might be possible
to obtain a QSH state through auxiliary heating of electrons even in regimes
of low plasma current.

##  5 Summary and discussion

In summary, for the first time, NIMROD simulations find the formation of QSH
state in KTX regimes with Lundqusit number $S=5\times 10^{4}$. The simulation
starts from a MH state, which transitions to a QSH state when secondary modes
decrease in amplitude. The simulation indicates that the dominate mode of QSH
state may develop from the linear mode that has the maximum growth rate. The
plasma $\beta$ appears to be the key parameter that governs the emergence and
duration of the QSH state. QSH form when $\beta$ is above certain threshold
and a higher $\beta$ leads to a longer duration of QSH.

An effective method to achieve and maintain the QSH state may be the auxiliary
heating that enhances the electron temperature, even in regimes of low plasma
current. However, in the transition between the MH phase and the QSH phase, it
is not clear what mechanism suppresses the magnetic fluctuations of the
secondary modes and maintain the dominate mode. In addition, how two-fluid
effects and anomalous viscosity in high plasma current regimes may affect the
QSH state should be also studied in future work.

##  6 Acknowledgment

This research was supported by the National Magnetic Confinement Fusion
Program of China under Grant Nos. 2014GB124002 and 2015GB101004, the National
Natural Science Foundation of China under Grant Nos. 11635008, the U.S.
Department of Energy Grant Nos. DE-FG02-86ER53218 and DE-FC02-08ER54795, and
the 100 Talent Program of the Chinese Academy of Science. We are grateful for
the support from the NIMROD and KTX teams. This research used resources from
the Supercomputing Center of University of Science and Technology of China,
and the National Energy Research Scientific Computing Center, a DOE Office of
Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231.

##  7 Reference

## References

  * [1] Bodin H and Newton A. 1980 Nucl. Fusion 20 1255
  * [2] P. Martin. et al 2000 Phys. Plasmas 7 1984
  * [3] L. Marrelli. et al 2002 Phys. Plasmas 9 2868
  * [4] L. Frassinetti. et al 2007 Phys. Plasmas 14 112510
  * [5] R. IKEZOE. et al 2008 Plasma Fusion Res. 3 029
  * [6] G. Rostagni 1995 Fusion Eng. Des 25 301
  * [7] Dexter R.N. et al 1991 Fusion Technol 19 131
  * [8] Brunsell P.R. et al 2001 Plasma Phys. Control. Fusion 43 1457
  * [9] S. Masamune et al 2007 J. Phys. Soc. Jpn 76 123501
  * [10] D.F. Escande. et al 2000 Plasma Phys. Control Fusion 42 B243
  * [11] D.F. Escande. et al 2000 Phys. Rev. Lett 85 1662
  * [12] D.F. Escande. et al 2000 Phys. Rev. Lett 85 3169
  * [13] R. Lorenzini. et al 2008 Phys. Rev. Lett 101 025005
  * [14] P. Piovesan. et al 2008 Nucl. Fusion 49 085036
  * [15] M. Gobbin. et al 2011 Phys.Rev.Lett 106 025001
  * [16] L. Carraro. et al 2009 Nucl. Fusion 49 055009
  * [17] R. Lorenzini. et al 2009 Nature Phys. 5 570
  * [18] D.Bonfiglio. et al 2013 Phys. Rev. Lett 111 085002
  * [19] E.J.Caramana. et al 1983 Phys. Fluids 26 1305
  * [20] D.D.Schnack. et al 1985 Phys. Fluids 28 321
  * [21] S.Cappello. et al 1992 Phys. Fluids B 4 611
  * [22] S.Cappello. et al 1996 Nucl. Fusion 36 571
  * [23] D. Montgomery. et al 1992 Plasma Phys. Control. Fusion 34 1157
  * [24] S.Cappello. et al 2000 Phys. Rev. Lett 85 3838
  * [25] S.Cappello. et al 2004 Plasma Phys. Control. Fusion 46 B313
  * [26] M.Veranda et al 2013 Plasma Phys. Control. Fusion 55 074015
  * [27] D.Bonfiglio et al 2015 Plasma Phys. Control. Fusion 57 044001
  * [28] Sovinec et al 2004 J. Comput. Phys. 195 355
  * [29] Wandong Liu. et al 2014 Plasma Phys. Control. Fusion 56 094009
  * [30] D.D.Schnack et al 1987 J. Comput. Phys 70 330
  * [31] Finn J. M et al 1992 Phys. Fluids B 4 1262
  * [32] M.Onofri et al 2008 Phys. Rev. Lett 101 255002
  * [33] H.P. Furth et al 1963 Phys. Fluids 6 459
  * [34] H.P. Furth et al 1973 Phys. Fluids 16 1054
  * [35] J.P Freidberg et al 1981 J. Plasma Phys. 26 177

Figure 1:  Radial profile of the safety factor q(r) of the initial equilibrium
for KTX.

Figure 2: Two kinds of NIMROD grids for a cylindrical configuration: (left)
“rectangular” grid and (right) “circular” grid.

Figure 3: Contours of the radial, poloidal, and toroidal components of the
perturbed magnetic field, i.e. $b_{r}$, $b_{\theta}$, and $b_{\phi}$ for the
(m=1,n=9) mode in the $\phi=0$ plane of the “circular” grid (a-c) and the
$\theta=0$ plane of the “rectangular” grid (d-f).

Figure 4:  (a) Linear growth rate as a function of toroidal mode number n for
different values of Lundquist number S. (b) Linear growth rate as a function
of Lundquist number S for the n=8 mode.

Figure 5:  Radial profiles as function of the minor radius r for (a) the
normalized $v_{r}$ ($\theta=90^{o}$) and (b) the normalized $b_{r}$
($\theta=180^{o}$) of the (m=1, n=8) mode with $S=10^{6}$. A singular layer
locates around r = 0.163 m.

Figure 6:  Contours of (a) $v_{r}$ and (b) $b_{r}$ of the (m=1, n=8) mode with
$S=10^{6}$ in poloidal plane.

Figure 7:  (a) The magnetic energies of different toroidal mode components as
functions of time, and (b) evolution of the spectral spread $N_{s}$ for
$S=5.0\times 10^{4}$ and $\beta=0.05$

Figure 8: Snapshots of the Poincare plots in poloidal plane for the simulation
case shown in Figure 7: $(a)-(f)$ correspond to t=0 ms, 1.5 ms, 2.3 ms, 3.1
ms, 3.9 ms, and 6.2 ms, respectively.

Figure 9:  The spectra in toroidal mode n for the m=1 mode at time t=0.71 ms,
1.5 ms and 6.2 ms in the MH state (a-c), and 3.1 ms, 3.9 ms, and 5.16 ms in
the QSH state (d-f).

Figure 10:  Spectral spread in $N_{s}$ as a function of time for different
$\beta$ with $S=5.0\times 10^{4}$.