The topology of a chaotic attractor in the Kuramoto–Sivashinsky equation

The importance of unstable periodic orbits (UPOs) to the understanding of chaotic dynamics has long been appreciated. It is believed, and in some cases proven, that UPOs are dense within strange attractors. Modern computational methods enable the detection and convergence of large numbers of such UPOs in high-dimensional discretizations of dynamical systems governed by partial differential equations (PDEs), including the Navier–Stokes equations from fluid dynamics.^1–4 Especially in the context of fluid turbulence, UPOs carry the hope that we may be able to describe statistical properties of the turbulent flow as weighted ensemble averages over these UPOs. This envisioned form of “periodic orbit theory” would transfer well-studied ergodic theory from low-dimensional chaotic ordinary differential equation (ODE) systems^5–7 to PDEs.

One alternative approach to use periodic orbits for understanding the behavior of chaotic attractors is through ideas from topology. In three dimensions, UPOs, which are closed loops in state space, form knots, and pairs of UPOs form links. Such knots and links possess topological invariants, which can be used to analyze the topological structure of the underlying chaotic attractor. In particular, the Birman–Williams theorem^8,9 shows that a chaotic attractor in a three-dimensional system can be identified with a two-dimensional branched manifold, a so-called template, on which all periodic orbits lie. The structure of this branched manifold—the number of branches and how they link together—can be deduced from the linking numbers of the periodic orbits. This formalizes the observation that chaotic attractors in canonical examples such as the Lorenz and Rössler systems seem to be very thin, almost two-dimensional structures, and the templates found for these systems match the structures visible in a long chaotic trajectory in the usual coordinates. The method of determining templates for all periodic orbits has been successfully applied to many three-dimensional ODE systems, including ones for which the topology is not obvious from observing trajectories.^8,10,11

In four and higher dimensions, all UPOs are, however, equivalent (i.e., ambient isotopic) to the unknot, and thus the specific methods discussed here cannot be directly applied. The solution spaces of PDEs are infinite dimensional. In practice, when solving PDE dynamical systems numerically, this infinite dimensional space is discretized to be finite-dimensional, but the number of dimensions is still on the order of at least tens, if not much higher, and so the ideas of Birman and Williams are not directly applicable (although more modern approaches exist which can be applied to dimensions higher than three, as discussed in Sec. V). Despite the high-dimensional solution spaces, the long-time dynamics of many systems often collapses onto an attracting set of much lower dimension.

The dissipative PDE systems that arise in fluid dynamics are typically associated with a parameter, such as the Rayleigh number or Reynolds number, whose increase is correlated with an increase in the complexity of the dynamics. When such a parameter is sufficiently high the system is fully turbulent, as characterized by a large number of positive Lyapunov exponents. In this case, the global attractor is a relatively high- (but still finite-) dimensional fractal structure in state space. In contrast, when the parameter is sufficiently small, the system collapses onto a single attracting fixed point. Between these regimes, it is often possible to find chaotic dynamics with only one positive Lyapunov exponent, which indicates a strange attractor with less than three dimensions (in the absence of symmetries), embedded within the infinite-dimensional state space of the PDE. If it were possible to construct an embedding of this low-dimensional attractor into three-dimensional space, this would preserve the topological properties, and we would be able to use the ideas of Birman and Williams to find a template for the system.

Since the topological properties of dynamical systems are independent of the choice of coordinates used for parameterizing state space, they provide a way to compare different representations of the same or similar systems. This is particularly relevant to PDE systems, in which the dynamics are necessarily discretized in order to be solved on a computer. Different choices of discretization (for example, finite differences or Galerkin methods) may lead to very different state space geometries, but the topology must be equivalent if the discretizations are accurate. Furthermore, it is often desirable to find a “reduced-order-model,” a small set of ordinary differential equations approximating the dynamics of a PDE system, either to increase computational efficiency or to aid the interpretation of physical mechanisms underlying the dynamics. One way to validate the accuracy of such a model would be to check that the model preserves the topology of the dynamics of the full system. However, in order for this we need a method to distill the topology of a chaotic PDE, which does not currently exist.

