11institutetext:  ASD/IMCCE, CNRS-UMR8028, Observatoire de Paris, PSL
University, Sorbonne Université, 77 Avenue Denfert-Rochereau, 75014 Paris,
France
11email<EMAIL_ADDRESS>
# High-order regularised symplectic integrator for collisional planetary
systems

Antoine C. Petit  Jacques Laskar  Gwenaël Boué  Mickaël Gastineau

(Accepted XXX. Received YYY; in original form ZZZ)

We present a new mixed variable symplectic (MVS) integrator for planetary
systems that fully resolves close encounters. The method is based on a time
regularisation that allows keeping the stability properties of the symplectic
integrators while also reducing the effective step size when two planets
encounter. We used a high-order MVS scheme so that it was possible to
integrate with large time-steps far away from close encounters. We show that
this algorithm is able to resolve almost exact collisions (i.e. with a mutual
separation of a fraction of the physical radius) while using the same time-
step as in a weakly perturbed problem such as the solar system. We demonstrate
the long-term behaviour in systems of six super-Earths that experience strong
scattering for 50 kyr. We compare our algorithm to hybrid methods such as
MERCURY and show that for an equivalent cost, we obtain better energy
conservation.

###### Key Words.:

Celestial mechanics, Planets and satellites: dynamical evolution and
stability, Methods: numerical

##  1 Introduction

Precise long-term integration of planetary systems is still a challenge today.
The numerical simulations must resolve the motion of the planets along their
orbits, but the lifetime of a system is typically billions of years, resulting
in computationally expensive simulations. In addition, because of the chaotic
nature of planetary dynamics, statistical studies are often necessary, which
require running multiple simulations with close initial conditions (Laskar &
Gastineau, 2009). This remark is particularly true for unstable systems that
can experience strong planet scattering caused by close encounters.

There is therefore considerable interest in developing fast and accurate
numerical integrators, and numerous integrators have been developed over the
years to fulfill this task. For long-term integrations, the most commonly used
are symplectic integrators. Symplectic schemes incorporate the symmetries of
Hamiltonian systems, and as a result, usually conserve the energy and angular
momentum better than non-symplectic integrators. In particular, the angular
momentum is usually conserved up to a roundoff error in symplectic
integrators.

Independently, Kinoshita et al. (1991) and Wisdom & Holman (1991) developed a
class of integrators that are often called mixed variable symplectic (MVS)
integrators. This method takes advantage of the hierarchy between the
Keplerian motion of the planets around the central star and the perturbations
induced by planet interactions. It is thus possible to make accurate
integrations using relatively large time-steps.

The initial implementation of Wisdom & Holman (1991) is a low-order
integration scheme that still necessitates small time-steps to reach machine
precision. Improvements to the method have since been implemented. The first
category is symplectic correctors (Wisdom et al., 1996; Wisdom, 2006). They
consist of a modification of the initial conditions to improve the scheme
accuracy. Because it is only necessary to apply them when an output is
desired, they do not affect the performance of the integrator. This approach
is for example used in WHFAST (Rein & Tamayo, 2015). The other approach is to
consider higher order schemes (McLachlan, 1995b; Laskar & Robutel, 2001;
Blanes et al., 2013). High-order schemes permit a very good control of the
numerical error by fully taking advantage of the hierarchical structure of the
problem. This has been used with success to carry out high-precision long-term
integrations of the solar system (Farrés et al., 2013).

The principal limitation of symplectic integrators is that they require a
fixed time-step (Gladman et al., 1991). If the time-step is modified between
each step, the integrator remains symplectic because each step is symplectic.
However, the change in time-step introduces a possible secular energy drift
that may reduce the interest of the method. As a consequence, classical
symplectic integrators are not very adapted to treat the case of systems that
experience occasional close encounters where very small time-steps are needed.

To resolve close encounters, Duncan et al. (1998) and Chambers (1999) provided
solutions in the form of hybrid symplectic integrators. Duncan et al. (1998)
developed a multiple time-step symplectic integrator, SYMBA, where the
smallest time-steps are only used when a close encounter occurs. The method is
limited to an order two scheme, however. The hybrid integrator MERCURY
(Chambers, 1999) moves terms from the perturbation step to the Keplerian step
when an interaction between planets becomes too large. The Keplerian step is
no longer integrable but can be solved at numerical precision using a non-
symplectic scheme such as Burlisch-Stoer or Gauss-Radau. However, the switch
of numerical method leads to a linear energy drift (Rein & Tamayo, 2015).
Another solution to the integration of the collisional $N$-body problem has
been proposed in (Hernandez & Bertschinger, 2015; Hernandez, 2016; Dehnen &
Hernandez, 2017). It consists of a fixed time-step second-order symplectic
integrator that treats every interaction between pairs of bodies as Keplerian
steps.

Another way to build a symplectic integrator that correctly regularises close
encounters is time renormalisation. Up to an extension of the phase space and
a modification of the Hamiltonian, it is indeed always possible to modify the
time that appears in the equations of motion. As a result, the real time
becomes a variable to integrate. Providing some constraints on the
renormalisation function, it is possible to integrate the motion with a fixed
fictitious time-step using an arbitrary splitting scheme. Here we show that
with adapted time renormalisation, it is possible to resolve close encounters
accurately.

While time renormalisation has not been applied in the context of close
encounters of planets, it has been successful in the case of perturbed highly
eccentric problems (Mikkola, 1997; Mikkola & Tanikawa, 1999; Preto & Tremaine,
1999; Blanes & Iserles, 2012), see (Mikkola, 2008) for a general review. We
here adapt a time renormalisation proposed independently by Mikkola & Tanikawa
(1999) and Preto & Tremaine (1999). We show that it is possible to use the
perturbation energy to monitor close encounters in the context of systems with
few planets with similar masses. We are able to define a MVS splitting that
can be integrated with any high-order scheme.

We start in section 2 by briefly recalling the basics of the symplectic
integrator formalism. In section 3 we present the time renormalisation that
regularises close encounters, and we then discuss the consequence of the
renormalisation on the hierarchical structure of the equations (section 4). In
section 5 we numerically demonstrate over short-term integrations the
behaviour of the integrator at close encounter. We then explain (section 6)
how our time regularisation can be combined with the perihelion regularisation
proposed by (Mikkola, 1997). Finally, we show the results of long-term
integration of six planet systems in the context of strong planet scattering
(section 7) and compare our method to a recent implementation of MERCURY
described in (Rein et al., 2019), to SYMBA, and to the non-symplectic high-
order integrator IAS15 (section 8).

##  2 Splitting symplectic integrators

We consider a Hamiltonian $H({\bf p},{\bf q})$ that can be written as a sum of
two integrable Hamiltonians,

$H({\bf p},{\bf q})=H_{0}({\bf p},{\bf q})+H_{1}({\bf p},{\bf q}).$ (1)

A classical example is given by $H_{0}=T({\bf p})$ and $H_{1}=U({\bf q}),$
where $T({\bf p})$ is the kinetic energy and $U({\bf q})$ is the potential
energy. In planetary dynamics, we can split the system as $H_{0}=K({\bf
p},{\bf q}),$ where $K$ is the sum of the Kepler problems in Jacobi
coordinates (_e.g._ Laskar, 1990) and ${H_{1}=H_{\mathrm{inter}}({\bf q})}$ is
the interaction between the planets.

Using the Lie formalism (_e.g._ Koseleff, 1993; Laskar & Robutel, 2001), the
equation of motion can be written

$\frac{\mathrm{d}z}{\mathrm{d}t}=\\{H,z\\}=L_{H}z,$ (2)

where $z=({\bf p},{\bf q})$, $\\{\cdot,\cdot\\}$ is the Poisson bracket111 We
use the convention ${\\{f,g\\}=\sum_{i}\frac{\partial f}{\partial
p_{i}}\frac{\partial g}{\partial q_{i}}-\frac{\partial f}{\partial
q_{i}}\frac{\partial g}{\partial p_{i}}}$. , and we note
$L_{f}=\\{f,\cdot\\}$, the Lie differential operator. The formal solution of
Eq. (2) at time $t=\tau+t_{0}$ from the initial condition $z(t_{0})$ is

$z(\tau+t_{0})=\exp(\tau
L_{H})z(t_{0})=\sum_{k=0}^{+\infty}\frac{\tau^{k}}{k!}L_{H}^{k}z(t_{0}).$ (3)

In general, the operators $L_{H_{0}}$ and $L_{H_{1}}$ do not commute, so that
$\exp(\tau L_{H})\neq\exp(\tau L_{H_{0}})\exp(\tau L_{H_{1}}).$ However, using
the Baker-Campbell-Hausdorff (BCH) formula, we can find the coefficients
$a_{i}$ and $b_{i}$ such that

$\exp(\tau(L_{H}+L_{H_{\mathrm{err}}}))=\mathcal{S}(\tau)=\prod_{i=1}^{N}\exp(a_{i}\tau
L_{H_{0}})\exp(b_{i}\tau L_{H_{1}}),$ (4)

where $H_{\mathrm{err}}=\mathrm{O}(\tau^{r})$ is an error Hamiltonian
depending on $H_{0}$, $H_{1}$, $\tau$ and the coefficients $a_{i}$ and
$b_{i}$.

Because $H_{0}$ and $H_{1}$ are integrable, we can explicitly compute the
evolution of the coordinates $z$ under the action of the maps $\exp(\tau
L_{H_{0}})$ and $\exp(\tau L_{H_{1}})$. The map $\mathcal{S}(\tau)$ is
symplectic because it is a composition of symplectic maps. Moreover,
$\mathcal{S}(\tau)$ exactly integrates the Hamiltonian $H+H_{\mathrm{err}}$.

If there is a hierarchy in the Hamiltonian $H$ in the sense that
$|H_{1}/H_{0}|\simeq\varepsilon\ll 1$, we can choose the coefficients such
that the error Hamiltonian is of the order of

$\sum_{i=1}^{n}\mathrm{O}(\tau^{r_{i}}\varepsilon^{i}),$ (5)

(see McLachlan, 1995b; Laskar & Robutel, 2001; Blanes et al., 2013; Farrés et
al., 2013). For small $\varepsilon$ and $\tau$, the solution of
$H+H_{\mathrm{err}}$ is very close to the solution of $H$. In particular, it
is thought that the energy error of a symplectic scheme is bounded. Because
$H_{\mathrm{err}}$ depends on $\tau$, a composition of steps
$\mathcal{S}(\tau)$ also has this property if the time-step is kept constant.
Otherwise, the exact integrated dynamics changes at each step, leading to a
secular drift of the energy error.

In planetary dynamics, we can split the Hamiltonian such that $H_{0}$ is the
sum of the Keplerian motions in Jacobi coordinates and $H_{1}$ is the
interaction Hamiltonian between planets, which only depends on positions and
thus is integrable (_e.g._ Laskar, 1990). This splitting naturally introduces
a scale separation $\varepsilon$ given by

$\varepsilon=\frac{\sum_{k=1}^{N}m_{k}}{m_{0}}$ (6)

where $N$ is the number of planets, $m_{k}$ is the mass of the $k$-th planet,
and $m_{0}$ is the mass of the star. If the planets remain far from each
other, $H_{1}$ is always $\varepsilon$ small with respect to $H_{0}$.

The perturbation term is of the order of $\varepsilon/\Delta,$ where $\Delta$
is the typical distance between the planets in units of a typical length of
the system. During close encounters, $\Delta$ can become very small, and the
step size needs to be adapted to $\varepsilon/\Delta_{\min}$ Here
$\Delta_{\min}$ is the smallest expected separation between planets normalised
by a typical length of the system.

##  3 Regularised Hamiltonian

###  3.1 General expression

In order to construct an adaptive symplectic scheme that regularises the
collisions, we extend the phase space and integrate the system with a
fictitious time. Let $s$ be such that

$\mathrm{d}t=g({\bf p},p_{t},{\bf q})\mathrm{d}s,$ (7)

where $g$ is a function to be determined and $p_{t}$ is the conjugated
momentum to the real time $t$ in the extended phase space. In order to have an
invertible function $t(s)$, we require $g$ to be positive. We consider the new
Hamiltonian $\Gamma$ defined as

$\Gamma({\bf p},p_{t},{\bf q},t)=g({\bf p},p_{t},{\bf q})\left(H({\bf p},{\bf
q})+p_{t}\right).$ (8)

$\Gamma$ does not depend on $t,$ therefore $p_{t}$ is a constant of motion.
The equations of motion of this Hamiltonian are

$\frac{\mathrm{d}t}{\mathrm{d}s}=\\{\Gamma,t\\}=g({\bf p},p_{t},{\bf
q})+\frac{\partial g}{\partial p_{t}}({\bf p},p_{t},{\bf q})\left(H({\bf
p},{\bf q})+p_{t}\right)$ (9)

and for all function $f(z)$

$\frac{\mathrm{d}f(z)}{\mathrm{d}s}=\\{\Gamma,f(z)\\}=g({\bf p},p_{t},{\bf
q})\\{H,f\\}+\left(H+p_{t}\right)\\{g,f(z)\\}.$ (10)

In general, $H$ is not a constant of motion of $\Gamma$. We have

$\frac{\mathrm{d}H}{\mathrm{d}s}=\\{\Gamma,H\\}=\left(H+p_{t}\right)\\{g,H\\}.$
(11)

If we choose initial conditions $z_{0}$ such that $p_{t}=-H(z_{0}),$ we have
$\left.\Gamma\right|_{t=0}=0$. Because $\Gamma$ is constant and $g$ is
positive, we deduce from equation (3) that we have at all times

