Predicting chaotic statistics with unstable invariant tori

Chaotic dynamics arise naturally from simple interactions in many physical systems, from fluid dynamics to electrical circuits and nonlinear optics. Studying the chaotic dynamics in terms of simple unstable, non-chaotic solutions to the underlying evolution equations, which are embedded within the stable chaotic attractor, provides key insights into the observed physics. In the absence of special symmetries, two types of special unstable solutions are generally studied: equilibria, zero-dimensional unstable fixed points in the state space of the system; and periodic orbits, time-periodic solutions corresponding to one-dimensional closed loops in state space. Though there is no a priori reason to expect large numbers of equilibria, analytical results^1,2 suggest that, for certain classes of chaotic systems, infinitely many periodic orbits are embedded within the attractor, densely covering it. Trajectories within chaotic attractors closely shadow these unstable periodic orbits (UPOs), and, consequently, periodic orbits are often described as the “backbone” of chaos. Dynamical averages of quantities of interest can be expressed in terms of dynamical zeta functions,³ and if the full collection of periodic orbits is known for a hyperbolic system, these can, in turn, be evaluated by summing over UPOs.^4–6 For systems exhibiting continuous symmetries, which are common in physics, extensions of these results imply summing over relative periodic orbits (RPOs) instead of periodic orbits.^7,8

Of course, for most realistic systems, it is never possible to know the full infinite collection of UPOs, and so these exact results of ergodic theory are of limited direct practical use. In practice, for the complicated spatiotemporal chaos with a large number of positive Lyapunov exponents, which is typical in fluid-dynamical systems, it is not usually possible to find very large numbers of unstable (relative) periodic orbits, or even to quantify how many have been missed, and so the formal periodic orbit theory cannot be applied. Nevertheless, many authors have successfully found moderate numbers of equilibria and periodic orbits in such flows in computations^9–13 and experiments,^14,15 and using ad-hoc choices of weights, reasonable predictions have been achieved for flow statistics by summing over these solutions.^9,13

On a periodic domain, relative periodic orbits are, in fact, a special case of invariant tori. Invariant two-tori are, after equilibria and periodic orbits, the next simplest solutions in continuous dynamical systems. They represent quasiperiodic behavior, in which two different fundamental frequencies interact. In a previous paper,¹⁶ it was shown that for a particular dissipative system of ordinary differential equations that exhibit hyperchaos with three positive Lyapunov exponents, invariant two-tori are generically found embedded within the attractor, appearing whenever the UPOs undergo Neimark–Sacker bifurcations. Suri et al.¹⁷ studied a hyperchaotic¹¹ fluid-dynamical system, and by searching the unstable manifolds of equilibria and UPOs, demonstrated the existence of an unstable invariant two-torus embedded within the attractor.

These pieces of evidence suggest the hypothesis that unstable invariant tori (UITs), as well as unstable periodic orbits, should be generic within hyperchaotic systems. If this is the case, we propose that both classes of solutions should be taken into account when constructing the ad-hoc sums, which were mentioned above. Generic invariant tori could allow the study of key phenomenology within certain fluid-dynamical systems, for which periodic orbits are rare or at least computationally problematic to detect. For example, in wall-bounded turbulence, it has proven difficult to find periodic orbit solutions which capture the interaction between different processes at different length-scales.¹⁸ In the absence of phase locking, we expect the different temporal frequencies associated with these length-scales to lead to the dynamics manifesting as invariant tori rather than periodic orbits. It is conjectured that as the complexity of spatiotemporal chaos increases, which is associated with a greater number of positive Lyapunov exponents, the likelihood of finding and relative importance of UITs increases. However, a method to find such tori generically remains elusive.

When compared with the computation of periodic orbits or relative periodic orbits, the convergence of invariant tori, even given a good initial guess, is a difficult computational problem. Various previous authors^19–26 have developed methods for the continuation of tori, and have successfully applied these to converge attracting tori and continue them as they become unstable, or to converge specific tori that arise from bifurcations. The majority of previous studies have concentrated on Hamiltonian dynamical systems, where invariant tori are known to be a generic feature. Nevertheless, in dissipative systems, (stable) invariant tori have been studied, as the breakdown of a two-torus is one of the possible routes to chaos.^27–29 

