A templex-based study of the Atlantic Meridional Overturning Circulation dynamics in idealized chaotic models

Understanding dynamics from time series is one of the most significant goals in nonlinear sciences. In the deterministic framework, this goal can be achieved seeking and finding invariants in phase space. Such invariants can be either metric or topological, but only the topological ones provide information on how the phase space is mapped onto itself under the action of the governing equations. Atlantic Meridional Overturning Circulation (AMOC) models have never been studied from this perspective. Here, we consider the time series produced by two idealized models: the first one describes the self-sustained chaotic oscillations characteristic of the glacial intervals, and the second one incorporates the time-varying position of the edge of sea ice driving the system in and out of this chaotic behavior. Using a templex, i.e., a mathematical object that describes the topology of a deterministic flow in phase space, this paper dissects the dynamics of the AMOC into the different components that decrypt its evolution. The approach is applied both in the autonomous and nonautonomous settings, and compared with previous studies using the concept of pullback attraction.

The term tipping point, popularized in the climate scientific community by Lenton et al.,¹ was originally introduced in a non-scientific context. Its early figurative use emerged in the 1950s, primarily referring to social dynamics.² The term was later repurposed in climate science to describe thresholds in environmental systems.³ This adaptation emphasizes its applicability to complex systems undergoing transitions, such as those observed in our planet’s dynamics.⁴ 

The Earth’s climate system has experienced intervals of glaciation and interglacial intervals. Paleoclimatic records from the last glaciation, including data from the Greenland ice cores and North Atlantic Ocean temperatures, reveal a chaotic variability corresponding to the Dansgaard–Oeschger (D–O) events.^5–7 These events involve slow cooling over several centuries followed by abrupt warming over a few decades. Recent studies have hypothesized that the chaotic variability of the Atlantic Meridional Overturning Circulation (AMOC) can be explained by deterministic self-sustained oscillations between two alternating phases: a slow on-phase during which the circulation strengthens, and a rapid off-phase during which the circulation weakens. These on and off phases are characteristic of relaxation oscillations,⁸ in contrast with the relative stability of the Holocene epoch.⁹ Temporal changes in the position of the edge of sea ice (ESI) have been proposed as a modification in the system’s driving term.¹⁰ 

Chaotic relaxation oscillations, described as an assemblage of slow and fast subsystems in increasingly higher dimensions,¹¹ have been studied from a topological perspective.¹² Understanding the topological structure of a flow in phase space allows researchers to identify the stretching, squeezing, folding, or tearing mechanisms that generate strange attractors.¹³ By analyzing the topology, one can discern how trajectories evolve, aiding in the prediction of transitions between different regimes.¹⁴ Additionally, topological techniques facilitate the interpretation of complex behavior in high-dimensional phase spaces.¹⁵ For three-dimensional systems, the topological structure of an attractor can be thoroughly described using knots,¹⁶ but knots unknot in four dimensions and beyond. To deal with higher-dimensional chaos, it is necessary to use topological tools without dimensional restrictions, such as cell complexes.¹⁷ 

A cell complex is a sort of skeleton made up of cells of different dimensions,¹⁸ allowing for the computation of layered invariants (homology and torsion groups), which can be used to distinguish different topologies. The exploration of phase space with cell complexes began in the nineties.^19,20 Today, we have a special class of cell complexes that approximate branched manifolds:¹² they are called BraMAH complexes—the acronym stands for “Branched Manifold Analysis through Homologies.”²¹ This approach was successfully used, for instance, to analyze speech data²⁰ or to detect Lagrangian Coherent Sets in fluid flows.²² The difficulty with homologies, however, is that they do not consider the fact that a cell complex in phase space is visited by the flow in a specific way, defining circuits that can be homologically equivalent while topologically distinct. Unveiling the non-equivalent paths of the flow upon the structure is not possible without equipping the cell complex with a directed graph (or digraph) that indicates how cells are connected by the flow. The mathematical object incorporating the flow into the picture was introduced in 2022 under the name “templex.”²³ It consists of a BraMAH complex endowed with a directed graph, whose nodes coincide with the highest local dimension cells of the complex. The edges of the graph contain the previously missing flow information. The templex is nothing more than the generalization of the cell complex for the algebraic extraction of the topological properties of a flow in phase space.²⁴ Its properties can be used to identify the fundamental processes that are responsible for the system’s behavior.

Understanding the AMOC dynamics is crucial because of its impact on the planet’s heat distribution.^25,26 There are a number of mechanisms in the literature identified as responsible for the AMOC dynamics. Idealized models coupling ocean, atmosphere, and ice suggest that millennial oscillations are mainly endogeneous and deterministic.²⁷ Sévellec and Fedorov⁸ formulate a low-order model of the AMOC that realistically reproduces the power spectra for different sets of initial conditions. This model is a modification of the Howard–Malkus loop,²⁸ formulated for three variables: two salinity gradients and the strength of the overturning circulation. It is derived from the equations of motion of the ocean circulation in its zonally averaged form in the depth-latitude plane.

