Jacobian-free variational method for computing connecting orbits in nonlinear dynamical systems

A broad spectrum of physical systems, from fluid flows¹ to nonlinear optics^2,3 and suspensions of motile micro-organisms,^4,5 exhibit spatiotemporally chaotic dynamics. In the framework of dynamical systems, the spatiotemporal chaos is viewed as the evolution of a chaotic trajectory in the state space of the governing equations. Embedded in the state space are non-chaotic, time-invariant solutions including equilibria, periodic orbits, and invariant tori. These invariant solutions are dynamically unstable so that the chaotic trajectory visits them transiently, yet recurringly. Spatiotemporal chaos can, thus, be viewed as a walk through a forest of invariant solutions that form the elementary building blocks of the chaotic solution.^6–8 Consequently, individual invariant solutions are able to capture essential properties of the observed spatiotemporal structures, and, collectively, they promise an avenue toward quantitatively predicting statistical properties of the chaotic dynamics. The increasing computational resources and algorithmic advances have enabled these concepts, originally developed in the context of low-dimensional chaotic dynamical systems, to be applied to high-dimensional problems including transitional turbulence where descriptions based on invariant solutions have proven to be particularly useful.^9–14 

While equilibria and periodic orbits form the building blocks of the dynamics, the chaotic evolution from the neighborhood of one unstable invariant solution to another is mediated by connecting orbits. These hetero- and homoclinic connections provide dynamic pathways between different periodic orbits or equilibria within the chaotic attractor. Therefore, a complete dynamical description of the chaotic dynamics in terms of state-space structures requires to both identify equilibria, periodic orbits, and invariant tori embedded in the chaotic attractor and compute connecting orbits between them. In the context of fluid dynamics, for example, Kawahara and Kida¹⁵ and van Veen and Kawahara¹⁶ use connecting orbits to explain the turbulent bursting in plane Couette flow; Suri et al. study the network of connecting orbits that underpins the transient dynamics in a quasi-two-dimensional Kolmogorov flow;¹⁷ Reetz and Schneider characterize the time-dependent dynamics of inclined layer convection using connecting orbits between coexisting invariant solutions;¹⁸ and De Jesús and Graham discuss a network of dynamically connected relative periodic orbits that organizes the attractor in a two-dimensional Kolmogorov flow.¹⁹ 

We specifically focus on connecting orbits between equilibrium solutions. Such connecting orbits have been identified as dynamically relevant in fluid systems,^17,18,20 and they are involved in global bifurcations when, for instance, a periodic orbit bifurcates off a homoclinic orbit or a heteroclinic cycle.^21,22 Connecting orbits are located within the intersection of the unstable manifold of one equilibrium with the stable manifold of another or the same equilibrium solution if they are of heteroclinic or homoclinic type, respectively. In the vicinity of an equilibrium solution, a trajectory approaches/departs the equilibrium along its stable/unstable manifold exponentially in time. Consequently, the time required to traverse the entire connecting orbit is not finite. This infinite passage time makes computing connecting orbits very challenging.

One approach to handle the computational challenge of the infinite passage time is to truncate the connecting orbit and compute part of the orbit that is traversed in finite time. Under favorable conditions, the truncated orbit can be computed using shooting methods. Geometrically, the truncation approach searches for an initial condition on the unstable manifold of the origin equilibrium whose resulting trajectory lands on the stable manifold of the destination equilibrium. Due to their curvature, parametrizations of stable and unstable manifolds are usually not accessible. Consequently, they need to be approximated locally by the corresponding tangent spaces associated with the origin and destination equilibrium. Practically, a connecting orbit is, thus, approximated by identifying an initial condition in the intersection of the unstable tangent space of and a hypersphere around the origin equilibrium, which after forward time integration reaches a distance below a chosen threshold from the destination equilibrium.²³ If the hypersphere is chosen small enough, the unstable tangent space accurately approximates the unstable manifold, and thus, the obtained trajectory accurately represents a connecting orbit.