$H+p_{t}=0.$ (12)

Because $p_{t}$ is also a constant of motion, $H$ is constant for all times.
We can simplify the equations of motion (9) and (10) into

$\displaystyle\frac{\mathrm{d}t}{\mathrm{d}s}$ $\displaystyle=g({\bf
p},p_{t},{\bf q})$ $\displaystyle\frac{\mathrm{d}f}{\mathrm{d}s}(z)$
$\displaystyle=g({\bf p},p_{t},{\bf q})\\{H,f(z)\\}.$ (13)

On the manifold $p_{t}=-H(t_{0})$, equations (13) describe the same motion as
equation (2). We call them the regularised equations.

We now wish to write $\Gamma$ as a sum of two integrable Hamiltonians such as
in section 2. Based on previous works (Preto & Tremaine, 1999; Mikkola &
Tanikawa, 1999; Blanes & Iserles, 2012), we write

$H+p_{t}=(H_{0}+p_{t})-(-H_{1}),$ (14)

for $H=H_{0}+H_{1}$ and we define $g$ as

$g({\bf p},p_{t},{\bf q})=\frac{f(H_{0}+p_{t})-f(-H_{1})}{H_{0}+p_{t}+H_{1}},$
(15)

where $f$ is a smooth function to be determined. $g$ is the difference
quotient of $f$ and is well defined when ${H_{0}+p_{t}+H_{1}\to 0}$. We have

$\left.g({\bf p},p_{t},{\bf
q})\right|_{H+p_{t}=0}=f^{\prime}(H_{0}+p_{t})=f^{\prime}(-H_{1}).$ (16)

With this choice of $g$, the Hamiltonian $\Gamma$ becomes

$\Gamma=f(H_{0}+p_{t})-f(-H_{1})=\Gamma_{0}+\Gamma_{1},$ (17)

where we note $\Gamma_{0}=f(H_{0}+p_{t})$ and $\Gamma_{1}=-f(-H_{1})$. We
remark that $\Gamma_{0}$ ($\Gamma_{1}$) is integrable because it is a function
of $H_{0}+p_{t}$ ($H_{1}$), which is integrable. Moreover, we have

$\displaystyle L_{\Gamma_{0}}$
$\displaystyle=f^{\prime}(H_{0}+p_{t})L_{H_{0}+p_{t}},$ $\displaystyle
L_{\Gamma_{1}}$ $\displaystyle=f^{\prime}(-H_{1})L_{H_{1}}.$ (18)

Because $H_{0}+p_{t}$ (resp. $H_{1}$) is a first integral of $\Gamma_{0}$
(resp. $\Gamma_{1}$), we have

$\displaystyle\exp(\sigma L_{\Gamma_{0}})$ $\displaystyle=\exp(\sigma
f^{\prime}(H_{0}+p_{t})L_{H_{0}+p_{t}})=\exp(\tau_{0}L_{H_{0}+p_{t}}),$
$\displaystyle\exp(\sigma L_{\Gamma_{1}})$ $\displaystyle=\exp(\sigma
f^{\prime}(-H_{1})L_{H_{1}})=\exp(\tau_{1}L_{H_{1}}),$ (19)

where

$\tau_{0}=\sigma f^{\prime}(H_{0}+p_{t})\text{ and }\tau_{1}=\sigma
f^{\prime}(-H_{1}).$ (20)

The operator $\exp(\sigma L_{\Gamma_{0}})$ ($\exp(\sigma L_{\Gamma_{1}})$) is
equivalent to the regular operator $\exp(\tau_{0}L_{H_{0}+p_{t}})$
($\exp(\tau_{1}L_{H_{1}})$) with a modified time step. The operator
$\exp(\sigma L_{\Gamma})$ can be approximated by a composition of operators
$\exp(\sigma L_{\Gamma_{k}})$

$\mathcal{S}_{\Gamma}(\sigma)=\prod_{i=1}^{N}\exp(a_{i}\sigma
L_{\Gamma_{0}})\exp(b_{i}\sigma L_{\Gamma_{1}}).$ (21)