This paper aims at advancing the methodology for generically finding and exploiting UITs in the description of chaos. To find generic tori inside a chaotic attractor, we need both a method that can converge very unstable tori from relatively poor initial guesses and a way to find such guesses. In order to address the first point, in this paper, we present a new convergence method, which exploits the quasi-periodicity of non-phase-locked UITs to parameterize their full surfaces. This allows an entirely local formulation of the problem, which is amenable to massive parallelization on GPUs, and thereby we are able to apply this to a partial differential equation (PDE) system. In order to address the issue of finding guesses for UITs, we consider a forced generalized Kuramoto–Sivashinsky equation (gKSE). Here, we can find invariant tori through their connection to relative periodic orbits in the unforced case, which avoids the difficult open problem of detecting UITs directly in chaos. We show that UITs are plentiful and can be summed over to predict various average quantities of the chaotic dynamics.

This paper is laid out as follows: in Sec. II, we present both the unforced and forced gKSE, and the Lyapunov exponents associated with our chosen parameters; in Sec. III, we present and apply methods to find UITs and their stability properties, which are used in Sec. IV to estimate dynamical averages for the system; and the results are discussed in Sec. V.

Figure 2 shows a projection of the chaotic attractors at the discussed parameter values yielding the classical KSE, the unforced gKSE, and the forced gKSE. The choice of projection is important here—since $ϵ≠0$ breaks the continuous phase-shift symmetry, we choose a projection which factors out this phase shift, namely, the absolute value of the complex Fourier coefficients.

We time-march these PDEs using an exponential time-differencing fourth-order Runge–Kutta scheme, following Kassam and Trefethen.³⁷ This is coded in Julia so that we may find Jacobians through automatic differentiation using the package Zygote.³⁸ Some additional calculations were performed using an identical algorithm in MATLAB. Throughout, we discretize the domain with $24$ points, or $N=16$ Fourier modes after 2/3 dealiasing, which is sufficient for this relatively low value of $L$⁠.

The average rate of growth or decay of perturbations to a trajectory within a chaotic attractor is measured by the Lyapunov exponents. The presence of at least one positive Lyapunov exponent is a necessary and sufficient condition for an attractor to be chaotic. A chaotic attractor is called hyperchaotic if its Lyapunov spectrum contains at least two positive exponents.

We compute the leading Lyapunov exponents using Algorithm 1, which is inspired by an approach initially proposed by Benettin et al.³⁹ In this algorithm, the state $u$ and the orthonormal basis $Q$ are mapped to $f τ(u)$ and $J u τ Q$⁠, respectively, where $f τ$ is the time evolution operator for time $τ$ and $J u τ= ∇ u f τ(u)$ is the Jacobian of the flow. The deformed basis is reorthonormalized via QR decomposition, which at the same time allows us to extract the finite-time evaluation of the Lyapunov exponents from the diagonal elements of the right triangular matrix.⁴⁰ This is repeated for $n≫1$ segments of fixed duration $τ$ along a long chaotic trajectory. The Lyapunov exponents $χ i$ are then computed by averaging the $n$ finite-time evaluations. Prior to this loop, the dynamics are integrated for a sufficiently long time $τ 0$ to ensure that the trajectory is confined to the chaotic attractor.

We employ a slight modification to the usual algorithm: instead of using arbitrary linearly independent vectors to construct the initial basis $Q ~$⁠, Algorithm 1 fills the first column of $Q ~$ with $∂ tu$⁠, and—if $ϵ=0$ and continuous translations in $x$ are allowed—the second column with $∂ xu$ (lines 3–6). Therefore, the first (two) column(s) of $Q ~$ span the neutral subspace of the Jacobian. The neutral action of the Jacobian on this subspace is enforced at each iteration by filling the first or the first two columns of $W$ in a similar manner (lines 12–15). As a result, $X$ starts with the expected zero values whose number equals the dimension of the neutral subspace, followed by the nonzero exponents in descending order. This enables us to distinguish exponents with very small absolute values from zeros associated with continuous symmetries.

