Heteroclinic and Homoclinic Connections in a Kolmogorov-Like Flow

TL;DR

This paper investigates the complex network of dynamical connections among unstable solutions in a weakly turbulent Kolmogorov flow, revealing heteroclinic and homoclinic structures that underpin turbulence.

Contribution

It provides the first extensive computation of heteroclinic and homoclinic connections in a Kolmogorov-like flow, linking various types of solutions and elucidating the chaotic dynamics.

Findings

Numerous heteroclinic connections between equilibria, periodic, and quasi-periodic orbits.

Identification of a homoclinic connection forming a chaotic repeller.

Evidence that the homoclinic tangle underpins transient turbulence.

Abstract

Recent studies suggest that unstable recurrent solutions of the Navier-Stokes equation provide new insights into dynamics of turbulent flows. In this study, we compute an extensive network of dynamical connections between such solutions in a weakly turbulent quasi-two-dimensional Kolmogorov flow that lies in the inversion-symmetric subspace. In particular, we find numerous isolated heteroclinic connections between different types of solutions -- equilibria, periodic, and quasi-periodic orbits -- as well as continua of connections forming higher-dimensional connecting manifolds. We also compute a homoclinic connection of a periodic orbit and provide strong evidence that the associated homoclinic tangle forms the chaotic repeller that underpins transient turbulence in the symmetric subspace.

Figures

Tables

Table 1. Table 1: The number of unstable directions N u subscript 𝑁 𝑢 N_{u} and the leading eigenvalue(s) (Floquet exponent(s)) λ 𝜆 \lambda for the recurrent solutions we computed in 𝒮 𝒮 \mathcal{S} at R ​ e = 22.05 𝑅 𝑒 22.05 Re=22.05 .

	$N_u$		$\lambda$
	$𝒮$	$ℳ\backslash 𝒮$	$𝒮$	$ℳ\backslash 𝒮$
EQ0	0	2	$-0.004\pm 0.020i$	$0.058\pm 0.047i$
EQ1	1	2	$0.017$	$0.029$
EQ2	2	5	$0.035\pm 0.114i$	$0.022$
EQ3	1	2	$0.068$	$0.029\pm 0.029i$
PO0	0	3	$-0.003\pm 0.059i$	0.151
PO1	1	2	0.036	$0.054\pm 0.033i$
PO2	1	3	0.044	0.049
QP1	1	-	-	-

Table 2. Table 2: Dynamical connections computed at R ​ e = 22.05 𝑅 𝑒 22.05 Re=22.05 .

	ECS-∞	ECS∞	d	$d_u^{-\infty }$	$k_s^{\infty }$
DC1, DC2	EQ2	EQ0	2	2	0
DC3	EQ2	EQ1	1	2	1
DC4, DC5	EQ1	EQ0	1	1	0
DC6–DC9	EQ2	PO1	1	2	1
DC10	EQ2	EQ3	1	2	1
DC11, DC12	EQ3	EQ0	1	1	0
DC13, DC14, DC17	PO1	EQ0	2	2	0
DC15, DC18	PO1	PO1	1	2	1
DC16, DC19	PO1	EQ3	1	2	1
DC20	PO2	EQ0	2	2	0
DC21, DC22	PO2	PO1	1	2	1
DC23	PO2	PO0	2	2	0
DC24, DC25	PO2	QP1	1	2	1
DC26	QP1	PO0	3	3	0

Equations

$${\partial }_t\mathbf{u}+\beta \mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\frac{1}{Re}\left({\nabla }^2\mathbf{u}-\gamma \mathbf{u}\right)+\mathbf{f},$$

$$\mathbf{u}(\theta ,0)={\mathbf{u}}_{eq}+ϵ\cos (\theta ){\hat{\mathbf{e}}}_1^{\prime }+ϵ\sin (\theta ){\hat{\mathbf{e}}}_2^{\prime }.$$

$$\mathbf{u}(\theta ,t)\approx {\mathbf{u}}_{eq}+ϵe^{\sigma t}\cos (\theta +\mu t){\hat{\mathbf{e}}}_1^{\prime }+ϵe^{\sigma t}\sin (\theta +\mu t){\hat{\mathbf{e}}}_2^{\prime }.$$

$$c_k(t)=\frac{(\mathbf{u}(\theta ,t)-{\mathbf{u}}_{eq})\cdot {\hat{\mathbf{e}}}_k^{\prime }}{D_c},$$

$$s(\theta ,t)=\frac{{\tau }_c}{\Vert \mathbf{u}(\theta ,t)\Vert }\Vert {\partial }_t\mathbf{u}(\theta ,t)\Vert .$$

$$D_0=\underset{t}{\min }\frac{\Vert \mathbf{u}(\theta ,t)-{\mathbf{u}}_{eq}\Vert }{D_c},$$

$${\mathbf{u}}^{\pm }(t)={\mathbf{u}}_{eq}\pm ϵe^{{\lambda }_{1}t}{\hat{\mathbf{e}}}_1,$$

$${\mathbf{u}}^{\pm }(\eta ,0)={\mathbf{u}}_{po}\pm ϵe^{{\lambda }_1\eta }{\hat{\mathbf{e}}}_1,$$

$$D_1=\underset{t,\eta }{\min }\frac{\Vert \mathbf{u}(\theta ,t)-{\mathbf{u}}_{po}(\eta )\Vert }{D_c},$$

$$\mathbf{u}(t)={\mathbf{u}}_{ecs}+\sum _i{a}_i(t){\hat{\mathbf{e}}}_i.$$

$$a_i(t)={\hat{\mathbf{e}}}_i^{†}\cdot (\mathbf{u}(t)-{\mathbf{u}}_{ecs}),$$

$${\hat{\mathbf{e}}}_i=\sum _j{T}_{ij}{\hat{\mathbf{e}}}_j^{\prime },$$

$$c_j=\sum _i{a}_i{T}_{ij}/D_c$$

$$r(t)=\underset{\tau \in [{\tau }_l,{\tau }_h]}{\min }\frac{\Vert \mathbf{u}(t)-\mathbf{u}(t+\tau )\Vert }{D_c},$$

Full text

Dynamical systems, Exact Coherent Structures, Invariant Solutions, Invariant manifolds

I Introduction

Turbulent fluid flows are ubiquitous; they can be found in the atmosphere and the oceans, water and oil pipelines, and even in the human aorta. Despite its great practical relevance, a tractable description of turbulent dynamics has remained elusive. However, recent numerical Kerswell (2005); Kawahara et al. (2012) and experimental Hof et al. (2004); de Lozar et al. (2012); Suri et al. (2017) studies have shown that unstable recurrent solutions of the Navier Stokes equation, which governs fluid flows, may prove pivotal in solving this longstanding problem. Often termed Exact Coherent States (ECSs), such solutions exist for the same parameters as turbulence but are more amenable to numerical analysis given their simple (e.g., steady or periodic) temporal behavior.

The state space description of turbulence best illustrates the dynamical role of ECSs Hopf (1948); Gibson et al. (2008). A turbulent flow in physical space maps to a winding trajectory in state space, with each point on it representing a flow field (see supplementary videos 1-7) 111. In contrast, ECSs such as steady and time-periodic flows are simpler objects (fixed points, closed loops), as shown in Fig. 1. Being unstable, each ECS is a saddle in state space; trajectories in its unstable manifold are repelled away, while those in the stable manifold converge to the ECS Duguet et al. (2008a); Viswanath and Cvitanovic (2009); Suri et al. (2017, 2018). Heteroclinic (homoclinic) trajectories – which originate in the unstable manifold of an ECS and terminate in the stable manifold of another (the same) ECS – connect different ECSs and create a chaotic saddle in state space Halcrow et al. (2009); Duguet et al. (2008a); Farano et al. (2018). Dynamics on ECSs, their stable manifolds, and on homo/heteroclinic connections are asymptotically non-chaotic. In this geometrical picture, turbulence represents a deterministic walk between neighborhoods of different ECSs, guided by the corresponding dynamical connections Kerswell and Tutty (2007); Suri et al. (2017), both homoclinic and heteroclinic.