With the BCH formula,
$\mathcal{S}_{\Gamma}(\sigma)=\exp(\sigma(L_{\Gamma}+L_{\Gamma_{\mathrm{err}}}),$
where $\Gamma_{\mathrm{err}}$ is an error Hamiltonian that depends on
$\sigma$. The symplectic map $\mathcal{S}_{\Gamma}(\sigma)$ exactly integrates
the modified Hamiltonian ${\Gamma+\Gamma_{\mathrm{err}}}$. The iteration of
$\mathcal{S}_{\Gamma}(\sigma)$ with fixed $\sigma$ is a symplectic integrator
algorithm for $\Gamma$.

When the timescale $\sigma$ is small enough, the numerical values of $H_{0}$
and $H_{1}$ do not change significantly between each step of the composition.
We have $\mathcal{S}_{\Gamma}(\sigma)\simeq\mathcal{S}(\tau)$ with
$\tau\simeq\sigma f^{\prime}(-H_{1})$. In other words, $\mathcal{S}_{\Gamma}$
behaves as $\mathcal{S}$ with an adaptive time-step while keeping the bounded
energy properties of a fixed time-step integrator.

###  3.2 Choice of the regularisation function

We wish the step sizes (20) to become smaller when planets experience close
encounters. These time-steps are determined by the derivative of $f$. For
nearly Keplerian systems, Mikkola & Tanikawa (1999) and Preto & Tremaine
(1999) studied renormalisation functions such that $f^{\prime}(x)\propto
x^{-\gamma}$, where $\gamma>0$ (this corresponds to power-law functions and
the important case of $f=\ln$).

However, these authors considered splitting of the type $H_{0}=T({\bf p})$ and
$H_{1}=U({\bf q})$. As pointed out in (Blanes & Iserles, 2012), when the
Hamiltonian is split as the Keplerian part plus an integrable perturbation, it
appears that both terms $K({\bf p},{\bf q})+p_{t}$ and $-H_{1}$ can change
signs, resulting in large errors in the integration.

We remark that the use of $f=\ln$ gives the best result when two planets
experience a close encounter, but it leads to large energy errors far away for
collision when $H_{1}$ is nearly 0. Based on these considerations, we require
$f$ to verify several properties to successfully regularise the perturbed
Keplerian problem in presence of close encounters:

  * •

$f^{\prime}$ should only depend on the magnitude of the perturbation $H_{1}$,
therefore we require $f^{\prime}$ even (and $f$ odd).

  * •

As already pointed out, $t$ must be an increasing function of $s,$ therefore
$f^{\prime}>0$

  * •

$f^{\prime}$ should be smooth, therefore we exclude piecewise renormalisation
functions.

  * •

We require that the regularisation vanishes in absence of perturbation _i.e_
$f^{\prime}(0)=1$. It results that for vanishing perturbation, $\sigma=\tau$.

  * •

Let $E_{1}$ be a typical value of $H_{1}$ far from close encounters that we
determine an expression for below. $f^{\prime}$ should only depend on
$H_{1}/E_{1}$ in order to only track relative changes in the magnitude of the
perturbation.

  * •

For high values of the perturbation (_e.g._ during close encounters), we
require $f^{\prime}(H_{1})\sim|E_{1}/H_{1}|$ such that we preserve the good
properties pointed out by previous studies (Preto & Tremaine, 1999; Mikkola &
Tanikawa, 1999).

  * •

To reduce the computational cost of the integrator, $f^{\prime}$ has to be
numerically inexpensive to evaluate.

Figure 1: Upper panel: $f^{\prime}$, Eq. (22) as a function of $h/E_{1}$. The
asymptotic value for $h/E_{1}\to+\infty$ is given in green.
Lower panel: $f/E_{1}$, Eq. (23) as a function of $h/E_{1}$.

These properties lead to a very natural choice for $f^{\prime}$ (see in figure
1). We choose

$f^{\prime}(h)=\frac{1}{\sqrt{1+\left(\frac{h}{E_{1}}\right)^{2}}}.$ (22)

Taking the odd primitive of Eq. (22), we find

$f(h)=E_{1}\operatorname{arcsinh}\left(\frac{h}{E_{1}}\right).$ (23)

For now on, $f$ always refers to definition (23). With this choice of the
function $f$, the Hamiltonian $\Gamma$ takes the form

$\Gamma=E_{1}\operatorname{arcsinh}\left(\frac{H_{0}+p_{t}}{E_{1}}\right)+E_{1}\operatorname{arcsinh}\left(\frac{H_{1}}{E_{1}}\right),$
(24)

where we used the oddity of $f$.

We need to define $E_{1}$ more explicitly. When planets are far from each
other, their mutual distance is of the same order as the typical distance
between the planets and the star. Using the same idea as in (Marchal & Bozis,
1982; Petit et al., 2018), we define a typical length unit of the system based
on the initial system energy $E_{0}$. We have

$a_{\mathrm{typical}}=-\frac{\mathcal{G}M^{*}}{2E_{0}},$ (25)

where $M^{*}=\sum_{0\leq i<j}m_{i}m_{j}$. The typical value for the
perturbation Hamiltonian far away from collision can be defined as

$E_{1}=\frac{\mathcal{G}m^{*}}{a_{\mathrm{typical}}}=\frac{2|E_{0}|m^{*}}{M^{*}},$
(26)

where $m^{*}=\sum_{1\leq i<j}m_{i}m_{j}$. We note that we have
$E_{1}/E_{0}=\mathrm{O}(\varepsilon)$.

The behaviour of the higher order derivative of $f$ is useful for the error
analysis, and in particular, their dependence in $\varepsilon$. The $k$-th
derivative of $f$ has for expression

$f^{(k)}(h)=E_{1}^{1-k}\operatorname{arcsinh}^{(k)}\left(\frac{h}{E_{1}}\right)=\mathrm{O}(\varepsilon^{1-k}).$
(27)

##  4 Order of the scheme

###  4.1 Analytical error estimates

As explained in section 2, most of the planet dynamics simulations are made
with a simple a second-order scheme such as the Wisdom-Holman leapfrog
integrator (Wisdom & Holman, 1991). It is indeed possible to take advantage of
the hierarchy between $H_{0}$ and $H_{1}$. This can be done by adding
symplectic correctors (Wisdom et al., 1996; Rein & Tamayo, 2015) or by
cancelling the term of the form $\varepsilon\tau^{k}$ up to a certain order
(Laskar & Robutel, 2001).

Hierarchical order schemes such as $\mathcal{SABA}(10,6,4)$ (Blanes et al.,
2013; Farrés et al., 2013) behave effectively as a tenth-order integrator if
$\varepsilon$ is small enough. Cancelling only selected terms reduces the
number of necessary steps of the scheme, which reduces the numerical error and
improves the performances.

Unfortunately, this property cannot be used for the regularised Hamiltonian
because $\Gamma_{0}$ and $\Gamma_{1}$ are almost equal in magnitude.
Nevertheless, it should be noted that the equations of motion (13) and the Lie
derivatives (18) keep their hierarchical structure. The Poisson bracket of
$\Gamma_{0}$ and $\Gamma_{1}$ gives

$\\{\Gamma_{0},\Gamma_{1}\\}=f^{\prime}(H_{0}+p_{t})f^{\prime}(H_{1})\\{H_{0},H_{1}\\}.$
(28)

Because $f^{\prime}$ does not depend directly on $\varepsilon$ (by choice,
$f^{\prime}$ only tracks relative variations of $H_{1}$),
$\\{\Gamma_{0},\Gamma_{1}\\}$ is of the order of $\varepsilon$. However, for
higher order terms in $\sigma$ in $H_{\mathrm{err}}$, it is not possible to
exploit the hierarchical structure in $\varepsilon$. We consider the terms of
the order of $\sigma^{2}$ in the error Hamiltonian for the integration of
$\Gamma$ using the leapfrog scheme. We have (_e.g._ Laskar & Robutel, 2001)

$\Gamma_{\mathrm{err}}=\frac{\sigma^{2}}{12}\\{\\{\Gamma_{0},\Gamma_{1}\\},\Gamma_{0}\\}+\frac{\sigma^{2}}{24}\\{\\{\Gamma_{0},\Gamma_{1}\\},\Gamma_{1}\\}+\mathrm{O}(\sigma^{4}).$
(29)

In order to show the dependence in $\varepsilon$, we develop the Poisson
brackets in Eq. (29)

$\displaystyle\Gamma_{\mathrm{err}}=$ $\displaystyle\ \ \
\frac{\sigma^{2}}{12}f^{\prime}(H_{0}+p_{t})^{2}f^{\prime}(H_{1})\\{\\{H_{0},H_{1}\\},H_{0}\\}$
$\displaystyle-\frac{\sigma^{2}}{12}f^{\prime}(H_{0}+p_{t})^{2}f^{\prime\prime}(H_{1})\left(\\{H_{0},H_{1}\\}\right)^{2}$
(30)
$\displaystyle+\frac{\sigma^{2}}{24}f^{\prime}(H_{1})^{2}f^{\prime}(H_{0}+p_{t})\\{\\{H_{0},H_{1}\\},H_{1}\\}$
$\displaystyle+\frac{\sigma^{2}}{24}f^{\prime}(H_{1})^{2}f^{\prime\prime}(H_{0}+p_{t})\left(\\{H_{0},H_{1}\\}\right)^{2}+\mathrm{O}(\sigma^{4}).$

The first and third terms only depend on $f^{\prime}$ and nested Poisson
brackets of $H_{0}$ and $H_{1}$. Their dependency on $\varepsilon$ is thus
determined by the Poisson brackets as in the fixed time-step case. The first
term is of the order of $\varepsilon\sigma^{2}$ and the third of the order of
$\varepsilon^{2}\sigma^{2}$. On the other hand, the second and last terms
introduce the second derivative of $f$ as well as a product of Poisson bracket
of $H_{0}$ and $H_{1}$. Equation (27) shows that they are of the order of
$\varepsilon\sigma^{2}$.

In order to cancel all terms of the order of $\varepsilon\sigma^{2}$, it is
necessary to cancel both terms in $\sigma^{2}$ in Eq. (29). Thus, the strategy
used in (Blanes et al., 2013) does not provide a scheme with a hierarchical
order because every Poisson bracket contributes with terms of the order of
$\varepsilon$ to the error Hamiltonian.

It is easy to extend the previous result to all orders in $\sigma$. We
consider a generic error term of the form

$\Gamma_{\mathrm{gen}}=\sigma^{n-1}\\{\\{\Gamma_{k_{0}},\Gamma_{k_{1}}\\},\dots,\Gamma_{k_{n}}\\},$
(31)

where $k_{j}$ is either 0 or 1 and $k_{0}=0$ and $k_{1}=1$. The development of
$\Gamma_{\mathrm{gen}}$ into Poisson brackets of $H_{0}$ and $H_{1}$ contains
a term that has for expression

$\sigma^{n-1}f^{(n_{0})}(H_{0}+p_{t})f^{(n_{1})}(H_{0}+p_{t})f_{1}^{\prime
n_{0}}f_{0}^{\prime n_{1}}(\\{H_{0},H_{1}\\})^{n-1},$ (32)

where $n_{j}$ is the number of $\Gamma_{j}$ in $\Gamma_{\mathrm{gen}}$, and
$f_{1}^{\prime}=f^{\prime}(H_{1})$. Because $n=n_{0}+n_{1}$, we deduce from
Eq. (27) that the term (32) is of the order of $\varepsilon\sigma^{n-1}$.
Thus, the effective order of the Hamiltonian is always
$\varepsilon\sigma^{r_{\min}}$ , where $r_{\min}$ is the smallest exponent
$r_{k}$ in Eq. (5). The error is still linear in $\varepsilon$. Hence, it is
still worth using the Keplerian splitting.

We use schemes that are not dependent on the hierarchy between $H_{0}$ and
$H_{1}$. McLachlan (1995a) provided an exhaustive list of the optimal methods
for fourth-, sixth-, and eighth-order integrators. Among the schemes he
presented, we select a sixth-order method consisting of a composition of
${n=7}$ leapfrog steps introduced by Yoshida (1990) and an eighth-order method
that is a combination of ${n=15}$ leapfrog steps. A complete review of these
schemes can be found in (Hairer et al., 2006). The coefficients of the schemes
are given in Appendix B.

In order to solve the Kepler step, we adopt the same approach as (Mikkola,
1997; Rein & Tamayo, 2015). The details on this particular solution as well as
other technical details are given in Appendix A.

###  4.2 Non-integrable perturbation Hamiltonian

When the classical splitting of the Hamiltonian written in canonical
heliocentric coordinates is used, $H_{1}$ depends on both positions and
momenta. We can write $H_{1}$ as a sum of two integrable Hamiltonian
${H_{1}=T_{1}+U_{1}}$, where $T_{1}$ is the indirect part that only depends on
the momenta and $U_{1}$ is the planet interaction potential, only depending on
positions (Farrés et al., 2013). We thus approximate the evolution operator
(19) by

$\exp(\sigma
L_{\Gamma_{1}})=\exp\left(\frac{\tau_{1}}{2}L_{T_{1}}\right)\exp\left(\tau_{1}L_{U_{1}}\right)\exp\left(\frac{\tau_{1}}{2}L_{T_{1}}\right)+\mathrm{O}(\varepsilon^{3}\tau_{1}^{3}),$
(33)

The numerical results suggest that heliocentric coordinates give slightly more
accurate results at constant cost. It is possible to use an alternative
splitting that is often called democratic heliocentric splitting (Laskar,
1990; Duncan et al., 1998). With this partition of the Hamiltonian, the
kinetic and the potential part of the perturbation Hamiltonian commute.
Therefore, the step $\exp(\sigma L_{\Gamma_{1}})$ is directly integrable using
the effective step size $\tau_{1}$ (20). In the following numerical tests, we
always use the classical definition when we refer to heliocentric coordinates.

##  5 Error analysis near the close encounter

###  5.1 Time-step and scheme comparison

Figure 2: Comparison between various schemes detailed in the body of the text
for the integration of a single synodic period for a system of equal-mass
planets initially on planar circular orbits with
$\varepsilon=2m/m_{0}=10^{-5}$ and $\alpha=0.8$. The solid lines represent the
integrations using a fixed real time-step and the dots the adaptive time-
steps. The closest planet approach is 0.19992 AU. Figure 3: Same as figure 2
for a system of two equal-mass planets initially on planar circular orbits
with ${\varepsilon=2m/m_{0}=10^{-5}}$ and $\alpha=0.97$. The closest planet
approach is ${\numprint{3.68E-5}}$ AU.

In this section, we test how a single close encounter affects the energy
conservation. To do so, we compare different integration schemes for a two-
planet system, initially on circular orbits, during an initial synodic period

$T_{\mathrm{syn}}=\frac{2\pi}{n_{1}-n_{2}}=\frac{2\pi}{n_{2}(\alpha^{-3/2}-1)},$
(34)

where $n_{i}$ is the mean motion of planet $i$ and $\alpha$ is the ratio of
the semi-major axis.

Because the time-step is renormalised, we need to introduce a cost function
that depends on the fictitious time $s$ as well as on the number of stages
involved in the scheme of the integrator. We define the cost of an integrator
as the number of evaluations of $\exp(a\sigma L_{\Gamma_{0}})\exp(b\sigma
L_{\Gamma_{1}})$ that are required to integrate for a given real time period
$T_{\mathrm{syn}}$

$\mathcal{C}=\frac{s_{\mathrm{syn}}n}{T_{\mathrm{syn}}\sigma},$ (35)

where $s_{\mathrm{syn}}$ is the fictitious time after $T_{\mathrm{syn}}$,
$\sigma$ is the fixed fictitious time-step, and $n$ is the number of stages of
the integrator. We also compare the renormalised integrators to the same
scheme with fixed time-step. For a fixed time-step, the cost function is
simply given by $\mathcal{C}_{\mathrm{fixed}}=n/\tau$, where $\tau$ is the
time-step (Farrés et al., 2013). We present different configurations in
figures 2, 3, and 4.

In the two first sets, we integrate the motion of two equal-mass planets on
circular orbits, starting in opposition with respect to the star. In both
simulations, we have $\varepsilon=(m_{1}+m_{2})/m_{0}=10^{-5}$, the stellar
mass is $1\ {\rm M}_{\odot}$ , and the outer planet semi-major axis is 1 AU.
In both figures, we represent the averaged relative energy variation

$\frac{\Gamma}{E_{0}}\simeq f^{\prime}(H_{1})\frac{\Delta E}{E_{0}}$ (36)

as a function of the inverse cost $\mathcal{C}^{-1}$ of the integration for
various schemes:

  * •

the classical order 2 leapfrog $\mathcal{ABA}(2,2)$,

  * •

the scheme $\mathcal{ABA}(8,2)$ from (Laskar & Robutel, 2001),

  * •

the scheme $\mathcal{ABA}(10,6,4)$ from (Blanes et al., 2013)

  * •

the sixth- and eighth-order schemes from (McLachlan, 1995a) that we introduced
in the previous section.

For each scheme, we plot with a solid line the result of the fixed time-step
algorithm and with dots the results of the adaptive time-step integrator. In
the results presented in figure 2 and 3, the systems are integrated in Jacobi
coordinates. In figure 4 we integrate in heliocentric coordinates. For the
three cases, the integrations were also carried out with the other set of
coordinates with almost no differences.

In Figure 2 the initial semi-major axis ratio $\alpha$ is 0.8 and the closest
approach of the two planets is 0.19992 AU. In this case, a fixed time-step
algorithm provides accurate results, and when a generalized order scheme such
as $\mathcal{ABA}(8,2)$ or $\mathcal{ABA}(10,6,4)$ is used, the performances
are better with a fixed time-step than with an adaptive one. However, we see
no sensitive differences in accuracy between the fixed and adaptive versions
of the $\mathcal{ABA}(6*)$ and $\mathcal{ABA}(8*)$ schemes. For the most
efficient schemes, machine precision is reached for an inverse cost of a few
$10^{-3}$. This case illustrates the point made in section 4 that there is no
advantage in using the hierarchical structure of the original equation in the
choice of the scheme.

In Figure 3 the initial semi-major axis ratio $\alpha$ is 0.97 and the closest
approach of the two planets is $\numprint{3.68E-5}$ AU. The radius of a planet
of mass 5×10-6 $\mathrm{M}_{\odot}$ with a density of ${5\
\mathrm{g.cm^{-3}}}$ , denser than Earth, is 5.21×10-4 AU, that is, ten times
greater. We can compute the density that corresponds to a body of mass 5×10-6
$\mathrm{M}_{\odot}$ and radius half the closest approach distance. We obtain
a density of 1.1×105 $\mathrm{g.cm^{-3}}$ , which is close to the white dwarf
density. In this case, the fixed time-step algorithm is widely inaccurate
while the performance of the adaptive time-step algorithm remains largely
unchanged. This demonstrates how powerful this new approach is for the
integration of systems with a few planets.

Figure 4: Same as figure 2 for a system of two equal-mass planets initially on
planar circular orbits with ${\varepsilon=2m/m_{0}=10^{-3}}$ and $\alpha=0.9$.
The closest planet approach is ${\numprint{1.7e-2}}$ AU. This system is
integrated in heliocentric coordinates to demonstrate that we obtain similar
results to the case of Jacobi coordinates.

We also give an example with higher mass planets in figure 4. We plot the
relative energy error as a function of the inverse cost for a system of two
planets of mass ${\numprint{5e-4}\leavevmode\nobreak\ \mathrm{M}_{\odot}}$
that correspond to the case of $\varepsilon=10^{-3}$. The initial semi-major
axis ratio is 0.9 to ensure a closest-approach distance of 1.7×10-2 AU.

For this particular close encounter, it is still possible to reach machine
precision with a fixed time-step at the price of a very small time-step. On
the other hand, the close encounter is perfectly resolved for the adaptive
time-step, even if $\varepsilon$ is larger by two orders of magnitude. Because
of the higher masses, it is necessary to take smaller time-steps to ensure
that the relative energy error remains at machine precision.

###  5.2 Behaviour at exact collision

Figure 5: Mutual distance of the planets as a function of the fictitious time
for 1000 initial conditions chosen as described in the body of the text. On
the upper axis, we give the real time for the initial condition that reached
the closest approach.

To demonstrate the power of our algorithm in resolving close encounters, we
considered the case of two planets that almost exactly collide (as if the
bodies were material points). We take two planets of mass $m=10^{-5}\
\mathrm{M}_{\odot}$, the first on an orbit of eccentricity 0.1 and semi-major
axis 0.95 AU, and the second on a circular orbit of semi-major axis 1 AU. We
then fine-tuned the mean longitudes of the two planets to ensure that we
approach the exact collision. We integrated 1000 initial conditions for which
we linearly varied $\lambda_{1}$ between -0.589 and -0.587, and we took
${\lambda_{2}=-1.582}$.

In figure 5 we represent for each initial condition the mutual distance of the
planets as a function of the fictitious time $s$. In the upper axis, we give
the real time of the particular initial condition that gave the closest
approach. In this simulation, all initial conditions reached a mutual distance
smaller than 1.1×10-4 AU and the closest approach is 5.5×10-9 AU. For the
chosen masses, we can compute the physical radius of such bodies assuming a
particular bulk density. We have

$R_{p}=\numprint{5.21e-2}\left(\frac{\rho}{1\mathrm{g.cm^{-3}}}\right)^{-1/3}\left(\frac{m_{p}}{1\mathrm{M_{\odot}}}\right)^{1/3},$
(37)

where $\rho$ is the bulk density of the planet. The radius is only a weak
function of the density over the range of planet bulk densities. When we
assume a very conservative value of $\rho=6\ \mathrm{g.cm^{-3}}$, we obtain a
planet radius of 6.17×10-4 AU, that is, well above the close encounters
considered here.

Figure 6: Relative energy error as a function of the closest mutual approach
reached during a close encounter of super-Earths in a two-planet system and in
a six-planet system. The particular initial conditions are described in the
text. In the upper part, we give the mutual distance in units of planet radii
for a planet of density $\rho=6\ \mathrm{g.cm^{-3}}$.

We plot the energy error as a function of the closest approach distance in
figure 6. We also considered a more practical case where four other planets
were added to the system. The additional planets were located on circular
orbits, had the same masses of $10^{-5}\ \mathrm{M}_{\odot}$ and were situated
outside of the orbit of the second planet, with an equal spacing of 0.15 AU.
Their initial angles are randomly drawn but were similar in every run.

In both set-ups, we used the heliocentric $\mathcal{ABA}(8*)$ scheme, and the
fictitious time-step was $\sigma=10^{-2}$. In the case of the two-planet
system, machine precision is reached even for an approach of approximately
$10^{-6}$ AU, which corresponds to a planet radius of about $10^{-3}$ . When
planets are added, it is more difficult to reach very close encounters because
the the other planets are perturbed. Nonetheless, approaches up to 6×10-7 AU
are reached, and the close encounter is resolved at machine precision up to a
closest distance of $\numprint{3e-5}$ AU or a planet radius of $5\times
10^{-2}$ .

##  6 Pericenter regularisation

When planet-planet scattering is studied, the closest distance between the
planet and the star may be reduced. As a result, the fixed time-step becomes
too large and the passage at pericenter is insufficiently resolved. In order
to address this problem, we can combine the close-encounter regularisation of
section 3 with a regularisation method introduced by (Mikkola, 1997).

We first detail Mikkola’s idea. When integrating a few-body problem
Hamiltonian, a time renormalisation $\mathrm{d}t=g({\bf q})\mathrm{d}u$ can be
introduced, where

$g({\bf q})=\left(\sum_{i=1}^{N}\frac{A_{i}}{q_{i}}\right)^{-1}.$ (38)

Here $A_{i}$ are coefficients that can be chosen arbitrarily. In this article,
we define them as

$A_{i}=\frac{m_{i}a_{\mathrm{typical}}}{\sum_{j=1}^{N}m_{j}},$ (39)

where $a_{\mathrm{typical}}$, defined in Eq. (25), is the typical length scale
of the system. The new Hamiltonian $\Upsilon$ in the extended phase space
$({\bf q},t,{\bf p},p_{t})$ is

$\Upsilon=\Upsilon_{0}+\Upsilon_{1}=g({\bf q})(H_{0}+p_{t})+g({\bf q})H_{1}.$
(40)

As in section 3, in the sub-manifold $\\{p_{t}=H(0)\\}$, $\Upsilon$ and $H$
have the same equations of motion up to the time transformation. If $H_{1}$
only depends on ${\bf q}$ (using Jacobi coordinates, e.g.), $\Upsilon_{1}$ is
trivial to integrate.

###  6.1 Kepler step

It is also possible to integrate $\Upsilon_{0}$ as a modified Kepler motion
through the expression of $g$ (Mikkola, 1997). We denote
$\upsilon_{0}=\Upsilon_{0}(0)$ and $\tilde{H}_{0}=g^{-1}({\bf
q})(\Upsilon_{0}-\upsilon_{0})$. $\tilde{H}_{0}$ has the same equations of
motion as $\Upsilon_{0}$ up to a time transformation
$\mathrm{d}\tilde{t}=g^{-1}({\bf q})\mathrm{d}u$. We have

$\displaystyle\tilde{H}_{0}$
$\displaystyle=H_{0}+p_{t}-\upsilon_{0}\sum_{i=1}^{N}\frac{A_{i}}{q_{i}}$
$\displaystyle=p_{t}+\sum_{i=1}^{N}\frac{\|{\bf
p}_{i}\|^{2}}{2m_{i}}-\frac{\mu_{i}m_{i}+\upsilon_{0}A_{i}}{q_{i}}=p_{t}+\sum_{i=1}^{N}\tilde{K}_{i},$
(41)

where $\tilde{K}_{i}$ is the Hamiltonian of a Keplerian motion of planet $i$
with a modified central mass

$\tilde{\mu}_{i}=\mu_{i}\left(1+\frac{\upsilon_{0}A_{i}}{\mu_{i}m_{i}}\right)=\mu_{i}\left(1+\frac{\upsilon_{0}a_{\mathrm{typical}}}{\mu_{i}\sum_{j=1}^{N}m_{j}}\right).$
(42)

However, the time equation must be solved as well. We integrate with a fixed
fictitious time-step $\Delta u$. The time $\Delta\tilde{t}(u)$ is related to
$u$ through the relation

$\Delta u=\int_{0}^{\Delta\tilde{t}}g^{-1}({\bf
q}(\tilde{t}))\mathrm{d}\tilde{t}=\sum_{i=1}^{N}A_{i}\int_{0}^{\Delta\tilde{t}}\frac{1}{r_{i}(\tilde{t}))}\mathrm{d}\tilde{t},$
(43)

