∎

11institutetext: Hanlun Lei (✉)22institutetext: School of Astronomy and Space
Science, Nanjing University, Nanjing 210023, China
Key Laboratory of Modern Astronomy and Astrophysics in Ministry of Education,
Nanjing University, Nanjing 210023, China
22email<EMAIL_ADDRESS>
# Structures of secular resonances for inner test particles in hierarchical
planetary systems

Hanlun Lei

(Received: date / Accepted: date)

###### Abstract

In this study, dynamics of secular resonances for inner test particles are
investigated under the octupole-level approximation by taking non-perturbative
approaches. In practice, webs of the major secular resonances are produced by
identifying families of stable periodic orbits and the associated stable
libration zones are obtained by analysing Poincaré surfaces of section. By
taking different values of the factor $\epsilon$ ($\epsilon$ measures the
contribution of octupole terms), the influences of the octupole-order terms
upon the dynamical structures are evaluated. Under the condition of
$\epsilon=0$ (no octupole-order contribution), the dynamical model is totally
integrable and there is only Kozai resonance arising in the phase space. When
the factor $\epsilon$ is different from zero, the dynamical structure in the
phase space becomes complicated due to varieties of secular resonances
appearing. Numerical results further indicate that (a) distributions of
libration centres and stable libration zones remain qualitatively similar with
different values of $\epsilon$, (b) Kozai resonance disappears due to the
chaotic motion in the low-eccentricity region, and (c) the chaotic area
arising in the low-eccentricity region increases with the factor $\epsilon$.
Secular resonances are the source of many important dynamical phenomena, such
as chaos, orbit alignment and orbit flipping, and thus the results presented
in this work could be useful to understand the secular dynamics for those
high-eccentricity and/or high-inclination objects in hierarchical planetary
systems.

###### Keywords:

Secular resonances Poincaré surfaces of section Families of periodic orbits

††journal: Celestial Mechanics and Dynamical Astronomy

##  1 Introduction

Hierarchical planetary systems are common in the Universe, covering a large
range of astrophysical systems from the satellite and planet scales to
supermassive black hole (Naoz 2016). In particular, when one of the triple
system is assumed as a massless particle, the system reduces to the so-called
hierarchical restricted three-body problem.

To explore the secular dynamics of test particles on timescales much longer
than their orbital periods, the short-period effects over the orbital periods
of the inner and outer binaries are averaged out by means of averaging
approach (Lidov 1962; Kozai 1962; Ford et al 2000; Naoz et al 2013; Luo et al
2016; Lei 2019, 2020). Under the standard double-averaged model, Lidov (1962)
studied the secular dynamics of artificial satellites perturbed by the Sun and
Moon and, at the same age, Kozai (1962) investigated the secular influences of
Jupiter’s gravitational perturbations upon those highly inclined and eccentric
asteroids in our Solar system. It is known that, when the inclination is
greater than $39.2^{\circ}$, the eccentricity and inclination of satellites or
asteroids would undergo a coupled oscillation (Lidov 1962; Kozai 1962; Thomas
and Morbidelli 1996). The coupled oscillation between eccentricity and
inclination is due to the existence of Kozai resonance, which happens between
the longitude of pericentre $\varpi$ and longitude of ascending node $\Omega$.
The resulting mechanism exciting eccentricity and/or inclination of test
particles is called von Zeipel–Lidov–Kozai mechanism (Ito and Ohtsuka 2019).
In both Lidov (1962) and Kozai (1962), the following two assumptions are made:
(a) the test-particle approximation for the satellites or asteroids and (b)
circular-orbit approximation for the orbit of perturber. The circular-orbit
approximation would lead to axisymmetric potential where the longitude of
ascending node $\Omega$ disappears from the Hamiltonian, so that the dynamical
model becomes an integrable system. Usually, the double-averaged disturbing
function is truncated at the lowest order in the semimajor axis ratio, which
is called the test particle quadrupole approximation (TPQ) by Naoz (2016).

However, if the circular-orbit approximation is relaxed, the next level
(octupole order) of approximation needs to be taken into consideration and the
resulting dynamical model is of two degrees of freedom with $\Omega$ and
$\omega$ as angular coordinates. Under such a relaxation, the standard von
Zeipel–Lidov–Kozai mechanism naturally becomes the eccentric version.

According to the octupole-level secular perturbation theory, Lee and Peale
(2003) studied the apsidal secular resonance with critical argument of
$\sigma=\varpi_{1}-\varpi_{2}$ under coplanar, hierarchical, two-planet
systems and they showed that this resonance causes large amplitude variations
of eccentricities. Under the octupole-order approximation, Lithwick and Naoz
(2011) proved that the particle’s vertical angular momentum is no longer a
constant and its Kozai oscillations are modulated on timescales much longer
than Kozai period. Particularly, the modulation of Kozai cycles could make the
particle’s orbit flip from prograde to retrograde and back again, and at the
flipping instant the orbit could reach arbitrarily high eccentricities (this
is called the eccentric von Zeipel–Lidov–Kozai mechanism). Such an eccentric
von Zeipel–Lidov–Kozai mechanism provides a possible clue for the formation of
retrograde ‘hot Jupiters’ (Naoz et al 2011): (a) the secular perturbation
coming from the perturber on an elliptic orbit causes orbit flipping, which is
accompanied by an excursion to high eccentricity and thus low altitude of
periastron, and (b) strong planet–star tidal interactions damp the
eccentricity of the inner orbit, forming a retrograde hot Jupiter. Considering
the fact that the octupole-order potential causes a slow and cyclic modulation
of the Kozai–Lidov cycles, Katz et al (2011) studied the very long-term
evolution analytically by averaging the secular equations of motion over the
Kozai–Lidov cycles (at the same time introducing a new constant of motion) and
derived an analytical criterion for orbit flipping. Li et al (2014) studied
the eccentric von Zeipel–Lidov–Kozai mechanism systematically by taking
advantage of Poincaré surfaces of section and the Lyapunov exponents and they
found that, when the mutual inclination is high, the resonances due to the
octupole potential are the main causes that lead to orbits flipping,
significant eccentricity excitation, and chaotic behaviors. Sidorenko (2018)
interpreted the eccentric von Zeipel–Lidov–Kozai effect that causes orbit
flipping under the octupole-order approximation as a resonance phenomenon with
the critical argument $\sigma=\Omega+{\rm sign}(\cos{i})\omega$ librating
around $0^{\circ}$ or $180^{\circ}$ ($\rm sign$ is the sign function, $\Omega$
the longitude of ascending node, $\omega$ the argument of pericentre and $i$
the inclination).

Regarding the dynamical model truncated at the octupole order, the previous
studies indicate that the secular resonances govern the long-term dynamics,
including the coupled excitation of eccentricity and inclination, orbit
flipping and chaotic behaviours. However, no systematic mapping of the secular
resonances has been conducted yet. In this context, we explore the global
dynamical structures (or webs of secular resonances) in the phase space with
an aim at understanding the secular dynamics for those high-eccentricity
and/or high-inclination objects in hierarchical restricted three-body systems.

The remaining part of this work is organised as follows. In Sect. 2, the
dynamical model for long-term evolutions of test particles is briefly
introduced and Kozai dynamics under the quadrupole-order approximation is
discussed and, in Sect. 3, the Poincaré surfaces of sections at different
levels of Hamiltonian are produced and analysed. In Sect. 4, families of
symmetric periodic orbits (or webs of secular resonances) are determined and,
in Sect. 5, the boundaries of stable libration zones are produced with
different system parameters. Finally, the conclusion is summarised in Sect. 6.

##  2 Hamiltonian model and Kozai dynamics

The dynamical model truncated at the octupole order is briefly presented, and
then under the quadrupole-level approximation the Kozai dynamics is recalled
for the purpose of comparison.

###  2.1 Hamiltonian function

We take the hierarchical restricted three-body problem composed of two
primaries (a central body and a distant planet) and an inner test particle as
the fundamental model. Under this model, the primaries move around their
barycentre in fixed elliptic orbits and the test particle moves around the
central body in a perturbed Keplerian orbit.

To describe the orbits of test particles and perturber, we choose an inertial
right-handed reference frame originated at the central body, denoted by
$o$-$xyz$, with the perturber’s orbit as the fundamental plane, the $z$-axis
along the total angular momentum vector and the $x$-axis along the
eccentricity vector of the perturber. Under the coordinate system, the orbits
of the test particle (and the perturber) are characterised by six orbital
elements, including the semimajor axis $a$ ($a_{p}$), eccentricity $e$
($e_{p}$), inclination $i$ ($i_{p}$), longitude of ascending node $\Omega$
($\Omega_{p}$), argument of pericentre $\omega$ ($\omega_{p}$) and the mean
anomaly $M$ ($M_{p}$). Unless otherwise specified, in the entire work we
utilise the variables with subscript ‘$p$’ for the perturber and the ones
without any subscript for test particles. Due to the choice of the reference
plane and the $x$-axis, both the inclination of perturber $i_{p}$ and the
longitude of pericentre $\varpi_{p}$ are zero.

In hierarchical configurations, it requires that the semimajor axis of the
test particle is much smaller than that of the secondary, i.e., the semimajor
axis ratio $\alpha=a/a_{p}$ is a small parameter. Following Harrington (1968,
1969), the Hamiltonian, describing the motion of test particles perturbed by a
distant body, can be written as

${\cal H}=-\frac{\mu}{{2a}}-\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\sum\limits_{n=2}^{\infty}{{\alpha^{n}}{{\left({\frac{r}{a}}\right)}^{n}}{{\left({\frac{{{a_{p}}}}{{{r_{p}}}}}\right)}^{n+1}}{P_{n}}\left({\cos\psi}\right)},$
(1)

where $\mu={\cal G}m_{0}$ is the gravitational parameter of the central body,
$m_{0}$ and $m_{p}$ are the mass of the central body and perturber, $r$ and
$r_{p}$ are the distances of the test particle and perturber relative to the
central body, $\psi$ stands for the angular separation between their position
vectors, and $P_{n}(\cos{\psi})$ are the Legendre polynomial of $\cos{\psi}$
with degree $n$. The second term in the right hand of equation (1) is called
the disturbing function.

To study the secular dynamics of test particles in hierarchical systems, it is
required to make a secular approximation for the Hamiltonian function (Ford et
al 2000; Naoz et al 2013; Naoz 2016). Usually, secular approximation is
achieved by means of phase averaging over the orbital periods of the inner and
outer binaries as follows:

${\cal
H}^{*}=\frac{1}{{4{\pi^{2}}}}\int\limits_{0}^{2\pi}{\int\limits_{0}^{2\pi}{{\cal
H}{\rm d}M{\rm d}{M_{p}}}}.$

To make sure that the double-averaged Hamiltonian could effectively predict
the long-term evolution of the test particle, it is required that the mass of
the perturber should be much smaller than that of the central body. Otherwise,
the effects filtered by double averaging need to be taken into account in
formulating the dynamical model of long-term evolution (Ćuk and Burns 2004;
Luo et al 2016; Lei et al 2018; Lei 2019; Hamers and Samsing 2019).

In the averaged model, the mean anomaly of test particle disappears from the
Hamiltonian, so that its semimajor axis remains unchanged in the long-term
evolution. Let us further truncate the Hamiltonian at the third order in
$\alpha$ and eliminate those constant terms, so that the octupole-level
approximation of secular Hamiltonian becomes (Ford et al 2000; Lithwick and
Naoz 2011; Naoz 2016)

${\cal H}^{*}=-{{\cal C}_{0}}\left({{F_{\rm quad}}+\epsilon{F_{\rm
oct}}}\right)$

where the coefficient ${{\cal C}_{0}}$ is

${{\cal C}_{0}}=\frac{3}{8}\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}{\left({\frac{a}{{{a_{p}}}}}\right)^{2}}\frac{1}{{{{\left({1-e_{p}^{2}}\right)}^{3/2}}}},$

the quadrupole- and octupole-order terms ${F_{\rm quad}}$ and ${F_{\rm oct}}$
are given by

${F_{\rm
quad}}=-\frac{1}{2}{e^{2}}+{\cos^{2}}i+\frac{3}{2}{e^{2}}{\cos^{2}}i+\frac{5}{2}{e^{2}}\left({1-{{\cos}^{2}}i}\right)\cos\left({2\omega}\right),$
$\displaystyle{F_{\rm oct}}$
$\displaystyle=\frac{5}{{16}}\left({e+\frac{3}{4}{e^{3}}}\right)\left({1-11\cos
i-5{{\cos}^{2}}i+15{{\cos}^{3}}i}\right)\cos\left({\omega-\Omega}\right)$
$\displaystyle+\frac{5}{{16}}\left({e+\frac{3}{4}{e^{3}}}\right)\left({1+11\cos
i-5{{\cos}^{2}}i-15{{\cos}^{3}}i}\right)\cos\left({\omega+\Omega}\right)$
$\displaystyle-\frac{{175}}{{64}}{e^{3}}\left({1-\cos
i-{{\cos}^{2}}i+{{\cos}^{3}}i}\right)\cos\left({3\omega-\Omega}\right)$
$\displaystyle-\frac{{175}}{{64}}{e^{3}}\left({1+\cos
i-{{\cos}^{2}}i-{{\cos}^{3}}i}\right)\cos\left({3\omega+\Omega}\right)$

