Review article: Dynamical systems, algebraic topology and the climate sciences

The definition of climate itself cannot be given without a proper understanding of the key ideas of long-term behavior of a system, as provided by dynamical systems theory. Hence, it is not surprising that concepts and methods of this theory have percolated into the climate sciences as early as the 1960s. The major increase in public awareness of the socio-economic threats and opportunities of climate change has led more recently to two major developments in the climate sciences: (i) the Intergovernmental Panel on Climate Change's successive Assessment Reports and (ii) an increasing understanding of the interplay between natural climate variability and anthropogenically driven climate change. Both of these developments have benefited from remarkable technological advances in computing resources, relating throughput as well as storage, and in observational capabilities, regarding both platforms and instruments.

Starting with the early contributions of nonlinear dynamics to the climate sciences, we review here the more recent contributions of (a) the theory of non-autonomous and random dynamical systems to an understanding of the interplay between natural variability and anthropogenic climate change and (b) the role of algebraic topology in shedding additional light on this interplay. The review is thus a trip leading from the applications of classical bifurcation theory to multiple possible climates to the tipping points associated with transitions from one type of climatic behavior to another in the presence of time-dependent forcing, deterministic as well as stochastic.

1 Introduction and motivation

This paper is based on the invited talks given by the two authors in an online series on “Perspectives on climate sciences: From historical developments to research frontiers”. The series had twice-monthly talks from July 2020 to July 2021 and its success led to the idea of having a special issue of Nonlinear Processes in Geophysics. The talks of the two co-authors are available at https://youtu.be/xjccOfptYII (last access: 27 September 2023) (Michael Ghil) and https://youtu.be/W1yndTsvR0g (last access: 27 September 2023) (Denisse Sciamarella). In the present paper, we go beyond the lively but more perishable video version to what we hope is a more coherent and permanent record of the convergence between two strains of Henri Poincaré's heritage – dynamical systems theory (Poincaré, 1892, 2017) and algebraic topology (Poincaré, 1895; Siersma, 2012) – and their joint applications to the climate sciences. This convergence resulted from the two authors meeting in November 2018 at the University of Buenos Aires, where Michael Ghil gave a series of six lectures on “Mathematical Problems in Climate Dynamics” at the invitation of Denisse Sciamarella; see Ghil (2021a, b, c, d, e).

1.1 Dynamical systems and climate dynamics

Many of the ideas and methods of dynamical systems theory were introduced into the climate sciences by a generation of pioneers in the 1960s. Stommel (1961) formulated a two-box and two-pipe model for the oceans' overturning circulation that had two stable, steady-state solutions with counter-rotating flows. Veronis (1963) used another form of reduced-order model, by projecting a one-layer, single-gyre wind-driven circulation model onto a small number of Fourier modes and likewise found multiple solutions, both steady and periodic. Most intriguingly, Lorenz (1963a) applied the same type of low-truncation Galerkin method to a Boussinesq model of flow between two horizontal plates, heated from below and cooled from above. In this setting, transitions from a quiescent fluid to two mutually symmetric flow patterns and to chaotic solutions occurred. The latter, particularly simple, convection model governed by only three ordinary differential equations (ODEs) provided inspiration for hundreds of papers on deterministically chaotic phenomena in the climate sciences and way beyond.

None of the pioneering papers mentioned above, though, nor any of the thousands of papers since, exhibits all the phenomena – mathematical and physical – of interest in this review paper. As we proceed, the illustrative examples will be taken from atmospheric, oceanographic and climate models that capture best one or a few of these phenomena.

It is important to realize that Poincaré had already seen the analogy between the chaos he found in the so-called reduced three-body problem of celestial mechanics (Poincaré, 1892; Gray, 2013; Poincaré, 2017) and the “sensitive dependence on initial conditions” that he realized occurred in the evolution of the weather. In fact, in Book I, Chap. IV of Poincaré (1908), called “Le Hasard”, he states that “it may happen that small differences in the initial conditions produce very great ones in the final phenomena”. And his second example of sensitive dependence is weather:

Our second example will be very analogous to the first and we shall take it from meteorology. Why have the meteorologists such difficulty in predicting the weather with any certainty? Why do the rains, the tempests themselves seem to us to come by chance, so that many persons find it quite natural to pray for rain or shine, when they would think it ridiculous to pray for an eclipse? We see that great perturbations generally happen in regions where the atmosphere is in unstable equilibrium. The meteorologists are aware that this equilibrium is unstable and that a cyclone is arising somewhere; but where they can not tell; one-tenth of a degree more or less at any point, and the cyclone bursts here and not there, and spreads its ravages over countries it would have spared. This we could have foreseen if we had known that tenth of a degree, but the observations were neither sufficiently close nor sufficiently precise, and for this reason all seems due to the agency of chance. Here again we find the same contrast between a very slight cause, unappreciable to the observer, and important effects, which are sometimes tremendous disasters.

The translations are both from Poincaré (2003, Book I, Chap. IV, “Chance”).

The work of Lorenz (1963a), while actually referring to a highly simplified model of thermal convection, illustrates perfectly Poincaré's insights about the role of what we now call deterministic chaos rather than pure chance. Ghil and Childress (1987) presented the applications of dynamical systems theory to large-scale atmospheric and climate dynamics as well as to dynamo theory and geomagnetism, in a systematic book form; they included a gradual introduction to the basic mathematical concepts and tools involved, and this book was reissued by Springer in 2012 as an e-book. Dijkstra (2013) provided a considerably expanded version of Ghil and Childress (1987), in terms of both the mathematical content and the areas of climatic applications, which include oceanographic and coupled ocean–atmosphere phenomena.

Aside from its applications to the climate sciences, the dynamical systems literature is quite extensive, in covering both mathematical fundamentals and applications to other areas. Holmes (2007) provides a fine historical overview of the field from 1885 to 1965, along with a fairly complete bibliography. Two important books are Arnol'd (2012) and Guckenheimer and Holmes (1983). Further references – including some that treat the subject for infinite-dimensional function spaces, like those that describe the solutions of the partial differential equations of fluid dynamics – are given in the subsequent list of noteworthy insights that dynamical systems theory has provided for the climate sciences.

Following Ghil et al. (1991) and Ghil (2019), we summarize herewith some key insights in the climate sciences from the theory of autonomous dynamical systems, in which time-dependent forcing or coefficients are absent.

1. The equations of continuum mechanics are nonlinear. Surprisingly many phenomena can be explained by linearization about a particular fixed basic state. Many more cannot.

2. Behavior of solutions to nonlinear equations – subject to some reasonable mathematical assumptions – changes qualitatively only at isolated points in phase-parameter space, called bifurcation points. Behavior along a single branch of solutions, between such points, is modified only quantitatively and can be explored by linearization about the basic state, which changes as the parameters change. That is, nonlinear dynamics are much like linear dynamics, only more so (Lorenz, 1963a, b; Ghil and Childress, 1987).

3. Bifurcation trees lead from the simplest, most symmetric states to highly complex and realistic ones, with much lower symmetry in either space or time or both. These trees can be explored partially by analytic methods (Jin and Ghil, 1990; Jordan and Smith, 2007) and more fully by numerical ones, such as pseudo-arclength continuation (Legras and Ghil, 1985; Dijkstra, 2005).

4. The truly nonlinear behavior near bifurcation points involves robust transitions, of great generality, between single and multiple fixed points (saddle–node, pitchfork and transcritical bifurcations), fixed points and limit cycles (Hopf bifurcation), and limit cycles and strange attractors (“routes to chaos”: Eckmann, 1981; Guckenheimer and Holmes, 1983). As the complexity of the behavior increases, its predictability decreases (Ghil, 2001).

5. Behavior in the most realistic, chaotic regime can be described by the ergodic theory of dynamical systems. In this regime, statistical information similar to but more detailed than for truly random behavior can be extracted and used for predictive purposes (Eckmann and Ruelle, 1985; Mo and Ghil, 1987; Ghil and Robertson, 2000).

6. Chaos and strange attractors are not restricted to low-order systems. They can be shown to exist for the full equations governing continuum mechanics (Constantin et al., 1989; Temam, 2000). The detailed exploration of finite- but high-dimensional attractors is in full swing (Legras and Ghil, 1985; Dijkstra, 2005; Simonnet et al., 2009; Doedel and Tuckerman, 2012; Dijkstra et al., 2014).

7. Single time series (Takens, 1981) and single numbers derived from them (e.g., Grassberger, 1983) have been used to describe chaotic behavior. This very simple and straightforward use of a nonlinear concept has attracted considerable attention to deterministically chaotic dynamics, including in the geosciences (Nicolis and Nicolis, 1984; Tsonis and Elsner, 1988). The use of single time series, while exciting in theory, is not very promising when the series are short and noisy (Ruelle, 1990; Smith, 1988). The increasing availability of a large number of similar series at different points in space, combined with physical insight, is compensating more and more for the shortcomings of each individual time series in describing the complexity of many phenomena in the geosciences, as well as advancing their prediction (Ghil et al., 2002).

Further details on the contributions of autonomous dynamical systems theory in general and the concepts and methods of bifurcation theory in particular, appear in Sect. 2.1. The recent contributions of the theory of non-autonomous and random dynamical systems (NDSs and RDSs) – with their generalization of bifurcations to tipping points – are reviewed in Sect. 2.2.

1.2 Algebraic topology and chaotic dynamics

What is the topology of chaos, and why is it important in the theory of dynamical systems and in the time series analysis for nonlinear and chaotic dynamics? We attempt here to provide answers to these questions, with an emphasis on applications to the climate sciences. Essentially, the concepts and tools of algebraic topology can be applied to the evolution of systems in both phase space and physical space as well as to the interesting back-and-forth trip between the two spaces. This complementary view of the way that dynamics and topology interact is a main motivation of the present article.

The emphasis on time dependence and dynamics here should not allow us to forget, though, the huge role that homologies have already been playing in the fields of image processing and visualization (e.g., Heine et al., 2016; Singh Bansal et al., 2022, and references therein). Significant advances in computational topology (Edelsbrunner and Harer, 2022) have helped substantially in these more static area of applications and will clearly do so in the more dynamic ones contemplated herein.

In Sect. 3.1, we present the rather novel approach of Branched Manifold Analysis through Homologies (BraMAH) (Sciamarella and Mindlin, 2001; Charó et al., 2021b) for approximating the branched manifolds (Birman and Williams, 1983a, b) of dynamical systems by a cell complex that allows one to characterize the manifold by its homology groups in phase space (Poincaré, 1895; Sciamarella and Mindlin, 1999). The detection and description of localized coherent sets (LCSs) in two-dimensional flows in physical space by BraMAH-based methods is reviewed in Sect. 3.2.

The most recent developments of the merging of the two strands of Poincaré's heritage – algebraic topology and dynamical systems – are covered in Sect. 3.3 and 3.4. In Sect. 3.3, we introduce the templex, a novel concept in algebraic topology (Charó et al., 2022, 2023), which complements the previously mentioned cell complexes of BraMAH by a directed graph (digraph), whose nodes are the cells and which approximates the flow on the branched manifold. The extension of this concept to the noise-perturbed chaotic attractors of RDS theory follows in Sect. 3.4.

In the rest of this section, we provide some quick historical background to the current interest in the ways in which algebraic topology can help one infer a system's chaotic dynamics from one or more time series of its observables. The first methods of time series analysis that associated geometric properties with experimental time series appeared in the early 1980s (e.g., Packard et al., 1980). These geometric methods continue to be used, for instance, to analyze datasets of Lagrangian trajectories and understand the geometry of transport (Banisch and Koltai, 2017).

But is geometry the best lens one can use to classify data according to underlying differences in dynamics? Classifying dynamics is possible thanks to invariants or quasi-invariants in phase space. Gilmore (1998) classified invariants as belonging to three distinct categories:

• metric invariants, such as dimensions of various types, e.g., correlation dimension (Grassberger and Procaccia, 1983) or multifractal scaling functions (Halsey et al., 1986);

• dynamic invariants, such as Lyapunov exponents (Oseledec, 1968; Wolf et al., 1985), further discussed by Eckmann and Ruelle (1985) and Abarbanel and Kennel (1993); and finally

• topological invariants, linking numbers, relative rotation rates, Conway polynomials and branched manifolds (Williams, 1974).

The first two kinds of invariants do not provide information on how to model the system's dynamics, while topological invariants actually do. Why is this so? Topology deals with the properties of a geometric object that do not change when continuous deformations are performed. Stretching, twisting, crumpling or bending preserve topology; while cutting or suturing holes, gluing separated pieces, or producing self-crossings do not. Volumes in phase space can be stretched or squeezed, folded or torn. The particular manner in which these processes are combined repetitively in phase space leads to a structure. The topology of such a structure is the signature of the mechanisms acting to build certain dynamics.

The “recipe” to “knead” the Lorenz (1963a) strange attractor is illustrated in Fig. 1 as a sequence of steps that are topological in nature. Quoting Gilmore and Lefranc (2003), “sets of initial conditions (cubes) are sliced, by running into an axis with a stable and unstable direction (the z axis for Lorenz-like systems), for example. The different parts flow off in different directions in the phase space, where they may encounter other sliced parts from different regions of the phase space. These are squeezed together and eventually return to regions they originated from (recursion).”

Figure 1 Sketch of the topological processes that intervene in obtaining the strange attractor of the Lorenz (1963a) convection model. From Letellier and Gilmore (2013) with permission by World Scientific Publishing Corp.

The advantage of using topology instead of geometry or fractality to describe chaos lies in the fact that topology provides information about the elementary stretching, folding, tearing or squeezing mechanisms that act in phase space to shape the flow. Geometric features may differ, but if the underlying dynamics obey certain equivalence principles, the topology should be the same. Topological equivalence between branched manifolds is defined by isotopy. In other words, two objects are isotopic if it is possible to mold one into the other without tearing or gluing it. It is in this sense that we speak of dynamical equivalence. Different geometric deformations of the Lorenz attractor that preserve its topology are sketched in Fig. 2. There is a two-way correspondence between topology and dynamics, in a sense that will be clarified in Sect. 3.

Figure 2 Point clouds associated with different geometrical representations of the Lorenz (1963a) attractor. They are obtained by integrating the model's governing equations using coordinate transformations for some of the variables. The butterfly is deformed, but the topological structure of the butterfly is maintained.

A good starting point for this quick historical perspective is the pioneering paper of Henri Poincaré (Poincaré, 1895), who first described the way in which a dynamical system's properties depend upon its topology; see also Gray (2013, Chap. 8). The concept of a branched manifold, introduced by Williams (1974), was anticipated in Edward N. Lorenz's famous convection paper: on Lorenz (1963a, p. 138), he remarks that “the [computed] trajectory is confined to a pair of surfaces which appear to merge in the lower portion of Fig. 3.” The paper's Fig. 3 is reproduced here, coincidentally, as Fig. 3 as well. Lorenz plots the isopleths of X as a function of Y and Z of the strange attractor, to approximate surfaces formed by all points on limiting trajectories. The etymology of “isopleth” combines “iso” with the ancient Greek word plêthos, “a great number”, as in the modern English word “plethora”. It is generically used to refer to a curve of points sharing the same value of some quantity. We will return to a stochastically perturbed version of the Lorenz (1963a) model in Sects. 2.2 and 3.4.

Figure 3 Isopleths of X (thin solid curves) as a function of Y and Z, based on a single trajectory of length 6000 time steps of the Lorenz (1963a) attractor, and isopleths of the lower of the two values of X where two values occur (dashed curves) for approximate surfaces formed by all points on nearby trajectories. The heavy solid curve and the extension of the dotted curves indicate natural boundaries of the surfaces. From Lorenz (1963a), published in 1963 by the American Meteorological Society.