where $r_{i}(\tilde{t})$ follows a Keplerian motion. We can rewrite equation
(43) with the Stumpff formulation of the Kepler equation (Mikkola, 1997; Rein
& Tamayo, 2015, see Appendix A.1). Because
$\int_{0}^{\Delta\tilde{t}}\frac{1}{r_{i}(\tilde{t}))}\mathrm{d}\tilde{t}=X_{i}$,
we have

$\Delta u=\sum_{i=1}^{N}A_{i}X_{i}.$ (44)

As a consequence, the $N$ Stumpff-Kepler equations

$\Delta\tilde{t}=r_{0i}X_{i}+\eta_{0i}G_{2}(\beta_{0i},X_{i})+\zeta_{0i}G_{3}(\beta_{0i},X_{i})=\kappa_{i}(X_{i})$
(45)

must be solved simultaneously with equation (44). To do so, we used a
multidimensional Newton-Raphson method on the system of $N+1$ equations
consisting of the $N$ Kepler equations (45) and equation (44) of unknowns
${Y=(X_{1},\dots,X_{N},\Delta\tilde{t})}$. The algorithm is almost as
efficient as the fixed-time Kepler evolution because it does not add up
computation of Stumpff’s series. At step $k$, we can obtain $Y^{(k+1)}$
through the equation

$Y^{(k+1)}=Y^{(k)}-\mathrm{d}F^{-1}(Y^{(k)})(F(Y^{(k)})),$ (46)

where

$F=\left(\begin{array}[]{c}\kappa_{1}(X_{1})-\Delta\tilde{t}\\\ \vdots\\\
\kappa_{N}(X_{N})-\Delta\tilde{t}\\\ \sum_{i=1}^{N}A_{i}X_{i}-\Delta
u\end{array}\right)$ (47)

and

$\mathrm{d}F=\left(\begin{array}[]{cccc}\kappa^{\prime}_{1}(X_{1})&0&\cdots&-1\\\
0&\ddots&0&-1\\\ \cdots&0&\kappa^{\prime}_{N}(X_{N})&-1\\\
A_{1}&\cdots&A_{N}&0\end{array}\right),$ (48)

with $\kappa_{i}^{\prime}$ being the derivative with respect to $X_{i}$ of
$\kappa_{i}$ (Eq. 45).

Equation (46) can be rewritten as a two-step process where a new estimate for
the time $\Delta\tilde{t}^{(k+1)}$ is computed and is then used to estimate
$X^{(k+1)}_{i}$. We have

$\Delta\tilde{t}^{(k+1)}=\left.\left(\Delta
u+\sum_{i=1}^{N}\frac{A_{i}\left(\kappa_{i}(X_{i}^{(k)})-\kappa^{\prime}_{i}(X_{i}^{(k)})X_{i}^{(k)}\right)}{\kappa^{\prime}_{i}(X_{i}^{(k)})}\right)\
\middle/\
\left(\sum_{i=1}^{N}\frac{A_{i}}{\kappa^{\prime}_{i}}\right)\right.,$ (49)

and

$X^{(k+1)}_{i}=\frac{\Delta\tilde{t}^{(k+1)}+\kappa^{\prime}_{i}(X_{i}^{(k)})X_{i}^{(k)}-\kappa_{i}(X_{i}^{(k)})}{\kappa^{\prime}_{i}(X_{i}^{(k)})}.$
(50)

###  6.2 Heliocentric coordinates

In heliocentric coordinates, $H_{1}^{H}$ depends on ${\bf p}$ as well. As a
result, $gH_{1}^{H}$ cannot be easily integrated, and it is not even possible
write it as a sum of integrable Hamiltonians. To circumvent this problem, we
can split $gH_{1}^{H}$ into $gU_{1}^{H}$ and $gT_{1}^{H}$. The potential part
$gU_{1}^{H}$ is integrable, but a priori, $gT_{1}^{H}$ is not integrable.

The integration of $gT_{1}^{H}$ can be integrated using a logarithmic method
as proposed in (Blanes & Iserles, 2012). The evolution of $gT_{1}^{H}$ during
a step $\Delta t$ is the same as the evolution of $\log gT_{1}^{H}=\log g+\log
T_{1}^{H}$, which is separable, for a step $\mathcal{T}_{1}\Delta t$ where
$\mathcal{T}_{1}=gT_{1}^{H}|_{t=0}$. Therefore, we can approximate $\log
gT_{1}^{H}$ using a leapfrog scheme, and the error is of order $\Delta
u^{2}\varepsilon^{3}$ as well.

Then, we can approximate the heliocentric step using the same method as in
Sect. 4.2. In this case, it is necessary to approximate the step even when
democratic heliocentric coordinates are used because $g({\bf q})$ does not
commute with $T_{1}^{H}$.

###  6.3 Combining the two regularisations

We showed that the Hamiltonian $\Upsilon$ is separable into two parts that are
integrable (or nearly integrable for heliocentric coordinates). Therefore we
can simply regularise the close encounters by integrating the Hamiltonian

$\tilde{\Gamma}=f(\Upsilon_{0}+p_{u})+f(\Upsilon_{1})=f(g({\bf
q})(H_{0}+p_{t})+p_{u})+f(g({\bf q})H_{1}),$ (51)

where $p_{u}$ is the momentum associated with the intermediate time $u$ used
to integrate $\Upsilon$ alone. The time equation is

$\frac{\mathrm{d}t}{\mathrm{d}s}=\frac{\partial\tilde{\Gamma}}{\partial
p_{t}}=f^{\prime}(\Upsilon_{0})g({\bf q}).$ (52)

We need to place ourselves on the sub-manifold such that
$\Upsilon_{0}+p_{u}+\Upsilon_{1}=0$ and $H_{0}+p_{t}+H_{1}=0$, in order to
have the same equations of motion for $\tilde{\Gamma}$, $\Upsilon$ and $H$ .
Both of these conditions are fulfilled by choosing $p_{t}=-E_{0}$ and
$p_{u}=0$.

##  7 Long-term integration performance

So far, we only presented the performance of the algorithm for very short
integrations. In this section, we present the long-term behaviour of the
integrator for systems with a very chaotic nature. We consider two different
configurations, a system composed of equal planet masses on initially
circular, coplanar and equally spaced orbits, used as a test model for
stability analysis since the work of Chambers et al. (1996). The second is a
similar system but with initially moderate eccentricities and inclinations.

###  7.1 Initially circular and coplanar systems

Figure 7: Evolution of the relative energy error for 100 systems composed of
six planets on initially circular orbits. The light curves represent
individual systems, and the thicker blue curves show the average over the 100
systems. The orange straight line is proportional to $\sqrt{t}$ in order to
show that the integrator follows Brouwer’s law. The orange line is the same in
all energy error figures for comparison.

We integrated 100 systems of six initially coplanar and circular planets. The
planet masses were taken to be equal to $10^{-5}\ \mathrm{M}_{\oplus}$, and
the semi-major axis of the outermost planet was fixed to 1 AU, while the semi-
major axis ratios of the adjacent planets were all equal to 0.88. These valus
were chosen to ensure that the lifetime of the system was of the order of 300
kyr before the first collision. The fixed fictitious time-step was
$\sigma=10^{-2}$ yr. Because of the renormalisation, this approximately
corresponds to a fixed time-step of 6.3×10-3 yr in terms of computational
cost. The initial period of the inner planet was of the order of 0.38 yr, that
is, we have about 50 steps per orbit.

The simulations were stopped when two planet centres approached each other by
less than half the planet radii, assuming a density of 6 $\mathrm{g.cm^{-3}}$.
This stopping criterion is voluntarily nonphysical as it allows for a longer
chaotic phase that leads to more close encounters. We also monitored any
encounters with an approach closer than 2 Hill radii at 1 AU (0.054 AU) and
recorded its time, the planets involved, and the closest distance between the
two planets.

For the majority of the integrations, we observe moderate semi-major axis
diffusion without close encounters. About 1 kyr before the final collision,
the system enters a true scattering phase with numerous close encounters. The
integrations lasted on average 353 kyr, the shortest was 129 kyr long and the
longest 824 kyr. On average, we recorded 557 close encounters, and 68%
occurred during the last 10000 years. On average, 6.4 approaches below 0.01 AU
per system were recorded. Of these very close encounters, 95% occurred during
the last 1000 years.

The relative energy error evolution is shown in figure 7. We observe that the
integrator follows the Brouwer (1937) law because the energy error behaves as
a random walk (the error grows as $\sqrt{t}$). The scattering phase does not
affect the energy conservation: no spikes of higher error values occur towards
the end of the integrations.

###  7.2 Planet scattering on inclined and eccentric orbits

Figure 8: Evolution of the semi-major axis, periapsis, and apoapsis for a typical initial condition. The bold curve is the semi-major axis, and the filled region represents the extent of the orbit. The mutual inclinations are not represented. Table 1: Summary of the main statistics and specifications of the spatial and eccentric runs Run | J1 | J2 | H1 | H2
---|---|---|---|---
Coordinates | Jacobi | Jacobi | Heliocentric | Heliocentric
Pericenter | No | Yes | No | Yes
Equivalent fixed time-step | $\numprint{2.6e-3}\ \mathrm{yr}^{-1}$ | $\numprint{2.6e-3}\ \mathrm{yr}^{-1}$ | $\numprint{2.9e-3}\ \mathrm{yr}^{-1}$ | $\numprint{2.9e-3}\ \mathrm{yr}^{-1}$
Time-step per shortest orbit | 112 | 112 | 102 | 102
Average lifetime | 53.3 kyr | 53.0 kyr | 52.6 kyr | 46.3 kyr
Average number of close encounters | $\numprint{1.26e4}$ | $\numprint{1.17e4}$ | $\numprint{1.20e4}$ | $\numprint{1.07e4}$
Average number of very close encounters | 458 | 415 | 426 | 384

