Noise-driven topological changes in chaotic dynamics

The use of topological concepts to describe complex physics goes back at least to Lord Kelvin, who proposed that atoms were knotted vortices swirling in the ether, an invisible medium believed at that time to fill the surrounding space.⁶ Though incorrect, Kelvin’s vision stimulated active research in knot theory^7,8 and eventually led to the awareness that knots and related tangled structures do actually form in various physical phenomena and can have a pivotal, albeit still poorly understood, influence on turbulent fluid dynamics,^9,10 quantum field theory,¹¹ and magnetic fields,^12,13 to cite but a few areas of physics. In his seminal work, Moffatt¹⁴ showed that helicity measures the total knottedness and linkage of a flow and that it is an invariant in ideal, viscosity-lacking fluids, like liquid helium. In viscous fluids, helicity fluctuates, and knots can experience transformations or unravel.¹⁵ 

In dissipative chaotic dynamical systems, applications of algebraic topology,¹ which studies algebraic invariants that classify topological spaces up to homeomorphism,¹⁶ have also attracted much attention over the past three decades.^17–19 Branched manifolds,^3,20 in particular, have played a central role in the characterization of the coarse-grained features of chaotic dynamics. The concept was anticipated in Lorenz’s famous 1963 paper,²¹ where he remarks that the flow on the attractor’s “surface” passes “back and forth from one spiral to the other without intersecting itself.” This surface is topologically equivalent to what is now called a branched manifold.²² 

A branched manifold has the property that each point in it has a neighborhood that is homoeomorphic to either a full ball or a half ball, thus allowing for the presence of boundaries and, therefore, of branches.¹⁹ Branched manifolds may thus have branches or not: they constitute a broader category, to which boundaryless manifolds belong as well. These objects provide a simple but powerful “cartoon” of the stretching and folding responsible for the creation of strange attractors and are also referred to as a template or a knot holder.

If the system is deterministic, the topology of its branched manifold is an invariant in phase space. Metric invariants, such as correlation dimension,²³ and dynamic invariants, such as Lyapunov exponents,²⁴ are very useful in describing and classifying dynamical systems. The topology of a branched manifold, though, provides more specific information on the mechanisms acting in phase space to shape the flow and, therefore, on how to model the system’s dynamics.¹⁸ In the present paper, we extend the applications of branched manifolds to chaotic systems with time-dependent forcing, in particular to noise-driven ones.

Interest in such nonautonomous and random dynamical systems (NDSs and RDSs) has been greatly stimulated recently—among many other areas of application, as discussed in Refs. 25–27—by the effects of anthropogenic forcing on the climate system, on the one hand, and by the multiscale character of the system’s intrinsic variability, on the other.^28,29 Several authors—see Ref. 30, and further references therein—have addressed the latter problem by stochastic parameterizations of unresolved scales in high-order models, and more realistic low-order models have thus been obtained as well.^31–33 

With this motivation in mind, Chekroun et al.⁵ studied interactions between noise and nonlinear effects from a qualitative dynamical viewpoint, through high-resolution numerical computations. The pullback attractors (PBAs) of RDS theory^34,35 provide a natural framework for doing so and are the mathematically rigorous counterpart of the heuristically defined snapshot attractors of nonlinear physics.³⁶ Further theoretical details on PBAs and RDSs appear in  Appendixes A and  B.

In the random, RDS setting, a PBA is computed by following an ensemble of trajectories, each driven by the same noise path. This way, the well-known smoothing effect of noise due to a single very long integration of the stochastic system disappears in the pullback approach and the fractal structure of the chaotic dynamics is fully captured when the unperturbed nonlinear system’s dynamics is itself chaotic.

PBAs thus provide a natural way to reveal the time-dependent stretching and folding caused by the interactions between noise and nonlinearities, when the unperturbed dynamics is chaotic.⁵ In their study of the stochastically perturbed Lorenz²¹ model’s PBA, dubbed LORA, Chekroun et al.⁵ pointed out that the stretch-and-fold mechanisms give rise to three distinct components of the flow in phase space: a smooth evolution tied to the Lorenz model’s deterministic convection; a pervasive, local “jiggling” due to the roughness of the driving multiplicative noise; and sudden deformations of the random PBA’s overall support. The static picture of the strange attractor generated by deterministic chaos^37,38 is thus replaced, in the presence of noise, by a dynamic version that is even stranger, and that we will refer to as noise-driven chaos in what follows.

Noise-driven chaos offers therewith a self-consistent framework for extending the concepts and tools of chaos topology toward more realistic situations in fluid dynamics, the climate sciences, and elsewhere, to advance the understanding of complex nonlinear dynamics in the presence of noise. The purpose of this article is to examine noise-driven chaos through the lens provided by algebraic-topology tools.

Topological data analysis studies the structure of large data sets by focusing on their connectivity.^39–41 These methods are used more and more in numerous and diverse research areas, including image processing,⁴² the spread of social and biological contagions on networks,⁴³ percolation theory,⁴⁴ genomics and evolutionary dynamics,^45,46 and protein structure.⁴⁷ Even though topological data analysis refers broadly to all methods using notions of shape and connectivity to find structure in data,⁴⁰ the term is nowadays being used more narrowly to refer to a particular method called persistent homology (PH).^48,49 This method has attained great popularity in mathematical imaging and vision.

