2.1 Newton descent

The main observation here is that discrete Newton iterations are an explicit Euler approximation of a continuous fictitious-time dynamical system (Saupe Reference Saupe1988). This well-known fact is rarely mentioned in the literature, prompting Smale (Reference Smale1981) to write ‘This geometry is based on an idea which is known, but not usually explicated. Namely, Newton’s method for solving $f(z)=0$ is an Euler approximation to the ordinary differential equation …’. Here, we first derive the continuous fictitious-time Newton method in the infinite-dimensional setting and then show how the discrete Newton iterations follow from there.

The norm $\Vert \boldsymbol{F}(\boldsymbol{u})\Vert _{L^{2}}$ evolves along the trajectories of (2.4) according to

where the linear operator $\mathbf{\mathscr{L}}(\boldsymbol{u};\cdot )$ is the Gâteaux derivative of $\boldsymbol{F}$ evaluated at the state $\boldsymbol{u}$ , and is defined by

Requiring $\boldsymbol{G}$ to satisfy

one obtains the continuous fictitious-time Newton descent in infinite dimensions. Substituting this expression in (2.6), we obtain

which has the exact solution $\Vert \boldsymbol{F}(\boldsymbol{u}({\it\tau}))\Vert _{L^{2}}=\Vert \boldsymbol{F}(\boldsymbol{u}(0))\Vert _{L^{2}}\exp (-{\it\tau})$ . This shows that, as long as (2.8) has a solution, $\Vert \boldsymbol{F}(\boldsymbol{u}({\it\tau}))\Vert _{L^{2}}$ decays to zero exponentially fast along the solutions of (2.4). In other words, the continuous fictitious-time Newton method converges to an equilibrium solution of (2.1) for almost all initial conditions (Saupe Reference Saupe1988; Cvitanović & Lan Reference Cvitanović, Lan and Antoniou2003; Lan & Cvitanović Reference Lan and Cvitanović2004).

For a given state $\boldsymbol{u}$ , (2.8) is a PDE that can in principle be solved for $\boldsymbol{G}$ . In practice, it is approximated through some finite truncation to take the form $\boldsymbol{{\rm\nabla}}\boldsymbol{F}(\boldsymbol{a})\boldsymbol{G}(\boldsymbol{a})=-\boldsymbol{F}(\boldsymbol{a})$ , with $\boldsymbol{a}$ being the Galerkin coefficients. The solution $\boldsymbol{G}$ of this large, but finite-dimensional, linear system is often approximated by some variant of GMRES iterations (Trefethen & Bau Reference Trefethen and Bau1997).

Furthermore, the PDE (2.4) is discretized in time to yield the explicit Euler step

with $0<{\it\delta}{\it\tau}\leqslant 1$ . The standard Newton iterations correspond to ${\it\delta}{\it\tau}=1$ , while damped Newton iterations adjust ${\it\delta}{\it\tau}$ in order to ensure that the $L^{2}$ norm decreases (Boyd & Vandenberghe Reference Boyd and Vandenberghe2004).

The above discrete iterations do not guarantee the global convergence that the continuous fictitious-time Newton descent does (Saupe Reference Saupe1988). As mentioned in the introduction, the convergence of Newton iterations is only guaranteed if the initial guess is sufficiently close to an equilibrium. In fact, if such a close initial guess is available, the convergence to the equilibrium is super-exponential (Deuflhard Reference Deuflhard2011).

Anecdotal evidence suggests that, in the context of fluid flow, the basin of attraction of Newton iterations is rather small (Waleffe Reference Waleffe1997; Viswanath Reference Viswanath2007). This domain can be enlarged by choosing a higher-order scheme for temporal discretization of (2.4). This is, however, computationally too demanding for large-dimensional systems, such as turbulent fluid flow.

2.2 Adjoint descent

The adjoint descent corresponds to an alternative choice of $\boldsymbol{G}$ that may be expressed analytically, thus avoiding the approximation involved in the Newton–GMRES iterations. Moreover, the adjoint direction can be evaluated at a relatively low computational cost, rendering a higher-order time discretization of (2.4) feasible.

To express the adjoint direction, note that (2.6) can be written as

where the adjoint $\mathbf{\mathscr{L}}^{\dagger }(\boldsymbol{u};\cdot )$ is the linear operator satisfying

Requiring

ensures the descent of the norm $\Vert \boldsymbol{F}(\boldsymbol{u}({\it\tau}))\Vert _{L^{2}}$ along the trajectories of (2.4), since

We refer to $\boldsymbol{G}$ given by (2.13) as the adjoint direction. The finite-dimensional analogue of this adjoint direction is $\boldsymbol{G}(\boldsymbol{a})=-[\boldsymbol{{\rm\nabla}}\boldsymbol{F}(\boldsymbol{a})]^{\text{T}}\boldsymbol{F}(\boldsymbol{a})$ , where $\text{T}$ denotes the usual matrix transposition.

The downside of using the adjoint descent is that, unlike the continuous fictitious-time Newton descent, (2.14) does not guarantee an exponential decay of the $L^{2}$ norm. In fact, we observe a relatively slow decay in § 4.3, and propose a hybrid adjoint–Newton algorithm to rectify this drawback.

In § 3, we derive the explicit form of the adjoint operator for the incompressible Navier–Stokes equations.

3.1 The adjoint descent equation for equilibria

Define

which is the right-hand side of (3.1a ). The reason for using the subscript $\mathbf{0}$ in $\boldsymbol{F}_{\mathbf{0}}$ becomes clear in § 3.2, where we derive the adjoint equation for travelling waves. We seek steady states $\boldsymbol{u}$ such that $\boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u})=\mathbf{0}$ .

To this end, we derive an adjoint PDE, such that, along its trajectories $\boldsymbol{u}({\it\tau})$ , the norm $\Vert \boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u}({\it\tau}))\Vert _{\mathscr{A}}$ decreases monotonically and converges to zero asymptotically. Our choice of the norm does not compromise the accuracy of the resulting equilibrium solutions of the Navier–Stokes equations: owing to the positive-definiteness of the operator $\mathscr{A}$ , $\Vert \boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u})\Vert _{\mathscr{A}}=0$ if and only if $\boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u})=\mathbf{0}$

Leaving the derivation to appendix A, the adjoint descent equations read

where $\boldsymbol{u}^{\prime \prime }=\boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u})$ , i.e.