and the factor $\epsilon$, measuring the significance of octupole-order term,
is given by (Lithwick and Naoz 2011)

$\epsilon=\frac{a}{{{a_{p}}}}\frac{{{e_{p}}}}{{1-e_{p}^{2}}}=\alpha\frac{{{e_{p}}}}{{1-e_{p}^{2}}}.$

For the same dynamical model, the double-averaged Hamiltonian can be
formulated up to an arbitrary order in the semimajor axis ratio $\alpha$, as
given by Lei (2021).

Evidently, the octupole-level contribution to the secular Hamiltonian is
determined by the magnitude of $\epsilon$, which is positively correlated to
the semimajor axis ratio $\alpha=a/a_{p}$ and the eccentricity of the
perturber’s orbit $e_{p}$. In Fig. 1, several level curves of $\epsilon$ are
given in the space $(\alpha,e_{p})$. Obviously, each value of $\epsilon$
stands for a series of dynamical models. In particular, when the perturber is
on a circular orbit (corresponding to $\epsilon=0$), the octupole-level terms
vanish. In this case, the dynamical model reduces to the well-known circular
restricted three-body problem (CRTBP), in which the octupole-order term has no
contribution.

Figure 1: Representative values of the system parameter $\epsilon$ adopted in
this study. As a special case, the dynamical model specified by $\epsilon=0$
will be used in Sect. 2.2. Dynamical models specified by
$\epsilon=0.03,0.05,0.07,0.1$ are to be used in Sects 3–5.

Considering the fact that the coefficient ${\cal C}_{0}$ is a constant in the
long-term evolution, we can scale the Hamiltonian ${\cal H}^{*}$ by ${\cal
C}_{0}$ as follows (for simplicity, the notation ${\cal H}$ is used to
represent the scaled Hamiltonian):

${\cal H}={{\cal H}^{*}}/{{\cal C}_{0}}=-\left({{F_{\rm quad}}+\epsilon{F_{\rm
oct}}}\right).$ (2)

The normalised Hamiltonian given by equation (2) has been widely used in
previous works (Lithwick and Naoz 2011; Li et al 2014; Naoz 2016). The scaling
of the Hamiltonian given by Eq. (2) is equivalent to changing the unit of
time.

To formulate the equations of motion, let’s adopt the scaled Delaunay
variables (the action variables are scaled by the constant factor $\sqrt{\mu
a}$) in the following manner (Lithwick and Naoz 2011):

$\displaystyle g=\omega,\quad G=\sqrt{1-{e^{2}}},$ $\displaystyle
h=\Omega,\quad H=G\cos i.$

Thus, the Hamiltonian ${\cal H}$ can be re-organised as a function of the
defined Delaunay variables by ${\cal H}(G,H,g,h)$. The dynamical model holds
the following symmetric properties:

${\cal H}(G,H,g,h)={\cal H}(G,-H,2\pi-g,h)={\cal H}(G,-H,g,2\pi-h)={\cal
H}(G,H,g+\pi,h+\pi)$

where the first two symmetries can be found in Sidorenko (2018). The (scaled)
Hamiltonian canonical relations hold

$\displaystyle\frac{{{\rm d}g}}{{{\rm d}t}}=$
$\displaystyle\frac{{\partial{\cal H}}}{{\partial G}},\quad\frac{{{\rm
d}G}}{{{\rm d}t}}=-\frac{{\partial{\cal H}}}{{\partial g}},$ (3)
$\displaystyle\frac{{{\rm d}h}}{{{\rm d}t}}=$
$\displaystyle\frac{{\partial{\cal H}}}{{\partial H}},\quad\frac{{{\rm
d}H}}{{{\rm d}t}}=-\frac{{\partial{\cal H}}}{{\partial h}},$

where the scaled time is (Lithwick and Naoz 2011)

$\displaystyle t$ $\displaystyle=\frac{{\cal C}_{0}}{\sqrt{\mu a}}t_{r}$
$\displaystyle=\frac{3}{8}\frac{{{m_{p}}}}{{{m_{0}}}}n{\left({\frac{a}{{{a_{p}}}}}\right)^{3}}\frac{1}{{{{\left({1-e_{p}^{2}}\right)}^{3/2}}}}{t_{r}}$

with $n$ as the mean motion frequency of the test particle moving around the
central body and $t_{r}$ as the true time.

To remove the singularity of variables at small eccentricities and
inclinations, the following variables are also used for describing the motion
of test particles:

$\displaystyle k$ $\displaystyle=e\cos\omega,\quad h=e\sin\omega,$
$\displaystyle q$ $\displaystyle=\sin{\frac{i}{2}}\cos\Omega,\quad
p=\sin{\frac{i}{2}}\sin\Omega.$

It is not difficult to realise the transformations among the following sets of
variables:

$(e,i,\Omega,\omega)\Leftrightarrow(G,H,g,h)\Leftrightarrow(k,h,q,p).$

###  2.2 Kozai dynamics

For completeness, here we make some discussions about Kozai dynamics in the
special case of $\epsilon=0$, which is equivalent to the quadrupole-level
approximation. When $\epsilon=0$, the (scaled) Hamiltonian given by Eq. (2)
can be organised as

${\cal
H}=\frac{1}{2}\left({1-{G^{2}}+3{H^{2}}-5\frac{{{H^{2}}}}{{{G^{2}}}}}\right)-\frac{5}{2}\left({1-{G^{2}}+{H^{2}}-\frac{{{H^{2}}}}{{{G^{2}}}}}\right)\cos\left({2g}\right).$
(4)

The angular variable $h$ is a cyclic variable, showing that its conjugate
momentum $H=\sqrt{1-e^{2}}\cos{i}$ becomes a motion integral under the
quadrupole-order approximation. Naturally, the dynamical model represented by
Eq. (4) is of one degree of freedom. Similar to Kozai (1962), let us denote
the motion integral by

$H=\sqrt{1-e^{2}}\cos{i}=\cos{i_{*}},$ (5)

where $i_{*}$ stands for the critical inclination (known as Kozai parameter)
when the eccentricity is assumed at zero.

The equilibrium points should satisfy the following stationary condition:

$\displaystyle\dot{g}$
$\displaystyle=-G+5\frac{{{H^{2}}}}{{{G^{3}}}}+5\left({G-\frac{{{H^{2}}}}{{{G^{3}}}}}\right)\cos\left({2g}\right)=0,$
(6) $\displaystyle\dot{G}$
$\displaystyle=-5\left({1-{G^{2}}+{H^{2}}-\frac{{{H^{2}}}}{{{G^{2}}}}}\right)\sin\left({2g}\right)=0.$

The second equality of Eq. (6) indicates that the equilibrium points are
located at $2g=0$ or $2g=\pi$. It is not difficult to demonstrate that the
equilibrium points at $2g=\pi$ are stable in dynamics (because the Hamiltonian
takes its maximum value when $2g=\pi$), showing that Kozai resonance occurs at
$2g=\pi$.

At the Kozai centre ($2g=\pi$), the first equality of Eq. (6) yields

$G^{4}=\frac{5}{3}H^{2},$ (7)

which leads to the following relation between eccentricity and inclination at
Kozai centres (Kozai 1962),

$\cos{i}=\pm\sqrt{\frac{3}{5}(1-e^{2})}.$ (8)

Eq. (8) determines two curves of Kozai centre in the eccentricity–inclination
space, one is in the prograde region and the other one is in the retrograde
region (see Fig. 4). Furthermore, the eccentricity of Kozai centre can be
expressed as a function of the Kozai parameter $i_{*}$,

$e=\sqrt{1-\sqrt{\frac{5}{3}\cos^{2}{i_{*}}}},$ (9)

and the inclination of Kozai centre is expressed as a function of $i_{*}$,

$\quad\cos{i}=\sqrt{\frac{3}{5}}\cos{i_{*}}.$ (10)

Physical solution of Eq. (9) implies that the Kozai parameter should satisfy
$\cos^{2}{i_{*}}<3/5$ (i.e., $\left|H\right|<\sqrt{3/5}$ or
$39.23^{\circ}<i_{*}<140.77^{\circ}$). Thus, the occurrence of Kozai resonance
requires $39.23^{\circ}<i_{*}<140.77^{\circ}$ (or
$\left|H\right|<\sqrt{3/5}$), meaning that Kozai resonance disappears when
$i_{*}$ is smaller than $39.23^{\circ}$ or greater than $140.77^{\circ}$.

In Fig. 2, we report the phase portraits (i.e., level curves of Hamiltonian)
in the space $(e\cos{\omega},e\sin{\omega})$ for the motion integral at
$\left|H\right|=\sqrt{3}/2$ (corresponding to the Kozai parameter
$i_{*}=30^{\circ}$ or $i_{*}=150^{\circ}$) and at $\left|H\right|=0.5$
(corresponding to the Kozai parameter $i_{*}=60^{\circ}$ or
$i_{*}=120^{\circ}$). In the case of $\left|H\right|=\sqrt{3}/2$ (see the left
panel of Fig. 2), Kozai resonance cannot happen and all the motions in the
phase space are of circulation. In the case of $\left|H\right|=0.5$ (see the
right panel of Fig. 2), the Kozai resonance takes place around $2\omega=\pi$
and the regions of libration and circulation are divided by the dynamical
separatrix.

Figure 2: Phase portraits of the dynamics at the quadrupole order specified by
the motion integral $\left|H\right|=\sqrt{3}/2$ corresponding to the Kozai
parameter $i_{*}=30^{\circ}$ or $i_{*}=150^{\circ}$ (_left panel_) and
$\left|H\right|=0.5$ corresponding to the Kozai parameter $i_{*}=60^{\circ}$
or $i_{*}=120^{\circ}$ (_right panel_). The zero-eccentricity point (i.e., the
coordinate origin) is marked by a black star. In the _left panel_ , all
motions are of circulation. In the _right panel_ , there are three equilibrium
points and the regions of libration and circulation are divided by the
dynamical separatrix (shown by the red line). The separatrix corresponds to
the level curve of ${\cal H}=-H^{2}=-\cos^{2}{i_{*}}$.

For the one-degree-of-freedom dynamical model at hand, the Hamiltonian ${\cal
H}$ and the motion integral $H$ (or the Kozai parameter $i_{*}$) are two
conserved quantities and their magnitudes determine the mode of motion
(libration or circulation). Next, the libration zones of Kozai resonance are
discussed.

On one hand, let us discuss the circulation and libration zones of Kozai
resonance in the space spanned by ${\cal H}$ and $H$ (or $i_{*}$). According
to the Hamiltonian given by Eq. (4), we can get that the minimum Hamiltonian
happens at $G=\left|H\right|$ (i.e., the inclination $i$ is equal to
$0^{\circ}$ or $180^{\circ}$) and $2g=0$,

${\cal H}=-2+H^{2}=-2+\cos^{2}{i_{*}},$ (11)

which provides a lower boundary of the allowed Hamiltonian ${\cal H}$ as a
function of $H$ (or $i_{*}$). According to resonant theory, the Hamiltonian
takes its maximum value at the resonant centre (stable equilibrium point).
Accordingly, we get that the maximum Hamiltonian happens at the Kozai centre
(i.e., $2g=\pi$ and $G^{4}=5/3H^{2}$),

${\cal
H}=3+4H^{2}-2\sqrt{15}\left|H\right|=3+4\cos^{2}{i_{*}}-2\sqrt{15}\left|\cos{i_{*}}\right|,$
(12)

which provides the upper boundary of the allowed Hamiltonian ${\cal H}$ as a
function of $H$ (or $i_{*}$). Additionally, it is observed from the phase
portrait (see the right panel of Fig. 2) that the separatrix passes through
the zero-eccentricity point. Thus, the dynamical separatrix holds the
Hamiltonian

${\cal H}=-H^{2}=-\cos^{2}{i_{*}},$ (13)

which provides the boundary between the libration and circulation zones in the
space $({\cal H},H)$ or $({\cal H},i_{*})$.

Figure 3: Distribution of libration and circulation domains in the space
spanned by Hamiltonian ${\cal H}$ and Kozai parameter $i_{*}$. The critical
lines of $i_{*}=39.23^{\circ},140.77^{\circ}$ are marked.