To reproduce the interglacial intervals, the model is made more complex.¹⁰ The latitudinal location of the AMOC with respect to the atmospheric precipitation is found to control the stability properties of the system, having a single stable steady state in the interglacials that switches to solutions wandering between three unstable steady states. Consistently with numerical studies, latitudinal shifts of the AMOC are assumed to follow the southward edge of the sea ice. To study this nonautonomous system, a sawtooth function with zero mean and a period of $100000$ years is imposed. Nonautonomous dynamical systems are often studied using the so-called pullback approach,^29–33 i.e., integrating the time-dependent equations of the nonautonomous model over a chosen time interval from a large set of initial conditions, and plotting the probability density function of the system. Sévellec and Fedorov¹⁰ show that computing the pullback attractors for different diagnostic times separated by $5000$ years is enough to unveil that the AMOC dynamics alternates between an interglacial interval and a chaotic burst during the glacial interval. In three-dimensional phase space, the pullback attractor has different shapes depending on the choice of the initialization time and contracts to a single point during the interglacials.

This work provides the first topological analysis of the flow in phase space associated with these deterministic models of the AMOC. The study is undertaken using the templex approach in both autonomous and nonautonomous settings, with the assumption that determinism and continuity imply topological invariance in the absence of false neighbors, i.e., in the absence of different states in phase space that are presented as coincident because of a projection. Adding phase space dimensions to eliminate self-crossings amounts to rewriting the nonautonomous system as a higher-dimensional bounded system in an autonomous form.³⁴ The analysis will be conducted by building the templex directly from the time series provided by the numerical solutions of both model equations, with as many dimensions per case as required, and by computing and comparing their templex properties.

The article is organized as follows. The first section is devoted to describing the templex approach, using the Rössler spiral chaotic attractor as a reference example. The following sections focus on the templex-based analysis of the AMOC models. The three-dimensional system is treated first to characterize the topology of the chaotic dynamics of the glacial intervals reproduced by the autonomous model.⁸ The same system is reconsidered when the sawtooth time-dependent forcing is introduced. The new templex is shown to be a compound structure that can be fragmented at certain joints. The concept of active templex properties is introduced at this point. The terminology, borrowed from Rinzel’s active and silent phases for bursting mechanisms in excitable systems,³⁵ enables the identification of tipping points within a deterministic templex. A comparison with the pullback approach is provided. Conclusions can be found in the last section.

In algebraic topology, a cell complex $K$ is a layered structure, made up of cells of dimension $k$ with $k∈ N 0$⁠. Thus, a $2$-cell is a polygon, whose faces are $1$-cells or segments, whose faces are $0$-cells or points. In a cell complex $K$ every face of the cell complex is in $K$⁠, and so, every pair of distinct cells in $K$ have disjoint interiors. Such cells are the basic building blocks of homology theory. Homology groups, which are a series of scaled invariants, can be calculated from cell complexes. These invariants encode the information on the existence, on the number, and on the location of holes of different dimensions: they are called $k$-generators or $k$-holes. The zero-order homology group $H 0(K)$ measures the connectivity of the complex, and its rank is associated with the number of connected components; the first homology group $H 1(K)$ identifies non-trivial loops of $1$-cells around the complex (these loops are called $1$-holes); the second homology group $H 2(K)$ accounts for the loops of $2$-cells enclosing empty cavities or $2$-holes. Torsion groups can be computed if the cell complex is uniformly oriented. Above the minimal number of cells (for each dimension) needed to correctly approximate a space,³⁶ the topological invariants (number of connected components, non-trivial loops, enclosed cavities or torsions) are independent of the number or the particular arrangement of cells in the cell complex.¹⁵ 

It has recently been shown that cell complexes, and current tools in algebraic topology, are insufficient to distinguish between attractors with the same number of holes and global torsions.²³ For instance, the spiral³⁷ and funnel³⁸ Rössler attractors are two variants of the Rössler system (with different control parameters) which are not dynamically equivalent in spite of having the same homology groups: they both have one $0$-hole, a single $1$-hole and no $2$-holes. As mentioned in the introduction, traditional topology aims at describing an object in any space, but does not describe the flow if the space is a phase space and the object is an attractor. To add the flow information, the cell complex must be complemented with a digraph. These two elements (the complex and the digraph) make up the twofold mathematical object called templex.

Let us introduce here some details regarding the mathematical construction of the templex from time series. The starting point is a numerical solution of a deterministic dynamical system—a set of time series—leading to a point cloud in phase space. For a deterministic system, the point cloud must be free of false neighbors, so as to study the topology of the solutions themselves and not of their projection onto a lower-dimensional subspace. The point cloud is segmented into $d$-cells where $d$ is the local dimension of the manifold on which the points lie. The set of points assembled in a $d$-cell must also satisfy the condition of being a good approximation to a $d$-dimensional Euclidean set.¹² A cell complex constructed in this manner is called a BraMAH complex—full definitions can be found in Sciamarella and Charó.²⁴ The number of cells forming this complex is neither too low as to risk losing details in the approximation of the underlying structure, nor too high as to needlessly increase the computational cost of calculating homology and torsion groups.