Even if the unstable manifold can be accurately approximated by the unstable tangent space, a systematic search for an initial condition that eventually reaches the destination equilibrium is a formidable task, especially for a chaotic system where nearby trajectories diverge exponentially with time. When the unstable manifold at the origin equilibrium solution is two-dimensional, an exhaustive search strategy can be employed.^11,17,24,25 In this case, the search space is a circle on the unstable tangent space with an angle being the only variable. However, when the unstable tangent space at the origin equilibrium has more than two dimensions, the search space is too large for an exhaustive search. To improve the dimensionality drawback, Farano et al.²⁶ propose an adjoint-based variational method for finding a state on an energy shell around the origin equilibrium whose trajectory reaches another energy shell around the destination equilibrium. They do not constrain the initial condition to be located on the unstable tangent space at the origin equilibrium, hence as a second step, the trajectory is confirmed to shadow a connecting orbit by matching the end points of the trajectory against the linearized dynamics around the two equilibria. In all these methods determining the size of the hypersphere around the origin equilibrium solution is not a trivial task: the hypersphere should be small enough in order for the tangent space to accurately approximate the manifold, and large enough to let the required time integration intervals be feasibly short.

An alternative to the shooting-based methods, which search for a single state on the connecting orbit, is to search in the space of connecting curves, i.e., all smooth curves in the state space which connect the two equilibria. Among all such curves, only connecting orbits are integral curves of the vector field induced by the governing equation. The idea is to start from a connecting curve pivoted on the two fixed points, then deform the curve until the tangent velocity coincides with the local field vector along the entire curve, and thus, a connecting orbit is achieved. This approach has several advantages over the reviewed shooting-based methods for computing connecting orbits: First, there is no limitation on the dimensionality of the unstable manifold at the origin equilibrium because no exhaustive search is needed; Second, the approach does not suffer from the exponential separation of trajectories with time since the connecting curve is deformed locally and no time integration is required; and finally, this approach yields the exact and the entire connecting orbit without requiring to truncate it.

Despite the conceptual advantages of searching in the space of connecting curves over the shooting-based alternatives, this approach is not extensively developed on the practical side. Liu et al.²⁷ use rational Chebyshev basis functions for the spectral representation of variables along the infinite temporal direction probably for the first time in this context. They formulate the problem as a system of nonlinear equations by setting the temporal derivative equal to the right-hand side of the governing equation for every state variable at every temporal collocation point. They then solve the resulting system of equations using standard Newton iterations. Dong and Lan²⁸ extend the variational method of Lan and Cvitanović,²⁹ originally developed for finding periodic orbits, to the problem of computing connecting orbits. They view the problem of deforming connecting curves toward a connecting orbit as a minimization problem: a connecting orbit is found by minimizing a cost function, which penalizes the deviation of a connecting curve from being an integral curve of the vector field. They employ an infinitesimal-step version of Newton iterations for continuously deforming the curve and use finite differences for calculating the tangent velocity vector. In his Ph.D. thesis, Pallantla³⁰ employs the same spectral representation of variables in the temporal direction as in Ref. 27 and deforms the curve in the direction of the steepest descent of the cost function. The common drawback of the aforementioned algorithms is that they all require explicit construction of the Jacobian matrix. In a system with $M$ temporal and $N$ spatial degrees of freedom, the size of the Jacobian matrix scales as $O( M 2 N 2)$⁠, which can be prohibitively large for high-dimensional dynamical systems such as three-dimensional fluid flows.

In order to transfer the advantages of searching in the space of connecting curves to high-dimensional dynamical systems, we propose a Jacobian-free variational method for computing connecting orbits between two equilibrium solutions. The method employs an adjoint-based optimization technique to minimize a cost function, which measures the deviation of a connecting curve between two equilibria from an integral curve of the vector field. We construct a globally contracting dynamical system in the space of connecting curves. Fixed points of this dynamical system are minima of the non-negative cost function. The global minima of the cost function, taking zero value, correspond to connecting orbits of the original dynamical system. Therefore, connecting orbits are found by integrating the dynamics in the space of connecting curves. Due to the explicit construction of the dynamical system in the space of connecting curves, the memory requirement scales as $O(MN)$ which allows the proposed method to be applied to high-dimensional dynamical systems.