Substantial numerical evidence has emerged for the dynamical relevance of ECSs in recent years, mostly from research on three-dimensional (3D) wall-bounded shear flows, such as plane-Couette Kawahara and Kida (2001); Kawahara (2005); Viswanath (2007); Gibson et al. (2009), pipe Faisst and Eckhardt (2003); Wedin and Kerswell (2004); Kerswell and Tutty (2007); Pringle and Kerswell (2007), and channel flows Waleffe (2001); Itano and Toh (2001); Toh and Itano (2003). Direct numerical simulation (DNS) of flows in small, spatially periodic domains Kim et al. (1987); Hamilton et al. (1995) suggests that turbulence at moderate Reynolds numbers ($Re$) is organized around unstable solutions such as equilibria (EQ) Wang et al. (2007); Gibson et al. (2008, 2009), traveling waves (TW) Waleffe (2001); Faisst and Eckhardt (2003); Wedin and Kerswell (2004); Pringle and Kerswell (2007), periodic orbits (PO) Kawahara and Kida (2001); Viswanath (2007); Duguet et al. (2008b), and relative periodic orbits (RPO) Cvitanović and Gibson (2010); Budanur et al. (2017). Here, TWs (RPOs) are solutions that correspond to steady (time-periodic) states in a reference frame moving in the direction of a continuous symmetry (e.g., along the axis of a pipe). Flow fields resembling ECSs were also observed in a few laboratory experiments Hof et al. (2004); de Lozar et al. (2012); Lemoult et al. (2014); Suri et al. (2017), which further validated their relevance in turbulence.

The geometry of chaotic saddle shaped by invariant (stable/unstable) manifolds and dynamical connections between ECSs, however, remains under-explored. In particular, the connectivity of different neighborhoods can be determined by generating a dense set of trajectories spanning the unstable manifold of an ECS and identifying which neighborhoods are subsequently visited by each trajectory Kerswell and Tutty (2007); Gibson et al. (2008). In general, dynamically relevant ECSs have several (three of more) unstable directions Gibson et al. (2009); Kerswell and Tutty (2007); Budanur et al. (2017); Farano et al. (2018); Suri et al. (2018), which renders this procedure computationally expensive. To circumvent this challenge numerical studies to date have analyzed dynamics confined to invariant subspaces (e.g., symmetry subspaces, laminar-turbulent boundary) which reduces the number of dynamically relevant ECSs as well as the dimensionality of their unstable manifolds. Using this technique, several previous studies Kawahara and Kida (2001); Toh and Itano (2003); Kerswell and Tutty (2007); Duguet et al. (2008b, a); Budanur and Hof (2017) identified trajectories that originate at ECSs with only one or two unstable directions and subsequently approach another ECS. However, such trajectories were not proven to asymptotically converge to an ECS. Hence, they can only be regarded as likely signatures of dynamical connections.

Even within invariant subspaces, very few dynamical connections between ECSs have been found. Gibson *et al.*aGibson et al. (2008) and Halcrow *et al.*aHalcrow et al. (2009) computed four heteroclinic connections between unstable EQs in plane-Couette flow (PCF). Using low-dimensional projections of state space, the authors showed that turbulent trajectories are transiently guided by these connections. The structural stability of these connections – their robustness to small changes in $Re$ – was also discussed using dimension counting arguments Smale (1967). Two homoclinic orbits of a PO in PCF were computed by van Veen *et al.*avan Veen and Kawahara (2011) using a multi-shooting algorithm van Veen et al. (2011). The authors suggested that dynamics along these connections resemble “bursting” phenomenon observed in turbulent boundary layers Kline et al. (1967). Riols *et al.*aRiols et al. (2013), in a study of subcritical transition to magnetorotational dynamo in Keplerian shear flows, computed both homoclinic and heteroclinic connections between unstable POs. Pershin *et al.*aPershin et al. (2019) found a heteroclinic connection from an EQ to a nearby PO in PCF. In both Riols *et al.*aand Pershin et al., connections were computed very close to the saddle-node bifurcations leading to the formation of these invariant solutions. Recently, Budanur *et al.*aBurak Budanur et al. (2018) computed a homoclinic orbit of a spatially localized RPO in pipe flow and suggested that the associated homoclinic tangle underlies transient turbulence in this flow. All the connections listed above originate (terminate) at ECSs with only one/two (one/zero) unstable directions, which facilitates the use of simple shooting and bisection algorithms. Recently, Farano *et al.*aFarano et al. (2018) showed that an adjoint-based method can be used to find connections between neighborhoods of unstable EQs in PCF where the originating (destination) EQs has more than two (one) unstable direction.

Despite these advances, an extensive exploration of state space to detect signatures of connections between dynamically dominant ECSs has not been carried out yet. For instance, some previous studies computed connections between ECSs of the same type (i.e., between EQs or between POs) Gibson et al. (2008); Halcrow et al. (2009); van Veen and Kawahara (2011) while others reported only isolated connections between POs and EQs Riols et al. (2013); Pershin et al. (2019). The chaotic saddle, however, is shaped by ECSs of different types and a complex network of dynamical connections between them. To address this shortcoming of previous studies, we report in this article a systematic and exhaustive exploration of low-dimensional unstable manifolds of ECSs to detect dynamical connections between various types on ECSs.

The system we numerically study is the quasi-two-dimensional (Q2D) Kolmogorov-like flow in a shallow electrolyte layer driven by a horizontal, (spatially) near-sinusoidal body force. Q2D flows are computationally more tractable than 3D flows since they can be accurately described using a 2D model Bondarenko et al. (1979); Suri et al. (2014). In fact, the relative simplicity of 2D DNS has prompted researchers in recent years to carry out the most systematic exploration of ECSs in 2D turbulence Chandler and Kerswell (2013); Lucas and Kerswell (2014); Farazmand (2016). For instance, Chandler *et al.*aChandler and Kerswell (2013) studied 2D Kolmogorov flow on a periodic domain to test whether turbulent statistics can be reproduced using suitable averages over time-periodic solutions Cvitanović (1988). More recently, Suri *et al.*aSuri et al. (2017, 2018) have validated the dynamical role of EQs and their unstable manifolds, for the first time in laboratory experiments, using Q2D Kolmogorov-like flow as the test bed. 2D DNS in these studies Suri et al. (2017, 2018) were performed on a numerical domain with lateral dimensions and boundary conditions identical to those in the experiment, facilitating quantitative comparison between DNS and experiment.

This article is structured as follows: In Sect. II we discuss the 2D model for Q2D flows and its numerical implementation. In Sections III.1 – III.4 we explore invariant manifolds of various EQs and POs and identify heteroclinic and homoclinic connections. The stability of dynamical connections to small changes in Reynolds number is discussed in Sect. III.5. In Sect. III.6 we discuss the relation between the dynamical connections and transient turbulence. Finally, we summarize our findings and discuss their significance in Sect. IV.

II Quasi-2D Kolmogorov-like Flow

The evolution of weakly turbulent flows in an electromagnetically driven shallow electrolyte layer is modeled here using a strictly 2D equation (Suri et al., 2014)

[TABLE]

which is derived by averaging 3D Navier Stokes equation in the vertical ($z$) direction Suri et al. (2014). Here, $\nabla$ is the gradient operator restricted to the horizontal dimensions $x$, $y$, and ${\bf u}(x,y,t)$ is the velocity field at the electrolyte-air interface which is assumed to be incompressible ($\nabla\cdot{\bf u}=0$). This is an accurate approximation for moderate Reynolds numbers ($Re\lesssim 40$) Suri et al. (2014). ${\bf f}$ is the nondimensionalized depth-averaged horizontal forcing profile, and $p$ plays the role of the 2D kinematic pressure. In Q2D flows, the solid boundary beneath the fluid layer causes a vertical gradient in the magnitude of the horizontal velocity. The prefactor $\beta\neq 1$ and linear friction $-\gamma{\bf u}$ model the effective change in inertia and the shear stress, respectively, due to this vertical gradient Bondarenko et al. (1979); Suri et al. (2014). For the flow in experiments detailed in Suri *et al.*aSuri et al. (2017), which we numerically study in this article, $\beta=0.83$ and $\gamma=3.22$. We note that these values differ significantly from $\beta=1$ and $\gamma=0$ corresponding to an unphysical strictly 2D flow typically studied in numerics Chandler and Kerswell (2013); Farazmand (2016). The Reynolds number $Re$ describes the strength of electromagnetic forcing and serves as the parameter that controls the complexity of flow. The dimensional form of equation (1) and analytical expressions for $\beta$, $\gamma$, and $Re$ as functions of experimental parameters are provided in references Suri et al. (2014) and Tithof et al. (2017).

DNS of the flow was performed using a finite-difference code previously employed in studies Suri et al. (2017, 2018) and Tithof et al. (2017). Velocity and pressure fields on a computational domain with lateral dimensions $(L_{x},L_{y})=(14,18)$ were spatially discretized using a 280$\times$360 staggered grid with spacing $\delta x=\delta y=0.05$ between grid points. No-slip boundary conditions were imposed on the velocity field, and spatial derivatives were approximated using second-order central finite difference formulas. Temporal integration of equation (1) was performed using the semi-implicit P2 projection scheme to enforce incompressibility of the velocity field at each time step (Armfield and Street, 1999). The linear (nonlinear) terms in equation (1) were discretized in time using second order implicit Crank-Nicolson (explicit Adams-Bashforth) method. A time step $\delta t\leq 1/100$ was used for temporal integration to ensure the CFL number ${\max\{{u_{x}}\delta t/\delta x,u_{y}\delta t/\delta y\}}\leq 0.5$.

