# Decomposing the Dynamics of the Lorenz 1963 model using Unstable Periodic
Orbits: Averages, Transitions, and Quasi-Invariant Sets

Chiara Cecilia Maiocchi<EMAIL_ADDRESS>Valerio Lucarini
<EMAIL_ADDRESS>Centre for the Mathematics of Planet Earth,
University of Reading  Department of Mathematics and Statistics, University of
Reading

###### Abstract

Unstable periodic orbits (UPOs) are a valuable tool for studying chaotic
dynamical systems. They allow one to extract information from a system and to
distill its dynamical structure. We consider here the Lorenz 1963 model with
the classic parameters’ value and decompose its dynamics in terms of UPOs. We
investigate how a chaotic orbit can be approximated in terms of UPOs. At each
instant, we rank the UPOs according to their proximity to the position of the
orbit in the phase space. We study this process from two different
perspectives. First, we find that, somewhat unexpectedly, longer period UPOs
overwhelmingly provide the best local approximation to the trajectory, even if
our UPO-detecting algorithm severely undersamples them. Second, we construct a
finite-state Markov chain by studying the scattering of the forward trajectory
between the neighbourhood of the various UPOs. Each UPO and its neighbourhood
are taken as a possible state of the system. We then study the transitions
between the different states. Through the analysis of the subdominant
eigenvectors of the corresponding stochastic matrix we provide a different
interpretation of the mixing processes occurring in the system by taking
advantage of the concept of quasi-invariant sets.

††preprint: APS/123-QED

##  I Introduction

Unstable periodic orbits (UPOs) play an important role in the analysis of
dynamical systems that exhibit chaotic behaviour and in some cases they
provide a powerful tool to describe the statistical properties of the system.
As early noticed by Poincaré, poincare_1893  unstable periodic orbits
represent a rigid skeleton hidden in the chaos of the dynamics and might be
the key to the understanding of chaotic systems cvitanovic_1991 . Since then,
UPOs have been at the centre of mathematical investigation. Unstable periodic
orbits are skeletal for the dynamics in the sense that they can be used to
reconstruct the statistical properties of a chaotic dynamical system
grebogi_1988 . In fact, since UPOs are dense in the attractor eckmann_1985 ,
they can approximate any trajectory in the system with an arbitrary accuracy
bowen_1975 . UPOs can be considered as islets of order in a landscape of
chaos. The trajectory moves from the neighbourhood of one UPO to another,
undergoing a mechanism of repulsion because of their instability. UPOs provide
a modal decomposition of the dynamics, enabling to identify the dynamical
properties of the trajectory with those of neighbouring UPOs. Within this
context it is possible to develop a theory that allows to write dynamical
averages as weighted sum over the full set of UPOs. Gutzweller gutzwiller_2013
first demonstrated that UPOs are the essential building blocks of chaotic
dynamics. Ruelle later derived the dynamical $\zeta$ function ruelle_2004  ,
that allows one to write the natural measure of the system as a weighted sum
over the infinite set of UPOs. It is well known that a complete understanding
of the full set of short UPOs would suffice in obtaining an accurate
approximation of ergodic averages eckhardt_1994 ; cvitanovic_1995 ;
artuso_1990 . These results are proven to be valid for dynamical systems
exhibiting strong chaoticity ruelle_1999 ; katok_1997 , such as hyperbolic and
Axiom A systems smale_1967 ; bowen_1972 . However, in complex models of fluid
flows, it is often difficult, if not impossible, to verify the hypothesis
required for the validity of periodic orbit expansion. When turbulent
conditions are considered, such systems live, after transients have died out,
in nonequilibrium steady state (NESS) gallavotti_2014 ; lucarini_2014 . This
state is in general characterised by generation of entropy, contraction of
phase space and finite-time predictability. The ’Chaotic hypothesis’ of
Gallavotti gallavotti_1998 ; gallavotti_1995  offers a possible solution to
the first problem, allowing to consider ’a turbulent fluid as a transitive
Axiom A system for the purpose of computing macroscopic properties of the
system.’

###  I.1 Unstable Periodic Orbits: Applications

A first application of periodic orbit expansion was performed by Auerbach et
al. in auerbach_1987 , by proving that UPOs are experimentally accessible and
capable of unfolding the structure of chaotic trajectories. In fact, by
extracting the complete set of UPOs of order $n$ and calculating their
stabilities, they approximated the fractal dimension and topological entropy
of the strange attractor of the paradigmadic Hénon map with very high
accuracy. Cvitanović cvitanovic_1988  argued that UPOs are the optimal
practical tool for measuring the invariant properties of a dynamical system
and provided a solid ground for the applications of periodic orbit theory
(POT). Artuso et al. tested this procedure through a series of applications
artuso_1990 ; artuso_1990_2  and demonstrated that cycle expansion of the
dynamical $\zeta$ function is instrumental for the analysis of deterministic
chaos, even in more generic settings than the ones required by cvitanovic_1988
, i.e. when the system is not uniformly hyperbolic nor the results depend on
the assumption of the existence of invariant measures or structural stability
of the dynamics. Eckhardt and Ott eckhardt_1994  presented one of the first
numerical applications of the periodic orbit formalism for studying the
statistical and the dynamical properties of the Lorenz 1963 (L63) system
lorenz_1963 . A subsequent analysis of the linear and nonlinear response of
the L63 to perturbations found that specific UPOs were responsible for
resonance mechanisms leading to an amplified response lucarini2009 .

Later on, periodic orbit theory found fruitful applications also within the
context of higher dimension forced and dissipative non-equilibrium systems,
specifically in the case of (geophysical) fluid dynamics. UPOs can be
considered as a mean to simplify and interpret qualitative behaviour of a
complex system lucarini_2020 , allowing to extract information and distill its
dynamical structure. This observation, together with the study of the
stability and thus predictability properties of the tangent space, allows to
associate relevant dynamical features of the flow to specific UPOs or classes
of UPOs. UPOs, true nonlinear modes of flow, can be interpreted as a
generalisation of the normal modes observed in a network of coupled linear
oscillators, that allow for a study of the system in its complexity, without
the necessity of considering a heavily truncated model.

Even though a periodic-orbit-theory-type expansion for turbulence is still a
far reaching goal, many steps have been made in this direction. The path
towards the description of a turbulent flow in terms of UPOs is paved with
promising results cvitanovic_2013 . Kawahara and Kida kawahara_2001 , who
found a UPO embedded in the attractor of a numerical simulation of plane
Couette flow, showed that one UPO only manages to capture in a surprisingly
accurate way the turbulence statistics. Chandler and Kerswell chandler_2013
identified $50$ UPOs of a turbulent fluid and used them to reproduce the
energy and dissipation probability density functions of the system as
dynamical averages over the orbit. These encouraging results suggested that
periodic orbit theory could represent a valid investigation tool also in the
realm of climate systems. Lucarini and Gristun lucarini_2020  considered POT
as a novel mathematical framework for investigating atmospheric, ocean and
climate dynamics. Gritsun gritsun_2013  numerically applied the theory of UPOs
expansion ruelle_1999 ; cvitanovic_2005  to reconstruct the probability
distribution function of a simple atmospheric model based on the barotropic
vorticity equation of the sphere. Gritsun and Lucarini gritsun_2017  used the
periodic orbit machinery for interpreting non trivial resonant responses to
forcing by associating it to a set of extremely unstable UPOs.