The remainder of the present article is organized as follows. In Sec. II, the problem of computing a connecting orbit is set up as a minimization problem, and in Sec. III, the adjoint-based minimization technique is formulated for a general autonomous dynamical system. In Sec. IV, a spectral representation suitable for the discretization along the unbounded temporal domain is discussed. To demonstrate the robustness of the proposed variational method, in Sec. V, we consider the one-dimensional Kuramoto–Sivashinsky equation in a spatiotemporally chaotic regime and show that several connecting orbits can be converged reliably. Finally, in Sec. VI, the paper is summarized and an outlook for future research is given.

The dynamical system $∂ u/∂τ= G( u)$ is constructed in the space of connecting curves $C g$ defined in Eq. (3). This imposes certain BCs on the residual $r( u)$ and the operator $G( u)$⁠. In the temporal direction, $lim s → ± ∞ r= 0$ since $u$ satisfies the correct asymptotic BCs for all $τ$⁠, and $lim s → ± ∞ G= 0$ since the correct asymptotic values of $u$ must be preserved. In space, $u$ satisfies the correct BCs at $∂Ω$ for all $τ$⁠; consequently, the spatial BCs of $r$ and $G$ are determined following similar arguments. For example, $r$ and $G$ will be periodic in directions where $u$ is periodic, will take zero value at the boundaries where $u$ satisfies Dirichlet boundary conditions and so forth. These properties must be taken into account while deriving the adjoint operator from the definition (14). Derivation of the adjoint operator for the Kuramoto–Sivashinsky system, introduced in Sec. V, is presented in  Appendix A where the zero asymptotic values of $r$ and $G$ in the temporal direction and their periodicity in space enable us to derive the adjoint operator as an explicit function of the space–time field $u$⁠.

Both heteroclinic and homoclinic connections can be numerically identified using the introduced variational method. In the case of a homoclinic connection to an equilibrium solution, zero variation in time, i.e., the equilibrium solution itself, is a trivial solution satisfying the definition (2). Therefore, depending on the initial connecting curve from which the integration starts, a trivial or a nontrivial solution with $J=0$ can be obtained. The definition of a heteroclinic connection does not have any trivial solution.

On an abstract level, we construct the operator $G$ following the same logic as that employed by Farazmand³¹ in computing equilibria and by Azimi et al.³² in computing periodic orbits. However, in the different contexts, the form of this operator differs as the resulting variational dynamics evolves a spatial or space–time field representing the specific sought-after invariant solution: In computing equilibria, the variational dynamics evolves a spatial field that geometrically corresponds to moving a point in the state space; in computing periodic orbits, the variational dynamics evolves a space–time field periodic in the temporal direction that corresponds to deforming a loop in the state space, and here the variational dynamics evolves a space–time field satisfying the asymptotic conditions (2) in the temporal direction that corresponds to deforming a curve in the state space that connects two fixed points.

An efficient implementation of the proposed adjoint-based variational method is aided by an accurate spectral representation of a space–time field $q( x,s)∈ C s$ in the $s$ direction such that the asymptotic conditions at $s→±∞$ are directly enforced by the chosen expansion. The spectral accuracy significantly reduces the number of time sections, and thereby memory, required for an accurate representation of connecting orbits. We use rational Chebyshev basis functions for the spectral representation in the temporal direction (see Chap. 17 of Ref. 33 for details).

The expansion in a rational Chebyshev basis allows us to represent the space–time objects in the unbounded temporal direction, and we can expect spectral accuracy with fast convergence as a function of the expansion’s truncation order. Rational Chebyshev functions form a generic basis for the spectral representation of functions over the entire real axis with constant asymptotic values and are, thus, a suitable expansion for connecting orbits for any studied physical system.

