New equilibrium solution branches of plane Couette flow discovered using a project-then-search method

TL;DR

This paper introduces a novel project-then-search method to discover 16 new equilibrium solutions in plane Couette flow, leveraging resolvent mode projections to efficiently explore the flow's solution space.

Contribution

The study presents a new approach combining resolvent mode projections with Newton-Krylov searches to find previously unknown equilibrium branches in plane Couette flow.

Findings

Discovered 16 new equilibrium solutions.

Achieved a 96% convergence rate with the new method.

Identified bifurcations and disconnected solution branches.

Abstract

New equilibrium solution branches for plane Couette flow are reported which add to the inventory of known solution branches. The exact solutions are found by projecting known equilibria onto the resolvent modes of McKeon & Sharma (J. Fl. Mech., vol. 658, 2010, pp. 336-382) to generate approximate solutions that are subsequently used as seeds in a Newton-Krylov-hookstep search. Searches initialised with these projections have a convergence rate of 96%, leading to the discovery of 16 new equilibria. The low-dimensional nature of the resulting equilibria is attributed to the projection generated by the resolvent model; resolvent modes for a given equilibrium solution span a low-dimensional space which the solution approximately inhabits. This property is exploited to generate new branches of equilibria in plane Couette flow, and to jump to known branches. The new branches include bifurcations from previously known bifurcation curves and a disconnected bifurcation curve which displays interesting behaviour when continued in Reynolds number. The unstable manifolds of the new exact coherent states are also computed to expand the known state space structure of plane Couette flow.

Figures

Tables

Table 1. Table 1 : Properties of the new equilibria and the equilibria they are derived from; a dividing line separates parent equilibria (previously known solutions) from their child equilibria (new solutions) from their grandchild equilibria (solutions found from projections of the child equilibria). The ‘Root’ column denotes the initial velocity field that lead to the new discovery, e.g. a rank-5 projection of EQ9 EQ9 \mathrm{EQ}9 lead to the discovery of EQ16 EQ16 \mathrm{EQ}16 . The ‘ R ​ e 𝑅 𝑒 Re ’ column indicates the distinct Reynolds number that the solution was discovered at. The ‘mean’ values are given for comparison, and are calculated from a spatially and temporally averaged turbulent flow. The ℒ 2 superscript ℒ 2 \mathscr{L}^{2} -norm of the velocity field is ‖ 𝒖 ~ ‖ norm ~ 𝒖 \|\tilde{\boldsymbol{u}}\| , E 𝐸 E is the kinetic energy density, D 𝐷 D is the dissipation rate, H 𝐻 H is the isotropy subgroup, dim( W u superscript 𝑊 𝑢 W^{u} ) is the dimensionality of the equilibrium’s unstable manifold, and dim( W H u subscript superscript 𝑊 𝑢 𝐻 W^{u}_{H} ) is the dimensionality of the unstable manifold within the H 𝐻 H -invariant subspace. The accuracy of the solution (Acc.) is calculated using ( 3 ). The return column (Ret.) denotes how often the solution was found using projections of itself. The occurrence column (Occ.) denotes how often the solution was found from searches initiated using projections of other equilibria.

Root	$Re$	$\mathrm{EQ}$	$\Vert \widetilde{𝒖}\Vert$	$E$	$D$	$H$	dim($W^u$)	dim($W_H^u$)	Acc.	Ret.	Occ.
	400	mean	0.2997	0.1016	2.6017	$\left\{e\right\}$
		0	0.0000	0.1667	1.0000	$\Gamma$	0	0		0	199
${\Pi }_3\left(\mathrm{EQ4}\right)$	400	3	0.1259	0.1382	1.3177	$\Sigma$	4	2	${10}^{-4}$	73	36
	400	4	0.1681	0.1243	1.4537	$\Sigma$	6	3	${10}^{-6}$	73	2
${\Pi }_1\left(\mathrm{EQ2}\right)$	400	1	0.2091	0.1363	1.4293	$\Sigma$	1	1	${10}^{-6}$	97	66
	400	5	0.2186	0.1073	2.0201	$\Sigma$	11	4	${10}^{-3}$	32	1
${\Pi }_5\left(\mathrm{EQ23}\right)$	400	2	0.3858	0.0780	3.0437	$\Sigma$	8	2	${10}^{-4}$	96	62
${\Pi }_2\left(\mathrm{EQ5}\right)$	400	7	0.0936	0.1469	1.2523	$\Theta$	3	1	${10}^{-4}$	66	44
${\Pi }_7\left(\mathrm{EQ4}\right)$	400	9	0.1565	0.1290	1.4048	${\Theta }_3$	5	3	${10}^{-4}$	74	6
	400	10	0.3285	0.1080	2.3721	${\Theta }_3$	10	7	${10}^{-4}$	48	1
	400	11	0.4049	0.0803	3.4322	${\Theta }_3$	13	10	${10}^{-3}$	1	1
${\Pi }_5\left(\mathrm{EQ9}\right)$	400	20	0.2405	0.1289	1.6034	${\Theta }_3$	3	2	${10}^{-5}$	97	49
${\Pi }_3\left(\mathrm{EQ10}\right)$	400	21	0.2683	0.1242	1.7630	${\Theta }_3$	4	3	${10}^{-6}$	43	79
${\Pi }_7\left(\mathrm{EQ10}\right)$	400	22	0.3037	0.1160	2.0713	${\Theta }_3$	8	6	${10}^{-5}$	40	25
${\Pi }_7\left(\mathrm{EQ11}\right)$	400	23	0.4014	0.0759	3.2474	${\Theta }_3$	10	4	${10}^{-5}$	33	40
${\Pi }_{10}\left(\mathrm{EQ11}\right)$	400	24	0.4049	0.0813	3.3612	${\Theta }_3$	15	9	${10}^{-5}$	2	11
	330	mean	0.2541	0.0959	1.6660	$\left\{e\right\}$
	330	6	0.2751	0.0972	2.8185	$\Sigma$	19	6	${10}^{-3}$	16	6
${\Pi }_6\left(\mathrm{EQ6}\right)$	330	13	0.2168	0.1337	1.4705	$\Sigma$	1	1	${10}^{-3}$	98	26
${\Pi }_{40}\left(\mathrm{EQ6}\right)$	330	14	0.2375	0.1052	2.2785	$\Sigma$	15	6	${10}^{-2}$	57	42
${\Pi }_2\left(\mathrm{EQ6}\right)$	330	12	0.1145	0.1410	1.3433	$\Theta$	3	1	${10}^{-3}$	67	92
${\Pi }_{67}\left(\mathrm{EQ6}\right)$	330	16	0.2348	0.1063	2.3047	${\Theta }_3$	15	8	${10}^{-2}$	17	5
${\Pi }_{42}\left(\mathrm{EQ6}\right)$	330	15	0.2674	0.0988	2.6947	${\Theta }_3$	18	9	${10}^{-2}$	27	37
${\Pi }_{31}\left(\mathrm{EQ15}\right)$	330	25	0.2331	0.1292	1.5650	${\Theta }_3$	3	2	${10}^{-3}$	93	7
${\Pi }_{66}\left(\mathrm{EQ15}\right)$	330	26	0.2707	0.0975	2.6274	${\Theta }_3$	17	9	${10}^{-2}$	40	1
	270	mean	0.2286	0.1089	1.4813	$\left\{e\right\}$
	270	8	0.3466	0.0853	3.6719	$\Theta$	15	2	${10}^{-3}$	42	1
${\Pi }_{11}\left(\mathrm{EQ8}\right)$	270	18	0.2292	0.1294	1.5415	$\Sigma$	1	1	${10}^{-2}$	34	9
${\Pi }_3\left(\mathrm{EQ8}\right)$	270	17	0.1546	0.1301	1.5530	$\Theta$	5	1	${10}^{-2}$	67	32
${\Pi }_{38}\left(\mathrm{EQ8}\right)$	270	19	0.3148	0.0904	3.1529	$K$	12	0	${10}^{-2}$	67	9
${\Pi }_2\left(\mathrm{EQ18}\right)$	270	27	0.2297	0.1292	1.5444	${\Theta }_3$	2	2	${10}^{-2}$	91	64

Equations

$$⟨\mathbf{a},\mathbf{b}⟩$$

$$\Vert \mathbf{a}{\Vert }^2$$

$$G(\tilde{\mathbf{u}},\sigma ,t)=\sigma f^t(\tilde{\mathbf{u}})-\tilde{\mathbf{u}}=0,$$