###  I.2 This paper

Our research aims to give a contribution towards the computational problems in
the application of periodic orbits theory. We believe that POT could indeed
represent a very powerful tool allowing to gain insight and simplify complex
systems such as the climate, but its implementation is often hindered by the
numerical difficulties associated with UPOs search gritsun_2008 . We consider
the L63 model as a test case. The use of UPOs for averaging in L63 has been
hugely debated in the literature (See discussion in section II.1.2 ). We will
not delve into this matter, but we rather want to shed light on the shadowing
mechanism and investigate its underpinning dynamics. In order to avoid
confusion, since the term ”shadowing” can be used with the meaning given in
bowen_1975 , we specify that by ”shadowing” we refer to the process of
approximation of the chaotic trajectory in terms of UPOs. In particular, at
each point in time we rank UPOs based on their distance with respect to the
given chaotic trajectory. We name this process ranked shadowing. Within this
context, our approach is twofold. On the one side we aim at numerically
understand how chaotic trajectories are approximated in terms of UPOs. On the
other side, we study transitions between neighbourhood of UPOs.

In the first part of our work we analyse the shadowing process by taking a
local point of view and studying the shadowing statistics of a chaotic
trajectory. Namely, we numerically implement the shadowing process, focusing
on the point-wise dynamics along the trajectory, by locating the closest UPOs
to each point, ordered accordingly to their distance from the chaotic
trajectory. Extending periodic search to higher dimensional systems still
represents an open challenge chandler_2013  and standard periodic averaging
with an incomplete set of UPOs does not provide accurate results. A possible
strategy to overcome the aforementioned difficulties has been proposed by
Lasagna lasagna_2018 ; lasagna_2020 , who found numerical evidence that long
period orbits could be used as proxies of chaotic trajectories. Following the
idea that long UPOs might be able to capture statistical properties of chaotic
trajectories, periodic search could be limited to few long UPOs, notably
reducing the computational complexity of the problem. In line with this work,
we find that longer UPOs have the lion’s share in reproducing the invariant
measure of the system.

In the second part of this paper, inspired by the work of Froyland and Padberg
on quasi-invariant sets froyland_2009 ; froyland_2005 , we assume a global
point of view on the dynamics by extracting the dynamical behaviour of the
system via a Markov model froyland_2001 . In doing this we abandon the
classical Ulam’s partition ulam_2004  of the phase space in virtue of a
possibly novel one: each UPO is seen as a unitary building block of the
system, a spatially extended state, and the shadowing process can be
interpreted as transitions between different states. We will show that this
new perspective on the dynamics allows for a different interpretation of
quasi-invariant sets froyland_2005  in terms of UPOs and provides a better
understanding of mass transport. Namely, by studying the spectral properties
of the discretised transfer operator, we obtain a partition of the phase space
in different bundles of UPOs, each one identifying a quasi-invariant region.
We prove that UPOs represent a valid tool to investigate diffusion properties
of the system, in fact, being exact solutions, they retain a memory of the
geometrical structure of the attractor.

The structure of the rest of the document is as follows. We describe the
ranked shadowing process in section II. We introduce the mathematical
framework, describing the characteristics of the selected database of UPOs. We
then discuss the statistics emerging from the shadowing of a chaotic
trajectory. We prove the robustness of the results in dependence of the
shadowing criteria. It emerges that longer period orbits play a major role in
reproducing the invariant measure of the system. In section III we study the
shadowing dynamics from a probabilistic point of view. We introduce the
discretised transfer operator and describe transitions in terms of UPOs. We
identify quasi-invariant sets through the study of the spectrum of the
transition matrix and recollect the mass diffusion dynamics. In section IV we
outline our conclusions and perspectives for future works. In the Appendix
appendix  we provide a more extensive description of the algorithms considered
for our analysis.

##  II Shadowing mechanism

###  II.1 Experimental Setting

In this section we numerically investigate the shadowing mechanism, with the
aim of understanding whether it is possible to identify a set of key UPOs that
capture the qualitative behaviour of the system.

####  II.1.1 Mathematical Framework

We consider a continuous-time autonomous dynamical system $\dot{x}=f(x)$ on a
compact manifold $M\subset\mathbb{R}^{n}$, where $x(t)=S^{t}x_{0}$, with
$x_{0}=x(0)$ initial condition and $S^{t}$ evolution operator defined for
$t\in\mathbb{R}_{t>0}$. We define $\Omega\subset M$ as the compact attracting
invariant set of the dynamical system. A differentiable function (observable)
$\varphi:M\rightarrow\mathbb{R}$ is also given, and it specifies the
probability measure $\rho$, invariant and ergodic with respect to $S^{t}$, in
the following manner:

$<\varphi>=\int\rho(dx)\varphi(x)=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\varphi(S^{t}x_{0})dt$
(1)

for almost all initial conditions $x_{0}$ belonging to the basin of attraction
of $\Omega$.

A periodic orbit of period $T$ is defined as

$S^{T}(x)=x.$ (2)

This representation is not unique. In fact, if equation 2 is satisfied,
$S^{nT}(x)=x$ is verified as well $\forall n\in\mathbb{N}$ . By the semigroup
property of the evolution operator, we also have that $S^{T}(y)=y$ if
$y=S^{s}(x)$ for any choice of $s$. From now onward we will considered a
periodic orbit to be identified by its prime period $T$ and initial condition
$x_{0}$.

####  II.1.2 The Model

We choose as testing ground the L63 model, which represents one of the most
paradigmatic chaotic systems, since for some parameters value it exhibits a
strange attractor tucker_1999  densely populated by an infinite number of UPOs
galias_1998 .

The set of equations of the model are given by:

$\displaystyle\dot{x}_{1}=-\sigma(x_{1}+x_{2})$ $\displaystyle\dot{x}_{2}=\rho
x_{1}-x_{2}-x_{3}x_{1}$ $\displaystyle\dot{x}_{3}=-\beta x_{3}+x_{1}x_{2}$

where the three parameters $\sigma,\rho,\beta$ are positive numbers
respectively proportional to the Prandtl number, Rayleigh number and geometry
of the considered region. In this work we consider the classical Saltzmann
parameter $\sigma=10$, $\rho=28$ and $\beta=8/3$. For such values, the
dynamics of the system is characterised by a chaotic behaviour on a singularly
hyperbolic attractor that supports an SRB measure tucker_2002 . With chaotic
we indicate the property of sensitive dependence on initial conditions of the
model. In particular, the first Lyapunov exponent, that gives information on
the asymptotic rate of divergence of the trajectories, is positive
pikovsky_2016 .