Previous authors have attempted to find reduced-order models with the same dynamics as infinite dimensional systems, for example, by using global modeling to model the dynamics with a three-dimensional ODE, from which a template can be derived.¹² Our approach is instead to consider accurate periodic trajectories of the full PDE system, and then to embed these trajectories in a three-dimensional reduced space.

Reducing the state space of a high-dimensional dynamical system to an accurate but approximate low-dimensional representation is an important and active area of research. Traditional methods, such as proper orthogonal decomposition¹³ (POD) and dynamic mode decomposition,¹⁴ use numerical linear algebra tools to determine the most important “modes” from a timeseries. This is often successful, although these methods have a fundamental limitation when the dynamics and state-space geometry are strongly nonlinear. For these reasons, researchers have applied the more flexible and nonlinear methods of deep learning to dimensionality reduction. In this paper, we use POD as well as the more modern approach of an autoencoder, a pair of deep neural networks which are trained to find an invertible (as close as possible) mapping from the full state space to a low-dimensional “latent” space and back. Due to their nonlinearity, autoencoders have been applied as a dimensionality reduction tool to a wide variety of problems in nonlinear dynamics, including chaotic systems. Questions on the interpretability of the latent space, such as the minimum number of latent dimensions and the physical meaningfulness of latent coordinates, have also been studied. Linot and Graham¹⁵ investigate how the errors of these networks vary depending on how close the latent dimension is to the manifold dimension of the chaotic attractor. Page et al.¹⁶ showed that by applying deep convolutional autoencoders to 2D Navier–Stokes, they are able to identify meaningful low-dimensional representation of two-dimensional turbulence. The apparent success of autoencoders at finding a low-dimensional representations for such complicated systems indicates the potential that the methods presented in this paper may be applicable to higher-dimensional systems including chaotic fluid flows.

The idea that autoencoders can in some cases preserve the topology of a chaotic attractor is supported by the studies of Uribarri and Mindlin¹⁷ and Fainstein et al.,¹⁸ who applied autoencoders to timeseries from synthetic and experimental data. Unlike in the present work, the aim there was to determine whether the known topological properties are preserved so that the autoencoders give an embedding. Indeed, both of these papers showed that this is not always the case, and so autoencoders should be used with caution.

We consider the Kuramoto–Sivashinsky equation (KSE), a dissipative PDE in one spatial dimension, which is frequently used as a model for the Navier–Stokes equations and other more complicated PDE systems. In this system, the domain size $L$ dictates the complexity of the dynamics. When constrained to its antisymmetric subspace, stable chaotic dynamics is first observed at $L≈18$⁠, whose existence was proven by Wilczak and Zgliczyński.¹⁹ This chaotic attractor is mapped into a three-dimensional space using both an autoencoder and also a traditional POD analysis. Using these two mappings to $R 3$⁠, we can compute linking numbers for the known periodic orbits and construct a template for them. If both mappings preserve the topology, as hypothesized, the template should be equivalent for the POD and the autoencoder-based mapping.

The topological properties of chaos within the KSE (in a very different parameter regime) have been studied by Siminos²⁰ using an entirely different and more general technique. In it, the dimensionality reduction technique of manifold learning was applied to Poincaré return maps, essentially learning the topology of a discrete-time dynamical system. By examining the linking of periodic orbits, our method gains a further understanding of the topology as a continuous-time system.

The remainder of this paper is laid out as follows. In Sec. II, the PDE system we study is described. In Sec. III, we explain the dimensionality reduction methods and their application to the system and its periodic orbits. In Sec. IV A, we give background on the topological approach to periodic orbits, and in Sec. IV B the periodic orbits are used to find a template for the system. Finally, concluding remarks are given in Sec. V.

Equation (1) is invariant under the transformation $x↦−x$⁠, $u↦−u$ and therefore admits an invariant subspace of solutions for which $u(−x)=−u(x)$⁠. Following previous studies,^26–28 we consider the dynamics constrained to this (usually unstable) subspace. This has the advantage of removing the continuous symmetry $x↦x+Δ$⁠, which otherwise would increase the dimension of any $x$-varying attractor by 1, significantly complicating topological analysis. In this subspace, the first known chaotic attractor occurs at $L≈18.05$⁠. Wilczak and Zgliczyński¹⁹ give a computer assisted proof for the existence of this attractor at $L=2π/ 0.1212=18.048…$⁠, which is the parameter value we use. Indeed, they proved the existence of a symbolic dynamics for the system, which in principle allows the enumeration and computation of a countable infinity of UPOs. We found, empirically, that this attractor has a very small basin of attraction: it is necessary to use an initial condition very close to the states depicted in Wilczak and Zgliczyński;¹⁹ otherwise, the system rapidly collapses onto a fixed point.