In Fig. 3, the lower and upper boundaries of the allowed regions and the
separatrix dividing libration and circulation zones are reported in the space
$({\cal H},i_{*})$. The allowed regions are bounded by the lines corresponding
to ${\cal H}=\cos^{2}{i_{*}}-2$ and ${\cal
H}=3+4\cos^{2}{i_{*}}-2\sqrt{15}\left|\cos{i_{*}}\right|$ and the white areas
in the space are forbidden regions. The circulation and libration zones are
separated by the critical line of ${\cal H}=-\cos^{2}{i_{*}}$ (called
separatrix). From Fig. 3, it is observed that the occurrence of Kozai
resonance requires that the Kozai parameter and Hamiltonian should satisfy the
following conditions: $39.23^{\circ}<i_{*}<140.77^{\circ}$ and $-0.6<{\cal
H}<3$.

About the distribution of libration and circulation zones, similar plot shown
in the space $({\cal H},\left|H\right|)$ has been obtained in Sidorenko
(2018), who derived the boundaries by introducing auxiliary functions.

Figure 4: Distribution of libration and circulation domains in the
eccentricity–inclination space. The level curves of the motion integral $H$
are plotted and the polar line is marked by a vertical blue dashed line. The
separatrices corresponding to $i=39.23^{\circ}$ and $i=140.77^{\circ}$ are
shown by two vertical black lines.

On the other hand, let us discuss the distribution of libration and
circulation zones of Kozai resonance in the eccentricity–inclination space. As
shown in the right panel of Fig. 2, the width of Kozai resonance can be
measured at $2\omega=180^{\circ}$ (location of resonant centre) and it is
denoted by the variation of eccentricity $\Delta e=e_{2}-e_{1}$ (the resonant
width in terms of inclination variation can be naturally denoted by $\Delta
i=i_{2}-i_{1}$).

Along a certain isoline specified by $H$ or $i_{*}$, we can see that the lower
boundary of eccentricity reaches $e_{1}=0$ where the inclination is equal to
the Kozai parameter $i_{1}=i_{*}$ (see Fig. 4). About the eccentricity $e_{2}$
at the upper boundary (the action variable is denoted by
$G_{2}=\sqrt{1-e_{2}^{2}}$), it satisfies the following equality (see the
phase portrait shown in the right panel of Fig. 2),

${\cal H}\left(H=\cos{i_{*}};2g=\pi,G=G_{2}\right)={\cal
H}\left(H=\cos{i_{*}};2g=0,G=1\right),$

which leads to the eccentricity at the upper boundary (Naoz 2016),

$e_{2}=\sqrt{1-\frac{5}{3}\cos^{2}{i_{*}}}.$ (14)

The motion integral $H=\sqrt{1-e_{2}^{2}}\cos{i_{2}}=\cos{i_{*}}$ leads to the
inclination at the upper boundary (Naoz 2016),

$\cos{i_{2}}=\pm\sqrt{\frac{3}{5}}$ (15)

which gives $i_{2}\approx 39.23^{\circ}$ or $i_{2}\approx 140.77^{\circ}$
(called Kozai angles). This means that the boundary of inclination ($i_{2}$)
for Kozai librations is independent on the Kozai parameter $i_{*}$. In Fig. 4,
the separatrices corresponding to $i=39.23^{\circ}$ and $i=140.77^{\circ}$ are
shown in black lines. In addition, the conservation of $H$ means that the
inclination of test particles cannot go over the limit of $i=90^{\circ}$, so
that the Kozai oscillations are confined to the prograde region with
inclination ranging from $39.23^{\circ}$ to $90^{\circ}$ or to retrograde
region with inclination ranging from $90^{\circ}$ to $140.77^{\circ}$.

Evidently, the resonant width in terms of the variation of eccentricity is
$\Delta e=e_{2}-e_{1}=\sqrt{1-\frac{5}{3}\cos^{2}{i_{*}}}$ and the one in
terms of the variation of inclination is $\Delta
i=i_{2}-i_{1}=i_{*}-39.23^{\circ}$ for the prograde case and $\Delta
i=i_{2}-i_{1}=140.77^{\circ}-i_{*}$ for the retrograde case. Please see Fig. 5
for $\Delta e$ and $\Delta i$ as functions of the Kozai parameter $i_{*}$. It
is observed that, when the parameter $i_{*}$ is closer to $90^{\circ}$, the
resonant width ($\Delta e$ or $\Delta i$) is larger.

Figure 5: Libration width in terms of $\Delta e$ and $\Delta i$ for Kozai
resonance under the dynamical model at the quadrupole-level approximation. The
vertical green dashed line stands for the location of $i_{*}=90^{\circ}$.

In summary, at the quadrupole-level approximation ($\epsilon=0$), the
dynamical model is totally integrable. In the phase space with Kozai parameter
$39.23^{\circ}<i_{*}<140.77^{\circ}$, there is only Kozai resonance and the
resonant centre is located at ($2g=\pi$,
$G^{2}=\sqrt{5/3}\left|\cos{i_{*}}\right|$). For such an integrable model, the
modes of motion are determined by the Hamiltonian ${\cal H}$ and motion
integral $H$ (or $i_{*}$) and the separatrix between libration and circulation
regions is given by the critical line of ${\cal H}=-H^{2}=-\cos^{2}{i_{*}}$
(see Fig. 3). Evaluating the resonant width of Kozai island, we can further
get that the libration zone of Kozai resonance covers the inclination space
ranging from $i=39.23^{\circ}$ to $90^{\circ}$ (oscillations in the prograde
region) and from $i=90^{\circ}$ to $140.77^{\circ}$ (oscillations in the
retrograde region), and the resonant width ($\Delta e$ or $\Delta i$) is
larger when $i_{*}$ is closer to $90^{\circ}$, as shown in Figs. 4 and 5.

##  3 Poincaré surfaces of section

When the parameter $\epsilon$ is different from zero, the dynamical model is a
two-degree-of-freedom system but with one conserved quantity (i.e., the
Hamiltonian), leading to the fact that, in general, the resulting dynamical
model is not integrable. For such a kind of two degree of freedom model,
Poincaré surface of section is a powerful tool to uncover the global dynamics
in the phase space (Lithwick and Naoz 2011; Li et al 2014). Thus, in this
section, we will analyse the surfaces of section at different levels of
Hamiltonian to get some preliminary dynamics.

In practice, we define two types of sections. In the first type, the section
is defined as

$\Omega=\pi/2,\quad{\dot{\Omega}}{{\dot{\Omega}}_{0}}>0,$

and, in the second type, the section is defined by

$\Omega=\pi,\quad{\dot{\Omega}}{{\dot{\Omega}}_{0}}>0.$

The condition ${\dot{\Omega}}{{\dot{\Omega}}_{0}}>0$ implies that the test
particle crosses the section of $\Omega=\pi/2$ (or $\Omega=\pi$) in the same
direction (${\dot{\Omega}}_{0}$ stands for the initial time derivative of the
longitude of ascending node).

It should be mentioned that our choice of section at $\Omega=\pi/2$ (the first
type) and $\Omega=\pi$ (the second type) is based on the consideration that
there is a correspondence between stable periodic orbits and centres of
libration island arising in the Poincaré sections and thus the initial
condition of stable periodic orbits can be provided by locating the centre of
libration islands. In addition, it is better to place the starting points of
symmetric periodic orbits at the aligned apsides with $\varpi_{0}=0$ or at the
anti-aligned apsides with $\varpi_{0}=0$. As a result, the first-type section
defined by $\Omega=\pi/2$ aims to produce periodic orbits associated with the
resonances centred at $\omega=\pi/2$ or $\omega=3\pi/2$ (for example, Kozai
resonance) and the second-type section defined by $\Omega=\pi$ aims to
generate those periodic orbits associated with the resonances centred at
$\omega=\pi$ or $\omega=0$ (for example, the apsidal resonance).

The points on the defined sections are presented in the
$(e\cos\omega,e\sin\omega)$ plane, in which it is clearly observed that the
regular regions are filled with smooth curves and chaotic regions are filled
by scattered points. In the regular regions, there are resonant and
circulatory trajectories and, in the regions filled with scattered points, the
motions are chaotic (Li et al 2014). On the Poincaré sections, at each point
specified by $(e\cos\omega,e\sin\omega)$, the variable $\dot{\Omega}$ can be
retrieved from the given Hamiltonian ${\cal H}$. Thus, the Poincaré sections
could provide the full state information for those trajectories crossing the
defined sections for a given Hamiltonian value. For the simulations performed
in this section, we take the factor $\epsilon=0.03$ as an example to
characterise the dynamical model. The Poincaré sections of the first type
(defined by $\Omega=\pi/2$) are reported in Fig. 6 for different levels of
Hamiltonian.

Figure 6: Poincaré surfaces of section (defined by $\Omega=\pi/2$) for
different levels of Hamiltonian. The red star represents the zero-eccentricity
point. Note that all points in the Poincaré section are prograde.

When the Hamiltonian is ${\cal H}=-0.43$ (see panel ‘a’), there are two
islands of libration centred at $\omega=90^{\circ}$ and $\omega=270^{\circ}$.
It is observed that these two islands are surrounded by chaotic regions, and
the chaotic regions are distributed along the separatrix of the Kozai
resonance shown in Fig. 2. This phenomenon can be understood as follows: the
perturbation transforms the separatrix into a chaotic layer, and such a
chaotic layer becomes larger and larger when the perturbation increases.

In the centre of each island, there is a (stable) periodic orbit corresponding
to the resonant centre with $\omega$ at $90^{\circ}$ and $270^{\circ}$ (these
two islands of libration are due to the well-known Kozai resonance), and
inside the island the motions are quasi-periodic trajectories (also known as
resonant trajectories). The size of the libration island, marked by $\Delta
e$, can measure the range of the stable libration zone of Kozai resonance,
which is to be discussed in Sect. 5 for more details.

When the Hamiltonian is decreased to ${\cal H}=-0.5$, the structures of
Poincaré sections are significantly changed (see panel ‘b’). Besides the
islands appearing in panel (a), an additional island of libration arises at
$\omega=90^{\circ}$ in the region with higher eccentricities. Totally, there
are three islands of libration in the considered section, two of them
correspond to the Kozai resonance with argument of $\sigma=\omega$ (the
resonant centres are located at $\sigma=90^{\circ}$ and $\sigma=270^{\circ}$),
and the third island with a higher eccentricity is due to the (anti-aligned)
apsidal secular resonance with critical argument $\sigma=\omega+\Omega$
librating around $\sigma=\pi$.

When the Hamiltonian is further decreased to ${\cal H}=-0.63$ (see panel ‘c’)
and ${\cal H}=-0.88$ (see panel ‘d’), the dynamical structures arising in the
Poincaré section become relatively simple. In both panels, there is only one
island of libration in the section, which corresponds to the (aligned) apsidal
secular resonance with $\sigma=\omega+\Omega$ librating around $\sigma=0$. The
size of the libration island is also marked by $\Delta e$, which could measure
the range of stable libration zone, as shown in panel (d).

From Fig. 6, we can see that there exists Kozai resonance in the Poincaré
section when the Hamiltonian is equal to ${\cal H}=-0.43,-0.5$ and there is no
Kozai resonance when the Hamiltonian is equal to ${\cal H}=-0.63,-0.88$. This
is in agreement with the discussion made in Section 2.2 that the occurrence of
Kozai resonance requires the Hamiltonian $-0.6<{\cal H}<3$. However, the
perturbation due to the octupole-order terms in the Hamiltonian introduces
significant changes about the dynamical structures. To show this point, here
we make a direct comparison between the phase portraits produced under the
quadrupole-level approximation (see Fig. 2) and the Poincaré surfaces of
section obtained from the octupole-order approximation (see Fig. 6). Firstly,
the comparison between panel (a) of Fig. 2 and panels (c,d) of Fig. 6 shows
that there are no any resonance under the quadrupole-order model but, under
the octupole-order model, the apsidal secular resonance with critical argument
$\sigma=\omega+\Omega$ appears. Secondly, the comparison between panel (b) of
Fig. 2 and panels (a,b) of Fig. 6 indicates that (a) the separatrix appearing
in the phase portrait under the quadrupole-order model is replaced by a
chaotic layer under the octupole-order model and (b) besides the Kozai
resonance the apsidal secular resonance with critical argument of
$\sigma=\omega+\Omega$ takes place under the octupole-order model.

Figure 7: Poincaré surfaces of section (defined by $\Omega=\pi$) for two
levels of Hamiltonian at ${\cal H}=-0.43$ and ${\cal H}=-0.5$. The red dot
stands for the zero-eccentricity point. Note that all points in the Poincaré
section are on prograde orbits.