By definition, the dimension of a cell complex is that of its highest dimensional cell(s), and by construction, the dimension of a BraMAH complex coincides with the largest local dimension on which the solution lies. Thus, if the structure supporting the point cloud is a surface, the highest dimensional cells will be $2$-cells and the complex will be said to be of dimension $2$⁠. Note that in a complex of dimension $κ$ there can exist sections of highest local dimension $κ l$ such that $κ l≤κ$⁠. The cells of largest local dimension are assigned a node in the digraph of the templex, and the edges of the digraph provide the flow-compatible connections between the nodes. The templex approach inherits from homology theory the advantage of canceling out extraneous information, and allowing for the algebraic extraction of the fundamental properties of the topological structure of the flow, regardless of the number and distribution of cells.

A templex can be studied homologically using its complex, while its digraph can be analyzed in terms of the distinct cycles than can be established through the junction(s) between different paths. The cells in a junction make up what we call the joining locus. Joining loci can be easily detected scanning the complex for sets of three or more $κ l$-cells sharing a $( κ l−1)$-cell. In other words, a joining locus is the location in a cell complex where there are three or more $κ l$-cells held by a common hinge—if $κ l=2$ this hinge is a $1$-cell (or a singly connected chain of $1$-cells) which may be called joining line. The flow traverses such joining line going from two or more “ingoing” $2$-cells into a unique “outgoing” $2$-cell: a templex fulfilling this condition concerning the cell structure at the joining locus is called a generating templex.

Henri Poincaré was the first to notice how phase portraits are structured around fixed points. Fixed points, however, do not suffice to determine the topology of a chaotic attractor.⁴¹ Let us build the templex for the spiral Rössler attractor from a solution in $R 3$⁠. The BraMAH complex is constructed following the strategy detailed by Sciamarella and Mindlin,¹² where the rationale is dividing the point cloud into subsets that locally approximate $d$-disks to form $d$-cells from them. The local dimension $d$ and the number of points in these $d$-disks is computed using singular value decomposition and linear regressions within a range of maximum and minimum numbers of points that are set beforehand for the entire computation. These $d$-cells form the complex, endowed with a directed graph or digraph every time a $d$-cell is connected with another one by the flow.

An algebraic code provided in https://community.wolfram.com/groups/-/m/t/3079776 is available to compute the templex properties from a generating templex of dimension $2$⁠. The code takes as input a set of $2$-cells spelled in terms of the indexed $0$-cells, as well as the digraph listing of the connections between the $2$-cells. It returns the indexed $k$-holes of the homology groups as well as the torsion groups if there are any, the location of the joining loci, if they exist, and the generatex/stripex set representing the non-equivalent pathways. A generatex is a representative cycle of the original digraph that starts in an outgoing node of a joining locus and comes back to it. Two cycles in the digraph of a templex are said to be equivalent if they share the same ingoing and outgoing nodes. In general, there will be several equivalent cycles: a generatex thus represents a family of equivalent closed paths around a joining locus. To introduce the concept of stripex, we need to define the order of a generatex. An order- $p$ generatex, where $p∈ N$ and $p≥1$⁠, has $p$ ingoing nodes. A generatex is a stripex if $p=1$⁠, and it is composed of $p$ stripexes if $p≠1$⁠. For a complete set of definitions, we refer the reader again to Sciamarella and Charó.²⁴ 

Let $T( R)$ be the templex shown in Fig. 1(a), where $R$ designates the point cloud obtained by integrating Eq. (1), $K( R)$ is the BraMAH cell complex and $D( R)$ is the digraph. The first panel displays a set of four $2$-cells labeled $γ i$⁠, $i=1,2,3,4$⁠. In the cell decomposition mentioned above, we typically obtain around twenty $2$-cells. The number of $2$-cells is purposefully reduced here for simplicity of display. The corresponding digraph is shown in (b), with one node per $2$-cell and directed edges following the flow direction.

Let us now see how the templex properties are extracted from $T(R)$⁠. We start with the homology groups, computed from $K( R)$⁠. The complex has one $0$-generator or $0$-hole since all $0$-cells are connected between them, they are all homologous. Any of the $0$-cells can play the role of generator of the zeroth order homology group. If the representative $0$-cell is labeled $⟨0⟩$⁠, this is written $H 0(K( R))= [ ⟨ 0 ⟩ ]≃ Z$⁠. A complex with one $0$-hole is simply a single-piece complex. $K( R)$ has only one $1$-generator or $1$-hole, shown in green in Fig. 1. All the $1$-holes, i.e., all the non-trivial loops encircling the superficial void in the center of the disk are homologous to the green chain of $1$-cells. We can, therefore, write $H 1(K( R))= [ ⟨ 1 , 3 ⟩ + ⟨ 3 , 6 ⟩ − ⟨ 1 , 6 ⟩ ]≃ Z$⁠. There are no $2$-holes (no enclosed cavities) and so $H 2(K( R))≃0$⁠.