We study here LORA’s changes in time by applying a topological data analysis methodology to describe the topology of branched manifolds using homologies (BraMAH).⁴ BraMAH was conceived three decades ago to replace the braids or knots made with reconstructed orbits, as done by Gilmore and co-workers¹⁸ up to that point. Braids or knots dissolve into trivial objects beyond three dimensions and, besides, knotted or braided orbits are difficult to reconstruct from data sets that are not noise-free. These shortcomings motivated the search for knot- or braidless and orbitless methods, as stated by Natiello and Solari.^50,51

Homologies, which had been used to analyze experimental dynamical systems involving boundaryless manifolds,¹⁹ provided an avenue for achieving this goal. The topological description provided by the Euler characteristic and Betti numbers associated with homology groups,¹⁹ however, did not suffice to correctly identify a branched manifold from a point cloud. BraMAH went further to ensure the extraction of all the pertinent information, including torsions, branch crossings, and weak boundaries.⁵² Applied at the beginning to noisy experimental data⁴ and to the analysis of four-dimensional chaotic solutions of a Shilnikov type,⁵³ BraMAH has been recently refined to handle data sets generated by Hamiltonian systems,⁵² as well as by nonautonomous ones.⁵⁴ 

In this work, we take a significant step forward, from the purely deterministic to the stochastic realm. Topology is no longer an invariant for RDSs but it is still expected to encode the stretching and squeezing mechanisms shaping the flow in phase space. We define here an RDS’s time-evolving branched manifold locally as an integer-dimensional set in phase space that provides a robust skeleton of the point cloud associated, at each instant, with the invariant measure supported by its random attractor; see  Appendix C for the theoretical underpinnings of the latter concepts.

The point cloud in BraMAH is first decomposed into subsets of points, called cells, which are glued together into a cell complex that forms the coarse skeleton of the structure upon which the point cloud lies. In algebraic topology, a cell complex is a finite set of cells of different dimensions such that the interiors of the cells do not intersect.⁵⁵ The BraMAH cell complex construction method is designed to approximate a point cloud lying on a branched manifold. This way, the cell complex’s topological features unveil the topological signature of the branched manifold. The BraMAH methodology is briefly summarized and further refined in Sec. II and  Appendix D provides further details.

The numerical results and their robustness are described in Sec. III. The key results are the following:

• the BraMAH methodology captures LORA’s time-evolving topological structure by constructing the branched manifold at each instant in time;

• LORA’s branched manifold differs from the static one obtained for the deterministic Lorenz model’s strange attractor; and

• the noise-driven model’s otherwise gradually evolving branched manifold exhibits sharp transitions at discrete points in time.

These results are further discussed in Sec. IV.

Pursuing our goal of examining noise-driven chaos through the lens of algebraic-topology tools, we first provide some background on PH in general and on BraMAH in particular.

The PH method measures connectivity in terms of intersections between $ϵ$-sized balls centered at each point in the point cloud. Connection rules to construct cell complexes can be defined in terms of these intersections as the value of $ϵ$ grows from zero to some value, which is related to the full size of the cloud. The set of rules establishing when and how to connect the points in the cell complex receives the name of filtration—such as the Čech or Vietoris–Rips filtrations^48,56—and $ϵ$ is called the filtration parameter. For a given $ϵ$-value, homologies can be used to identify holes (regions in the point cloud that are void of connections).

Examining how the cell complexes change with $ϵ$⁠, holes can be described as having a lifetime, i.e., an $ϵ$-range within which they “persist.” The result of a PH study is often expressed using $ϵ$ as the horizontal axis of a barcode, which furnishes a connectivity fingerprint of that particular point set. When applied to point clouds corresponding, for instance, to a data set generated by the equations of the Lorenz attractor,⁵⁷ one chooses a filtration method and obtains an $ϵ$-dependent portrait of the point set. For a certain $ϵ$-range, the PH study identifies the two holes of the butterfly’s wings, along with many other holes, which are less persistent. But is any point cloud with two persistent holes in a certain filtration range identical to the Lorenz attractor’s branched manifold? Of course not. Counting the holes and pondering their persistence may be indicative of the attractor’s identity but not conclusive. The information retrieved with a PH study contains valuable scale-dependent details within a point cloud but the approach is not at all tailored to obtain the geometry-independent invariant of a deterministic dynamical system in phase space.

Recently, Strommen et al.⁵⁸ have adapted PHs to yield a computationally tractable identification and classification of multiple weather regimes^38,59–61 on the basis of the topological structure of the atmospheric flow in a low-dimensional, high-variance subspace of a given data set. Given the large number of classification schemes—based on phase space clustering and temporal-persistence criteria, among others^62,63—such a topological approach could provide additional insights into the still poorly understood mechanisms of the complex organization of the atmosphere’s “regime diagram.”⁶⁴ The weather regime classification of Ref. 58 resorts to bifiltrations in order to incorporate density as an extra filtration parameter into their PH study.