In a 2D Kolmogorov flow, the forcing profile is strictly sinusoidal, i.e., ${\bf f}\propto\sin(\pi y){\,\bf\hat{x}}$. In experiments detailed in references (Suri et al., 2017; Tithof et al., 2017), however, the electromagnetic forcing is sinusoidal only near the center of the domain and decays to zero at the boundaries. To replicate this forcing, we used a dipole lattice approximation of the magnet array in the experiment and computed the resulting electromagnetic forcing. Comparison between experimentally measured forcing profile and the numerically estimated one was provided in Tithof *et al.*aTithof et al. (2017). The 2D forcing profile ${\bf f}$ computed from the dipole lattice model is anti-symmetric under the coordinate transformation $\mathcal{I}:(x,y)\rightarrow(-x,-y)$, i.e., $\mathcal{I}{\bf f}=-{\bf f}$. This 2-fold inversion symmetry ($\mathcal{I}^{2}=1$) is equivalent to rotation by $\pi$ about the $z-$axis passing through the lateral center of the computational domain. Under $\mathcal{I}$, the velocity field transforms as $\mathcal{I}{\bf u}(x,y,t)\rightarrow-{\bf u}(-x,-y,t)$, which makes equation (1) equivariant under ${\mathcal{I}}$.

A consequence of this equivariance is that the symmetry of a rotationally invariant flow is preserved during its time evolution, i.e., the symmetry subspace $\mathcal{S}=\{{\bf u}\in\mathcal{M}\,|\,\mathcal{I}{\bf u}={\bf u}\}$ of $\mathcal{M}$ is invariant. Here, $\mathcal{M}$ represents the full state space. All ECSs and connections between them we report in this article lie in $\mathcal{S}$. Since trajectories in $\mathcal{S}$ are generally unstable, numerical errors will accumulate such that ${\bf u}(t)$ will eventually leave $\mathcal{S}$ even if ${\bf u}(0)\in\mathcal{S}$. To prevent this, we augmented the numerical integrator by projecting ${\bf u}(t)$ back into $\mathcal{S}$ after every time step Suri (2017).

III Results

As $Re$ is increased, the Kolmogorov-like flow undergoes a number of bifurcations before transitioning to weak turbulence at $Re\approx 18$ Tithof et al. (2017). At $Re\approx 22$ the flow in full state space $\mathcal{M}$ is chaotic, which was confirmed by computing the Lyapunov spectrum of a long turbulent trajectory using continuous Gram-Schmidt orthogonalization Farmer et al. (1983); Suri et al. (2018). In this dynamical regime we previously identified 31 unstable equilibria of the 2D model Suri et al. (2018). Twenty eight of these were outside of $\mathcal{S}$. Furthermore, their unstable manifolds were relatively high-dimensional in $\mathcal{M}$ (with the number of unstable directions as low as three and as high as twenty) and thus computationally forbidding to map out. Hence, we instead started our analysis with the equilibria in $\mathcal{S}$, labeled EQ0, EQ1, and EQ2.

The stability properties of these equilibria and other recurrent solutions we found in $\mathcal{S}$ are summarized in Table 1. We have listed the number of unstable directions $N_{u}$ and the leading eigenvalues $\lambda$ for EQs (Floquet exponents for POs) both in the symmetry invariant subspace $\mathcal{S}$ and in the rest of state space, i.e., $\mathcal{M}\backslash\mathcal{S}$. All ECSs we computed have two or fewer unstable directions in $\mathcal{S}$, which allows their unstable manifolds to be well-approximated with a dense set of trajectories Suri et al. (2018). Note that EQ0 and PO0 are stable in $\mathcal{S}$. Hence, turbulence in $\mathcal{S}$ can be transient, with trajectories converging to one of these two ECSs asymptotically in time. This is indeed what numerical simulations show, as we discuss in Sect. III.6.

Connection(s) from an origin (ECS*-∞) to a destination (ECS∞) lie at the intersection of the unstable manifold of ECS-∞* and the stable manifold of ECS*∞*. Since the dimensionality of stable manifolds is very high ($O(10^{5})$ in our case), we will search the low-dimensional unstable manifolds of ECSs for dynamical connections. Such unstable manifolds can be conveniently approximated by families of trajectories ${\bf u}({\bf p},t)$ that start near an ECS, where the unstable manifold is locally tangent to the linear subspace parametrized by the unstable and marginal eigenvectors of that ECS. Here, ${\bf p}$ parametrizes the family of manifold trajectories. The specifics of constructing a linear subspace, and consequently its parametrization, depend on the type of ECS (EQ or PO) as well as the dimensionality of its unstable manifold, as discussed below.

Table 1 shows that several solutions we identified have a single unstable direction in $\mathcal{S}$. The unstable manifold of such ECS is naturally divided into two halves, which correspond to the positive and negative halves of the corresponding linear subspace (see Fig. 1). We shall hereafter refer to these halves as ${\bf u}^{\pm}({\bf p},t)$. Furthermore, ${\bf u}^{\pm}({\bf p},t)$ lie on the opposite sides of the stable manifold of the ECS, which serves as a local separatrix and repels trajectories near the ECS along ${\bf u}^{+}({\bf p},t)$ or ${\bf u}^{-}({\bf p},t)$. This allows us to compute connections terminating at ECSs with one unstable direction using a simple bisection method Itano and Toh (2001). In the following sections we provide a detailed discussion of connections between various ECSs listed in Table 1.

III.1 Connections originating at EQ2

We start our analysis with equilibrium EQ2, which has the highest dimensional unstable manifold of all the EQs in $\mathcal{S}$. The vorticity field ($\omega=\nabla\times{\bf u}$) associated with EQ2 is shown in Fig. 2(a). The 2D unstable manifold of EQ2 is locally tangent to the plane spanned by the complex conjugate pair of unstable eigenvectors ${\,\bf\hat{e}}_{1},{{\,\bf\hat{e}}}_{2}={\,\bf\hat{e}}_{1}^{\ast}$ associated with $\lambda_{1,2}=\sigma\pm i\mu$, where $\sigma=0.035$ and $\mu=0.114$. To construct this 2D surface, which lies entirely in $\mathcal{S}$, we generated 360 initial conditions ${\bf u}(\theta,0)$ on a circle around EQ2 in the plane spanned by ${{\,\bf\hat{e}}}_{1}$, ${{\,\bf\hat{e}}}_{2}$:

[TABLE]

Here, ${\bf u}_{eq}$ corresponds to EQ2 and ${\,\bf\hat{e}}_{1}^{\prime}$, ${\,\bf\hat{e}}_{2}^{\prime}$ are real orthonormal vectors constructed from the real and imaginary parts of ${{\,\bf\hat{e}}}_{1}$, ${\,\bf\hat{e}}_{2}$. For numerical convenience we chose $\epsilon=0.002\cdot\|{\bf u}_{eq}\|$ and positioned the initial conditions ${\bf u}(\theta,0)$ at equal angular intervals $\Delta\theta=2\pi/360$ on the circle. Numerical integration of initial conditions ${\bf u}(\theta,0)$ generates trajectories ${\bf u}(\theta,t)$ that approximate the 2D manifold and parametrize it by $\theta$, $t$. Each manifold trajectory ${\bf u}(\theta,t)$ was computed on a temporal interval of length 25$\tau_{c}$, where $\tau_{c}=12.5$ (nondimensional time units) is the average temporal auto-correlation time Suri et al. (2018).

In the neighborhood of EQ2, the various trajectories ${\bf u}(\theta,t)$ initially evolve spiraling outward, i.e.,

[TABLE]

To illustrate this, we plotted a portion of the 2D unstable manifold in Fig. 3, projected onto vectors ${\,\bf\hat{e}}_{1}^{\prime}$ and ${\,\bf\hat{e}}_{2}^{\prime}$. Only one in every ten trajectories generated is shown and the segment of each trajectory ${\bf u}(\theta,t)$ plotted corresponds to $0\leq t\leq 11\tau_{c}$. Coordinates $c_{1},c_{2}$ in Fig. 3 are the normalized inner products

[TABLE]

where, $D_{c}=\max_{t,t^{\prime}}\|{\bf u}(t)-{\bf u}(t^{\prime})\|=300$ is the empirically estimated maximum separation between two points on an 800$\tau_{c}$-long turbulent trajectory in $\mathcal{S}$, which defines the “diameter” of the chaotic repeller. Normalizing distances with $D_{c}$ ensures that points separated by nearly unit distance in the low-dimensional projection are very far apart in full state space. Notice that, farther away from EQ2 ($|c_{1}|,|c_{2}|>0.2$) the shape of ${\bf u}(\theta,t)$ becomes fairly complicated due to the nonlinearity of the governing equation (1).