and we use the shorthand notation $\tilde{\boldsymbol{u}}^{\prime \prime }:=\mathscr{A}\boldsymbol{u}^{\prime \prime }$ .

Two pressure-type terms $p$ and $p^{\prime \prime }$ appear in the adjoint equation. As in the case of the Navier–Stokes equation, the pressure terms play the role of the Lagrange multipliers, enforcing the divergence-free constraints on $\boldsymbol{u}$ and $\boldsymbol{u}^{\prime \prime }$ (cf. appendix A). Taking the divergence of (3.4) we get

enforcing $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime }=0$ . Similarly, taking the divergence of (3.3a ) yields

ensuring $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ .

Owing to the periodic boundary conditions, these divergence-free constraints can be easily enforced by a projection operator in the spectral space. In appendix C, we derive the Fourier-space representation of the adjoint equations (3.3). This spectral representation closely resembles that of the Navier–Stokes equations. Therefore, an existing pseudo-spectral code for the Navier–Stokes equations can also be used, with straightforward modifications, to solve the adjoint equations. Furthermore, the adjoint equations (3.3) enjoy the same degree of parallelizability as the Navier–Stokes equations.

3.2 Adjoint descent for the travelling waves

In the absence of forcing, the Navier–Stokes equation (3.1) is invariant under Galilean translations. This symmetry allows for the existence of travelling-wave (or relative equilibrium) solutions. A travelling wave has the general form $\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}(\boldsymbol{x}-\boldsymbol{c}t)$ , where $\boldsymbol{c}=(c_{1},c_{2},c_{3})\in \mathbb{R}^{3}$ is a constant wave velocity. The forcing term $\boldsymbol{f}$ may break the translational symmetry in one or more directions, in which case, the wave speed $c_{i}$ corresponding to the symmetry-broken direction $x_{i}$ is identically zero.

Substituting the travelling-wave solution $\boldsymbol{u}(\boldsymbol{x}-\boldsymbol{c}t)$ into the Navier–Stokes equation (3.1) gives the time-independent equations

Therefore, finding travelling-wave solutions to the Navier–Stokes equation amounts to finding the kernel of the operator

If the wave velocity $\boldsymbol{c}$ is known, the solutions to $\boldsymbol{F}_{\boldsymbol{c}}(\boldsymbol{u})=\mathbf{0}$ can be found through an adjoint descent equation similar to (3.3). In general, however, the wave velocity $\boldsymbol{c}$ is unknown. To address this general case, we allow $\boldsymbol{c}$ to be a function of the fictitious time ${\it\tau}$ , and enforce its derivative $\dot{\boldsymbol{c}}=\text{d}\boldsymbol{c}/\text{d}{\it\tau}$ to change in such a way that $\Vert \boldsymbol{F}_{\boldsymbol{c}({\it\tau})}(\boldsymbol{u}({\it\tau}))\Vert _{\mathscr{A}}$ decreases monotonically to zero as ${\it\tau}$ increases.

Leaving the details to appendix B, we obtain the set of adjoint descent equations

where $\boldsymbol{u}^{\prime \prime }=\boldsymbol{F}_{\boldsymbol{c}}(\boldsymbol{u})$ , i.e.

and $\tilde{\boldsymbol{u}}^{\prime \prime }=\mathscr{A}\boldsymbol{u}^{\prime \prime }$ . The pressure-type functions $p$ and $p^{\prime \prime }$ satisfy (3.5) and (3.6), respectively, with $\boldsymbol{u}^{\prime \prime }$ given by (3.10). The spectral representation of the adjoint equations (3.9) is given in appendix C.

The adjoint descent for travelling waves includes the evolution equation (3.9b ) for the unknown wave velocity $\boldsymbol{c}$ . If the symmetry in the $x_{i}$ direction is broken due to the forcing term $\boldsymbol{f}$ , the corresponding wave speed $c_{i}$ is identically zero and the equation for ${\dot{c}}_{i}$ is eliminated from (3.9b ).

3.3 Choice of the norm

The conventional choice for the norm $\Vert \cdot \Vert _{\mathscr{A}}$ is the $L^{2}$ norm, corresponding to $\mathscr{A}$ being the identity operation. We find, however, that the resulting equations are stiff, requiring very small time steps for their numerical integration. This is in line with a similar observation made by Yang & Lakoba (Reference Yang and Lakoba2007) in the context of nonlinear wave equations. Here we choose the $H^{-1}$ norm that renders the adjoint PDEs less stiff by damping the higher Fourier modes.

This choice is somewhat arbitrary, but we find it to work very well for the Kolmogorov flow presented in § 4.1. Moreover, the implementation of the resulting equations in the Fourier space is straightforward.

The $H^{-1}$ can be computed as follows. Let $\boldsymbol{q}=(\boldsymbol{u},p)$ be the velocity–pressure pair and define the operator $\mathscr{A}$ through its action in the Fourier space,

where $\widehat{\boldsymbol{q}}$ denotes the Fourier transform of $\boldsymbol{q}$ and $\boldsymbol{k}=(2{\rm\pi}k_{1}/L_{1},2{\rm\pi}k_{2}/L_{2},2{\rm\pi}k_{3}/L_{3})$ , with $k_{i}\in \mathbb{Z}$ being the wavenumber. Note that the operator $\mathscr{A}$ commutes with the divergence operation, i.e.  $\mathscr{A}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\mathscr{A}\boldsymbol{u}=0$ .

With this choice and using Parseval’s identity, the inner product (2.2) can be written explicitly in the Fourier space as

where $^{\ast }$ denotes complex conjugation. This inner product induces the $H^{-1}$ norm $\Vert \boldsymbol{q}\Vert _{H^{-1}}^{2}=\langle \boldsymbol{q},\boldsymbol{q}\rangle _{H^{-1}}$ .

The operator $\mathscr{A}$ appears in the adjoint equations (3.3) and (3.9) through $\tilde{\boldsymbol{u}}^{\prime \prime }=\mathscr{A}\boldsymbol{u}^{\prime \prime }$ , which can be readily computed in the Fourier space through (3.11).

4.1 Kolmogorov flow