Many studies on UPOs of the Lorenz system have been carried out. It is known
that a good understanding of the UPOs of the model plays a fundamental role in
the characterisation of the system. In fact, the invariant measure of the
system can be reconstructed through the use of trace formulas cvitanovic_2005
by considering the following expression for the average of any measurable
observable $\varphi$:

$\displaystyle<\varphi>=\lim_{t\to\infty}\frac{\sum_{U^{p},p\leq
t}w^{U^{p}}\bar{\varphi}^{U^{p}}}{\sum_{U^{p},p\leq t}w^{U^{p}}}$ (3)

where $U^{p}$ is a UPO of prime period $p$, $w^{U^{p}}$ is its weight and
$\bar{\varphi}^{U^{p}}$ is the average of the observable along the orbit. For
uniformly hyperbolic dynamical systems this result is exact and the weight can
be obtained, to a first approximation, by
$w^{U^{p}}\propto\exp(-ph_{ks}^{U^{p}})$ grebogi_1988  , with $h_{ks}$ being
the Kolmogorov-Sinai entropy of the system. This quantity provides information
on the rate of creation of information due to the chaoticity of the system.
From the knowledge of the Lyapunov exponents of the system $\Lambda_{i}$
pikovsky_2016 , which describe how linearised neighborhoods of the trajectory
evolve in time, we can recollect an explicit expression for $h_{ks}$ via Pesin
theorem ott_2002

$\\\ h_{ks}\leq\sum_{\Lambda_{i}>0}\Lambda_{i},$ (4)

where the left and right hand sides are equal if the invariant measure is of
the SRB type (as in the case of L63 for the chosen parameters’ value).
Eckhardt and Ott eckhardt_1994  presented one of the first numerical
applications of the periodic orbit formalism by considering an approximate
symbolic coding cvitanovic_1988  (UPOs with period up to $9$) to calculate
Hausdorff dimensions and Lyapunov exponents. Franceschini, Giberti and Zheng
franceschini_1993  calculated a number of UPOs of the Lorenz attractor at both
standard and non standard parameter values and used them to approximate the
topological entropy and Hausdorff dimension. Zoldi zoldi_1998  investigated to
what extent trace formulas can can predict the structure of the histogram of
chaotic time series data extracted from the run of the L63 model with
different parameter values. The use of a correct weighting in the trace
formula has been vastly investigated, see for a debate saiki_2010 ; saiki_2009
; zaks_2010 .

In this paper we do not want to improve the way UPOs are used to compute
sample averages, nor enlarging the already extensive database UPOs of the
system, instead we wish to investigate how well an imperfect subset of UPOs
can be used to reconstruct a long chaotic trajectory and interpret the
transitions from one neighborhood of the orbit to another.

####  II.1.3 The Database

In this section we describe some general characteristics of the collected set
of UPOs. The extrapolation of UPOs represents a complex numerical challenge.
One of the main obstacles in applying periodic orbits theory resides in fact
in the intrinsic difficulty of numerically extracting the orbits from the
model. Via Newton’s method (see Appendix appendix  for more details) we obtain
a set of $\bar{M}=34191$ orbits, whose statistics of prime periods is shown in
Fig. 1. The periods span from $T_{min}=1.5587$ to $T_{max}=17.1398$. L63,
being a chaotic system for the current choice of parameters’ value, is
characterised by the presence of a positive Lyapunov exponent in its SRB
measure, e.g. $\Lambda_{1}>0$ pikovsky_2016 . The distribution of values of
$\Lambda_{1}>0$ among the UPOs included in the database is presented in Figure
2, calculated as per benettin_1980 . As we see, the distribution of values is
peaked with a small spread of just few percents, and, in particular, no UPO
has a vanishing or negative value of $\Lambda_{1}$ (which would go against the
chaotic nature of the flow). Note that, as well known, the local instability
of the L63 model varies wildly within its attractor, where regions with very
high instability alternate with regions where one observes return-of-skill for
finite-time forecast Palmer1993 . Hence in this case, as opposed to what
observed in lucarini_2020 , the heterogeneity of the attractor in terms of
instability cannot be explained using UPOs, possibly because we are
considering here a very low-dimensional flow.

Our UPO sample does not present the characteristic exponential growth with the
period bowen_1970  and is biased towards shorter UPOs, which is known are
easier to be found gritsun_2013 ; lucarini_2020 . This issue could be resolved
by improving the search algorithm and allocating more computational resources.
Many numerical algorithms have been proposed so far. Saiki saiki_2007
reviewed the Newton-Raphson-Mees method, proposing a value for the damping
coefficient related to the stability exponent of the orbit, while Barrio et
al. barrio_2015  carried out an extensive high-precision numerical simulation
in order to gather a benchmark database of UPOs for L63. Galias and Tucker
galias_2009 , computed all the UPOs of symbolic dynamics period $T\leq 14$.
However, we want to remark that our analysis is not intended as an extensive
calculation of UPOs in L63, but aims at developing tools that would enable us
to take advantage of the periodic orbits machinery even when the symbolic
dynamic is not known or when higher dimensional systems are considered. In
virtue of this motivation we thus decided to settle with the aforementioned
database.

Figure 1: Number of detected UPOs vs their prime period Figure 2: Distribution
of the maximal Lyapunov exponent

####  II.1.4 Shadowing the Chaotic Trajectory

In this section we present results on how the UPOs shadow a chaotic
trajectory. Since the system is ergodic and we consider a long trajectory
compared to the timescale of the system, we assume that the statistics
presented here do not depend on the choice of the chaotic trajectory.

Let us denote the set of UPOs of the database as
$\mathcal{U}=\\{U_{k}\\}_{k=1}^{\bar{M}}$ where the UPO $U_{k}$ is intended as
a set of points in the system phase space
$U_{k}=\\{u_{k}(s)\\}_{s=1}^{dt*T_{k}}$, with $T_{k}$ being its period and
$dt$ the time step. We consider a numerical chaotic trajectory
$\mathcal{X}_{chaotic}$ consisting of the set of points
$\mathcal{X}_{chaotic}=\\{x_{t}\\}_{t=1}^{T_{max}}$ where
$T_{max}=N_{chaotic}*dt$ , being $N_{chaotic}$ the number of points of the
chaotic trajectory. We say that the UPO $U_{\bar{k}}$ shadows the chaotic
trajectory $\mathcal{X}_{chaotic}$ at time t if

$\displaystyle\min_{s}|u_{\bar{k}}(s)-x(t)|=\min_{k}(\min_{s}|u_{k}(s)-x(t)|)$
(5)

This calculation was carried out using all available periodic orbits, with an
output time-step $dt=0.01$ (See Appendix appendix ). for more details on the
algorithm specification).