What does BraMAH do? BraMAH builds a cell complex that serves as an algebraic representation of the branched manifold underlying a point cloud in phase space, and it computes the full set of topological properties that describe it. In other terms, BraMAH uses homologies as a means to reconstruct a branched manifold. The BraMAH methodology requires approximating the structure of the branched manifold that supports the point cloud by a cell complex, as defined in algebraic topology.⁵⁵ The structure of this invariant manifold can be approximated with “building blocks.” In algebraic topology, these building blocks are Euclidean closed sets (segments, disks, etc.), called $n$-cells, where $n∈N0$⁠. A point is a $0$-cell, a line segment joining two points is a $1$-cell, a polygon is a $2$-cell, and so forth. Constructing and assembling these cells according to the BraMAH rules detailed in  Appendix D builds a cell complex that is the branched manifold’s skeleton.

The topological structure of this cell complex is studied by using homology groups⁵⁵ and orientability chains.⁴ The homology groups $Hk$ identify the $k$-dimensional holes (⁠$k$-holes) of a topological space of dimension $n$⁠, where $k∈N0$ and $k≤n$⁠. The group $H0$ identifies the connected components (⁠$0$-holes), $H1$ the cycles (⁠$1$-holes), $H2$ the cavities (⁠$2$-holes), and if $k≥2$⁠, $Hk$ the hypercavities (⁠$k$-holes). By construction, a BraMAH cell complex is uniformly oriented. The orientability chain indicates the number and location of torsions in the cell complex. For further details of the BraMAH method, see  Appendix D.

As an example, we provide here a BraMAH analysis of the classical Lorenz²¹ attractor. The strange attractor and the BraMAH cell complex associated with its branched manifold are shown in Fig. 1. The cell complex presents one connected component, i.e., $H0∼Z$⁠, and two homologically independent $1$-holes associated with each wing of the butterfly, i.e., $H1∼Z2$⁠; there are no torsions, i.e., $O1∼∅$⁠, and no enclosed cavities, i.e., $H2∼∅$⁠. The manner in which the two wings meet at the regular saddle on the $z$ axis can be assessed from the location of the $1$-holes, taking into account how they interweave without twisting.

In contrast to the PH approach, BraMAH’s target is not counting holes and examining their persistence in terms of filtration rules that are established in terms of one or more free parameters. The manner in which BraMAH constructs a cell complex relies neither on $ϵ$-sized balls around every point of the point cloud nor on how such balls intersect. Knowing that the point cloud being analyzed lies on a branched manifold, one builds a single cell complex and not a nested family of parameter-dependent cell complexes, as done in the PH approach. There is no analog of $ϵ$ or of a filtration rule in BraMAH since a branched manifold is not parameter-dependent in the PH filtration sense. This property of the dynamics generating the point cloud is utilized to produce a cell complex in which each cell represents a large subset of points that jointly approximate a locally Euclidean set, whose dimension respects the local dimension of the branched manifold. Some of the branches may have torsions that are important for the correct identification of an attractor; hence BraMAH also distinguishes, for instance, a standard strip from a Möbius strip or a Möbius strip from a Klein bottle. The method’s specificities hinge on the cell complex construction process as well as on the way it extracts the topological properties associated with the complexes so constructed.

BraMAH has been applied so far to speech data,⁴ to slow-fast systems,^20,53,65 and to Lagrangian dynamics in fluid flows.⁵⁴ Here, it is being applied for the first time to data sets generated by noise-driven chaos. When evaluating a time-dependent structure like LORA, one has to ask whether its topology is time-dependent. In order to answer this question, we take two steps: (1) generate the point clouds that approximate LORA at successive time instants, called “snapshots”^5,,66 and (2) compute the topological properties associated with each snapshot using BraMAH.

In step (1), the point cloud approximating LORA is generated by solutions of the stochastic Lorenz model introduced by Chekroun et al.⁵ In (1), the equations of the classical Lorenz²¹ deterministic model are perturbed by a multiplicative noise in the Itô sense,³⁵ with $Wt$ a Wiener process and $σ>0$ the noise intensity,

Here, $r=28$⁠, $s=10$⁠, and $b=8/3$ are the standard parameter values for deterministically chaotic behavior.

The time-dependent sample measures $μt$ associated with the system (1) are probability measures for the population density of any ensemble of initial data driven by the same noise realization until time $t$⁠, after removal of the transient behavior. Mathematically, these measures are of Sinai–Ruelle–Bowen (SRB) type,^67–69 i.e., they are supported by the foliation of unstable manifolds that structure the random attractor.⁵ 

A numerical estimation $μt^$ of such a measure can be computed at any time instant $t$ by a pullback approach, i.e., by letting a large set of $N0$ initial points ${xj(0)=(xj,yj,zj)(t=0):j=1,…,N0}$ “flow” in phase space from the remote past until time $t$⁠, for a fixed noise realization $ω$⁠. The convergence of the sample measure’s approximation $μt^=μt^(N0)$ is studied as the number $N0$ of initial points increases; it is observed herein for $N0≃108$⁠.