Considered the forced Navier–Stokes equation (3.1) on the two-dimensional torus $\boldsymbol{x}=(x_{1},x_{2})\in \mathbb{T}^{2}=[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ with the forcing $\boldsymbol{f}(\boldsymbol{x})=\sin (nx_{2})\boldsymbol{e}_{1}$ , where $\boldsymbol{e}_{1}=(1,0)^{\text{T}}$ and $n$ is a positive integer. This choice of forcing yields the two-dimensional Kolmogorov equation,

where $\boldsymbol{u}=(u_{1},u_{2})$ and ${\it\nu}=1/Re$ with $Re$ being the Reynolds number.

Here, we review some of the relevant properties of this equation that can also be found in Platt, Sirovich & Fitzmaurice (Reference Platt, Sirovich and Fitzmaurice1991) and Chandler & Kerswell (Reference Chandler and Kerswell2013). The Kolmogorov equation has an equilibrium solution that can be expressed analytically as

at any Reynolds number $Re$ . Following Chandler & Kerswell (Reference Chandler and Kerswell2013), we refer to this solution as the laminar state $E_{0}$ .

The instantaneous energy $E$ , energy dissipation $D$ and energy input $I$ of a state $\boldsymbol{u}$ are defined, respectively, by

where $|\cdot |$ denotes the usual Euclidean norm and ${\it\omega}$ is the non-zero component of the vorticity field $\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}$ . One can readily show that

along an arbitrary solution of the Navier–Stokes equation. For equilibrium and travelling-wave solutions, the energy is time-independent and therefore $I=D=\text{const}$ . For the laminar solution $E_{0}$ , for instance, we have

For $n=1$ , the laminar state is the global attractor for all Reynolds numbers, and therefore all solutions converge to it asymptotically (see e.g. Marchioro Reference Marchioro1986; Foias et al. Reference Foias, Manley, Rosa and Temam2001, § 3.1). For $n>1$ and large enough Reynolds numbers, the laminar state becomes unstable. The numerical study of Platt et al. (Reference Platt, Sirovich and Fitzmaurice1991) suggests that, for $n=4$ and large enough Reynolds numbers, all invariant solutions of the Kolmogorov equation become unstable, resulting in a chaotic attractor. More recently, Chandler & Kerswell (Reference Chandler and Kerswell2013) confirmed this observation. To be in agreement with these studies, we also choose the forcing wavenumber $n=4$ . Our numerical results are carried out at $Re=40$ , unless otherwise mentioned explicitly. Chandler & Kerswell (Reference Chandler and Kerswell2013) also carry out their most exhaustive search for invariant solutions at $Re=40$ .

4.2 Temporal and spatial discretization

To evaluate the right-hand side of the Navier–Stokes equation (4.1) and the adjoint equations (3.3) and (3.9), we use a standard pseudo-spectral method with $2/3$ dealiasing. Chandler & Kerswell (Reference Chandler and Kerswell2013) consider $Re=40$ , $60$ , $80$ and $100$ using uniform $256\times 256$ Fourier modes for all Reynolds numbers. For $Re=40$ , this resolution is unnecessarily high. Instead, we use $128\times 128$ Fourier modes for $Re=40$ . The resulting energy spectra (cf. figure 1), and the fact that we have been able to reproduce the (relative) equilibria found by Chandler & Kerswell (Reference Chandler and Kerswell2013), reassures that $128^{2}$ modes are sufficient.

Figure 1.

Typical energy spectrum of a turbulent trajectory at $Re=40$ .

For the time integration of the adjoint equations (3.3) and (3.9), one can take an explicit adaptive Euler time step. Denoting the right-hand side of either adjoint equation by $\boldsymbol{G}$ , the $i$ th Euler step reads

where the length ${\it\delta}{\it\tau}_{i}$ of the time step is chosen small enough such that the residue $\Vert \boldsymbol{F}(\boldsymbol{u}_{i+1})\Vert _{\mathscr{A}}$ is less than $\Vert \boldsymbol{F}(\boldsymbol{u}_{i})\Vert _{\mathscr{A}}$ .

Yang & Lakoba (Reference Yang and Lakoba2007) derive an upper bound for the admissible time step ${\it\delta}{\it\tau}_{i}$ . Computing this upper bound for the Navier–Stokes equation is, however, not straightforward. In practice, therefore, one starts with a large value ${\it\delta}{\it\tau}_{i}$ and decreases it incrementally until the residue $\Vert \boldsymbol{F}(\boldsymbol{u}_{i+1})\Vert _{\mathscr{A}}$ reduces compared to the previous iteration.

To gain accuracy, however, we use a higher-order numerical scheme for our time integrations. This is computationally feasible due to the low cost of evaluating the right-hand sides of (3.3) and (3.9) by the pseudo-spectral method. More specifically, we use the embedded Runge–Kutta scheme RK5(4) of Dormand & Prince (Reference Dormand and Prince1980). This scheme allows for an automatic adaptive time step size. Roughly speaking, RK5(4) takes fourth- and fifth-order Runge–Kutta (RK) steps. The fifth order is eventually used for the time stepping. The fourth-order prediction is used to adjust the step size as follows. Let err denote the absolute difference between the fourth- and fifth-order predictions. Then the step size is chosen such that $\mathtt{err}<\mathtt{atol}+|\boldsymbol{u}_{i}|\mathtt{rtol}$ . The prescribed constants atol and rtol are the absolute and relative errors, respectively. We refer the reader to Press et al. (Reference Press, Teukolsky, Vetterling and Flannery2007, § 17.2) for further details and an implementation of the RK5(4) scheme. This integrator is also implemented in the MATLAB function ode45. For integrating the adjoint equations, we choose $\mathtt{atol}=\mathtt{rtol}=10^{-10}$ . Even with this conservative choice, step sizes as large as $10$ fictitious-time units were taken by RK5(4).

We also use RK5(4) for the temporal discretization of the Navier–Stokes equation (4.1). However, as the time step sizes were much smaller for integrating this equation, we used the less stringent error tolerances $\mathtt{atol}=\mathtt{rtol}=10^{-5}$ .

4.3 Hybrid adjoint–Newton iterations

We find that the adjoint equation (3.3) does take arbitrary initial conditions to Navier–Stokes equilibria. The convergence is, however, rather slow. To demonstrate this, we take the initial condition $\boldsymbol{u}(\boldsymbol{x},0)=(\cos (2x_{2}),\cos (x_{1}))$ (see figure 2 a) and evolve it under the adjoint equation (3.3) to time ${\it\tau}=500$ . The result is shown in figure 2(b). This integration took 54 s (on an iMac with a single processor: Intel Core i5, 2.9 GHz). The $L^{2}$ norm of the residue $\boldsymbol{F}_{\mathbf{0}}$ , however, decreases to $\simeq 5\times 10^{-2}$ approximately, indicating that much longer integration times are required to reduce the residue sufficiently enough, say to $10^{-10}$ .