It is clear from the definition that it is possible that not all the UPOs of
the database would be considered in the shadowing process.

Figure 3: Number of shadowing UPOs as a function of the length of the shadowed
chaotic trajectory

We can see from Fig. 3 that a chaotic orbit of length $T_{max}=10^{5}$ almost
exhausts the database without saturating the UPOs (98.6% of the UPOs of the
database are shadowing UPOs). We believe that this is the ideal length to
produce reliable statistics. From now onward all the statistics produced are
derived from data obtained by shadowing a chaotic trajectory of length
$T_{max}=10^{5}$ with the full database.

###  II.2 Statistics of the Shadowing

####  II.2.1 Shadowing time and Occupancy ratio

Our main goal is to shed light on the shadowing mechanism. It is crucial to
develop some tools that are capable of measuring to what extent each UPO
contributes in the approximation of the chaotic trajectory. With this purpose
in mind, we define the shadowing time and occupancy ratio of a UPO. The
shadowing time of a UPO is defined as the amount of time that the UPO spends
shadowing the chaotic trajectory. More precisely, if the UPO $\mathcal{U}_{k}$
is selected as shadowing orbit $t_{k}$ times, its shadowing time will be
$r_{k}=t_{k}*dt$ . This quantity is a good indicator for the absolute
shadowing time, but it does not take into account the length of the UPO. It is
sensible that a longer period UPO corresponds to a longer trajectory and
therefore, since it explores more extensively the phase space, it might
shadows the trajectory for longer time. To measure the absolute shadowing time
we introduce the occupancy ratio for the UPO $\mathcal{U}_{k}$, defined as
$o_{k}=\frac{t_{k}}{T_{k}/dt}$ with $T_{k}$ being the period of the UPO. In
this way we are able to measure the shadowing time normalised over the period
of the UPO.

####  II.2.2 Longer Period Orbits Shadow the Trajectory for a Longer Time

Figure 4:  UPOs with period larger than $5.5$ (which represent 20% of the
shadowing UPOs) shadows the trajectory $80\%$ of the total time.

(a)

(b)

Figure 5:  Shadowing time (panel LABEL:sub@error_bar_shadowing) and occupancy
ratio (panel LABEL:sub@error_bar_occupancy) increase with the period of the
UPO ($2\sigma$ error bars).

(a)

(b)

Figure 6: Distribution of the shadowing time (panel
LABEL:sub@distribution_shadowing) and occupancy ratio (panel
LABEL:sub@distribution_occupancy). Dashed lines: UPOs with period smaller than
5.5. Solid lines: UPOs with period greater than 5.5.

We find that higher period orbits have a major role in reproducing the
invariant measure of the system. In fact, not only the trajectory is shadowed
for most of the time by UPOs with period greater than $5.5$ (Fig. 4), but we
also observe a positive correlation between occupancy ratio and UPOs period
(Fig. 5b). This is even more interesting if we note that longer period UPOs
are less represented in the database we considered (Fig. 1). We also remark
that for $T<5.5$ the average occupancy ratio has a magnitude of $10^{-1}$,
while for $T>5.5$ the average is of order 1 (See Fig. 6b).

In order to assess the robustness of our results and test the sensitivity with
respect to numerical rounding, we conduct a similar analysis, allowing for the
shadowing algorithm to select $N$ shadowing UPOs at each step instead of
limiting the selection to the minimum only. We want to test whether even
allowing for a looser definition of shadowing UPOs, the positive correlation
between higher period and longer occupancy ratio remains consistent.

More precisely, for each point $x_{t}$ along the chaotic trajectory
$\mathcal{X}_{t}$ we select a set of $N$ ordered UPOs consisting of the first
$N$ UPOs that minimise the distance with $x_{t}$. In this context, we are
interested in average quantities over all tiers. Namely, we define the average
occupancy ratio at time $t$ as

$\displaystyle\bar{o}_{t}=\frac{1}{N}{\sum_{k=1}^{N}o_{k}}$ (6)

where $N$ is the number of tiers and $o_{k}$ is the UPO that shadows the
trajectory at time t. Similarly we define the average period at time $t$ in
tier $k$. Note that in each tier the shadowing UPO is chosen among the full
set of orbits, meaning that a cycle could shadow the trajectory at point $t$
in more than one tier.

Figure 7:  Solid line: Average period of the shadowing UPOs calculated over
different tiers. Each UPO is weighted proportionally to how many times it is
chosen as shadowing UPO. On the x-axis we indicate the percentage of UPO of
the full database in each tier.

Figure 7 confirms the importance of higher period UPOs. The average period of
the shadowing UPOs over all tiers is sensibly higher than the average period
of the set of shadowing UPO and it monotonically decreases as orbits in
further tiers are selected. This reflects the fact that longer period UPOs are
indeed prioritised in the shadowing process and this choice is consistent over
the full length of the chaotic trajectory.

The robustness of the analysis is confirmed when reproducing the statistics
presented in Fig. 5 with $N$ shadowing UPOs. The results are presented in Fig.
8. Allowing for more shadowing UPOs does not affect the correlation found in
the previous section when considering average quantities.

(a)

(b)

Figure 8: Shadowing time (panel LABEL:sub@shad_moretiers) and occupancy ratio
(panel LABEL:sub@occ_moretiers)averaged over first 10 (black), 30 (red) and
100(blue) tiers ($2\sigma$ error bar).

These results are robust when under sampling the set of UPOs down to 2 orders
of magnitude. They provide an insight into the shadowing dynamics and are
particularly useful in the applications of periodic orbits theory. In fact,
knowing a priori that longer UPOs play a major role in reproducing the
invariant measure of the system could give a contribution to the optimisation
of periodic search algorithms by addressing computational resources to a
limited subset of UPOs. Hence our findings seem to provide support to the
approach proposed by Lasagna in lasagna_2018 ; lasagna_2020 , nonetheless it
is not clear to us why this might hold true.

##  III TRANSITIONS

In this section we use UPOs as a tool to investigate the mixing properties of
the system. The ranked shadowing will be used to define the Markov process
that describes transitions between neighbourhood of UPOs.

###  III.1 Extracting a Markov Chain from the Dynamics

A very valuable tool to study transitions and evolution of measures in
dynamical systems is offered by the transfer operator. Let
$\textbf{F}:\Omega\in\mathbb{R}^{d}\to\mathbb{R}^{d}$ be a smooth vector
field, generating the dynamical system or flow
$\\{\Phi^{t}\\}_{t\in\mathbb{R}}$, with $\Phi_{t}:\Omega\to\Omega$. We define
froyland_2009  the Perron-Frobenius operator or transfer operator
$\mathcal{P}_{t}:L^{1}(\Omega)\to L^{1}(\Omega)$ by

$\displaystyle\mathcal{P}_{t}f(x)=f(\Phi_{-t}(x))|det(D\Phi_{-t}(x))|,$ (7)