Each point within a given cloud at time $t$ is mapped to a value of $μt^$ that is obtained by averaging over a volume surrounding that point: higher or lower $μt^$-values correspond to more or less populated regions of the random attractor. We wish to characterize the topology of the point cloud’s most populated regions but also to ascertain this topology’s robustness. In order to do this, a threshold for $μt^$ must be selected, as discussed below.

Solutions of the stochastic Lorenz model (1) generate data sets in the form of three-dimensional point clouds ${xi∈R3:1≤i≤N0}$ with $N0=108$⁠. In order to obtain the robust skeleton of the point cloud associated, at each instant, with the invariant measure supported by LORA, these point clouds are sieved using a threshold value $n¯$ in the sample measure $μt^$⁠. This yields point clouds that reveal the most populated regions of the cloud.

The $n¯$-values are chosen within limits that guarantee that the connectivity of the point cloud is not lost so that there is only one connected component. To this end, we define $n¯max$ as the $n¯$ value for which there is still one connected component, i.e., $H0∼Z$⁠, and normalize all $n¯$ values by the value of $n¯max$⁠; for notational simplicity, we keep the symbol $n¯$⁠, though, for the normalized value as well.

Figure 2 shows how LORA at $t=40.27$ flakes off as the point cloud is sieved. We use the term “sieve” to avoid confusing this process with that of a filtration, as defined in PH. Notice that when $n¯>1$⁠, the point cloud is broken into three sufficiently populated connected components (⁠$H0∼Z3$⁠).

In order to study LORA’s changing topology at time $t$⁠, a criterion must be adopted to select a value $n¯t∗$ of $n¯$ that adequately represents the branched manifold for that value of $t$⁠. For each $n¯$-value, the BraMAH cell complex built from ${xi:μt^(xi)≥n¯}$ can be analyzed, yielding a set of 1-holes, among other results. As one can see in Fig. 2 for $t=40.27$⁠, some of the 1-holes are robust as the value of $n¯$ increases. Based on this observation, we build a robustness plot, using $n¯$ as the horizontal axis and the number of $H1$ generators—i.e., of 1-holes—on the vertical axis, which piles up in the order of their appearance. Normalizing the lengths $ℓ¯$ of the bars against the longest bar $ℓ¯max$ of the plot, we define the most robust 1-holes as those for which $ℓ¯/ℓ¯max≥0.5$⁠. With this criterion of robustness in mind, we finally require $n¯t∗$ at each time $t$ to be a value of $n¯$ at which the $H1$ generators of the BraMAH cell complex coincide with the most robust $1$-holes.

Note that the density here is not a filtration parameter in the PH sense, i.e., the value of $n¯$ is not involved in the rules used to construct the BraMAH cell complex. For this reason, the robustness plot should not be read as a PH barcode: at each $n¯$-value, BraMAH uses a point cloud that differs in its size and elements to construct the cell complex, whereas in PH calculations, the single point cloud is under consideration and the horizontal axis sweeps parameter values that regulate the construction of nested cell complexes for the same point cloud.

Figure 3(a) shows the robustness plot for $t=40.27$ with noise intensity $σ=0.3$⁠. Fourteen $H1$ generators (1-holes) are represented and colored. The most robust holes are labeled with numbers (1)–(4). In panel (b), we show an $(y,z)$ projection of the sieved point cloud with $n¯=0.167$⁠, in which holes labeled (1)–(3) can be identified. In (c), we show a projection onto the plane $−2.25x−20y+6z=0$ of the same point cloud where hole labeled (4) can be visualized. The colors of the holes correspond to those of the bars in (a).

Following the criteria formulated in Secs. II and III A and illustrated in Figs. 2 and 3, we analyze the topology of three successive LORA snapshots for a given, fixed noise intensity of $σ=0.3$⁠. These three points in time—$t=40.09,40.18,$ and $40.27$—have been chosen visually from the LORA video provided in the supplementary material of Ref. 5. The video clearly suggests a gradual, smooth evolution of the invariant measure, except for sudden changes at discrete times. Two such points in time appear to lie in the open intervals $(40.09,40.18)$ and $(40.18,40.27)$⁠, both of which are quite short with respect to the characteristic nondimensional time of $1.0$ for the deterministic Lorenz model.²¹ 

At these three points, we find that the local dimension of the branched manifold is the same as in the deterministic case: locally, LORA remains a two-branched manifold and no torsions are ever observed, i.e., $O1∼∅$⁠. Noise adds neither connected components nor cavities to the manifold so that the groups $H0∼Z$ and $H2∼∅$ are the same as in Fig. 1. The topological distinctions between LORA and the classical Lorenz²¹ attractor are thus solely present in $H1$⁠. One notices immediately, though, that the number of $1$-holes changes from one snapshot to the next, with holes created or destroyed by the noise in the course of time.