(a) J1: Jacobi coordinates, no pericenter regularisation. (b) J2: Jacobi
coordinates, pericenter regularisation.

(c) H1: heliocentric coordinates, no pericenter regularisation. (d) H2:
heliocentric coordinates, pericenter regularisation.

Figure 9: Evolution of the energy error for the four runs detailed in the text
and Table 1. The light curves represent individual systems, and the thicker
blue curves show the medians over the 100 systems. The orange filled lines are
proportional to $\sqrt{t}$ and are the same in all energy error figures for
comparison. The dashed orange lines represents a linear error in time to show
the deviation from the Brouwer law. We discuss the outliers in the main text.

Our second long-term test problem was a set of planetary systems with moderate
eccentricities and inclinations. The goal is here to test the integrator in a
regime where strong scattering occurs on long timescales with a much higher
frequency in very close encounters (shorter than 0.01 AU). In addition, this
case helps to demonstrate the advantage of the pericenter regularisation. The
initial conditions were chosen as follows:

  * •

Similarly to the previous case, we considered systems of six planets with an
equal mass of $10^{-5}\ \mathrm{M}_{\odot}$.

  * •

The semi-major axis ratio of adjacent pairs was 0.85.

  * •

The Cartesian components of the eccentricities $(e\cos(\varpi),e\sin(\varpi))$
and inclinations $(i\cos(\Omega),i\sin(\Omega))$ were drawn from a centred
Gaussian distribution with a standard deviation of 0.08. The average
eccentricities and inclinations were about 0.11.

  * •

The mean longitudes were chosen randomly.

We also kept the same stopping criterion and recorded the close encounters
(smaller than 2 Hill radii).

The integrations were performed with two sets of coordinates, Jacobi and
heliocentric; with and without pericenter regularisation. The main statistics
from the runs are summarised in table 1. For all integrations, we used the
eighth-order scheme and a fictitious time-step $\sigma=4\times 10^{-3}$. The
systems on average experience a close encounter every 5 yr. A typical
evolution is presented in figure 8. In this particular run, we see numerous
planet exchanges, almost 35000 close encounters occur, and the system
experiences more than 1200 very close encounters. Nevertheless, the final
relative energy error is $\numprint{6.8e-15}$.

We plot in Figs. 9(a), 9(b), 9(c), and 9(d) the relative energy errors of runs
J1, J2, H1, and H2, respectively. We plot in light colours the relative energy
error of individual systems and in thicker blue the median of the relative
energy error. The choice to use the median was made because some initial
conditions lead to a worse energy conservation, which makes the average less
informative. For the longer times, the median is no longer reliable because
the majority of the simulations had already stopped.

Despite the extreme scattering that occurs, the energy is conserved in most
systems within the numerical round-off error prescription. We also remark that
the heliocentric coordinates seem to provide a more stable integrator even
though the interaction step is approximated. In the two cases without the
pericenter regularisation, J1 and H1, it seems that the energy error is no
longer proportional to $\sqrt{t}$ but linear in $t$.

As explained above, because of the high eccentricities, the innermost
pericenter passage is not very well resolved, and therefore the step-to-step
energy variation is no longer close to machine precision. This effect has been
confirmed by the detailed analysis of the time-series close to energy spikes
by restarts of the integration from a binary file. We observed that the energy
did not change during the close encounter, but shortly after, when the inner
planet passed closer to the star because of the scattering.

We also observe an initial condition well above the machine-precision level in
all four runs. In this particular system, we have $e_{5}(t=0)=0.81,$ which
gives an initial pericenter at 0.15 AU. In this extreme case, the pericenter
regularisation leads to an improvement by two orders of magnitude of the error
in Jacobi coordinates. The heliocentric coordinates remain more efficient,
however, and we do not see an improvement in this particular case.

In all runs, the time distribution of the close encounters appears to be
rather uniform. Moreover, the probability distribution function (PDF) of the
closest approach during a close encounter is linear in the distance, as shown
in figure 10. This proves that the systems are in a scattering regime where
the impact parameters for the close encounters are largely random (see also
Laskar et al., 2011, figure 5).

Figure 10: PDF of the closest approach during close encounters. The PDFs of
the four runs are almost identical.

##  8 Comparison with existing integrators

Figure 11: Comparison of the different adaptive integrator runs presented in
Table 1 to other codes: IAS15, MERCURIUS, and SYMBA. The same initial
conditions were integrated. We plot the median of the energy error.

We ran the same eccentric and inclined initial conditions using the REBOUND
(Rein & Liu, 2012) package integrators IAS15 (Rein & Spiegel, 2015) and
MERCURIUS (Rein et al., 2019). IAS15 is an adaptive high-order Gauss–Radau
integrator, and MERCURIUS is a hybrid symplectic integrator similar to MERCURY
(Chambers, 1999). This choice is motivated because the REBOUND team made
recent and precise comparisons with other available integrators (Rein &
Spiegel, 2015; Rein & Tamayo, 2015). Moreover, we used a very similar
implementation of the Kepler step. We also ran the same tests with the multi-
time-step symplectic integrator SYMBA (Duncan et al., 1998).

We plot the energy error median that results from the integration. Our
integrator performs similarly to the non-symplectic integrator IAS15 and
better than MERCURIUS and SYMBA. Our integrator runs at the same speed as
MERCURIUS and SYMBA with a time-step of $\Delta t=10^{-4}$. However, IAS15 is
much faster on this particular problem (by a factor 10 on average). MERCURIUS
and SYMBA are not designed to work efficiently at such low-energy errors.
However, our motivation for this work was to find a precise symplectic
integrator for close encounters. We thus tried to obtain the best precision
from MERCURIUS and SYMBA to make the comparison relevant.

##  9 Discussion

We showed that time renormalisation can lead to very good results for the
integration of highly unstable planetary systems using a symplectic integrator
(section 3). The algorithm can use a scheme of arbitrary order, however, we
have not been able to take advantage of the hierarchical structure beyond the
first order in $\varepsilon$. This is due to the non linearity of the time
renormalisation as explained in section 4. It is possible, however, to cancel
every term up to arbitrary order in the time-step using non-hierarchical
schemes.

Our time renormalisation uses the perturbation energy to monitor when two
planets encounter each other. The algorithm is thus efficient if the two-
planet interaction contributes to a significant part of the perturbation
energy. As a result, this integrator is particularly well adapted for systems
of a few planets (up to a few tens) with similar masses (within an order of
magnitude). Moreover, it behaves extremely well in the two-planet problem even
when one planet mass is significantly lower than the other.

However, the algorithm is not yet able to spot close encounters between
terrestrial planets in a system that contains giant planets, such as the solar
system. In the case of the solar system, the perturbation energy is dominated
by the interaction between Jupiter and Saturn. We plan to address this
particular problem in the future.

## References

  * Blanes et al. (2013) Blanes, S., Casas, F., Farrés, A., et al. 2013, Applied Numerical Mathematics, 68, 58
  * Blanes & Iserles (2012) Blanes, S. & Iserles, A. 2012, Celestial Mechanics and Dynamical Astronomy, 114, 297
  * Brouwer (1937) Brouwer, D. 1937, The Astronomical Journal, 46, 149
  * Chambers et al. (1996) Chambers, J., Wetherill, G., & Boss, A. 1996, Icarus, 119, 261
  * Chambers (1999) Chambers, J. E. 1999, Monthly Notices of the Royal Astronomical Society, 304, 793
  * Dehnen & Hernandez (2017) Dehnen, W. & Hernandez, D. M. 2017, Monthly Notices of the Royal Astronomical Society, 465, 1201
  * Duncan et al. (1998) Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, The Astronomical Journal, 116, 2067
  * Farrés et al. (2013) Farrés, A., Laskar, J., Blanes, S., et al. 2013, Celestial Mechanics and Dynamical Astronomy, 116, 141
  * Gladman et al. (1991) Gladman, B., Duncan, M., & Candy, J. 1991, Celestial Mechanics and Dynamical Astronomy, 52, 221
  * Hairer et al. (2006) Hairer, E., Lubich, C., & Wanner, G. 2006, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Springer Science & Business Media)
  * Hernandez (2016) Hernandez, D. M. 2016, Monthly Notices of the Royal Astronomical Society, 458, 4285
  * Hernandez & Bertschinger (2015) Hernandez, D. M. & Bertschinger, E. 2015, Monthly Notices of the Royal Astronomical Society, 452, 1934
  * Kahan (1965) Kahan, W. 1965, Commun. ACM, 8, 40
  * Kinoshita et al. (1991) Kinoshita, H., Yoshida, H., & Nakai, H. 1991, Celestial Mechanics and Dynamical Astronomy, 50, 59
  * Koseleff (1993) Koseleff, P. V. 1993, in Applied Algebra, Algebraic Algorithms and Error-Correcting Codes, ed. G. Cohen, T. Mora, & O. Moreno, Lecture Notes in Computer Science (Springer Berlin Heidelberg), 213–230
  * Laskar (1990) Laskar, J. 1990, in Les Méthodes Modernes de La Mécanique Céleste. Modern Methods in Celestial Mechanics, 89–107
  * Laskar & Gastineau (2009) Laskar, J. & Gastineau, M. 2009, Nature, 459, 817
  * Laskar et al. (2011) Laskar, J., Gastineau, M., Delisle, J.-B., Farrés, A., & Fienga, A. 2011, Astronomy and Astrophysics, 532, L4
  * Laskar & Robutel (2001) Laskar, J. & Robutel, P. 2001, Celestial Mechanics and Dynamical Astronomy, 80, 39
  * Marchal & Bozis (1982) Marchal, C. & Bozis, G. 1982, Celestial Mechanics, 26, 311
  * McLachlan (1995a) McLachlan, R. 1995a, SIAM Journal on Scientific Computing, 16, 151
  * McLachlan (1995b) McLachlan, R. I. 1995b, BIT Numerical Mathematics, 35, 258
  * Mikkola (1997) Mikkola, S. 1997, Celestial Mechanics and Dynamical Astronomy, 67, 145
  * Mikkola (2008) Mikkola, S. 2008, in Dynamical Evolution of Dense Stellar Systems, Vol. 246, 218–227
  * Mikkola & Innanen (1999) Mikkola, S. & Innanen, K. 1999, Celestial Mechanics and Dynamical Astronomy, 74, 59
  * Mikkola & Tanikawa (1999) Mikkola, S. & Tanikawa, K. 1999, Celestial Mechanics and Dynamical Astronomy, 74, 287
  * Petit et al. (2018) Petit, A. C., Laskar, J., & Boué, G. 2018, Astronomy & Astrophysics, 617, A93
  * Preto & Tremaine (1999) Preto, M. & Tremaine, S. 1999, The Astronomical Journal, 118, 2532
  * Rein et al. (2019) Rein, H., Hernandez, D. M., Tamayo, D., et al. 2019, arXiv:1903.04972 [astro-ph]
  * Rein & Liu (2012) Rein, H. & Liu, S.-F. 2012, Astronomy and Astrophysics, 537, A128
  * Rein & Spiegel (2015) Rein, H. & Spiegel, D. S. 2015, Monthly Notices of the Royal Astronomical Society, 446, 1424
  * Rein & Tamayo (2015) Rein, H. & Tamayo, D. 2015, Monthly Notices of the Royal Astronomical Society, 452, 376
  * Stiefel & Scheifele (1971) Stiefel, E. L. & Scheifele, G. 1971, Linear and Regular Celestial Mechanics. (Springer, Berlin)
  * Stumpff (1962) Stumpff, K. 1962, Himmelsmechanik, Band I (VEB Deutscher Verlag der Wissenschaften, Berlin)
  * Wisdom (2006) Wisdom, J. 2006, The Astronomical Journal, 131, 2294
  * Wisdom & Holman (1991) Wisdom, J. & Holman, M. 1991, The Astronomical Journal, 102, 1528
  * Wisdom et al. (1996) Wisdom, J., Holman, M., & Touma, J. 1996, Fields Institute Communications, 10, 217
  * Yoshida (1990) Yoshida, H. 1990, Physics Letters A, 150, 262

Table B.1: Coefficients of the schemes used in this article. The values are computed from McLachlan (1995a). Scheme | Order | Stages | $a_{i}$ |  $b_{i}$
---|---|---|---|---
$\mathcal{ABA}(6*)$ | 6 | 7 | $a_{1}=0.39225680523877863191$ |  $b_{1}=0.78451361047755726382$
|  |  | $a_{2}=0.51004341191845769875$ | $b_{2}=0.23557321335935813368$
|  |  | $a_{3}=-0.471053385409756436635$ | $b_{3}=-1.17767998417887100695$
|  |  | $a_{4}=0.068753168252520105975$ | $b_{4}=1.3151863206839112189$
$\mathcal{ABA}(8*)$ | 8 | 15 | $a_{1}=0.370835182175306476725$ |  $b_{1}=0.74167036435061295345$
|  |  | $a_{2}=0.166284769275290679725$ | $b_{2}=-0.409100825800031594$
|  |  | $a_{3}=-0.109173057751896607025$ | $b_{3}=0.19075471029623837995$
|  |  | $a_{4}=-0.191553880409921943355$ | $b_{4}=-0.57386247111608226666$
|  |  | $a_{5}=-0.13739914490621317141$ | $b_{5}=0.29906418130365592384$
|  |  | $a_{6}=0.31684454977447705381$ | $b_{6}=0.33462491824529818378$
|  |  | $a_{7}=0.324959005321032390205$ | $b_{7}=0.31529309239676659663$
|  |  | $a_{8}=-0.240797423478074878675$ |  $b_{8}=-0.79688793935291635398$