To have any hope of applying methods based on knot theory to such a system, the fractal dimension of the attracting chaotic set must be less than three. An upper bound is given by the Lyapunov dimension,²⁹ which we calculate to be approximately $2.227$ based on the four leading exponents in the Lyapunov spectrum being $0.0149$⁠, $0$⁠, $−0.0658$⁠, and $0.0989$⁠. Note that no positive Lyapunov exponents were found at this parameter value by Edson et al.²⁵ (in that work corresponding to $L≈9$ with odd-periodic boundary conditions) but this is unsurprising given the localized nature of the chaotic attractor and the very limited range of $L$ for which it exists.

In order to be able to embed the chaotic attractor in a space of dimension 3, it is necessary but not sufficient for the fractal dimension to be less than three. Consider for example the Klein bottle, whose dimension is only two but cannot be embedded within a Euclidean space of dimension less than four. Experiments at larger values of $L$⁠, for example, the well-studied^27, $L≈39$⁠, and the parameter values studied by Dong²⁸ and Siminos²⁰ of $L≈36$⁠, found attractors whose fractal dimension was apparently less than three, but for which we could only produce a non-injective immersion into a 3D space. An injective embedding is necessary for our topological approach to be valid.

We discretize the spatial dimension of the anti-symmetric system with $N x=32$ points, turning the function $u$ into a 32-dimensional vector $u$⁠. We observe a satisfactory drop in the amplitudes of the Fourier modes for this discretization, suggesting a well-resolved discretization of the PDE. For comparison, Lan and Cvitanović³⁰ also use 32 modes in the $L=38.5$ system to find periodic orbits, which exhibits stronger chaos²⁵ and smaller-scale spatial features than our system. To calculate the POD modes in Sec. III A and train the autoencoder in Sec. III B, we generate a long direct numerical simulation (DNS) of the system using the ETDRK4 scheme,³¹ with $t max=20000,dt=0.1$⁠. Due to the small basin of attraction of the chaotic attractor, we pick a point that we know for certain lies on the attractor as initial condition for the DNS.

A section of this timeseries is shown in Fig. 1, in the typical manner of plotting space along the vertical axis and using a color scale to show the value of $u$ as time varies. Unlike the usual projection of trajectories in the Lorenz system, for example, no topology for the dynamics is evident in this figure. Indeed, for this value of $L$⁠, it is not even obvious that the dynamics is chaotic.

Using a long timeseries from a simulation of the KSE as described in Sec. II A, we find the three leading eigenvalues to be $λ 1≈19.49, λ 2≈1.92$⁠, and $λ 3≈1.10$⁠, with the corresponding modes plotted in Fig. 2. The projection of the attractor defined by the first three POD modes is plotted in Fig. 3 (together with one UPO to emphasize the shape). We obtain the POD latent coordinates $( ξ 1, ξ 2, ξ 3)$ by projecting onto the POD modes. Formally, this is done by evaluating the integral $ξ i(t)= ∫ − L / 2 L / 2u(x,t) ϕ i(x)dx$⁠. However, since the spatial dimension is discretized, we obtain the latent POD coordinates from the vector dot product $ξ i(t)=u(x,t)⋅ ϕ i$ Note that the terminology of “latent space” for the expansion into POD modes is chosen in analogy to the alternative autoencoder-based dimensionality reduction approach discussed below.

The three leading POD modes capture 99.3% of the variance, which is sufficiently high that we hope the map from the full space to the coefficients of these three modes preserves the topology and, thus, allows consistent application of the methods of Sec. IV B.