In Fig. 7, the Poincaré sections of the second type (defined by $\Omega=\pi$)
are reported. The Hamiltonian is taken as ${\cal H}=-0.43$ in the left-hand
panel and ${\cal H}=-0.5$ in the right-hand panel. In each panel, it is
observed that there are three islands of libration. The first two islands
centred at $\omega=90^{\circ}$ and $\omega=270^{\circ}$ corresponds to the
Kozai resonance, which have been appeared in the first two panels of Fig. 6.
The third island of libration centred at $\omega=0^{\circ}$ is a new one in
comparison to the ones appearing in Fig. 6 and it is due to the (anti-aligned)
apsidal secular resonance with $\sigma=\omega+\Omega$ librating around
$\sigma=\pi$.

At the centre of each island of libration, there is a (stable) periodic orbit.
According to dynamical system theory, the (stable) periodic orbit projected
onto the Poincaré section corresponds to the resonant centre, and the quasi-
periodic orbits correspond to the associated librating trajectories inside the
islands (Lithwick and Naoz 2011; Li et al 2014; Naoz 2016). Thus, we could
compute stable periodic orbits in order to identify the location of resonant
centres and then evaluate the size of libration islands filled with quasi-
periodic orbits in order to specify the boundaries of stable libration zones.
This is the key idea of the current work.

In the current study, we concentrate on the symmetric periodic orbits with
starting points on the aligned or anti-aligned apsides (i.e., $\varpi_{0}=0$
or $\varpi_{0}=\pi$). Those asymmetric periodic orbits and their associated
resonances are out of the scope of the current study. According to the
Poincaré sections shown in Fig. 6, we can summarise that the initial angles of
the first two types of symmetric periodic orbits are
$\omega_{0}=3\pi/2,\Omega_{0}=\pi/2$ (case I) and
$\omega_{0}=\pi/2,\Omega_{0}=\pi/2$ (case II), and the sections shown in Fig.
7 indicate that the initial angles of the third type of symmetric periodic
orbits are $\omega_{0}=0$ and $\Omega_{0}=\pi$ (case III).

It should be mentioned that this study is restricted to those major secular
resonances that are visible in the Poincaré surfaces of section, and thus they
only represent a subset of varieties of secular resonances. This choice is
justified by the fact that these resonances considered are the strongest ones,
dominating the secular dynamics of test particles.

In the coming section, we determine the initial eccentricity and inclination
for symmetric periodic orbits of cases I, II and III (their initial angles
$\Omega_{0}$ and $\omega_{0}$ are known).

##  4 Webs of secular resonances (or families of stable periodic orbits)

As discussed in the previous section, the appearance of libration islands in
the Poincaré sections is due to secular resonances and there is a
correspondence between stable periodic orbits and resonant centres. Thus, we
could determine the exact location of resonance by computing the associated
(stable) periodic orbit.

In the two-degree-of-freedom system, a periodic orbit should satisfy the
following constraints:

${\left\\{{k,h,q,p}\right\\}_{t=0}}={\left\\{{k,h,q,p}\right\\}_{t=T}},$

where $T$ is the period. Obviously, the topic of solving periodic orbits is a
typical two-point boundary-value problem (TPBVP). In general, for such a
TPBVP, we need to determine five parameters including the initial states
$(k_{0},h_{0},q_{0},p_{0})$ and period $T$.

Let us take case I as an example. The initial states of periodic orbits (for
case I) can be expressed as follows (it is known that the initial angles are
$\omega_{0}=3\pi/2$ and $\Omega_{0}=\pi/2$):

$k_{0}=0,\quad h_{0}=e_{0},\quad q_{0}=0,\quad p_{0}=\sin{\frac{i_{0}}{2}},$

where only the initial eccentricity and inclination ($e_{0}$ and $i_{0}$) are
undetermined variables. In addition, for a given Hamiltonian ${\cal H}$, only
one of the variables $e_{0}$ and $i_{0}$ is independent, i.e., $e_{0}$ (or
$i_{0}$) can be known from the expression of Hamiltonian ${\cal
H}(e_{0},i_{0},\Omega_{0}=3\pi/2,\omega_{0}=\pi/2)$. In summary, for a given
Hamiltonian, we only need to identify two variables including $e_{0}$ (or
$i_{0}$) and $T$ by means of correction algorithm. The Hamiltonian ${\cal H}$
varies in a certain interval, then the entire family of periodic orbits can be
obtained by means of continuation algorithm. Please refer to Howell (1984) and
Lara and Peláez (2002) for the correction and continuation algorithms.

To judge the stability of a periodic orbit, we need to know its Monodromy
matrix, which corresponds to the state transition matrix (STM) evaluated at
one period. For the two-degree of freedom dynamical model considered in the
present study, there are two pairs of eigenvalues for the Monodromy matrix.
One pair of them is equal to unity, and the other non-unity pair of
eigenvalues determine the stability. Let us denote the non-unity pair of
eigenvalues as $(\lambda_{1},\lambda_{2})$ and define the stability index as
$S=\lambda_{1}+\lambda_{2}$. The stable condition of periodic orbits states
that the stability index $S$ should be smaller than 2, i.e., $S<2$ (Hénon
1969).

In Fig. 8, the characteristic curves of stable periodic orbits ($S<2$) are
presented in the action space $(e_{0},i_{0})$ for cases I, II and III. In case
I (see panel a) and case II (see panel b), there are eight families of
periodic orbits in the whole action space, and they are denoted by $f_{i}$ for
case I and $g_{i}$ for case II ($i=1,2,...,8$). In case III (see panel c),
there are two families of stable periodic orbits, which are denoted by $h_{1}$
for the prograde case and $h_{2}$ for the retrograde case. In a certain family
of periodic orbits, all the members hold the same resonant argument. Please
refer to Table 1 for the resonant arguments and their libration centres.

Figure 8: Characteristic curves of families of stable periodic orbits (or webs of secular resonances) in the case of $\epsilon=0.03$ for cases I, II and III. Table 1: Resonant arguments for the families of stable periodic orbits shown in Fig. 8. $\sigma_{i}$ stands for the critical argument and $\sigma_{i,c}$ is the associated resonant centre. Family | Resonances | $\sigma_{i}$ | $\sigma_{i,c}$
---|---|---|---
$f_{1}(g_{1})$ | Aligned apsidal secular resonance |  $\omega+\Omega$ ($\omega-\Omega$) | $0^{\circ}$
$f_{2}(g_{2})$ | Aligned apsidal secular resonance |  $\omega+\Omega$ ($\omega-\Omega$) | $0^{\circ}$
$f_{3}(g_{3})$ | Kozai resonance | $\omega$ |  $270^{\circ}$ ($90^{\circ}$)
$f_{4}(g_{4})$ | Kozai resonance | $\omega$ |  $270^{\circ}$ ($90^{\circ}$)
$f_{5}(g_{5})$ | Kozai resonance | $\omega$ |  $270^{\circ}$ ($90^{\circ}$)
$f_{6}(g_{5})$ | Kozai resonance | $\omega$ |  $270^{\circ}$ ($90^{\circ}$)
$f_{7}(g_{7})$ | High-order secular resonance |  $3\omega-\Omega$ ($3\omega+\Omega$) | $0^{\circ}$
$f_{8}(g_{8})$ | Anti-aligned apsidal secular resonance |  $\omega-\Omega$ ($\omega+\Omega$) | $180^{\circ}$
$h_{1}(h_{2})$ | Anti-aligned apsidal secular resonance |  $\omega+\Omega$ ($\omega-\Omega$) | $180^{\circ}$

Observing the characteristic curves shown in Fig. 8, we can find that (a) the
members in $f_{i}$ and the ones in $g_{i}$ (with the same $i$) are symmetric
with respect to the polar line and (b) the members in family $h_{1}$ and the
ones in $h_{2}$ are symmetric relative to the polar line. These symmetric
properties of solution are due to the symmetry inherent in the Hamiltonian, as
discussed in Sect. 2.

In panel (a) of Fig. 8, one representative member in each family is marked by
a green star. For convenience, these representatives are denoted by letters
from ‘A’ in family $f_{1}$ to ‘H’ in family $f_{8}$. Their time histories of
resonant argument during one resonant period are shown in Fig. 9. It is
observed from Fig. 9 that (a) for members ‘A’ in family $f_{1}$ and ‘B’ in
$f_{2}$, the resonant angle $\sigma_{1,2}=\omega+\Omega$ librates around
$0^{\circ}$, (b) for members ‘C’, ‘D’, ‘E’ and ‘F’ in families $f_{i}$ with
$i=3,4,5,6$, the resonant angle $\sigma_{i}=\omega$ librates around
$270^{\circ}$, (c) for member ‘G’ in family $f_{7}$, the resonant angle
$\sigma_{7}=3\omega-\Omega$ librates around $0^{\circ}$ and (d) for member ‘H’
in family $f_{8}$, the argument $\sigma_{8}=\omega-\Omega$ librates around
$180^{\circ}$.

Figure 9: Time histories of the resonant angle for the example periodic orbits
marked in panel (a) of Fig. 8.

From the viewpoint of geometrical configurations, we can see that the families
of stable periodic orbits shown in Fig. 8 belong to several types of secular
resonances (see Table 1): (a) families $f_{1,2}(g_{1,2})$ correspond to the
aligned apsidal secular resonances with critical argument librating around
zero, (b) families $f_{i}(g_{i})$ with $i=3,4,5,6$ correspond to Kozai
resonances with $\omega$ as the critical argument librating around
$90^{\circ}$ or $270^{\circ}$, (c) families $f_{7}(g_{7})$ correspond to high-
order secular resonances and (d) families $f_{8}(g_{8})$ and $h_{1}(h_{2})$
correspond to the anti-aligned apsidal secular resonances with critical
argument librating around $180^{\circ}$.

Figure 10: Characteristic curves of stable periodic orbits (i.e., webs of
secular resonances) obtained under dynamical models specified by different
values of $\epsilon$. The upper panel is for case I and the bottom panel is
for case III. The results in case II can be easily obtained by symmetry.

Following the similar approach, we can further compute families of stable
periodic orbits under the dynamical models characterised by different values
of $\epsilon$. The characteristic curves (or webs of secular resonances) of
those primary families are reported in Fig. 10. It is noted that families of
stable periodic orbits in case II are not considered due to symmetry and
families $f_{3}$, $f_{6}$ and $f_{7}$ are not presented considering their weak
strength of resonance (their islands are not visible from Poincaré sections).
For convenience of comparison, characteristic curves corresponding to
$\epsilon=0$ (see Fig. 3) and the ones corresponding to $\epsilon=0.03$ (see
Fig. 8) are also presented in Fig. 10.

Observing Fig. 10, we can see that (a) in the case of $\epsilon=0$
(quadrupole-level approximation) there are only two families of stable
periodic orbits corresponding to Kozai resonances (one is in the prograde
region and the other one is in the retrograde region) and other families begin
to appear when $\epsilon$ is different from zero, (b) characteristic curves
are changed when $\epsilon$ is different, meaning that the location of
resonant centres is dependent on the magnitude of $\epsilon$, (c) when
$\epsilon$ is different from zero the characteristic curves of the families of
Kozai resonance (i.e., families $f_{4}$ and $f_{5}$) can not extend to the
region around the zero-eccentricity point because the low-eccentricity regions
are chaotic in dynamics, (d) in the case of $\epsilon=0$ (quadrupole-level
approximation) the curves of Kozai centres are symmetric with respect to
$i_{0}=90^{\circ}$, while such a symmetry breaks when $\epsilon$ is different
from zero, (e) the characteristic curve of $f_{8}$ is insensitive to the
eccentricity, (f) characteristic curves of $h_{1}$ and $h_{2}$ are always
symmetric with respect to $i_{0}=90^{\circ}$ at arbitrary $\epsilon$ and (g)
secular resonances in families $h_{1}$ and $h_{2}$ occur in the region with
eccentricities smaller than $\sim$0.35.

In the next section, we will measure the boundaries of stable libration zones
in the action space by analysing Poincaré sections.

##  5 Boundaries of stable libration zones

In the previous section, we have produced the distribution of exact resonant
centres in the action space $(e_{0},i_{0})$ by determining families of
symmetric periodic orbits. In addition, we know that the islands of libration
appearing in the Poincaré sections are due to the secular resonances, thus the
boundaries of stable libration zones can be measured by analysing Poincaré
sections.