where $D\Phi_{-t}(x)$ is the Jacobian matrix corresponding to the spatial
derivatives of $\Phi_{t}$. The transfer operator is a statistical tool that
describes how the densities globally evolve with time. From the study of its
spectral properties we can deduce significant statistical information about
the system, such as mixing properties, invariant densities and decay of
correlations baladi_2000 ; froyland_2001 . For instance, fixed points of
$\mathcal{P}_{t}$ represent invariant densities for the dynamics, that remain
unaltered by the flow.
We need to define an appropriate numerical estimate of the transfer operator
$\mathcal{P}_{t}$. In fact, in order to tackle the problem from a numerical
standpoint, we have to consider the transfer operator within a finite
dimensional setting, where the phase space is not interpreted as a continuum,
but it is appropriately discretised into a finite collection of regions, with
mass moving from one region to the other at each iteration of the transfer
operator. It is important to notice that at this stage the dynamics occurring
within each set of the partition is ignored, and we are just interested in the
macroscopic movement of mass. Different methods for defining this
approximation have been developed. For the well-known Ulam’s method ulam_2004
the approximation takes the form of a regular lattice covering the phase
space. More precisely, let $\\{B_{i}\\}_{i=1}^{M}$ be a partition of the phase
space in $M$ boxes of identical size, so that $\Omega=\bigcup_{i=1}^{M}B_{i}$.
Let $\chi_{B_{i}}$ be the characteristic function of $B_{i}\subset\Omega$. We
consider the projection $\pi_{M}:L^{1}(\Omega)\to
sp\\{\chi_{B_{1}},...,\chi_{B_{M}}\\}$ given by

$\displaystyle\pi_{M}f=\sum_{i=1}^{M}\left(\frac{1}{m(B_{i})}\int_{B_{i}}f\quad
dm\right)\chi_{B_{i}},$ (8)

with $m$ being the normalised Lebesgue measure on $\Omega$. The action of
$\mathcal{P}_{t}$ on $sp\\{\chi_{B_{1}},...,\chi_{B_{M}}\\}$ admits the
following matrix representation

$\displaystyle P_{i,j}^{\tau}=\frac{m(B_{i}\cap\Phi^{-\tau}B_{j})}{m(B_{i})}.$
(9)

By construction we have that $P^{\tau}$ is a stochastic matrix. It follows
that $P_{i,j}^{\tau}\geq 0$ and $\sum_{i=1}^{M}P_{i,j}^{\tau}=1$ $\forall
j=1,...,M$, since $P_{i,j}^{\tau}$ gives the probability of being in state $i$
at time $t+\tau$ being the system in state $j$ at time $\tau$.

We propose here a different way to discretise the dynamics of the system. We
select $M$ numerical UPOs $U_{1},...,U_{M}$ and we associate the states
$A_{1},...,A_{M}$ obtained by considering the UPOs together with their
neighbourhoods. Each $A_{i}$ represents one of the possible discrete states of
the system. We implement the shadowing algorithm: at each time step $t$ the
UPO $U_{k}$ that minimises the distance with the chaotic trajectory is
selected (See section II.1.4 for more details on the algorithm). Hence we say
that the system is in the state $A_{k}$ at time $t$. The stochastic variable
$s:\\{1,...,N_{chaotic}\cdot dt\\}\subset\mathbb{N}\to\mathcal{A}$ describes
the shadowing process just outlined as follows:

$\displaystyle s(t)=A_{k}$ (10)

with $A_{k}$ being the shadowing UPO at time t and corresponding
neighbourhood. We then construct the stochastic matrix along the lines of
equation 9 as