This situation is illustrated in Fig. 4 for the three successive snapshots, at $t=40.09,40.18,$ and $40.27$⁠; for each snapshot, we use the value of $n¯∗$ specified in the caption, which leads to LORA’s branched manifold at time $t$⁠. The sieved point clouds are given in the three top panels, while the corresponding BraMAH cell complexes are presented in the three bottom panels. The branched manifold’s structure undergoes large changes not only in the number of $1$-holes but also in the organization of the branches from one snapshot to the next.

The BraMAH analysis herein thus fully supports the intuition provided by the LORA video of Chekroun et al.⁵ on sudden topological changes associated with those visible in the invariant measure. Further work clearly needs to complete the present analysis in terms of narrowing down even further the exact location of the topological tipping points and being more exhaustive in cataloging the possible changes in the branched manifolds, for LORA and for other stochastically perturbed chaotic systems.

We have topologically analyzed here a paradigmatic random attractor⁵ associated with the Lorenz²¹ convection model. To accomplish this, we have extended the deterministic concept of a branched manifold² by defining 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 instant $t$⁠. This new definition of a branched manifold is used to perform what is—to the best of the authors’ knowledge—the first topological analysis of a time-evolving chaotic attractor.

Our topological analysis is based on a topological data analysis methodology called BraMAH,^4,52–54 originally designed to study the topological structure of chaotic data in a deterministic setting.⁴ This homology computation tool leads from a point cloud to the full set of topological features that characterize a branched manifold. BraMAH differs from other generic topological data analysis methods by its using explicitly the fact that the point cloud under study lies on a branched manifold and by detecting relevant features such as branch torsions or weak boundaries if they exist.

In the presence of noise-driven chaos, the branched manifold is time-evolving and must, therefore, be analyzed time-wise, i.e., using point clouds that approximate the random attractor by successive snapshots. The branched manifolds so constructed by BraMAH provide topological portraits at different stages in the life of LORA, the Lorenz model’s random attractor.

Our study shows a marked difference between the deterministic case and the noise-driven one. The stochastically perturbed system’s random attractor, dubbed LORA, presents a much richer structure than the deterministic strange attractor, with a topology that also changes drastically in time. The framework introduced in this article to characterize such changes in topological features appears to hold promise for the understanding of topological tipping points in general.

A fairly straightforward BraMAH application to the climate sciences might clarify the following quandary of so-called subseasonal-to-seasonal prediction.⁶³ The quandary deals with the role of intermittent vs oscillatory low-frequency variability in the atmosphere.^60,62 Low-frequency variability refers to phenomena whose characteristic time of 10–100 days exceeds that of typical midlatitude storms of 5–7 days (see, for instance, Chap. 6 of Ref. 38).

Such phenomena include the so-called blocking of the westerlies and intraseasonal oscillations with periodicities of 40–50 days; the former is intermittent, or particle-like, the latter oscillatory or wave-like.^29,60 The quandary is which type of phenomenon dominates low-frequency variability and thus could contribute most to subseasonal-to-seasonal prediction (cf. Fig. 1 of Ref. 70).

Blocking has been recently studied in Lucarini and Gritsun⁷¹ by the unstable periodic orbits methodology of Gilmore.¹⁸ It would appear that the BraMAH methodology proposed herein could address, more generally and efficiently, the waves-vs-particles quandary of low-frequency variability and how it might be affected by global change.⁷⁰ 

More broadly, tipping points have been given a precise definition in the climate sciences as a generalization to NDSs and RDSs of classical bifurcations in autonomous systems,^28,72 and they are being actively pursued.^73–78 Topological tipping points seem to be a further generalization of the concept that could help us apprehend sudden and drastic changes in time of model behavior as well as drastic changes due to mean forcing intensity.

The gradual change of atmospheric concentrations of greenhouse gases and aerosols may continue to modify global temperatures in a fairly smooth way,⁷⁹ as the case has been so far, although it might also lead to dynamical tipping points, as suggested by the above-cited papers. On the other hand, the intrinsic noise associated with cloud processes on small space and time scales affects the entire climate system through their interaction with dynamic and radiative processes on larger scales and may lead to hitherto unsuspected topological tipping points.²⁹ 

Topological tipping points seem to be a further generalization of the concept that could help us apprehend sudden and drastic changes in time of model behavior as well as drastic changes due to mean forcing intensity, whether deterministic or stochastic.

We have concentrated throughout much of this paper on problems related to the climate sciences, as had Lorenz in his pioneering paper.²¹ With all due modesty, it is not unlikely—considering the great generality of topological methods^1,18—to expect the results obtained herein to have some applicability to other areas of the physical, life, and socio-economic sciences.

We thank many colleagues who asked good questions at several talks and posters given on this work (see, for instance, Refs. 80 and 81). It is a pleasure to thank D. Faranda and two anonymous reviewers for their constructive comments that helped improve the presentation of our results. This work is supported by the French National program LEFE-NOISE (Les Enveloppes Fluides et l’Environnement) and by the CLIMAT-AMSUD 21-CLIMAT-05 project (D.S.). G.D.C. gratefully acknowledges her postdoctoral scholarship from CONICET and wishes to thank Juan Ruiz and the CIMA computing staff for support with the computations. This work was also partially supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 810370), by a Ben May Center Grant for theoretical and/or computational research and by the Israeli Council for Higher Education (CHE) via the Weizmann Data Science Research Center (M.D.C). The article is TiPES Contribution No. 46; this project has received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 820970 (M.G.).