The four leading Lyapunov exponents of the chaotic attractor over $2 − 7< 10 3ϵ< 2 4$ (with $β=0.01$ and $L=22$ being fixed) are shown in Fig. 3. The parameters $τ 0=2.5× 10 3$⁠, $τ=2.0$⁠, and $n=7.5× 10 5$ are fixed in the calculations. Outside the plotted range of $ϵ$⁠, the attractor behaves as follows: for large values of $ϵ$⁠, the global attractor is a stable fixed point; as $ϵ$ decreases, the system becomes chaotic via the Ruelle–Takens–Newhouse route to chaos: at $ϵ≈0.1202$⁠, the fixed point loses its stability and the attractor becomes a limit cycle; at $ϵ≈0.1021$⁠, the attractor becomes a stable torus; and at $ϵ≈0.0998$⁠, the attractor turns chaotic and remains chaotic until $ϵ≈0.0819$⁠. Below this value, the attractor—or sometimes coexisting attractors—consist(s) of stable periodic orbits or stable tori. At $ϵ≈0.0190$⁠, the attractor again becomes chaotic and remains so for all smaller values of $ϵ$⁠.

Due to the presence of the continuous translational symmetry in the case of unforced $ϵ=0$ system, we do not expect to find many periodic orbits but rather relative periodic orbits (RPOs). For $β=0$⁠, the combination of the continuous and discrete symmetries leads to a dense class of periodic orbits called preperiodic orbits by Cvitanović et al.,⁴¹ but these are not present when $β≠0$⁠. As the continuous symmetry is broken by $ϵ>0$⁠, the remaining RPOs should generically become invariant two-tori, as discussed in Sec. III B. Intuitively, for two-tori to be embedded within a hyperbolic chaotic attractor, we require more than one positive Lyapunov exponent—with only one, the chaotic attractor locally looks like a very thin sheet, which cannot geometrically include an invariant two-dimensional non-chaotic manifold. Indeed, we see two positive Lyapunov exponents for $ϵ>0$⁠, so this is consistent. Note that the Kaplan–Yorke dimension $D K Y$ of the attractor slightly decreases with increasing $ϵ$⁠, while the first and fourth Lyapunov exponents become more negative (see Table I). In all cases, $D K Y≈4.2$—the system we are studying exhibits low-dimensional chaos, in contrast with the spatiotemporal chaos seen at higher $L$ in the KSE, for example.

The fact that $R$ and $S$ are independent of $ρ$ and $σ$ partially constrains the choice of parameterization, which improves numerical convergence, but also assumes quasiperiodic dynamics. This is a strong assumption, which is certainly not valid in Arnold tongues, regions of parameter space in which phase locking implies the existence of periodic orbits on an invariant torus. In such regions, we may instead directly converge the periodic orbits.¹⁶ In practice, when the phase locking of a UIT is such that the periodic orbits on it are of very long period, it is numerically indistinguishable from a quasiperiodic invariant torus, and our algorithm will converge.

We start with an initial guess, which geometrically describes a two-torus in state space, but not an invariant manifold. We then iteratively deform this torus until (6) is satisfied at every point on the surface. The torus is discretized with $N=16$ Fourier modes in $x$⁠, as for (2), so the linear terms are entirely local differential operators in Fourier space. In $ρ$ and $σ$⁠, we discretize on an evenly spaced grid of $M×M=64×64$ points, and the derivatives are calculated with dense trigonometric differentiation matrices. This means that we can compute a sparse Jacobian for the left hand side of (6), of size $N M 2×(N M 2+2)$⁠, where the two additional columns come from derivatives with respect to $R$ and $S$⁠. The matrix consists of dense $N×N$ blocks on the main diagonal, corresponding to the nonlinear terms, and diagonal $N×N$ blocks elsewhere, corresponding to the derivatives with respect to $ρ$ and $σ$⁠. This gives a total of $M 2 N 2+ M 4N+ M 2N$ nonzero entries. Following Cvitanović et al.,⁴¹ we can solve this system using a Levenberg–Marquardt algorithm. As a consequence, we do not need to introduce additional constraints for this underconstrained optimization problem. At each iteration of Levenberg–Marquardt, the sparse linear system is solved using the conjugate-gradient-like least squares minimum residual.⁴² Since this requires only sparse matrix multiplications, this can be performed on a GPU in Julia at great speed despite the size of the system.