It is not known a priori which ${\bf u}(\theta,t)$ approach an ECS (EQ or PO) closely after leaving the neighborhood of EQ2. Previously, Gibson *et al.*aGibson et al. (2008) and Halcrow *et al.*aHalcrow et al. (2009) inferred possible dynamical connections by inspecting low-dimensional state space projections. In contrast, Riols *et al.*aRiols et al. (2013) analyzed time-series of (magnetic) energy and identified close passes to POs using intervals with “periodic” behavior. Both these techniques, however, require visual inspection and detecting signatures of dynamical connections cannot be automated. In this study we tested two methods to detect signatures of dynamical connections that do not require laborious visual inspection. The method employed for detecting connections that terminate at EQs proved very effective and is discussed next. A method for detecting connections to POs using recurrence analysis Duguet et al. (2008a) is discussed in Appendix B.

For each trajectory ${\bf u}(\theta,t)$ in the unstable manifold of EQ2, we computed the normalized instantaneous state space speed $s(\theta,t)$, which we define as Suri et al. (2017, 2018):

[TABLE]

Since $\partial_{t}{\bf u}=0$ ($s=0$) for any EQ, $s(\theta,t)\ll 1$ indicates that a trajectory ${\bf u}(\theta,t)$ in state space is near an EQ Suri et al. (2017); Neelavara et al. (2017). In the physical space, this corresponds to the evolution of the flow dramatically slowing down. Figure 4 shows, as examples, the speed plots for two different manifold trajectories ${\bf u}(\pi/2,t)$ and ${\bf u}(3\pi/2,t)$. Clearly, ${\bf u}(\theta,t)$ lie in the linear neighborhood of EQ2 for $0<t\lesssim 10\tau_{c}$ ($s\ll 1$) and $s(\theta,t)$ grows exponentially according to equation (3). After this initial transient, however, dynamics described by various trajectories can be qualitatively very different. For instance, the shape of $s(3\pi/2,t)$ (gray curve) suggests that ${\bf u}(3\pi/2,t)$ displays turbulent evolution. In contrast, $s(\pi/2,t)$ (black curve) suggests that ${\bf u}({\pi/2},t)$ approaches an EQ closely after a brief turbulent excursion.

Hence, to test if a trajectory ${\bf u}(\theta,t)$ approaches an EQ after it departs from the neighborhood of EQ2, we computed the minimum speed $s_{m}(\theta)=\min_{t}s(\theta,t)$ for $10\tau_{c}\leq t\leq 25\tau_{c}$, which is shown in Fig. 5. The low values of $s_{m}(\theta)$ for $\theta\in\left(\theta_{1},\theta_{2}\right)$, $\theta\in\left(\theta_{3},\theta_{4}\right)$, and $\theta={\theta_{5}}$ strongly suggest that the corresponding trajectories ${\bf u}(\theta,t)$ closely shadow either heteroclinic or homoclinic connections from EQ2 to other EQs. In Fig. 5, the uncertainty in $\theta$ is limited by the angular resolution $\Delta\theta$ in the initial conditions ${\bf u}(\theta,0)$. Consequently, to compute exact dynamical connections, the estimates for $\theta$ may require refinement in some cases. The labels $\theta_{1}$ through $\theta_{5}$ in Fig. 5 correspond to these refined values. Lastly, $s_{m}(\theta)>0$ even for trajectories ${\bf u}(\theta,t)$ that correspond exactly to dynamical connections. This is because, minima in the $s(\theta,t)$ are computed on a finite temporal interval, while dynamical connections converge ($s_{m}=0$) to the destination EQs only in infinite time.

The destination EQs that different dynamical connections approach were computed using a Newton solver Saad and Schultz (1986); Kelley (2003); Mitchell (2013) initialized with the flow field ${\bf u}(\theta,t)$ at the instant when $s(\theta,t)=s_{m}(\theta)$. For $\theta\in\left(\theta_{1},\theta_{2}\right)$ and $\theta\in(\theta_{3},\theta_{4})$ the solver converged to EQ0, which is a stable node in $\mathcal{S}$. For $\theta=\theta_{5}$ the solver converged to EQ1, which is a saddle with one unstable direction in $\mathcal{S}$. The flow fields corresponding to EQ0 and EQ1 are shown in Figs. 2(b) and 2(c). Dimensional arguments (discussed in Sect. III.5) suggest that connections from EQ2 to EQ0 should comprise a one-parameter family forming a two-dimensional manifold, just like those in the interval $(\theta_{3},\theta_{4})$. However, the difference $\theta_{2}-\theta_{1}$ is smaller than the resolution $\Delta\theta$ in case of the interval $(\theta_{1},\theta_{2})$. Consequently, $\theta_{1}$, $\theta_{2}$ cannot be distinguished in Fig. 5. Hereafter, we will refer to $\theta$ in this narrow interval $(\theta_{1},\theta_{2})$ collectively as $\theta_{n}$. The wide interval $(\theta_{3},\theta_{4})$ will be referred to as $\theta_{w}$. Lastly, to confirm that the dynamical connections originating at EQ2 indeed terminate at either EQ0 or EQ1, we need to make sure that the distance

[TABLE]

vanishes in the limit $t\to\infty$. Here, ${\bf u}_{eq}$ is the velocity field corresponding to the destination EQ. For all connections reported in this study, we chose $D_{0}<0.005$ as the criterion for convergence.

Let us start our analysis with the trajectories ${\bf u}(\theta_{n},t)$ and ${\bf u}(\theta_{w},t)$ that approach the equilibrium EQ0 which is stable in $\mathcal{S}$. The eigenvalues of weakly contracting modes of EQ0 are $\lambda_{1,2}=-0.004\pm 0.020i$, $\lambda_{3,4}=-0.012\pm 0.42i$. Consequently, most nearby trajectories in $\mathcal{S}$ converge to EQ0 at a rate determined by the real part of $\lambda_{1}$ (${\lambda_{1}}$). Hence, we computed trajectories ${\bf u}(\theta_{n},t)$ and ${\bf u}(\theta_{w},t)$ for an interval of duration $350\tau_{c}\gg 1/\lambda_{1}\approx 10\tau_{c}$. At the end of numerical integration both ${\bf u}(\theta_{n},t)$ and ${\bf u}(\theta_{w},t)$ approached EQ0 to within a distance $D_{0}<10^{-8}$, confirming that we have indeed found heteroclinic connections from EQ2 to EQ0. Hereafter, we shall refer to these dynamical (heteroclinic) connections corresponding to ${\bf u}(\theta_{n},t)$ and ${\bf u}(\theta_{w},t)$ as DC1 and DC2, respectively.

A connection from the family DC1 is shown (gray curve) in Fig. 6(a) and thirty connections from DC2 that are equally spaced on the interval $(\theta_{3},\theta_{4})$ are shown in Fig. 6(b). Both figures show the projection of state space onto an orthogonal basis constructed from the stable eigenvectors ${\,\bf\hat{e}}_{1}$, ${\,\bf\hat{e}}_{2}$, and ${\,\bf\hat{e}}_{5}$ of EQ0. The vectors ${\,\bf\hat{e}}_{1}$, ${\,\bf\hat{e}}_{2}$ were chosen to visualize asymptotic dynamics along the connections terminating at EQ0. The vector ${\,\bf\hat{e}}_{5}$ was chosen because trajectories ${\bf u}(\theta,t)$ far away from EQ0 have large components along this direction. We note that, unless mentioned otherwise, these modes are used for all state space projections in this article.

In Fig. 6, both DC1 and DC2 converge to EQ0 spiraling inwards, almost entirely confined to the $c_{1}-c_{2}$ plane. This is a consequence of the large separation between the real values of $\lambda_{1,2}$ and $\lambda_{3,4}$ which makes ${\,\bf\hat{e}}_{1}$, ${\,\bf\hat{e}}_{2}$ the only dynamically relevant eigenvector pair near EQ0. Another interesting feature of Fig. 6(b) is that DC2 initially forms a flat strip bounded by manifold trajectories ${\bf u}(\theta_{3},t)$ and ${\bf u}(\theta_{4},t)$ (black curves). However, farther away from EQ2, this strip widens and folds such that ${\bf u}(\theta_{3},t)$ and ${\bf u}(\theta_{4},t)$ come close to each other, trace loops which are strikingly similar in shape, and eventually merge (see inset in Fig. 8). This is not an artifact of low-dimensional projection of the state space. As we explain in the following section, the black trajectories in Fig. 6(b) correspond to heteroclinic connections from EQ2 to a different ECS (PO1).