The Rössler attractor is hence a structure with a single $1$-hole contouring the fixed point in the main disk of the attractor. The second fixed point is relatively far from the attractor, and leads to the existence of a fold, not captured by homology groups, which only capture voids. The fold is however related to the existence of a joining locus in the templex. The joining locus is formed here by three $2$-cells: $γ 1$ (outgoing), $γ 3$ (ingoing), and $γ 4$ (ingoing). The three $2$-cells share the joining line $⟨0,1⟩$⁠, the red $1$-cell holding the three joining $2$-cells together, as shown in Fig. 1(a).

A stripex is said to have a local twist if its free edges change their relative positions with respect to the orientation from the center of the attractor to the periphery. In our case study, the second stripex $S 2(R)$ has a local twist. If $S 2(R)$ was isolated from the rest of the templex, it would have a torsion as the Möbius strip, while $S 1(R)$ would have no torsions at all, as in the cylinder or standard strip¹⁸ (see the  Appendix). All the templex properties, namely, the existence and location of (i) a single $1$-hole in the complex, (ii) a single joining line, (iii) two stripexes glued to each other, and (iv) a local twist in the most external of them, provide a complete topological description of the attractor. For a templex-based study of different well-known chaotic attractors, the reader is referred to Charó et al.²³ 

Nonautonomous systems can, in principle, be rewritten in an autonomous form at the price of scaling up in phase space dimension, as observed in Charó et al.⁴³ Topological methods based on knot-theoretical tools, which are limited to three dimensions, are unable to handle higher-dimensional solutions, but the templex approach is a “knotless” method and can therefore deal with systems having more than three dependent variables. When the forcing corresponds to a periodic function, the nonautonomous system can be seen as driven by a one-degree-of-freedom oscillator. When the forcing is periodic, this involves, in principle, two additional phase space dimensions. As the relationship between the two extra variables $F B T$ and $F N S$ is fixed here [see Eq. (4)], only one of the two variables will be required to obtain a false-neighbor-free solution. This is in agreement with false nearest neighbor (FNN) tests.⁴⁴ FNN is an algorithm usually used for estimating the embedding dimension. The FNN method is based on the principle that, in the step from dimension $d$ to dimension $d+1$⁠, one can differentiate between true and false neighbors. A false neighbor is a point in the dataset that appears to be a neighbor only because the points are being projected onto a lower-dimensional space. Here, we apply the FNN technique to the normalized time series and identify the false neighbors in the point cloud. Figure 3 highlights the regions where points that appear to be close together in three dimensions are, in fact, separated (thickening the structure) by the missing dimension. The algorithm used to identify these points, as shown in Fig. 3, is available in Python at https://git.cima.fcen.uba.ar/caterina.mosto/fnn. Notice that the four-dimensional thickening concerns mainly the base disk and a relatively small fragment of the detaching strips.

The point set $N$ corresponds to a two-stripex attractor with a twisted path in the peripheral part of the cell complex, as in the Rössler case. However, the homological description differs: the twisted stripex is detached from the main disk in $N$⁠, whereas this is not true for $R$⁠. This detachment leads to the formation of an additional $1$-hole, which causes $S 1(N)$ and $S 2(N)$ to separate at node $3$⁠, i.e., after the flow passes through $γ 3$⁠. This tearing mechanism is absent in $R$⁠, making the two systems distinct in this aspect. In $R$⁠, $γ 2$ separates $S 1(R)$ and $S 2(R)$ without tearing them apart, keeping their lateral boundaries attached. In this context, homologies can be used to differentiate between tearing and folding. Due to the additional $1$-hole separating the two stripexes, the weakening regime of the AMOC is separated from the strengthening regime by a void of states that the solution cannot access.

In studies using knots, where the topology of the flow is considered independently of the structure serving as support, branched manifolds composed of a normal and a twisted strip (with the strip in a knot-holder being the precursor to the stripex) may be regarded as topologically equivalent to the spiral Rössler attractor. However, it is possible to distinguish different mechanisms leading to this topology by examining the first-return maps to a Poincaré section of the flow.^42,45 Each monotone branch of the map corresponds to a strip. If the two monotone branches are separated by a critical point where the map is differentiable, as in the case of the spiral Rössler attractor, folding occurs. If the map is non-differentiable, tearing occurs (as in a cusp). An example of a two-strip attractor with tearing was proposed by Rössler and Ortoleva,⁴⁶ where tearing can also be identified by the presence of an additional $1$-hole.