While POD projections capture a great amount of information, they are known to generalize less well to highly nonlinear systems and are outperformed by autoencoders.³⁴ Autoencoders are neural networks that consist of two parameterized functions, namely, the encoder $E: R N x→ R N h$⁠, where $N h$ is the dimension of the latent space, and the decoder $D: R N h→ R N x$⁠. Typically, $N h≪ N x$⁠. The parameters are trained so that the autoencoder approximates the identity. That is, for an input vector $u$⁠, we have $( D ° E )(u)≈u$⁠. Thus, the encoder reduces the dimension of the input, while the decoder does the inverse operation and increases it again. The structure of our autoencoder was chosen so that $E$ and $D$ are smooth functions. The goal is to choose $N h$ as small as possible while still preserving the topology of the system. If it is possible to use $N h=3$ such that $E$ is an embedding, invertible for points on the attractor, the methods of Sec. IV A can be applied.

Since we define the latent dimension to be $N h=3$⁠, this gives us the desired three-dimensional representation of the KSE, in which we would like to test the methods described in Sec. IV A. Taking datapoints from a long physical timeseries and applying the encoder, we obtain a 3D latent timeseries. Figure 6 shows a 3D plot of the image of the attractor in latent space, illustrating its band-like shape, and 2D projections.

In Secs. III A and III B, we described how proper orthogonal decomposition and autoencoders allow us to obtain three-dimensional approximate coordinates for the KSE, denoted by $( ξ 1, ξ 2, ξ 3)$ and $( h 1, h 2, h 3)$⁠, respectively. In this section, for each of the state space reduction techniques, we will define a Poincaré section and a return map. These will help us with identifying UPOs and allow us to obtain a set of names for them.

In the POD latent coordinates $( ξ 1, ξ 2, ξ 3)$⁠, we define a Poincaré section $P POD$ given by the plane $ξ 1=1.3$ and flow direction $∂ t ξ 1<0$⁠. The value $ξ 1=1.3$ is chosen so that the flow is approximately orthogonal to the Poincaré section. This defines a unique crossing section along the band (see Fig. 7) and gives us a sequence ${ ξ 2 ( i ) }$ of the $ξ 2$ coordinate of the crossings through $P POD$⁠. The return map $f 1$ is defined by $f 1: ξ 2 ( i )↦ ξ 2 ( i + 1 )$⁠. A plot of such subsequent crossings is shown in the top left panel of Fig. 8. The points approximately fall on a smooth unimodal curve, which we fit with a polynomial of order three. This simplified first return map appears to only have one intersection with the identity $ξ 2 ( i )= ξ 2 ( i + 1 )$⁠, suggesting that there is only one periodic orbit with $p=1$⁠, where $p$ is the number of intersections of an orbit with $P POD$⁠. By looking at further iterations/compositions of the return map (see other panels of Fig. 8), we also obtain an indication that there are 1, 0, and 1 periodic orbits with $p=$ 2, 3, and 4, respectively (accounting for crossings that correspond to periodic orbits with smaller $p$—e.g., the orbit with $p=1$ will also account for one intersection in all subsequent graphs). Finally, we split the band domain $I=[−3.229,−2.952]$ at the minimum value of the first return map, giving approximately $I 1=[−3.229,−3.055]$ and $I 2=[−3.055,−2.952]$⁠. From this, we obtain an indication of the qualitative dynamics: points starting in $I 1$ are mapped to the full interval $I$⁠, while points starting in $I 2$ are mapped to a subset of $I 1$⁠. This is shown in Fig. 9.

We proceed in the same way with the autoencoder’s latent coordinates $( h 1, h 2, h 3)$ . We define a Poincaré section $P AE$ given by the plane $h 2=3/4$ and flow direction $∂ t h 2>0$⁠, so that the flow is approximately orthogonal to it. Figure 10 shows that this defines a unique crossing section along the band, giving us a sequence ${ h 1 ( i ) }$⁠. The return map $g 1$ is defined by $g 1: h 1 ( i )↦ h 1 ( i + 1 )$ (see the top left panel of Fig. 11). We again fit the resulting smooth curve, with a polynomial of order three. We make the same observations as for the POD return map: $g 1$ seems to have one intersection with the identity $h 1 ( i )= h 1 ( i + 1 )$⁠, implying one periodic orbit with $p=1$⁠. Further iterations/compositions of the return map (see Fig. 11) indicate that there are 1, 0, and 1 periodic orbits with $p=$ 2, 3, and 4, respectively. We again split the band domain $J=[0.198,0.247]$ at the minimum value of $g 1$⁠, giving approximately $J 1=[0.198,0.23]$ and $J 2=[0.23,0.247]$⁠. We observe similar qualitative dynamics to the POD return map: points starting in $J 1$ are mapped to $J$⁠, while points starting in $J 2$ are mapped to a subset of $J 1$ (see Fig. 12).