Unlike DC1 and DC2, the trajectory ${\bf u}(\theta_{5},t)$ approaches the solution EQ1 which has one unstable direction (cf. Table 1). Since the stable/unstable manifolds of EQ1 guide the evolution of trajectories in its neighborhood, we explored their geometry to compute the connection from EQ2 to EQ1. Trajectories ${\bf u}(\theta,t)$ for $\theta\approx\theta_{5}$ that approach EQ1 should subsequently depart following its unstable manifold, which coincides with a pair of trajectories

[TABLE]

in the linear neighborhood of EQ1 (cf. Fig. 1). Here, ${\,\bf\hat{e}}_{1}$ is the unstable eigenvector of EQ1 and $\epsilon>0$ is some small constant. In particular, we found a pair of adjacent trajectories ${\bf u}(\theta,t)$ and ${\bf u}(\theta+\Delta\theta,t)$ with $\theta\approx\theta_{5}$ approach EQ1 and depart its neighborhood in opposite directions shadowing ${\bf u}^{+}(t)$ and ${\bf u}^{-}(t)$, respectively. Hence, a heteroclinic connection DC3 from EQ2 to EQ1 lies between these two trajectories, i.e., in the stable manifold of EQ1.

The projection of DC3 (black curve) is shown in Fig. 6(a). It was computed by refining the estimate for $\theta_{5}$ to within $\delta\theta=2^{-12}\Delta\theta$ using bisection Itano and Toh (2001); Gibson et al. (2008); Halcrow et al. (2009). This refinement reduced the separation between ${\bf u}(\theta_{5},t)$ and EQ1 to $D_{0}=0.002$ (cf. equation (6)). Further refinement in $\theta$ did not significantly reduce $D_{0}$ since numerical noise on ${\bf u}(\theta_{5},t)$ amplifies rapidly in the direction of unstable eigenvector $\pm{\,\bf\hat{e}}_{1}$ of EQ1. This behavior stems from the strong asymmetry in the real parts of unstable ($\lambda_{1}=0.017$) and weakly stable ($\lambda_{2,3}=-0.003\pm 0.109i$) eigenvalues of EQ1. Nevertheless, a better accuracy can be achieved by employing multi-shooting van Veen et al. (2011) or approximating DC3 as a piece-wise continuous solution Toh and Itano (2003). Using the latter technique allowed us to decrease $D_{0}$ to less than $10^{-6}$. Lastly, equation (7) governs the evolution of trajectories in the unstable manifold of EQ1 only in its linear neighborhood. Farther away, both ${\bf u}^{\pm}(t)$ display turbulent evolution for $t\gg 100\tau_{c}$ and eventually approach EQ0. Hence these two trajectories define (long) heteroclinic connections (DC4, DC5) from EQ1 to EQ0.

III.2 Dynamics Near PO1 and EQ3

While trajectories ${\bf u}(\theta,t)$ for $\theta\in(\theta_{3},\theta_{4})$ quickly converge to EQ0, trajectories just outside this interval exhibit qualitatively different dynamics; this can be inferred from $s_{m}$ changing abruptly at $\theta_{3}$ and $\theta_{4}$. This qualitative difference suggests that trajectories ${\bf u}(\theta_{3},t)$ and ${\bf u}(\theta_{4},t)$ play the role of separatrices on the unstable manifold of EQ2 and hence lie in the stable manifold of another ECS. In case of ${\bf u}(\theta_{4},t)$, for example, we found this ECS by inspecting the two trajectories ${\bf u}(\theta_{4}^{-},t)$ and ${\bf u}(\theta_{4}^{+},t)$ obtained as a result of successive bisections; here $\theta_{4}^{\pm}=\theta_{4}\pm\delta\theta$ with $\delta\theta=2^{-12}\Delta\theta$. The corresponding state space speed plots are shown in Fig. 7, which suggest that ${\bf u}(\theta_{4}^{\pm},t)$ evolve almost indistinguishably for about $55\tau_{c}$ and subsequently separate. For $t/\tau_{c}\in[30,55]$ the state space speed oscillates with a period of approximately $1.4\tau_{c}$, suggesting that ${\bf u}(\theta_{4}^{\pm},t)$ shadow a periodic orbit during this interval.

Using the flow field corresponding to ${\bf u}(\theta_{4}^{+},t)$ at $t/\tau_{c}=55$ as an initial condition into a Newton solver, we indeed found an unstable periodic orbit PO1 with a period $T=1.32\tau_{c}$. A similar refinement using bisection showed that trajectories at $\theta_{3}^{\pm}=\theta_{3}\pm\delta\theta$ also approach PO1 and separate after shadowing it for an extended period of time (the corresponding speed plots are not shown). Figure 8 shows the state space projection of ${\bf u}(\theta_{4}^{\pm},t)$ and ${\bf u}(\theta_{3}^{\pm},t)$ approaching PO1, shadowing it closely (inset), and subsequently leaving its neighborhood. The projection coordinates here are the same as in Fig. 6, but the viewing angle is different. Lastly, the result that ${\bf u}(\theta_{3}^{\pm},t)$ and ${\bf u}(\theta_{4}^{\pm},t)$ approach PO1 is consistent with the folding of DC2 shown in Fig. 6(b).

PO1 has a single real unstable direction in $\mathcal{S}$ with an associated Floquet exponent $\lambda_{1}=0.036$ (cf. Table 1). Hence, the departure of ${\bf u}(\theta_{3}^{\pm},t)$ and ${\bf u}(\theta_{4}^{\pm},t)$ from the neighborhood of PO1 is guided by its unstable manifold, which is two-dimensional since it is associated with a PO which is itself one-dimensional. This 2D unstable manifold can be constructed by evolving initial conditions Riols et al. (2013)

[TABLE]

where ${\bf u}_{po}$ is a reference point on PO1, ${\,\bf\hat{e}}_{1}$ is the Floquet vector at ${\bf u}_{po}$, $\epsilon>0$ is sufficiently small (we chose $\epsilon=0.001\cdot\|{\bf u}_{po}\|$), and $\eta\in(0,T]$ parametrizes different initial conditions. A total of $720$ initial conditions, equally spaced in $\eta$, were generated and the corresponding trajectories ${\bf u}^{\pm}(\eta,t)$ were computed on an interval of length $35\tau_{c}$ to approximate the manifold. We found that trajectories ${\bf u}^{+}(\eta,t)$ uneventfully converge to EQ0 while ${\bf u}^{-}(\eta,t)$ display turbulent evolution after leaving the neighborhood of PO1. To illustrate this, we plotted in Fig. 8 a pair of trajectories ${\bf u}^{\pm}(\eta,t)$ (red dashed curves) evolving in opposite directions.

Since PO1 has a single unstable direction, its stable manifold divides the state space in the neighborhood of PO1 into two halves Avila et al. (2013). We found that trajectories ${\bf u}(\theta_{3}^{+},t)$ and ${\bf u}(\theta_{4}^{-},t)$ smoothly converge to EQ0 after approaching PO1 and hence lie on one (the same) side of the stable manifold. Meanwhile ${\bf u}(\theta_{3}^{-},t)$ and ${\bf u}(\theta_{4}^{+},t)$ display turbulent excursions and hence lie on the opposite side of the stable manifold, as shown in Fig. 8. Therefore, there exist two heteroclinic connections from EQ2 to PO1 sandwiched between ${\bf u}(\theta_{3}^{\pm},t)$ (DC6) and ${\bf u}(\theta_{4}^{\pm},t)$ (DC7), that lie on the stable manifold and asymptotically converge to PO1. Just like in the case of DC3, simple shooting cannot be employed to compute DC6 and DC7 in their entirety since they are very unstable. Nevertheless, for sufficient refinement in $\theta$, the trajectories ${\bf u}(\theta_{3}^{\pm},t)$ and ${\bf u}(\theta_{4}^{\pm},t)$ approximate DC6 and DC7, respectively, reasonably accurately (cf. Fig. 6(b)). This was tested by computing the smallest distance between ${\bf u}(\theta,t)$ and PO1

[TABLE]

where ${\theta}=\theta_{3}^{\pm}$ or $\theta_{4}^{\pm}$ and ${\bf u}_{po}(\eta)$ are states along PO1 parametrized by $\eta\in(0,T]$. We found that $D_{1}<0.005$ in both cases.

Fig. 8 shows a toroidal structure traced out by ${\bf u}(\theta_{4}^{\pm},t)$ and ${\bf u}(\theta_{3}^{\pm},t)$ as these trajectories approach PO1 along its stable manifold. A magnified view of this region shown in the inset illustrates the complicated shape of the manifold. The numerous small loops correspond to the fast constant amplitude oscillation along PO1. In contrast, the large spiral corresponds to a slowly decaying oscillation described by the stable Floquet vectors ${\,\bf\hat{e}}_{2}$, ${\,\bf\hat{e}}_{3}$ of PO1 with exponents $\lambda_{2,3}=-0.0044\pm 0.0264i$. The characteristic time associated with this slow oscillation is $2\pi/0.0264\approx 20\tau_{c}$ and manifests itself in the weak modulation of the state space speed during the interval $t\in[30\tau_{c},50\tau_{c}]$ in Fig. 7.