Another characteristic of the AMOC three-variable model is the joining line, which does not traverse the entire main disk as in the spiral Rössler case. In $N$⁠, the joining line is a narrow “slit” near the center of the main disk, as shown in the zoom of Fig. 4(a). The Rössler–Ortoleva attractor shares this feature with $N$ (see Fig. 2 in Rössler and Ortoleva⁴⁶), as do systems with a Shilnikov’s saddle-focus homoclinic orbit,^11,12 where the flow is expelled and reinjected near the disk’s center. Another example can be found in the case discussed by Rössler with an artificially composed Lorentzian equation.⁴⁷ 

Figure 5 provides further illustration, offering a schematic diagram of the stripexes in $N$⁠.

In this section, we present the templex analysis of the model corresponding to the original nonautonomous case, where chaotic variability alternates with nearly steady ocean conditions. We work with a single glacial–interglacial cycle yielding point cloud $M∈ R 4$⁠. The additional state variable is $F B T$⁠.

Figure 6 shows the three-dimensional phase portrait of $M$ spanned by state variables $ω$⁠, $F N S$ and $S N S$⁠. Both the glacial and the interglacial intervals are represented, as well as the transitions between them. The point cloud $M$ is dissected into four parts. Let $M= ⋃ i = 1 4 M i$⁠, where $M i∈ R 4$ (⁠ $i=1,…,4$⁠) and $M 1$⁠, $M 2$⁠, $M 3$⁠, and $M 4$ correspond to point sets related to the four sections in the time series. Section $M 2$ corresponds to the glacial interval, characterized by its chaotic dynamics. Sections $M 1$ and $M 3$ correspond to the transitions in and out of the glacial intervals, and $M 4$ corresponds to the interglacial interval.

We will, therefore, construct and analyze $T(M)=(K(M),D(M))$⁠, the templex of the complete four-dimensional point cloud $M$⁠, taking into account the dissection, $K(M)= ⋃ i = 1 4K( M i)$⁠, so that $K( M i)$ is the cell complex obtained from subset $M i$⁠. The location of the holes and stripexes of the dissected templex is accessible through the $0$-cells, whose coordinates in the cell complex are known by construction. Indeed, the BraMAH cell complex provides a link between the topological and geometrical information, allowing for a direct connection between the topological structure and the flow in phase space.

One important caveat regarding the dissection of the templex is that the topological properties of the subtemplexes should be considered with caution. A property such as a $1$-hole in a subcomplex may not correspond to a real void but simply indicate a missing section. Just as with the figure of speech known as synecdoche, the whole should not be mistaken for the part.

The latitudinal shift of the sea ice edge, absent in $N$⁠, has clearly altered the topological structure of the glacial chaotic regime, introducing a third stripex, $S 3( M 2)$⁠, separated from $S 2( M 2)$ by $h 3$⁠. This stripex is “new” in the sense that its cells visit regions of phase space which remained unexplored in $N$⁠. The sketch in Fig. 8 helps recognize the $T( M 2)$ stripexes in the projection spanned by $ω$⁠, $S B T$⁠, and $S N S$⁠.

Two regimes can be identified within the glacial interval of the time series: the first regime, occurring in $[5000yr,12000yr]$ and $[32000yr,35000yr]$⁠, presents an isolated extreme excursion with $ω$ exceeding $0.04$⁠. In contrast, the second regime, operating in $[12000yr,32000yr]$⁠, features twin extremes: a higher peak followed by a lower one, both above $0.04$⁠. Figure 6 (right panel) illustrates how these two regimes correlate with the templex structure. The first regime in the time series of the glacial interval corresponds to the trajectory visiting stripex $S 1( M 2)$ (in red) and $S 3( M 2)$ (in violet), where there is only one excursion far from the main disk, resulting in a single extreme value (with $ω$ above $0.04$⁠). The second regime, occurring in the middle of the glacial interval, corresponds to the trajectory visiting $S 1( M 2)$ (in red) and stripex $S 2( M 2)$ (in blue). This twisted stripex features a double-loop excursion far from the base disk, as shown in Fig. 8. The stripex structure thus explains the twin peaks in the middle regime and the single peak in the regime at the beginning and the end of the glacial interval.

These results can be compared with the time series of the lower-dimensional system discussed in Sec. IV, where only two stripexes were available: the main disk responsible for the low-amplitude, slowly increasing peaks in $ω$—stripex $S 1(N)$—and the double-peak extreme excursions through the twisted stripex $S 2(N)$⁠. The forcing thus introduces the single-peak extreme excursions, which occur when the ESI values are either below $−0.01$ or above $0.05$⁠.

In light of these results, tipping points between different regimes can be interpreted based on the sets of stripexes being visited. As the ESI values increase, the dynamics switches or “tips” from visiting $S 1 , 3( M 2)$ to visiting $S 1 , 2( M 2)$⁠, resulting in qualitative differences in extreme excursions and typical peak values.