In this work, we evaluate boundaries of stable libration zones when the angles
are fixed at the resonant centre, i.e. fixing $\Omega_{0}=\pi/2$ and
$\omega_{0}=3\pi/2$ for case I (corresponding to families $f_{i}$ with
$i=1,2,...,8$), $\Omega_{0}=\pi/2$ and $\omega_{0}=\pi/2$ for case II
(corresponding to families $g_{i}$ with $i=1,2,...,8$), and $\Omega_{0}=\pi$
and $\omega_{0}=0$ for case III (corresponding to families $h_{1}$ and
$h_{2}$). It should be noted that the results (boundaries of stable libration
zones) will be changed if we adopt different definitions of Poincaré sections.

To describe the process of evaluating libration zones in the space
$(e_{0},i_{0})$, let us take case I ($\Omega_{0}=\pi/2$ and
$\omega_{0}=3\pi/2$) as an example. A procedure is designed to identify the
boundaries of stable libration zones by searching for the upper and lower
limits that enclose the islands of libration as follows:

  1. 1.

For a given Hamiltonian ${\cal H}_{0}$, the eccentricity and inclination
($e_{0}$ and $i_{0}$) of the exact resonant centre (i.e. the periodic orbit)
can be determined, as shown in Figs 8 and 10.

  2. 2.

Starting from $e_{0}$, increase the eccentricity and keep the Hamiltonian as a
constant equal to ${\cal H}_{0}$. When the trajectory switches between
libration and circulation, we halve the step of eccentricity until the step is
smaller than a given tolerance. Let us record the maximum eccentricity at
which the trajectory is of libration as $e_{\max}$.

  3. 3.

Similar to the previous step, decreasing the eccentricity from $e_{0}$ and
keeping the Hamiltonian as a constant equal to ${\cal H}_{0}$, we can
determine the minimum eccentricity at which the trajectory is of libration.
Record the minimum value of eccentricity at which the trajectory is of
libration as $e_{\min}$.

  4. 4.

According to the Hamiltonian ${\cal H}_{0}$, it is possible to reproduce the
boundary points in the action space as $(e_{\min},i_{1})$ and
$(e_{\max},i_{2})$.

  5. 5.

When the Hamiltonian varies, the points $(e_{\min},i_{1})$ and
$(e_{\max},i_{2})$ constitute the boundaries of stable libration zones.

In practical simulations, we limit the duration of numerical integration as
1000 in normalised units, which is much longer than the timescales of secular
resonances considered in the present work. In the regular regions, there is no
problem to distinguish the points with libration and circulation. However, in
the chaotic region, it becomes complicated and the following rules are adopted
to distinguish orbits of libration and circulation: (a) those trajectories are
recorded as librating orbits if their critical arguments are always of
libration during the maximum period of time, and (b) those trajectories are
recorded as non-libration if their critical arguments enter into the
circulation mode in the process of numerical integration.

Figure 11: Stable libration zones under dynamical model characterised by
$\epsilon=0.03$. The dashed lines specified by $i_{0}=39.2^{\circ}$ and
$i_{0}=140.8^{\circ}$ stand for the separatrices of Kozai zones in the case of
$\epsilon=0$.

In the case of $\epsilon=0.03$, the results are reported in Fig. 11, where the
stable libration zones are shown in the space of $(e,i)$ as shaded areas for
three cases of symmetric periodic orbits with angles at $\Omega_{0}=\pi/2$ and
$\omega_{0}=3\pi/2$ (see panel a), $\Omega_{0}=\pi/2$ and $\omega_{0}=\pi/2$
(see panel b) and $\Omega_{0}=\pi$ and $\omega_{0}=0$ (see panel c). The
characteristic curves standing for the location of resonant centre are also
shown in each panel. For convenience of comparison, the boundaries of Kozai
zones in the case of $\epsilon=0$ are shown in dashed lines, which are
specified by $i_{0}=39.2^{\circ}$ in the prograde region and
$i_{0}=140.8^{\circ}$ in the retrograde region.

From Fig. 11, we can observe that (a) the secular resonances associated with
$f_{3}(g_{3})$, $f_{6}(g_{6})$ and $f_{7}(g_{7})$ happen in the vicinities of
$i_{0}=39.2^{\circ}$ and $i_{0}=140.8^{\circ}$ and they have very small
domains in the considered space, (b) the secular resonances associated with
$f_{1,2}$ ($g_{1,2}$) occur in the low-eccentricity region with inclinations
smaller than $i_{0}\approx 39^{\circ}$ (greater than $i_{0}\approx
141^{\circ}$), (c) there is a large Kozai zone but it cannot cover the low-
eccentricity region because the area near the zero-eccentricity point is
filled with chaotic motion, as shown by the Poincaré sections in Figs. 6 and
7, (d) in the high-eccentricity region the boundaries of Kozai zones move
rightward (leftward) for case I (case II) in comparison to the ones in the
case of $\epsilon=0$, (e) the resonances associated with $h_{1}$ and $h_{2}$
occurs in the mid- and low-eccentricity region (these two resonances do not
exist in the high-eccentricity space), and (f) the libration zones associated
with $h_{1}$ and $h_{2}$ are symmetric with respect to the polar line.

Figure 12: Stable libration zones under dynamical models characterised by
$\epsilon=0.05,0.07,0.1$. The dashed lines specified by $i_{0}=39.2^{\circ}$
and $i_{0}=140.8^{\circ}$ stand for the separatrices of Kozai zones in the
case of $\epsilon=0$. The panels in the left column correspond to case I and
the ones in the rigth column stands for case III (results about secular
resonances in case II can be easily produced by means of symmetry).

Extending to higher values of $\epsilon$, we produce the boundaries of stable
libration zones in the $(i_{0},e_{0})$ space, as shown in Fig. 12 where the
dashed lines stand for the boundaries of Kozai zone in the case of
$\epsilon=0$. In practice, the cases of $\epsilon=0.05,0.07,0.1$ are taken
into account, and families $f_{3}$, $f_{6}$ and $f_{7}$ are not considered due
to their small strength. From Fig. 12, we can see that (a) the topological
structures of libration zones for the primary resonances remain similar to the
ones shown in Fig. 10, (b) when $\epsilon$ increases, the area of Kozai
libration zone becomes smaller and smaller, showing that the gap in low-
eccentricity region becomes larger and larger, (c) libration zones associated
with $f_{1}$ and $f_{2}$ have larger area with higher $\epsilon$, and (d) the
libration zones associated with $f_{8}$, $h_{1}$ and $h_{2}$ are insensitive
to the value of $\epsilon$.

Regarding the the apsidal secular resonances in families $f_{1}$ and $f_{2}$,
Lei (2021) analytically investigated their dynamics using the perturbation
theory based on Lie-series transformation and concluded that the considered
resonances govern long-term evolutions of test particles with inclinations
smaller than $39^{\circ}$. Comparing to the numerical results in the current
work, we can see a good agreement between the stable libration zones
associated with $f_{1}$ and $f_{2}$ (see Figs. 11 and 12) and the analytical
resonant width presented by Lei (2021) (see Fig. 6 in his work).

##  6 Conclusions

In this work, we systematically investigated varieties of secular resonances
for inner test particles in hierarchical planetary systems by means of non-
perturbative approaches. To formulate the dynamical model for long-term
evolution, the double-averaged Hamiltonian truncated at the third order in the
semimajor axis ratio (the octupole-level approximation) is adopted. The
resulting dynamical model is of two degrees of freedom with $\omega$ and
$\Omega$ as the angular coordinates.

In order to obtain preliminary dynamics, we produced Poincaré surfaces of
sections at different levels of Hamiltonian. In the surfaces of section, some
structures are observed: some areas are filled with chaotic motion, some areas
are filled with resonant trajectories (quasi-periodic orbits) and some areas
are filled with circulatory trajectories. According to dynamical system
theory, at the centre of each libration island (corresponding to the resonant
centre), there is a stable periodic orbit. In general, we summarised three
cases of symmetric periodic orbits as follows: case I with $\omega_{0}=3\pi/2$
and $\Omega_{0}=\pi/2$, case II with $\omega_{0}=\pi/2$ and $\Omega_{0}=\pi/2$
and case III with $\omega_{0}=0$ and $\Omega_{0}=\pi$.

Using numerical correction and continuation techniques, we computed families
of stable and symmetric periodic orbits. The characteristic curves of families
of stable periodic orbits are shown in the space $(e_{0},i_{0})$ and these
curves are also named webs of secular resonances. To measure the boundaries of
stable libration zones, we numerically determined the upper and lower limits
in the action space by analysing Poincaré surfaces of section. The parameter
spaces enclosed by the boundaries are called stable libration zones, inside
which the trajectories are of libration during the period of numerical
integration. The resonances considered in this study can be classified as
aligned or anti-aligned apsidal secular resonances, Kozai resonances and high-
order secular resonances. It should be mentioned that in this work only the
resonances which are visible in the Poincaré surfaces of section are taken
into account (those weak resonances are not considered). This is based on the
consideration that they are the strongest from the viewpoint of strength. In
addition, the boundaries of stable libration zones are obtained by analysing
Poincaré sections, so results will be changed if different definitions of
Poincaré section are adopted.

The effects of the octupole potential are evaluated by taking different values
of $\epsilon$ ($\epsilon$ specifies the contribution of octupole terms). Under
the quadrupole-order approximation (corresponding to $\epsilon=0$), the
dynamical model is totally integrable and there is only Kozai resonance, which
occurs in the inclination space ranging from $39.2^{\circ}$ to
$140.8^{\circ}$. Under the octupole-level approximation with non-zero
$\epsilon$, besides (perturbed) Kozai resonance, new secular resonances (such
as apsidal resonances and high-order resonances) appear in the phase space and
all these resonances constitute a web (called web of resonances). Due to the
perturbation of octupole-order potential, dynamics of Kozai resonance is
changed: Kozai resonance disappears in the low-eccentricity region because
this zone is filled with chaos. Results show that the chaotic area increases
with the factor $\epsilon$.

The results obtained in this work are helpful to understand the secular
dynamics of hierarchical exoplanetary systems in terms of the following
concerns: (a) how many kinds of secular resonances existing in the action
space, (b) the place where the exact resonance occurs, (c) the region where a
certain secular resonance can occur or not, and (d) the size of libration zone
changing with eccentricity and/or inclination for a particular resonance.

###### Acknowledgements.

The author thanks two anonymous reviewers for their insightful and rigorous
comments that help to improve the quality of this paper. This work is
performed with the financial support of the National Natural Science
Foundation of China (No. 12073011).

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * Ćuk and Burns (2004) Ćuk M, Burns JA (2004) On the secular behavior of irregular satellites. AJ 128(5):2518
  * Ford et al (2000) Ford EB, Kozinsky B, Rasio FA (2000) Secular evolution of hierarchical triple star systems. ApJ 535(1):385
  * Hamers and Samsing (2019) Hamers AS, Samsing J (2019) Analytic computation of the secular effects of encounters on a binary: features arising from second-order perturbation theory. MNRAS 487(4):5630–5648
  * Harrington (1968) Harrington RS (1968) Dynamical evolution of triple stars. AJ 73:190–194
  * Harrington (1969) Harrington RS (1969) The stellar three-body problem. Celest Mech 1(2):200–209
  * Hénon (1969) Hénon M (1969) Numerical exploration of the restricted problem, v. A&A 1:223–238
  * Howell (1984) Howell KC (1984) Three-dimensional, periodic,‘halo’orbits. Celest Mech 32(1):53–71
  * Ito and Ohtsuka (2019) Ito T, Ohtsuka K (2019) The lidov–kozai oscillation and hugo von zeipel. Monogr Environ Earth Planets 7:1–113
  * Katz et al (2011) Katz B, Dong S, Malhotra R (2011) Long-term cycling of kozai-lidov cycles: extreme eccentricities and inclinations excited by a distant eccentric perturber. Phys Rev Lett 107(18):181,101
  * Kozai (1962) Kozai Y (1962) Secular perturbations of asteroids with high inclination and eccentricity. AJ 67:591
  * Lara and Peláez (2002) Lara M, Peláez J (2002) On the numerical continuation of periodic orbits-an intrinsic, 3-dimensional, differential, predictor-corrector algorithm. A&A 389(2):692–701
  * Lee and Peale (2003) Lee MH, Peale S (2003) Secular evolution of hierarchical planetary systems. ApJ 592(2):1201
  * Lei (2019) Lei H (2019) A semi-analytical model for secular dynamics of test particles in hierarchical triple systems. MNRAS 490(4):4756–4769
  * Lei (2020) Lei H (2020) Dynamical models for secular evolution of navigation satellites. Astrodynamics 4:57–73
  * Lei (2021) Lei H (2021) Secular resonance of inner test particles in hierarchical planetary systems. MNRAS 506(2):1879
  * Lei et al (2018) Lei H, Circi C, Ortore E (2018) Modified double-averaged hamiltonian in hierarchical triple systems. MNRAS 481(4):4602–4620
  * Li et al (2014) Li G, Naoz S, Holman M, Loeb A (2014) Chaos in the test particle eccentric kozai-lidov mechanism. ApJ 791(2):86
  * Lidov (1962) Lidov M (1962) The evolution of orbits of artificial satellites of planets under the action of gravitational perturbations of external bodies. P&SS 9(10):719–759
  * Lithwick and Naoz (2011) Lithwick Y, Naoz S (2011) The eccentric kozai mechanism for a test particle. ApJ 742(2):94
  * Luo et al (2016) Luo L, Katz B, Dong S (2016) Double-averaging can fail to characterize the long-term evolution of lidov–kozai cycles and derivation of an analytical correction. MNRAS 458(3):3060–3074
  * Naoz (2016) Naoz S (2016) The eccentric kozai-lidov effect and its applications. ARA&A 54:441–489
  * Naoz et al (2011) Naoz S, Farr WM, Lithwick Y, Rasio FA, Teyssandier J (2011) Hot jupiters from secular planet–planet interactions. Nature 473(7346):187
  * Naoz et al (2013) Naoz S, Farr WM, Lithwick Y, Rasio FA, Teyssandier J (2013) Secular dynamics in hierarchical three-body systems. MNRAS 431(3):2155–2171
  * Sidorenko (2018) Sidorenko VV (2018) The eccentric kozai–lidov effect as a resonance phenomenon. Celest Mech Dynam Astron 130(1):4
  * Thomas and Morbidelli (1996) Thomas F, Morbidelli A (1996) The kozai resonance in the outer solar system and the dynamics of long-period comets. Celest Mech Dynam Astron 64(3):209–229