As mentioned earlier, based on dimensional arguments, DC1 should be a 2D manifold which corresponds to a one-parameter family of trajectories connecting EQ2 to EQ0. Since the corresponding interval $\theta_{n}$ is very narrow, $\theta_{2}-\theta_{1}<\Delta\theta$, we resampled a wider interval of width $2\Delta\theta$ that includes $\theta_{n}$ using 100 equally spaced initial conditions. This refinement showed that $\theta_{2}-\theta_{1}\approx 0.5\Delta\theta$ and all trajectories inside this interval converge to EQ0, i.e., DC1 is indeed a 2D manifold. Thirty trajectories from DC1 are shown (as gray curves) in Fig. 9. The projection coordinates and the viewing angle are the same as in Fig. 6. As in the case of DC2, the trajectories ${\bf u}(\theta_{1},t)$ and ${\bf u}(\theta_{2},t)$ at the left and right edges of DC1 are separatrices on the 2D unstable manifold of EQ2. Using bisection, we identified that both ${\bf u}(\theta_{1},t)$ and ${\bf u}(\theta_{2},t)$ approach PO1 closely ($D_{1}<0.005$) and well-approximate heteroclinic connections DC8, DC9 from EQ2 to PO1 shown (as black curves) in Fig. 9.

The refinement in initial conditions around $\theta_{n}$ also revealed that a few trajectories ${\bf u}(\theta,t)$ for $\theta_{1}-\Delta\theta<\theta<\theta_{1}$ approach an unstable equilibrium EQ3. The flow field corresponding to EQ3 is shown in Fig. 2(d). EQ3 has a single unstable direction in $\mathcal{S}$ and its unstable manifold is one-dimensional, similar to EQ1 (see Fig. 1). Also, we found that a pair of adjacent trajectories from EQ2 approach EQ3 and subsequently depart in opposite directions. Hence, using bisection we computed the heteroclinic connection DC10 from EQ2 to EQ3 (solid red curve in Fig. 9). Additionally, we also found that one of the manifold trajectories of EQ3 uneventfully converges to EQ0. This connection (DC11) is also shown (as dashed red curve) in Fig. 9. The manifold trajectory from EQ3 evolving in the opposite direction, however, also converged to EQ0 after a long turbulent excursion and hence corresponds to another connection (DC12). This qualitative difference in dynamics on the two sides of the stable manifold of EQ3 is not observed in case of EQ1. As mentioned previously, both trajectories in the 1D unstable manifold of EQ1 display long turbulent excursions.

III.3 Connections originating at PO1

As we have established, all trajectories ${\bf u}^{+}(\eta,t)$ from PO1 quickly converge to EQ0 (cf. Fig. 8). These trajectories constitute a one-parameter family of heteroclinic connections from PO1 to EQ0 which form a 2D manifold DC13. The union of DC6 and DC7 forms a 1D boundary of the union of the 2D manifolds DC2 and DC13. Similarly, the union of DC8 and DC9 forms a 1D boundary of the union of 2D manifolds DC1 and DC13.

Unlike ${\bf u}^{+}(\eta,t)$, the trajectories ${\bf u}^{-}(\eta,t)$ display turbulent excursions after departing from the neighborhood of PO1. Hence, we tested whether for some $\eta$ they approach other ECSs, as in the case of trajectories ${\bf u}(\theta,t)$ in the 2D unstable manifold of EQ2 (cf. Sect. III.1). Since ${\bf u}^{-}(\eta,t)$ is also a one-parameter family of trajectories, we followed the procedure discussed in Sect. III.1 to search for signatures of dynamical connections originating at PO1.

To begin, we computed the state space speed $s(\eta,t)$ for each trajectory ${\bf u}^{-}(\eta,t)$ (e.g., Fig. 16 in Appendix B) and calculated $s_{m}(\eta)=\min_{t}s(\eta,t)$ for $t/\tau_{c}\in[{15,35}]$. From the plot of $s_{m}$ versus $\eta$, shown in Fig. 10, we identified about 20 trajectories with $s_{m}\ll 1$ as possible dynamical connections. Surprisingly, all these trajectories converged to EQ0. This was tested by extending the corresponding ${\bf u}^{-}(\eta,t)$ until the separation from EQ0 decreased to $D_{0}<0.01$. Since a detailed analysis of dynamics along all of these trajectories is not feasible, we limit the discussion to the families of connections DC14 that correspond to the interval $(\eta_{1},\eta_{2})$ marked in Fig. 10.

About thirty trajectories (gray curves) from the family DC14 are shown in Fig. 11 which, unlike all other figures, employs a projection of the state space onto an orthonormal basis constructed from Floquet vectors ${\,\bf\hat{e}}_{1}$ (unstable, real) and ${\,\bf\hat{e}}_{2},\,{\,\bf\hat{e}}_{3}$ (stable, complex conjugate pair) of PO1. Trajectories in DC14 leave the neighborhood of PO1 along ${\bf e}_{1}$. After a brief turbulent excursion, they visit the neighborhood of EQ2 and finally converge to EQ0. An interesting feature of Fig. 11 is that several trajectories in DC14 appear to visit the neighborhood of PO1 enroute to EQ0. This raises the question of whether a homoclinic orbit ${\bf u}^{-}(\eta,t)$ of PO1 exists near the interval $(\eta_{1},\eta_{2})$. Since PO1 has one unstable direction, a homoclinic orbit should be sandwiched between trajectories that show qualitatively different dynamics Burak Budanur et al. (2018). Such a behavior is indeed observed for trajectories near $\eta_{1}$ and $\eta_{2}$ where $s_{m}$ changes abruptly. Using bisection, we refined the estimate for $\eta_{2}$ and found that ${\bf u}^{-}(\eta_{2},t)$ indeed approaches PO1 very closely (within $D_{1}<0.005$) and therefore well-approximates a homoclinic connection DC15 shown (black curve) in Fig. 11. The projection used in this figure allows visualization of the asymptotic dynamics along DC15 for both early times (on the unstable manifold of PO1) and late times (on the stable manifold of PO1). Refining the estimate for $\eta_{1}$ using bisection revealed that ${\bf u}^{-}(\eta_{1},t)$ instead converges to EQ3 (with $D_{0}=5\times 10^{-4}$), approximating a heteroclinic connection DC16 (blue curve in Fig. 11).

We found that trajectories ${\bf u}^{-}({\eta,t})$ in the interval marked $(\eta_{3},\eta_{4})$ in Fig. 10 also converge to EQ0, after a brief turbulent excursion. These trajectories constitute a one-parameter family of connections DC17 from PO1 to EQ0. Using bisection, we identified that ${\bf u}^{-}(\eta_{3},t)$ is another homoclinic orbit (DC18) of PO1 while ${\bf u}^{-}(\eta_{4},t)$ is a heteroclinic connection (DC19) from PO1 to EQ3. The shapes of DC17, DC18, and DC19 are very similar to that of DC14, DC16, and DC15, respectively, and hence are not shown. Note that we have so far inspected state space along manifold trajectories to detect signatures of dynamical connections to EQs. An alternative metric, which allows one to identify signatures of close passes to both EQs and POs is discussed in more detail in Appendix B. Trajectories ${\bf u}^{-}({\eta,t})$ originating at PO1 are analyzed using both metrics and the results are compared in Fig. 16.

III.4 Connections originating at PO2

Visual inspection of state space speed and recurrence plots (e.g., Fig. 16) revealed that some turbulent trajectories ${\bf u}^{-}(\eta,t)$ from PO1 shadow an unstable periodic orbit PO2. This orbit, which we computed using a Newton solver, has a period $T=4.61\tau_{c}$ and is only moderately repelling in $\mathcal{S}$; its leading Floquet exponents are $\lambda_{1}=0.044$, $\lambda_{2,3}=-0.019\pm 0.016i$. However, we did not find evidence of a short heteroclinic connection from PO1 to PO2, i.e., ${\bf u}^{-}(\eta,t)$ for $t/\tau_{c}\in(15,35)$ did not approach PO2 very closely. Hence, we tested whether a connection instead exists from PO2 to PO1.