The data that support the findings of this study are available from the corresponding author upon request.

We summarize herein some pertinent facts on nonautonomous systems of ordinary differential equations,

considered in the framework of nonautonomous dynamical systems (NDSs); $F$ denotes a smooth time-dependent vector field that governs the time evolution of the state $x$ in a phase space $X$⁠, taken for simplicity to be the Euclidean space $RN$⁠.

Once existence and uniqueness are guaranteed, one can assign to this NDS a solution map $Φ(t,s)$⁠, which provides a two-time description of the motion: the time $s$ when the system was initialized, and the time $t≥s$ of the system’s current state. Thus,

denotes the solution of Eq. (A1) at time $t$⁠, when initialized at $x(t)=x0$ for time $t=s$⁠. In the autonomous case, only the time interval $t′=t−s$ separating $s$ and $t$ matters and $Φ(t,s)$ reduces to a standard flow $Φ(t′)$⁠.

In the case of forced and dissipative systems, such as the climate system,^29,38 one can define a collection of subsets called a pullback attractor (PBA).^{5,35,76,82,83}

The subsets ${A(t)}$ are often called snapshots.^5,36,84,85 According to (1) and (2) above, the indexed family of snapshots is invariant under the system’s dynamics and it attracts at each time $t$ all compact initial subsets $B$ from the remote past (see also Fig. A.1 of Ref. 83 for a simple illustration).

In physical systems, such as those encountered in the climate sciences, random time-dependent forcing is often present.²⁹ When that is so, it becomes necessary to model this type of system using stochastic differential equations (SDEs).³⁵ 

In the theory of RDSs, random PBAs are known as random attractors and they can be constructed in an extended phase space composed of the phase space $X$ and a probability space associated with the paths of the driving noise. A probability space $(Ω,F,P)$ is a three-tuple, where $Ω$ is the sample space; $F$ is the event space, formulated as a $σ$-algebra; and $P$ is a probability measure on $F$ (see Appendix A of Ref. 35).

The probability space is then endowed with a time-dependent shift $θt$⁠. In the case of a stochastic–dynamic system driven by a Wiener process, as is the case here, this shift is defined on $Ω$ according to $Ws(θtω)=Wt+s(ω)−Ws(ω)$⁠.³⁵ With this mapping $θt$ in hand, the noise realization $ω$ evolves in time, and one can define a cocycle to describe the evolution of the state $x$⁠.

A mapping $φ:R×Ω×X→X$ has the cocycle property when $φ(t,ω)=φ(t,ω,⋅):X→X$ satisfies the following conditions:³⁵ 

• $φ(0,ω)x=x$⁠, for all $x∈X$ and $ω∈Ω$⁠, and

• $φ(t+s,ω)=φ(t,θs(ω))°φ(s,ω)$⁠, for all $s,t∈R$ and $ω∈Ω$⁠, where $°$ denotes the composition operation for mappings of $X$⁠.

Property (i) just sets the initial state of the cocycle, while property (ii) states that a cocycle is an expression of the existence and uniqueness of solutions, in the sense that going from a copy of $X$ at time $0$ to one at time $s$ and from there on to one at time $t+s$ is the same as going directly from $0$ to $t+s$ (see also Fig. A.2 of Ref. 83). This cocycle property is satisfied for a broad class of SDEs like those of interest here (see Ref. 35).

Mathematically, given an SDE with the right properties, the probability space $(Ω,F,P)$ equipped with the collections of shifts $θ={θt}t∈R$⁠, and its associated cocycle $φ$⁠, form what is called an RDS $(φ,θ)$⁠, also called sometimes an RDS $φ$ over $θ$⁠.

The evolution of a stochastic–dynamic system can thus be modeled by an RDS, while its associated random attractor ${A(t;ω)}t∈R$ provides the natural extension of a PBA (as defined in  Appendix A above) to the random setting in which each individual PBA depends on the specific noise realization $ω∈Ω$⁠. The resulting family $⋃ω∈ΩA(t;ω)$ of random compact sets provides a complete description of all the system’s possible states that are likely to be observed at time $t$⁠.

A number of interesting properties follow from the fact that the RDS $(θ,φ)$ has a random attractor. One of these is the existence of invariant measures of $(θ,φ)$⁠, in the sense of RDS theory. In this Appendix, we briefly clarify this notion and discuss the properties of these measures.

To do so, recall that any RDS $(θ,φ)$ generates a skew-product semiflow ${Θ(t)}t≥0$ on $Ω×X$ by the formula