As a proof of concept, we apply the introduced method for identifying connecting orbits to the one-dimensional Kuramoto–Sivashinsky equation (KSE).^34,35 The KSE is a nonlinear partial differential equation, which emerges in various physical contexts such as flame propagation,³⁵ plasma physics,³⁶ or interfacial fluids instability.³⁷ The KSE is also commonly used as a model system for examining new methods developed for chaotic fluid flows and transitional turbulence since it exhibits spatiotemporally chaotic behavior and displays some similar features to the Navier–Stokes equations.

For $L<2π$⁠, the trivial equilibrium solution $u(x,t)= const .$ is linearly stable and is the global attractor of the dynamics. By increasing $L$⁠, solutions of the KSE undergo a series of bifurcations before the dynamics can exhibit spatiotemporally chaotic behavior for sufficiently large domain sizes. For certain ranges of the control parameter, heteroclinic cycles between equilibria emerge as the most attracting invariant sets. For instance, for $13.03≲L≲14.92$⁠, a network of isolated heteroclinic orbits between two symmetry-related equilibrium solutions attracts the nearby trajectories, and the dynamics exhibits regular transitions between the two equilibria.^38,39 For $18.05≲L≲20.89$⁠, we again observe the network of connecting orbits between two symmetry-related equilibria to be attracting. In this range of $L$⁠, the two equilibria are connected by a one-parameter continuous family of connecting orbits which provides infinitely many heteroclinic cycles and, therefore, transitions take place via irregular bursts.

We demonstrate the application of the proposed method by computing connecting orbits between equilibrium solutions of the KSE for $L=22$⁠. This domain size is large enough for the KSE to sustain the spatiotemporally chaotic dynamics. Hence, there is no attracting heteroclinic cycle, and connecting orbits cannot be obtained by time-marching the KSE. However, this domain size is small enough to have low-dimensional unstable manifolds at the equilibria found, over which an exhaustive search for possible connecting orbits is practical. The state space geometry of the KSE for this parameter value has previously been explored in detail by Cvitanović and collaborators.²⁵ They identified several connecting orbits using the shooting method described in Sec. I. Here, at least one connecting orbit between any two equilibrium solutions is computed, or it is confirmed by the exhaustive search in Ref. 25 that no connecting orbit exists between the two equilibria.

The KSE (22) preserves the spatial mean value of the evolving field. Consequently, the spatial mean along a connecting orbit is constant and the same as the end point equilibrium solutions. We consider the dynamics of the KSE in the subspace of fields with zero spatial mean. The zero mean value is not enforced during the evolution of a connecting curve toward a connecting orbit. However, since the two end point equilibria do have zero spatial mean, a converged connecting orbit with $J=0$ takes zero mean value as well.

In spectral space, the connecting curve $u$ is represented by an $(M+2)×N$ matrix of complex numbers $u ^ m , j= u ^ j( s m)$⁠. The derivative of order $q∈ W$ of this space–time field with respect to $x$ is obtained by the Hadamard product $D x ( q )⊙ u ^$⁠, where $D x ( q ) m , j= ( 2 π j i / L ) q$⁠, and its derivative of order $q∈ W$ with respect to $s$ is obtained by multiplying $u ^$ from the left by $D t q$⁠, where the temporal differentiation matrix $D t$ is defined in Eq. (21). The residual $r$ and the descent direction $G$ are discretized in the same way with the only difference that their time sections at $s 0$ and $s M + 1$ (corresponding to $s→+∞$ and $s→−∞$⁠, respectively) are identically zero (see Sec. III). The nonlinear terms are calculated in physical space where products are of elementwise Hadamard type. Transforming back and forward between physical and spectral representations of the space–time fields requires one-dimensional forward or backward discrete Fourier transformation of each time section.