##  Appendix A Implementation

We give technical details on our implementation choices in this appendix.

###  A.1 Kepler equation

The key step in any Wisdom-Holman algorithm is the numerical resolution of the
Kepler problem,

$H_{\mathrm{Kepler}}=\frac{\mathbf{v}^{2}}{2}-\frac{\mu}{r},$ (53)

where $\mu=\mathcal{G}M$ and $M$ is the central mass in the set of coordinates
used in the integration. Because this is the most expensive step from a
computational point of view, it is particularly important to optimise it.
Here, we closely followed the works by Mikkola (1997) and Mikkola & Innanen
(1999) and refer to them for more details. Rein & Tamayo (2015) presented an
unbiased numerical implementation that can be found in the package
REBOUND222https://rebound.readthedocs.io/en/latest/.

We consider a planet with initial conditions ${\bf r}_{0}$ and
$\mathbf{v}_{0}$. The goal is to determine the position of the planet ${\bf
r}$ and its velocity $\mathbf{v}$ along the Keplerian orbit after a time $t$.
In order to avoid conversions from Cartesian coordinates into elliptical
elements, we use the Gauss $f$\- and $g$-function333These functions are
different from the time-renormalisation functions used in this article.
formalism (_e.g._ Wisdom & Holman 1991). We have

$\displaystyle{\bf r}$ $\displaystyle=f{\bf r}_{0}+g\mathbf{v}_{0},$
$\displaystyle\mathbf{v}$ $\displaystyle=\dot{f}{\bf
r}_{0}+\dot{g}\mathbf{v}_{0},$ (54)

where the values of $f$, $g$, $\dot{f}$ and $\dot{g}$ depend on $t$, ${\bf
r}_{0}$, and $\mathbf{v}_{0}$ and are given in equations (59). When two
planets encounter each other, their orbits may become hyperbolic. To be able
to resolve these events as well as ejection trajectories, we use a formulation
of the Kepler problem that allows hyperbolic orbits. In order to do so,
Stumpff (1962) developed a general formalism that contains the hyperbolic and
elliptical case in the same equations. Moreover, this approach avoids the
singularity for an eccentricity close to 1. Stumpff introduced special
functions

$c_{n}(z)=\sum_{j=0}^{+\infty}\frac{(-z)^{j}}{(n+2j)!}.$ (55)

The $c$-functions allows us to compute the so-called $G$-functions (Stiefel &
Scheifele 1971), which are defined as

$G_{n}(\beta,X)=X^{n}c_{n}(\beta X^{2}).$ (56)

In this formalism, the Kepler equation takes the form (Stumpff 1962)

$t=r_{0}X+\eta_{0}G_{2}(\beta,X)+\zeta_{0}G_{3}(\beta,X),$ (57)

of unknown $X=\int_{0}^{t}\mathrm{d}t/r$ and where

$\displaystyle\beta=\frac{2\mu}{r_{0}}-\mathbf{v}_{0}^{2},$
$\displaystyle\eta_{0}={\bf r}_{0}\cdot\mathbf{v}_{0},$ (58)
$\displaystyle\zeta_{0}=\mu-\beta r_{0}.$

In Eq. (57), $X$ plays a similar role to the eccentric anomaly in the
classical form. Equation (57) can be solved by the Newton method (Rein &
Tamayo 2015). The new position and velocity are then obtained with Eq. (54)
and

$\displaystyle f$ $\displaystyle=1-\mu\frac{G_{2}(\beta,X)}{r_{0}},$
$\displaystyle\dot{f}$ $\displaystyle=-\frac{\mu G_{1}(\beta,X)}{r_{0}r},$
$\displaystyle g$ $\displaystyle=t-\mu G_{3}(\beta,X),$ $\displaystyle\dot{g}$
$\displaystyle=1-\mu\frac{G_{2}(\beta,X)}{r},$ (59)

where $r=r_{0}+\eta_{0}G_{1}+\zeta_{0}G_{2}$.

###  A.2 Computing the effective time-step

In both the Kepler and the perturbation step, it is important to precisely
compute the effective time-step (20). For the Kepler step in particular,
$\tau_{0}$ depends on the difference $H_{0}-E_{0}$ , where $H_{0}$ and $E_{0}$
are similar. To avoid numerical errors, we compute the initial energy with
compensated summation (Kahan 1965). We save the value and the associated
error. We then evaluate the Keplerian energy $H_{0}$ using compensated
summation and then derive the difference. Compensated summation is also used
to update the positions and velocities, and to integrate the real time
equation.

During the numerical tests, we realised that in Jacobi coordinates, the
perturbation energy $H_{1}$ is most of the time lower by almost an order of
magnitude than the typical planet energy interaction $E_{1}$. On the other
hand, $H_{1}$ and $E_{1}$ are roughly of the same order for heliocentric
coordinates. This shows that the algorithm is less efficient in the detection
of an increase in interaction energy (that monitors the close encounters). To
circumvent this problem, we slightly modify expression (20) for Jacobi
coordinates. The effective step sizes are computed with

$\tau_{0}=\sigma f^{\prime}(H_{0}-E_{0}+E_{1})\text{ and }\tau_{1}=\sigma
f^{\prime}(-(H_{1}-E_{1})).$ (60)

The total energy sum is still zero,

$(H_{0}-E_{0}+E_{1})+(H_{1}-E_{1})=0,$ (61)

which preserves the equation of motion according to Eq. (13). With this
modification, the results between Jacobi and heliocentric coordinates are
comparable.

##  Appendix B MacLachlan high-order schemes

In McLachlan (1995a), the scheme coefficients $w_{i}$ are the coefficients for
a symmetric composition of second-order steps. We list in Table B.1 the
corresponding coefficients $a_{i}$ and $b_{i}$ for the schemes that we used,
$\mathcal{ABA}(6*)$ and $\mathcal{ABA}(8*)$. McLachlan provided 20 significant
digits, therefore we made the computation in quadruple precision and truncated
to the appropriate precision.


Let $\\{\lambda_{i}\left({\bm{S}}\right):i=1,\ldots,P\\}$ be eigenvalues of
the matrix ${\bm{S}}$ in descending order and
$\\{\lambda_{i}\left({\bm{J}}\right):i=1,\ldots,P\\}$ be eigenvalues of the
matrix ${\bm{J}}$ in descending order. Then

$\displaystyle\left|\lambda_{i}\left({\bm{S}}\right)-\lambda_{i}\left({\bm{J}}\right)\right|\leq\left|\left|{\bm{S}}-{\bm{J}}\right|\right|.$
(B.1)

######  Lemma 2 (Lemma S.1 of Lam and Yao (2012)).

Let ${\bm{F}}$ be a $k\times k$ symmetric matrix such that

$\displaystyle{\bm{F}}=\left(\begin{matrix}{\bm{G}}&{\bm{H}}\\\
{\bm{H}}^{\top}&{\bm{D}}\\\ \end{matrix}\right)$ (B.2)

with ${\bm{G}}:k_{1}\times k_{1}$, ${\bm{D}}:k_{2}\times k_{2}$ and
$\lambda_{k_{1}}\left({\bm{G}}\right)>\lambda_{1}\left({\bm{D}}\right)$. Note
that $k_{1}+k_{2}=k$. Then for $1\leq j\leq k_{2}$ ,

$\displaystyle
0\leq\lambda_{j}\left({\bm{D}}\right)-\lambda_{k_{1}+j}\left({\bm{F}}\right)\leq\frac{\lambda_{1}\left({\bm{H}}{\bm{H}}^{\top}\right)}{\lambda_{k_{1}}\left({\bm{G}}\right)-\lambda_{j}\left({\bm{D}}\right)}.$
(B.3)

######  Lemma 3 (Lemma 3 of Lam et al. (2011)).

Suppose ${\bm{F}}$ and ${\bm{F}}+{\bm{E}}$ are $P\times P$ symmetric matrices
and that ${\bm{Q}}=\left({\bm{Q}}_{1},{\bm{Q}}_{2}\right)$, where
${\bm{Q}}_{1}$ has size $P\times k$ and ${\bm{Q}}_{2}$ has size
$P\times(P-k)$, is an orthogonal matrix such that
span$\left({\bm{Q}}_{1}\right)$ is an invariant subspace for the matrix
${\bm{F}}$, that is, ${\bm{F}}\times span\left({\bm{Q}}_{1}\right)\subset
span\left({\bm{F}}\right)$. Partition the matrices
${\bm{Q}}^{\top}{\bm{F}}{\bm{Q}}$ and ${\bm{Q}}^{T}{\bm{E}}{\bm{Q}}$ as
follows.

$\displaystyle{\bm{Q}}^{\top}{\bm{F}}{\bm{Q}}=\left(\begin{matrix}{\bm{D}}_{1}&\textbf{0}\\\
\textbf{0}&{\bm{D}}_{2}\\\ \end{matrix}\right),\ \
{\bm{Q}}^{\top}{\bm{E}}{\bm{Q}}=\left(\begin{matrix}{\bm{E}}_{11}&{\bm{E}}_{21}^{\top}\\\
{\bm{E}}_{21}&{\bm{E}}_{22}\\\ \end{matrix}\right).$ (B.4)

If
$sep\left({\bm{D}}_{1},{\bm{D}}_{2}\right):=\min_{\lambda\in\Lambda\left({\bm{D}}_{1}\right),\mu\in\Lambda\left({\bm{D}}_{2}\right)}\left|\lambda-\mu\right|>0$,
where $\Lambda\left({\bm{D}}_{1}\right)$ denotes the set of eigenvalues of the
matrix ${\bm{D}}_{1}$ and $\left|\left|{\bm{E}}\right|\right|\leq
sep\left({\bm{D}}_{1},{\bm{D}}_{2}\right)/5$, then there exists a matrix
${\bm{P}}:(P-k)\times k$ with

$\displaystyle\left|\left|{\bm{P}}\right|\right|\leq\frac{4\left|\left|{\bm{E}}_{21}\right|\right|}{sep\left({\bm{D}}_{1},{\bm{D}}_{2}\right)}$
(B.5)

such that the columns of the matrix
$\widehat{{\bm{Q}}}_{1}=\left({\bm{Q}}_{1}+{\bm{Q}}_{2}{\bm{P}}\right)\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{-1/2}$
define an orthogonal basis for a subspace that is invariant for the matrix
${\bm{F}}+{\bm{E}}$.

######  Lemma 4.

Under Assumptions 1-8, we have

$\displaystyle\lambda_{j}\left({\bm{L}}_{1}\right)\asymp P^{2-2\delta_{1}},\ \
j=1,\ldots,r_{1}.$ (B.6)
$\displaystyle\lambda_{r_{1}+j}\left({\bm{L}}_{1}\right)\asymp
P^{2-2\delta_{2}},\ \ j=1,\ldots,r_{2};$ (B.7)
$\displaystyle\lambda_{r_{1}+r_{2}+i}\left({\bm{L}}_{1}\right)=o_{p}\left(P^{1-\delta_{1}}\right),\
\ i=1,\ldots,P-(r_{1}+r_{2}).$ (B.8)
$\displaystyle\lambda_{i}\left({\bm{L}}_{2}\right)\asymp P^{2},\ \
i=1,\ldots,r_{2}.$ (B.9)

###### proof of Lemma 4.

Recall the FHFM

$\displaystyle{\bm{y}}_{t}={\bm{B}}{\bm{k}}_{t}^{(1)}+{\bm{A}}{\bm{k}}_{t}^{(2)}+\bm{\varepsilon}_{t}.$
(B.10)

From the expression (B.10), the population covariance matrix of ${\bm{y}}_{t}$
has the following decomposition

$\displaystyle\bm{\Sigma}_{y}(1)={\bm{B}}{\bm{M}}_{1}+{\bm{A}}{\bm{M}}_{2}+\bm{\Sigma}_{\varepsilon}(1),$
(B.11)

where

$\displaystyle{\bm{M}}_{1}=\bm{\Sigma}_{k}^{(1)}(1){\bm{B}}^{\top}+\bm{\Sigma}_{k}^{(12)}(1){\bm{A}}^{\top},\
\ \ \
{\bm{M}}_{2}=\bm{\Sigma}_{k}^{(2)}(1){\bm{A}}^{\top}+\bm{\Sigma}_{k}^{(21)}(1){\bm{B}}^{\top}.$

Based on Lemma 1, we can evaluate the $j$-th eigenvalue of ${\bm{L}}_{1}$
below, $j=1,\ldots,r_{1}$,