# A Transition System Abstraction Framework for Neural Network Dynamical
System Models

Yejiang Yang, Zihao Mo, Hoang-Dung Tran, and Weiming Xiang This research was
supported by the National Science Foundation, under NSF CAREER Award no.
2143351, NSF CNS Award no. 2223035, and NSF IIS Award nos. 2331937,
2331938.Yejiang Yang, Zihao Mo, and Weiming Xiang are with the School of
Computer and Cyber Sciences, Augusta University, Augusta GA 30912, USA. Hoang-
Dung Tran is with the School of Computing, University of Nebraska-Lincoln,
Lincoln, NE 68588, USA.Yejiang Yang is also with the School of Electrical
Engineering, Southwest Jiaotong University, China.

###### Abstract

This paper proposes a transition system abstraction framework for neural
network dynamical system models to enhance the model interpretability, with
applications to complex dynamical systems such as human behavior learning and
verification. To begin with, the localized working zone will be segmented into
multiple localized partitions under the data-driven Maximum Entropy (ME)
partitioning method. Then, the transition matrix will be obtained based on the
set-valued reachability analysis of neural networks. Finally, applications to
human handwriting dynamics learning and verification are given to validate our
proposed abstraction framework, which demonstrates the advantages of enhancing
the interpretability of the black-box model, i.e., our proposed framework is
able to abstract a data-driven neural network model into a transition system,
making the neural network model interpretable through verifying specifications
described in Computational Tree Logic (CTL) languages.

##  I Introduction

In recent years, the popularity of machine learning models has surged
dramatically, primarily owing to their remarkable capacity to offer valuable
insights and accurate predictions, such as the remaining useful life
prediction of lithium-ion batteries in [1], controlling physical human-robot
interaction in [2], modeling the thermal conductivity of water-based
nanofluids containing magnetic copper nanoparticles in [3]. Nonetheless,
conventional machine learning models are often regarded as “black boxes”,
making the control and verification of neural network models relying on real-
time trajectory computation with tools such as those proposed in [4], [5, 6,
7] and [8], etc., which are often considered computationally expensive and
challenging for real-time computation. Additionally, the dilemma of black-box
verification limits the applications of neural networks in dynamical systems
modeling, thus, the need for an interpretable abstraction for neural network
dynamical system that unveils the black box is critical.

Thanks to recent works on hybrid system learning in [9], [10], [11] and
research on learning frameworks via distributed learning structures in [12],
[13] and [14], etc., it is possible to break down one general system into
multiple sections and the relationships between them. However, most of the
existing works tend to concentrate on modeling hybrid dynamics using hybrid
models, inadvertently overlooking the essential benefits of employing
distributed learning models, which can significantly enhance the
interpretability [15] of the learning model while also broadening the scope of
model verification, e.g., verifying the learning model with user-specified
properties via Linear Temporal Logic [16].

Diverging from the approach of employing multiple neural networks to learn
hybrid systems, which primarily emphasizes subsystem learning, the challenge
of investigating inter-subsystem relationships is rooted in the abstraction of
the model. To begin with, since prior knowledge of a system is often limited,
especially in the context of data-driven modeling, we introduce a transition
system abstraction approach. This method involves partitioning the model’s
working zone using techniques founded on Maximum Entropy (ME) partitioning, as
discussed in [17]. Following this, the transitions between distinct partitions
will be approximated via reachability analysis, which enables us to abstract
the neural network dynamical system and the partition relationship into a
transition system. Consequently, this method provides an intuitive means of
discovering the relationship between subsections of a neural network and
contributes to its improved interpretability.

The main contributions of this paper can be summarized as follows.

  * •

The proposed framework involves the introduction of a novel transition system
abstraction for neural network dynamical systems. This model serves the dual
purpose of verifying the behaviors of localized subsystems and shedding light
on their interrelationships by estimating transitions within the neural
network dynamical system’s working zone through the application of
reachability analysis, a technique uniquely suited to analyzing localized
subsystem interactions.

  * •

The proposed neural transition system can be established on an existing neural
network dynamical system to make the black-box model interpretable. The
specifications can be verified through Computational Tree Logic (CTL), thereby
broadening the possibilities for neural-network-based model verification. The
approach is validated by human handwriting dynamics learning and verification.

The paper is structured as follows: Preliminaries and problem formulations are
given in Section II. The main result, a data-driven neural transition system
abstraction-based framework, is given in Section III. In Section IV, the
applications for abstracting the learning model of the LASA handwriting data
set with our proposed framework are presented. Conclusions are given in
Section V.

_Notations_ : In the rest of the paper, $\mathbb{N}$ denotes the natural
number sets, $\mathbb{R}$ is the field of real numbers; $\mathbb{R}^{n}$
stands for the vector space of $n$-tuples of real number; $\underline{X}$ and
$\overline{X}$ are the lower bound and upper bound of an interval $X$,
respectively.

##  II Preliminaries and Problem Formulation

A neural network dynamical system model describes the dynamical system
behaviors through an input-output point of view in the form of

$\displaystyle x(k+1)=\Phi(x(k),u(k))$ (1)

in which $\Phi:\mathbb{R}^{n_{x}+n_{u}}\to\mathbb{R}^{n_{x}}$ is the trained
neural network. A neural network dynamical system model is trained based on
samples, which can be well-trained under an adequate sample set. Therefore, we
assume that the $\Phi$ applies only to a localized working zone $\mathcal{X}$
which can be obtained based on samples.

Neural network $\Phi$ can approximate dynamical systems with high accuracy,
but the lack of interpretability burdens the processes of verification and
control synthesis. In this paper, we aim to develop a transition system
abstraction for a neural network dynamical system model defined as

###### Definition 1

A transition system abstraction is a tuple
$\mathcal{T}\triangleq\left\langle\mathcal{X},P,\mathcal{E}\ \right\rangle$
consisting of the following components:

  * •

_Working zone_ $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is the working zone that
users are paying attention to.

  * •

_Partitions_ : $P=\\{P_{1},\cdots,P_{n}\\}$ is a finite set of partitions,
which are subspaces of $\mathcal{X}$ where: 1) $P_{i}\subset\mathcal{X}$; 2)
$\mathcal{X}=\bigcup{P_{i}}$; 3) $P_{i}\bigcap P_{j}=\emptyset$.

  * •

_Transitions_ : $\mathcal{E}\subset\mathbb{B}^{n\times n}$ is the transition
matrix where $e_{ij}$ is the $i$th line, $j$th row boolean variable of
transition matrix that describes whether there exists a transition from $i$th
partition to $j$th partition. For example, a transition from $i$th partition
to $j$th partition is taking place when $x(k)\in P_{i}$ and $x(k+1)\in P_{j}$,
$e_{ij}=1\in\mathcal{E}$, if there is no transition, $e_{ij}=0$.

In this section, we will be presenting the preliminaries and problem
formulation in this paper, which includes the variation of the Shannon
Entropy, feed-forward neural network, and its reachability analysis.

###  II-A Variation of the Shannon Entropy

Learning models are usually established under limited prior knowledge, which
means the system information is provided mostly by samples. The Shannon
Entropy measures the uncertainty or disorder in a given data distribution
while only relying on samples. Specifically, according to [18] the Shannon
Entropy of the dynamics with $l\in\mathbb{R}$ partitions (denoted as $H(l)$)
is

$\displaystyle H(l)=-\sum_{i=1}^{l}p_{i}\log_{2}{p_{i}},$ (2)

where $p_{i}$ denotes the probability of the $i$th partition of the dynamics.
If we partition the dynamics into $k\geq l$ local partitions, we can compute
the variation of the Shannon Entropy based on (2) in

$\displaystyle\Delta H=H(k)-H(l)$ (3)

where $\Delta H$ reviles the variation of the Shannon Entropy, which helps in
identifying the difference between different partitioning sets. We denote the
probability of one partition in the form of

$\displaystyle p_{i}=\frac{N_{i}}{N},$ (4)

in which $p_{i}$ is the occurrence probability of partition $\mathcal{P}_{i}$
in (2), $N_{i}$ is the number of samples that are within $\mathcal{P}_{i}$,
and $N$ is the number of all samples.

###  II-B Neural Network Dynamical System Model and Reachability Analysis

Our proposed framework focuses on developing an interpretable neural
transition model on an existing trained neural network dynamical system model
in the form of (1). Among all the learning structures of neural networks,
feed-forward neural networks have received particular attention for their
simple structure and varied training approaches. In this section, We will take
a feed-forward neural network as an example to investigate the internal
propagation mechanism of a neural network dynamical system model.

For a feed-forward neural network
$\Phi:\mathbb{R}^{n_{0}}\to\mathbb{R}^{n_{L}}$ with $L$ layers, the process of
inter-layer propagation is in the form of

$\displaystyle x_{k+1}=h_{k}(W_{k}x_{k}+b_{k})$ (5)

in which $x_{k}\in\mathbb{R}^{n_{k}}$ is the output of $k$th layer,
$W_{k}\in\mathbb{R}^{n_{k+1}}\times\mathbb{R}^{n_{k}}$ is the weight matrix
that connects $k$th layer to $k+1$th layer, $b_{k}\in\mathbb{R}^{n_{k}}$,
$h_{k}:\mathbb{R}^{n_{k}}\to\mathbb{R}^{n_{k}}$ is the activation function of
$k$th layer.

Based on (5), namely the weight and bias are given, we will be able to compute
the reachable set output defined as follows.

###### Definition 2

For a neural network dynamical system model
$\Phi:\mathbb{R}^{m}\to\mathbb{R}^{n}$, when given a set input
$\mathcal{X}\subset\mathbb{R}^{m}$, we define

$\displaystyle\mathcal{Y}=\\{y\in\mathbb{R}^{n}\mid y=\Phi(x),\forall
x\in\mathcal{X}\\},$ (6)

as the reachable set output of the neural network $\Phi$ when given input set
$\mathcal{X}$.

The reachable set defined in Definition 2 facilitates verification through the
input-output set, which presents a unique angle to tackle the challenge of
black-box verification.

###### Definition 3

Given initial state set $\mathcal{X}_{(0)}$ and input set $\mathcal{U}$, the
reachable set at $k$th time step $\mathcal{X}_{(k)}$ of a neural network
dynamical system model $\Phi$ is defined as

$\displaystyle\mathcal{X}_{(k)}=\left\\{x(k;x_{0},u(\cdot))\in\mathbb{R}^{n}\mid
x_{0}\in\mathcal{X}_{(0)},~{}u(k)\in\mathcal{U}\right\\},$ (7)

and the union of $\mathcal{X}_{(k)}$ over $[0,K]$ defined by

$\mathcal{R}_{(K)}=\bigcup_{k=0}^{K}\nolimits\mathcal{X}_{(k)},$ (8)

is the reachable set over time interval $[0,K]$.

###### Remark 1

Instead of relying on the input-output point of view, reachability analysis
makes the set-valued computation possible so that we can handle the
perturbations and inaccuracies in safety-critical scenarios. However, the
reachability analysis heavily relies on iterations of the reachable set
computation which can be very computationally expensive and time-consuming.