Joan Birman and Robert F. Williams used branched manifolds to classify chaotic attractors in terms of the way periodic orbits are “knotted” in dynamical systems (Birman and Williams, 1983a, b). These authors discovered that systems whose branched manifolds have the same topology, are dynamically equivalent. From this discovery, a topologist's dream blooms: can one classify types of dynamics as one classifies the elements in Mendeleev's table?

In the late 1990s, it became possible to determine whether or not two three-dimensional (3-D) dissipative dynamical systems are equivalent by using knot theory (Gilmore, 1998; Gilmore and Lefranc, 2003; Natiello et al., 2007; Letellier and Gilmore, 2013). In Sect. 3.1, we address the question of how these authors and many more worked with knots, the difficulties that arose with knot theory, and how the latter were solved, at least in part, using the homology groups of algebraic topology. Section 3.2–3.4 describe the applications of BraMAH to the Lagrangian analysis of fluid flows in physical space, the introduction of digraphs to complement cell complexes in describing the flow on a branched manifold in phase space and the extension of the templexes that describe the latter to noise-driven chaotic systems.

2 Dynamical systems theory for the climate sciences

As we indicated in Sect. 1.1, dynamical systems theory entered the evolution of the climate sciences – at that time consisting mainly of meteorology and oceanography – in the 1960s, in the pioneering papers of Lorenz (1963a, b), Veronis (1963), Stommel (1961) and others. These and other papers of the 1960s and early 1970s did not necessarily include explicit references to bifurcation theory, although awareness of the fundamental concepts and methods was clearly present in one form or another. In the 1970s, another set of papers, on energy balance models (EBMs), reported on the possibility of alternative stable steady states (warm and cold) of Earth's climate system (Held and Suarez, 1974a; North, 1975). Ghil (1976a) specifically introduced the saddle–node bifurcation into this climate setting as well as numerical methods needed to deal with it in the context of a full partial differential equation model of climate rather than of low-order or otherwise simplified models.

The contributions of “nonlinear dynamics”, as dynamical systems theory tended to be referred to by physicists and other non-mathematicians by training, were presented for the first time in a quadrennial report (1987–1991) of the US geosciences community to the International Union of Geodesy and Geophysics (IUGG) by Ghil et al. (1991). The presentation of elementary bifurcations below is for a broad audience and is based on Boers et al. (2022). It focuses on multistability and the possible transitions between different regimes of behavior: in Sect. 2.1 for systems with time-independent forcing and coefficients and in Sect. 2.2 for systems in which time dependence is present in either the forcing or the coefficients or both.

2.1 Autonomous dynamical systems

Assume that the state of a system of interest can be described by a vector x∈ℝ^d and that the time evolution of x(t) is governed by the following equation of motion, namely a first-order autonomous ODE:

$$\begin{array}{}\text{(1)}& \dot{\mathit{x}}=\mathit{f}(\mathit{x};p).\end{array}$$

Here $(\cdot {)}^{\cdot }=\mathrm{d}(\cdot )/\mathrm{d}t$, f denotes a generally nonlinear, smooth – i.e., continuously differentiable, up to some order – vector field and p a scalar parameter or, in more general cases, a small set of parameters. For clarity, one separates the variables x from the parameter p by a semicolon. The term “autonomous” refers to the fact that in Eq. (1) both the coefficients and the forcing are constant in time. This means that changes in p are assumed to be infinitely slow or, at least, very slow compared to the characteristic internal-variability times of the system being modeled.

Points x^* for which $\mathit{f}({\mathit{x}}^{*};p)=\mathrm{0}$ are called fixed points. Linearizing the equation of motion around a given fixed point x^* yields, for a small perturbation $\stackrel{\mathrm{̃}}{\mathit{x}}=\mathit{x}-{\mathit{x}}^{*}$,

$$\begin{array}{}\text{(2)}& \dot{\stackrel{\mathrm{̃}}{\mathit{x}}}={\mathit{f}}^{\prime }({\mathit{x}}^{*};p)\stackrel{\mathrm{̃}}{\mathit{x}};\end{array}$$

here ${\mathit{f}}^{\prime }({\mathit{x}}^{*};p)$ is the Jacobian matrix comprised of the elements $\partial f_i/\partial x_j$. For an initial condition ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{0}}$, the solution to this linearized equation is given by

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{̃}}{\mathit{x}}\left(t\right)=e^{f^{\prime }t}{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{0}}.\end{array}$$

We call a fixed point x^* linearly stable if all eigenvalues of f^′ have negative real part and linearly unstable otherwise. A scalar example will be given in Sect. 2.1.1.

The bifurcations of a dynamical system that we deal with in this subsection describe the creation and annihilation of fixed points as well as changes in their linear stability. Further types of bifurcations are considered in the next subsection.

Typically, bifurcations lead to abrupt qualitative changes in the dynamics, explaining why they are often invoked as a mathematical model for abrupt regime shifts or state transitions in real-world systems. Until fairly recently, bifurcations were studied mostly in the context of autonomous dynamical systems. The more realistic situations in which the forcing is allowed to depend explicitly on time are addressed in Sect. 2.2. In this broader context, bifurcations have been called “tippings” in the climate sciences (Lenton et al., 2008; Ashwin et al., 2012; Kuehn, 2011; Ghil, 2019) and elsewhere.

There are at least two different interpretations of “tipping” and “tipping points” in the literature. One of these, emanating from Gladwell (2000) and Lenton et al. (2008), interprets tipping merely as a sudden change, whether due to a well-defined bifurcation or not. In this interpretation, a tipping point is merely a threshold. The other interpretation sees a tipping point as a generalization to non-autonomous systems of a bifurcation point (Kuehn, 2011; Ghil, 2019).

In the latter case, tipping is necessarily related to a tipping point in phase-parameter space as opposed to just a threshold in some parameter value; thus, not every jump or critical transition arises from a such a point. Both points of view – pun intended, of course – have their merits, but confusion should be avoided to the extent possible. Clearly, in this review article, we follow the more unambiguously defined mathematical version.

2.1.1 The double-well potential as a source of bistability

As an instructive and widely used example, we briefly introduce a prototype model to describe scalar dynamical systems than can occupy either one of two stable fixed points, separated by an unstable one, as plotted in Fig. 4a. The double-well potential $U(x;p)=x^{\mathrm{4}}/\mathrm{4}-x^{\mathrm{2}}/\mathrm{2}-px$ leads to the equation of motion

$$\begin{array}{}\text{(4)}& \dot{x}=U^{\prime }(x;p)=-x^{\mathrm{3}}+x+p,\end{array}$$

where $U^{\prime }=\partial U/\partial x$. For $p<-p^{*}$ this dynamical system only has the stable fixed point $x_{-}^{*}$ and for $p>p^{*}$ it only has the stable fixed point $x_{+}^{*}$, while for $-p^{*}<p<p^{*}$ the two stable fixed points $x_{\pm }^{*}$ coexist and there is a third, unstable fixed point $x_{\mathrm{0}}^{*}$ in between these two.

Figure 4 Bifurcation diagrams for (a) the double-fold, (b) supercritical pitchfork and (c) supercritical Hopf bifurcation. Stable fixed-point branches are indicated by solid lines and unstable ones by dashed lines. The colored lines in panel (a) correspond to a hysteresis cycle. See text for details. After Boers et al. (2022) under CC-BY license.

The two stable fixed points correspond to the two minima of the potential U above, whereas $\pm p^{*}$ represents the two critical thresholds of the system. In this scalar case, the basins of attraction of the two minima are the intervals $-\mathrm{\infty }<x<x_{\mathrm{0}}^{*}$ and $x_{\mathrm{0}}^{*}<x<+\mathrm{\infty }$, respectively. They are separated by the unstable fixed point $x_{\mathrm{0}}^{*}$, which is a local maximum of the potential U(x).

Changing p slowly from, say, $p=-\mathrm{1}$ to $p=+\mathrm{1}$ will lead to a bifurcation-induced critical transition from $x_{-}^{*}$ to $x_{+}^{*}$ at the critical value $p=p^{*}$. When p is subsequently changed back from p=1 to $p=-\mathrm{1}$, the transition from $x_{+}^{*}$ back to $x_{-}^{*}$ will only occur at $p=-p^{*}$. This phenomenon of jumps from one fixed point to the other occurring at distinct parameter values is called hysteresis, and it is highly relevant to the practical reversibility of abrupt transitions. It was studied, for instance, in electromagnetic systems by James Clerk Maxwell and by Pierre Curie, and it is important in the physical, biomedical, engineering and socio-economic sciences. In the context of the climate sciences, a hysteresis loop like the one seen in Fig. 4a has been described in detail for EBMs by Ghil and Childress (1987, Chap. 10) and by Ghil (1994), using solar insolation as the parameter p.

The bifurcation introduced above is called a double-fold bifurcation, since it is obtained by combining a supercritical fold (with the stable branch reaching forward to p→∞) with a subcritical one (with the stable branch reaching backward to $p\to -\mathrm{\infty }$). A more recent version of such a double-fold bifurcation is plotted in Fig. 2 of Von der Heydt et al. (2016) for an energy balance model with respect to carbon dioxide concentration as the parameter p. The single-fold bifurcation is often called a saddle–node bifurcation (super- or subcritical) since in two dimensions it corresponds to the merging of a node that is stable in both directions on one branch with a saddle that is stable in one direction and unstable in the perpendicular direction on the other branch; see, for instance, Ghil and Childress (1987, Fig. 12.3) for sketches of the stability of fixed points for a linear autonomous ODE system in two dimensions.

2.1.2 Bistability in the presence of symmetry: the pitchfork bifurcation

Another example of bistability is given by a pitchfork bifurcation (Fig. 4b). Its so-called normal form, i.e., the simplest ODE that exhibits the change in behavior of interest, is

$$\begin{array}{}\text{(5)}& \dot{x}=x(p-x^{\mathrm{2}}).\end{array}$$

This bifurcation captures bistable behavior in systems in which spatial mirror symmetry prevails for low p values. A well-known example in the climate sciences is symmetry in a meridional plane for an idealized Atlantic Meridional Overturning Circulation (AMOC) at low buoyancy forcing by a weak pole-to-Equator temperature and precipitation gradient (Quon and Ghil, 1992; Ghil, 2001; Dijkstra, 2005). In this case, water sinks at both poles and rises on either side of the Equator, forming two overturning cells that are symmetric with respect to the equatorial plane.

The solutions of Eq. (5) are x₀=0 and $x_{\pm }=\pm \sqrt{p}$. In this normal form, the scalar symmetry of the latter two solutions with respect to 0 stands for the mirror symmetry of the AMOC's overturning cells with respect to the Equator.

The bifurcation occurs as the parameter p, which is a normalized form of the thermal and salinity forcing in the AMOC case and crosses over from negative to positive values. It is easy to check that, for p≤0, x₀=0 is the unique fixed point, while for p>0, the three fixed points coexist. Their linear stability is given by considering infinitesimal perturbations around a given steady-state solution $x=x_{\ast }+\mathit{\xi }$.

With the scalar version $f=f(x;p)$ of the notation in Eq. (1), we have the scalar version of Eq. (2) in the specific case at hand given by

$$\begin{array}{rl}(x_{\ast }+\mathit{\xi }{)}^{\cdot }& =f(x_{\ast }+\mathit{\xi };p)\\ & =\mathrm{0}+{\displaystyle \frac{\partial f(x;p)}{\partial x}}{|}_{x=x_{\ast }}\mathit{\xi }+\mathcal{O}\left({\mathit{\xi }}^{\mathrm{2}}\right)\simeq p-\mathrm{3}{x}_{\ast }^{\mathrm{2}}.\end{array}$$

Since ${\dot{x}}_{\ast }=\mathrm{0},$ this leaves

$$\begin{array}{}\text{(6)}& \dot{\mathit{\xi }}=p-\mathrm{3}{x}_{\ast }^{\mathrm{2}}\end{array}$$

to determine linear stability for small ξ (Ghil and Childress, 1987; Dijkstra and Ghil, 2005). Thus it is clear that, for p<0, the unique solution $x_{\ast }=\mathrm{0}$ is linearly stable; but, for p>0, this null solution becomes linearly unstable, while the two mutually symmetric solutions $x_{\ast }=x_{\pm }=\pm \sqrt{p}$ are stable, since $p-\mathrm{3}p<\mathrm{0}$. We thus suspect that, for sufficiently strong buoyancy forcing, the two-cell AMOC will lose its stability and yield the approximately single-cell AMOC that is currently observed; see Stocker and Wright (1991) or Quon and Ghil (1992), for instance.

2.1.3 Beyond bistability: Hopf bifurcation and limit cycles

Bistability is only the first step up the bifurcation tree that leads from system behavior with the highest degree of symmetry in space and time – possibly as simple as uniform in both – to behavior that has greater and greater complexity (Eckmann, 1981; Ghil and Childress, 1987; Strogatz, 2018). We outline now one further step up this tree, the one leading from fixed points to stable periodic solutions, called limit cycles in dynamical systems parlance.

In polar coordinates, this normal form is given by

$$\begin{array}{}\text{(7)}& \dot{\mathit{\rho }}=\mathit{\rho }(\mathit{\mu }-{\mathit{\rho }}^{\mathrm{2}}),\dot{\mathit{\theta }}=\mathit{\upsilon },\end{array}$$

with $\mathit{\upsilon }=-\mathrm{1}$ for uniform counterclockwise rotation around the origin. The two equations above are decoupled and the one for the radial variable $\mathit{\rho }=x^{\mathrm{2}}+y^{\mathrm{2}}$ has exactly the same form as the pitchfork normal form for x in Eq. (5). Note, though, that here ρ is necessarily nonnegative and the mirror symmetry of Fig. 4b is replaced by the rotational symmetry of Fig. 4c.

The version shown in Fig. 4 is the supercritical one, which leads to a smooth increase with μ in the amplitude of an oscillation generated by the Hopf bifurcation. For plots of the subcritical Hopf bifurcation, please see Ghil and Childress (1987, Figs. 12.8 and 12.9). In particular, in the presence of higher-order terms, as in Fig. 12.9b of Ghil and Childress (1987), a sharp jump from no oscillation to a finite-amplitude one occurs as μ passes a critical threshold, and one can have a hysteresis cycle between no oscillation and a large-amplitude oscillation. For instance, it is a matter of some debate whether the Mid-Pleistocene Transition – during which both the amplitude and the dominant periodicity of climatic variability changed – might be associated with a sub- or a supercritical Hopf bifurcation; see Riechers et al. (2022, and references therein).

2.1.4 Successive bifurcations and routes to chaos

A further step on the route to chaos for deterministic systems with no explicit time dependence (Eckmann, 1981; Ghil and Childress, 1987; Strogatz, 2018) involves the transition from a one-dimensional limit cycle to a two-dimensional torus in phase space. In the latter case, the motion on the torus is quasi-periodic – i.e., the coordinates of the point on the torus are of the form $(x,y)=\left(x\right(t)=f({\mathit{\upsilon }}_{\mathrm{1}}t,{\mathit{\upsilon }}_{\mathrm{2}}t),y(t)=g({\mathit{\upsilon }}_{\mathrm{1}}t,{\mathit{\upsilon }}_{\mathrm{2}}t\left)\right)$, where the functions f(s,t) and g(s,t) are arbitrary and the two angular frequencies υ₁ and υ₂ are incommensurable; i.e., ${\mathit{\upsilon }}_{\mathrm{1}}/{\mathit{\upsilon }}_{\mathrm{2}}$ is not a rational number.

This kind of motion is typical in celestial mechanics (Arnold et al., 2007; Ghil and Childress, 1987), and, in fact, the periodicities that are associated with the orbital forcing of the glacial–interglacial cycles are of this type, although one usually refers to them by truncated values – such as those in Table 12.1 of Ghil and Childress (1987) – that could suggest that the ratios between these periodicities, like 41 kyr for the obliquity and 19 kyr for the precessional parameter, do have a common denominator. The latter view is clearly an oversimplification but this is not the place to discuss chaos in the solar system, whether Hamiltonian or, more recently, dissipative.