This subsection considers the cell structure corresponding to the interglacial time window. The system abandons the glacial section through stripex $S 3( M 2)$⁠, specifically through the cell labeled $γ 13$⁠. The interglacial interval exhibits regular oscillations as the system transitions in and out of the glacial interval. Between $M 3$ and $M 1$⁠, there is $M 4$⁠, with nearly steady conditions.

Topologically, the interglacial interval is represented by a union of subcomplexes, namely, $K(I)=K( M 1)∪K( M 3)∪K( M 4)$⁠. The cell structure involves $2$-cells in $M 1$ and $M 3$⁠, and a $1$-cell in $M 4$⁠. These cells are shown in Fig. 9 using the projection $ω$⁠, $S N S$⁠, and $F B T$⁠, along with the underlying point sets, all colored consistently with Fig. 6.

Transitions between the subsections of $K(I)$ are signaled by topological changes in the cell structure. Let us recall that the local dimension $κ l$ rises from $2$ to $3$ when entering the solid torus—visited by stripex $S 1( M 2)$—and drops from $2$ to $1$ when the orange surface shrinks to the green line covered by the $1$-cell labeled $σ$ in $M 4$⁠.

We are now in a position to analyze the finite-time templex of the glacial and interglacial dynamics from the four-dimensional point cloud, re-establishing the union of the subsets described separately thus far.

From a topological perspective, the trajectory exits the glacial phase due to the existence of a splitting locus in cell $γ 13$ (belonging to $S 3$ in $M 2$⁠), which allows the flow to “tip” to $γ 26$⁠, as shown in Fig. 10. Since the glacial interval is exited via $S 3( M 2)$⁠, the $ω$ time series exhibits single-peak excursions with $ω$ values exceeding $0.04$ before entering the interglacial interval. In other words, it is the regime $S 1 , 3$ in $M 2$⁠, followed by the quasi-periodic oscillations in $M 3$ underlying $S 4$ (cells $γ 13 , 26 … 33$⁠), that suggest the existence of a topological sign of the incoming interglacial interval. On the other hand, the trajectory switches or “tips” from the interglacial interval to the glacial one through cells or nodes $γ 20 … 25$⁠, which lead inexorably to $γ 1$⁠. Here, the oscillations of increasing amplitude in $M 1$⁠, represented by the $2$-cells following the $1$-cell $σ$ in $S 4$⁠, seem to constitute a topological early warning signal (EWS) of the incoming glacial interval. The changes in the local maximal dimension of the cell complex should correspond to an increase in variance, which is consistent with the statistical EWS approach to AMOC collapse.^48,49

To assemble the results of the topological analysis of $M$ into a single picture that illustrates the topological transformations or tipping points occurring within $M$⁠, we introduce the concept of active templex properties. In simple terms, when a templex property, such as a stripex, is not visited for a certain time, it is considered inactive or silent. The same applies to a $1$-hole that ceases to be surrounded by the trajectory, thereby becoming silent or inactive. Hereafter, we will refer to active holes or active stripexes concerning a specific time window and its associated subtemplex. This concept is summarized in the table provided in Fig. 11, using the $ω$ time series as reference.

Point clouds $M 1$⁠, $M 3$⁠, and $M 4$ do not possess any active $1$-holes; the three $1$-holes, $h 1$⁠, $h 2$⁠, and $h 3$⁠, are only active for cloud $M 2$⁠. The stripexes $S1$⁠, $S2$⁠, and $S3$ are active in $M 2$⁠. We can discriminate between the two regimes within the glacial interval in terms of active stripexes: $S 1 , 3$ at the beginning and end of the $M 2$ interval, and $S 1 , 2$ in the middle. During the interglacial interval, which covers $I$ (⁠ $M 1$⁠, $M 3$⁠, and $M 4$⁠), the only active stripex is $S4$⁠. It is noteworthy that $S 4$ remains bonded to the glacial part of the templex. Note also that the part of $S4⊂I$ is composed of cells with progressively lower $κ l$ (from $3$-cells to $2$-cells and $1$-cells). Subsequently, the local dimension of the subcomplexes increases again to reach a new chaotic burst.

It is important to emphasize that our templex is constructed from a single glacial–interglacial cycle, meaning that the non-bonded part of stripex $S 4$ is visited only once. The templex structure remains consistent across all the cycles we have considered, provided that each cycle is analyzed separately. However, if the solution used to build the templex encompasses multiple interglacial-glacial cycles, the resulting structure will display features not studied here. A multi-cycle templex for the non-autonomous case is left for future work.

Nonautonomous systems are often studied using the concept of pullback attraction, which generalizes global attractors for nonautonomous systems.^31,33 Pullback attractors describe the long-term behavior of a system as it is “pulled back” in time, offering insight into how trajectories evolve in varying environments.³⁰ 

The nonautonomous idealized model of the AMOC has been analyzed by its authors through the lens of pullback attraction. Figure 12 presents four snapshots in the pullback sense,¹⁰ juxtaposed with different subsets of $K(M)$⁠. In gray, we plot the points representing the highest density regions, along with the subcomplexes in the $(ω, S N S, F B T)$ projection of the four-dimensional phase space.