We find periodic orbits (such as the one plotted in Figs. 3 and 6) by generating guesses using the methodology from Beck et al.³⁵ and converging them using the algorithms from Azimi et al.³⁷ This is similar to the variational method used by Dong²⁸ to find UPOs and symbolic dynamics at a higher $L$ value in the KSE. We obtain the guesses by defining arbitrary closed curves (loops) within the latent space, and decoding them to the physical space. This gives us a timeseries that is already time-periodic and lies close to the chaotic attractor; however, it does not necessarily satisfy the PDE. The convergence algorithm then deforms the loop in an attempt to turn it into a solution, by minimizing a cost function $J$ that quantifies the misalignment between the tangent vectors and the flow vectors at each point of the loop. Thus, a root of the cost function corresponds to a closed curve in the state space that satisfies the flow equations everywhere and, hence, is a periodic orbit.

Following this methodology, we find many periodic orbits of the system, not all of which necessarily lie on the chaotic attractor that we are considering in this paper. To characterize the attractor, we need to ensure that only those periodic orbits that lie on the attractor are considered in the analysis. Since the latent space is three-dimensional, we can plot the attractor with the found periodic orbits. Visual inspection allows us to discard those orbits that clearly do not lie on the attractor. In our analysis, we thus only consider the periodic orbits that appear to lie on the chaotic attractor in the latent space (of which we find 15 with $p<12$⁠). A list of periods $T$ and intersections $p$ with $P POD$ and $P AE$ (as defined in Sec. III C) is given in Table I.

Comparing the extracted list of 15 periodic orbits with the return map and its iterates in Fig. 11 suggests that we have successfully identified all relevant periodic orbits with $p<7$⁠. For example, we initially find periodic orbits with $p=3$⁠, but upon visual inspection, we discarded all of those as not lying on the attractor. This agrees with our return maps. Since the only intersection of the third iterate with the identity (in the bottom left panel of Fig. 11) corresponds to the $p=1$ orbit, the iterated return map indeed suggests that no orbit with $p=3$ on the attractor exists. This is consistent with the symbolic dynamics derived from the approximate return map as shown in Fig. 12, since period 3 symbolic sequences are not admissible.³⁸ Recall that the existence of a period 3 orbit implies chaos, but the converse is not true. Likewise, all identified orbits exactly correspond to those expected from the return maps. For $p≥7$⁠, it is difficult to be certain about the number of periodic orbits by purely looking at the return maps, and so we do not claim that we have found all of them with certainty. All orbits that appear to lie on the attractor were included in our analysis, and none were discarded purely on the basis of a return map argument. For the $15$ orbits, we arbitrarily chose names based on $p$ and a running index, as indicated in Table I.

The Birman–Wiliams theorem also states that periodic orbits of the strange attractor are in one-to-one correspondence with the periodic orbits on the template, with at most two extraneous orbits.^9,40 Moreover, the template has the property that all periodic orbits can be projected onto it without altering any of their topological invariants, crucially including their (self-)linking numbers. Certain preconditions of the Birman–Williams theorem can be relaxed; see Chapter 13 in Gilmore and Lefranc⁴¹ for more details. Given a higher dimensional flow whose Lyapunov dimension is less than 3, one can first project the flow along its stable directions to an intermediate manifold where we can apply the Birman–Williams theorem. In the following, we attempt to construct the template of the chaotic attractor from the identified set of periodic orbits and thereby identify the structure of the spectrum of periodic orbits on the attractor of the PDE.