The cocycle property for $φ$ is equivalent to the semigroup property for $Θ(t)$⁠, namely, $Θ(t+s)=Θ(t)Θ(s)$⁠. In what follows, we denote by $B$ the $σ$-algebra of Borel sets in $X$ (see Ref. 35). We have then the following definition.

It is known that any probability measure $μ$ on $(Ω×X,F×B)$ possesses a disintegration or factorization,³⁵ given by a function $(ω,B)↦μω(B)$ from $Ω×B$ into the interval $[0,1]$ such that:

• for any $B∈B$⁠, $μω$ is $F$-measurable;

• there exists a measurable set $Ω′$ in $Ω$ such that $P(Ω′)=1$ and $μω$ is a probability measure on $(X,B)$ for all $ω$ in $Ω′$⁠; and

• for all $f$ in $Lμ1(Ω×X)$⁠, we have

  $∫Ω×Xf(ω,x)μ(dω,dx)=∫Ω(∫Xf(ω,x)μω(dx))P(dω).$

  (C2)

The disintegration $μω$ is unique $P$-almost surely and it is also called a sample measure (see Ref. 68). The invariance property (i) of Definition 2 translates into $φ(t,ω)μω=μθtω$⁠, in terms of sample measures.

When an RDS $(θ,φ)$ possesses a random compact attractor, then it supports every invariant measure, i.e., $μω(A(ω))=1$ for almost all $ω$ in $Ω$⁠. In this case, the sample measure possesses a useful interpretation. To understand it, recall that, roughly speaking, the random attractor $A(t;ω)$ determines the portions of the phase space $X$ onto which any bounded set $B$ is mapped at time $t$⁠, when $B$ is propagated by the cocycle $φ$ from a remote past, for a given noise realization $ω$⁠. The sample measure $μθtω$ supported by the random attractor $A(t;ω)$ provides, therewith, the spatiotemporal probability distributions of the portions of the phase space $X$ occupied by the RDS, at time $t$ and for the noise realization $ω$⁠.

For a given stochastic–dynamic system, the set of possible invariant measures is rather large. This raises the question whether a particular class of invariant measures is “naturally chosen” by the dynamics. Physical measures are sample measures of special interest in this respect.^5,29 A probability measure $μ$ is physical if, for any continuous observable $ψ:X→R$⁠, the time average equals the ensemble average for almost all initial data $x0$ that lie in a Lebesgue-positive set $Bμ$⁠, called the basin of attraction of $μ$ [see Eq. (5) of Ref. 5].

Sinai–Ruelle–Bowen (SRB) measures^67–69 form a closely related class of probability measures. A probability measure $μ$ is an SRB measure if its conditional measures on the unstable manifolds are absolutely continuous with respect to Lebesgue measure. For many dynamical systems, the class of physical measures coincides with that of SRB measures; small differences may exist, however, between the two concepts for certain systems. Regarding the stochastic Lorenz system (1), Chekroun et al.⁵ have rigorously shown in their Appendix C that LORA supports a random SRB measure when the Kolmogorov operator associated with (1) is hypoelliptic and its leading Lyapunov exponent is positive. Both conditions are met for the parameter values and stochastic forcing used in Eq. (1). Strong numerical evidence provided in Ref. 5 suggests that the measure to which the point clouds used herein converge must be physical, too.

Note that the existence of an SRB measure $μ$ does not guarantee its uniqueness and that two such measures $μ≠ν$ may also have different basins of attraction, $Bμ≠Bν$⁠. The extensive numerical calculations in Ref. 5 and herein have given no indication, though, of such nonuniqueness. Still, the uniqueness of LORA’s random SRB measure or of its physical measure has not been proven rigorously, to the best of the authors’ knowledge.

BraMAH computes the topology of point clouds whose points are locally distributed on a branched manifold. An $m$-manifold is a topological space with the property that each point has a neighborhood that is homeomorphic to either a full $m$-ball or a half $m$-ball, with $m$ in $N$⁠.¹⁹ In nonlinear dynamics, an $m$-branched manifold is a mathematical object embedded in a phase space of dimension $d$ (⁠$m≤d$⁠) that is a manifold everywhere but at the tear points of the branches.¹⁸ We refer to $m$ as the local dimension of the branched manifold.

Subsets of points of such a point cloud can be locally approximated by $m$-disks. This decomposition into point subsets or “patches” is called “patch decomposition.”⁴ A BraMAH complex is a cell complex constructed so that each patch is associated with a cell in the complex.

To define cell complexes in general, we start with a single cell. A cell of dimension $k$ (⁠$k$ in $Z+$⁠) or $k$-cell is a set that can be mapped through a homeomorphism into the interior of a $k$-disk so that the boundaries of the image are divided into a finite number of lower-dimensional cells, called faces. A cell complex $K$ is a finite set of cells, such that (a) their faces are also elements of the complex and (b) the interiors of two cells never intersect. $K$ is said to be an $h$-complex or to have dimension $h$⁠, if its highest dimensional cell is an $h$-cell.