As a consequence of using the Levenberg–Marquardt algorithm, our method is able to converge from relatively poor initial guesses, since it is effectively gradient-descent when far from a solution, and so shares the convergence properties of recent improvements in methods for the computation of periodic orbits.^12,13,43 Figure 5 shows the torus-shaped but not invariant initial guess at $ϵ=0.002$ (generated by continuation, see Sec. III B), and the drastically different converged solution.

Our method is also able to converge some UPOs which exist when strong phase locking is present, as UPOs are a degenerate case of the solution given by (6) in which either $∂ ρu≡0$ or $∂ σu≡0$⁠. Indeed, if they were both zero, the algorithm would have converged to a steady state solution, but this never occurred.

A relative periodic orbit is a particular case of an invariant torus, where one of the two dimensions of the invariant manifold exactly corresponds to a continuous symmetry of the system. In the presence of such a symmetry (and the absence of discrete symmetries), RPOs are expected to be the generic structure in chaotic attractors, with pure periodic orbits rare special cases, for example, related to bifurcations of solutions outside the chaos. The unforced system with $ϵ=0$ has a continuous symmetry, and we can find RPOs in it using the now-routine method of recurrent flow analysis, followed by a Newton–Krylov based shooting method.⁹ Our implementation of recurrent flow analysis exactly matches that of Cvitanović et al.⁴¹ Since $β≠0$⁠, the system has no discrete symmetry, and so we do not anticipate the existence of “pre-periodic orbits,” a special case of (non-relative) periodic orbits.⁴¹ 

It was found to be sufficient to use the RPO as an initial condition for the algorithm at $ϵ=0.001$ and then use this converged UIT and the $ϵ=0$ RPO to linearly extrapolate an initial guess for $ϵ=0.002$⁠. Several of the UITs were continued to $ϵ=0.005$ and beyond—Fig. 6 shows an RPO at $ϵ=0$ and its continuation to a UIT at $ϵ=0.01$—but for the more geometrically complex tori this proved difficult and higher resolution discretizations of $ρ$ and $σ$ are likely to be required. Indeed, we only searched for RPOs for periods up to $T=100$ for this reason. In other cases, strong phase lockings were detected, at which point it is no longer possible to satisfy (6) and the continuation algorithm breaks down. Such phase lockings are generic when invariant tori are continued, as Arnold tongues are encountered.⁴⁴ It is nevertheless still possible to continue the UITs past these tongues with considerable additional effort,¹⁶ but this was not performed here.

Figure 7 shows two different projections of a simple RPO at $ϵ=0$⁠, which is continued to a UIT at $ϵ=0.01$⁠. With a projection showing the absolute values of the Fourier coefficients, which are invariant under shifts in the $x$ direction, the RPO appears as a simple loop but at non-zero $ϵ$⁠, the structure is clearly a full two-dimensional structure. A more naïve projection showing the real parts of the Fourier coefficients, which are not invariant under phase shifts, shows that the RPO is indeed topologically a torus at $ϵ=0$⁠, but then the difference between the UIT and the RPO is obscured.

Twelve distinct UITs were successfully converged at $ϵ=0.001$ from the continuation of 27 distinct RPOs. A projection of all 12 of these is shown in Fig. 8. Of these, eight were successfully continued to $ϵ=0.002$⁠. In addition to these, a common traveling wave solution at $ϵ=0$ was continued to give a periodic orbit, but this was omitted from our calculations in Sec. IV, as were any UPOs resulting from phase lockings of the UITs.

Algorithms for finding the stability properties of invariant tori are notoriously difficult to use.^45,46 Here, we propose a simple iterative method which accurately finds the leading Lyapunov exponents for a UIT, i.e., the real parts of the stability eigenvalues. It is sufficient to know these to calculate the weights for the statistical averages discussed in Sec. IV. Our method does not give the imaginary parts of the eigenvalues or the linear manifolds associated with the eigenvectors.