The four snapshots of the pullback attractor are obtained for different values of $( t 0,t)$⁠, where $t 0$ is the initialization time and $t$ is the diagnostic time. Computations start $5$ kyr before $t 0$⁠. In the first frame (A), the pullback point set lands along the cells in $S 1 , 2(M)$ in accordance with the value of $t 0=−175$ kyr. In the second snapshot (B), the highest density regions visit the cells in $S 1(M)$ and the periphery of $S 3(M)$⁠, as expected for $t 0=−150$ kyr. In (C), the pullback attractor is concentrated in the interglacial part of $S 4(M)$ corresponding to $M 3 , 4⊂I$⁠, in agreement with $t 0=−100$ kyr. In (D), with $t 0=−125$ kyr, the central parts of $S 1 , 3$ take the lead.

The time-varying point sets of the evolving pullback attractor can hence be translated into the varying active properties of our templex. $T(M)$ supports the highest density regions (the most visited parts of the three-variable phase space) obtained from a very large number of numerical experiments (two million), each with different initial conditions. Let us note that calculating the pullback attractor for the idealized nonautonomous model requires as much numerical effort as simulating several decades with a state-of-the-art ocean general circulation model. This result suggests that the most populated stripexes of the idealized model are contained within our finite-time templex, which was constructed from a single glacial–interglacial cycle and a single numerical experiment (one initial condition). Taking into account the juxtaposition of the snapshots upon $T(M)$⁠, we can conclude that the not-yet-constructed multi-cycle templex may have extra stripexes, but these should correspond to seldom-visited ones. Using first-return maps to the Poincaré section of the flow should help confirm this.

This work shows that tipping points in a deterministic dynamical system’s behavior can be effectively analyzed through changes in the topological active properties of its phase space flow, which can be accurately detected using the templex approach. The methodology is applied to a recently proposed idealized model of the Atlantic Meridional Overturning Circulation in autonomous and nonautonomous settings. In the autonomous case, the model mimics the chaotic variability of the glacial intervals using three variables: the anomalies in the overturning rate, the bottom–top, and the north–south salinity gradient. In the nonautonomous case, the model takes into account the glacial–interglacial variations in the sea ice cover. The fourth variable is the Fourier projection of the bottom–top surface-salt flux.

The templex is a mathematical object based on the definition of a cell complex. Composed of cells of different dimensions that are glued together, a cell complex allows for the computation of topological properties that do not depend on the specific number or arrangement of cells. These properties may include torsions and voids, referred to as holes, which are quantified through homology groups derived from the cell structure. These groups reveal the presence of connected components (at level $0$⁠), non-trivial loops (at level $1$⁠), empty cavities (at level $2$⁠) and higher-dimensional voids, thus providing a layered description of the topology of the structure. In chaos topology, this structure often takes the form of a manifold with branches. A special type of cell complex called BraMAH can be constructed to describe features such as joining (splitting) loci, where branches merge (separate). The main difference between a templex and a cell complex is that the templex incorporates the “arrows” of a flow on the cell structure. A templex is therefore a BraMAH cell complex endowed with a directed graph, which carries information about how the cells are connected by the flow. As a result, templex properties not only describe the topology of a structured point cloud in a phase space of arbitrary dimensions, but also the topologically distinct cyclic paths allowed by the flow around that structure. The structure itself is characterized by holes, torsions, and joining or splitting loci, while the flow is represented by the generatex/stripex set, which corresponds to a set of non-equivalent pathways around the joining loci.

In the autonomous case, the chaotic oscillations represent the irregular strengthening and weakening of the AMOC associated with the Dansgaard–Oeschger events. The rapid switches between these on- and off-phases correspond to a branched manifold with two $1$-holes and two stripexes, one of which is twisted. One of the $1$-holes is of the focus type: the empty space is circled by the flow. The other $1$-hole is associated with a tearing of the flow, which creates an empty gap between the two stripexes. The circulation weakening is directly related to the twisted stripex.

In the nonautonomous scenario, when the driving force is active, the dynamics can be embedded in four dimensions without false neighbors. A third stripex and a third $1$-hole appear in this context, located in regions of phase space that remained unexplored in the projection spanned by the three variables of the autonomous model. Stripexes have parts in common that keep them bonded, but the interesting part of the dynamics is given by their non-overlapping parts. Two regimes can be distinguished based on how the stripexes are visited: one of these, associated with the two stripexes already present in the autonomous case, exhibits a double peak during the weakening phase, while the other shows a single peak.