$\displaystyle\lambda_{j}\left({\bm{L}}_{1}\right)$
$\displaystyle=\sigma^{2}_{j}\left(\bm{\Sigma}_{y}(1)\right)\geq\left[\sigma_{j}\left({\bm{B}}{\bm{M}}_{1}\right)-\sigma_{1}\left({\bm{A}}{\bm{M}}_{2}+\bm{\Sigma}_{\varepsilon}(1)\right)\right]^{2}$
$\displaystyle\geq\left[\sigma_{j}\left({\bm{B}}{\bm{M}}_{1}\right)-\sigma_{1}\left({\bm{A}}{\bm{M}}_{2}\right)-\sigma_{1}\left(\bm{\Sigma}_{\varepsilon}(1)\right)\right]^{2}$
$\displaystyle=\left[\sigma_{j}\left({\bm{M}}_{1}\right)-\sigma_{1}\left({\bm{M}}_{2}\right)-\sigma_{1}\left(\bm{\Sigma}_{\varepsilon}(1)\right)\right]^{2}$
$\displaystyle\asymp\sigma^{2}\left({\bm{M}}_{1}\right)\geq\left[\sigma_{j}\left(\bm{\Sigma}_{k}^{(1)}(1){\bm{B}}^{\top}\right)-\sigma_{1}\left(\bm{\Sigma}_{k}^{(12)}(1){\bm{A}}^{\top}\right)\right]^{2}$
$\displaystyle\asymp\sigma^{2}_{j}\left(\bm{\Sigma}_{k}^{(1)}(1){\bm{B}}^{\top}\right)=\sigma^{2}_{j}\left(\bm{\Sigma}_{k}^{(1)}(1)\right)\geq\left|\left|\bm{\Sigma}_{k}^{(1)}(1)\right|\right|^{2}_{\min}=P^{2-2\delta_{1}},$

where the first and second inequalities use Lemma 1; the second equality uses
the matrices ${\bm{B}}$ and ${\bm{A}}$ being orthonormal assumed in Assumption
1; and the last inequality and equality both utilize Assumption 2.

Hence, the first $r_{1}$ largest eigenvalues of the matrix ${\bm{L}}_{1}$ have
the order of $P^{2-2\delta_{1}}$.

Now we consider the order of the left $p-r_{1}$ eigenvalues of the matrix
${\bm{L}}_{1}$. In terms of Weyl’s inequality in Lemma 1, we use the
eigenvalues of the matrix
$\widetilde{{\bm{L}}}_{1}=\widetilde{\bm{\Sigma}}_{y}(1)\widetilde{\bm{\Sigma}}_{y}(1)$
to approximate the eigenvalues of ${\bm{L}}_{1}$, where
$\widetilde{\bm{\Sigma}}_{y}(1)={\bm{B}}{\bm{M}}_{1}+{\bm{A}}{\bm{M}}_{2}$. In
fact,

$\displaystyle\quad\left|\lambda_{r_{1}+j}\left({\bm{L}}_{1}\right)-\lambda_{r_{1}+j}\left(\widetilde{{\bm{L}}}_{1}\right)\right|\leq\left|\left|{\bm{L}}_{1}-\widetilde{{\bm{L}}}_{1}\right|\right|$
$\displaystyle\leq\left|\left|{\bm{B}}{\bm{M}}_{1}+{\bm{A}}{\bm{M}}_{2}+\bm{\Sigma}_{\varepsilon}(1)\right|\right|\cdot\left|\left|\bm{\Sigma}_{\varepsilon}(1)\right|\right|+\left|\left|{\bm{B}}{\bm{M}}_{1}+{\bm{A}}{\bm{M}}_{2}\right|\right|\cdot\left|\left|\bm{\Sigma}_{\varepsilon}(1)\right|\right|$
$\displaystyle=o\left(P^{1-\delta_{1}}\right),$ (B.12)

where the last equality uses Assumption 2.

Now we evaluate the order of
$\lambda_{r_{1}+j}\left(\widetilde{{\bm{L}}}_{1}\right)$. Note that the rank
of $\widetilde{{\bm{L}}}_{1}$ is no larger than $r_{1}+r_{2}$. So, when
$j>r_{1}+r_{2}$, $\lambda_{r_{1}+j}\left(\widetilde{{\bm{L}}}_{1}\right)=0$.
Hence, next we investigate the case of $j=1,\ldots,r_{2}$.

Decompose $\widetilde{{\bm{L}}}_{1}$ in the following way.

$\displaystyle\widetilde{{\bm{L}}}_{1}$
$\displaystyle=\left(\begin{matrix}{\bm{B}}&{\bm{A}}\\\
\end{matrix}\right)\left(\begin{matrix}{\bm{M}}_{1}\\\ {\bm{M}}_{2}\\\
\end{matrix}\right)\left(\begin{matrix}{\bm{M}}_{1}^{\top}&{\bm{M}}_{2}^{\top}\end{matrix}\right)\left(\begin{matrix}{\bm{B}}^{\top}\\\
{\bm{A}}^{\top}\\\ \end{matrix}\right)$
$\displaystyle=\left(\begin{matrix}{\bm{B}}&{\bm{A}}\end{matrix}\right)\left(\begin{matrix}{\bm{M}}_{1}{\bm{M}}_{1}^{\top}&{\bm{M}}_{1}{\bm{M}}_{2}^{\top}\\\
{\bm{M}}_{2}{\bm{M}}_{1}^{\top}&{\bm{M}}_{2}{\bm{M}}_{2}^{\top}\\\
\end{matrix}\right)\left(\begin{matrix}{\bm{B}}^{\top}\\\ {\bm{A}}^{\top}\\\
\end{matrix}\right).$ (B.13)

Because

$\displaystyle\left(\begin{matrix}{\bm{B}}^{\top}\\\ {\bm{C}}^{\top}\\\
\end{matrix}\right)\left(\begin{matrix}{\bm{B}}&{\bm{C}}\end{matrix}\right)={\bm{I}},$
(B.14)

we have
$\lambda_{j}\left({\bm{L}}_{1}\right)=\lambda_{j}\left({\bm{M}}\right)$, where

$\displaystyle{\bm{M}}=\left(\begin{matrix}{\bm{M}}_{1}{\bm{M}}_{1}^{\top}&{\bm{M}}_{1}{\bm{M}}_{2}^{\top}\\\
{\bm{M}}_{2}{\bm{M}}_{1}^{\top}&{\bm{M}}_{2}{\bm{M}}_{2}^{\top}\\\
\end{matrix}\right).$ (B.15)

It follows from Lemma 3 and Assumption 2 that

$\displaystyle\lambda_{r_{1}+j}\left({\bm{M}}\right)\leq\lambda_{j}\left({\bm{M}}_{2}{\bm{M}}_{2}^{\top}\right)\asymp
P^{2-2\delta_{2}},$ (B.16)

and

$\displaystyle\lambda_{r_{1}+j}\left({\bm{M}}\right)\geq\sigma_{j}^{2}\left({\bm{M}}_{2}\right)-\frac{\sigma^{2}_{1}\left({\bm{M}}_{1}{\bm{M}}_{2}^{\top}\right)}{\sigma^{2}_{r_{1}}\left({\bm{M}}_{1}\right)-\sigma^{2}_{j}\left({\bm{M}}_{2}\right)}\asymp
P^{2-2\delta_{2}},$ (B.17)

where the last $\asymp$ above uses the fact that
$\sigma_{r_{1}}^{2}\left({\bm{M}}_{1}\right)\asymp P^{2-2\delta_{1}}$,
$\sigma_{j}^{2}\left({\bm{M}}_{2}\right)\asymp P^{2-2\delta_{2}}$ and
$\sigma_{1}^{2}\left({\bm{M}}_{1}{\bm{M}}_{2}^{\top}\right)=O\left(P^{4-2(\delta_{1}+\delta_{2})}\right)$.

Combining (B.16) and (B.17), we can get

$\displaystyle\lambda_{r_{1}+j}\left({\bm{M}}\right)\asymp P^{2-2\delta_{2}},\
\ j=1,\ldots,r_{2}.$ (B.18)

Then it follows from (B), (B.18) and Assumption 2 that

$\displaystyle\lambda_{r_{1}+j}\left({\bm{L}}_{1}\right)\asymp
P^{2-2\delta_{2}},\ \ j=1,\ldots,r_{2};$ (B.19)
$\displaystyle\lambda_{r_{1}+r_{2}+i}\left({\bm{L}}_{1}\right)=o_{p}\left(P^{1-\delta_{1}}\right),\
\ i=1,\ldots,P-(r_{1}+r_{2}).$ (B.20)

Finally, the order of $\lambda_{i}\left({\bm{L}}_{2}\right)$ can be derived
from Proposition 2.1 of Fan et al. (2013) directly. ∎

### Proof of Theorem 1

###### Proof.

Let ${\bm{E}}_{L}^{(1)}=\widehat{{\bm{L}}}_{1}-{\bm{L}}_{1}$ with
$\widehat{{\bm{L}}}_{1}=\widehat{\bm{\Sigma}}_{y}(1)\widehat{\bm{\Sigma}}_{y}(1)^{\top}$
and ${\bm{L}}_{1}=\bm{\Sigma}_{y}(1)\bm{\Sigma}_{y}(1)^{\top}$. First we
evaluate the order of $\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|$. In terms
of simple calculations, we have

$\displaystyle\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|\leq\left|\left|\widehat{\bm{\Sigma}}_{y}(1)-\bm{\Sigma}_{y}(1)\right|\right|^{2}+2\left|\left|\widehat{\bm{\Sigma}}_{y}(1)\right|\right|\cdot\left|\left|\widehat{\bm{\Sigma}}_{y}(1)-\bm{\Sigma}_{y}(1)\right|\right|.$
(B.21)

In terms of (B.6) in Lemma 4, we have
$\left|\left|\bm{\Sigma}_{y}(1)\right|\right|\asymp P^{1-\delta_{1}}$. From
(B.11), we can get

$\displaystyle\left|\left|\widehat{\bm{\Sigma}}_{y}(1)-\bm{\Sigma}_{y}(1)\right|\right|$
$\displaystyle\leq\left|\left|\widehat{{\bm{M}}}_{1}-{\bm{M}}_{1}\right|\right|+\left|\left|\widehat{{\bm{M}}}_{2}-{\bm{M}}_{2}\right|\right|$
$\displaystyle\quad+\left|\left|\widehat{\bm{\Sigma}}_{\varepsilon}(1)-\bm{\Sigma}_{\varepsilon}(1)\right|\right|,$
(B.22)

where

$\displaystyle\widehat{{\bm{M}}}_{1}=\widehat{\bm{\Sigma}}_{k}^{(1)}(1){\bm{B}}^{\top}+\widehat{\bm{\Sigma}}_{k}^{(12)}(1){\bm{A}}^{\top},\
\ \ \
\widehat{{\bm{M}}}_{2}=\widehat{\bm{\Sigma}}_{k}^{(2)}(1){\bm{A}}^{\top}+\widehat{\bm{\Sigma}}_{k}^{(21)}(1){\bm{B}}^{\top},$

with $\widehat{\bm{\Sigma}}_{k}^{(1)}(1)$,
$\widehat{\bm{\Sigma}}_{k}^{(12)}(1)$, $\widehat{\bm{\Sigma}}_{k}^{(2)}(1)$
and $\widehat{\bm{\Sigma}}_{k}^{(21)}(1)$ are the sample covariances
corresponding to the population covariances $\bm{\Sigma}_{k}^{(1)}(1)$,
$\bm{\Sigma}_{k}^{(12)}(1)$, $\bm{\Sigma}_{k}^{(2)}(1)$ and
$\bm{\Sigma}_{k}^{(21)}(1)$, respectively.

Hence, we evaluate (B) further

$\displaystyle\left|\left|\widehat{\bm{\Sigma}}_{y}(1)-\bm{\Sigma}_{y}(1)\right|\right|$
$\displaystyle\leq\left|\left|\widehat{\bm{\Sigma}}_{k}^{(1)}(1)-\bm{\Sigma}_{k}^{(1)}(1)\right|\right|+\left|\left|\widehat{\bm{\Sigma}}_{k}^{(12)}(1)-\bm{\Sigma}_{k}^{(12)}(1)\right|\right|$
$\displaystyle\quad+\left|\left|\widehat{\bm{\Sigma}}_{k}^{(2)}(1)-\bm{\Sigma}_{k}^{(2)}(1)\right|\right|+\left|\left|\widehat{\bm{\Sigma}}_{k}^{(21)}(1)-\bm{\Sigma}_{k}^{(21)}(1)\right|\right|$
$\displaystyle\quad+\left|\left|\widehat{\bm{\Sigma}}_{\varepsilon}(1)-\bm{\Sigma}_{\varepsilon}(1)\right|\right|$
$\displaystyle=O_{p}\left(\frac{P^{1-\delta_{1}}}{T^{1/2}}\right)+O_{p}\left(\frac{P^{1-\delta_{2}}}{T^{1/2}}\right)+O_{p}\left(\frac{P}{T}\right)$
$\displaystyle=O_{p}\left(\max\left(\frac{P}{T},\frac{P^{1-\delta_{1}}}{T^{1/2}}\right)\right),$
(B.23)

where the last second equality uses (A8) of Lam et al. (2011) which
demonstrates
$\left|\left|\widehat{\bm{\Sigma}}_{\varepsilon}(1)-\bm{\Sigma}_{\varepsilon}(1)\right|\right|=O_{p}\left(\frac{P}{T}\right)$.

Then it follows from (B.21) and (B) that

$\displaystyle\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|$ $\displaystyle=$
$\displaystyle
O_{p}\left(\max\left(\frac{P^{2}}{T^{2}},\frac{P^{2-2\delta_{1}}}{T},\frac{P^{2-\delta_{1}}}{T},\frac{P^{2-2\delta_{1}}}{T^{1/2}}\right)\right)$
(B.24) $\displaystyle=$ $\displaystyle
O_{p}\left(\frac{P^{2-2\delta_{1}}}{T^{1/2}}\right),$ (B.25)