The algorithm is an extension of that discussed in Sec. II A to calculate Lyapunov exponents for the attractor. In theory, if we start with a point exactly on an invariant set, the Lyapunov exponents calculated using the given algorithm should give us exactly the Lyapunov exponent of the invariant set. However, as the sets of interest are all unstable, a trajectory started from a calculated point of the invariant set will quickly drift away. Consequently, we augment the algorithm and in each iteration, we ensure that the trajectory starts from the correct point on the invariant manifold, assuming the dynamics are given locally by (5), as described in Algorithm 2.

All the UITs found had either two unstable eigenvalues (distinct values or repeated values, the latter suggesting a complex conjugate pair) or only one unstable eigenvalue. Any solution embedded within the chaotic attractor can have at most two unstable eigenvalues, since this is the number of positive Lyapunov exponents for the attractor itself, as discussed in Sec. II. Any UIT should have, by definition, two zero eigenvalues, and in each case, we numerically found two Lyapunov exponents with magnitude very close to zero. The same method applied to UPOs gives just one (approximately) zero eigenvalue and one or two positive ones. The full results are listed in Table II.

Tables II and III list the UITs that were successfully converged at $ϵ=0$ and continued to $ϵ=0.001$ and $ϵ=0.002$⁠, respectively. The calculated Lyapunov exponents $χ 1$ and $χ 2$ for each of them are used to compute weights, and these weights are used to give predictions for the long-time averages of the energy production/dissipation and the first three Fourier coefficients in the chaotic attractor. The formal periodic orbit theory suggests that it is possible to calculate the Lyapunov exponents for the attractor by summing over UPOs,⁵ and, thus, we have also attempted this here using UITs.

For the Lyapunov exponents, which are properties of a full trajectory rather than instantaneous quantities, the standard deviation is not well-defined and so we do not give a normalized error. Nevertheless, from the values in the table, we immediately see that the predicted values are rather poor estimates of the true Lyapunov exponents, when compared against the difference between the true values.

We have demonstrated that, for this contrived system at carefully chosen parameter values, UITs are common and readily found. Based on the projection shown in Fig. 8, all identified UITs are located within the attractor and collectively give a convincing visual representation of it. This suggests that the UITs may indeed be spanning a significant part of the attractor.

For these particular collections of UITs, we have predicted the energy production/dissipation to within less than 5% of one standard deviation, despite the measured values associated with individual tori varying by much more than the standard deviation. We interpret this observation as evidence that our particular collections of tori, which were generated from random chaotic timeseries, represent a reasonably unbiased sample of the dynamics.

The predicted values of the Lyapunov exponents of the chaotic attractor are significantly poorer, which is to be expected as these are notoriously difficult to calculate and sensitive to which parts of the attractor are accounted for. Whether the predicted values become more precise as more UITs are included—or whether instead there is some systematic reason that our predictions of the Lyapunov exponents are poor—remains an open question.

By definition, if UITs exist within an attractor, chaotic trajectories must pass arbitrarily close to them. We might well expect, therefore, that given a number of diverse UITs that cover much of the attractor we should be able to average quantities over the tori to approximate averages over the attractor itself. This system was deliberately chosen to give a method to find large numbers of UITs, and so it is no great surprise that our method is successful here.

What remains an open question is whether UITs exist in large numbers in realistic systems for which it has been difficult to detect UPOs, such as wall-bounded turbulence. If they do, they could capture physical processes that are currently poorly understood.¹⁸ Even if such UITs exist in general high-dimensional, dissipative, chaotic systems, their detection and convergence remains a computational challenge. In this work, we have specifically constructed a system in which we can exploit the breaking of a continuous symmetry to find UITs; this method will not be applicable in more generic systems. The detection of UITs directly from a chaotic timeseries is an open problem and will certainly require more sophisticated data analysis methodologies than the traditional recurrent flow analysis used to find UPOs.

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.

Jeremy P. Parker: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Omid Ashtari: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Tobias M. Schneider: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

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