###  II-C Problem Formulation

Neural networks are often considered black-box models, primarily due to their
limited interpretability. Consequently, this presents challenges for real-time
verification, as available methods are constrained and computational costs are
high.

###### Problem 1

Based on a neural network dynamical system $\Phi$, how can we abstract one
general $\Phi$ into a transition system $\mathcal{T}$ that makes it possible
to avoid computational-expensive real-time verification while enhancing the
interpretability of a neural network dynamical system?

In the rest of the paper, we will focus on solving Problem 1 in detail.

##  III Transition System Abstraction via ME Partitioning and Reachability
Analysis

We aim to solve Problem 1 with our proposed transition system abstraction
while computing its transitions with set-valued reachability analysis to
generate the abstracted transition graph to reveal its inner transitions.

To begin with, the neural network dynamical system will be abstracted within
the working zone. The assumption of the working zone is given as follows.

###### Assumption 1

The working zone of neural network dynamical system $\Phi$ is within the
localized state space $x\in\mathcal{X}$, given the external input bound where
$u\in[\underline{u},\overline{u}]$.

###### Remark 2

A neural network dynamical system is trained based on samples, which can be
well-trained under an adequate sample set. Therefore, we assume that the
learning model applies only to a localized working zone based on samples. This
localized approach allows for better abstraction of the dynamics within that
specific region.

###  III-A ME Partitioning

Based on the working zone of the neural dynamical system, we are able to
generate numerous traces that are included in the set of samples when given
random initial starting points with (1).

###### Definition 4

Set of samples $\mathcal{W}=\\{w_{1},w_{2},\cdots,w_{L}\\}$ of neural network
dynamical system (1) is a collection of sampled $L$ traces obtained by
measurement, where for each trace $w_{i}$, $i=1,\ldots,L$, is a finite
sequence of time steps and data
$(k_{0,i},d_{0,i}),(k_{1,i},d_{1,i}),\cdots,(k_{M_{i},i},d_{M_{i},i})$ in
which

  * •

$k_{\ell,i}\in(0,\infty)$ and $k_{\ell+1,i}=k_{\ell,i}+1$,
$\forall\ell=0,1,\ldots,M_{i},\ \forall i=1,2,\ldots,L$.

  * •

$d_{\ell,i}=[x^{\top}_{i}(k_{\ell,i}),~{}u^{\top}_{i}(k_{\ell,i})]^{\top}\in\mathbb{R}^{n_{x}+n_{u}}$,
$\forall\ell=0,1,\ldots,M_{i}$, $\forall i=1,2,\ldots,L$, where
$x_{i}(k_{\ell,i}),u_{i}(k_{\ell,i})$ denote the state and input of the system
at $\ell$th step for $i$th trace, respectively.

In this section, we will be presenting the ME partitioning technique that
divides one general working zone $\mathcal{X}$ into a set of partitions
$\\{P_{1},\ldots,P_{n}\\}$ based on the set of samples, while computing the
transitions between partitions with reachable set computation.

First, ME partitioning will partition the working zone into multiple
partitions based on the Shannon Entropy of the abstraction. Based on (8), it
takes the samples within one partition for Shannon Entropy computation. We
define the segment data set, which is subsequently defined under this process,
as follows.

###### Definition 5

Given a set of samples $\mathcal{W}$ of neural network dynamical system $\Phi$
and a set of partitions $P=\\{P_{1},\ldots,P_{N}\\}$, the segment data set is

$\displaystyle X=\\{X_{1},X_{2},\ldots,X_{N}\\}$ (9)

in which

$\displaystyle X_{q}\triangleq\\{x_{i}(k_{\ell,i})\mid x_{i}(k_{\ell,i})\in
P_{q}\\}$ (10)

in which $\forall i=1,2,\ldots,L$, $\forall\ell=0,1,\ldots,M_{i}$,
$x_{i}(k_{\ell,i})$ is any sampled state $\forall q=1,\ldots,N$ in traces
$w_{i}$. Therefore, $X_{q}$ contains all the sampled state traces evolving
within partition $P_{q}$.

###### Remark 3

In this paper, the working zone $\mathcal{X}$ is defined by hyper-rectangles
which are given as follows: For any bounded state set
$\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$, we define
$\mathcal{X}\subseteq\bar{\mathcal{X}}$, where
$\bar{\mathcal{X}}=\\{x\in\mathbb{R}^{n_{x}}\mid\underline{x}\leq
x\leq\bar{x}\\}$, in which $\underline{x}$ and $\bar{x}$ are defined as the
lower and upper bounds of $x$ in $\mathcal{X}$ as
$\underline{x}=[\mathrm{inf}_{x\in\mathcal{X}}(x_{1}),\ldots,\mathrm{inf}_{x\in\mathcal{X}}(x_{n_{x}})]^{\top}$
and
$\bar{x}=[\mathrm{sup}_{x\in\mathcal{X}}(x_{1}),\ldots,\mathrm{sup}_{x\in\mathcal{X}}(x_{n_{x}})]^{\top}$,
respectively. We denote the distance of $\mathcal{X}$ at $k$th dimension as
$d(\mathcal{X}_{k})\triangleq(\sup_{x\in\mathcal{X}}(x_{k})-\mathrm{inf}_{x\in\mathcal{X}}(x_{k}))$.

The goal of the Maximum Entropy Partitioning Algorithm is to determine the
optimal set of partitions. Specifically, we bisect one partition into two once
at a time, e.g., $P_{i}\to\\{P_{i,1},P_{i,2}\\}$ will lead to the variation of
the Shannon Entropy

$\displaystyle\Delta
H_{i}=\frac{N_{i,1}\log_{2}{\frac{N_{i,1}+N_{i,2}}{N_{i,1}}}+N_{i,2}\log_{2}{\frac{N_{i,1}+N_{i,2}}{N_{i,2}}}}{N_{i}}$
(11)

where $N_{i,1}=\left|X_{i,1}\right|$ is the sample number of $P_{i,1}$ under
(10), etc.

Based on (11), the Shannon entropy will always increase due to $\Delta
H_{i}>0$ in (11) always holds. To provide a proper partition set, we will set
a threshold $entropy>0$ for $\Delta H_{i}$ as the stop condition, namely the
bisection process will stop if $\Delta H<entropy$. The detailed ME
partitioning is given in pseudo-code in Algorithm 1 in which we set a minimum
threshold $L_{min}$ for the length of the partitions.

Input : Set of samples $\mathcal{W}$; Working Zone $\mathcal{X}$; Threshold
for Variation of The Shannon Entropy $entropy$; Threshold for minimum length
of the partitions $L_{min}$.

Output : Set of Partitions $\\{P_{1},P_{2},\ldots,P_{M}\\}$.

1 $P_{save}\leftarrow\emptyset$; $\mathcal{W}_{save}\leftarrow\emptyset$;

2 $P_{1}\leftarrow\mathcal{X}$; $r\leftarrow\dim{\mathcal{X}}$;
$Length\leftarrow\max(d(\mathcal{X}))$;

3 Obtain $X_{1}$ from $\mathcal{W}$; $q\leftarrow|X_{1}|$;

4

5while  _$Length\geq L_{min}$_ do

6 $[i,Length,k]\leftarrow$ $\max(d(P))$

7 Obtain $P_{temp1}$ and $P_{temp2}$ by bisecting $P_{i}$ on dimension $k$

8 Obtain $X_{temp1}$ and $X_{temp2}$ based on $P_{temp1}$ and $P_{temp2}$

9 if  _the variation of the Shannon Entropy $\Delta H_{i}\geq entropy$_ then

10 $P_{i}\leftarrow\\{P_{temp1},P_{temp2}\\}$

11 $X_{i}\leftarrow\\{X_{temp1},X_{temp2}\\}$

12 else

13 Add $P_{i}$ to $P_{save}$ and delete $P_{i}$

14 Add $\mathcal{W}_{i}$ to $\mathcal{W}_{save}$ and delete $\mathcal{W}_{i}$

15 end if

16

17 end while

18$P\leftarrow\\{P,P_{save}\\}$

return _$\\{P_{1},P_{2},\ldots,P_{M}\\}$._

Algorithm 1 Pseudo Code for ME Partitioning

###  III-B Transition Computation via Reachability Analysis

After ME partitioning, we will be able to compute the transition relationship
between partitions through the reachability analysis. According to (7), we
apply the reachable set computation technique in the form of

$\displaystyle\mathcal{Y}=[\Phi](\mathcal{X})$ (12)

in which $[\Phi]$ is a reachable set computation tool such as the Neural
Network Verification (NNV) toolbox in [6], the set-valued reachability
analysis in [19], Beta-crown in [5], etc.

When given the input of one partition $P_{i}$, we will be able to detect
whether the reachable set output of the neural dynamical system intersects
with any partitions (included itself), namely, given input set $\mathcal{U}$
if $[\Phi](P_{i},\mathcal{U})\cap P_{j}\neq\emptyset$, then $e_{ij}=1$, based
on which the transitions can be obtained. The process of transition
computation is given in Algorithm 2.

Input : Neural Network Dynamical System $\Phi(\cdot)$; Set of Partitions $P$.

Output : Set of Transitions $\mathcal{E}$.

1 Get the number of partitions $n$ from $P$

2 $i\leftarrow 1$; $j\leftarrow 1$;

3 while  _$i\leq n$_ do

4 Computing $[\Phi](P_{i})$ using (12)

5 while  _$j\leq N$_ do

6 if  _$[\Phi](P_{i})\cap P_{j}\neq\emptyset$_ then

7 $e_{ij}=1$

8 else

9 $e_{ij}=0$

10 end if

11 $j\leftarrow j+1$

12 end while

13 $i\leftarrow i+1$

14 end while

return _$\mathcal{E}=\begin{bmatrix}e_{11} &\cdots&e_{1n}\\\
\vdots&\ddots&\vdots\\\ e_{n1}&\cdots&e_{nn}\end{bmatrix}$._

Algorithm 2 Pseudo Code for Transition Computation

According to the transition system in Definition 1, $e_{ii}=1$, if
$[\Phi](P_{i},\mathcal{U})\cap P_{i}\neq\emptyset$, this will result in the
self-loop of transition system abstraction. From the verification point of
view, $e_{ii}=1$ indicates two scenarios:

  1. 1.

$P_{i}$ contains invariant set for neural network dynamical system $\Phi$
given the external set $\mathcal{U}$.

  2. 2.

$P_{i}$ does not contain an invariant set, $e_{ii}=1$ only because sampling
time being too short, i.e., $P_{i}$ does not contain an invariant set when
$\exists x(k_{l,i})\in P_{i},\ x(k_{l+1,i})\in P_{i}$.

Scenario 2 may lead to potential misjudgments when we apply the verification
technique such as CTL in verification [20]. Therefore, we should take
$e_{ii}=0$ under certain conditions. In practice, we reduce the self-loop
under Scenario 2 with the help of the traces that are segmented by partitions.

###### Definition 6

Given trace $w_{i}$, partition $P_{q}$, the segment of $w_{i}$ is

$\displaystyle X_{i,q}=\\{x_{i}(k_{\ell,i})\mid x_{i}(k_{\ell,i})\in
P_{q}\\},\ \forall\ell=0,1,\ldots,M_{i}$ (13)

in which $P_{q}$ and $w_{i}$ are the given partition and trace.

Under Definition 6, we can reduce the self-loop on Scenario 2 with the
following assumption.

###### Assumption 2

Given threshold $n\in\mathbb{N}$, if there is no invariant set within $P_{q}$
when given external input set $\mathcal{U}$, then $|X_{i,q}|\leq n,\forall
i=1,\ldots,L$.

Assumption 2 can be very helpful in simplifying the self-loop under Scenario 2
when given adequate neural network dynamical system information.

Neural network dynamical systems $\Phi$ are often trained as the approximation
of the system’s ideal dynamical system description $f$ in the form of

$\displaystyle x(k+1)=f(x(k),u(k)).$ (14)

To make sure we have an accurate abstraction for the dynamical system, we can
provide guaranteed transition computation, which means our transition
abstraction can still be applied for $f$ with the following proposition.

###### Proposition 1

Given an output reachable set estimation in
$\mathcal{X}^{\prime}=[\Phi](P_{i})$ and $\mathcal{X}^{\prime}\subset P_{j}$,
$\dim(\mathcal{X}^{\prime})=m$, the estimation for transitions of the neural
transition system will be the same as the dynamics $f$, there is a transition
in neural transition system as is in the ideal dynamical system description
$e_{ij=1}$, if: (1)
$\underline{x}^{\prime}_{k}-\epsilon\geq\underline{p}_{j,k}$, (2)
$\overline{x}^{\prime}_{k}+\epsilon\leq\overline{p}_{j,k}$, $\forall k\leq m$,
where $\epsilon\geq 0$ is the guaranteed distance between dynamics $f$ and its
neural network approximation $\Phi$ in [8], $\underline{x}^{\prime}_{k}$ is
the lower bound of $\mathcal{X}^{\prime}$ at $k$th dimension, etc.