Figure 2.

(a) The vorticity field for the initial condition $\boldsymbol{u}(\boldsymbol{x},0)=(\cos (2x_{2}),\cos (x_{1}))^{\text{T}}$ . (b) The initial condition evolved to time ${\it\tau}=500$ under the adjoint equation (3.3) to obtain $\boldsymbol{u}(\boldsymbol{x},500)$ . (c) The result $\boldsymbol{u}(\boldsymbol{x},500)$ shown in panel (b) is used as the initial guess for the Newton–GMRES–hook iterations, which converged after seven iterations with residue $5.98\times 10^{-11}$ . This equilibrium solution is labelled in table 1 as $E_{4}$ . (d) The evolution of the $L^{2}$ residue $\Vert \boldsymbol{F}_{\mathbf{0}}(\boldsymbol{u}(\cdot ,{\it\tau}))\Vert _{L^{2}}$ as the initial condition in panel (a) evolves under the adjoint equation from ${\it\tau}=0$ to ${\it\tau}=500$ .

This slow convergence is due to nearly marginal, a priori unknown eigenmodes of the adjoint operator (Lakoba & Yang Reference Lakoba and Yang2007). For particular wave equations, ad hoc methods have been proposed to eliminate these modes, and hence speed up the convergence (Yang & Lakoba Reference Yang and Lakoba2007). Owing to the complexity of the Navier–Stokes equations, it is unclear how such mode elimination techniques could be employed here.

Using the state $\boldsymbol{u}(\boldsymbol{x},0)=(\cos (2x_{2}),\cos (x_{1}))$ directly as the initial guess for the NGh iterations does not converge to an equilibrium either: after the first $20$ iterations, the residue plateaued around $2\times 10^{-2}$ and remained so for the $75$ iterations that were carried out.

Instead, when we use the result of the adjoint integration, i.e.  $\boldsymbol{u}(\boldsymbol{x},{\it\tau}=500)$ , as the initial guess for NGh iterations, it converges after seven iterations with residue $5.98\times 10^{-11}$ (figure 2 c). These seven iterations took 156.84 s.

This suggests that, although the convergence of the adjoint equation to the equilibrium is slow, it takes generic initial guesses to the vicinity of the equilibrium at a relatively low computational cost. This can be seen by comparing panels (b) and (c) of figure 2, which shows that the state $\boldsymbol{u}(\boldsymbol{x},{\it\tau}=500)$ is quite similar to the equilibrium $E_{4}$ , and thus likely to lie in its NGh domain of attraction. Indeed, switching to NGh took this state to the equilibrium solution within a few iterations.

Figure 3.

State space depiction of the hybrid adjoint–Newton method. An initial guess may be too far away from the desired equilibrium solution EQ for the Newton–GMRES–hook iterations (NGh) to converge to it. The adjoint descent may eventually converge to EQ, but alone it would take a long integration time to do so. In the hybrid adjoint–Newton method the adjoint descent takes the initial guess to the domain of attraction of NGh. Once there, NGh converges to EQ in a few iterations.

Based on this observation, we propose the following hybrid adjoint–Newton iterations, sketched in figure 3. We take a prescribed positive real number ${\it\tau}_{0}$ and a positive integer N. An initial condition $\boldsymbol{u}(\boldsymbol{x},0)$ is then integrated under the adjoint equation (3.3) to obtain $\boldsymbol{u}(\boldsymbol{x},{\it\tau}_{0})$ . This state is then fed into the NGh algorithm and N iterations of NGh are performed. The output is then fed back into the adjoint equation as the initial condition. This hybrid loop is repeated until the residue $\Vert \boldsymbol{F}_{\mathbf{0}}\Vert _{L^{2}}$ decreases below some prescribed tolerance tol. This procedure is summarized in Algorithm 1. The hybrid adjoint–Newton iterations for finding travelling waves are similar except that the adjoint equations (3.9) are solved at each iteration (see Algorithm 2).

For the computations reported below, we set ${\it\tau}_{0}=100$ , $\mathtt{tol}=10^{-10}$ and $\mathtt{N}=1$ . Since NGh steps are relatively expensive, we only take one NGh step ( $\mathtt{N}=1$ ) per iteration of adjoint–Newton. At Reynolds number $Re=40$ , for instance, the integration of the adjoint equations to ${\it\tau}=100$ took approximately 10 s on average while each NGh step took approximately 50 s on average.

4.4 Equilibrium solutions

In this section, we report the equilibria found by the hybrid adjoint–Newton iterations (Algorithm 1). We also compare its performance with that of pure NGh iterations without the adjoint step.

To this end, we consider the set of initial guesses

for a range of integers $m_{1}$ and $m_{2}$ . We find that for $m_{2}>4$ , both Newton–GMRES–hook and the hybrid adjoint–Newton iterations converge to the laminar equilibrium $E_{0}$ . Therefore, we restrict the range of the integers to $1\leqslant m_{1},m_{2}\leqslant 4$ . This leads to $16$ distinct initial guesses.

Figure 4.

The decay of the residue (i.e. the $L^{2}$ -norm of the right-hand side) using Newton–GMRES–hook iterations (a) and the hybrid adjoint–Newton iterations (b).

Figure 5.

The convergence diagram. Double indices $(m_{1},m_{2})$ label the initial guesses (4.7). A dashed black arrow indicates that both the hybrid adjoint–Newton and the NGh iterations converge to the same equilibrium. A dotted blue arrow indicates that the NGh iterations converge to a different equilibrium from the one reached by the hybrid adjoint–Newton (red arrow). A solid red arrow (with no blue arrow originating from the same initial guess) indicates that only the hybrid adjoint–Newton converged to the equilibrium, i.e. the NGh iterations alone failed to converge.

These initial conditions are divergence-free by construction, consistent with our incompressible Navier–Stokes setting. They are also explicit, rendering the following results reproducible. Furthermore, the generic nature of the initial conditions illustrates the main advantage of the adjoint descent method, namely that its convergence does not require good initial guesses. For brevity, we shall refer to the states $\boldsymbol{u}_{0}^{(m_{1},m_{2})}$ by their indices $(m_{1},m_{2})$ .