$\displaystyle
P_{i,j}^{dt}\approx\frac{\\#\\{k:(s(k)=A_{j})\wedge(s(k+dt)=A_{i})\\}}{N_{chaotic}}$
(11)

where $\\#$ defines the cardinality of the set.

###  III.2 Spectral Properties of the Transfer Operator

In this section we use the spectrum of the stochastic matrix $P^{dt}$ to study
the mixing properties of the system. We focus on the process of scattering
that the forward trajectory undergoes by being repelled continuously between
the neighbourhood of the various UPOs.

Let us recall a few basic properties of the spectrum of a general stochastic
matrix. Its leading eigenvalue is $\lambda=1$, and its corresponding
eigenvector $\textbf{w}^{(1)}$ , in the case of an ergodic Markov chain,
determines the unique invariant measure. The other eigenvalues, which can be
proven to be inside the unit circle, fulfill the condition
$\sum_{j}w_{j}^{(\lambda)}=0$. The subdominant eigenvalues , ordered
accordingly to
$1>\Re(\lambda_{2})\geq\Re(\lambda_{3})\geq...\geq\Re(\lambda_{M})$, can be
thought of as modes of decay, as they determine the time scale of convergence
to the stationary probability measure. We can quantify these time scales by
defining the corresponding decay rate as
$\tau_{k}=-\frac{dt}{log\left(\Re\left(\lambda_{k}\right)\right)}$, where $dt$
takes into account how we have discretised the dynamics in the time domain. In
particular, $\tau_{2}$ identifies the mixing time scale pikovsky_2016 .

We consider a set of $M$ UPOs and corresponding neighbourhoods chosen at
random from our database. We present results for $M=300$. The corresponding
plots for $M=100,1000$ are displayed in the supplementary material appendix
and indicate the robustness of our results. Choosing a sufficiently large $M$
makes sure that, with overwhelming probability, we have a selection of UPOs
that has a distribution of periods that is comparable with the one of the full
dataset. We remark that we are unable to use the full set of $\bar{M}$ UPOs
because it leads to substantial computational difficulties. We derive the
matrix $P^{dt}$ following the procedure outlined in Section III.1, by
considering the shadowing of a chaotic trajectory with length
$T_{max}=10^{5}$. $P^{dt}$ is a stochastic matrix by construction, its first
eigenvalues are
$\lambda_{1}=1,\lambda_{2}=0.9931,\lambda_{3}=0.9921,\lambda_{4}=0.9907$ and
the corresponding decay rates are
$\tau_{2}=1.4443,\tau_{3}=1.2608,\tau_{4}=1.0154$. We also verified that the
associated stochastic process is Markovian and that for each matrix
(corresponding to $M=100,300,100$ UPOs) there exist a value $\hat{N}$ so that
$P_{i,j}^{\hat{N}}\neq 0$ $\forall$ $i,j$, implying that the process is
ergodic.

We want to portray these modes in $\mathbb{R}^{3}$, in order to understand the
geometrical structure associated with the UPOs. We will proceed as follows. It
is known sparrow_1982  that the attractor of the L63 model is contained in a
compact subset of $\mathbb{R}^{3}$, being it the Cartesian product
$\mathcal{P}=[-20,20]\times[-27.5,27,5]\times[1,48]$. We cover this region
with $103400$ boxes $\mathcal{D}=\\{D_{i}\\}_{i=1}^{103400}$ of equal size
$1$. The boxes are built having adjacent sides, so that $\mathcal{D}$
constitutes a partition of $\mathcal{P}$. Each UPO and corresponding
neighbourhood intersects a certain number of boxes and each box might contain
contributions from different orbits. We would like to define a numerical
quantity (mass) that weighs the contribution given by UPOs of different types
within each box. We set a fixed number of points $\bar{N}$ to be represented
in the phase space a priori. We then assign the points to each box, depending
on the weight given by the different UPOs and corresponding neighbourhood from
$v^{(2)}$. We consider the transformation
$\hat{v}^{(2)}=|v^{(2)}|/||v^{(2)}||$. Each state $A_{i}$ will be awarded of a
number of points $\bar{n}_{i}$ given by
$\bar{n}_{i}=\hat{v}^{(2)}_{i}\cdot\bar{N}$. We will sample these points
uniformly along the UPO trajectory. Thus, each states of the system (i.e. each
periodic orbit and its neighbourhood) will be characterised by two quantities:
$s(v^{(2)}_{i})$, its ”type”, and $\bar{n}_{i}$, its weight. Let $\theta$ be a
map $\theta:\mathcal{D}\to\mathbb{R}$ that associates a numerical value to
each box, with the aim of quantifying the proportion of orbits of each kind.
We consider the box $D_{j}$. Let us call $pos_{j}$ and $neg_{j}$ the number of
positive and negative points contained in the box. We define
$\theta(D_{j})=pos_{j}+neg_{j}$, to quantify the mass contained in each cube.
The results of the aforementioned decomposition are presented in Figure 9.
Panel 9a describes the invariant measure, panels 9b, 9c, 9d describe the
subdominant eigenvalues. The latter show the geometrical structures describing
how an initial probability measure converges to the invariant one.

We wish to provide a somewhat alternative interpretation of the plots
corresponding to the subdominant eigenvalues as follows, keeping in mind that
UPOs are invariant sets and provide exact solutions of the evolution
equations. We propose that regions characterised by the same colour (red and
blue in our figures) are associated with separate bundles of UPOs. The forward
trajectory undergoes transitions between the neighbourhood of the UPOs
belonging to a bundle, and is repelled with low probability towards the
neighbourhood of an UPO belong to the other bundle. The closer to one the real
part of an eigenvalue, the less efficient is the exchange between regions of
different colours in the corresponding mode. We observe that the smaller the
eigenvalue, thus associated to faster decay rate, the finer the geometrical
structure associated with the mode. This agrees with our intuition on how
diffusion works.

(a)

(b)

(c)

(d)

Figure 9: Invariant Measure of the system obtained by projection of
$\lambda_{1}$ LABEL:sub@fig:_subfig_11a. Distribution of mass in the phase
space as a projection of: LABEL:sub@fig:_subfig_11b: $\lambda_{2}$,
LABEL:sub@fig:_subfig_11c: $\lambda_{3}$, LABEL:sub@fig:_subfig_11d:
$\lambda_{4}$.

###  III.3 Quasi-Invariant Sets

It is possible to provide further support to the interpretation given above to
the geometrical structures determined by the subdominant modes. In fact, these
resemble the quasi-invariant sets studied by Froyland and Padberg for L-63
froyland_2009 . Loosely speaking, quasi-invariant sets are macroscopic
dynamical structure such that the probability of individual trajectories
beginning in the subset would leave it in short time is very little (see
Appendix appendix  for more details). Let us explore the similarity more in
details by adapting the heuristic approach implemented in froyland_2005 ;
froyland_2008 ; froyland_2003  and originally proposed by Dellnitz and Junge
dellnitz_1999 ; dellnitz_1997  to detect quasi-invariant sets.

The adapted algorithm is set up as follows. We calculate the transition matrix
$P^{dt}$ with the methodology outlined in the previous paragraph. Let
$\textbf{w}^{(2)}$ be the eigenvector associated with $\lambda_{2}$. We
consider the mapping $s:\\{w^{(2)}_{1},...,w^{(2)}_{M}\\}\to\\{1,-1\\}$
provided by $\textbf{w}^{(2)}$:

$s(w^{(2)}_{i})=\begin{cases}-1&w^{(2)}_{i}\leq 0\\\
1&w^{(2)}_{i}>0\end{cases}$ (12)

where $w^{(2)}_{i}$ indicates the $i^{th}$ component of the eigenvector
$w^{(2)}$ corresponding to the state $i$. Note that this mapping closely
resembles what we have done in Section III.2 except that the weighting of the
mass is removed, whereas the sign of the component is kept. This allows us to
define two sets $B_{1}$ and $B_{2}$:

$\displaystyle B_{1}=\bigcup_{i\in\mathcal{I}_{1}}A_{i}\quad
where\quad\mathcal{I}_{1}=\\{i:s(w_{i}^{(2)})=1\\}$ (13) $\displaystyle
B_{2}=\bigcup_{i\in\mathcal{I}_{2}}A_{i}\quad
where\quad\mathcal{I}_{2}=\\{i:s(w_{i}^{(2)})=-1\\}$ (14)

We thus grouped the UPOs and corresponding neighbourhoods into two different
bundles as discussed in froyland_2009 . These two sets heuristically identify
quasi-invariant regions between which little communication is expected.

We test the heuristic procedure for three sets of states of different size.
Namely, we choose at random from our database a set of 100, 300 and 1000 UPOs.
We found robustness with respect to the size of the sample and thus present
the results only for 300 states. The other two cases can be found in the
supplementary material appendix .
For each set of UPOs and corresponding neighbourhoods we calculate the
transition matrix, its leading eigenvalues and the associated decay time. The
result are presented in Figure 10. The second eigenvector $\textbf{w}^{(2)}$
provides information on slowest rate at which the system is mixing, defining
the largest quasi-invariant macroscopic structures. We see a clear qualitative
correspondence between the patterns observed in Fig. 9b and 10a.

(a)

(b)

(c)

Figure 10:  Quasi-invariant bundles of UPOs obtained with the method outlined
in Section III.3. LABEL:sub@fig:_subfig_9a: $\lambda_{2}$=0.9931,
$\tau_{2}=1.4443$; LABEL:sub@fig:_subfig_9b: $\lambda_{3}$=0.9921,
$\tau_{3}=1.2608$; LABEL:sub@fig:_subfig_9c: $\lambda_{4}$=0.9902,
$\tau_{4}=1.0154$.

$B_{1}$ and $B_{2}$ \- corresponding to the alternating red and blue regions -
identify two quasi-invariant structures in the phase space between which we
expect little diffusion. We now take a step further and extract information
from the other subdominant eigenvalues $\lambda_{k}$ of the transfer operator.
More precisely, the subdominant eigenvectors $\textbf{w}^{(k)}$ provide an
ordering of the quasi-invariant structures in terms of ”leakiness” (See e.g.
froyland_2003 ; huisinga_2006 ; deuflhard_2005 ) . In fact, we observe that
also in this case, finer quasi-invariant structures are associated to
eigenvalues with smaller real part, which corresponds to faster diffusion
(Fig.10b, 10c) . This correspondence is interesting but not surprising,
because the physical process responsible for the slow decay described in Fig.
10 is indeed the slow mixing occurring in phase space between the regions
described by the quasi-invariant sets depicted in Fig. 9.

##  IV Conclusion

The theory of UPOs has found extensive applications in the study of low-
dimensional chaotic systems, in particular as a mean to calculate dynamical
averages through the use of trace formulas eckhardt_1994 ; franceschini_1993 ;
zoldi_1998 . In recent times promising developments has been made in the
context of higher dimension dynamical systems cvitanovic_2013 ; chandler_2013
; kawahara_2001 . Our work wishes to provide a contribution towards the goal
of using UPOs to interpret the dynamics of complex systems. The ranked
shadowing, intended as the representation of approximation of a chaotic
trajectory through the nearest UPOs, was at the center of our investigation
which has been performed on L63. We found that longer UPOs, despite being
underrepresented in our sample, are the most effective in shadowing the orbit
of the system. We then proposed a discrete Markovian interpretation of the
dynamics, defining the stochastic matrix controlling the shadowing of UPOs
along the trajectory. In this context each UPO represents a possible state of
the system. We found that this perspective on the dynamics allows for a new
interpretation of quasi-invariant structures in terms of UPOs. By looking at
the eigenvectors corresponding to subdominant eigenvalues we are able to
explain on the one side the geometrical pattern associated with the decay of
perturbations with respect to the invariant measure and on the other side link
this behaviour to the structure of the invariant sets. This work prepares the
path for different future lines of research. We would like to extend this
analysis to higher dimensional system of practical relevance, characterised by
high variability in the number of unstable dimension. In particular we would
like to extend the work of Lucarini and Gritsun lucarini_2020  on blocking
events, investigating transitions between zonal flow and blocking by applying
the methodology developed in this paper on the Lucarini-Gritsun UPOs dataset,
by investigating the dynamical properties between zonal and blocked states.

###### Acknowledgements.

The authors have benefitted from scientific exchanges with P. Cvitanović, J.
Dorrington, G. Ducci, G. Froyland, A. Gritsun, C. Nesbitt, M. Santos, N.
Zagli, M. Zaks. VL acknowledges the support received from the EPSRC project
EP/T018178/1 and from the EU Horizon 2020 project TiPES (grant no. 820970).
CCM has been supported by an EPSRC studentship as part of the Centre for
Doctoral Training in Mathematics of Planet Earth (grant number EP/L016613/1).

## References

  * [1] Henri Poincaré.  Les méthodes nouvelles de la mécanique céleste: Méthodes de MM. Newcomb, Glydén, Lindstedt et Bohlin. 1893, volume 2.  Gauthier-Villars it fils, 1893.
  * [2] Predrag Cvitanović.  Periodic orbits as the skeleton of classical and quantum chaos.  Physica D: Nonlinear Phenomena, 51(1-3):138–151, 1991.
  * [3] Celso Grebogi, Edward Ott, and James A Yorke.  Unstable periodic orbits and the dimensions of multifractal chaotic attractors.  Physical Review A, 37(5):1711, 1988.
  * [4] J-P Eckmann and David Ruelle.  Ergodic theory of chaos and strange attractors.  The theory of chaotic attractors, pages 273–312, 1985.
  * [5] Rufus Bowen.  $\omega$-limit sets for axiom a diffeomorphisms.  Journal of differential equations, 18(2):333–339, 1975.
  * [6] Martin C Gutzwiller.  Chaos in classical and quantum mechanics, volume 1.  Springer Science & Business Media, 2013.
  * [7] David Ruelle.  Thermodynamic formalism: the mathematical structure of equilibrium statistical mechanics.  Cambridge University Press, 2004.
  * [8] Bruno Eckhardt and Gerolf Ott.  Periodic orbit analysis of the lorenz attractor.  Zeitschrift für Physik B Condensed Matter, 93(2):259–266, 1994\.
  * [9] Predrag Cvitanović.  Dynamical averaging in terms of periodic orbits.  Physica D: Nonlinear Phenomena, 83(1-3):109–123, 1995.
  * [10] Roberto Artuso, Erik Aurell, and Predrag Cvitanovic.  Recycling of strange sets: I. cycle expansions.  Nonlinearity, 3(2):325, 1990.
  * [11] David Ruelle.  Smooth dynamics and new theoretical ideas in nonequilibrium statistical mechanics.  Journal of Statistical Physics, 95(1):393–468, 1999.
  * [12] Anatole Katok and Boris Hasselblatt.  Introduction to the modern theory of dynamical systems.  Number 54. Cambridge university press, 1997.
  * [13] Stephen Smale et al.  Differentiable dynamical systems.  Bulletin of the American mathematical Society, 73(6):747–817, 1967\.
  * [14] Rufus Bowen.  Periodic orbits for hyperbolic flows.  American Journal of Mathematics, 94(1):1–30, 1972.
  * [15] Giovanni Gallavotti.  Nonequilibrium and irreversibility.  Springer, 2014.
  * [16] Valerio Lucarini, Richard Blender, Corentin Herbert, Francesco Ragone, Salvatore Pascale, and Jeroen Wouters.  Mathematical and physical ideas for climate science.  Reviews of Geophysics, 52(4):809–859, 2014.
  * [17] Giovanni Gallavotti.  Chaotic dynamics, fluctuations, nonequilibrium ensembles.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 8(2):384–392, 1998.
  * [18] Giovanni Gallavotti and Ezechiel Godert David Cohen.  Dynamical ensembles in nonequilibrium statistical mechanics.  Physical review letters, 74(14):2694, 1995.
  * [19] Ditza Auerbach, Predrag Cvitanović, Jean-Pierre Eckmann, Gemunu Gunaratne, and Itamar Procaccia.  Exploring chaotic motion through periodic orbits.  Physical Review Letters, 58(23):2387, 1987.
  * [20] Predrag Cvitanović.  Invariant measurement of strange sets in terms of cycles.  Physical Review Letters, 61(24):2729, 1988.
  * [21] Roberto Artuso, Erik Aurell, and Predrag Cvitanovic.  Recycling of strange sets: Ii. applications.  Nonlinearity, 3(2):361, 1990.
  * [22] Edward N Lorenz.  Deterministic nonperiodic flow.  Journal of atmospheric sciences, 20(2):130–141, 1963.
  * [23] Valerio Lucarini.  Evidence of dispersion relations for the nonlinear response of the lorenz 63 system.  Journal of Statistical Physics, 134(2):381–400, 2009.
  * [24] Valerio Lucarini and Andrey Gritsun.  A new mathematical framework for atmospheric blocking events.  Climate Dynamics, 54(1):575–598, 2020.
  * [25] Predrag Cvitanović.  Recurrent flows: the clockwork behind turbulence.  Journal of Fluid Mechanics, 726:1–4, 2013.
  * [26] Genta Kawahara and Shigeo Kida.  Periodic motion embedded in plane couette turbulence: regeneration cycle and burst.  Journal of Fluid Mechanics, 449:291, 2001.
  * [27] Gary J Chandler and Rich R Kerswell.  Invariant recurrent solutions embedded in a turbulent two-dimensional kolmogorov flow.  Journal of Fluid Mechanics, 722:554–595, 2013.
  * [28] A Gritsun.  Statistical characteristics, circulation regimes and unstable periodic orbits of a barotropic atmospheric model.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1991):20120336, 2013.
  * [29] Predrag Cvitanovic, Roberto Artuso, Ronnie Mainieri, Gregor Tanner, Gábor Vattay, Niall Whelan, and Andreas Wirzba.  Chaos: classical and quantum.  ChaosBook. org (Niels Bohr Institute, Copenhagen 2005), 69:25, 2005\.
  * [30] Andrey Gritsun and Valerio Lucarini.  Fluctuations, response, and resonances in a simple atmospheric model.  Physica D: Nonlinear Phenomena, 349:62–76, 2017.
  * [31] AS Gritsun.  Unstable periodic trajectories of a barotropic model of the atmosphere.  Russian Journal of Numerical Analysis & Mathematical Modelling, 23(4), 2008.
  * [32] Davide Lasagna.  Sensitivity analysis of chaotic systems using unstable periodic orbits.  SIAM Journal on Applied Dynamical Systems, 17(1):547–580, 2018\.
  * [33] Davide Lasagna.  Sensitivity of long periodic orbits of chaotic systems.  Physical Review E, 102(5):052220, 2020.
  * [34] Chiara Cecilia Maiocchi and Valerio Lucarini.  https://figshare.com/projects/Decomposing_the_Dynamics_of_the_Lorenz_1963_model_using_Unstable_Periodic_Orbits_Averages_Transitions_and_Quasi-Invariant_Sets/119907
  * [35] Gary Froyland and Kathrin Padberg.  Almost-invariant sets and invariant manifolds—connecting probabilistic and geometric descriptions of coherent structures in flows.  Physica D: Nonlinear Phenomena, 238(16):1507–1523, 2009.
  * [36] Gary Froyland.  Statistically optimal almost-invariant sets.  Physica D: Nonlinear Phenomena, 200(3-4):205–219, 2005.
  * [37] Gary Froyland.  Extracting dynamical behavior via markov models.  In Nonlinear dynamics and statistics, pages 281–321. Springer, 2001\.
  * [38] Stanislaw M Ulam.  Problems in modern mathematics.  Courier Corporation, 2004.
  * [39] Warwick Tucker.  The lorenz attractor exists.  Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 328(12):1197–1202, 1999.
  * [40] Z Galias and P Zgliczyński.  Computer assisted proof of chaos in the lorenz equations.  Physica D: Nonlinear Phenomena, 115(3-4):165–188, 1998.
  * [41] Warwick Tucker.  A rigorous ode solver and smale’s 14th problem.  Foundations of Computational Mathematics, 2(1):53–117, 2002.
  * [42] Arkady Pikovsky and Antonio Politi.  Lyapunov exponents: a tool to explore complex dynamics.  Cambridge University Press, 2016.
  * [43] Edward Ott.  Chaos in dynamical systems.  Cambridge university press, 2002.
  * [44] Valter Franceschini, Claudio Giberti, and Zhiming Zheng.  Characterization of the lorentz attractor by unstable periodic orbits.  Nonlinearity, 6(2):251, 1993.
  * [45] Scott M Zoldi.  Unstable periodic orbit analysis of histograms of chaotic time series.  Physical review letters, 81(16):3375, 1998.
  * [46] Yoshitaka Saiki and Michio Yamada.  Reply to “comment on ‘time-averaged properties of unstable periodic orbits and chaotic orbits in ordinary differential equation systems’”.  Physical Review E, 81(1):018202, 2010.
  * [47] Yoshitaka Saiki and Michio Yamada.  Time-averaged properties of unstable periodic orbits and chaotic orbits in ordinary differential equation systems.  Physical Review E, 79(1):015201, 2009.
  * [48] Michael A Zaks and Denis S Goldobin.  Comment on “time-averaged properties of unstable periodic orbits and chaotic orbits in ordinary differential equation systems”.  Physical Review E, 81(1):018201, 2010.
  * [49] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn.  Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory.  Meccanica, 15(1):9–20, 1980.
  * [50] T.N. Palmer.  Extended-range atmospheric prediction and the lorenz model.  Bulletin of the American Meteorological Society, 74(1):49–65, 1993\.
  * [51] Rufus Bowen.  Topological entropy and axiom a.  In Proc. Sympos. Pure Math, volume 14, pages 23–41, 1970.
  * [52] Y Saiki.  Numerical detection of unstable periodic orbits in continuous-time dynamical systems with chaotic behaviors.  Nonlinear Processes in Geophysics, 14(5):615–620, 2007.
  * [53] Roberto Barrio, Angeles Dena, and Warwick Tucker.  A database of rigorous and high-precision periodic orbits of the lorenz model.  Computer Physics Communications, 194:76–83, 2015.
  * [54] Zbigniew Galias and Warwick Tucker.  Symbolic dynamics based method for rigorous study of the existence of short cycles for chaotic systems.  In 2009 IEEE International Symposium on Circuits and Systems, pages 1907–1910. IEEE, 2009.
  * [55] Viviane Baladi.  Positive transfer operators and decay of correlations, volume 16.  World scientific, 2000.
  * [56] C Sparrow.  The lorenz equations.  Edited by Arun V. Holden, page 111, 1982.
  * [57] Gary Froyland.  Unwrapping eigenfunctions to discover the geometry of almost-invariant sets in hyperbolic maps.  Physica D: Nonlinear Phenomena, 237(6):840–853, 2008.
  * [58] Gary Froyland and Michael Dellnitz.  Detecting and locating near-optimal almost-invariant sets and cycles.  SIAM Journal on Scientific Computing, 24(6):1839–1863, 2003.
  * [59] Michael Dellnitz and Oliver Junge.  On the approximation of complicated dynamical behavior.  SIAM Journal on Numerical Analysis, 36(2):491–515, 1999.
  * [60] Michael Dellnitz and Oliver Junge.  Almost invariant sets in chua’s circuit.  International Journal of Bifurcation and Chaos, 7(11):2475–2485, 1997.
  * [61] Wilhelm Huisinga and Bernd Schmidt.  Metastability and dominant eigenvalues of transfer operators.  In New algorithms for macromolecular simulation, pages 167–182. Springer, 2006.
  * [62] Peter Deuflhard and Marcus Weber.  Robust perron cluster analysis in conformation dynamics.  Linear algebra and its applications, 398:161–184, 2005.