The defined dynamical system $∂u/∂τ=G$ is globally contracting and we are only concerned about the asymptotic state $u= u 0+ ∫ 0 ∞G dτ$⁠. Consequently, we select the numerical integration scheme based on simplicity and stability rather than accuracy. We use semi-implicit forward Euler time-stepping scheme, which has first-order accuracy in $τ$ and treats the linear terms of $G$ in $u$ implicitly and the nonlinear terms explicitly. The code was developed in C++ with OpenMP parallelization of local calculations.

In the subspace of fields with zero spatial mean, the KSE with $L=22$ has four known equilibrium solutions including the trivial solution $u=0$⁠. Hereafter, we denote the trivial equilibrium solution by $E 0$⁠, and the nontrivial ones by $E 1$⁠, $E 2$⁠, and $E 3$ as shown in Fig. 2. We compute these equilibrium solutions following the adjoint-based variational method of Farazmand.³¹  $E 1$⁠, $E 2$⁠, and $E 3$ are invariant under inversion about the origin, $σ$⁠. $E 2$ and $E 3$ are also symmetric under discrete shifts $γ(1/2)$ and $γ(1/3)$⁠, respectively. Therefore, in addition to the inversion about $x=0$ and $L/2$⁠, $E 2$ is symmetric under inversion about $x=L/4$ and $3L/4$⁠, and $E 3$ is symmetric under inversion about $x=L/6$⁠, $L/3$⁠, $2L/3$⁠, and $5L/6$ as well. The repelling eigenvalues of all four equilibrium solutions are listed in Table I, and their associated eigenvectors are shown in Figs. 15–18 in  Appendix B.

Due to the continuous translational symmetry of the KSE, $E i$ with $i=1,2,3$ represent their so-called group orbits of all symmetry-related states, i.e., $γ(α) E i$ where $α∈[0,1)$⁠. Every connecting orbit, therefore, has infinite dynamically equivalent copies corresponding to similar translations of the origin and the destination equilibrium solutions. We search for connecting orbits of a certain relative phase between the two end points by fixing the origin equilibrium and shifting the destination equilibrium solution when constructing the initial connecting curve using Eq. (29). In the following, we first demonstrate the application of the introduced method by computing a connecting orbit from $E 1$ to $E 2$⁠. We then present converged connecting orbits between other equilibrium solutions and compare them to the same orbits obtained from other methods reported in the literature if applicable.

The search for a heteroclinic connection from $E 1$ to $E 2$ is initialized by a connecting curve constructed using Eq. (29) in the inversion-symmetric subspace of $M$ (⁠ $a=0$⁠). We discretize the space–time domain by $N=64$ Fourier modes in space and $M=550$ rational Chebyshev grid points in time. The scaling of the temporal discretization is set to $S=55$ and the center of the distribution to $s 0=0$⁠. For this system, the integration scheme described in Sec. V C is stable for $Δτ=0.01$⁠.

After a sharp initial decrease, the arc length cost function decays exponentially with the fictitious time, as shown in Fig. 3, and reaches the convergence criterion, $J arc= 10 − 12$⁠, at $τ≈1.25× 10 4$⁠. In the vector field induced by $G$⁠, heteroclinic connections are attracting fixed points. The exponential decay of the cost function suggests that when the evolving connecting curve gets close enough to the connecting orbit, the dynamics is dominated by the leading, i.e., the slowest, eigendirection of the linearized dynamics in the vicinity of the fixed point of $∂ u/∂τ=G$⁠.

Figure 4 shows six snapshots of the continuous deformation of the connecting curve from $E 1$ to $E 2$ governed by the dynamics in the space of connecting curves (25) toward a heteroclinic connection. A substantial deformation toward the final shape of the connecting orbit takes place at the beginning of the evolution. The major remaining part of the integration time is spent on the slight remaining deviation from the final orbit. The space–time field corresponding to the initial connecting curve [snapshot (i) in Fig. 4] and the converged connecting orbit [snapshot (vi) in Fig. 4] are displayed in panels (a) and (b) of Fig. 5, respectively.