$$\frac{\Vert f^{T=1}(\tilde{\mathbf{u}})-\tilde{\mathbf{u}}\Vert }{\Vert \tilde{\mathbf{u}}\Vert }.$$

$$\mathbf{u}(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{i\omega t}\hat{\mathbf{u}}(\omega )d\omega +\mathbf{T}(t)$$

$$\mathbf{u}$$

$$\frac{\partial \tilde{\mathbf{u}}}{\partial t}$$

$$=\mathcal{L}\tilde{\mathbf{u}}(t)+\tilde{\mathbf{g}}(t)$$

$$\hat{\mathbf{u}}(\omega )={\left(\omega I-\mathcal{L}\right)}^{-1}\hat{\mathbf{g}}(\omega ),$$

$${\mathcal{H}}_K(\omega )\mathbf{a}$$

$$=\sum _{m=1}^{\infty }{\mathbf{\psi }}_K^m(\omega ){\sigma }_K^m(\omega ){⟨{\mathbf{\varphi }}_K^m(\omega ),\mathbf{a}⟩}_y.$$

$$⟨{\mathbf{\varphi }}_K^m(\omega ),{\mathbf{\varphi }}_K^{m^{\prime }}(\omega )⟩$$

$$⟨{\mathbf{\psi }}_K^m(\omega ),{\mathbf{\psi }}_K^{m^{\prime }}(\omega )⟩$$

$${\sigma }_1\ge {\sigma }_2\ge \dots$$

$$\mathbf{u}(x,y,z,t)=\frac{1}{2\pi }\sum _{\alpha ,\beta }{\int }_{-\infty }^{\infty }{e}^{i(\omega t+\alpha x+\beta z)}\sum _{m=1}^{\infty }{\mathbf{\psi }}_K^m(\omega ,y)c_K^{m}d\omega .$$

$$\mathbf{u}={\Pi }_M\mathbf{u}+{\mathbf{u}}^{\perp }$$

$${\Pi }_M\mathbf{u}$$

$${\mathbf{u}}^{\perp }$$

$$\sigma \dot{\tilde{\mathbf{u}}}=\sigma f(\tilde{\mathbf{u}})=f(\sigma \tilde{\mathbf{u}}).$$

$${\sigma }_z$$

$${\sigma }_x$$

$${\sigma }_{xz}$$

$$\tau (\delta x,\delta z)$$

$${\sigma }_x=$$

$${\tau }_z{\sigma }_z=$$

$${\tau }_z{\sigma }_{xz}=$$

$${\tau }_x{\sigma }_z=$$

$${\tau }_x{\sigma }_x=$$

$${\tau }_{xz}=$$

$${\tau }_{xz}{\sigma }_x=$$

$$\Theta$$

$$K$$

$$\Sigma$$

$${\Theta }_3$$

$$\mathbb{U}=\{\mathbf{u}\in L^2(\Omega )|$$

$$\mathbf{u}(x,\pm 1,z)=0,$$

$$\mathbf{u}(x,y,z)=\mathbf{u}(x+L_x,y,z)=\mathbf{u}(x,y,z+L_z)\}.$$

$${\mathbb{U}}_H=\{\mathbf{u}\in \mathbb{U}|h\mathbf{u}=\mathbf{u},\forall h\in H\},$$

$$a_n(t)=⟨\mathbf{u}(t),{\mathbf{e}}_n⟩,$$

$$\begin{array}{ccccc}& & {\tau }_x& {\tau }_z& {\tau }_{xz}\\ {\mathbf{e}}_1={\gamma }_1(1+{\tau }_x+{\tau }_z+{\tau }_{xz}){\mathbf{u}}_{EQ2}& & S& S& S,\\ {\mathbf{e}}_2={\gamma }_2(1+{\tau }_x-{\tau }_z-{\tau }_{xz}){\mathbf{u}}_{EQ2}& & S& A& A,\\ {\mathbf{e}}_3={\gamma }_3(1-{\tau }_x+{\tau }_z-{\tau }_{xz}){\mathbf{u}}_{EQ2}& & A& S& A,\\ {\mathbf{e}}_4={\gamma }_4(1-{\tau }_x-{\tau }_z+{\tau }_{xz}){\mathbf{u}}_{EQ2}& & A& A& S,\end{array}$$

$$\mathbf{u}(0)={\mathbf{u}}_{EQ}\pm ϵ{\mathbf{v}}_{EQ},$$

Full text

1 Introduction

A complete understanding of turbulent flows is one of the most important and fascinating problems in classical physics. The problem has been looked at from various perspectives; here we adopt the nonlinear dynamical systems standpoint. Wall-bounded turbulent flows contain persistent macroscopic structure that recurs in space and time known as coherent structures. Kline et al. (1967) performed experiments which showed spatially organised streaks in the turbulent boundary layer. Further work has revealed that there is order within apparently random turbulence (Adrian, 2007; Cantwell, 1981; Holmes et al., 1998; Robinson, 1991; Sharma & McKeon, 2013). Hamilton et al. (1995) identified a well-defined quasi-cyclic process in plane Couette flow at low Reynolds number. This quasi-cyclic process eventually manifested in the ‘Self-Sustaining Process’ theory (SSP) (Waleffe, 1997) that describes the roll-streak behaviour in terms of forced response of streaks to rolls, growth of streak instabilities, and nonlinear feedback from streak instabilities to rolls. It is widely thought that similar processes are important in turbulent flow.

The dynamical systems picture can be traced back to Hopf (1948), who proposed that turbulence is described by a finite-dimensional manifold invariant under the action of the state space flow. In this picture, equilibria are fixed points that give structure to the invariant manifold, and their bifurcation structure determines behaviour across Reynolds number. In a sense, all work on turbulence is ultimately aimed at producing a concise description of this invariant manifold.

The idea is that turbulent flow can be described through invariant solutions of the NSE and their connecting orbits. There are four classes of such solution that we are concerned with. Equilibria, which are constant in time; relative equilibria (also called travelling wave solutions), which are constant in a co-moving frame of reference; periodic orbits, which repeat over a fixed time period; and relative periodic orbits, which are periodic orbits in a co-moving frame of reference.

For these reasons, Kawahara et al. (2012) pose the finding of sufficiently many invariant solutions to fully describe a turbulent flow as one of the outstanding problems in the field. Unfortunately, with current methods, computing invariant solutions involves trial and error and computational searches often fail. In what follows, we focus on the simplest objects, equilibria, and show how many new equilibria may be generated cheaply from those already known.

Trivial homogeneous equilibrium solutions of the NSE, such as the stable laminar flow state in plane Couette and Poiseuille flow, can be found with very little effort; these solutions are easy to find analytically as well as numerically. Less symmetric solutions on the other hand are more difficult to find. The most commonly used method for finding invariant solutions is the Newton search because it converges quadratically. However, as typically implemented the convergence is only local; to guarantee convergence an initial guess that is close to a solution must be provided. Finding suitable initial states for the algorithm is the tricky task. During the course of this investigation Farazmand (2016) developed a hybrid adjoint-Newton algorithm that provides global convergence from any given initial condition to find equilibria of two-dimensional Kolmogorov flow.

There are various methods that can be used to generate initial guesses for searches, the most widely used of which we list below.

1. (a)

Bisection:

This approach involves varying the initial amplitude of a velocity field depending on its temporal evolution. The time-evolved flow field which sits on the boundary between the laminar and turbulent states (a flow field that neither laminarises nor becomes fully turbulent, a so-called edge state) is used as an initial guess for a Newton search. This method is used by Duguet et al. (2008, 2009); Hwang et al. (2016) and Avila et al. (2013). A related but more sophisticated method, which modifies the Reynolds number with on-line feedback control to automatically achieve the same effect, is used by Willis et al. (2017). 2. (b)

Homotopy:

In this method a convergence control parameter is introduced and solutions are found as the parameter changes via continuation. Waleffe (2003) describes it as smoothly deforming the base flow into the desired flow while tracking the solutions with Newtonfls method. Several studies have successfully used this approach to build solution connections between Taylor-Couette, plane Couette and plane Poiseuille flow, see Nagata (1990, 1997, 1998), Faisst & Eckhardt (2003) and Waleffe (2001, 2003). 3. (c)

Recurrence plots:

This technique involves studying the time-evolution of a perturbed flow field for periodic patterns. The most dynamically important fields are then used as initial guesses in a Newton search. The method has been successfully used to find equilibria, travelling waves and periodic orbits of the NSE, see Chandler & Kerswell (2013); Cvitanović & Gibson (2010); Gibson et al. (2009). For a comprehensive review of the technique refer to the work of Marwan et al. (2007). 4. (d)

Projections:

Introduced here, low-rank projections of known equilibrium solutions are used as initial guesses in the search for other solutions. The projected velocity field resembles its parent and maintains its dominant flow characteristics (in the sense of the norm used to define the projection). The projected fields have structures derived from known solutions, making them good initial guesses for the search algorithm.

The first to discover non-trivial high-dimensional nonlinear solutions of the NSE at moderate-Reynolds number was Nagata (1990), who computed a pair of unstable three-dimensional equilibria in plane Couette flow using homotopy and bifurcations from a wavy vortex solution of Taylor-Couette flow. Continuing the two solutions downwards in Reynolds number reveals that they originate from a saddle-node bifurcation and contain a wavy low-velocity streak flanked by counter-rotating vortices; the upper branch solutions consist of weak streaks with strong vortices whereas the lower branch solutions consist of weak vortices but stronger streaks. Their structure persists at higher Reynolds numbers and resembles coherent structures observed in the near-wall region of wall-bounded turbulent flows (Jeong et al., 1997; Stretch, 1990). The same solutions were found independently by Clever & Busse (1997) and Waleffe (1998, 2003). Itano & Generalis (2009), Nagata (1997), Gibson et al. (2008, 2009) and Schmiegel (1999) found a multitude of equilibria in plane Couette flow that are not related to Nagata (1990)’s solutions but exhibit similar flow structures. Relative equilibria have been found using the homotopy approach in plane Poiseuille flow (Ehrenstein & Koch, 1991; Itano & Toh, 2001), pipe flow (Faisst & Eckhardt, 2003; Wedin & Kerswell, 2004) and square ducts (Okino et al., 2010; Uhlmann et al., 2010; Wedin et al., 2009). For a comprehensive review refer to Kawahara et al. (2012).

Gibson et al. (2008) developed a representation-independent method of visualising invariant solutions in state space, the details of which are given in §3.1.

Figure 1 depicts the state space portrait that Gibson et al. constructed. The portrait allows us to follow a transitional or turbulent trajectory as it wanders near fixed points and is directed by the unstable manifolds of the equilibria. This state space representation organises solutions relative to the laminar solution thus revealing the structure of the laminar attractor. Additionally, the portrait helps in seeing heteroclinic connections between fixed points as shown by Gibson et al. (2008), who were the first to find a heteroclinic connection between two non-trivial equilibria. Halcrow et al. (2009) found further heteroclinic connections between equilibria and demonstrated that these connections are pertinent to obtaining a global picture of turbulent flow dynamics in state space.

It is claimed that the connecting orbits between invariant solutions of the NSE can be used to describe low-dimensional processes embedded within fully developed turbulent flows (Halcrow, 2008). The high-dimensional dynamics of turbulence contains macroscopic organised structure that looks like and shares mechanisms with invariant solutions. For this reason, Waleffe (2001) referred to the exact solutions of the NSE as exact coherent structures. The solutions closely resemble coherent structures seen in experiments and direct numerical simulations (DNS), but they lack the complex spatio-temporal intermittency characteristics of coherent structures observed in turbulent flows. Therefore, understanding exact coherent structures and their connections can clarify the fundamental dynamics that form the backbone of fully turbulent flows. Jiménez & Simens (2001) showed that the underlying dynamics in the near-wall region reveal themselves in intermittent and deterministic processes that are low-dimensional and Holmes et al. (1998) argue that low-dimensional processes are embedded in global high-dimensional fully developed turbulence. From the nonlinear dynamical systems viewpoint, connecting orbits can describe these low-dimensional processes. With this in mind, we provide structure to the laminar attractor first.

The structure of the paper is as follows. The project-then-search methodology is described in §2. The classification of symmetries in plane Couette flow that support equilibria are given in §3. The new equilibria are presented and discussed in §4 and §5. Conclusions are given in §6.

2 Methodology

2.1 Plane Couette flow

The flow geometry (see figure 2) and all parameters are chosen to be consistent with the work of Gibson et al. (2009). The nondimensional NSE for an incompressible fluid are {subeqnarray} ∂\boldsymbol{}{u}$∂t +$\boldsymbol{}{u}\cdot\nabla\boldsymbol{u}&=-\nabla p+Re^{-1}\nabla^{2}\boldsymbol{u},\\ \nabla\cdot\boldsymbol{u}=0,where$$\boldsymbol{}{u}(x,y,z,t)=[u\ v\ w]^{T}$is the velocity vector in the streamwise$x$, wall-normal$y$and spanwise$z$directions,$t$is time,$p(x,y,z,t)$is the pressure,$\nabla$is the gradient operator,$\nabla^{2}$is the Laplace operator and$Re$is the Reynolds number. The Reynolds number is defined as$Re=Uh/\nu$, where$U$is half the relative velocity of the plates,$h$is the channel half-height and$\nu$ is the kinematic viscosity.

The domain has periodic boundary conditions in the streamwise and spanwise directions. Spatial periodicity is specified in terms of fundamental streamwise and spanwise Fourier wavenumbers, $\alpha$ and $\beta$, respectively. The relation between the spatial wavenumbers and the domain is given as $L_{x}=2\pi m/\alpha$ and $L_{z}=2\pi n/\beta$ where $m,n\in\mathbb{Z}^{+}$. The equilibria are computed on the domain $\Omega=[\frac{2\pi}{1.14},2,\frac{2\pi}{2.5}]$ which is discretised onto a $[32,35,32]$ grid. In this study the pressure gradient is fixed at zero and integration forward in time is performed with time step $\Delta t=0.03125$.

The unit vectors in the $x,y,z$ directions are denoted $\hat{\boldsymbol{x}},\hat{\boldsymbol{y}},\hat{\boldsymbol{z}}$ and the plane Couette base flow is defined as $y\hat{\boldsymbol{x}}$. Hence, the total velocity is defined as $\boldsymbol{u}=\tilde{\boldsymbol{u}}+y\hat{\boldsymbol{x}}$, where $\tilde{\boldsymbol{u}}$ is the velocity difference from laminar. In the present study the $\mathscr{L}^{2}$-inner product and norm are defined as

[TABLE]

The energy is defined as $\|\tilde{\boldsymbol{u}}\|=\sqrt{\langle\tilde{\boldsymbol{u}},\tilde{\boldsymbol{u}}\rangle}$, the kinetic energy density as $E=\frac{1}{2}\|\boldsymbol{u}\|^{2}$ and the dissipation rate as $D=\|\nabla\times\boldsymbol{u}\|^{2}$. The role that certain solutions play in turbulent dynamics can be interpreted from the dissipation rate, e.g. a lower dissipation rate means that the equilibrium is far from the turbulent attractor and closer to laminar flow.

2.2 Computational method for search

The new equilibria found in this report are derived from equilibria previously discovered by Nagata (1990), Gibson et al. (2008, 2009) and Halcrow et al. (2009), which can be found at channelflow.org. The open source library Channelflow was used to find, continue and analyse all solutions in our investigation (Gibson, 2014). All computational settings are the same as in the work of Gibson et al. (2009).

In the present work, the Newton-Krylov-hookstep (NKH) search algorithm developed by Viswanath (2007) is used to find invariant solutions of the NSE, for a detailed explanation of the algorithm refer to the work of Halcrow (2008) and Viswanath (2007). The algorithm finds approximate solutions to the following equation,

[TABLE]

where $f^{t}(\tilde{\boldsymbol{u}})$ is a time-mapped instance of the initial flow field at time $t$, $\tilde{\boldsymbol{u}}$ is the initial flow field and $\sigma$ is the symmetry of the flow field (see §3). The convergence criteria for the searches is $\|G\|\leq 10^{-13}$. An equilibrium solution is defined as $\tilde{\boldsymbol{u}}(x,y,z,t)=\tilde{\boldsymbol{u}}_{EQ}(x,y,z)$, therefore when there is close to no difference between the time-evolved state and the initial flow field an equilibrium solution is deemed to have been found. For the searches performed in the current work, the NKH algorithm is unconstrained with respect to symmetries, only the final solution is inspected for symmetries.

Following Gibson et al. (2009) and Viswanath (2007), to determine the accuracy of a solution it is interpolated from a grid resolution of $[32,35,32]$ onto a $[48,49,48]$ grid and then time integrated for $T=1$ with $dt=0.02$. Then the accuracy is determined by calculating the residual, defined as

[TABLE]

All searches for new equilibria and their bifurcation curves were also performed at a higher grid resolution of $[48,65,48]$ to ensure that the searches were well-resolved; all results were found to agree with those found at the lower grid resolution of $[32,35,32]$ down to the fifth significant figure. Invariant solutions have been found in the context of minimal flow units (small periodic domains just large enough to sustain turbulent flow or contain a single coherent structure, Jiménez & Moin (1991)).

2.3 Low-rank projections

The practical difficulty with using Newton’s method to find equilibria is that the search may fail or return the trivial laminar solution. The aim here is to generate initial guesses that are close enough to new equilibria to converge. As such, we generate such initial guesses by projecting known equilibria onto resolvent modes. The resolvent model is used to provide a physically relevant, ordered basis in which the velocity field can be expanded. In principle, projections onto other bases could be used.

To this end, a short summary of the pertinent aspects of the resolvent model and associated projection is provided below. For a detailed derivation of the model in the context of turbulence, the reader may refer to McKeon & Sharma (2010), and for the projection step, to Sharma et al. (2016b).

It was previously shown in Sharma et al. (2016b) that these projections are often close to the original solutions, even when the projection is very low rank. We will see that the projections are also often close to other, nearby, solutions.

Consider a long-time solution to (2.1), $\boldsymbol{u}$. This solution may be expanded into its harmonic and transient parts,

[TABLE]

where the dependence on $x,y,z$ has been suppressed. For the case where flow has already decayed onto the attractor, or for any recurrent flow, $\boldsymbol{T}=0$. Notice that the temporal mean is associated with $\hat{\boldsymbol{u}}(0)$.

Writing (2.1) as $\partial\boldsymbol{u}/\partial t=\boldsymbol{f}(\boldsymbol{u})$, an expansion about the temporal mean $\bar{\boldsymbol{u}}:=\hat{\boldsymbol{u}}(\omega=0)$ gives

[TABLE]

where $\tilde{\boldsymbol{g}}$ represents the second-order terms in the expansion of $\boldsymbol{f}$ (the Reynolds stress gradients). Note that we are performing an expansion around the mean of a specific equilibrium solution, rather than the turbulent mean, as in the earlier work. Similarly expanding $\boldsymbol{g}$ in its Fourier coefficients and rearranging gives

[TABLE]

at any $\omega\neq 0$. In this formulation, the second-order terms act to excite the state rather than being truncated.

The operator $\mathcal{H}(\omega)=\left(\omega I-\mathcal{L}\right)^{-1}$ is the resolvent of $\mathcal{L}$ and is a linear mapping from the Reynolds stress gradients to the velocity field. The idea of the projection step is to find the optimal projection $\Pi(\omega)$ of rank $M$ that approximates $\mathcal{H}(\omega)$. Since the resolvent operator is linear, the optimal projection is provided by the singular value decomposition (SVD) of $\mathcal{H}$.

Noting that the SVD of $\mathcal{H}$ induces a Fourier decomposition in the spatially invariant directions ($x$ and $z$, see Sharma et al. (2016a)), it is profitable to perform the SVD separately at each Fourier frequency-wavenumber combination $K=(\alpha,\ \beta)$ for the Fourier representation of (6),

[TABLE]

The SVD has the useful orthogonality properties

[TABLE]

In this case, the $\bar{\boldsymbol{u}}$ to be used in forming $\mathcal{L}$ is the $\alpha=0,\ \beta=0$ component of the temporal mean.

The velocity field may then be expressed as an expansion in resolvent modes,

[TABLE]

Applying the optimal rank-$M$ projection $\Pi_{M}$ gives

[TABLE]

with the sum over $m$ in (9) being split as

[TABLE]

The idea of the projection is that the flow resides mostly (in an energy sense) in the subspace that $\Pi_{M}$ projects onto. For an explanation in terms of Koopman modes, refer to Sharma et al. (2016a). The projection $\Pi_{M}$ is calculated for each known solution and then used as an initial guess for the NKH search.

3 Symmetries in plane Couette flow

A symmetry operation $\sigma$ is a linear transformation of the state of a dynamical system which commutes with integration forward in time,

[TABLE]

We define an isotropy group of $\tilde{\boldsymbol{u}}$ as a group that contains all symmetries that satisfy $\sigma\tilde{\boldsymbol{u}}=\tilde{\boldsymbol{u}}$.

Solutions in an equivariant system (such as plane Couette flow, which is highly symmetric) can satisfy all of the system’s symmetries, a subgroup of the symmetries or none of the symmetries. Typically, a turbulent trajectory has no symmetries, i.e. its isotropy group consists of the identity operation $\{e\}$ only. The laminar solution in plane Couette flow obeys every continuous symmetry that the geometry allows (see figure 2), this symmetry group is defined as $\Gamma$ (see Halcrow (2008) for full derivation).

For plane Couette flow, the NSE are invariant under reflection in the $yx$-plane, rotation about the $z$-axis by $\pi$, pointwise-inversion through the origin and continuous translations in the $x$ and $z$ axes. Accordingly isotropy groups of the equilibria reported here contain combinations of reflection, rotation, pointwise-inversion and translation, defined as

[TABLE]

respectively.

The symmetries that a particular solution adheres to dictate the type of solution it is. The reflection symmetry reverses the spanwise velocity, $w$, therefore any solution that is invariant under $\sigma_{z}$ cannot be a spanwise travelling wave. Similarly, the rotation symmetry reverses the streamwise velocity, $u$, therefore for any solution that is invariant under $\sigma_{x}$ does not permit a streamwise travelling wave. Consequently, if a solution is invariant under $\sigma_{xz}$ it cannot be a travelling wave solution in $x$ or $z$. This implies that the solutions that obey $\sigma_{xz}$ are equilibria since they are spatially static. All of the equilibria in this study satisfy $\sigma_{xz}$, since the equilibria are derived from the solutions found by Nagata (1990) and Gibson et al. (2009), who sought equilibria that obeyed $\sigma_{xz}$.

The periodic boundary conditions impose discrete translation symmetries on the equilibria. If a field is fixed under a discrete shift $\tau(L_{x}/n,0)$, it is periodic on the smaller spatial domain $x\in[0,L_{x}/n],n\in\mathbb{Z}^{+}$, similarly for $z$. Using half-cell shifts in the streamwise ($\varDelta x=L_{x}/2$) and spanwise ($\varDelta z=L_{z}/2$) directions only, the following symmetry operations can be defined

[TABLE]

These operations are used to define the isotropy groups to which all solutions in our work belong. For more details on the other isotropy groups of plane Couette flow see Halcrow (2008). All equilibrium solutions found in the present study belong to the following one of the following isotropy subgroups,

[TABLE]

Here $\Theta_{3}\subset\Sigma,K\subset\Theta\subset\Gamma$. The isotropy subgroup $K$ only admits translations in the spanwise direction; an equilibrium solution that belongs to this symmetry subgroup was discovered and is detailed in Section 4.1. The isotropy subgroup $\Sigma$ is called $S$ in the works of Nagata (1990), Gibson et al. (2009) and Waleffe (2003). $\Theta$ can also be expressed as $\Theta=\Sigma\times\{e,\tau_{xz}\}$. These are only some of the symmetry subgroups of $\Gamma$. Other subgroups might also play an important role in turbulent dynamics; there may be other equilibria that obey different symmetries, or none at all.

The full space of velocity fields for plane Couette flow can be written as

[TABLE]

The isotropy subgroups can be used to define invariant subspaces of velocity fields that are invariant under symmetry operations contained within the isotropy subgroups. For an arbitrary isotropy subgroup $H$, the $H$-invariant subspace can be defined as

[TABLE]

where $h$ is a symmetry operation. $H$ can be replaced with any of the isotropy subgroups defined in (17).

3.1 Geometry of plane Couette state space

Plane Couette flow dynamics can be represented as state space portraits. Gibson et al. (2008) developed a basis that defines a representation independent state space which shows the relationships between equilibria and allows us to chart turbulent trajectiories. A brief summary of the state space visualisation method developed by Gibson et al. is given below.

The premise of the method is that velocity fields may be projected onto an orthonormal basis,

[TABLE]

where $\boldsymbol{e}_{n}$ is a unit basis vector and $a_{n}$ is the low-dimensional projection of a given velocity field $\boldsymbol{u}(t)$. Therefore, a state space trajectory is projected onto the $\{\boldsymbol{e}_{n}\}$ coordinate frame.

Following the work of Gibson et al. (2008), an orthonormal translational basis is constructed based on streamwise and spanwise half-domain shifts of Nagata (1990)’s upper branch solution ($\boldsymbol{u}_{\mathrm{EQ2}}$), defined as

[TABLE]

where $\gamma_{n}$ is a normalisation constant such that $\|\boldsymbol{e}_{n}\|=1$. Here, $\tau_{i}$ represents a half-domain shift in the direction specified by the subscript $i$. On the right hand side the last three columns denote the symmetry of each basis function under the appropriate translation, e.g $A$ in the $\tau_{x}$ column means that $\tau_{x}\boldsymbol{e}_{n}=-\boldsymbol{e}_{n}$ (anti-symmetric) and $S$ means $\tau_{x}\boldsymbol{e}_{n}=\boldsymbol{e}_{n}$ (symmetric). The origin in this state space is the laminar solution $\boldsymbol{u}_{\mathrm{EQ0}}$ since it is invariant under all symmetries and all solutions presented here are expressed as are differences from laminar. As emphasised by Gibson et al. (2008), this orthonormal basis definition is one of many.

It should be noted that a basis can be constructed from any field as there is no pre-determined method of selecting a fluid state, being at the author’s discretion. Following Gibson et al. (2008), we chose to select an orthonormal basis formed from $\boldsymbol{u}_{\mathrm{EQ2}}$.

4 Results

Firstly, the energy and symmetry properties of all equilibria collated in table 1 are described and compared in §4.1. Secondly, we visualise the new equilibria and their invariant manifolds in state space (defined in (21)) in §4.2. Thirdly, the equilibria are continued in $Re$ and the resulting bifurcation curves are analysed in §4.3. For comparison with the the work of Gibson et al. (2009), projections were only performed at three distinct Reynolds numbers of $400$, $330$ and $270$. The results are labelled accordingly.

4.1 Equilibria catalogue

The properties of all previously known and new equilibria are given in table 1, where the solutions are organised by isotropy subgroup and sorted by descending dissipation rate, $D$, within. In the following sections equilibria flow fields will be labelled as EQX rather than $\boldsymbol{u}_{EQX}$ for clarity. The equilibria labels follow the naming convention of Gibson et al. (2009); EQ0 is the laminar solution, EQ1 and EQ2 are Nagata (1990)’s lower- and upper-branch equilibria, respectively, and EQ3-EQ11 are the equilibria discovered by Gibson et al. (2009) and Halcrow et al. (2009). The rest of the equilibria reported here are new and the numeric labelling denotes the chronological order in which they were found. In the following, a projection of integer rank $M$ is denoted $\Pi_{M}(\cdot)$, where $1\leq M\leq 99$.

The unstable parts of the new equilibria and their parents’ spectra are shown in figure 3. From this figure we can get an idea of the nature of their unstable manifolds computed using the leading unstable eigenvalues and eigenvectors. For example, the leading unstable eigenvalues of EQ11 are a complex conjugate pair, so we expect the trajectories, perturbed along the leading unstable eigenvector of EQ11, to travel on a surface away from the fixed point. Whereas, for EQ24 the leading unstable eigenvalue is real only, hence we expect all perturbed trajectories to travel along a one-dimensional curve. For figures depicting these unstable manifolds, see §4.2.

4.1.1 Child equilibria found at $Re=400$

EQ20, EQ21, EQ23.

These equilibria were discovered from projections of EQ9 and EQ10. The new solutions obey the same symmetries as the parent equilibria and all three of them resemble EQ10; the fields contain a streamwise streak flanked by a pair of counter-rotating asymmetric vortices. The main difference between the flow fields is that the child equilibria have smoother structures with lower velocity magnitudes. The unstable manifolds of the child solutions have reduced dimensionality in $\mathbb{U}$ and within the $\Theta_{3}$-invariant subspace $\mathbb{U}_{\Theta 3}$, as compared to the parent solutions. The magnitudes of the unstable eigenvalues of the children are reduced compared to those of their parents, see figure 3(a).

EQ22, EQ24

These equilibria were discovered from projections of EQ11. These equilibria were also found by Gibson et al. (2009) by continuation of EQ10 with varying spanwise wavenumber. EQ22 resembles EQ10 and EQ24 is similar in appearance to EQ11 and both children and parents adhere to the $\Theta_{3}$ isotropy subgroup. Both EQ22 and EQ24 have reduced dimensionality in $\mathbb{U}_{\Theta 3}$. However, EQ24 is more unstable than EQ11 in $\mathbb{U}$ and further from the laminar state, so presumably it is more likely to be approached by a turbulent trajectory. The unstable eigenvalues of the new equilibria and EQ11 have similar real parts in both $\mathbb{U}$ and $\mathbb{U}_{\Theta 3}$, see figure 3(a).

4.1.2 Child equilibria found at $Re=330$

All of the following equilibria were discovered from projections of EQ6.

$\boldsymbol{\mathrm{EQ12}}$.

This solution sits on the lower branch of the EQ7-8 bifurcation curve at $Re=330$, see figure 7. The solution looks like EQ7 in terms of velocity field structure and also obeys the same symmetries as the solutions on the EQ7-8 branch (all solutions on the curve belong adhere to the $\Theta$-isotropy subgroup). Of all the equilibria in table 1, it has one of the lowest energy norms and dissipation rates, both properties being close to the laminar solution. This is not surprising as the solution is on the EQ7 branch, which is the closest solution to the laminar state.

$\boldsymbol{\mathrm{EQ13}}$.

This solution sits on Nagata’s lower branch at $Re=330$ (see figure 8), accordingly it resembles and has the same symmetries as EQ1. Similar to EQ1, the new solution has a one dimensional unstable manifold within in the $\Sigma$-invariant subspace and the full space, and the magnitude of their unstable eigenvalues are similar.

$\boldsymbol{\mathrm{EQ14}}$.

This equilibrium solution sits on the lower branch of the EQ5-6 bifurcation curve at $Re=330$, see figure 9. EQ14 forms a lower-upper branch pair with EQ6 at $Re=330$, consequently the solution looks similar to and belongs to the same isotropy subgroup as its upper-branch partner. EQ14’s unstable manifold has the same dimensionality in the $\Sigma$-invariant subspace as EQ6.

$\boldsymbol{\mathrm{EQ16}},\boldsymbol{\mathrm{EQ26}}$.

These solutions are located on a new branch shown in figure 7. The flow field structures of the new solutions are similar to EQ6 with two fast-travelling streamwise streaks and two slow-travelling streamwise streaks above, arranged in an alternating fashion in the spanwise direction. The child equilibria belong to a different isotropy subgroup than their parent solution; the structures within EQ16 and EQ26 are asymmetric with respect to rotation and reflection symmetries, hence they belong to $\Theta_{3}$. Both solutions have high dimensional unstable manifolds in $\mathbb{U}_{\Theta 3}$, and their unstable manifolds have even greater dimensionality than their parent in $\mathbb{U}$. The magnitudes of the unstable eigenvalues are also higher than other solutions discovered at $Re=330$, see figure 3(b).

$\boldsymbol{\mathrm{EQ15}}$.

This equilibrium solution inhabits a newly discovered branch. It is similar to EQ6 in terms of flow field structures, however its symmetry subgroup is different; the solution is part of the $\Theta_{3}$-isotropy subgroup, hence it is less symmetric than its parent. EQ15 is highly dissipative and very unstable; its unstable manifold has greater dimensionality than EQ6 in the full space and the symmetry-invariant subspace. It is the most unstable solution found at $Re=330$, both its energy norms and dissipation rates are close to those of EQ6.

$\boldsymbol{\mathrm{EQ25}}$.

This solution sits on the lower branch of the EQ20-21-23 curve at $Re=330$. In terms of field structure the solution is identical to EQ20, except it is reflected in the $yx$-plane, i.e. EQ25 $=\sigma_{z}$EQ20. It also has the same number of unstable eigenvalues as EQ20 in the full space and symmetry-invariant subspace, and also belongs to the same isotropy subgroup as EQ20.

4.1.3 Child equilibria found at $Re=270$

All of the following equilibria were discovered from projections of EQ8.

$\boldsymbol{\mathrm{EQ17}}$.

This solution sits on the lower branch of the EQ7-8 bifurcation curve at $Re=270$ and shares much of its properties with EQ7. The solution looks like EQ7 in terms of velocity field structure and obeys the same symmetries as all solutions on the EQ7-8 bifurcation curve.

$\boldsymbol{\mathrm{EQ18}}$.

This solution is situated on the lower branch of EQ1-2 bifurcation curve at $Re=270$, it has the same symmetries as and looks similar to EQ1. The new equilibrium also has the same number and similar magnitude of unstable eigenvalues as EQ1.

$\boldsymbol{\mathrm{EQ19}}$.

This solution is located on a new branch shown in figure 7. It looks like an asymmetric version of EQ2 and is very dissipative. The solution does not have an unstable manifold in the symmetry-invariant subspace, which has not been seen before for any type of solution, and in the full space it has many unstable eigenvalues and hence a high dimensional unstable manifold. EQ19 is part of a symmetry subgroup previous solutions have not been observed to belong to. It is labelled as $K$ in table 1, defined as $K=\{e,\theta_{1},\theta_{2},\theta_{3}\}$ (see Section 3 for more details). The isotropy subgroup $K$ is unlike other isotropy subgroups as it allows the shift-reflect and shift-rotate symmetries given in (16) but only permits translations in the spanwise direction.

$\boldsymbol{\mathrm{EQ27}}$.

This solution, like EQ18 sits on the lower branch of the EQ1-2 bifurcation curve, however when continued in $Re$ it follows the path of the EQ20-21-23 bifurcation curve and is very similar to EQ25. Like EQ18, this is an interesting solution being close to the bifurcation points of many difference branches in figure 7.

4.2 Visualising equilibria in state space

In order to understand the inter-relations between equilibria and to study their manifolds we must visualise the solutions in an appropriately reduced state space. To this end, we use the representation independent state space projections described in §3.1. Figures 4, 5 and 6 show the equilibria in state space at the distinct Reynolds numbers of $400$, $330$ and $270$, respectively. The different regions of the figures have labels in the corners to denote where translated siblings of the equilibria are projected. For example, if we shift EQ3 in the streamwise direction it is reflected in the $a_{4}$ axis (see figure 4(a)); this behaviour is due to the symmetry properties of $\boldsymbol{e}_{4}$ as it is anti-symmetric with respect to streamwise shifts, see (21). The distances between equilibria in the state space portraits below relate to the differences in the magnitude of fluctuations, i.e. the $\mathscr{L}^{2}$-norms of the fixed points. Therefore, the further a fixed point is from the laminar solution (the origin) the greater its deviation from the laminar flow field.

The trajectories seen leaving the equilibria in figures 4, 5 and 6 depict the unstable manifolds of the equilibria. We use the method outlined by Gibson et al. (2008) to compute the unstable manifolds. In short, we perturb an equilibrium solution in the direction of its leading unstable eigenvector and integrate the perturbed field forward in time to determine the unstable manifold. Mathematically the perturbed field is defined as

[TABLE]

where $\boldsymbol{u}_{EQ}$ is an equilibrium solution, $\epsilon$ is a scaling constant and $\boldsymbol{v}_{EQ}$ is the unstable eigenvector. If the leading unstable eigenvalues of a solution are a complex conjugate pair, $\epsilon$ is varied in the range $10^{-2}\geq\epsilon\geq 10^{-4}$, whereas if the leading unstable eigenvalue is positive real only $\epsilon=10^{-4}$. Note that the unstable eigenvector selected in (22) belongs to the solution’s symmetry-invariant subspace, hence all trajectories eventually decay onto the laminar solution.

Focussing on figure 4 first, we can see that EQ20, EQ21 and EQ22 have filled in the region between EQ1 and EQ10. Those three new equilibria are close to the $a_{1}$ axis, which means that they do not vary as much when they are translated in the streamwise direction. They do, however, vary significantly in the spanwise direction as is evident from figure 4(e). The unstable manifolds of EQ20 and EQ21 form surfaces since the leading unstable eigenvalues are a complex conjugate pair (see figure 3) and their perturbed trajectories spiral outwards from the fixed points. However, there are some trajectories of EQ20 that are directed towards the laminar solution — those trajectories closely follow the manifold of EQ1 towards the laminar solution. EQ23 and EQ24 sit near EQ2 and EQ11 respectively, and their unstable manifolds are not as locally complicated as expected when looking at the dimensionality of their unstable manifolds in table 1; the two new manifolds closely follow the shape of EQ11’s manifold. Both EQ23 and EQ24’s unstable manifolds should form lines, however when perturbed in the negative eigendirection the trajectories swing back around to follow the dominant positive eigendirection.

Figure 5 shows that all solutions found from EQ6 are closer to the laminar solution than EQ6. This matches the trends seen in table 1, where all child equilibria of EQ6 have a lower $\|\tilde{\boldsymbol{u}}\|$ and lower dissipation rate. All solutions near EQ6 (including EQ6) are highly unstable with complicated manifold shapes that closely follow each other at certain points and generally resemble each other’s manifold shape. EQ12, EQ13 and EQ25 belong to the EQ7-8, EQ1-2 and EQ20-21-23 solution branches (see figure 7). Interestingly, the unstable manifolds of EQ13 and EQ25 are more unstable than the solutions that belong to their respective branches at $Re=400$, i.e. EQ1 and EQ20 in figure 4. The unstable manifolds push trajectories away from the laminar solution at $Re=330$, which is in contrast to the behaviour seen at $Re=400$. The manifold of EQ12 encases all of the other manifolds, but because of EQ12’s symmetry subgroup, its manifold is only visible in the plots with $a_{4}$ as one of the axes. The unstable manifolds of EQ6, EQ14, EQ15, EQ16 and EQ26 are quite complicated in shape and seem to be entangled with each other. EQ6 and EQ15’s manifolds initially follow each other very closely, eventually going their separate ways but maintaining a loose resemblence.

Figure 6 shows that EQ8’s child equilibria are closer to the laminar solution. EQ8’s leading unstable eigenvalues are a complex conjugate pair indicating that a surface will be formed in state space. However, the perturbed trajectories converge to almost one line that loops around and takes aim at the laminar solution. Hence any turbulent trajectory near EQ8 would quite quickly be guided towards the laminar solution. EQ17’s unstable manifold is similar to EQ8’s, except it has a larger turning radius. EQ18 is hidden behind EQ27, both solutions belong to a closely related pair of families, EQ18 to the EQ1-2 and EQ27 to the EQ20-21-23 branch, which originates from a pitchfork bifurcation from the EQ1-2 lower branch at approximately $Re=270$, see §4.3. The unstable manifolds of EQ18 and EQ27 are more unstable than other solutions in their respective families at higher Reynolds numbers since the perturbed trajectories stray quite far from their respective fixed points. The manifolds of EQ18 and EQ27 follow each other closely until they diverge after crossing into opposing quadrants. As noted in table 1, EQ19 does not have any unstable eigenvalues in the symmetry-invariant subspace, thus no manifold could be computed.

4.3 Continuation under $Re$

The solutions were parametrically continued with $Re$ as the bifurcation parameter and subsequently grouped into pairs by observing their bifurcation curves and identifying which part of the curve each solution inhabits. The naming convention for branch pairs is EQX-Y, where X indicates the ‘lower’ branch solution (lower dissipation rate) and Y signifies the ‘upper’ branch solution (higher dissipation rate). There is one case where three equilibria sit on one curve, in that case the naming convention is in ascending order of dissipation rate of each equilibrium solution.

The bifurcation diagrams in figures 7, 8, 9 and 10 are plotted with the dissipation rate $D$ against the bifurcation parameter $Re$. The curves are independent of each other and any apparent intersections between curves are not bifurcations, except in the case of the pitchfork bifurcations which are highlighted in the description below. Each curve represents a family of solutions with an upper and lower branch, originating from bifurcations at certain Reynolds numbers.

Figure 11 shows how the bifurcation curves vary with respect to kinetic energy density $E$ as well as $D$ and $Re$. The figure shows that solutions with high $D$ have low $E$ and the bifurcations seem to roughly form a pseudo-hyperplane in their $Re-D-E$ projection. We have already seen in the previous section that equilibria with large dissipation rates are typically further from the laminar solution. Since turbulence is highly dissipative, the less dissipative equilibria would seem to be more relevant in understanding transitional trajectories.

4.3.1 Bifurcation curves of equilibria found at $Re=400$

In ascending $Re$, the first saddle-node bifurcation encountered is at $Re=218.1$ which produces the Nagata (1990) EQ1-2 branches, respectively. The EQ20 branch originates from a pitchfork bifurcation at $Re=268.5$ from Nagata’s lower branch, see figure 8. There is a saddle-node bifurcation at $Re=270.2$ which forms the EQ21-23 curve that reverse bifurcates with the EQ20 branch at $Re=434.9$, forming the EQ20-21-23 curve, which merges with the EQ2 branch at $Re=980.4$.

The EQ10-11 branches merge with another bifurcation curve at $Re=609$ and almost form a loop, the curves remain disconnected between $Re\approx 635-655$, as shown in figure 8. When EQ22 and EQ24 are continued in $Re$, it is found that they inhabit the lower and upper branches of the curve that EQ10-11 merges with, respectively. It is also found that when continued, EQ22 and EQ24 close the aforementioned incomplete loop. There are no symmetry changes observed for any solution in the bifurcation loop and the flow field structure is similar across both sets of branches. The EQ10-11 branches come from a saddle-node bifurcation at $Re=348.2$, and the EQ22-24 branches originate from a saddle-node bifurcation at $Re=324.3$.

The EQ9 branch originates from a pitchfork bifurcation at $Re=369.3$ from the EQ3-4 curve as shown in figure 8. When continued up in Reynolds number, the EQ9 branch undergoes a reverse saddle-node bifurcation at $Re=712.8$ ($D=1.4$) followed by a saddle-node bifurcation at $Re=565.1$ ($D=1.54$), after which it smoothly continues to higher Reynolds numbers.

4.3.2 Bifurcation curves of equilibria found at $Re=330$

Some of the new solutions are situated on branches that were previously known, hence when the child equilibria are continued in Reynolds number it leads to the rediscovery of a branch that was already known. This can be seen in figures 8, 9 and 10 when symbols signifying different Reynolds numbers appear on the same curve; different symbols on the same curve are indicative of a branch-jump.

When EQ14 is continued down in Reynolds number, the bifurcation curve follows the path to EQ6 and continues further to form another branch. The curve increases to a dissipation rate of 3.2 at $Re=335.4$ and then rapidly but smoothly decreases in dissipation rate with decreasing Reynolds number towards a saddle-node bifurcation at $Re=312.2$. The two bifurcation curves merge at $Re=335.4$, note that there is no bifurcation here as there are no symmetry or stability changes at this point.

The EQ5-6 curve also has a pitchfork bifurcation that takes place at $Re=328$, leading to EQ15. The bifurcation is accompanied by a symmetry change from $\Sigma$ to $\Omega_{3}$ for all solutions that inhabit the new branch.

Lastly, the EQ16-26 bifurcation curve starts from a saddle-node bifurcation at $Re=298.1$. It consists of a lower-upper branch pair of solutions that cross in dissipation at $Re=453.7$, so that the lower and upper branches become the upper and lower branches, respectively. The upper branch then smoothly curves back on itself; the dissipation rate decreases followed by a slight increase as it loops back on itself while $Re$ increases. After that both $D$ and $Re$ decrease almost linearly, until the curve stops at $Re=478.1$, $D=2.1$. Focussing on the lower branch, the bifurcation curve then merges with another bifurcation curve at $Re=587.9$, which originates from a saddle-node at $Re=517.9$. The merged curve then continues upwards in Reynolds number until $Re=822.7$, at this point the curve loops and continues downwards in Reynolds number. There are three further loops where the symmetry of the solutions change from $\Theta$ to $K$ indicating that a bifurcation has taken place. The final part of the curve goes through a saddle-node bifurcation at $Re=436.1$, $D=3.61$ and continues upwards in $Re$ whilst maintaining a very low dissipation rate. There is a further merging with another bifurcation curve at $Re=1010.9$ that is spawned from a saddle-node bifurcation at $Re=943.3$. This curve does not exhibit any twists and turns beyond $Re=1200$ and note that in figure 9 it seems as though this curve joins with the EQ5-6 lower branch, however this is not the case as can be seen in figure 11. It should be noted that at $Re=478,D=2.1$ and $Re=450,D=3.4$ there are unstable fixed points that, when approached, cause the continuation algorithm to reverse direction and retrace the path it took to those points.

4.3.3 Bifurcation curves of equilibria found at $Re=270$

The EQ7-8 branches are formed from a saddle-node bifurcation at $Re=234.2$, see figure 10. At $D=5.2$ the curve merges with another bifurcation curve; this is evident from a rapid and smooth continuation at $Re=364.2$. A lower-upper branch pair specifically at $Re=270$ can be formed with EQ17-8, and EQ7 can be paired with the corresponding solution at $Re=400$ on the upper branch.

The EQ19 curve originates from a pitchfork bifurcation at $Re=259.6$. The symmetry of the equilibria change at the bifurcation point, going from EQ19 to EQ7 the symmetry subgroup changes from $K$ to $\Theta$. The curve reaches a maximum dissipation at $D=3.71$ at $Re=433$, then experiences a slight decrease in dissipation to $D=3.66$ at $Re=435.2$.

As mentioned before, EQ18 and EQ27 are points on the bifurcation curves for the EQ1-2 and EQ20-21-23 branches, respectively. Hence, when the new equilibria are continued in Reynolds number their curves are identical to the aforementioned parent branches. EQ12 (discovered at $Re=330$) is the result of a branch-jump from the EQ5-6 curve to the EQ7-8 curve, see figure 10.

5 Discussion

The results in §4 reveal that, in most cases, equilibria with less fluctuation energy and lower dissipation (compared to their parent solutions) can be found using the project-then-search method. We speculate that the project-then-search technique allows us to take known equilibria on higher-dimensional unstable manifolds and use their projections (which are close to lower-dimensional manifolds) to find equilibria that sit on lower-dimensional unstable manifolds. The low-dimensional nature of the new equilibria is attributed to the projection generated by the resolvent model; resolvent modes for a given equilibrium solution span a low-dimensional space which the solution approximately inhabits. This is evident when comparing the dimensionality of the unstable manifolds of child equilibria against their parent solutions; almost all second/third generation solutions have an unstable manifold with reduced dimensionality in both the symmetry-invariant subspace and the full space. Also, it is clear from table 1 that there is a relationship between a solutions’ dissipation rate, energy and the dimensionality of its unstable manifold in both the symmetry-invariant subspace and the full space. In general, the higher the dissipation rate and fluctuation magnitude are, the more unstable the solution’s manifold is.

To illustrate the mechanics of the project-then-search method, the process can be shown in state space. For the purposes of demonstration, successful searches from EQ10’s low-rank projections are shown in figure 12. The search algorithm is quick to locate the region where a potential solution exists. Once the search is very close to a solution, it takes just a few Newton steps (not shown) to converge onto the solution. It is clear from figure 12 that small differences in the initial state used for the NKH algorithm can lead to very different results, indicating that the initial conditions are close to a basin boundary. For example $\Pi_{3}$(EQ10) and $\Pi_{4}$(EQ10) are extremely close to each other in state space, yet they lead to two different equilibria.

Visualising the project-then-search process in the reduced state space allows us to see the distance and direction of the projections relative to the original solution. Figure 12 shows that lowering the rank produces projections that are simultaneously further form the original solution and closer to the laminar state, in an energy sense. This is to be expected as the resolvent modes used to generate the projections are derived with respect to an energy norm. As the resolvent modes are truncated during the projection process, the small length-scale content of the original field is removed and energetic large-scale motions remain. The resulting fields of varying rank are ‘smoothed’ versions of the original solution, thus, the projected fields tend to be less dissipative. The projected fields physically resemble the original fields but with lower fluctuation magnitude (see figure 13). Consequently, less dissipative equilibria are discovered, with EQ20 and EQ27 being exceptions; they have a higher dissipation rate and energy magnitude than their parent. By way of example, the flow fields of the projections and final solutions from figure 12 are shown in figure 13. The plots clearly show that the projections remove small length-scales near the wall while maintaining the larger structures in the central region of the domain.

For every solution all available ranks were used to generate projections. The ‘Root’ column in table 1 denotes the lowest ranked-projection that leads to the discovery of a solution, but this does not mean that that is the only projection that yields that solution. All ranked-projections of the equilibria lead to the discovery of other non-trivial unstable solutions, although there are some equilibria whose projections only converge to the laminar state or return to the original fixed point (this fact is denoted in table 1 under the ‘Ret.’ column). Examples of such equilibria include EQ1, EQ3, EQ4, EQ7 and EQ20. There are, however, some projections that fail to initialise any successful search at all; there were a total of 107 failed searches, which translates to a $4\%$ failure rate. For these cases the search was repeated with eight times as many Newton steps ($160$ rather than $20$), yet the solutions failed to converge; Viswanath (2007) found that for a system with $10^{5}-10^{6}$ unknowns, the NKH algorithm takes less than $100$ steps to find an exact solution. In relation to (2), the required tolerance for successful searches is $\|G\|\sim 10^{-15}$, in contrast we find that for unsuccessful searches $10^{-6}\lesssim\|G\|\lesssim 10^{-3}$.

To demonstrate what unsuccessful searches look like in state space we have plotted the first 40 Newton steps of searches initiated from $\Pi_{38,40,43}($EQ10$)$ in figure 14. The figure shows that the Newton steps oscillate as they near a point, suggesting that the search is stuck near or between manifolds. Another potential reason for failure to converge may be the prescribed time horizon of the DNS phase when checking (2); in our study a moderate time horizon of $T=20$ was chosen in keeping with Gibson et al. (2009) who chose it to keep the computation time relatively short. Increasing the time horizon to a much larger value could help in converging the unsuccessful searches and potentially lead to the discovery of other solutions. Alternatively, the Newton step size could be reduced to prevent overshooting.

Focussing on the symmetries of the solutions reveals how spatially organised they are. That is to say, a highly symmetric solution has a very systematic recurrent structure, whereas a field which lacks symmetry, such as a turbulent field, lacks exactly repeating structure. Therefore, highly symmetric solutions are suitable for describing the low-dimensional nature of coherent structures. The EQ7-8 family of solutions adhere to the $\Theta$ isotropy subgroup: they are the most symmetric set of equilibria found and as shown in §4.2 their unstable manifolds are relatively simple and quickly aim trajectories towards the highly symmetric laminar state. Flow fields with fewer symmetry constraints can contain more complicated repeating structure. Equilibria found previously by Nagata (1990), Gibson et al. (2009), Clever & Busse (1997) and Waleffe (1998, 2003) in plane Couette flow conform to $\Sigma$ or $\Theta_{3}$. Most of the equilibria discovered in our investigation also adhere to those two isotropy subgroups. However, EQ19 adheres to the subgroup $K$ which exhibits periodicity in the spanwise direction only. To our knowledge this is the first time an equilibrium solution in plane Couette flow belonging to this symmetry subgroup has been found. Interestingly, $K$ is found from a pitchfork bifurcation off the EQ7-8 branch. As pointed out by Crawford & Knobloch (1991), when an equilibrium solution branch undergoes a bifurcation, states on the new branch have less symmetry and generally have more complicated dynamics. In the present work, all bifurcations lead to symmetry-breaking and the loss of symmetry is manifested by the appearance of a new pattern.

Furthermore, performing the project-then-search method at other Reynolds numbers may find more unknown branches. Highly symmetric equilibria along the upper-branches of bifurcation curves have the potential to lead to other solution branches because there are symmetries to break. For example, solutions along the upper-branch of the EQ1-2 bifurcation curve adhere to the isotropy subgroup $\Sigma$ which contains more than one symmetry operation. Therefore, conducting the project-then-search process at other Reynolds numbers on this branch may lead to the discovery of other symmetry-breaking bifurcations. The project-then-search method can also be applied to branches of solutions that adhere to the $\Theta_{3}$ isotropy subgroup of order 1. Doing so may lead to branches of chaotic solutions that do not obey any symmetries.

Presenting the equilibria in the representation-independent state space in §4.2 reveals that the spanwise-shifted equilibria are projected further from the original solution. Figure 4 shows EQ1’s unstable manifold and it is clear that a spanwise-shift leads to a larger deviation from its original state: EQ1 and $\tau_{x}$EQ1 are quite close in state space, whereas EQ1 and $\tau_{z}$EQ1 are quite far apart. This suggests that the structure within the field varies more in the spanwise direction. This property is observed in solutions that are far from the laminar state, and EQ1, EQ20, EQ21 and EQ22. Hence, transitional trajectories that come near the equilibria will be influenced more by the spanwise structure of the equilibria, rather than the streamwise structure. We argue that the spanwise-shifts play an important role in guiding and shaping transitional trajectories because the state space portraits show that the spanwise-shifts lead to larger deviations in state space. This is reminiscent of the SSP in which streamwise rolls, which redistribute the mean streamwise-momentum to sustain streaks, experience instability due to ‘spanwise inflections’ that lead to streamwise undulations, ultimately leading back to streamwise rolls via nonlinear self-interaction of the streamwise undulations (Waleffe, 1997). It is the spanwise inflections that lead to an instability that causes the streamwise-streaks to breakdown to recreate streamwise rolls. The breakdown and reformation of coherent structures in the turbulent boundary layer is referred to as bursting in the work of Viswanath (2007); this is an important process that is linked to the onset and maintenance of turbulent flow. Figure 6 shows that EQ18 and EQ27, which belong to the EQ1-2 and EQ20-21-23 solution families respectively, exhibit the same large deviation in the spanwise direction. Therefore, we suppose that these equilibria may also play an important role in transition and low-Reynolds number dynamics.

The state space portraits also reveal the shape of the laminar attractor. To our knowledge there are no solutions that have been found to inhabit an unstable manifold that encases EQ1’s manifold - this is the most spanwise varying manifold near the laminar solution, closely followed by EQ20.

6 Conclusions

A systematic exploration of the exact invariant solutions of the NSE was undertaken. New sets of equilibria were added to the collection of known exact coherent states found by Nagata (1990); Gibson et al. (2008, 2009) and Halcrow et al. (2009). The equilibria add more structure to the state space of plane Couette flow, and the computation of their unstable manifolds illuminates the pathways that exist in state space. These pathways can be used to find heteroclinic connections between equilibria and guide turbulent trajectories to desired physical states.

The key methodological advance reported here is the computationally cheap method by which initial guesses are generated for the NKH search. The resolvent model was used to generate low-rank projections of known equilibria which were then used as seeds for the NKH algorithm with a convergence rate of $96\%$. If the search succeeded in finding new solutions, then the project-then-search method was applied to the new equilibria, and so on until we obtained a closed set and no new equilibria were discovered. It is our belief that the project-then-search method’s success in finding new solutions is due to the combination of the projections’ properties and their location relative to other equilibria in state space. The energetically dominant physical characteristics of known solutions were maintained in the projections due to the nature of the projection. Searches initiated with the projections produced new solutions in the state space neighbourhoods of known solutions since the projections sit near the unstable manifolds of other solutions and the DNS phase of the NKH algorithm allows them to follow the directions of nearby manifolds. The low computational cost of the projection process makes it effective in obtaining numerous initial states, and thus solutions, from a few previously known ones. There are countless equilibria to be found; while we restricted ourselves to projections at three previously determined Reynolds numbers we fully expect that the same process at other Reynolds numbers would yield new branches.

In total, $4\%$ of searches initialised by projected flow fields failed to converge to the imposed tolerance despite extending the number of Newton steps. The failed searches seemed to be trapped in small regions of state space. This may be attributed to the small time horizon used for the DNS phase of the searches. This time horizon can be extended in future studies to avoid failed searches.

The current approach requires a known solution to find other solutions — no doubt this a significant limitation. Nonetheless, many avenues of further work suggest themselves. The projection method can be modified to truncate the Fourier modes in the homogeneous directions instead of (or as well as) the resolvent modes in the wall-normal direction; projections of equilibria at higher Reynolds numbers can be used to initiate searches; or a different projection basis can be used (for example, resolvent modes defined via a different norm).

The project-then-search approach can trivially be extended to the other classes of invariant solution. The case of periodic orbits will be presented in a later paper.

This work has been supported by the Air Force Office of Scientific Research (European Office of Aerospace Research and Development) under the award FA9550-14-1-0042. We would like to thank John Gibson for providing his code, equilibrium solutions and helpful comments. We would also like to acknowledge Edgar Knobloch for his guidance and Davide Lasagna for his probing questions and useful recommendations. This work was completed in part at the Kavli Institute for Theoretical Physics, with the support of the National Science Foundation under Grant No. NSF PHY11-25915.