Using the hybrid adjoint–Newton iterations for equilibria (see Algorithm 1), all $16$ runs converged, resulting in $10$ distinct equilibria. These equilibria are listed as $E_{0}$ to $E_{9}$ in table 1, where $E_{0}$ is the laminar state (4.2). All initial guesses converged to an equilibrium within the first $10$ iterations of hybrid adjoint–Newton iterations, except initial guess $(2,1)$ which took $17$ iterations (see figure 4). Figure 5 shows the convergence diagram, connecting each initial guess $(m_{1},m_{2})$ to the resulting equilibrium solution.

Using NGh alone, only nine out of $16$ runs converged. They converged to seven distinct equilibria: $E_{0}$ , $E_{1}$ , $E_{6}$ , $E_{8}$ , $E_{9}$ , $E_{10}$ and $E_{11}$ (see figure 5). This comes initially as a surprise since the exhaustive search carried out by Chandler & Kerswell (Reference Chandler and Kerswell2013) only returned a single equilibrium (i.e.  $E_{1}$ ). This can be explained, however, in terms of the method used for initiating the NGh iterations. Chandler & Kerswell (Reference Chandler and Kerswell2013) use recurrences to find initial guesses for their searches, as opposed to the generic initial guesses used here. Recurrences only happen in a subset of the state space where a generic turbulent trajectory spends most of its lifetime. As a result, using recurrences for initiating the NGh searches might preferentially yield the equilibria ‘close’ to this region. We will return to this subject in § 5, where the temporal intermittency is discussed.

We also sought further equilibria using states of the form $(\sin (m_{2}x_{2}),\cos (m_{1}x_{1}))$ as the initial guesses for the hybrid adjoint–Newton iterations. This resulted in five more equilibria: $E_{12}$ to $E_{16}$ in table 1. While all of these additional searches did converge, most of them re-converged to previously discovered equilibria, including $E_{10}$ and $E_{11}$ that were only found by NGh iterations when initial guesses (4.7) were used. Figure 6 shows the vorticity field of three select equilibrium solutions.

Figure 6.