Proof: To begin with, the guaranteed distance (Definition in [8]) is

$\displaystyle\left\|\Phi(x)-f(x)\right\|\leq\epsilon,$

as for the neural transition systems’ successive $x^{\prime}\in P_{j}$, it
yields $\left\|x^{\prime}-f(x)\right\|\leq\epsilon$. It leads to

$\displaystyle\begin{cases}\overline{x}^{\prime}_{k}\geq\overline{f}_{k}(x)-\epsilon\\\
\underline{x}^{\prime}_{k}\leq\underline{f}_{k}(x)+\epsilon\end{cases}\forall
k\leq m,$ (15)

in which $\overline{f}_{k}(x)$ is the $k$th dimension upper bound of the
output set for $f$, etc. If the successive set $\mathcal{X}^{\prime}\subset
P_{j}$, then $\Phi(x)\in P_{j},\ \forall x\in P_{i}$, which leads to

$\displaystyle\begin{cases}\underline{x}^{\prime}_{k}-\epsilon\geq\underline{p}_{j,k}\\\
\overline{x}^{\prime}_{k}+\epsilon\leq\overline{p}_{j,k}\end{cases}\forall
k\leq m.$ (16)

The proof is complete by combing (15) with (16), which shows $f(x)\in P_{j},\
\forall x\in P_{i}$, indicating there is a transition in $\mathcal{T}$ and the
ideal dynamical system description $f$. $\hfill\hfill\blacksquare$

Transition system abstraction for a neural network dynamical system can be
summarized as follows.

  * •

The localized working zone, i.e., partitions can be obtained based on an ME
partitioning process, which is completely data-driven and can be easily tuned
by adjusting the threshold.

  * •

The transitions can be off-line computed by reachability analysis. To provide
accurate transition estimations between the $\Phi$ and $f$, the guaranteed
distance is introduced in the proposed Proposition 1 for our abstraction
system.

After ME partitioning and transition computation, we will be able to abstract
the neural network dynamical system $\Phi$ into a transition system
abstraction $\mathcal{T}=\left\langle\mathcal{X},P,\mathcal{E}\right\rangle$.

##  IV Application to Learning and Verification on Human Handwriting Dynamics

Learning-based methods have been promoted as an effective way to model motions
[21, 22]. Consider the complex dynamical system such as human behaviors
borrowed from the LASA data set [23] which contains handwriting motion
demonstrations of human users of $30$ shapes. The challenges of modeling this
kind of complex dynamical system can be summarized as follows.

  * •

The lack of understanding of human complex behavior. Limited samples may lead
to a deep neural network dynamical system model not getting well-trained.

  * •

Human demonstrations contain abrupt changes, which means the behaviors of the
trained neural network dynamical system model will have differences in local
areas of the working zone.

  * •

We have limited methods to verify the neural network dynamical system model
for its lack of interpretability.

Due to the specifics above, we will demonstrate our proposed abstraction-based
framework based on the neural network dynamical system approximation of the
Human demonstrations. This process can be summarized as follows111The code for
the human handwriting application and transition system abstraction tool is
available at: https://github.com/aicpslab/Abstraction-Toolbox-for-Neural-
Transition-Systems.

  * •

We train $20$-ReLU activated neurons ELMs as the neural dynamical system model
for learning Angle and P shape handwriting demonstrations. By setting the
maximum and minimum values of the sample data, we obtain the working zone of
neural network dynamical systems.

  * •

A threshold $entropy=4\times e^{-2}$ is set as the variation of maximum
entropy in Algorithm 1, which results in $28$ partitions in transition system
abstraction for Angle shape Fig. 1 (a) and $27$ partitions in learning model
for P shape Fig. 1 (b). We partition the working zone based on the randomly
generated $20$ trajectories.

  * •

Based on the partitions of the neural transition systems, we employ the neural
verification toolbox for transition computations in Algorithm 2, and reduce
the self-loop with $n=100$ in Assumption 2. Then we will be able to provide
the transition graphs in 2 (a) and (b). From Fig. 1 to Fig. 2, the transition
relationships between local working zones are able to be interpreted through
our transition system abstraction.

  * •

We verify our transition system abstraction with the Computation Tree Logic
(CTL) formulae [24], in which $\diamond$ or $\Box$ denote the for $some$ or
$all$ traces, $\bigcirc$ denotes for next step. The results are given in Tab
I, in which $\phi_{1}$ is true means there exists a trace from $P_{28}$ to
$P_{12}$, $\phi_{2}$ is true for $\mathcal{T}_{Angle}$ means for all traces
that start from $P_{28}$, it will go to $P_{25}$ as the next partitions,
$\phi_{3}$ means for $\mathcal{T}_{Pshape}$ there exists a trace in $P_{4}$
then will be in $P_{16}$ next.

(a) $|P_{Angle}|=28$. (b) $|P_{Pshape}|=27$.

Figure 1: There are $28$ partitions from (a) while $27$ partitions from (b),
information-rich areas are allocated more partitions.

(a) Transition system abstraction $\mathcal{T}_{Angle}$. (b) Transition system
abstraction $\mathcal{T}_{Pshape}$.

Figure 2: The transition relationships are denoted by red arrows. Based on (a) and (b) we can clearly see the model’s working flow, which means the learning model’s interpretability is enhanced, such as results in Table I where complex CTL formulae can be verified. TABLE I: Verification results of $CTL$ formula: $\mathcal{T}_{Angle}$ with $P_{27}$ as the initial partition, and $\mathcal{T}_{Pshape}$ with $P_{26}$ as the initial partition.  $CTL$ formula | $\mathcal{T}_{Angle}$ | $\mathcal{T}_{Pshape}$
---|---|---
$\phi_{1}=\diamond\exists P_{12}$ | $true$ | $-$
$\phi_{2}=\Box\bigcirc P_{25}$ | $true$ | $-$
$\phi_{3}=\diamond\exists(P_{4}UP_{16})$ | $-$ | $true$

##  V Conclusions

In this paper, a transition system abstraction-based framework for the neural
network dynamical system model is given to promote interpretability. First,
the working zone of the transition system will be obtained based on the
samples. Then, we employ the ME partitioning to bisect the working zone of the
neural network dynamical system model into multiple partitions, while the
transitions of the learning model will be computed under reachable set
computations. Transition system-based abstractions and verification for neural
dynamical systems of Angle shape and P-shape are presented to illustrate the
effectiveness, which indicates the interpretability is enhanced while
promoting a new way of neural network dynamical system model verification.

## References

  * [1] S. Wang, Y. Fan, S. Jin, P. Takyi-Aninakwa, and C. Fernandez, “Improved anti-noise adaptive long short-term memory neural network modeling for the robust remaining useful life prediction of lithium-ion batteries,” Reliability Engineering & System Safety, vol. 230, p. 108920, 2023.
  * [2] A. Liu, T. Chen, H. Zhu, M. Fu, and J. Xu, “Fuzzy variable impedance-based adaptive neural network control in physical human–robot interaction,” Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 237, no. 2, pp. 220–230, 2023.
  * [3] M. Ghazvini, H. Maddah, R. Peymanfar, M. H. Ahmadi, and R. Kumar, “Experimental evaluation and artificial neural network modeling of thermal conductivity of water based nanofluid containing magnetic copper nanoparticles,” Physica A: Statistical Mechanics and its Applications, vol. 551, p. 124127, 2020.
  * [4] H. Zhang, T.-W. Weng, P.-Y. Chen, C.-J. Hsieh, and L. Daniel, “Efficient neural network robustness certification with general activation functions,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [5] S. Wang, H. Zhang, K. Xu, X. Lin, S. Jana, C.-J. Hsieh, and J. Z. Kolter, “Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification,” Advances in Neural Information Processing Systems, vol. 34, 2021.
  * [6] H.-D. Tran, X. Yang, D. Manzanas Lopez, P. Musau, L. V. Nguyen, W. Xiang, S. Bak, and T. T. Johnson, “NNV: the neural network verification tool for deep neural networks and learning-enabled cyber-physical systems,” in International Conference on Computer Aided Verification, pp. 3–17, Springer, 2020.
  * [7] M. Fazlyab, M. Morari, and G. J. Pappas, “Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming,” IEEE Transactions on Automatic Control, vol. 67, no. 1, pp. 1–15, 2020.
  * [8] Y. Yang, T. Wang, J. P. Woolard, and W. Xiang, “Guaranteed approximation error estimation of neural networks and model modification,” Neural Networks, vol. 151, pp. 61–69, 2022.
  * [9] M. Poli, S. Massaroli, L. Scimeca, S. Chun, S. J. Oh, A. Yamashita, H. Asama, J. Park, and A. Garg, “Neural hybrid automata: Learning dynamics with multiple modes and stochastic transitions,” Advances in Neural Information Processing Systems, vol. 34, pp. 9977–9989, 2021.
  * [10] C. Sprague, Efficient and Trustworthy Artificial Intelligence for Critical Robotic Systems.  PhD thesis, Kungliga Tekniska högskolan, 2022.
  * [11] O. A. Beg, H. Abbas, T. T. Johnson, and A. Davoudi, “Model validation of PWM DC-DC converters,” IEEE Transactions on Industrial Electronics, vol. 64, no. 9, pp. 7049–7059, 2017.
  * [12] A. Chowdhery, S. Narang, J. Devlin, M. Bosma, G. Mishra, A. Roberts, P. Barham, H. W. Chung, C. Sutton, S. Gehrmann, et al., “Palm: Scaling language modeling with pathways,” arXiv preprint arXiv:2204.02311, 2022.
  * [13] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, “Adaptive mixtures of local experts,” Neural Computation, vol. 3, no. 1, pp. 79–87, 1991.
  * [14] S. Horváth, D. Kovalev, K. Mishchenko, P. Richtárik, and S. Stich, “Stochastic distributed learning with gradient quantization and double-variance reduction,” Optimization Methods and Software, vol. 38, no. 1, pp. 91–106, 2023.
  * [15] X. Li, H. Xiong, X. Li, X. Wu, X. Zhang, J. Liu, J. Bian, and D. Dou, “Interpretable deep learning: Interpretation, interpretability, trustworthiness, and beyond,” Knowledge and Information Systems, vol. 64, no. 12, pp. 3197–3234, 2022.
  * [16] J. Liu, N. Ozay, U. Topcu, and R. M. Murray, “Synthesis of reactive switching protocols from temporal logic specifications,” IEEE Transactions on Automatic Control, vol. 58, no. 7, pp. 1771–1785, 2013.
  * [17] Y. Yang, Z. Mo, and W. Xiang, “A data-driven hybrid automaton framework to modeling complex dynamical systems,” in 2023 IEEE International Conference on Industrial Technology (ICIT), pp. 1–6, IEEE, 2023.
  * [18] V. Rajagopalan and A. Ray, “Symbolic time series analysis via wavelet-based partitioning,” Signal processing, vol. 86, no. 11, pp. 3309–3320, 2006.
  * [19] W. Xiang, H.-D. Tran, and T. T. Johnson, “Output reachable set estimation and verification for multilayer neural networks,” IEEE transactions on neural networks and learning systems, vol. 29, no. 11, pp. 5777–5783, 2018.
  * [20] M. Lahijanian, S. B. Andersson, and C. Belta, “Formal verification and synthesis for discrete-time stochastic systems,” IEEE Transactions on Automatic Control, vol. 60, no. 8, pp. 2031–2045, 2015.
  * [21] R. F. Reinhart, Z. Shareef, and J. J. Steil, “Hybrid analytical and data-driven modeling for feed-forward robot control,” Sensors, vol. 17, no. 2, p. 311, 2017.
  * [22] R. F. Reinhart and J. J. Steil, “Neural learning and dynamical selection of redundant solutions for inverse kinematic control,” in 2011 11th IEEE-RAS International Conference on Humanoid Robots, pp. 564–569, IEEE, 2011.
  * [23] M. Khansari, E. Klingbeil, and O. Khatib, “Adaptive human-inspired compliant contact primitives to perform surface–surface contact under uncertainty,” The International Journal of Robotics Research, vol. 35, no. 13, pp. 1651–1675, 2016.
  * [24] H. Pan, Y. Li, Y. Cao, and Z. Ma, “Model checking computation tree logic over finite lattices,” Theoretical computer science, vol. 612, pp. 45–62, 2016.