Quasi-periodic motion already looks much more irregular than purely periodic motion. Thus, for instance, the intervals between lunar or solar eclipses are highly irregular. Still, the 14th century scholar Nicole Oresme was already aware of the kinematic consequences of quasi-periodicity for celestial motions (Grant, 1961). He realized that a periodic and a quasi-periodic motion cannot be distinguished from each other during a finite observation interval. Oresme also knew that the motion of a point on a torus will describe a simple closed loop if the two angular velocities are commensurable, while the point's orbit will never close but densely cover the surface of the torus if the two velocities are incommensurable (Arnol'd, 2012), i.e., in a way that is visually indistinguishable from painting the whole torus a uniform color.

From quasi-periodic motion to a deterministically chaotic one there are several routes (Eckmann, 1981), as already mentioned in Sect. 1.1. Some of these routes to chaos were explored numerically in the climate sciences by Lorenz (1963a, b) and described more didactically in Chaps. V and VI of Ghil and Childress (1987) for atmospheric motions. For such routes in the paleoclimatic context, see Ghil and Childress (1987, Chap. XII) as well as Ghil (1994). We shall not go into greater detail herein but pass instead to the more recent insights from the theory of dynamical systems subject to time-dependent forcing.

2.2 Non-autonomous and random dynamical systems

Realistically, the natural systems that we want to describe in terms of dynamical systems theory are non-autonomous, meaning that f in Eq. (1) above has an explicit time dependence: $\partial \mathit{f}/\partial t\not\equiv \mathbf{0}$. The Earth system as a whole, as well as all its components, is clearly non-autonomous, being affected by time-dependent forcing, such as quasi-periodic variations in solar insolation due to gravitational perturbations in Earth's orbit (Milankovitch, 1920), along with anthropogenic forcing due to rising greenhouse gas concentrations (Arrhenius, 1896; Houghton et al., 1990; Solomon, 2007; IPCC, 2014, 2021).

Moreover, there is typically high-frequency forcing, such as cloud processes or weather variability. In a drastic simplification, this type of forcing is often represented by white noise (Hasselmann, 1976). Including both deterministic and stochastic time dependence requires a description of the dynamics in terms of stochastic differential equations of the form

$$\begin{array}{}\text{(8)}& \mathrm{d}\mathbf{X}=\mathbf{F}(\mathbf{X},t;p)\mathrm{d}t+\mathit{\sigma }\left(\mathbf{X}\right)\mathrm{d}\mathit{\eta },\end{array}$$

where dη denotes the infinitesimal increments of a Wiener process, which are stationary and independently distributed according to a normal distribution with mean μ=0 and variance $\mathbb{E}\left(\right|\mathrm{d}\mathit{\eta }{|}^{\mathrm{2}})=\mathrm{d}t$. Often, a further simplification is made in assuming that the noise is additive or state-independent, and thus σ=const. above. The possibly time-dependent but still deterministic term $\mathbf{F}(\mathbf{X},t;p)$ is called the drift.

Interest in autonomous dynamical systems and their bifurcations started over 2 centuries ago and can be traced back to Leonhard Euler and the Bernoullis, while that in non-autonomous and random dynamical systems (NDSs and RDSs) only goes back a few decades. We describe some key differences between the two cases next and justify the need for considering pullback attractors (PBAs) in the latter case.

2.2.1 NDSs, RDSs and pullback attraction

For the sake of simplicity, we assume that the physical system under consideration is described by a set of ODEs. In the autonomous case, such a set of ODEs can be formally written as

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}}=\mathbf{F}\left(\mathbf{X}\right),\mathbf{X}\left(t_{\mathrm{0}}\right)={\mathbf{X}}_{\mathrm{0}};\end{array}$$

here t∈ℝ, X∈ℝ^d, $\mathbf{F}:{\mathbb{R}}^d\to {\mathbb{R}}^d$, and d is the number of the system's dependent variables.

For the non-autonomous case, the brief presentation here follows Caraballo and Han (2017) and the paradigmatic formulation of the initial-value problem is

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}}=\mathbf{G}(t,\mathit{x}),\mathbf{X}\left(t_{\mathrm{0}}\right)={\mathbf{X}}_{\mathrm{0}}.\end{array}$$

As in Eq. (9), t∈ℝ, X∈ℝ^d and $\mathbf{G}:\mathbb{R}\times {\mathbb{R}}^d\to {\mathbb{R}}^d$ in Eq. (10), and one still assumes that G has “nice” properties that guarantee the existence, uniqueness and continuous dependence on initial states and on parameters for the solutions of Eq. (10). Furthermore, Caraballo and Han (2017) show that, provided the vector field G(t,X) is dissipative, solutions of Eq. (10) exist and satisfy the two other properties globally, i.e., for all t∈ℝ. We call such a global solution $\mathit{\phi }(t,t_{\mathrm{0}},{\mathbf{X}}_{\mathrm{0}})$.

There are two key distinctions between the autonomous case and the non-autonomous one:

• a.

  In the autonomous setting, solutions cannot intersect, since there is only one trajectory through a given point X₀∈ℝ^d due to uniqueness. Hence, for d=2, the only possible (forward) attracting sets are fixed points and limit cycles; i.e., chaotic behavior and strange attractors can only occur for d≥3. The NDS setting is different in these respects; i.e., intersections are possible at 2 times t₁ and t₂≠t₁, and thus chaos can occur for d=2 and periodic forcing, as is the case, for instance, in the Van der Pol oscillator (e.g., Guckenheimer and Holmes, 1983).

• b.

  In the autonomous setting, solutions depend only on the time t−t₀ elapsed since initial time, while in the NDS setting, they depend separately on the initial time t₀ and the current time t, at which we observe the system. In the former setting, it suffices to consider forward-in-time attraction, which results in attractors that are fixed; time-independent objects, such as fixed points; limit cycles; tori; and strange attractors. In the latter case, we need to define pullback attraction and the PBAs that it leads to.

Before proceeding with a more rigorous justification for and definition of a PBA, here it is, in the simplest possible terms: a pullback attractor is a possibly time-dependent object in a system's phase space that exhibits attraction in the sense of convergence at each time t to a set, called a snapshot, to which the system's initial state at time s tends as s tends to $-\mathrm{\infty }.$ This is distinct from the forward attractors that can be defined for autonomous systems started at a fixed time t₀.

Given the uniqueness and the continuous dependence of the global solutions to Eq. (10) on initial states and on parameters, it is straightforward to verify that a global solution φ of Eq. (10) satisfies

• i.

  the initial value property at t=t₀, namely $\mathit{\phi }(t_{\mathrm{0}},t_{\mathrm{0}},{\mathbf{X}}_{\mathrm{0}})={\mathbf{X}}_{\mathrm{0}}$, and

• ii.

  the two-parameter semigroup evolution property,

  $$\mathit{\phi }(t_{\mathrm{2}},t_{\mathrm{0}},{\mathbf{X}}_{\mathrm{0}})=\mathit{\phi }(t_{\mathrm{2}},t_{\mathrm{1}},\mathit{\phi }(t_{\mathrm{1}},t_{\mathrm{0}},{\mathbf{X}}_{\mathrm{0}}\left)\right)\mathrm{for}{t}_{\mathrm{0}}\le t_{\mathrm{1}}\le t_{\mathrm{2}},$$

  which corresponds to the concatenation of solutions; i.e., in order to go from t₀ to t₂ one can go first from t₀ to t₁ and then from t₁ to t₂.

One can then provide the following definition of a process.

Definition 1. Let ${\mathbb{R}}_{+}^{\mathrm{2}}=\mathit{\left\{}\right(t,t_{\mathrm{0}})\in {\mathbb{R}}^{\mathrm{2}}:t\ge t_{\mathrm{0}}\mathit{\}}$. A process on ℝ^d is a family of mappings

$$\mathit{\phi }(t,t_{\mathrm{0}},\cdot )):{\mathbb{R}}^d\to {\mathbb{R}}^d,(t,t_{\mathrm{0}})\in {\mathbb{R}}_{+}^{\mathrm{2}},$$

which satisfy

• i.

  the initial value property $\mathit{\phi }(t_{\mathrm{0}},t_{\mathrm{0}},\mathbf{X})=\mathbf{X}$ for all X∈ℝ^d and any t₀∈ℝ,

• ii.

  the two-parameter semigroup property for all X∈ℝ^d and both $(t_{\mathrm{2}},t_{\mathrm{1}})\in {\mathbb{R}}_{+}^{\mathrm{2}}$ and $(t_{\mathrm{1}},t_{\mathrm{0}})\in {\mathbb{R}}_{+}^{\mathrm{2}}$, and

• iii.

  the continuity property that the mapping $(t,t_{\mathrm{0}},\mathbf{X})\mapsto \mathit{\phi }(t,t_{\mathrm{0}},\mathbf{X})$ be continuous on ${\mathbb{R}}_{+}^{\mathrm{2}}\times {\mathbb{R}}^d$.

An alternative NDS formulation for this process formulation is the so-called skew-product formulation, which goes back to the work of George R. Sell, as reviewed in Sell (1971); see also Kloeden and Yang (2020). A process as defined above is also called a two-parameter semigroup on ℝ^d, in contrast with the one-parameter semigroup of an autonomous dynamical system, since the former depends not just on the initial time t₀, as in the latter case, but also on the current time t.

This difference matters, in particular, in determining the asymptotic behavior of the solutions. In the autonomous case, a global solution is invariant with respect to translation in time: $\mathit{\phi }(t,t_{\mathrm{0}},{\mathbf{X}}_{\mathrm{0}})=\mathit{\phi }(t-t_{\mathrm{0}},\mathrm{0},{\mathbf{X}}_{\mathrm{0}})$. Hence, the usual forward asymptotic behavior for $t\to +\mathrm{\infty }$ and fixed t₀ is the same as the behavior for fixed t and $t_{\mathrm{0}}\to -\mathrm{\infty }$. This equivalence may no longer hold when the translation invariance is lost, as it is in the NDS case.

To illustrate the effect of this lost invariance, consider the following simple scalar ODE (Caraballo and Han, 2017, Sect. 3.2.1):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}=-ax+b\sin t,x\left(t_{\mathrm{0}}\right)=x_{\mathrm{0}},t\ge t_{\mathrm{0}},\end{array}$$

for which analytical computations can be carried out explicitly. Individual solutions do not have a forward limit as $t\to +\mathrm{\infty }$ for fixed t₀, but the difference between any two solutions vanishes in this limit. The particular solution

$$\begin{array}{}\text{(12)}& \mathcal{A}\left(t\right)={\displaystyle \frac{b(a\sin t-\cos t)}{a^{\mathrm{2}}+\mathrm{1}}}\end{array}$$

provides the long-term information on the behavior of all the solutions of Eq. (11). This result is best captured by recognizing that the pullback limit,

$$\begin{array}{}\text{(13)}& \underset{t_{\mathrm{0}}\to -\mathrm{\infty }}{lim}\mathit{\phi }(t,t_{\mathrm{0}},x_{\mathrm{0}})=\mathcal{A}\left(t\right)\mathrm{for}\mathrm{all}t\mathrm{and}{x}_{\mathrm{0}}\in \mathbb{R},\end{array}$$

yields 𝒜(t) as the PBA of all the solutions of Eq. (11).

One is thus led to the following rigorous definition of a PBA for a forced dissipative dynamical system subject to a time-dependent forcing, where we have generalized ℝ^d to a finite-dimensional metric space 𝒳 and have replaced t₀ by s, for greater symmetry.

Definition 2. A PBA 𝒜 is an indexed family of invariant sets (𝒜(t))_t∈R that depend on time and satisfy the following conditions:

1. For all t, 𝒜(t) is a compact subset in 𝒳 that is invariant with respect to the two-parameter semigroup ℱ(t,s),

   $$\begin{array}{}\text{(14)}& \mathcal{F}(t,s)\mathcal{A}\left(s\right)=\mathcal{A}\left(t\right)\mathrm{for}\mathrm{every}s\le t,\mathrm{and}\end{array}$$

2. for all t, pullback attraction is reached when

   $$\begin{array}{}\text{(15)}& \underset{s\to -\mathrm{\infty }}{lim}{D}_{\mathrm{H}}\left(\mathcal{F}\right(t,s)\mathcal{B},\mathcal{A}(t\left)\right)=\mathrm{0}\text{for all}\mathcal{B}\in \mathcal{C},\end{array}$$

   where D_H(E,D) is the Hausdorff semi-distance between two sets and 𝒞 is a collection of bounded sets in 𝒳.

In the physical literature, the invariant sets 𝒜(t) at a given t∈ℝ have been called snapshots (Romeiras et al., 1990) and this terminology has been used also in the recent climate literature on the applications of NDSs, RDSs and PBAs (Ghil et al., 2008; Chekroun et al., 2011; Tél et al., 2020).

The finite-dimensional definition above follows Charó et al. (2021b, Appendix A and references therein). In fact, both deterministic and stochastic versions of forcing have been applied, for instance, by Chekroun et al. (2018) in the study of an infinite-dimensional, delay-differential equation model of the El Niño–Southern Oscillation (ENSO). The deterministic forcing corresponded to the purely periodic, seasonal changes in insolation, while the stochastic component represented the westerly wind bursts appearing in various ENSO models by Fei-Fei Jin and Axel Timmermann (e.g., Timmermann and Jin, 2002); see also Chekroun et al. (2011, Sect. 4.3). This ENSO example, among many others, shows that there is great flexibility in the application of the concepts and methods of non-autonomous dynamical systems (NDS and RDS) theory to climate problems.

2.3 Simple examples of pullback and random attractors

A straight-line PBA. An even simpler example of a PBA than the one of Eqs. (11) and (12) above is given by

$$\begin{array}{}\text{(16)}& \dot{x}=-\mathit{\alpha }x+\mathit{\sigma }t,\end{array}$$

with both α and σ being positive. The example was provided by Mickaël D. Chekroun (personal communication, 2011). The autonomous part of this ODE, $\dot{x}=-\mathit{\alpha }x$, is dissipative, and all solutions $x(t;x_{\mathrm{0}})=x(t;x(\mathrm{0})=x_{\mathrm{0}})$ converge to 0 as $t\to +\mathrm{\infty }$. What about the non-autonomous, forced ODE?

Here, the time-dependent forcing σt and the state-dependent dissipation −αx will tend to balance. But, again, as in the example of Sect. 2.2.1, there is no forward limit as $t\to +\mathrm{\infty }$, and one has to use the pullback limit, i.e., replace x(t;x₀) by $x(s,t;x_{\mathrm{0}})=x(t;x(s)=x_{\mathrm{0}})$ and let $s\to -\mathrm{\infty }$. Doing so yields the snapshots

$$\begin{array}{}\text{(17)}& \mathcal{A}\left(t\right)={\displaystyle \frac{\mathit{\sigma }}{\mathit{\alpha }}}\left(t-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\right).\end{array}$$

These snapshots are, in the extremely simple case at hand, just the points along the straight line illustrated in Fig. 5, which is the graph of the PBA (𝒜(t))_t∈ℜ.

Figure 5 The graph of the PBA for the simple NDS example governed by Eq. (16) and given by the indexed family ${\left(\frac{\mathit{\sigma }}{\mathit{\alpha }}(t-\frac{\mathrm{1}}{\mathit{\alpha }})\right)}_{t\in \mathrm{\Re }}$, along with several trajectories that converge to it from times $s_{\mathrm{2}}<s_{\mathrm{1}}<t_{\mathrm{1}}<t_{\mathrm{2}}$. Here t₁<t₂ are the increasing times at which we observe the system, while s₂<s₁ are the decreasing times to which we have to pull back in order to get the convergence. Figure courtesy of Mickaël D. Chekroun.

A PBA with periodic forcing. To further improve the reader's intuition for PBAs, we provide a second illustrative example here. It was worked out in detail by Riechers et al. (2022).