The spatial resolution is chosen by monitoring the energy spectrum of spatial Fourier modes in a direct numerical simulation of the KSE for $L=22$⁠. The spatial resolution $N=64$ ensures at least six orders of magnitude drop in the modulus of spatial Fourier coefficients at all times. The converged connecting orbit from $E 1$ to $E 2$⁠, as an equilibrium solution to Eq. (25), is structurally stable for a wide range of temporal resolutions $M$⁠. However, the accuracy of the spectral representation in time, and therefore the minimum achieved value of the cost function, $J arc , min:= lim τ → ∞ J arc(τ)$⁠, varies with $M$⁠. Figure 6 show the spectral convergence of $J arc , min$ with $M$⁠. Notice that $J arc , min$ can be considerably higher than the convergence criterion when $M$ is not large enough. If a local minimum of the cost function is reached, in contrast, $J arc , min$ does not improve as the temporal resolution is increased. As an example of a failing search, we try to converge a connecting orbit between $E 2$ and $γ(1/4) E 2$ from an initial connecting curve constructed using Eq. (29) with $a=0$ (see Sec. V C 3 why such connection cannot exist). The integration from this initial connecting curve does not reach a global minimum but approaches a local minimum with $J arc , min=5.1× 10 − 2$⁠. As shown on Fig. 6, the minimum value does not decrease as the temporal discretization is refined, confirming that a converged local minimum has been identified and no connecting orbit was found.

We converge a heteroclinic connection from $E 0$ to $E 1$⁠, $E 2$ and $E 3$ from an initial connecting curve constructed using Eq. (29) with $a=0$⁠. A three-dimensional state space projection and the space–time contour of heteroclinic connections from $E 0$ to the other three equilibrium solutions are exhibited in Figs. 7 and 8, respectively. The algorithm settings are presented in  Appendix C.

The unstable manifold of $E 0$ is six-dimensional. Each of the repeated unstable eigenvalues of $E 0$⁠, Table I, is associated with one eigenvector symmetric under reflection across $x=0$ and another one symmetric under inversion about the origin (see Fig. 15). An exhaustive search in the unstable tangent space at $E 0$ is not practical even in the inversion-symmetric subspace of the KSE where the reflection-symmetric eigenvectors do not exist, and the unstable manifold is three-dimensional. Dong and Lan²⁸ have computed a heteroclinic connection from $E 0$ to $E 1$ using their variational method, which employs finite differences for calculating tangent velocity vectors. They have used 6000 sections to discretize this connecting orbit in time and obtain residuals of order $O( 10 − 6)$⁠. To achieve this value of $J arc$ (and similarly, the supremum norm of the residual $r$⁠), $M=25$ interior time sections suffice for the proposed variational method.

We demonstrated the details of converging a heteroclinic connection from $E 1$ to $E 2$ at the beginning of this section (see Figs. 3–6). We also converge a heteroclinic connection from $E 1$ to $E 3$ from an initial connecting curve constructed using Eq. (29) with $a=0$⁠. Figures 9 and 10 show a three-dimensional state space projection and the space–time contour plot of the converged heteroclinic connection from $E 1$ to $E 3$⁠, respectively. The algorithm settings are presented in  Appendix C.

The unstable manifold of $E 1$ is four-dimensional. One pair of complex conjugate unstable eigenvalues of $E 1$⁠, Table I, is associated with eigenvectors invariant under reflection across $x=0$⁠, while the other pair is associated with eigenvectors invariant under inversion about the origin (see Fig. 16). An exhaustive search in the four-dimensional unstable tangent space at $E 1$ is not practical. Cvitanović et al.²⁵ perform an exhaustive search in the two-dimensional plane spanned by the reflection-symmetric eigenvectors at $E 1$ and show that all trajectories starting from that plane are chaotic and do not reach any of the equilibrium solutions. They perform another exhaustive search in the two-dimensional plane spanned by the inversion-symmetric eigenvectors and show that trajectories starting from that plane form a one-parameter family of heteroclinic connections from $E 1$ to $E 2$⁠, except one bordering orbit that converges to $E 3$⁠.