Vorticity fields of the equilibrium solutions $E_{5}$ (a), $E_{12}$  (b) and $E_{15}$  (c). All panels show the entire domain $[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ .

Table 1.

List of equilibrium solutions at $Re=40$ . Energy $E$ , energy dissipation $D$ and energy input $I$ are defined in (4.3). The leading unstable eigenvalue of the equilibrium is ${\it\mu}_{1}+\text{i}{\it\omega}_{1}$ . The dimension of the linear unstable manifold of the equilibrium is denoted by $\dim E^{u}$ .

Figure 7.

Vorticity fields of the equilibrium solutions at Reynolds number $Re=60$ (a), $Re=80$ (b) and $Re=100$ (c). All panels show the entire $[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ domain.

Some of the equilibria (e.g.  $E_{15}$ shown in figure 6) exhibit vertical bands of localized vorticity that are separated by an almost zero-vorticity background. Such spatially localized equilibria of the Kolmogorov flow were only observed previously on domains with small aspect ratio ${\it\alpha}=L_{2}/L_{1}$ (Lucas & Kerswell Reference Lucas and Kerswell2014). The fact that they also exist on a domain with aspect ratio ${\it\alpha}=1$ comes as a surprise.

Although our focus here is on Reynolds number $Re=40$ , figure 7 showcases select equilibria at $Re=60$ , $80$ and $100$ . These equilibria are computed using the higher resolution of $256\times 256$ Fourier modes. They are found by our hybrid adjoint–Newton method while the previous studies using Newton–GMRES–hook iterations had not been able to discover them (Chandler & Kerswell Reference Chandler and Kerswell2013).

4.5 Travelling-wave solutions

The forcing term in the Kolmogorov equation (4.1) breaks the continuous symmetry in the $x_{2}$ direction. Therefore, only travelling-wave solutions of the type $\boldsymbol{u}(\boldsymbol{x}-\boldsymbol{c}t)$ with $\boldsymbol{c}=(c,0)$ are permitted. This reduces the wave velocity equations (3.9b ) to the scalar equation

Similarly, the term $\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\tilde{\boldsymbol{u}}^{\prime \prime }$ in (3.9a ) reduces to $c\,\partial _{x_{1}}\tilde{\boldsymbol{u}}^{\prime \prime }$ and the term $\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ in (3.10) reduces to $c\,\partial _{x_{1}}\boldsymbol{u}$ .

We search for travelling waves using Algorithm 2 and the generic initial guesses discussed in the previous section. For the initial wave speed, we used $c(0)=1$ . Other values of $c(0)$ yielded similar results.

Figure 8.

Vorticity field for the travelling-wave solutions $T_{4}$  (a),  $T_{7}$  (b) and $T_{8}$  (c). All panels show the entire domain $[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ .

Table 2.

List of travelling-wave solutions at $Re=40$ . The constant $c$ denotes the wave speed. Energy $E$ , energy dissipation $D$ and energy input $I$ are defined in (4.3). The leading unstable eigenvalue of the travelling wave is ${\it\mu}_{1}+\text{i}{\it\omega}_{1}$ . The dimension of the linear unstable manifold of the travelling wave is denoted by $\dim E^{u}$ .

Some of our searches for travelling waves converged to equilibrium solutions. This is to be expected, as equilibria are degenerate travelling waves with wave speed $c=0$ . In fact, the adjoint equation (3.9) admits such solutions. This can be readily verified by letting $\boldsymbol{u}(\boldsymbol{x})$ be an equilibrium solution of the Navier–Stokes equation and setting $\boldsymbol{c}\equiv 0$ . Then $\boldsymbol{u}^{\prime \prime }$ in (3.10) is zero, resulting in vanishing right-hand sides in (3.9).

Nonetheless, our hybrid adjoint–Newton searches led to nine distinct travelling waves listed in table 2. Only travelling waves $T_{1}$ and $T_{2}$ had been discovered previously (Chandler & Kerswell Reference Chandler and Kerswell2013). Figure 8 shows the vorticity field for three select travelling waves. As in the case of equilibria, we find that some of the travelling waves (e.g.  $T_{8}$ in figure 8) exhibit localized spatial structures, although the domain aspect ratio is one.

Figure 9.

Results for $Re=40$ . Grey curve: turbulent trajectory spanning $10^{3}$ time units. Red circles: equilibria. Blue squares: travelling waves. The green square marks the region where $I/I_{lam}<0.12$ and $D/D_{lam}<0.12$ . The turbulent trajectory spends $86.62\,\%$ of the total $10^{3}$ time units inside this region. The diagonal $I=D$ is marked by the dashed black line. Equilibria and travelling waves with $I/I_{lam}=D/D_{lam}>0.32$ are not shown.

5.1 A dynamical systems perspective on temporal intermittency

Although the Navier–Stokes equations generate an infinite-dimensional dynamical system, it is believed that, owing to the dissipative term ${\it\nu}{\rm\Delta}\boldsymbol{u}$ , its solutions converge exponentially fast to a finite-dimensional, invariant set, usually referred to as the inertial manifold (see Constantin et al. (Reference Constantin, Foias, Nicolaenko and Temam1989) for the rigorous definition). The existence of the inertial manifold for the Navier–Stokes equation, in its most general form, is an open mathematical problem. In practice, however, its existence is often assumed. In fact, this assumption underlies the finite Galerkin truncations used in computations (Foias et al. Reference Foias, Jolly, Kevrekidis, Sell and Titi1988).

Some (relative) equilibria and (relative) periodic orbits and portions of their stable and unstable manifolds belong to the inertial manifold. For large enough Reynolds numbers, these invariant solutions are typically unstable. A generic turbulent trajectory visits the neighbourhood of an invariant solution for a finite time before it is repelled along its unstable manifold towards the neighbourhood of another invariant solution (Ruelle Reference Ruelle and Sirovich1991; Halcrow et al. Reference Halcrow, Gibson, Cvitanović and Viswanath2009). This process continues indefinitely in a somewhat unpredictable fashion, thereby causing the complex temporal behaviour of turbulent trajectories (see figure 11 for an illustration).

Figure 11.

An initial condition (black square) decays rapidly to the inertial manifold $\mathscr{M}$ where the dynamics is governed by the (relative) equilibria, (relative) periodic orbits and portions of their stable/unstable manifolds that lie within the inertial manifold. The green dots represent equilibria. Highly unstable invariant solutions (e.g. the equilibrium $E$ ) are rarely visited by a generic trajectory.

From this perspective, intermittent episodes are viewed as close passages of the turbulent trajectory to invariant solutions that reside in a ‘less accessible’ region of the inertial manifold or, more precisely, the attractor (Holmes Reference Holmes1993). As depicted in figure 11, such passages are viable along the heteroclinic connections between the invariant solutions (Holmes & Stone Reference Holmes, Stone, Gatski, Speziale and Sarkar1992; Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012).

To characterize the less accessible regions of the attractor, one naturally needs to answer the following question: how frequently is an invariant solution visited by generic turbulent trajectories? There is no straightforward, a priori answer to this question (Cvitanović et al. Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2014, ch. 23). There are, however, some characteristics of the invariant solutions that are relevant. In tables 1 and 2, for instance, we report the dimension of the linear unstable manifold (i.e.  $\dim E^{u}$ ) of each invariant solution. The invariant solutions $E_{1}$ , $E_{4}$ , $E_{7}$ , $T_{1}$ , $T_{2}$ , $T_{3}$ , $T_{7}$ , $T_{8}$ and $T_{9}$ that reside close to the ergodic sea have at most $10$ linearly unstable eigenmodes. The remaining invariant solutions have at least $13$ unstable eigenmodes and seem to reside further away from the ergodic sea.

In tables 1 and 2, we also report the stability exponent ${\it\mu}_{1}+\text{i}{\it\omega}_{1}$ of the most unstable eigenmode of each invariant solution. It is tempting to assert that the invariant solutions with larger ${\it\mu}_{1}$ are less likely to be visited by a generic turbulent trajectory. This is, however, not the case. For instance, we have ${\it\mu}_{1}=0.62697$ for equilibrium $E_{4}$ and ${\it\mu}_{1}=0.61189$ for equilibrium $E_{8}$ . Equilibrium $E_{4}$ is located in the heart of the ergodic sea and is, in fact, visited by the turbulent trajectory quite often. Equilibrium $E_{8}$ , in spite of having a similar stability exponent, is rarely visited by the turbulent trajectory. The crucial difference between these two equilibria is the dimension of the linear unstable manifold, which is $\dim E^{u}=5$ for $E_{4}$ and $\dim E^{u}=17$ for $E_{8}$ .

In retrospect, the lack of correlation between the stability exponent and the frequency of visitations by the turbulent trajectory is to be expected. The stability exponent of an invariant solution is a local quantity. As such, for it to be meaningful, the turbulent trajectory should already be in the vicinity of the invariant solution. Once there, the stability exponent ${\it\mu}_{1}$ determines how quickly the trajectory will leave the neighbourhood.

At any rate, neither the dimension of the unstable manifold nor the stability exponents of an invariant solution decisively determine the frequency at which its neighbourhood is visited by a turbulent trajectory. Therefore, we take a more direct approach to quantify intermittency. Namely, we measure the $L^{2}$ distance between the turbulent trajectory and the computed invariant solutions.

The Kolmogorov equation is equivariant under a one-parameter family of continuous symmetries and a number of discrete symmetries (Sirovich Reference Sirovich1987). This implies that each solution $\boldsymbol{u}(t)$ has infinitely many equivalent copies. Therefore, when measuring the ‘distance’ between two states, one needs to make an informed choice among the equivalent copies of each state. This necessitates a discussion on the symmetries of the Kolmogorov flow.

5.2 Symmetries of Kolmogorov flow

The Kolmogorov equation (4.1) is equivariant with respect to $4n$ discrete symmetries and a continuous translational symmetry (Sirovich Reference Sirovich1987). We denote the complete set of such symmetries by

where $m\in \{0,1,\ldots ,2n-1\}$ . Here, $\mathscr{T}_{\ell }$ denotes $\ell$ shift in the $x_{1}$ direction, $\mathscr{R}$ denotes rotation through ${\rm\pi}$ and $\mathscr{S}$ denotes a simultaneous $({\rm\pi}/n)$ shift in the $x_{2}$ direction and a reflection in the $x_{1}$ direction.

One can readily verify that the above symmetry operations act on the vorticity field according to the following rules:

The glide reflection (or shift–reflect operation) $\mathscr{S}$ generates a cyclic group of order $2n$ ,

where $e$ denotes the identity $e=\mathscr{S}^{0}$ . The rotation-through- ${\rm\pi}$ operation $\mathscr{R}$ generates a cyclic group of order two, $R_{2}=\{e,\mathscr{R}\}$ . The complete set of discrete symmetries of the Kolmogorov equation, therefore, is the dihedral group of order $4n$ , i.e.

Note that the operations $\mathscr{R}$ and $\mathscr{S}$ do not commute, $\mathscr{R}\mathscr{S}\neq \mathscr{S}\mathscr{R}$ . Instead, we have $\mathscr{S}\mathscr{R}\mathscr{S}=\mathscr{R}$ .

Therefore, the solutions of the Kolmogorov equation have up to $4n$ distinct but equivalent copies due to its equivariance under the discrete symmetries $D_{4n}$ . They also have infinitely many equivalent copies due to the continuous symmetry $\mathscr{T}_{\ell }$ for any $\ell \in [0,2{\rm\pi}]$ .

An invariant solution may itself have some, all or none of the symmetries of the equations. The laminar solution $E_{0}$ , for instance, has the complete set of symmetries, i.e.  $gE_{0}=E_{0}$ for all $g\in D_{4n}$ and $\mathscr{T}_{\ell }E_{0}=E_{0}$ for all $\ell \in [0,2{\rm\pi}]$ . The laminar state, therefore, has only one copy. The travelling wave $T_{7}$ (see figure 8), on the other hand, has no symmetries and therefore possesses infinitely many equivalent copies. Incidentally, travelling wave $T_{7}$ happens to have the lowest-dimensional unstable manifold, $\dim E^{u}=2$ , among the solutions found here.

These symmetry-related copies greatly complicate the analysis of the state space of the Kolmogorov flow. When comparing the $L^{2}$ distance between two states $\boldsymbol{u}^{1}$ and $\boldsymbol{u}^{2}$ , one needs to take the minimum $L^{2}$ distance between $\boldsymbol{u}^{1}$ and $\boldsymbol{u}^{2}$ and all their symmetry-related copies. For example, let $\boldsymbol{u}^{2}(t)$ be a symmetry copy of a solution $\boldsymbol{u}^{1}(t)$ such that $\boldsymbol{u}^{2}(t)=(\mathscr{T}_{\ell }\boldsymbol{u}^{1})(t)$ for some $\ell \in (0,2{\rm\pi})$ . These two states are equivalent and both solve the Kolmogorov equation. However, the norm $\Vert \mathscr{T}_{\ell }\boldsymbol{u}^{1}-\boldsymbol{u}^{1}\Vert _{L^{2}}$ is generally non-zero. Therefore, an appropriate norm on the state space of the Kolmogorov flow is

where the minimum is taken over all $\ell \in [0,2{\rm\pi}]$ and $g\in D_{4n}$ . Evaluating this norm can be somewhat cumbersome.

Another approach is to map each state into a symmetry-invariant polynomial basis. For low-dimensional dynamical systems with simple discrete symmetries, such coordinates are available analytically (Gilmore & Letellier Reference Gilmore and Letellier2007). As the dimension of the system (and/or the order of the group) increases, the determination of the invariant coordinates becomes quickly prohibitive (Siminos & Cvitanović Reference Siminos and Cvitanović2011). As a result, and to the best of our knowledge, symmetry-invariant polynomial coordinates for the Kolmogorov flow are not known.

Here, we take an alternative approach, which also proves to be insightful in analysis of the temporal intermittency. We define the projection operator

which is the average over all copies of $\boldsymbol{u}$ given by the discrete symmetries $D_{4n}$ (Cvitanović et al. Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2014). The projection $\mathscr{P}\boldsymbol{u}$ is invariant under all discrete symmetries $g\in D_{4n}$ and therefore we refer to it as the symmetric part of the state $\boldsymbol{u}$ . The symmetric part of vorticity ${\it\omega}$ is defined analogously.

All $4n$ symmetry copies of the state $\boldsymbol{u}$ have a unique projection $\mathscr{P}\boldsymbol{u}$ . Working with the symmetric part of the states, therefore, eliminates the complications arising from the discrete symmetries.

Considering this symmetric part also has a physical motivation. The energy $E$ , dissipation $D$ and energy input $I$ defined in (4.3) are invariant under symmetry operations. Denoting the energy input of a state $\boldsymbol{u}$ by $I[\boldsymbol{u}]$ , we have $I[\boldsymbol{u}]=I[g\boldsymbol{u}]$ for $g$ being any symmetry of the Kolmogorov flow. The same holds for energy $E$ and dissipation $D$ . The particular linear form of the energy input $I[\boldsymbol{u}]$ implies that it is also invariant under the projection operation $\mathscr{P}$ , i.e.  $I[\boldsymbol{u}]=I[\mathscr{P}\boldsymbol{u}]$ .

Furthermore, we have

that is, the remainder $\boldsymbol{u}-\mathscr{P}\boldsymbol{u}$ does not contribute to the energy input. As discussed in § 5.1, the intermittent episodes of the flow are triggered by high energy inputs. Therefore, to study the temporal intermittency of the Kolmogorov flow, it is sufficient to consider the symmetric part $\mathscr{P}\boldsymbol{u}$ . Note, however, that, owing to the quadratic form of $E$ and $D$ , the energy and dissipation of the remainder $\boldsymbol{u}-\mathscr{P}\boldsymbol{u}$ are generally non-zero.

Figure 12.

The symmetric part of the vorticity field $\mathscr{P}{\it\omega}$ of the turbulent trajectory at $t=84$ (a), $t=94$ (b), $t=104$ (c), $t=124.2$ (d), $t=134.2$ (e) and $t=144.2$ (f). The energy input–dissipation pairs $(I,D)$ for the states are $(0.109,0.107)$ , $(0.115,0.114)$ , $(0.104,0.104)$ , $(0.095,0.094)$ , $(0.314,0.269)$ and $(0.099,0.093)$ , respectively.

The remaining continuous symmetry $\mathscr{T}_{\ell }$ is handled here by the method of slices (Cartan Reference Cartan1935; Field Reference Field1980; Rowley et al. Reference Rowley, Kevrekidis, Marsden and Lust2003). This method replaces all continuous symmetry copies of a state with a copy that belongs to a given hypersurface called the slice. This hypersurface is such that each group orbit $\mathscr{T}_{\ell }\boldsymbol{u}$ in a neighbourhood of a given state intersects the slice transversely at a unique point. Each group orbit $\mathscr{T}_{\ell }\boldsymbol{u}$ is then replaced by its unique intersection with the slice. The method of slices has only recently been used in the context of fluid dynamics (Willis et al. Reference Willis, Cvitanović and Avila2013; Willis, Short & Cvitanović Reference Willis, Short and Cvitanović2015). We use the first-Fourier-mode implementation of this method developed by Budanur et al. (Reference Budanur, Cvitanović, Davidchack and Siminos2015) (cf. appendix D for more detail).

5.3 Temporal intermittency

Figure 12 shows the symmetric part of the vorticity field $\mathscr{P}{\it\omega}$ for six select times along the turbulent trajectory. Since the forcing wavenumber is $n=4$ , there are $4n=16$ discrete symmetries. As a result, the symmetric part of each state exhibits recurring patterns that are related through discrete symmetry operations $\mathscr{R}$ and $\mathscr{S}$ . In other words, knowing the symmetric part $\mathscr{P}{\it\omega}$ on one-sixteenth of the domain $[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ , one can reproduce $\mathscr{P}{\it\omega}$ on the entire domain.

Figure 13.

The vorticity field ${\it\omega}$ for the snapshots of the turbulent trajectory shown in figure 12. The panels correspond to times $t=84$ (a), $t=94$  (b), $t=104$  (c),  $t=124.2$  (d),  $t=134.2$ (e) and $t=144.2$ (f).

Figure 12(e), showing the state at time $t=134.2$ , corresponds to the first intermittent peak in figure 10(a). A distinct change of topology occurs in the symmetric part of the vorticity field as the turbulent trajectory undergoes an intermittent episode. Before and after the episode, $\mathscr{P}{\it\omega}$ has at least two distinct co-rotating vortices in each positive (or negative) vorticity band. As the trajectory gets closer to the intermittent episode, the co-rotating vortices merge, resulting in horizontal bands of alternating positive and negative vorticity. After the episode (see figure 12 f), these bands split again into two distinct co-rotating vortices. This sharp distinction is not immediate from the full vorticity field ${\it\omega}$ (cf. figure 13). The same trend (i.e. the merger of co-rotating vortices) was observed during the intermittent episodes occurring at later times.

Figure 14.

(a) The vorticity field ${\it\omega}$ for the equilibrium solution $E_{13}$ . (b) The symmetric part of the vorticity $\mathscr{P}{\it\omega}$ for the equilibrium solution $E_{13}$ . (c) The symmetric part of the vorticity $\mathscr{P}{\it\omega}$ for the equilibrium solution $E_{6}$ .

We found by inspection that the symmetric part of equilibrium solution $E_{13}$ (see figure 14) is strikingly similar to that of the turbulent trajectory as it undergoes intermittency. This is visually appreciable from comparing figure 12(e) with figure 14(b). This observation suggests that close passages to equilibrium $E_{13}$ trigger the intermittent behaviour.

In figure 14(c), we also show the symmetric part $\mathscr{P}{\it\omega}$ for equilibrium $E_{6}$ . While the energy input–dissipation of equilibria $E_{6}$ and $E_{13}$ are very close (cf. table 1), their vorticity fields are very different. This demonstrates the fact that closeness of the energy input–dissipation of two states does not imply their closeness in the state space.

To quantitatively examine the role of equilibrium $E_{13}$ on intermittency, we consider the $L^{2}$ distance between the symmetric part of each turbulent state and the symmetric part of equilibrium $E_{13}$ , i.e.  $\Vert \mathscr{P}{\it\omega}(t)-\mathscr{P}E_{13}\Vert _{L^{2}}$ . To account for the continuous translational symmetry, each symmetric part is first brought to the first-Fourier-mode slice, as explained in appendix D, before the distance is computed.

Figure 15.

The $L^{2}$ distance $\Vert \mathscr{P}{\it\omega}(t)-\mathscr{P}E\Vert _{L^{2}}$ between the symmetric part of the turbulent trajectory $\mathscr{P}{\it\omega}$ and the symmetric parts of the equilibria $E_{13}$  (a),  $E_{11}$  (b) and $E_{1}$  (c). (d) The normalized energy input $I/I_{lam}$ of the turbulent trajectory, shown for reference (same as figure 10 a). Each distinguished peak of the energy input (i.e. intermittency) coincides with a significant decrease in the distance from the equilibrium $E_{13}$ .

Figure 15 shows the evolution of the $L^{2}$ distance as a function of time. The evolution of the energy input $I$ is also shown in the same figure. Every intermittent episode, i.e. significant peaks in the energy input, is concurrent with a significant dip in the $L^{2}$ distance. This demonstrates that the intermittent episodes correspond to close passages of the turbulent trajectory to equilibrium $E_{13}$ .

The distance between the turbulent trajectory and some of the other invariant solutions (e.g.  $E_{11}$ and $T_{4}$ ) also decreases during the intermittent episodes. Their decrease is, however, not as dramatic as the one observed for equilibrium $E_{13}$ . For instance, the evolution of the distance from $E_{11}$ is also shown in figure 15. The distance from $E_{11}$ fluctuates mostly around $14.0$ , which is significantly larger than the average distance from $E_{13}$ , which fluctuates around $8.0$ . During the intermittent episodes, the distance from both equilibria decreases. The distance from $E_{11}$ reaches $4.64$ at its minimum, which is approximately one-third of its average distance of $14.0$ . The decrease in the distance from $E_{13}$ is more significant: it attains a minimum as low as $1.21$ , which is approximately $6.6$ times smaller than its average distance of $8.0$ . In contrast, the distance from equilibria residing in the ergodic sea increases during each intermittent episode. The distance from the equilibrium $E_{1}$ is shown in figure 15 as an example.

If the Kolmogorov flow at Reynolds number $Re=40$ possesses an ergodic attractor, as the numerical evidence suggests (see e.g. Platt et al. Reference Platt, Sirovich and Fitzmaurice1991), every invariant solution embedded in the attractor will eventually be visited by a generic turbulent trajectory. Nonetheless, for the finite-time (and relatively short) turbulent trajectory computed here, equilibrium $E_{13}$ appears to explain its observed intermittencies.

Finding the invariant equilibrium and travelling-wave solutions is only the first step. A complete understanding of the Kolmogorov flow in terms of its invariant solutions will require a detailed state space analysis of the type carried out by Gibson, Halcrow & Cvitanović (Reference Gibson, Halcrow and Cvitanović2008) for the plane Couette flow. Such a thorough study deserves its own treatment, which will be presented elsewhere.