A system defined in polar coordinates by

$$\begin{array}{}\text{(18)}& \dot{\mathit{\rho }}=\mathit{\alpha }(\mathit{\mu }-\mathit{\rho }),\dot{\mathit{\varphi }}=\mathit{\upsilon },\text{with}\mathit{\rho },\mathit{\mu }\in {\mathbb{R}}^{+}\text{and}\mathit{\varphi }\in \mathbb{R}/\mathrm{2}\mathit{\pi },\end{array}$$

can easily be seen to exhibit a limit cycle in the (x,y) plane with $(x=\mathit{\rho }\cos \mathit{\varphi },y=\mathit{\rho }\sin \mathit{\varphi })$. An initial deviation of ρ from μ will decay exponentially, and the system converges to an oscillation of radius μ with the angular velocity υ. Here, we transform this autonomous dynamical system into a non-autonomous one by modulating the target radius μ with a sinusoidal forcing

$$\begin{array}{}\text{(19)}& \mathit{\mu }\to \mathit{\mu }\left(t\right)=\mathit{\mu }+\mathit{\beta }\sin \left(\mathit{\nu }t\right),\end{array}$$

where the modulation is moderate, so as to guarantee that $\mathit{\mu }+\mathit{\beta }\sin \left(\mathit{\nu }t\right)>\mathrm{0}$ for all t.

Since the dynamics of the phase ϕ and of the radius ρ are decoupled, the corresponding equations can be solved and analyzed separately. While the temporal development of the phase is trivial, the pullback invariant attracting set of the radius for the initial condition ρ(s)=ρ₀ is given by

$$\begin{array}{}\text{(20a)}& {\mathcal{A}}^{\left(\mathit{\rho }\right)}(t;{\mathit{\rho }}_{\mathrm{0}})=\underset{s\to -\mathrm{\infty }}{lim}\mathit{\rho }(t,s;{\mathit{\rho }}_{\mathrm{0}})=\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+\mathit{\mu },\end{array}$$

with

$$\begin{array}{}\text{(20b)}& \mathit{\vartheta }=\arctan (-\mathit{\nu }/\mathit{\alpha }),\end{array}$$

as shown in Riechers et al. (2022, Appendix B). Note that, in the limit $s\to -\mathrm{\infty }$, the dependence on the initial value ρ₀ vanishes and the attracting set ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ performs an oscillation of the same frequency as the forcing. It lags the phase of the time-dependent fixed point by the constant ϑ, while its amplitude is amplified by the factor α. Since ρ is restricted to positive values, this solution requires αβ<μ.

The PBA with respect to the coordinate ρ is comprised of the family of all the sets ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ as defined in Eqs. (20a) and (20b) and thus reads

$$\begin{array}{}\text{(21)}& {\mathcal{A}}^{\left(\mathit{\rho }\right)}=\mathit{\left\{}\mathit{\alpha }\mathit{\beta }\sin \right(\mathit{\nu }t+\mathit{\vartheta })+\mathit{\mu }{\mathit{\}}}_{t\in R}.\end{array}$$

Since the pullback limit for the phase ϕ does not exist, no constraints on it other than $\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi })$ are imposed by the dynamics. Hence, for the system (18) comprised of radius and phase, we find that

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \underset{t_{\mathrm{0}}\to -\mathrm{\infty }}{lim}{d}_{\text{H}}\left(\right(\mathit{\rho }(t;t_{\mathrm{0}},{\mathit{\rho }}_{\mathrm{0}}),\mathit{\varphi }(t;t_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})),\\ & \underset{{\mathcal{A}}_t}{\underbrace{\mathit{\left\{}\right(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+\mathit{\mu },\mathit{\varphi }):\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}}})=\mathrm{0},\end{array}\end{array}$$

where d_H denotes the Hausdorff semi-distance. The pullback-attracting sets 𝒜_t at time t are circles in the (x,y) plane with oscillating radius, and the system's PBA is given by the family of these circles:

$$\begin{array}{}\text{(23)}& \mathcal{A}=\mathit{\left\{}\right(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+\mathit{\mu },\mathit{\varphi }):\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }){\mathit{\}}}_{t\in R}.\end{array}$$

Figure 6 Trajectories and PBA of the system defined by Eqs. (18)–(20b). (a) Trajectories (ρ(t),ϕ(t)) of the system starting from different times in the past in the 3-D space spanned by the two Cartesian coordinates (x,y) and time t; the system's PBA lies on the red-shaded surface. (b) Heat map of the three trajectories' projection onto the (x,y) plane. A video of the heat map filling up, as more and more trajectories with different initial conditions are added, is provided in the Supplementary Material to Riechers et al. (2022). (c) Temporal evolution of the phase. (d) Temporal evolution of the radius (solid colored lines) together with its PBA (dashed red line). From Riechers et al. (2022) with thanks to the coauthors Keno Riechers, Takashita Mitsui and Niklas Boers.

Figure 6 shows trajectories of the system starting from different points in the past. In panel a the trajectories are depicted in the 3-D space spanned by the two Cartesian coordinates (x,y) and the time t, where the usual transformation from polar to Cartesian coordinates was applied. The shaded surface in this panel represents the PBA of the system. Panel b shows a heat map (Wilkinson and Friendly, 2009) that approximates a portion of the PBA's invariant measure projected onto the (x,y) plane. For a clean definition of such a measure in NDSs and RDSs, there are several references (e.g., Ghil et al., 2008; Chekroun et al., 2011; Caraballo and Han, 2017; Kloeden and Yang, 2020). Essentially, the heat map here counts the number of times that the trajectories in panel a cross small pixels in the (x,y) plane.

Note that the structure of the system's trajectories depends on the ratio $\mathit{\upsilon }/\mathit{\nu }$, and three different cases must be distinguished. If the radius is modulated with the same frequency as the oscillation itself, i.e.,  υ=ν, after one period the system practically repeats its orbit. More precisely, the radius of the oscillation does differ from one “round trip” to the next, but this difference tends to zero as ρ(t) asymptotically approaches ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}.$ If υ and ν are rationally related, mυ=nν with $n,m\in \mathbb{N}$, then the same quasi-repetition of the orbit occurs after n periods of the radial modulation and m periods of the system's oscillation. Such a trajectory will appear as an n-fold quasi-closed loop. Finally, if $\mathit{\upsilon }/\mathit{\nu }\notin \mathbf{Z}$, then the trajectory does not repeat itself but instead densely covers the annular disk $\mathcal{D}=\mathit{\left\{}\right(\mathit{\rho },\mathit{\varphi }):\mathit{\rho }\in [\mathit{\mu }-\mathit{\alpha }\mathit{\beta },\mathit{\mu }+\mathit{\alpha }\mathit{\beta }]$ and $\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi })\mathit{\}}$. The trivial evolution of the phase is depicted in panel c, while the trajectories of ρ(t) and their convergence to ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ are shown in panel d.

2.3.1 Random attractors (RAs)

Let us return now to the more general, nonlinear and stochastic case of Eq. (8) that includes not only deterministic time dependence F(X,t) but also random forcing:

$$\begin{array}{}\text{(24)}& \text{d}\mathbf{X}=\mathbf{F}(\mathbf{X},t)\text{d}t+\mathbf{G}\left(\mathbf{X}\right)\text{d}\mathit{\eta };\end{array}$$

here $\mathit{\eta }=\mathit{\eta }(t,\mathit{\omega })$ represents a Wiener process, which is taken to be scalar; its independent increments dη are commonly referred to as “white noise”, and ω labels the particular realization of this random process. More generally, one can also deal with a vector Wiener process, as in Eq. (8). See, for instance, Wax (1954) for early references on these matters.

When G= const. the noise is additive, while for $\partial \mathbf{G}/\partial \mathbf{X}\ne \mathbf{0}$, we speak of multiplicative noise. Intuitively, the distinction between dt and dη in the stochastic differential Eq. (24) is necessary since, roughly speaking and following the Einstein (1905) paper on Brownian motion, it is the variance of a Wiener process that is proportional to time and thus $\text{d}\mathit{\eta }\propto (\text{d}t{)}^{\mathrm{1}/\mathrm{2}}$. In Eq. (24), for the sake of simplicity, we also dropped the dependence on a parameter p that we had introduced, for the sake of generality, in Eq. (8).

The noise processes may include “weather” and volcanic eruptions when X(t) is “climate,” thus generalizing the linear model of Hasselmann (1976), or cloud processes when we are dealing with the weather itself: one person's signal is another person's noise, as the saying goes. In the case of random forcing, the concepts introduced by the simple example of Eq. (8) above can be illustrated by the random attractor 𝒜(ω) in Fig. 7.

Figure 7 Schematic diagram of a random attractor 𝒜(ω) and of the pullback attraction to it; here ω labels the particular realization of the random process θ(t)ω that drives the system. We illustrate the evolution in time t of the random process θ(t)ω (light solid black line at the bottom), the random attractor 𝒜(ω) itself (yellow band in the middle) with the snapshots ${\mathcal{A}}_{\mathrm{0}}\left(\mathit{\omega }\right)=\mathcal{A}(\mathit{\omega };t=\mathrm{0})$ and 𝒜(ω;t) (the two vertical sections, heavy solid), and the flow of an arbitrary compact set ℬ from “pullback times” $s=-{\mathit{\tau }}_{\mathrm{2}}$ and $s=-{\mathit{\tau }}_{\mathrm{1}}$ onto the attractor (heavy blue bands). See Appendix A in Ghil et al. (2008) for the requisite properties of the random process θ(t)ω that drives the RDS. After Ghil et al. (2008) with permission from Elsevier.

Chekroun et al. (2011) studied a specific case of such a random attractor for the paradigmatic, climate-related Lorenz (1963a) convection model. The authors introduced multiplicative noise into each of the ODEs of the original, deterministically chaotic system, as shown below:

$$\begin{array}{}\text{(25a)}& {\displaystyle }\text{d}X=P_r(Y-X)\text{d}t+\mathit{\sigma }X\text{d}\mathit{\eta },\text{(25b)}& {\displaystyle }\text{d}Y=(rX-Y-XZ)\text{d}t+\mathit{\sigma }Y\text{d}\mathit{\eta },\text{(25c)}& {\displaystyle }\text{d}Z=(-bZ+XY)\text{d}t+\mathit{\sigma }Z\text{d}\mathit{\eta };\end{array}$$

here r=28, P_r=10 and $b=\mathrm{8}/\mathrm{3}$ are the standard parameter values for chaotic behavior in the absence of noise and σ is a constant variance of the Wiener process that is not necessarily small. The well-known strange attractor of the deterministic case is replaced by the Lorenz model's random attractor, dubbed LORA by the authors by the authors of Chekroun et al. (2011).

Four snapshots 𝒜_t(ω) of LORA are plotted in Fig. 8 here, and a video of its evolution in time $\mathcal{A}\left(\mathit{\omega }\right)=\mathit{\left\{}{\mathcal{A}}_t\right(\mathit{\omega }){\mathit{\}}}_{t\in R}$ is available as Supplementary Material in Chekroun et al. (2011). What is actually plotted, in both the figure reproduced here and in the video, is the approximation of the time-dependent invariant measure ν_t(ω) supported by the attractor. A full definition of the sample measures of random attractors would occupy too much space in an already rather long review paper; please see Chekroun et al. (2011, Appendix A) and Charó et al. (2023, Appendix C), along with the references therein.

Figure 8 Heatmaps of the time-dependent invariant measure ν_t(ω) supported by four snapshots 𝒜_t(ω) of LORA. The values of the parameters r,s and b are the classical ones, while the variance of the noise is σ=0.5. The color bar, shown in Chekroun et al. (2011, Fig. 2) for a single snapshot, is on a log scale, and it quantifies the probability of landing in a particular region of phase space; shown is a projection of the 3-D phase space $(X,Y,Z)$ onto the (X,Z) plane. Note the complex, interlaced filament structures between highly populated regions (in yellow) and moderately populated ones (in red); the less populated a small patch, the darker its color. The time interval between the snapshots shown (left to right and top to bottom) is Δt=0.0875 in the nondimensional time units of the deterministic Lorenz (1963a) model. Reproduced from Chekroun et al. (2011, Fig. 3) with permission from Elsevier.

The striking effects of the noise on the nonlinear dynamics that are visible in Fig. 8 here and in the video of Chekroun et al. (2011) motivated much of the work reviewed in Sect. 3 below, starting with LORA's topological study by Charó et al. (2021b). The latter study gathered further insights into the abrupt changes in the snapshots' topology at critical points in time, changes that suggested the possibility of random processes giving rise to qualitative jumps in climate variability.

2.3.2 Abrupt transitions in non-autonomous systems

Ashwin et al. (2012) have proposed three classes of abrupt transitions in systems that can be described by Eq. (24): (i) bifurcation-induced transitions, (ii) noise-induced transitions and (iii) rate-induced transitions. An example of the first class has already been given in Sect. 2.1 and Fig. 4a above.

For an example of the second class, assume that the control parameter p remains constant in the drift term of Eq. (8), which is taken again to correspond to a double-well potential, as in Eq. (4). Noise-induced transitions occur when the noise amplitude is sufficiently high for the system to switch occasionally, and unpredictably, from one potential well to the other. Moreover, when p varies so as to push the system toward a bifurcation point, the noise will cause it to transition before – and, in certain cases, long before – the critical parameter value p^* of the corresponding deterministic system is reached.

Finally, the third class of rate-induced transitions arises when there is no strong separation between the system's intrinsic timescales and those at which the control parameter changes. So far, we implicitly assumed that, for each change in p, the system has sufficient time to adapt to the new equilibrium position; this type of slow change in p is sometimes called quasi-adiabatic. If this is not the case, the fixed point attracting the system may change its position so quickly that the system cannot follow and eventually loses track of the basin of attraction in which it started and falls into the other one (Ashwin et al., 2012).

Figure 9 Sketch of a double-fold bifurcation and how it leads to abrupt transitions and hysteresis in the temporal evolution of a system in a double-well potential with slowly changing parameter p=p(ϵt), where ϵ≪1, driven by additive white noise. The stable branch of fixed points is indicated by the solid part of the red line and the unstable one by the dashed part of the red line. Compare with Fig. 4a. After Boers et al. (2022).

Ashwin et al. (2012) have called these transitions tippings and refer to the three types described above as B tipping, N tipping and R tipping. Thus, aside from the rhetorically striking character of tipping points, tippings are the mathematically well-defined generalization of the bifurcations treated in the autonomous dynamical systems of Sect. 2.1 to the non-autonomous and random setting addressed herein. In fact, the first two types, B and N tipping, are not totally novel inasmuch as they only add deeper insight to what happens when a parameter p changes at a slow but finite (rather than infinitely slow rate). The biggest surprises occur for R tipping (Wieczorek et al., 2011; Feudel et al., 2018; Ghil, 2019; Pierini and Ghil, 2021), but we will not deal explicitly with this form of tipping herein.

We illustrate in Fig. 9 the N tipping of a system governed by $\mathrm{d}x=U^{\prime }(x;p)\mathrm{d}t+\mathit{\sigma }\mathrm{d}W$, with U as in Eq. (4) and dW as in Eq. (24), but σ= const. For example, simple EBMs (Ghil and Lucarini, 2020) exhibit a double-fold bifurcation of this kind, as described already in Sect. 2.1 above. The upper stable branch corresponds in this case to the current climate state, while the lower one corresponds to the Snowball Earth state (Held and Suarez, 1974b; Ghil, 1976b; Ghil and Lucarini, 2020).

To simulate the system's trajectory, the control parameter p is varied slowly from +1 to −1 and back to +1, causing the system to transition first from the upper stable branch to the lower one and then, at a considerably higher p value, back to the upper stable branch. Note that due to the noise driving the system, transitions typically occur earlier than expected from the corresponding deterministic dynamics governed by Eq. (4).