We converge two heteroclinic connections from $E 2$ to $E 3$ and $γ(1/4) E 2$⁠. The initial conditions are constructed using Eq. (29) by setting $a=0$ for the connecting orbit between $E 2$ and $E 3$ and setting $a=−1$ and $v=ℜ{ v 1 , 2}$ for the connecting orbit between $E 2$ and $γ(1/4) E 2$⁠, where $ℜ{ v 1 , 2}$ is the real part of the complex conjugate unstable eigenvectors at $E 2$ (see Fig. 17). In the latter, adding the symmetry breaking term (⁠ $a≠0$⁠) is necessary because $E 2$ and $γ(1/4) E 2$ are both symmetric under inversion about $x=kL/4$ with $k=0,1,2,3$⁠, thus, an initial connecting curve constructed by setting $a=0$ is symmetric under inversion about all these points. The dynamics (25) preserves all the four inversion symmetries while no connecting orbit can exist in such subspace of $M$ because the unstable eigenvectors of $E 2$ are symmetric only about $x=0$ and $L/2$⁠, meaning that as soon as a trajectory of the KSE leaves $E 2$⁠, the inversion symmetries about $x=L/4$ and $3L/4$ are broken. Consequently, as shown in Fig. 6, $a=0$ results in getting stuck in a local minimum of the cost function as the dynamics (25) is integrated. A three-dimensional state space projection and the space–time contour plot of the connecting orbits from $E 2$ to $E 3$ and $γ(1/4) E 2$ are shown in Figs. 11 and 12, respectively. The algorithm settings are presented in  Appendix C.

By an exhaustive search in the two-dimensional unstable tangent space at $E 2$⁠, Cvitanović et al.²⁵ show that the unstable manifold of $E 2$ is a one-parameter family of connecting orbits that converge to $γ(1/4) E 2$⁠, except one orbit that connects $E 2$ to $E 3$⁠.

We converge two heteroclinic connections from $E 3$ to $E 2$⁠. The initial conditions are constructed using Eq. (29) by setting $a=0$ in one and $a=−1$ and $v=sin⁡ ( x )$ in the other. A three-dimensional state space projection and the space–time contour of these connecting orbits are shown in Figs. 13 and 14, respectively. The algorithm settings are presented in  Appendix C.

The unstable manifold of $E 3$ is two-dimensional. The repeated positive eigenvalue of $E 3$⁠, Table I, is associated with one eigenvector symmetric under reflection across $x=0$⁠, and another eigenvector symmetric under inversion about the origin (see Fig. 18). Cvitanović et al.²⁵ conduct an exhaustive search in the two-dimensional unstable tangent space at $E 3$ and identify two heteroclinic connections from $E 3$ to $E 2$ corresponding to the perturbation of $E 3$ along the inversion-symmetric eigenvector and its opposite direction. Fixing $E 3$ and shifting $E 2$ in space by $L/3$ and $2L/3$ puts the translated copy of $E 2$ in the same relative phase to $E 3$ as the original configuration. Therefore, the exhaustive search identifies two other pairs of heteroclinic connections from $E 3$ to the group orbit of $E 2$⁠, which are copies of the first pair of connecting orbits shifted by $L/3$ and $2L/3$ in the $x$-direction.