where the last equality uses the assumption that
$P^{\delta_{1}}=o\left(T^{1/2}\right)$.

Now we use Lemma 3 to get the order of estimated factor loadings. In Lemma 3,
let ${\bm{F}}$ and ${\bm{E}}$ be ${\bm{L}}_{1}$ and
$\widehat{{\bm{L}}}_{1}-{\bm{L}}_{1}$, respectively. Let $k$ in Lemma 3 equal
to $r_{1}$. Then we have, from (B.6),

$\displaystyle sep({\bm{D}}_{1},{\bm{D}}_{2})\asymp P^{2-2\delta_{1}},$ (B.26)

where the defition of $sep(\cdot,\cdot)$ is provided in Lemma 3. Then
${\bm{E}}_{L}^{(1)}$ and
$sep({\bm{D}}_{1},{\bm{D}}_{2})$ satisfies

$\displaystyle\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|=o_{p}\left(sep({\bm{D}}_{1},{\bm{D}}_{2})\right)\leq\frac{sep({\bm{D}}_{1},{\bm{D}}_{2})}{5}.$
(B.27)

Hence Lemma 3 tells us that, there exists a matrix ${\bm{P}}:(P-r_{1})\times
r_{1}$ such that

$\displaystyle\left|\left|{\bm{P}}\right|\right|\leq\frac{4}{sep\left({\bm{D}}_{1},{\bm{D}}_{2}\right)}\cdot\left|\left|\left({\bm{E}}_{L}^{(1)}\right)_{21}\right|\right|\leq\frac{4\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|}{sep({\bm{D}}_{1},{\bm{D}}_{2})}$
(B.28)

and then
$\widehat{{\bm{B}}}=\left({\bm{B}}+{\bm{B}}^{c}{\bm{P}}\right)\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{-1/2}$
is an estimator of ${\bm{B}}$ with ${\bm{B}}^{c}$ being ${\bm{Q}}_{2}$ in
Lemma 3. In view of this, the rate of convergence for $\widehat{{\bm{B}}}$ can
be calculated as

$\displaystyle\left|\left|\widehat{{\bm{B}}}-{\bm{B}}\right|\right|$
$\displaystyle=\left|\left|\left({\bm{B}}+{\bm{B}}^{c}{\bm{P}}\right)\left({\bm{I}}+{\bm{P}}^{T}{\bm{P}}\right)^{-1/2}-{\bm{B}}\right|\right|$
$\displaystyle=\left|\left|\left[\left({\bm{B}}+{\bm{B}}{\bm{P}}\right)-{\bm{B}}\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{1/2}\right]\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{-1/2}\right|\right|$
$\displaystyle=\left|\left|\left({\bm{B}}\left[{\bm{I}}-\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{1/2}\right]+{\bm{B}}^{c}{\bm{P}}\right)\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{-1/2}\right|\right|$
$\displaystyle\leq\left|\left|{\bm{I}}-\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{1/2}\right|\right|+\left|\left|{\bm{P}}\right|\right|\leq
2\left|\left|{\bm{P}}\right|\right|,$ (B.29)

where the last second equality uses the fact that ${\bm{B}}$ and
${\bm{B}}^{c}$ are orthonormal; and the last equality uses the fact that

$\displaystyle\left|\left|{\bm{I}}-\left({\bm{I}}+{\bm{P}}^{\top}{\bm{P}}\right)^{1/2}\right|\right|$
$\displaystyle=$ $\displaystyle
1-\left(1+\lambda_{\min}\left({\bm{P}}^{\top}{\bm{P}}\right)\right)^{1/2}$
(B.30) $\displaystyle\leq$
$\displaystyle\lambda^{1/2}_{\max}\left({\bm{P}}^{\top}{\bm{P}}\right).$
(B.31)

Therefore, by (B) and (B.28), we obtain

$\displaystyle\left|\left|\widehat{{\bm{B}}}-{\bm{B}}\right|\right|=O_{P}\left(\frac{\left|\left|{\bm{E}}_{L}^{(1)}\right|\right|}{sep({\bm{D}}_{1},{\bm{D}}_{2})}\right)=O_{p}\left(\frac{1}{\sqrt{T}}\right).$
(B.32)

For the second factor model part, the estimation is to conduct principal
component analysis on the residual of the first step, i.e. estimating the
factor model

$\displaystyle\widehat{{\bm{u}}}_{t}={\bm{A}}{\bm{k}}_{t}^{(2)}+\bm{\eta}_{t},\
\ t=1,2,\ldots,T,$ (B.33)

where
$\widehat{{\bm{u}}}_{t}={\bm{y}}_{t}-\widehat{{\bm{B}}}\widehat{{\bm{k}}}_{t}^{(1)}$,
$\bm{\eta}_{t}$ is the new error component in the estimation at the second
step.

In order to derive the rate of convergence for $\widehat{{\bm{A}}}$, we also
utilize Lemma 3. Now let ${\bm{F}}$ and ${\bm{E}}$ in Lemma 3 are
${\bm{L}}_{2}$ and ${\bm{E}}_{L}^{(2)}:=\widehat{{\bm{L}}}_{2}-{\bm{L}}_{2}$,
respectively. Let $k$ in Lemma 3 equal to $r_{2}$.

First, we evaluate $\left|\left|{\bm{E}}_{L}^{(2)}\right|\right|$. Based on
(B.33), we have

$\displaystyle\left|\left|{\bm{E}}_{L}^{(2)}\right|\right|\leq\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\bm{\Sigma}_{u}(0)\right|\right|^{2}+2\left|\left|\bm{\Sigma}_{u}(0)\right|\right|\cdot\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\bm{\Sigma}_{u}(0)\right|\right|,$
(B.34)

where $\bm{\Sigma}_{u}(0)$ is the population covariance matrix of
${\bm{u}}_{t}$ and $\widehat{\bm{\Sigma}}_{\widehat{u}}(0)$ is the sample
covariance matrix of $\widehat{{\bm{u}}}_{t}$. Based on Assumption 3 and
Proposition 2.1 of Fan et al. (2013), we know that
$\left|\left|\bm{\Sigma}_{u}(0)\right|\right|\asymp p$. For the term
$\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\bm{\Sigma}_{u}(0)\right|\right|$,
we evaluate its order as follows.

$\displaystyle\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\bm{\Sigma}_{u}(0)\right|\right|$
$\displaystyle\leq\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\widehat{\bm{\Sigma}}_{u}(0)\right|\right|+\left|\left|\widehat{\bm{\Sigma}}_{u}(0)-\bm{\Sigma}_{u}(0)\right|\right|$
$\displaystyle\leq\frac{1}{T}\left|\left|\widehat{{\bm{B}}}\widehat{{\bm{K}}}^{(1)}-{\bm{B}}{\bm{K}}^{(1)}\right|\right|^{2}+\frac{2}{T}\left|\left|{\bm{B}}{\bm{K}}^{(1)}\right|\right|\cdot\left|\left|\widehat{{\bm{B}}}\widehat{{\bm{K}}}^{(1)}-{\bm{B}}{\bm{K}}^{(1)}\right|\right|$
$\displaystyle\quad+\left|\left|\widehat{\bm{\Sigma}}_{u}(0)-\bm{\Sigma}_{u}(0)\right|\right|,$
(B.35)

where
$\widehat{{\bm{K}}}^{(1)}=\left(\widehat{{\bm{k}}}_{1}^{(1)},\widehat{{\bm{k}}}_{2}^{(1)},\ldots,\widehat{{\bm{k}}}_{T}^{(1)}\right)$
and
${\bm{K}}^{(1)}=\left({\bm{k}}_{1}^{(1)},{\bm{k}}_{2}^{(1)},\ldots,{\bm{k}}_{T}^{(1)}\right)$.

From Assumption 2 and (B.32), it can be derived that

$\displaystyle\frac{1}{\sqrt{T}}\left|\left|\widehat{{\bm{B}}}\widehat{{\bm{K}}}^{(1)}-{\bm{B}}{\bm{K}}^{(1)}\right|\right|$
$\displaystyle=$ $\displaystyle
O_{p}\left(\frac{P^{1/2-\delta_{1}/2}}{\sqrt{T}}\right),$
$\displaystyle\frac{1}{\sqrt{T}}\left|\left|{\bm{B}}{\bm{K}}^{(1)}\right|\right|$
$\displaystyle=$ $\displaystyle O_{p}\left(P^{1/2-\delta_{1}/2}\right).$
(B.36)

Similar to (B), we can also get

$\displaystyle\left|\left|\widehat{\bm{\Sigma}}_{u}(0)-\bm{\Sigma}_{u}(0)\right|\right|$
$\displaystyle\leq\left|\left|\widehat{\bm{\Sigma}}_{k}^{(2)}(0)-\bm{\Sigma}_{k}^{(2)}(0)\right|\right|+\left|\left|\widehat{\bm{\Sigma}}_{\varepsilon}(0)-\bm{\Sigma}_{\varepsilon}(0)\right|\right|$
$\displaystyle=O_{p}\left(\frac{P}{T^{1/2}}\right)+O_{p}\left(\frac{P}{T}\right).$
(B.37)

In view of (B), (B) and (B), we can get

$\displaystyle\left|\left|\widehat{\bm{\Sigma}}_{\widehat{u}}(0)-\bm{\Sigma}_{u}(0)\right|\right|$
$\displaystyle=O_{P}\left(\frac{p^{1-\delta_{1}}}{\sqrt{T}}\right)+O_{p}\left(\frac{P}{\sqrt{T}}\right)+O_{p}\left(\frac{P}{T}\right)$
(B.38) $\displaystyle=O_{p}\left(\frac{P}{\sqrt{T}}\right).$

The order of $\left|\left|{\bm{E}}_{L}^{(2)}\right|\right|$ is obtained from
(B.34) and (B.38), i.e.

$\displaystyle\left|\left|{\bm{E}}_{L}^{(2)}\right|\right|=O_{p}\left(\frac{P^{2}}{\sqrt{T}}\right).$
(B.39)

Moreover, it follows from Proposition 2.1 of Fan et al. (2013) that

$\displaystyle sep({\bm{D}}_{1},{\bm{D}}_{2})\asymp P^{2}.$ (B.40)

Here ${\bm{D}}_{1}$ in Lemma 3 is the diagonal matrix corresponding to the
orthogonal matrix ${\bm{A}}$. Then we can get from Lemma 3 that

$\displaystyle\left|\left|\widehat{{\bm{A}}}-{\bm{A}}\right|\right|=O_{p}\left(\frac{\left|\left|{\bm{E}}_{L}^{(2)}\right|\right|}{sep\left({\bm{D}}_{2},{\bm{D}}_{1}\right)}\right)=O_{p}\left(\frac{1}{\sqrt{T}}\right).$
(B.41)

∎

### Proof of Theorem 2

###### Proof.

It follows from the model (2.1) and (2.2) that

$\displaystyle\frac{1}{P}\left|\left|\widehat{{\bm{y}}}_{t-1}(1)-{\bm{y}}_{t}\right|\right|^{2}=\frac{1}{P}\sum^{P}_{i=1}\left(\widehat{y}_{i,t-1}(1)-y_{it}\right)^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{P}\sum^{P}_{i=1}\left(\widehat{{\bm{b}}}^{\top}_{i}\widehat{{\bm{k}}}_{t-1}^{(1)}(1)+\widehat{{\bm{a}}}_{i}^{\top}\widehat{{\bm{k}}}^{(2)}_{t-1}(1)-{\bm{b}}_{i}{\bm{k}}_{t}^{(1)}-{\bm{a}}^{\top}_{i}{\bm{k}}_{t}^{(2)}-\varepsilon_{it}\right)^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{P}\sum^{P}_{i=1}\left(\left(\widehat{{\bm{b}}}_{i}-{\bm{b}}_{i}\right)^{\top}\widehat{{\bm{k}}}^{(1)}_{t-1}(1)\right)^{2}+\frac{1}{P}\sum^{P}_{i=1}\left(\left(\widehat{{\bm{a}}}_{i}-{\bm{a}}_{i}\right)^{\top}\widehat{{\bm{k}}}^{(2)}_{t-1}(1)\right)^{2}$
$\displaystyle+\frac{1}{P}\sum^{P}_{i=1}\left[{\bm{b}}_{i}^{\top}\left(\widehat{{\bm{k}}}_{t-1}^{(1)}(1)-{\bm{k}}^{(1)}_{t}\right)\right]^{2}+\frac{1}{P}\sum^{P}_{i=1}\left[{\bm{a}}_{i}^{\top}\left(\widehat{{\bm{k}}}_{t-1}^{(2)}(1)-{\bm{k}}^{(2)}_{t}\right)\right]^{2}+\frac{1}{P}\sum^{P}_{i=1}\varepsilon_{it}^{2}$
$\displaystyle=:$ $\displaystyle\sum^{5}_{d=1}C_{d}.$

From Theorem 1, we have the orders of $C_{1}$ and $C_{2}$ are $\frac{1}{T}$.
The orders of $C_{3}$ and $C_{4}$ reply on the forecasting error of two kinds
of factors $\widehat{{\bm{k}}}^{(j)}_{t-1}(1)$ with $j=1,2$. The last term
$\varepsilon_{it}$ comes from the second factor model error in the second
step. ∎