Note also that in the generalization from autonomous bifurcations to non-autonomous tippings, the phrase “tipping point” – aside from its threatening implication – is somewhat misleading: a bifurcation point is a point in phase-parameter space, like $(\pm x^{\ast },\pm p^{\ast })$ for the double well of Eq. (4) and Fig. 4a. The meaning attached to it by Gladwell (2000) in general and by Lenton et al. (2008) in the climate sciences refers only to the value of the forcing, like $\pm p^{\ast }$ in the case above.

3 Topological structure of flows in phase space and in physical space

At the end of Sect. 1.2, we mentioned that knot theory provided a first approach towards unveiling the topological structure of a flow in a 3-D phase space. In this case, the term “flow” does not refer to a fluid flow in physical space but to a family of solution curves of ODEs or other evolution equations (Arnol'd, 2012; Coddington and Levinson, 1955; Guckenheimer and Holmes, 1983). Of course, a flow in phase space may – as we will see later in Sect. 3.2 – refer to a particle in the Lagrangian description of a fluid flow in physical space. There is a strong link between the two situations, but the keywords refer to different motivations and objectives.

Clarifying the difference between these two kinds of flow, in physical space and in phase space, is relevant here because, in the community involved in the work been reviewed here, the phrase “topological chaos” is used when studying how fluid–particle trajectories are entangled in physical space during a mixing experiment. A noteworthy example is the motion induced by spatially periodic obstacles in a two-dimensional flow in order to form nontrivial braids (Gouillart et al., 2006; Thiffeault and Finn, 2006b), as shown in Fig. 10. Such motion generates exponential stretching of material lines and hence efficient mixing.

Figure 10 Topological chaos emerges in stirring or mixing experiments. Here we see stylized streamlines induced by pairs of rods on a periodic lattice and we see how these streamlines are stretched in physical space. From Thiffeault and Finn (2006a), under CC-BY license.

On the other hand, “topology of chaos” or “chaos topology”, for short, considers the problem of how multi-dimensional point clouds or trajectories are topologically structured in phase space. Such a study in phase space is not equivalent to the type of study illustrated in Fig. 10. Working with the topology of real fluid-flow trajectories in physical space requires working in no more than three dimensions, for example. The topological structure we will always be referring to in the present work is defined in phase space, even when studying how such a topological structure is related to the motion of fluid particles in physical space. In 3-D phase space, deterministic flows can be characterized by topological invariants and, therefore, in terms of knots.

Mathematically, a knot is an embedding of a circle in 3-D Euclidean space ℝ³. We can imagine a knot as a thin tangled rope in 3-D space whose ends are glued together (Prasolov and Sossinsky, 1997). Two mathematical knots are equivalent if one can be transformed into the other via a deformation of ℝ³ into itself, known as an ambient isotopy; these transformations correspond to manipulations of a knotted string that do not involve cutting it or passing it through itself. The knot approach – i.e., extracting the knot content of hyperbolic attractors – is based on a geometrical construction that was named template or knot holder. The first way of applying this approach consisted in computing certain knot invariants – such as linking numbers or Conway polynomials – by starting from a set of trajectories (Gilmore, 1998; Gilmore and Lefranc, 2003; Natiello et al., 2007; Letellier and Gilmore, 2013).

There are in fact three steps in this knot-theoretical approach, and the aim of each one is achieved in a particular way:

1. approximate the neighboring unstable periodic orbits (UPOs) around which the flow is evolving with an orbit or closed curve,

2. find a topological representation of the orbit structure and

3. obtain an algebraic description of the topological representation.

The first step is rooted in Henri Poincaré's observation that one can always choose a model's periodic solution as a first approximation of an aperiodic one (Poincaré, 1892). To achieve the first step, one thus applies a close returns method (Mindlin and Gilmore, 1992; Boyd et al., 1994). If the trajectories being studied have been obtained from a data-driven method rather than a model simulation – using, for instance, time-delay embeddings (Takens, 1981) – this step requires long, well-sampled time series that are noise free for orbits to be reconstructed accurately.

Knot theory comes in the procedure's second step and computing the identified knot invariants closes the procedure. Another possibility, instead of using knots, is resorting to braids, as discussed by Natiello et al. (2007). A braid is a collection of strands crossing over or under each other. The braid approach is based on results from Thurston on the classification of two-dimensional diffeomorphisms and on the braid content of a given diffeomorphism (Fathi, 1979). The spirit of the procedure is the same because when connecting the ends of a braid, one ends up with a knot. In the Letellier and Gilmore (2013) Festschrift for Robert Gilmore's 70th birthday, Mario Natiello's Chap. 7 is entitled “A braided view of a knotty story”. The reason is that knots dissolve into trivial objects in dimensions higher than three.

In “How topology came to chaos”, Gilmore (2013b, p. 175) explains that metric and dynamical invariants do not provide a way to distinguish among the different types of chaotic attractors and that a tool of a different nature was needed to create a dictionary of processes and mechanisms underlying a chaotic system. While Gilmore, Lefranc and co-workers were “mulling over implementing a program [based on building tables of linking numbers and relative rotation rates between trajectories, a] better solution became available”. Joan Birman and Robert Williams had shown that the dissipative nature of a flow in phase space allows projecting the points along the direction of the stable manifold by identifying all the points with the same future.

Gilmore continues as follows:

Suppose we have a dissipative [chaotic] flow in three dimensions: There is one positive Lyapunov exponent λ₁>0 [for the unstable direction,] one negative Lyapunov exponent λ₃<0 [for the stable direction], and one zero exponent λ₂=0 “along the direction of the flow”. The dissipative nature of the flow requires ${\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}=\mathrm{0}$. Then it is possible to project points in the phase space “down” along the direction of the stable manifold. This is done by identifying all the points with the same future:

$$x\simeq y\text{ if}\underset{t\to +\mathrm{\infty }}{lim}\left|x\right(t)-y(t\left)\right|=\mathrm{0},$$

[where] x(t) is the future in phase space of the point x=x(0) under the flow. This Birman-Williams identification effectively projects the [3-D] flow down to a two-dimensional set that is a manifold almost everywhere,

except at the points where the flow splits into branches heading towards distinct parts of phase space or at the points where two branches are squeezed together. These mathematical structures were called branched manifolds.

A branched manifold, in the strict sense of the two words that make up the term, can in fact be defined mathematically without reference either to a flow or to the Birman–Williams projection mentioned above. Following Kinsey (1993, p. 64), an n-dimensional manifold is a topological space such that every point has a neighborhood topologically equivalent to an n-dimensional open disk with center x and radius r. Such a manifold is said to be Hausdorff if and only if any two distinct points have disjoint neighborhoods. The second condition is not satisfied precisely at the junction between branches, i.e., at the locations that describe stretching and squeezing of a flow in phase space.

A branched manifold is, therefore, a manifold that is not required to fulfill the Hausdorff property. We prefer this more general definition, instead of the one related to the Birman–Williams projection, for several reasons, including the possibility of extending the concept of a branched manifold to the structure of instantaneous snapshots of random attractors, as we shall see in Sect. 3.4. This mathematical definition of a branched manifold will also let us extend the procedure to cases in which the hypotheses of the Birman–Williams theorem – in which the dynamical system must be hyperbolic, 3-D and dissipative – are not valid. In most geoscientific applications, for instance, uniform hyperbolicity does not apply.

As the topological structure of a branched manifold is closely related to the stretching and squeezing mechanisms that constitute the fingerprint of a certain chaotic attractor, its properties can be used to distinguish among different attractors. This is how one can justify the two-way correspondence between topology and dynamics. This correspondence remains valid in the case of four-dimensional semi-conservative systems (Charó et al., 2019, 2021a), for which the hypotheses of the Birman–Williams theorem do not hold.

The terms “branched manifold” and “template” have often been used interchangeably. We do not regard them as synonyms, for technical reasons that will be important in the development of the concept of templex in Sect. 3.3. A branched manifold is just a particular type of manifold that can be reconstructed from a set of points in ℝⁿ, by approximating subsets of points by disks of local dimension d≤n. Now, to describe branched manifolds immersed in n≥4, we still need a different tool. This tool is in fact provided by homology theory.

3.1 Branched manifold analysis through homologies

Homologies provide an algebraization of topology by building compressed representations of a certain object through cell complexes and by computing essential signatures of the object's shape through homology groups that do not depend on the particular representation used to compute them. Homology groups enable the analysis of n-dimensional manifolds or point clouds, with n as high as desired. This procedure can handle time-delay embeddings produced with shorter and reasonably noisy time series, since the method no longer relies on orbit reconstruction in phase space. In Natiello's terms (Letellier and Gilmore, 2013, Chap. 7), homologies are knotless and orbit-less, and the topological program can be extended to deal with higher-dimensional systems and with real, noisy data.

Other approaches that characterize aspects of dynamical chaos in arbitrary dimensions (e.g., Lefranc, 2006) are somewhat similar to cell complexes. These approaches so far only address estimating the entropy of the flow, which is still an important issue in and of itself.

To illustrate how homologies work, let us take as an example a point cloud obtained by the integration of the deterministic (Lorenz, 1963a) model. Here too the methodology has three steps, but they differ in their tasks and their objectives:

1. approximate the points as lying on a branched manifold,

2. find a topological approximation of the branched manifold and

3. obtain an algebraic description of the topological structure.

Essentially, the passage through the closed orbits is replaced by passing through the branched manifold.

A branched manifold is a generalization of a differentiable manifold that may have singularities of a very restricted type, which correspond to the branching, and it admits a well-defined tangent space at each point. In other words, such a manifold has the property that each point has a neighborhood that is homeomorphic to either a full 2-ball or a half 2-ball, and which is locally homeomorphic to Euclidean space or locally metrizable but not globally so because of the branching (Williams, 1974). A typical branching line is one that joins the “pair of surfaces which appear to merge in the lower portion of Fig. 3”.

As points in our cloud are assumed to lie on a branched manifold, we can classify the points into subsets that constitute a good local approximation of a d disk, where d is the local dimension of the branched manifold and n is the dimension of phase space (d≤n). In the case of the Lorenz attractor, d=2 and n=3. The topological representation is obtained if we convert each subset of points into an individual cell of a cell complex. This complex is sort of a skeleton of the object of interest, namely the Lorenz (1963a) attractor in the case at hand.

Here we use polygons for the cells that pave the attractor's branched manifold. These cells must be correctly glued to each other in order to retain the topological features of the original point cloud. Once the cell complex is constructed, homologies can be computed to yield an algebraic description of the approximating structure. In this review paper, we will not go into the mathematical definitions and theorems required to fully and correctly understand cell complexes and homology theory but only give a taste of the theoretical framework via challenging applications. The reader is referred to Kinsey (1993) for the full mathematics at a comfortable level and to Sciamarella (2019) for a more detailed explanation of the geoscientific applications.

The key point here is that the homology groups represent essential information about the branched manifold, while being independent of the number of cells used to construct the complex (Poincaré, 1895; Siersma, 2012). The topological structure describing the manifold can thus be identified and higher dimensions can be handled, and relatively short and noisy data can be sufficient for this purpose, too.

When Michael Ghil visited the University of Buenos Aires in fall 2018 and got acquainted with this methodology, whose first results were published 2 decades ago (Sciamarella and Mindlin, 1999, 2001), he suggested one should give it a name that identifies and distinguishes it from other methods that had become popular in the meantime in topological data analysis, in particular that of persistent homologies (PHs: Zomorodian and Carlsson, 2004; Edelsbrunner and Harer, 2008). The PH methodology has been enormously successful in problems of shape recognition and classification from large but incomplete datasets.

In dynamic problems, and especially in chaotic dynamics, the PH approach has to contend with the difficulty of finding robust criteria for the degree to which a cell complex represents a manifold that underlies a point cloud (Carlsson and Zomorodian, 2007). Instead of insisting on the improved approximation of such a manifold, PH chooses to display and evaluate the properties of a sequence of cell complexes constructed with a cell creation rule, called a filtration, which depends on a filtration parameter, such as the size of the balls used to approximate the original space around each point of the point cloud. The problem with filtrations is that it is perfectly possible that none of the complexes created by a dynamics-independent rule correctly approximates the branched manifold whose topology is to be described.

For this reason, the Buenos Aires group chose to establish special rules for the construction of a complex, namely rules that take into account that the objective of the reconstruction is not just any arbitrary shape but a branched manifold in phase space. Michael Ghil's suggestion led to the use of Branched Manifold Analysis through Homologies (BraMAH) for this method, a name that says it all and simultaneously recalls the Hindu god of creation and knowledge, which seems very auspicious. The precursors of this technique are four researchers of the Nonlinear Systems Laboratory of the Mathematics Institute at the University of Warwick, who extracted Betti numbers from time series (Muldoon et al., 1993). Betti numbers define the rank of the homology groups, and they can be seen as the number of “holes” in a point cloud. This method served as a guide to construct a cell complex from a point cloud, using singular value decomposition.

We review here briefly the improvements that Sciamarella and Mindlin (1999, 2001) brought to the Warwick approach. The information that was obtained as output by Muldoon et al. (1993) is useful but incomplete if one wishes to identify a branched manifold. As observed in the concluding remarks of the latter paper, the examples used therein involve boundaryless manifolds traversed by a dense orbit, but they suggest potential applications to a wider class of objects including branched manifolds. In order to identify a branched manifold from a point cloud through homologies, it is important to realize that there is much more information contained in a cell complex than just the Betti numbers and that much of this information is relevant to describing the underlying topology.

Sciamarella and Mindlin (1999) were able to show that the branched manifold could be reconstructed with all its features, including torsions and branch locations, from a noisy dataset. The example used was a time series associated with a voice signal of a Spanish speaker articulating the word casa. The topological analysis was carried out on the first vowel, showing that a 3-D time-delay embedding of the acoustic pressure yielded a point cloud with an organization that is typical of a branched manifold. The authors used this dataset to show that the BraMAH method could be applied to reconstruction from a noisy time series, where identifying unstable periodic orbits would have been very difficult or even impossible. They succeeded in characterizing the topology of this dataset but also in showing that their approach and its underlying principles had been fruitful.

In their follow-up paper, Sciamarella and Mindlin (2001) described the algorithm in detail, coded in Wolfram Mathematica, and presented an example of a four-dimensional dynamical system having chaotic solutions of the Shilnikov type. The flow generated by the set of ODEs considered therein was such that any 3-D projection contained self-intersections, stressing the truly four-dimensional nature of the dataset. Sciamarella and Mindlin (1999, 2001) thus showed that their approach could overcome the two main obstacles in the topological analysis of dynamical systems, namely the limitations of dimensionality imposed by the knot-theoretical approach and the noise.

BraMAH can also detect the presence of a Klein bottle in the data, like the one discovered by Mindlin and Solari (1997). Recall that a Klein bottle is a one-sided surface that is formed by passing the narrow end of a tapered tube through the side of the tube and flaring this end out to join the other end. Immersed in three dimensions – as usually shown in the drawings we are used to – a Klein bottle presents self-intersections, and this is why it is a paradigmatic example of a structure that is inherently four-dimensional. In phase space, self-intersections violate uniqueness, and this is why projections may be not only inconvenient but also misleading. Returning to the Muldoon et al. (1993) algorithm, the Betti numbers alone that it computes do not distinguish a Klein bottle from a Möbius strip. The moral is that the topological description of nonlinear dynamical systems in phase space should not only count the holes – as done today by many available topological toolkits – but should be carried out more fully, as in BraMAH. The method is illustrated in Fig. 11, where it is applied to the strange attractor of the deterministic Lorenz (1963a) model, according to Charó et al. (2022) and Charó et al. (2023).

The topological-analysis program has been applied to many fields of science: voice production (Sciamarella and Mindlin, 1999), ocean color (Tufillaro, 2013), biological motor patterns (Mindlin, 2013), financial economics (Gilmore, 2013a), nano-oscillators (Gilmore and Gilmore, 2013) and so on. What is the purpose? To quote Robert Gilmore: “Topological methods can be used to determine whether or not two dynamical systems are equivalent; in particular, they can determine whether a model developed from time-series data is an accurate representation of a physical system. Conversely, it can be used to provide a model for the dynamical mechanisms that generate chaotic data”. The topological program can hence be harnessed for multiple purposes, including but not restricted to

1. validating or refuting models (simulations vs. observations),

2. comparing models (time series generated by different models),

3. comparing datasets (e.g., in situ versus satellite data),

4. characterizing and labeling chaotic behaviors (towards a systematic classification), and

5. classifying sets of time series according to their main dynamical traits (e.g., in Lagrangian flow analysis).

Figure 11 BraMAH analysis of the Lorenz (1963a) attractor. (a) Cell complex with the 0-cells (vertices), 1-cells (line segments) and 2-cells (polygons) constructed as an approximation to subsets of points in a point cloud with N₀=25 000 points; reproduced from Charó et al. (2023) under CC-BY license. (b) A diagram showing the labeled 0-cells. Curved arrows indicate the orientation of the 2-cells $\mathit{\{}{\mathit{\gamma }}_i:i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{14}\mathit{\}}$. The heavy horizontal line in panel (b) indicates the singular line that unites the two branches; reproduced from Charó et al. (2022) with the permission of AIP Publishing.

3.2 Lagrangian coherence in fluid flows

In fluid mechanics, two viewpoints are possible. In the Eulerian viewpoint, fluid motion is observed at specific locations in space, as time passes. In the Lagrangian viewpoint, instead, the observer follows individual fluid particles as they move through the fluid domain. The Eulerian description is more often used for prediction and other purposes. Lagrangian analysis, though, is a powerful way to analyze fluid flows when tracking and understanding the origins and fates of individual particles are important (Bennett, 2006). The fluid envelopes of the Earth system, for instance, exhibit a wide variety of dynamical motions that can act quite differently on mixing and transport. In the ocean, for instance, fluid particles carry tracers such as nutrients, plankton, heat, salt or marine debris (Van Sebille et al., 2018). Hence, in the climate sciences, we are often interested in how particles in the ocean or the atmosphere move and how this motion affects tracer transport.

The oft observed formation of ordered patterns in fluids with complex behavior has led to the search for a theory that could explain Lagrangian coherence in terms of an underlying skeleton responsible for structuring the pathways of sets of fluid particles. These structures may have a finite lifetime, and so one refers to them as finite-time coherent sets (e.g., Williams et al., 2015). Sensitivity to initial conditions makes Lagrangian fluid motion inherently unstable, calling for methods from nonlinear dynamical systems theory (Haller, 2015). In this section, we show how algebraic and chaos topology can help one understand transport in fluid flows (Charó et al., 2020, 2021a) and, more specifically, we demonstrate BraMAH's potential in this setting.

The unsteady or driven double gyre (DDG) system is an analytic model, often used to show how much Lagrangian patterns may differ from patterns in Eulerian fields. Shadden et al. (2005) introduced the DDG model to mimic the motion of two adjacent oceanic gyres enclosed by land, and, since the work of Sulalitha Priyankara et al. (2017), it has been known to present chaotic transport between the two counter-rotating laterally oscillating vortices. The Lagrangian model is defined by the following set of ODEs:

$$\begin{array}{}\text{(26a)}& {\displaystyle }\dot{\mathit{x}}(t;t_{\mathrm{0}},{\mathit{x}}_{\mathrm{0}})=\mathit{v}(t;\mathit{x}(t;t_{\mathrm{0}},{\mathit{x}}_{\mathrm{0}}\left)\right),\text{(26b)}& {\displaystyle }\mathit{x}(t_{\mathrm{0}};t_{\mathrm{0}},{\mathit{x}}_{\mathrm{0}})={\mathit{x}}_{\mathrm{0}}.\end{array}$$

Here the initial conditions x₀ lie in a rectangular domain $\mathrm{\Omega }=[\mathrm{0},\mathrm{2}]\times [\mathrm{0},\mathrm{1}]$ and $\mathit{v}=(u,v)$ is the Eulerian velocity field, which is derived from the streamfunction $\mathit{\psi }=\mathit{\psi }(x,y,t)$ given by

$$\begin{array}{}\text{(27a)}& {\displaystyle }\mathit{\psi }(x_{\mathrm{1}},x_{\mathrm{2}},t)=A\sin \left(\mathit{\pi }f\right(x_{\mathrm{1}},t\left)\right)\sin \left(\mathit{\pi }{x}_{\mathrm{2}}\right),\text{(27b)}& {\displaystyle }\begin{array}{rl}& f(x_{\mathrm{1}},t)=a\left(t\right){x_{\mathrm{1}}}^{\mathrm{2}}+b\left(t\right)x_{\mathrm{1}},a\left(t\right)=\mathit{\eta }\sin \left(\mathit{\omega }t\right),\\ & b\left(t\right)=\mathrm{1}-\mathrm{2}a\left(t\right).\end{array}\end{array}$$

The usual parameter values are A=0.1, η=0.1 and $\mathit{\omega }=\mathit{\pi }/\mathrm{5}$. Note that $u=-\partial \mathit{\psi }/\partial y$ and $v=\partial \mathit{\psi }/\partial x$ and hence the flow is non-divergent at all times ($\partial u/\partial x+\partial v/\partial y=\mathrm{0}$).

Clearly, this DDG model is non-autonomous for η≠0, since the coefficients a,b are periodic in time. Note, however, that the streamfunction $\mathit{\psi }(x_{\mathrm{1}},x_{\mathrm{2}},t)$ given by Eqs. (27a) and (27b) would not correspond to a solution of the Navier–Stokes equations in two dimensions: it is a synthetic example that (i) exhibits somewhat familiar oceanic flow patterns; and (ii) chaotic behavior within certain subsets of the induced particle motion (Shadden et al., 2005). In fact, more realistic Eulerian flows that are solutions of the so-called quasi-geostrophic equations governing the wind-driven oceanic circulation subject to rotation (Pedlosky, 1987; Dijkstra, 2005) are themselves chaotic, rather than periodic in time, for realistic parameter values (Jiang et al., 1995; Dijkstra and Ghil, 2005).

From the Eulerian perspective, the DDG has a time-periodic and simple behavior, a snapshot of which is shown in Fig. 12a. What happens, though, if there is an “oil spill” in the middle of the DDG domain? When injecting a passive tracer, as in Fig. 12b, blank regions appear, i.e., zones of particles in motion that are never reached by the oil spill, and present circular or triangular shapes. The system being conservative, particle behavior depends on the initial particle position being integrated. The oil spill spreads in a chaotic sea surrounding regular islands containing particles where behavior is quasi-periodic. Between the regular islands and the chaotic sea, there are hermetic transport barriers, inhibiting particles to move from one region to another one. This simple, synthetic example demonstrates therewith that flow patterns can effectively differ depending on whether the system is observed in Eulerian or Lagrangian terms. The transport barriers are not even visible in the Eulerian perspective. For further details on particle behavior, the reader is referred to Charó et al. (2019).

Figure 12 Eulerian and Lagrangian perspectives for a fluid flow in the case of the Shadden et al. (2005) driven double gyre (DDG). (a) Vorticity field and (b) passive tracer injected in the middle of the domain (1,0.5). The blank regions are not static: they describe closed orbits within the right and left hemispheres of the domain and deform as they move. The video can be seen at https://youtu.be/W1yndTsvR0g (last access: 27 September 2023). Particles in these four blank regions are trapped within them and exhibit a regular, non-chaotic behavior, while particles in the region visited by the tracer do exhibit Lagrangian chaos. The DDG is known to have an embedded horseshoe near the point (1,1) (Sulalitha Priyankara et al., 2017). Subsets of particles behaving alike – and therefore sharing the same topology – swarm together robustly.

How can BraMAH help us in Lagrangian analysis? The interesting cases, as shown by the DDG example, correspond to dynamical systems that are non-autonomous. But in such systems, some processes involved in the particle dynamics derived from the Eulerian streamfunction are not explicitly described in the two-dimensional space spanned by the particle positions' coordinates. Many authors choose to work in an “extended phase space”, in which time is added as a phase space coordinate.

But such an extended phase space is in fact deceptive, since it assigns a double status to the time variable, which should not play the role of both an independent and a dependent variable. Due to this double status, some tools from autonomous dynamical systems theory do not apply (Charó et al., 2019). The importance of this point in the topology of chaos should not be neglected. In fact, one of the fundamental hypotheses in writing a dynamical system as a set of ODEs is that time is the only independent variable, while all state variables are time dependent.

Working in a space whose dimension is increased by 1 due to introducing the extra ODE $\dot{t}=\mathrm{1}$ leads to certain difficulties in using the tools borrowed from nonlinear dynamical systems theory – for instance, the state space is no longer bounded. In this extended phase space, a periodic orbit is no longer a closed curve, simply because when the system returns to the same state, it does not return to the same point. The very definition of phase space in which a point represents one-to-one a state of the system is no longer valid in the extended phase space.

Many of the properties that are valid in a well-defined phase space are altered in an extended phase space, and topology is one of them. In the case of the DDG model discussed by Charó et al. (2019), the starting point is a non-autonomous system of two ODEs. The extended phase space – with a third ODE written as $\dot{t}=\mathrm{1}$ – is three-dimensional. But the paper shows that a fourth dimension is needed to rewrite the system as an autonomous set of ODEs without using the standard extension trick. The genuine phase space of the autonomously written driven double gyre has four ODEs. Two additional variables are required: u and v.

Such a transformation gets rid of the explicit time dependence with a legitimate procedure that does not run into the previously explained inconsistency. In this four-dimensional phase space, and for certain initial conditions, the topological structure that is obtained is a Klein bottle. A Klein bottle cannot be immersed into a 3-D space without self-intersections: the role of the fourth dimension that is required to rewrite the system in an autonomous form is, therefore, highly relevant here. Thus, to use topological tools self-consistently, one must be prepared to work in a well-defined phase space, and with as many dimensions as required.

In the fluid-flow problem, the four-dimensional phase space complements the Lagrangian variables by an indirect representation of the Eulerian variables. A knotless approach like BraMAH does allow one to work in such a space, which was previously out of reach for a topological analysis. As we shall see, though, in Sect. 3.4, a more general approach to the topological study of NDS and RDS problems is to extend the time-independent BraMAH of Sect. 3.1 and the associated templexes of Sect. 3.3 to the corresponding time-dependent cases.

When applied to time series describing particle trajectories in fluid flows, BraMAH falls within a family of methods that measure the complexity of individual trajectories to identify coherent regions, i.e., regions with qualitatively different dynamical trajectory behavior. Rypina et al. (2011), for instance, use correlation dimension as a measure of complexity. Correlation dimension, though, is a metric invariant, which does not provide information on how to model the system's dynamics. Charó et al. (2020) applied BraMAH to Lagrangian trajectories $\mathit{x}(t;t_{\mathrm{0}},{\mathit{x}}_{\mathrm{0}})$ and obtained the topology of the associated branched manifold in the full four-dimensional phase space of the DDG equations in their Lagrangian form (26a) and (26b). This result is achieved by deriving the recipes that knead the DDG model's dynamical behavior in phase space, without having to look into the geometrical complexity of individual particle trajectories.

Returning now to the oil spill in the middle of the DDG system's domain Ω, Charó et al. (2020) applied BraMAH to 8528 fluid particles in a four-dimensional reconstructed phase space. Only five distinct topological classes emerge, and their characteristic cell complexes are plotted in Fig. 13.

Figure 13 The analysis of a set of 8528 particles advected by the DDG flow field yields five topological classes. These five classes are obtained by applying BraMAH to four-dimensional point clouds; the plots in the figure are three-dimensional projections of representative cell complexes for each of the five classes. Four of them involve quasi-periodic particle motion, and only one of them, which is represented by the third cell complex, points to a branched manifold that refers to the so-called chaotic sea (colored in blue in Fig. 14 below). The 1-cells, i.e., the lines that are highlighted in color, indicate the generators of the homology group H₁ of holes in each cell complex. From Charó et al. (2021a) with permission from Cambridge University Press.

From left to right, Class I corresponds to a strip, Class II to a torus and Class III to a branched manifold with three 1-holes – i.e., with a Betti number β₁=3 – and a torsion that is indicated by the orientability chain. The remaining complexes are of Class IV, with the topology of a Klein bottle, and of Class V, which is a very peculiar kind of torus that involves a torsion and a weak boundary. Each topological class is assigned a color (class I: green; class II: magenta; class III: blue; class IV: red; class V: orange) used in Fig. 14 to tag the particles in motion and thus identify distinct particle sets that stay coherent while moving and being distorted. The frontiers between differently colored regions will be called separators. Such flow separators are associated with LCSs that are known to separate dynamically distinct regions in fluid flows (Kelley et al., 2013).

Figure 14 Coloring of 8528 particles in motion in a DDG field, with colors corresponding to the topological structure of the particle trajectories in phase space. The boundaries between distinct colors are fairly well defined, displaying the existence of transport barriers that separate non-mixing regions, like the green, orange, red or magenta, vs. the chaotic sea (blue). A direct correspondence is found between the regions identified by the topological BraMAH analysis and those observed dynamically using a Poincaré section, as in Charó et al. (2019, Fig. 5), or a finite-time Lyapunov exponent study, as in You and Leung (2014, Fig. 12). From Charó et al. (2021a) with permission from Cambridge University Press.

The presence of the Klein bottle as Class IV among the five classes in Fig. 13 stresses the importance of being able to work in a sufficiently high-dimensional phase space that guarantees an autonomous setting: as mentioned in Sect. 3.1 before, the Klein bottle cannot be immersed in three dimensions without self-intersections.

Charó et al. (2021a) further emphasized that BraMAH can identify and describe LCSs in a fluid flow from a sparse set of particles and achieve this without inspecting relative particle positions. The method differs from previous ones because it describes transport by how particles behave without looking at where they go. Such a dynamical analysis ends up pointing to finite-time coherent sets, thanks to the property that particles sharing equivalent dynamics tend to stay together. The same authors have also successfully used BraMAH to study numerically generated fluid particle behavior in the wake behind a rotary oscillating cylinder (Charó et al., 2021a).

The BraMAH applications reviewed in this subsection demonstrate substantial progress in Lagrangian analysis, by providing a method that enables one to identify coherent sets without previous knowledge of the flow field. This particular set of results also shows methodological progress in chaos topology, since it appears that BraMAH can help describe the topological structure of non-dissipative, Hamiltonian systems. Recall, as a stepping stone in this direction, the analogy between the non-divergence of a fluid flow in physical space, like the DDG model, and the Hamiltonian character of a dynamical system's flow conserving volume in phase space, like the equations of celestial mechanics (Poincaré, 1892; Arnold et al., 2007).

3.3 Templexes for dynamical systems

Structures in phase space are special because they are not just spatial objects: they are associated with a semi-flow on them, which is sometimes represented by arrows. A cell complex can effectively encapsulate the properties of a branched manifold in standard space, but it will not convey the fact that, when the cells in a complex represent a semi-flow on a spatial object, they can be traversed in an arbitrary order only at the expense of forgetting about the semi-flow. In other words, time is absent from the description. Including the arrow of time in the description calls for a more refined mathematical object, in which the topological properties of a flow in phase space come to light through the combined analysis of both the spatial structure of the underlying branched manifold and of the semi-flow upon it.

Figure 15 Solution trajectories for (a) the spiral-type (a=0.343295, b=2, c=4) and (b) the funnel-type (a=0.492, b=2, c=4) attractors of the Rössler (1976) model.

Charó et al. (2022) introduced such a novel type of mathematical object and called it a templex, a word obtained from the contraction between “template” and “complex”. A template in dynamical systems theory is a synonym for a knot holder (Birman and Williams, 1983a; Tufillaro et al., 1992; Ghrist et al., 1997). Since Mindlin and Gilmore (1992), templates have been used to describe three-dimensional flows from experimental data in many fields: to study a three-species food chain model in ecology (Letellier and Aziz-Alaoui, 2002), to forecast the time series of sunspot numbers (Aguirre et al., 2008), or to better understand delayed interactions between cancer cells and the micro-environment (Ghosh et al., 2017). Albeit limited to three dimensions, a template provides a description of an attractor at a level of detail that homologies alone cannot achieve.

The Rössler (1976) model,

$$\begin{array}{}\text{(28a)}& {\displaystyle }\dot{x}=-y-z,\text{(28b)}& {\displaystyle }\dot{y}=x+ay,\text{(28c)}& {\displaystyle }\dot{z}=b+z(x-c),\end{array}$$

provides a simple example. Changing two of the parameter values in the governing Eqs. (28b) and (28c), one can produce two distinct chaotic attractors, shown in Fig. 15: in panel a, the spiral case, with $(a=\mathrm{0.343295},b=\mathrm{2},c=\mathrm{4})$, and in panel b, the funnel case, with $(a=\mathrm{0.492},b=\mathrm{2},c-\mathrm{4})$. As discussed by Charó et al. (2022), these two structures can be approximated by cell complexes that are homologically equivalent. But templates are able to discriminate between the two cases using the concept of strip.

For strongly dissipative systems, like the Rössler attractor, the number of monotone branches of the first-return map provides the number of strips required to construct the corresponding template. Strips are cylinders, in topological terms, but one must beware that the meaning of strip in a template is not introduced to refer to a topological class but to discriminate between the different paths followed by the flow along the branched manifold. A strip is typically defined between a splitting chart and a joining chart, in which the strips are split and joined, respectively. Thus, in the template terminology, the spiral attractor has two strips, while the funnel attractor has three strips, as shown in Fig. 16a and b, respectively.

Figure 16 Templates for the (a) spiral-type and (b) funnel-type attractors for the Rössler (1976) model. Reproduced from Letellier et al. (1995) with the permission of AIP Publishing.

Strips in a template are associated with a tearing of the flow. They are sometimes split in a fictitious manner, introducing false holes into the branched manifold, even if these strips are not necessarily delimited by boundaries or associated with holes in the sense of homologies. Their number can be obtained, for strongly dissipative systems, by computing the number of monotone branches of the first return map. But where are these strips in a cell complex? As mentioned above, they cannot be directly identified with holes in the latter. Can they be identified all the same from some other properties of the cell complex? The short answer is yes but not without the information that is contained in the flow on the cell complex rather than just in the cell complex itself.

The templex thus combines all the essential information that is relevant to the topology of the branched manifold and to the flow on it. The flow on the cell complex is represented by a directed graph (digraph) (e.g., Bang-Jensen and Gutin, 2008), whose nodes are the highest-dimensional cells and whose edges, or arcs, are provided by the cell connections that are consistent with the flow. In a templex, the cell complex and the digraph are interrelated. Computations carried out on the two complementary objects yield a description of the branched manifold and of the permitted nonequivalent paths around it.

Algebraic computations on a templex provide, on the one hand, the already known properties of the cell complex – such as the homology groups, torsion groups and weak boundaries – that describe the branched manifold; on the other hand, they provide the properties of the flow on this structure. The topology of a templex is described in terms of a set of sub-templexes that will be called stripexes, since they play the same role as strips in a template. This is no longer done at the price of introducing false holes or boundaries to separate the strips. It is achieved through a set of well-defined operations that include flow-orienting the cell complex; minimizing the cell structure at the joining loci, where the tearing of the flow takes place, to obtain a generating templex; calculating the cycles of the digraph; and checking for local twists, since uneven torsions in a strip correspond to a local twist in a stripex. The reader is referred to the steps in Charó et al. (2022) for further details. This dissection of the cell complex into stripexes provides the information that enables one to distinguish the topological properties of the two Rössler attractors from each other. In order to see how, consider Fig. 17, which illustrates the templexes for the two types of Rössler (1976) attractor.

Figure 17 Templexes for (a) the spiral and for (b) the funnel attractors corresponding to the structures in phase space shown in Fig. 15a and b, respectively. Cell complexes (above) are shown as planar diagrams, with the convention that points (0-cells) and segments (1-cells) with identical labels must be glued to each other. The digraphs provide the allowed connections between the polygons (2-cells) labeled γ_i for the spiral case and ${\mathit{\gamma }}_i^{\prime }$ for the funnel case, with i∈ℕ. The two cell complexes are homologically equivalent. Reproduced from Charó et al. (2022) with the permission of AIP Publishing.

The cell complex of a templex can be seen as a dynamic kirigami or cutout paper model, made of pieces that fit together; in this case, the pieces are polygons. Note that points or segments with the same label must be glued together when constructing the paper model. The digraph can be seen as a map of the flow-compatible connections between the pieces. Combining the cell complex and the digraph, we can define and algebraically compute the stripexes. For details on this procedure, the reader is again referred to Charó et al. (2022). The stripexes for the spiral attractor are given by two paths along the cell complex, indicated by the two cycles below, the first of which is twisted.

$$\begin{array}{}\text{(29a)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}\to {\mathit{\gamma }}_{\mathrm{2}}\to {\mathit{\gamma }}_{\mathrm{4}}\to {\mathit{\gamma }}_{\mathrm{6}}\to {\mathit{\gamma }}_{\mathrm{1}},\text{(29b)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}\to {\mathit{\gamma }}_{\mathrm{2}}\to {\mathit{\gamma }}_{\mathrm{3}}\to {\mathit{\gamma }}_{\mathrm{5}}\to {\mathit{\gamma }}_{\mathrm{1}}.\end{array}$$

There are three stripexes for the funnel attractor, as shown below, and only the middle one presents a local twist:

$$\begin{array}{}\text{(30a)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{2}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{3}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{6}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{1}}^{\prime },\text{(30b)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{2}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{3}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{5}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{1}}^{\prime },\text{(30c)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{2}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{3}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{6}}^{\prime }\to {\mathit{\gamma }}_{\mathrm{1}}^{\prime }.\end{array}$$

The description in terms of stripexes provided by the two templexes in Eqs. (29a), (29b), and (30a)–(30c) is equivalent to the strips in the templates of the spiral and the funnel case of the Rössler (1976) attractor, as shown in Fig. 16. Let us recall that templates are knot holders and can therefore only be obtained for three-dimensional flows, while templexes can be computed for four- or higher-dimensional dynamical systems, as shown in Charó et al. (2022, Sect. IV).

3.4 Algebraic topology and noise-driven chaos

BraMAH and the associated templexes, as presented so far, provide a topological description that holds within an autonomous and deterministic framework. As discussed in Sect. 2.2 regarding dynamical systems theory for the climate sciences, the question that naturally arises, though, is whether we can take one step beyond, namely extend the topological perspective to NDSs and RDSs, which provide the appropriate mathematical framework to tackle the effects of time-dependent forcing on intrinsic climate variability (Ghil, 2019; Ghil and Lucarini, 2020; Tél et al., 2020). Of the two forms of time-dependent forcing, it is the random one that is more challenging. Moreover, the topological characterization of noise-driven chaos is crucial in the understanding of complex systems in general, where part of the dynamics remains unresolved and is modeled as noise.

An example involving not only deterministic time dependence but also random forcing was presented in Eq. (25a)–(25c) and Fig. 8 of Sect. 2.2. In the stochastically perturbed (Lorenz, 1963a) model's random attractor, termed LORA (Chekroun et al., 2011), the stretching and folding mechanisms shape the flow in phase space yielding a time-evolving branched manifold, which must be analyzed accordingly. Nothing prevents one from applying BraMAH to successive point clouds, each of which corresponds to a single snapshot, and comparing the topological properties of these instantaneous cell complexes, as done for the first time by Charó et al. (2021b).

Such an analysis was performed by Charó et al. (2021b) for a fixed realization of the driving noise dη at different instants in time. In order to construct the cell complexes, these authors first sieved the LORA point clouds to retain the most populated regions in phase space. The deterministic concept of branched manifold (Williams, 1974) was extended to the stochastic framework by redefining it locally as an integer-dimensional set in phase space that robustly supports the point cloud associated with the system's invariant measure at each time instant. The numerical results show that BraMAH captures LORA's time-evolving homologies (Charó et al., 2021b), as shown here in Fig. 18. The topologies differ from the deterministic Lorenz model's strange attractor, and the noise-driven model's branched manifold exhibits sharp topological changes in time.

Figure 18 Three LORA snapshots with the noise variance σ=0.3 and cloud size N₀=10⁸. Sieved point clouds (a–c) and cell complexes (d–f): (a, d) t=40.09, (b, e) t=40.18 and (c, f) t=40.27. The cell complexes are not homologically equivalent from one snapshot to another: their Betti numbers are ${\mathit{\beta }}_{\mathrm{1}}=\mathrm{3},\mathrm{10}$ and 4 for $t=\mathrm{40.09},\mathrm{40.18}$ and 40.27, while the Betti number for the deterministic strange attractor in Fig. 11 is β₁=2, stemming from the two holes around the two convective fixed points on either “wing” of the butterfly. Reproduced from Charó et al. (2021b) with the permission of AIP Publishing.

The stochastic branched manifold, characterized by a single-cell complex for each snapshot, does not contain any information about the future or the past of the invariant measure. The flow in a cell complex representing the invariant measure on a random attractor can no longer be represented within that cell complex, as done when using a deterministic templex, like the one described in Sect. 3.3 for the Rössler (1976) model. Incorporating time into this formalism requires establishing a link between the cell complexes of distinct snapshots.

But how can one track changes between different cell complexes without using specific individual cells? Let us recall that the number of cells and their distribution in a cell complex are arbitrary and that homology groups are conceived so as to cancel out the extraneous information in the cells and to only retain the essential properties of the topological space. Homologies will thus provide the key to connect a cell complex of a random attractor at a given instant to a cell complex corresponding to another instant. For a random attractor, we will endow a set of cell complexes with a digraph that does not connect cells within a single complex, as in Fig. 17, but holes of cell complexes at distinct instants of time. This is the key idea that led Charó et al. (2023) to construct their random templexes.

Tracking holes requires some caveats, though. Homology groups and the associated Betti numbers are independent of the particular set of cells forming a cell complex. Hence, the holes or generators of a homology group can be expressed in terms of one of several representative cycles that need not strictly follow the boundary of the holes, as shown in Fig. 19. A representative cycle may wander around a hole, without tightly encircling the empty space. Still, the boundaries of the holes can be retrieved algebraically, from the cell complex itself, as shown by Charó et al. (2023). We can thus define a random templex as an indexed family of BraMAH cell complexes hanging together by a digraph. In this digraph, each node is a minimal hole of a given cell complex and the edges denote the connections between minimal holes occurring at successive time instants.

Figure 19 Cell complex of the deterministic Lorenz (1963a) model's attractor, as shown in Fig. 11a here. Emphasized in color in this figure are (a) the holes obtained in the homology computation and (b) the tight or minimal holes. From Charó et al. (2023) under CC-BY license.

What does the random templex, thus defined, encode? In the life of a random attractor, there may be time intervals within which the branched manifold evolves geometrically but maintains its homological properties. Topology can be said to change when the holes that are being tracked from one snapshot to the next are created or destroyed. Some of them can be found to split or merge. Such changes are associated with what we call hereafter a topological tipping point (TTP) (Charó et al., 2021b, 2023). Since the Betti numbers are integers, any changes in them must be sudden. In fact, these sudden changes could already be noticed visually in the LORA video published by Chekroun et al. (2011) at https://vimeo.com/240039610 (last access: 27 September 2023).

To confirm this further, Charó et al. (2021b, Fig. 4) showed that the time intervals over which the Betti numbers changed drastically were quite short, i.e., no longer than δt=0.09, as reproduced in Fig. 18 herein. This time interval is very short indeed, compared to the characteristic time to switch wings for a trajectory of the deterministic Lorenz (1963a) model, which is of the order of units.

Charó et al. (2023, Fig. 5) further showed that the numerically observed intervals over which the set of minimal 1-holes change can be even shorter, with δt≤0.065. More interestingly, these authors demonstrated that TTPs can be identified and classified using the digraph of a random templex.

Figure 20 here shows the “story” of two holes in a finite time window $T_{\mathrm{w}}=\mathrm{40.065}\le t\le \mathrm{40.110}$ of LORA's life in the form of two tree plots; the two holes, 73 and 74, lie on opposite wings of the LORA butterfly at the window's initial time. For the sake of simplicity, we kept only two of the 15 connected components of the complete finite-time random templex of LORA for T_w; see Charó et al. (2023, Fig. 6) for the complete picture. Square nodes correspond either to an initial or to a final node for a given time window. A splitting TTP occurs where two or more edges emerge and a merging node receives two or more edges. Similarly, there is a creation or annihilation TTP where an initial or a terminal node in a connected component of the digraph does not correspond to the boundaries of the time window: square nodes cannot be TTPs since the preceding or following instant in time is outside the inspected time window.

Figure 20 A “day in the life” of two mutually symmetric holes of LORA for a fixed noise realization and noise intensity σ=0.3. The time window is $T_{\mathrm{w}}=[\mathrm{40.065},\mathrm{40.11}]$. The nodes are highlighted in different colors according to the type of event: creation in green, destruction in black, splitting in red and merging in blue.

The indices in Fig. 20 label a hole at a certain instant. Tracking enables one to connect, for instance, hole 73 with 91, which will split into holes 108 and 109; this is why hole 91 is colored in red. A symmetric splitting event can be found on the other wing of the animated butterfly, where hole 74 becomes hole 96, which splits into holes 112 and 116. All these holes can be located in phase space using the coordinates of each hole's barycenter. Plotting the position of the barycenters of all the holes present in the analysis in phase space, we obtain a constellation set, as shown in Fig. 21. Each constellation contains the immersed nodes and edges forming a connected component in the digraph and transforms the tree plots into actual paths in phase space. In other words, embedding the digraph of the random templex into phase space, one can represent parsimoniously the evolution of LORA's topology over a given time interval. Such a representation might provide access to a more detailed description of the flow dynamics in a random attractor.

Figure 21 A particular constellation out of the set that represents the essence of the evolution of LORA's 1-holes within a time window $T_{\mathrm{w}}=[\mathrm{40.065},\mathrm{40.11}]$. The plot shows the embedding into the Lorenz (1963a) model's phase space of a connected component in the digraph of LORA's random templex by using the coordinates of the barycenters of the nodes. Regular nodes in a constellation are marked by open stars, while nodes associated with TTPs – such as splitting or merging of holes – are marked by filled stars. From Charó et al. (2023) under CC-BY license.

4 Concluding remarks

The purpose of this paper was to provide an account of the convergence between two strains of Henri Poincaré's heritage – dynamical systems theory (Poincaré, 1892, 2017) and algebraic topology (Poincaré, 1895; Siersma, 2012) – and their joint applications to the climate sciences.

4.1 Summary

In Sect. 1, we provided a bird's eye view of the evolution of these two strains of research since the mid-20th century and how they started to be applied to issues related to fluid flows at both engineering and planetary scales. Sections 2 and 3 developed next in greater detail (a) the concepts and methods associated with dynamical systems and their applications to the climate sciences and (b) those associated with algebraic topology and their applications first to engineering fluid dynamics and then to the climate sciences. Note that the pioneering references mentioned in Sects. 1.1 and 2 date back to the early 1960s, while those of Sects. 1.2 and 3 start in the early 1980s. It is clear that, on the whole, algebraic topology started playing a noticeable role in the climate sciences about 2 decades later than dynamical systems theory.

Section 2.1 covered autonomous dynamical systems, in which neither the forcing nor the coefficients depend explicitly on time. A very extensive and thorough mathematical theory exists and certain aspects of it are well known to a substantial fraction of climate scientists; see, for instance, Ghil and Childress (1987) and Dijkstra (2013). The contents of this section emphasized elementary bifurcations – saddle–node and fold, pitchfork, and Hopf, summarized in Fig. 4 – ending with bifurcation trees and routes to deterministic chaos.

The material in Sect. 2.2 refers to systems with explicit time dependence in the forcing or the coefficients, and it is much newer. The theory of NDSs and RDSs only started in the 1960s – with George Sell, followed by Ed Ott and colleagues and by Ludwig Arnold, Hans Crauel and Franco Einaudi (Sell, 1971; Romeiras et al., 1990; Crauel and Flandoli, 1994; Arnold, 1998; Caraballo and Han, 2017; Kloeden and Yang, 2020) – and its applications to the climate sciences started merely 15 years ago (Ghil et al., 2008; Chekroun et al., 2011; Tél et al., 2020).

We first explained in this section the essential difference between forward and pullback attraction, i.e., between convergence in time of single-parameter and two-parameter semigroups of solutions to the governing equations. Simple examples of pullback attractors (PBAs) were given to familiarize newcomers with the appropriate concepts and methods; see again Figs. 5 and 6. The sequence of examples was concluded with the striking random attractor of the stochastically perturbed Lorenz model, as introduced and studied by Chekroun et al. (2011); see Figs. 7 and 8. Finally, tipping points were introduced as the proper generalization to NDSs and RDSs of the elementary bifurcations for autonomous systems described in Sect. 2.1 (Fig. 9).

In Sect. 3, we presented topological methods in a dynamical systems perspective. We reviewed the advantages of working with homology theory in order to overcome the limitation of a three-dimensional space imposed by using knot theory, since knots simply disentangle in higher dimensions. In Sect. 3.1, we showed that homologies provide a knotless method and that a BraMAH cell complex can be used to describe the spatial structure of a flow in phase space by using homology group generators, weak boundaries and torsion groups (Fig. 11). We described an application of these concepts and methods to Lagrangian analysis in Sect. 3.2, by showing how to define and detect localized coherent sets (LCSs) for fluid flows in physical space; see again Figs. 12–14.

In Sect. 3.3, we dealt with the fact that BraMAH alone does not provide a robust skeleton of the flow in phase space on its branched manifold, even for an autonomous, deterministic system. To obtain such a robust and parsimonious flow description in phase space, we introduced a directed graph (digraph), whose nodes are the cells, while the edges point from one cell to another, in a way that is consistent with the flow on the branched manifold. The mathematical object that combines such a digraph with the underlying cell complex is called a templex. Homologically equivalent attractors – such as the spiral and funnel versions of the Rössler (1976) attractor – can be distinguished using a templex, given its digraph's properties; see Figs. 15–17.

Finally, in Sect. 3.4, we discussed how a digraph, and hence a templex, can be generalized from the autonomous and deterministic version of Sect. 3.3 to non-autonomous and random dynamical systems; see Figs. 18–20. To define a random templex, one needs to shift the perspective from defining a digraph on the single-cell complex of an autonomous system to an indexed family of cell complexes at successive instants in time and the vertices pointing from one cell at time t=t_j to the corresponding one at time $t=t_{j+\mathrm{1}}$.

The fact that the change in the set of minimal holes of a cell complex at t=t_j to the next one is sudden allowed us to rigorously define topological tipping points (TTPs) as happening at an instant at which such a sudden change occurs. It is these TTPs that are a matter of particular interest for future work in the overlap of the two fields that we considered in this review, namely dynamical systems and algebraic topology.

4.2 Perspectives

As usual, when stumbling upon some striking findings, there are two kinds of paths that one might wish to pursue: (i) more general or stronger theoretical results and (ii) interesting applications. Clearly we have some rather striking findings, and we will outline some intriguing paths to pursue, of both kinds, as well as connections between the two kinds of paths.

TTPs in the templex of an NDS or RDS are obviously connected with a lot more detailed information in phase space about the system under investigation than one might suspect from the usual kinds of bifurcation-induced, noise-induced and rate-induced transitions – or B tipping, N tipping and R tipping – discussed by Ashwin et al. (2012) and mentioned here toward the end of Sect. 2.2. But what does that say about changes in the flow in physical space? Could this localization in phase space say something about the association with localized sudden changes in the flow in physical space, i.e., with the “tipping elements” of Lenton et al. (2008)? The use of systematically derived reduced-order models (Kondrashov et al., 2015, 2018; Gutiérrez et al., 2021), for which both TTPs and the better understood dynamical tipping points can be computed fairly easily, could help clarify such relationships and the associated precursors of critical transitions.

An interesting example, among many, of localized changes in Earth's physical space is that of persistent anomalies (Dole and Gordon, 1983) or flow regimes (Legras and Ghil, 1985; Ghil and Childress, 1987, Chap. 6) or weather regimes (Hannachi et al., 2017). Strommen et al. (2023) have recently applied a multiparameter PH method (Carlsson and Zomorodian, 2007; Vipond et al., 2021) to decide more objectively the much debated existence of distinct regimes in the large-scale atmosphere's phase space (Ghil and Robertson, 2002; Hannachi et al., 2017; Robertson and Vitart, 2018). Their findings certainly strengthen the affirmative reply to the quandary. Before proceeding to the next quandary, though, let us consider briefly the issues that are still open in applying this approach to regime identification.

In applying the multiparameter PH method to the classical (Lorenz, 1963a) convection model, Strommen et al. (2023) essentially equate the existence of distinct regimes to the existence of two holes in its branched manifold. Doing so, however, is not quite enough. To explain why, we revisit the example of the Rössler attractor discussed in our Sect. 3.3. This attractor's spiral form, shown in Fig. 15a, has a single hole, but it has two strips: one strip related to the system's slow branch and the other strip to its fast branch. These two branches, however, are not separated by a hole in the sense of homology groups, and this is why homologies alone cannot distinguish between them. Instead, the templex introduced in the same subsection captures these two ways of circulating around the attractor in terms of two stripexes, despite the fact that the branched manifold is single-holed.

The existence of multiple regimes in a dynamical system is certainly associated with its attractor's nontrivial topological structure, as Strommen et al. (2023) state, but this nontrivial topology is not necessarily captured by homologies alone. As explained throughout this work, the templex – with its stripexes and digraph – contributes additional tools to accurately describe the phase-space topology of a flow and of its single or multiple regimes.

Ghil and Robertson (2002), however, asked a more subtle question: is the so-called low-frequency variability (LFV) of the atmosphere – which is closely related to the rapidly growing interest in subseasonal-to-seasonal (S2S) predictability (e.g., Robertson and Vitart, 2018) – oscillatory, i.e., wavelike, or episodic and intermittent, i.e., particle-like. These authors, with an obvious nod to the classical problem of quantum mechanics, formulated the question as “waves” vs. “particles”. Two decades later, this question is still far from settled, as discussed quite recently by Ghil et al. (2018) and by Ghil and Lucarini (2020).

Figure 22 Schematic overview of atmospheric LFV mechanisms. From Ghil et al. (2018) with the permission of Elsevier.

Which type of phenomena dominate atmospheric LFV? There are two apparently contradictory descriptions: oscillatory, wavelike flow features or geographically fixed, particle-like, episodic flow features; e.g., blocking of the westerlies (particle-like) or intraseasonal oscillations (wavelike), with periodicities of 40–50 d (Ghil et al., 2018). In fact, these two are by now accompanied by several more key dynamical mechanisms of midlatitude LFV variability, summarized in Fig. 22.

The simplest approach to persistent anomalies in midlatitude atmospheric flows on 10–100 d timescales is to regard them as due to the slowing down of Rossby waves or to their linear interference (Lindzen et al., 1982; Lindzen, 1986). This approach is illustrated in the sketch labeled c within the figure: zonal flow Z and blocked flow B are simply slow phases of a harmonic oscillation, like the neighborhood of $t=\mathit{\pi }/\mathrm{2}$ or $t=\mathrm{3}\mathit{\pi }/\mathrm{2}$ for a sine wave sin (t); or else they are due to an interference like that occurring for a sum Asin (t)+Bsin (3t) near $t=(\mathrm{2}k+\mathrm{1})\mathit{\pi }/\mathrm{2}$. A more versatile, quasi-linear version of this approach is to study long-lived resonant wave triads between a topographic Rossby wave and two free Rossby waves (Egger, 1978; Trevisan and Buzzi, 1980; Ghil and Childress, 1987, Sect. 6.2). Neither version of this approach, though, explains the anomalies' organizing into distinct flow regimes.

Rossby et al. (1939) initiated a different, genuinely nonlinear approach by raising the possibility of multiple equilibria as an explanation of preferred atmospheric flow patterns. These authors drew an analogy between such equilibria and hydraulic jumps and formulated simple models in which similar transitions between faster and slower atmospheric flows could occur. This multiple-equilibrium approach was then pursued quite aggressively in the 1980s (Charney and DeVore, 1979; Charney et al., 1981; Legras and Ghil, 1985; Ghil and Childress, 1987, Sect. 6.3–6.6), and it is illustrated in Fig. 22 by the sketch labeled a: one version of the sketch illustrates models that concentrated on the B–Z dichotomy (Charney and DeVore, 1979; Charney et al., 1981; Benzi et al., 1986) and the other on models (e.g., Legras and Ghil, 1985) that allowed for the presence of additional clusters, like those found by Kimoto and Ghil (1993a) or Smyth et al. (1999), viz. opposite phases of the North Atlantic Oscillation (NAO) and the Pacific North American (PNA) anomalies – dubbed RNA for Reverse PNA and BNAO for Blocked NAO in sketch a of Fig. 22. The LFV dynamics in this approach are given by the preferred transition paths of a Markov chain between two or more regimes.

A third approach is associated with the idea of oscillatory instabilities of one or more of the multiple fixed points that can play the role of regime centroids. Thus, Legras and Ghil (1985) found a 40 d oscillation arising by Hopf bifurcation off their blocked regime B, as illustrated in sketch b of the figure. An ambiguity arises, though, between this point of view and a complementary possibility, namely that the regimes are just slow phases of such an oscillation, caused itself by the interaction of the midlatitude jet with topography. Thus, Kimoto and Ghil (1993b) found, in their observational data, closed paths within a Markov chain whose states resemble well-known phases of an intraseasonal oscillation. Kondrashov et al. (2004) confirmed the likelihood of such a scenario in the intermediate-complexity model of Marshall and Molteni (1993). Furthermore, multiple regimes and intraseasonal oscillations can coexist in a two-layer model on the sphere within the scenario of “chaotic itinerancy” (Itoh and Kimoto, 1996, 1997).

Finally, sketch d in the figure refers to the role of stochastic processes in LFV variability and S2S prediction, whether it be noise that is white in time, as in Hasselmann (1976) or in linear inverse models (Penland, 1989, 1996; Penland and Ghil, 1993; Penland and Sardeshmukh, 1995), or red in time, as in empirical model reduction and multilayer stochastic models (Kravtsov et al., 2005, 2009; Kondrashov et al., 2013, 2015; Gutiérrez et al., 2021), or even non-Gaussian (Sardeshmukh and Penland, 2015). Stochastic processes may enter into models situated on various rungs of the modeling hierarchy, from the simplest conceptual models to high-resolution global climate models. In the latter, they may enter via stochastic parametrizations of subgrid-scale processes (e.g., Palmer and Williams, 2009, and references therein), while in the former they may enter via stochastic forcing, whether additive or multiplicative, Gaussian or not (e.g., Kondrashov et al., 2015; Gutiérrez et al., 2021, and references therein). Dorrington and Palmer (2023) recently drew attention to yet another mechanism of interaction between stochastic forcing and nonlinear regime dynamics that might modify the picture.

How might topological data analysis contribute to clarify this thicket of apparently contradictory descriptions of LFV? One hint is found in the work of Lucarini and Gritsun (2020), who showed that blocking can be studied by extracting from the complex high-dimensional dynamics of a model its essential building blocks, given by truly nonlinear modes. In this work, they abandoned the classic identification of weather regimes with fixed points, as in Charney and DeVore (1979), and directly considered the chaotic nature of the atmosphere, using the unstable periodic orbits (UPOs) that are a key component of the Gilmore (1998) topological analysis of the chaos program.

This UPO-based approach did confirm certain theoretical results of Legras and Ghil (1985) and the laboratory findings of Weeks et al. (1997) – about the relative stability and persistence of blocked and zonal flows – as well as providing further insights into the waves-versus-particles quandary (Ghil et al., 2018). UPOs can be very useful in characterizing a chaotic system, since the information about them can be obtained in a finite time, which is particularly useful in nonstationary systems, and because a single UPO can already provide substantial information (Amon and Lefranc, 2004). Still, Lucarini and Gritsun (2020) found it quite hard to carry out the necessary computations of very numerous UPOs even for the relatively simple Marshall and Molteni (1993) model.

As explained here in Sects. 1.2 and 3.1, BraMAH is crucially inspired by the Gilmore (1998) program. Yet it is more powerful than the knots-and-braids methodology, which is limited by the dimensionality of the phase spaces that it can be applied to. Likewise, it is more computationally efficient than the UPO methodology, and, as shown at the beginning of this subsection, it provides considerably more information than the PH methodology for chaotic dynamics.

It is thus conceivable, although it remains to be demonstrated, that the additional tools brought to the table by the mathematical object we called templex – namely the digraph and stripexes – could help explore, in a highly simplified setting, issues like the existence and multiplicity of regimes, as well as of the presence of oscillatory features in the dynamics. As explained in the Rössler (1976) attractor context, stripexes can greatly help, beyond counting holes, to determine regime multiplicity.

Finally, as stated at the end of Sect. 3.4, minimal or tight holes in the cell complex of a random templex can be located in phase space using the coordinates of each hole's barycenter. Plotting the position of the barycenters of all the holes present in the analysis in phase space yields a constellation set as shown in Fig. 21. The topology of this constellation set deserves further exploration. It might lead, quite conceivably, to the generalization of a stripex for a random templex and therefore to the extraction of the nonequivalent paths that a nonlinear system follows when driven by multiplicative noise. Random stripexes should provide us with the stretching, squeezing, folding and tearing mechanisms that knead, mold and alter the topological structure of a noise-driven flow in phase space.

The extension of the templex from autonomous and deterministic systems (Charó et al., 2022) to non-autonomous and stochastic ones (Charó et al., 2023) opens the way to the exploration of key aspects of the LFV quandaries associated with Fig. 22. More broadly, it can facilitate exploring a plethora of climate problems that are strongly affected by time-dependent forcing, such as anthropogenic greenhouse gas and aerosol emissions, and stochastic components, such as cloud microprocesses. One can imagine, for instance, applying methods from network theory (Bang-Jensen and Gutin, 2008; Coluzzi et al., 2011; Colon and Ghil, 2017) to investigate the presence of cyclicity in a given model's or dataset's digraph as well as issues of multimodality or multistability.

More broadly, complex networks (Zou et al., 2019) have found numerous applications in the climate sciences in recent years and could provide other links between topology and the multivariate time series analysis of nonlinear phenomena. The field of complex networks shares many of the challenges that are faced by the topology of chaos. Algebraic topology is not mentioned in the Zou et al. (2019) paper, but there have been some papers applying PH methods to complex networks (De Silva and Ghrist, 2007; Horak et al., 2009; Petri et al., 2013). The network approach is used to reconstruct the phase space, which is a preliminary and certainly necessary step for the analysis of the topological structure of flows from data.

The PH framework to obtain families of nested cell complexes from point clouds has only been mentioned in passing in this review article for the sake of brevity; it should be taken into account, though, as an important of branch of computational topology that is continuously providing us with solutions to algorithmic problems being faced in chaos topology and the climate sciences. So far, the complex network community seems to be lacking a dual object such as the templex to deal with nonstationarity. Finding such an object that captures the spatial structure and is, in addition, endowed with another object that captures the flow structure on the spatial object appears to be a worthwhile challenge.