How does BraMAH construct cells from patches? One first partitions the point cloud ${xi∈Rd:1≤i≤N}$ into overlapping patches that are homeomorphic to the interior of an $m$-disk. A patch ${xi=(xi,1,xi,2,…,xi,d),i=1,…,Nc}$ is built around a center $x0=(x0,1,x0,2,…,x0,d)$ by searching for the largest number $Nc$ of points around $x0$ that constitute approximations to a Euclidean set of dimension $m≤d$⁠. In order to calculate $Nc$⁠, candidate sets ${xi,i=1,…,nc}$ are computed with its elements sorted by distance from $x0$⁠, and with $N{min}≤nc≤N{max}$⁠. We use $N{min}$ typically $2%$ the total size of the cloud, and $N{max}$ an order of magnitude higher than $N{min}$⁠.

Given a candidate set, the points $xi$ within the ball with center $x0$ and radius $r$ are represented by the neighborhood matrix $X∈Rnc×d$⁠,

A local coordinate system centered at $x0$ is provided by the singular vectors of $X$⁠, with the singular values describing the distribution of the points inside the ball centered at $x0$⁠. For a patch that is approximately lying on an $m$-disk in $Rd$⁠, the local singular spectrum of $X$ has $m$ singular values that scale linearly with $nc$ as $r$ is increased. This property holds while the effects of the manifold’s curvature are negligible.^19,86 The remaining $(d−m)$ singular values, which measure the deviation from the tangent space, will scale as $rℓ$ with $ℓ≥2$⁠.

Using this rule, the $m$ relevant singular values and vectors that span the tangent space approximating the patch under consideration can be identified. The value of $Nc$ is obtained when the $m$ relevant singular values—as functions of $nc$ subject to $N{min}<nc<N{max}$—exhibit the best linear regression coefficient.

This procedure is applied repetitively until every point of the point cloud belongs to at least one patch. The patch axes are chosen so as to favor patch overlapping, in order to keep track of the gluing prescriptions between them. Each patch is transformed into a cell by using convex hulls, and the singular vectors are used to orient each cell so that neighboring cells have the same orientation, thus resulting in the complex having a uniform orientation. A BraMAH complex is a cell complex that results from applying these steps to a $d$-dimensional point cloud.

The next step in the BraMAH methodology is to compute the topological properties of the cell complex obtained from the algorithm above. For any given complex, a $k$-chain is a linear combination of $k$-cells with integer coefficients. The algebra of these chains allows one to describe the connectivity of the cells at each $k$-level. A $k$-hole is a closed chain, called a $k$-cycle, that is not the border of any higher-dimensional cell.

Our approach uses the labeled list of $0$-cells of the BraMAH complex in order to build a boundary matrix and extracts the linearly independent rows of it. Then, it computes the null spaces of the transpose of the boundary matrix and expresses the $k$-borders in terms of the $k$-cycles to determine which $k$-cycles are homological to others. The $k$-cycles that are homologically independent are appended to $Hk$⁠. The integer multiples found in the chain that sums up all the $k$-borders of the complex are used to form the orientability chain, and the orientability chain is in turn used to compute and locate torsions, as well as weak boundaries.

The results of this final phase of the method are the homology groups of the BraMAH complex expressed in terms of their generators ${Hk:k=0,…,h}$⁠, spelled out in terms of the $0$-cells. The advantages of the BraMAH methodology over other homology computation methods are the following:

• A BraMAH cell complex is uniformly oriented—thus allowing for the computation of orientability properties—and it has a dimension $h$ that agrees with the local dimension $m$ of the branched manifold; the latter is computed locally using the neighborhood matrix $X$ so that no cells with unnecessarily high dimensions are created.

• This cell complex is non-simplicial by construction, i.e., the number of sides of a cell is free, and each $h$-cell is constructed from a large subset of points; doing so yields a complex with a number of $h$-cells that is significantly lower than the number of points in the original point cloud, while PH-type methods—e.g., the Vietoris–Rips complex—tend to produce complexes with a number of cells that are vastly larger.^48,56

• The $0$-cells of the BraMAH complex are associated with a set of points whose coordinates are identified in the original point cloud; these coordinate values can be used to embed the complex in phase space and to ascertain the mutual organization of the branched manifold’s branches.

• The $k$-holes and orientability chains of the complex can be superimposed on the point cloud or—if $d≥4$⁠, on a projection of the point cloud—allowing one to identify the branches, as well as the torsions or twists along them.

As discussed in Sec. II, BraMAH is a topological data analysis approach that does not fall into the category of PH methods. The PH approach focuses on how connectivity properties of a point cloud change with respect to a parameter $ϵ$ that is involved in the filtration rule, i.e., in the rule used to build the cell complex. PH cell complexes are parameter-dependent or multiscale, in the sense that the scale at which points are connected in the cell complex is given by the value $ϵ$ of the filtration parameter: thus, at an $ϵ=0$-scale, all points are disconnected, and they start connecting as the balls grow and intersect.

The PH aim is to describe how the homologies of a cell complex so constructed are affected by the value of this parameter. BraMAH is, on the contrary, filtration-free and parameterless, in agreement with its aim, which is unambiguously identifying a branched manifold’s topology from a point cloud.