Connecting orbits are of significant importance for studying spatiotemporally chaotic dynamical systems in terms of their invariant state space structures. We introduce a variational method for computing connecting orbits between two equilibrium solutions by searching in the space of all smooth curves in the state space that connect the two equilibria. In this method, the deviation of a connecting curve from an integral curve of the vector field is penalized by a non-negative cost function. A dynamical system in the space of connecting curves is set up such that along its trajectories the cost function is guaranteed to decrease monotonically. All trajectories of this dynamical system eventually converge to an equilibrium, which corresponds to a minimum of the cost function. Global minima of the cost function, taking zero value, correspond to the connecting orbits of the original dynamics. This method is not limited by the dimensionality of the unstable manifold at the origin equilibrium solution, does not suffer from exponential separation of trajectories, and does not require any domain truncation. The introduced method is Jacobian-free, and its memory requirement scales linearly with the number of degrees of freedom, which allows this method to be applied to high-dimensional dynamical systems including three-dimensional fluid dynamics problems.

As a proof of concept, we apply the introduced variational method to the one-dimensional KSE and compute several connecting orbits between known equilibrium solutions of the system with domain size $L=22$⁠. The set of converged solutions contains at least one connecting orbit between any two equilibrium solutions unless it is known from an exhaustive search in the unstable manifold of the origin equilibrium solution that they are not connected. The stability analysis of the adjoint-descent dynamics (25) can determine if the identified solution is an isolated connecting orbit or belongs to a one-parameter continuous family of orbits. In the former case, all eigenvalues are negative except only one zero associated with continuous translations in the temporal direction $s$⁠. In the latter, the stability matrix has an additional zero associated with variation within the continuous family of connecting orbits. In that case, the connecting surface can be computed by tracing the neutral eigendirection associated with the second zero eigenvalue.

After demonstrating the feasibility of the introduced method for computing connecting orbits between equilibrium solutions of the one-dimensional KSE, we are extending the present work in two directions: One is applying this method to the three-dimensional wall-bounded fluid flows governed by the Navier–Stokes equations (NSEs). The challenge in applying this method to the wall-bounded NSEs lies not only in dealing with a dynamical system of considerably larger size but also in handling the incompressibility constraint and the pressure field. Pressure is not governed by an explicit evolution equation but by the so-called pressure Poisson equation to adapt itself to the velocity such that the velocity field remains divergence-free. An accurate computation of the pressure field associated with an instantaneous divergence-free velocity field in a wall-bounded domain is not a trivial task,⁴⁰ let alone the derivation of the adjoint operator in the presence of this nonlocal, nonlinear operator. The second direction is developing methods following a similar idea for computing connecting orbits between invariant solutions of other types, including between two periodic orbits and eventually between invariant tori. This requires to reformulate the variational dynamics within modified search spaces appropriate for the different asymptotic boundary conditions in the temporal direction. Together with improved methods for computing other invariant solutions,^32,41,42 the proposed methodology for computing connecting orbits represents a step toward a more complete description of the state-space structures supporting spatiotemporally chaotic dynamics. Eventually, the identification of connecting orbits mediating transitions between invariant solutions may allow for efficient forecasting of chaos even in high-dimensional systems including fluid turbulence.

This research has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 865677). The authors would like to thank Sajjad Azimi and Jeremy P. Parker for their helpful discussions.

The authors have no conflicts to disclose.

Omid Ashtari: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (equal). Tobias M. Schneider: Conceptualization (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); 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.

The KSE with $L=22$ has four known equilibrium solutions including the trivial solution $E 0=0$⁠, and three nontrivial solutions $E 1$⁠, $E 2$⁠, and $E 3$ as shown in Fig. 2. The repelling eigenvalues of these equilibrium solutions are listed in Table I. The corresponding eigenvectors of $E 0$⁠, $E 1$⁠, $E 2$⁠, and $E 3$ are shown in Figs. 15–18, respectively.

In all calculations presented in Sec. V, we have used $N=64$ Fourier modes in space, have set the center of the temporal distribution at the origin $s 0=0$⁠, and have used time step size $Δτ=0.01$⁠. The temporal resolution $M$ and the scaling $S$ are listed in Table II. The temporal resolution is set high enough so that the convergence criterion $J arc , min< 10 − 12$ is achieved (see Sec. V D).