The 2D unstable manifold of PO2 is composed of two sets of trajectories ${\bf u}^{\pm}(\xi,t)$ that start from initial conditions ${\bf u}^{\pm}(\xi,0)$ constructed using equation (8). Here, $\xi$ (instead of $\eta$) parametrizes states along PO2 as well as the trajectories ${\bf u}^{\pm}(\xi,t)$. Unlike the unstable manifold of PO1, we found that trajectories in both ${\bf u}^{+}(\xi,t)$ and ${\bf u}^{-}(\xi,t)$ display turbulent evolution. Hence, we analyzed ${\bf u}^{-}(\xi,t)$ as well as ${\bf u}^{+}(\xi,t)$ for signatures of dynamical connections; the plot of $s_{m}(\xi)$ for each set is included as Fig. S1 in the supplementary material. Inspecting state space speed (as well as recurrence) plots for each trajectory and following the procedure outlined in the previous sections, we identified a family of heteroclinic connections from PO2 to EQ0 (DC20) and two isolated connections to PO1 (DC21, DC22) that lie at the boundary of DC20. A connection from the family DC20 (black curve) and the isolated connection DC21 (gray curve) are shown in Fig. 12. The projection coordinates here are the same as in Fig. 8 and the viewing angle is similar. We note that DC20, DC21, and DC22 (not shown) approach their respective destinations quickly after leaving the neighborhood of PO2.

We also found that a narrow band of trajectories ${\bf u}^{+}(\xi,t)$ for $\xi\in(\xi_{1},\xi_{2})$ from PO2 approach a periodic orbit PO0 shown in Fig. 12. PO0 has a period $T=1.22\tau_{c}$ and is stable in $\mathcal{S}$, with its leading Floquet exponent being $\lambda_{1,2}=-0.0026\pm 0.059i$. Hence, a one-parameter family of connections DC23 from PO2 to PO0 was computed by evolving ${\bf u}^{+}(\xi,t)$ for $\xi\in(\xi_{1},\xi_{2})$ for $t\approx 400\tau_{c}$. The state space projection of a trajectory from the family DC23 is shown in Fig. 13. Since PO0 is stable in $\mathcal{S}$, DC23 is a 2D manifold and trajectories that correspond to the left ($\xi_{1}$) and right ($\xi_{2}$) edges separate trajectories converging to PO0 from those which become turbulent. Using bisection, we found that both ${\bf u}^{+}(\xi_{1},t)$, ${\bf u}^{+}(\xi_{2},t)$ approach a 2-torus representing an unstable quasi-periodic orbit (QPO) QP1 shown in Fig. 13 (inset). Hence, ${\bf u}^{+}(\xi_{1},t)$, ${\bf u}^{+}(\xi_{2},t)$ correspond to connections DC24, DC25 from PO2 to QP1.

The state space speed for ${\bf u}^{+}(\xi_{1},t)$ clearly shows dynamics with two different frequencies over a time interval $70\tau_{c}$ (Fig. S2 in supplementary material). Computing recurrence plots for this segment of ${\bf u}^{+}(\xi_{1},t)$, we estimated that the periods associated with large and small loops of QP1 are $T_{1}/\tau_{c}=12.21\pm 0.018$ and $T_{2}/\tau_{c}=1.21\pm 0.05$, respectively. Since the ratio $T_{1}/T_{2}\approx 10.05$ is close to an integer, we also used the Newton-Krylov solver to test whether QP1 is instead a periodic orbit with period $T\approx T_{1},\,2T_{1},$ or $3T_{1}$; in all cases the solver failed to converge, suggesting that QP1 is indeed a QPO. The procedure we used to compute QP1 also suggests that this solution (just like PO1) possesses just one unstable direction in $\mathcal{S}$ and its stable manifold separates initial conditions that quickly approach PO0 from those exhibiting transient turbulence. Computing the stability exponents of QP1, however, is beyond the scope of this study. Lastly, the trajectories from QP1 that uneventfully converge to PO0 constitute a two-parameter family of connections forming a 3D manifold (DC26).

III.5 Structural Stability of Dynamical Connections

All dynamical connections we reported so far were computed at $Re=22.05$. One may ask if these connections are robust to small changes in $Re$, i.e., whether these connections are structurally stable Smale (1967). Recall that a connection from ECS*-∞* to ECS*∞* is formed by an intersection of the unstable manifold of ECS*-∞* and the stable manifold of ECS∞, with respective dimensions $d_{u}^{-\infty}$ and $d_{s}^{\infty}$. If these two manifolds intersect (in our case we have shown that they do), their intersection will generically be of dimension $d=d_{u}^{-\infty}+d_{s}^{\infty}-N$, where $N$ is the dimension of the state space (in our case $\mathcal{S}$). We should have $d>0$ for the connection to be structurally stable.

Since both $d_{s}^{\infty}$ and $N$ are typically very large (or infinite), $d$ can be expressed in a more convenient form in terms of the co-dimension $k_{s}^{\infty}=N-d_{s}^{\infty}$ of the stable manifold of ECS∞ Halcrow et al. (2009). If $N_{u}^{\infty}$, $N_{s}^{\infty}$, and $N_{m}^{\infty}$ are the number of unstable, stable, and marginal directions of ECS∞, then $N_{u}^{\infty}+N_{s}^{\infty}+N_{m}^{\infty}=N$ and $d_{s}^{\infty}=N_{s}^{\infty}+N_{m}^{\infty}$. Hence, $k_{s}^{\infty}=N_{u}^{\infty}$ which yields $d=d_{u}^{-\infty}-N_{u}^{\infty}$. The criterion ($d>0$) for the structural stability of a connection is then simply $d_{u}^{-\infty}>N_{u}^{\infty}$. Table 2 lists the various connections we computed and the corresponding dimensions for each connection. The entire network of connections is also shown in schematic form in Fig. 14. Clearly, all of our connections satisfy the structural stability criterion. We also numerically validated the structural stability of the connections computed at $Re=22.05$ by continuing the ECSs for $Re\in[21.85,22.2]$ and analyzing their unstable manifolds. For $Re\in[21.85,22.2]$ the number of unstable directions of all the ECSs we computed (except PO2) remained unchanged. PO2 exists only for $Re\geq 22.0$ and its unstable manifold was analyzed only at $Re=22.2$.

III.6 Transient Turbulence

As we discussed previously, turbulence in the symmetry subspace $\mathcal{S}$ is transient, with most initial conditions eventually “relaminarizing” by converging to the stable equilibrium EQ0. While some states relaminarize relatively quickly, others stay turbulent for a significant interval of time. It is natural to ask what geometric structures are responsible for maintaining turbulent flow and for relaminarization. Our results show that the periodic orbit PO1 plays a key role in both processes. As we have demonstrated, PO1 possesses (at least) two distinct homoclinic connections (DC15, DC18). The presence of homoclinic connections, which lie at the intersection of the stable and unstable manifold of PO1, suggests that these manifolds intersect and form a homoclinic tangle. The tangle explains the fractal nature of the minimal state space speed shown in Fig. 10 which results from stretching and folding of the unstable manifold. More importantly, it implies the presence of a chaotic set (a chaotic repeller in our case) anchored by PO1 as well as the presence of arbitrarily long periodic orbits that visit the neighborhood of PO1 Smale (1967). These are precisely the ingredients required for transient turbulence.

To show that it is indeed the homoclinic tangle associated with PO1 that underlies transient turbulence in our system, we computed the distance $D_{1}$ ($D_{0}$) from PO1 (EQ0) to a particular long turbulent trajectory. As Fig. 15 demonstrates, the trajectory returns to the vicinity of PO1 many times before finally relaminarizing i.e., converging to EQ0 ($D_{0}\rightarrow 0$). Furthermore, just before relaminarization, this trajectory comes very close to PO1, which suggests PO1 plays an important role in this process. It should be pointed out that not all trajectories approach PO1 closely just before relaminarization. For instance, some trajectories pass through the neighborhood of EQ3 instead, as Fig. 11 illustrates. Both PO1 and EQ3 have stable manifolds with co-dimension one; states on one side of these manifolds relaminarize almost immediately and those on the other side exhibit transient turbulence. Hence, these two stable manifolds form portions of a local boundary between “laminar” and “turbulent” flows. This analogy, however, is not perfect since the chaotic set underlying the turbulent transient is not an attractor.

It should be mentioned that non-uniqueness of edge states lying on a “laminar-turbulent” boundary has been previously reported for other turbulent flows as well. For instance, Kerswell Kerswell and Tutty (2007) and Duguet *et al.*aDuguet et al. (2008b) have identified several edge states corresponding to different traveling wave solutions in a short periodic pipe. Our results provide further evidence that not only can multiple edge states coexist, they can be of different types (e.g., an EQ and a PO, in our case).

The relationship between transient turbulence and chaotic repellers has been suggested previously, mainly based on indirect evidence – the power law decay of the relaminarization times characterizing a memoryless process Hof et al. (2006); Eckhardt et al. (2007); Schneider et al. (2010); Borrero-Echeverry et al. (2010). Direct evidence, such as the presence of homoclinic tangles van Veen and Kawahara (2011); Burak Budanur et al. (2018) or a period-doubling cascade Kreilos and Eckhardt (2012), is more recent. Moreover, while the dynamics and stability in systems with heteroclinic cycles have been explored previously Kirk and Silber (1994); Krupa and Melbourne (1995), there is currently very little understanding of dynamics in the presence of both heteroclinic and homoclinic connections. Due to the relative simplicity of the numerical model and ease of experimental access, the system considered here is particularly attractive for studying the relation between transient dynamics, relaminarization, and the structure of connections.