The term active templex properties is introduced in the nonautonomous case to characterize the topological changes induced by external forcing. These properties, which may or may not be observed in a subtemplex during a given time interval, are referred to as active or silent. For instance, in the nonautonomous model, the three $1$-holes are active only during the chaotic burst and remain silent during the interglacial interval. Another key templex property is the local dimension of the cell complexes, which marks topological jumps in the time series. The local dimension of the cell structure is $1$ during nearly steady intervals, $2$ during quasi-periodic transitions, and rises to $3$ during chaotic bursts. Tipping points can therefore be identified as the cells or nodes where the trajectory switches from one set of active properties to another. Early warning signals can also be formulated using the topological terms of this approach, as both the glacial and interglacial intervals are preceded by specific topological indicators. In contrast, there are no tipping points in the autonomous case, as all templex properties remain continuously active.

An important comment is necessary regarding how the templexes in this work were constructed. They are all finite-time objects in the sense that they are built from finite time series, but there are significant differences in terms of the number of cycles in the time series considered for building the three-variable and the four-variable templexes. The two stripexes of the autonomous case correspond to approximately twenty strengthening and weakening cycles; thus, the two stripexes are visited around twenty times. In contrast, the interglacial interval, which is only reproduced by the nonautonomous model, is traveled only once in the point cloud used for the topological study. As far as we could test, the nonautonomous templex properties seem to remain unchanged for different glacial-interglacial cycles, as long as each cycle is investigated on its own. Of course, a nonautonomous multi-cycle templex structure is expected to present more complex interglacial properties related to the repeated visits of different interglacial intervals. This aspect has been left out of the scope of this study but will be considered in a separate work.

When the nonautonomous templex properties are compared with the analysis of the nonautonomous model conducted using the pullback approach, the templex from our study serves as support for the points corresponding to high probability densities of the pullback attractor. The most populated regions of snapshots computed for different initializations and diagnostic times from two million numerical experiments¹⁰ land within our nonautonomous templex. This suggests that our topological results are robust enough to describe the most visited stripexes of the idealized nonautonomous model, despite being constructed from a single glacial-interglacial time window based on only one numerical integration.

In conclusion, the templex approach has proven useful for detecting early warning signals and tipping points of AMOC collapse,^48,49 formulated in terms of topological active properties. This opens up the possibility of using the templex method to detect AMOC tipping points under different forcings. Such analyses could be conducted through the examination of state-of-the-art climate model simulations from both the historical period and future projections provided by the Coupled Model Intercomparison Project.⁵⁰ This will also be the focus of future work.

This work has received funding from the ANR project TeMPlex ANR-23-CE56-0002 (D.S.). Insightful comments and suggestions received during the review process are sincerely appreciated. G.D.C. gratefully acknowledges her postdoctoral scholarship at LSCE within the ANR project. C.M. gratefully acknowledges her doctoral scholarship from CONICET.

The authors have no conflicts to disclose.

Caterina Mosto: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (lead); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Gisela D. Charó: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Florian Sévellec: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Pierre Tandeo: Data curation (equal); Validation (equal); Writing – original draft (equal). Juan Ruiz: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Denisse Sciamarella: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

In algebraic topology, planar diagrams are used as a visual and combinatorial tool to represent cell complexes corresponding to higher-dimensional objects in a flat form, with gluing instructions indicated by repeating the labels of points or $0$-cells to be joined, just as one would when assembling a paper model.

Figure 13 shows a planar diagram for a cell complex with a single $2$-cell labeled $α$ and $0$-cells labeled $P$⁠, $Q$⁠, $U$⁠, and $V$⁠. By taking a rectangular piece of paper and gluing the left and right edges together, making the repeated labels coincide, a cylinder is formed. The cylinder can then be deformed into a standard band: a disk with a hole in the middle.

Gluing instructions are crucial: if the repeated labels are positioned differently on the rectangular paper, a Möbius band can be formed instead of the standard band, as shown in Fig. 14. In both cases, the associated cell complex has a single $2$-cell $α$⁠, but the $1$-cells and $0$-cells are assembled differently.

We can, of course, construct planar diagrams for the cell complexes throughout our work. Let us start with point set $N$⁠. The planar diagram shown in Fig. 15 is drawn using the same labels as those in the schematic diagram of Fig. 5. The simplified cell complex $K ′(N)$ has two $2$-cells labeled $γ 1 ′$ and $γ 2 ′$⁠, corresponding to $γ 1 ⋯ 8$ and $γ 9 ⋯ 29$ in Fig. 4. Thick lines are used to indicate the joining locus in the cell complex. Five $0$-cells are being used: $P$⁠, $Q$⁠, $R$⁠, $S$⁠, and $W$⁠.

Let us consider the point set in Fig. 7 along with its schematic diagram in Fig. 8. The corresponding planar diagram is shown in Fig. 16(a). The simplified cell complex is labeled $K ′( M 2)$⁠. Note that the $2$-cells in the main disk actually represent $3$-cells that have been flattened.

A planar diagram of a simplified cell complex $K ′(M)$ is given in Fig. 16(b). The rectangle at the top corresponds to the main disk. The labels correspond to the schematic diagrams in Fig. 8. The $1$-cell $σ$ in Fig. 9 is shown here with its label in $K(I)$⁠.