A branched manifold associated with a template of $n$ branches can be given an algebraic description consisting of a small number of integers: torsion and layering numbers.⁴¹ The torsion terms are encapsulated in a $n×n$ symmetric matrix of integers describing how the branches are knotted together. The entry $t i , j$ is the signed number of crossings between branch $i$ and branch $j$⁠, and $t i , i$ is the self-torsion of branch $i$⁠, i.e., signed number of crossings between the two edges of the branch. In addition, the layering of the branches can be encapsulated in a vector of integers, consisting of entries $l i j$ for $i<j$⁠. The entry $l i , j$ for $i<j$ is $1$ if branch $j$ is closer to the reader than branch $i$ and $−1$ otherwise. While this algebraic description is not topologically invariant (it depends on the choice of projection of the template in 2D), it determines the arrangement of the periodic orbits, which is the topological invariant of interest in this work.

These integers describing the template can themselves be determined by computing one of two topological invariants of the periodic orbits: their (self-)linking numbers or relative rotation rates.

In summary, a template can be determined from an input set of periodic orbits via the following algorithm:⁴¹ 

1. Compute the linking numbers of the input set of orbits.

2. If possible, determine the number of branches in the template from the return map. The number of branches is given by the stretching that occurs during one period. It equals the number of monotonic branches of the first return map.⁴¹ 

3. Iterate over the possible symbolic sequences for the orbits (and/or the number of branches) and solve the system of Eq. (7). If there is a unique valid solution, one can include additional orbits to check the validity of a candidate. If, on the other hand, multiple valid names and template candidates are found, the problem is undetermined and more orbits must be included.

Hence, the algorithm returns valid set(s) of possible names together with a template. The procedure is described precisely in Appendix A in Gilmore and Lefranc.⁴¹ An overview of the topological analysis program is given in Lefranc.⁴⁵ The previous algorithm does not immediately extend to flows in dimensions higher than three. However, in the considered case of the KSE, it appears that the attractor can be embedded in $R 3$ where the computations of linking numbers can be carried out.

Strange attractors can be classified by the topological organization of periodic orbits.⁴⁶ If the dynamics can be embedded in three dimensions, then the periodic orbits form knots and links and their topological organization is described by their linking numbers. Then, by the Birman–Williams theorem, one can project the flow onto a template while preserving the organization of the periodic orbits. A template not only supports all the periodic orbits and describes their topological organization but also offers a natural way for enumerating these periodic orbits. For $L≈18.05$⁠, the Lyapunov dimension associated with the KSE system is lower than three. We use two distinct tools: proper orthogonal decomposition and an autoencoder neural network, described earlier, to embed the dynamic in 3D latent spaces. Then, we compare the templates resulting from each method.

We use the POD described earlier to embed the dynamics of the Kuramoto–Sivashinsky PDE from the original phase space into an intermediary manifold having dimension three. Since the latent space is three-dimensional, the Birman–Williams theorem can be applied and the procedure described earlier can be carried out to determine the template for the KSE.

Note that after projecting the periodic orbits found in Sec. III D onto the POD space, we compute their linking numbers using the algorithm described by Qu and James.⁴⁷ The linking numbers are given in Table II. We did not calculate self-linking numbers for the periodic orbits, as this requires the choice of a framing, and it is not clear how to do this consistently given the dimensionality reduction. It is not necessary to compute self-linking numbers to determine a template, but without them it requires more periodic orbits to find a unique template.⁴⁸ 

Some of the UPOs considered have very long periods. We found that numerically calculating the linking numbers of these requires discretization of each loop using a very high number of points, given the necessity of accurately counting the crossings of periodic orbits in very close proximity. To resolve the most challenging cases, it was necessary to discretize the orbits with $2 13$ points.

A template provides a natural way to enumerate the periodic orbits: to each periodic orbit, one assigns a word whose letters represent the successive branches traversed. A collection of periodic orbits with unique symbolic sequences are a key element of a method for constructing a generating partition by interpolating the encoding map.^43,50

During the template identification, we obtained several sets of possible names so, according to the algorithm described in Sec. IV A, we should add more orbits in the input set. However, in our particular case of the KSE, we are able to determine the symbolic sequences of the input set of periodic orbits using the POD return map as shown in Fig. 9. The set of symbolic sequences for the input set of orbits is presented in Table I. The choice of the starting location of a symbolic sequence for a given orbit is arbitrary; we follow the convention given by Hao and Zheng,³⁸ i.e., choosing the sequence whose initial string has even parity.