IV Summary and Conclusions

Several recent studies on a dynamical description of fluid turbulence focused on ECSs and how they shape the state space geometry in their neighborhoods. However, complete understanding of turbulence requires a global picture which explains how the flow moves between neighborhoods of ECSs. Such a picture can be considered as a coarse description of the dynamics, in the spirit of symbolic dynamics, analogous to a route network where ECSs serve as nodes and dynamical (homo/heteroclinic) connections as links connecting the nodes. This study describes the first rigorous attempt to construct such a network for a turbulent fluid flow. We identified eight nodes and several tens of connections between them, far more than any study to date. Moreover, while most previous studies computed connections between ECSs of the same type (primarily EQs), we have identified connections between three different types of ECSs: EQ, PO, and QPO. Indeed, this is the first study to compute connections involving QPOs. We have also demonstrated, for the first time, the existence of higher (two/three) dimensional connections between ECSs, i.e., continua of trajectories from one ECS to another.

Despite the limited attention they have received, dynamical connections can play a very important role in turbulent evolution Gibson et al. (2008). For instance, dynamically dominant ECSs in the Kolmogorov-like flow are equilibria Suri et al. (2018). Being fixed points in the state space, EQs cannot guide turbulent trajectories in their neighborhoods in the same way POs or QPOs do. Therefore, connections between EQs (as well as other types of ECSs) become the dynamically dominant solutions that guide turbulent trajectories, shaping their evolution locally. Even in systems where the dominant solutions are POs or QPOs, ECSs constrain the dynamics only locally in state space and over short intervals of time. The network of connections, on the other hand, constrains the dynamics globally in state space and over arbitrarily long time intervals.

Identifying the connection network has potential applications such as forecasting Suri et al. (2017) and control of turbulent flows Hof et al. (2010); Kühnen et al. (2018). Even though quantitative prediction has a time horizon set by the leading Lyapunov exponents that characterize the sensitivity to initial conditions, qualitative predictions do not have this limitation. In principle, prediction of extreme events can be made based on the connectivity of different ECSs. Identifying the connection network can also facilitate “low-energy” control of turbulent flows, where small perturbations to the flow result in its subsequent (natural) evolution towards a particular ECS or region of state space with desired behavior Pringle et al. (2012). Connections can also provide new insight into laminar-turbulent transition in wall-bounded 3D shear flows Budanur and Hof (2017); Burak Budanur et al. (2018).

Lastly, we point out that constraining the dynamics to a symmetry invariant subspace lowered the dimensionality of unstable manifolds and dramatically simplified the procedure of computing connections between ECSs. Whether these connections are dynamically relevant in the full state space requires further exploration. Currently, we are lacking robust numerical methods for computing connections between ECSs with more than one or two unstable directions. Some approaches, such as adjoint looping have shown promise Farano et al. (2018), but whether they present a viable option for computing connections between different types of ECS remains an open problem.

Acknowledgements.

Suri thanks Burak Budanur for useful discussions and Björn Hof for providing financial support over the duration of this study. Suri acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 754411. MS and RG acknowledge funding from the National Science Foundation (CMMI-1234436, DMS-1125302, CMMI-1725587) and Defense Advanced Research Projects Agency (HR0011-16-2-0033).

Appendix A State Space Projections

In this appendix we briefly describe the procedure used to project the state space onto a low-dimensional subspace spanned by the eigenvectors (Floquet vectors) ${\,\bf\hat{e}}_{i}$ of an equilibrium (periodic orbit). Each trajectory ${\bf u}(t)$ is expressed as a linear combination of the vectors ${{\,\bf\hat{e}}_{i}}$ as follows:

[TABLE]

Here, ${\bf u}_{ecs}$ corresponds to an equilibrium (e.g., EQ0) or a point on a periodic orbit (e.g., PO1). The coefficients $a_{i}(t)$ are computed using the scalar product

[TABLE]

where ${\,\bf\hat{e}}_{i}^{\dagger}$ is the adjoint eigenvector (Floquet vector) such that ${\,\bf\hat{e}}_{i}^{\dagger}\cdot{\,\bf\hat{e}}_{j}=\delta_{ij}$ (Kronecker delta). Typically, the vectors ${\,\bf\hat{e}}_{i}$ are not orthonormal (${\,\bf\hat{e}}_{i}\cdot{\,\bf\hat{e}}_{j}\neq\delta_{ij}$). Hence, we construct orthonormalized vectors ${\,\bf\hat{e}}_{j}^{\prime}$ such that

[TABLE]

where, the matrix elements $T_{ij}={\,\bf\hat{e}}_{i}\cdot{\,\bf\hat{e}}_{j}^{\prime}$ can be computed using the orthonormality condition ${\,\bf\hat{e}}_{i}^{\prime}\cdot{\,\bf\hat{e}}_{j}^{\prime}=\delta_{ij}$. The normalized components

[TABLE]

along vectors ${\,\bf\hat{e}}_{i}^{\prime}$ are plotted to generate state space projections. Here, $D_{c}$ is the empirically estimated largest separation between two states on a long turbulent trajectory (cf. Sect. III.1).

Appendix B Recurrence analysis

In sections III.1-III.4 we identified close passes to EQs by computing state space speed $s(t)$ along turbulent trajectories (cf. Fig. 4). However, $s(t)$ is not zero (or constant) for POs, so detecting close passes to POs requires visual inspection of speed plots to find intervals of oscillatory behavior. To address this shortcoming, we tested recurrence analysis Eckmann et al. (1987) in a form similar to that discussed in Duguet *et al.*aDuguet et al. (2008a). For each trajectory ${\bf u}(t)$ we computed the normalized recurrence function

[TABLE]

where $0<\tau_{l}<\tau_{h}$ are appropriately chosen constants. Low values of $r$ indicate that the flow field at an instant $t$ nearly recurs during a later interval $[t+\tau_{l},t+\tau_{h}]$. Unlike state space speed, $r(t)=0$ for ${\bf u}(t)$ representing both EQs and POs. In the former case $\tau_{l}$ and $\tau_{h}$ are arbitrary, while in the latter case the period $T$ of the orbit should lie inside the interval $[\tau_{l},\tau_{h}]$.

Since a turbulent trajectory shadowing an ECS mimics its spatiotemporal behavior, intervals where $r(t)\ll 1$ correspond to the trajectory visiting neighborhoods of EQs or POs. To identify such visits, however, we should restrict how near (or far) into the future we search for recurrence. Choosing $\tau_{l}$ far smaller than the correlation time $\tau_{c}$ leads to spurious self-recurrence since ${\bf u}(t)$ and ${\bf u}(t+\tau)$ do not differ appreciably; the extreme case being $r(t)=0$ when $\tau_{l}=0$. The upper bound $\tau_{h}$ restricts the search to short periodic orbits (which tend to be dynamically relevant Cvitanović and Gibson (2010)) with period $T<\tau_{h}$. Besides, it also limits the overhead associated with computing $r(t)$. In our analysis we chose $\tau_{l}=\tau_{c}$ and $\tau_{h}=5\tau_{c}$.

Fig. 16(a) shows recurrence (black curve) and speed (gray curve) plots for a turbulent trajectory ${\bf u}(t)$ that originates at PO1. Initially, ${\bf u}(t)$ shadows PO1 and consequently $s(t)$ displays steady oscillations about a finite value; in contrast, $r(t)$ almost vanishes. After a brief turbulent excursion characterized by $r$ increasing to $O(1)$, ${\bf u}(t)$ visits the neighborhoods of EQ0 at $t/\tau_{c}\approx 30$ and PO2 for $33\lesssim t/\tau_{c}\lesssim 40$. Both these close passes correspond to $r$ decreasing to well below unity. Hence, to identify signatures of dynamical connections from PO1 to both EQs and POs, we computed $r(\eta,t)$ for each ${\bf u}^{-}(\eta,t)$ in the unstable manifold of PO1 (cf. Sect. III.2). We then computed $r_{m}(\eta)=\min_{t}r(\eta,t)$ for $t>15\tau_{c}$, i.e., after each trajectory initially leaves the neighborhood of PO1. The results are shown in Fig. 16(b) which compares $r_{m}(\eta)$ (black curve) with $s_{m}(\eta)$ (gray curve) for each trajectory ${\bf u}^{-}(\eta,t)$. Clearly, the prominent minima of the two metrics align, which suggests that recurrence-based analysis is capable of successfully identifying signatures of close passes to both EQs and POs. However, it is a slightly more expensive method to identify connections, compared with the minimal state space speed.