We proceed similarly with the autoencoder method. Again, as shown in Fig. 12, from the autoencoder return map we obtain the same set of symbolic words for the input set of periodic orbits. The linking numbers obtained from the symbolic sequences and the template identification (9) are consistent with the linking numbers obtained from the embedded periodic orbits in the autoencoder latent space.

Thus, consistent symbolic sequences for low-period UPOs are derived from two different dimensionality reduction methods and the symbolic sequences are validated by the template identification.

In this paper, we have found a candidate template for the topology of a chaotic attractor in the Kuramoto–Sivashinsky equation. To our knowledge, this is the first attempt to do this in a PDE. While it is possible that our template would be shown to be invalid when other periodic orbits are included, it was found using only five UPOs and validated with a further ten, so we are confident it gives an accurate representation of the topology of the dynamics. Furthermore, equivalent topologies were found using two different dimensionality-reduction techniques: proper orthogonal decomposition and an autoencoder. Though neither method is guaranteed to preserve the topology in general,¹⁷ the agreement between the two increases our confidence of the accuracy of the topology and, furthermore, provides strong evidence that these dimensionality reduction methods are sufficiently accurate to capture the behavior of the system in this case.

The main result of this paper is Fig. 13 (or equivalently Fig. 14). We now have a qualitative interpretation of the dynamics: as trajectories move around the chaotic attractor, they are divided into two branches, which intertwine and eventually merge. The bounding torus for this attractor has one hole, unlike the famous Lorenz attractor, but like the classical ODE systems of the Rössler and Sprott D attractors.⁴⁶ As with the templates for those latter ODE systems, the KSE attractor exhibits two branches, one of which has an additional twist, which demonstrates stretching-and-folding chaos (as opposed to “tearing-and-squeezing,” the other standard mechanism).

Topological descriptions of chaotic dynamics give a way of quantifying the stretching and folding processes which give rise to chaos. In a three-dimensional system like the Rössler system, this folding is directly visible from looking at the state space of the attractor. In higher dimensions, projections and dimensionality reduction methods do not immediately provide this information, but a diagram of the template, as given in Figs. 13 and 14, provide an immediate and intuitive interpretation of the dynamics. A topological description of a chaotic dissipative PDE is particularly useful, as it can be used to compare between different dimensionality reduction techniques, and different discretizations of the system, to confirm that the essential dynamics are preserved, even if the exact quantitive details of the system may have been distorted.

The particular system and parameters studied here were chosen so that this work would be possible, and it is not directly generalizable to other dissipative PDEs: we stress that not only is an attractor with fractal dimension less than three necessary, we also must find an embedding of this into three-dimensional Euclidean space. Indeed, in this case, it was proven that a symbolic dynamics exists, which gave further expectation that our approach would be successful.

However, hope for applying topological methods to other PDEs is given by more modern techniques: assigning an $n$-dimensional cell complex to a cloud of points on the attractor allows one to calculate invariants such homology groups, Betti number and Euler characteristic which provide information about the topological structure.^51–54 BraMAH (branched manifold analysis through homologies) can be applied to identify torsions, branch crossings, and weak boundaries.^55,56 Recent work has improved this through the notion of templex⁵⁷ by endowing it with an oriented graph encoding the direction of the flow. A templex associated with a four-dimensional system is constructed in Charó et al.⁵⁷ Though unlikely to be successful in the full state space of the discretized PDE, given the computational complexity of such a high-dimensional system, this is certainly possible in four or five dimensions, and so our approach of using an autoencoder to reduce the dimension can be combined with templexes to study the topology of more complicated PDE systems.

The authors are particularly indebted to an anonymous referee for finding an error in an initial version of the graphical representation of the template at the heart of this paper. The authors thank Omid Ashtari for fruitful discussions. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant No. 865677).

The authors have no conflicts to disclose.

Marie Abadie: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Pierre Beck: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jeremy P. Parker: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Tobias M. Schneider: Conceptualization (equal); Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

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