2.1 Problem definition

Figure 1. Geometry of our study. At the inlet ( $BC_{inlet}$ ) a uniform unit velocity ( $u=1;v=0$ ) is imposed. Dashed red line: free-slip condition ( $BC_{free\text{-}slip}$ ). Thick blue line: no-slip boundary condition ( $BC_{no\text{-}slip}$ ). A free-outflow boundary condition is imposed at the outlet.

The configuration considered is the two-dimensional incompressible viscous shear-driven flow of a Newtonian fluid over an open cavity with equal length and depth shown in figure 1. This configuration is the same as that considered by Sipp & Lebedev (Reference Sipp and Lebedev2007) and Barbagallo et al. (Reference Barbagallo, Sipp and Schmid2009), or more recently by Meliga (Reference Meliga2017). We use the unperturbed upstream velocity $U_{\infty }$ , the cavity length $L$ and the resulting advective time $L/U_{\infty }$ to non-dimensionalize the variables. The dynamics of the flow is governed by the incompressible Navier–Stokes equations

(2.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}P+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{U}\otimes \boldsymbol{U}),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{U}(\boldsymbol{x},t)=(U,V)^{\text{T}}$ and $P$ are the velocity and pressure fields. The Reynolds number $Re$ is defined as

(2.2) $$\begin{eqnarray}\displaystyle Re=\frac{U_{\infty }L}{\unicode[STIX]{x1D708}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the fluid, and will range between $Re=4000$ and $5000$ . The boundary conditions, illustrated in figure 1, are

(2.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{U}=\boldsymbol{e}_{\boldsymbol{x}}\quad \text{on }BC_{inlet},\\ \boldsymbol{U}=\mathbf{0}\quad \text{on }BC_{no\text{-}slip},\\ \unicode[STIX]{x2202}_{y}U=V=0\quad \text{on }BC_{free\text{-}slip},\\ \unicode[STIX]{x2202}_{x}\boldsymbol{U}=\mathbf{0}\quad \text{on }BC_{outlet}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The boundary conditions at the inlet and along the wall are crucial. The flow is given a uniform profile at the inlet and develops a boundary layer structure as it advances downstream. The instability occurs where the boundary layer reaches the upstream corner of the cavity and it is the thickness of the boundary layer at this point that controls the details of the transition. When free-slip conditions are imposed on the wall close to the inlet, then a boundary layer of an appropriate thickness develops over a shorter distance than would be the case if no-slip conditions were used over the entire wall. A shorter domain can be used, making the calculation more economical.

The Navier–Stokes equations are solved using the incompressible flow solver Nek5000 (Fischer, Lottes & Kerkemeir Reference Fischer, Lottes and Kerkemeir2008) which is based on the spectral element method. A $\mathbb{P}_{N}-\mathbb{P}_{N-2}$ formulation has been used: the velocity field is discretized using $N$ th-order Lagrange interpolants defined on the Gauss–Legendre–Lobatto quadrature points as basis and trial functions while the pressure field is discretized using Lagrange interpolants of degree $N-2$ defined on the Gauss–Legendre quadrature points. Finally, time integration is performed using the BDF3/EXT3 scheme: integration of the viscous term relies on backward differentiation while the convective terms are integrated explicitly using a third-order accurate extrapolation. In practice, the polynomial degree was set to $N=6$ while the computational domain was discretized using 4000 spectral elements. The resulting mesh refinement is thus similar to that used in Sipp & Lebedev (Reference Sipp and Lebedev2007).

2.2 Base flow and linearization

A base flow $\boldsymbol{U}_{b}(\boldsymbol{x})$ is a solution of the stationary Navier–Stokes equations

(2.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle 0=-\unicode[STIX]{x1D735}P_{b}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}_{b}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{U}_{b}\otimes \boldsymbol{U}_{b}),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}_{b},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with the boundary conditions again given by (2.3). Various techniques can be used to compute the base flow $\boldsymbol{U}_{b}(\boldsymbol{x})$ . Because of its simplicity, the selective frequency damping (SFD) technique initially proposed by Åkervik et al. (Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006) has been used; see also Jordi, Cotter & Sherwin (Reference Jordi, Cotter and Sherwin2014, Reference Jordi, Cotter and Sherwin2015), Cunha, Passaggia & Lazareff (Reference Cunha, Passaggia and Lazareff2015).

Once the equilibrium $\boldsymbol{U}_{b}(\boldsymbol{x})$ has been computed, we determine its linear stability. To do so, we consider an infinitesimal perturbation $\boldsymbol{u}(\boldsymbol{x},t)$ to the base flow $\boldsymbol{U}_{b}$ , whose dynamics is governed by the linearized Navier–Stokes equations

(2.5) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{U}_{b}\otimes \boldsymbol{u}+\boldsymbol{u}\otimes \boldsymbol{U}_{b}),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The boundary conditions are the homogeneous version of (2.3), i.e. we now prescribe a zero velocity profile at the inlet.

(2.6) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{u}=\mathbf{0}\quad \text{on }BC_{inlet},\\ \boldsymbol{u}=\mathbf{0}\quad \text{on }BC_{no\text{-}slip},\\ \unicode[STIX]{x2202}_{y}u=v=0\quad \text{on }BC_{free\text{-}slip},\\ \unicode[STIX]{x2202}_{x}\boldsymbol{u}=\mathbf{0}\quad \text{on }BC_{outlet}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Solutions to (2.5)–(2.6) are of the form $\boldsymbol{u}(\boldsymbol{x},t)=\hat{\boldsymbol{u}}(\boldsymbol{x})\text{e}^{(\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714})t}+\text{c.c.}$ , $p(\boldsymbol{x},t)=\hat{p}(\boldsymbol{x})\text{e}^{(\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714})t}+\text{c.c.}$ , from which we obtain the eigenvalue problem

(2.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}(\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714})\hat{\boldsymbol{u}}={\mathcal{L}}_{\boldsymbol{U}_{b}}\hat{\boldsymbol{u}},\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}_{\boldsymbol{U}_{b}}$ is the Jacobian of the Navier–Stokes equations linearized around $\boldsymbol{U}_{b}$ ,

(2.8) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{\boldsymbol{U}_{b}}\hat{\boldsymbol{u}}\equiv -\unicode[STIX]{x1D735}\hat{p}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\hat{\boldsymbol{u}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{U}_{b}\otimes \hat{\boldsymbol{u}}+\hat{\boldsymbol{u}}\otimes \boldsymbol{U}_{b}). & & \displaystyle\end{eqnarray}$$

The stability of the base flow is determined by the sign of the real part $\unicode[STIX]{x1D70E}$ of the leading eigenvalue, which is the growth rate of the perturbation. If $\unicode[STIX]{x1D70E}$ crosses zero for an eigenvalue with non-zero imaginary part $\unicode[STIX]{x1D714}$ , then a Hopf bifurcation leads to a limit cycle whose frequency at onset is $\unicode[STIX]{x1D714}$ . In our case, the base flow undergoes a first Hopf bifurcation at $Re_{2}=4126$ , leading to a limit cycle $LC_{2}$ and a second Hopf bifurcation at $Re_{3}=4348$ gives rise to $LC_{3}$ .

We computed the leading eigenvalues and eigenvectors using a time-stepper approach; see, e.g. Edwards et al. (Reference Edwards, Tuckerman, Friesner and Sorensen1994). Our stability calculation typically used a Krylov subspace of dimension $K=256$ and a sampling period $\unicode[STIX]{x0394}T=10^{-3}$ non-dimensional time units. Eigenvalues were considered to be converged if the residual obtained from the Arnoldi decomposition was below $10^{-6}$ .

2.3 Mean flow and linearization

At the threshold of a Hopf bifurcation, linearization about the base flow leads to an eigenvalue whose real part is zero and whose imaginary part is the frequency of the limit cycle which is produced. As the Reynolds number is increased and the limit cycle develops nonlinearly and deviates from the base flow, eigenvalues obtained by linearization about the base flow no longer correspond to the properties of the limit cycle. However, linearization about the mean flow often leads to an eigenvalue whose imaginary part is closer to the nonlinear frequency.

We consider a Reynolds decomposition of the instantaneous flow field, i.e. 

(2.9) $$\begin{eqnarray}\displaystyle \boldsymbol{U}(\boldsymbol{x},t)=\overline{\boldsymbol{U}}(\boldsymbol{x})+\boldsymbol{u}(\boldsymbol{x},t), & & \displaystyle\end{eqnarray}$$

where $\overline{\boldsymbol{U}}(\boldsymbol{x})$ is the mean flow and $\boldsymbol{u}(\boldsymbol{x},t)$ is the zero-mean fluctuation. Introducing this decomposition into the Navier–Stokes equations and averaging shows that $\overline{\boldsymbol{U}}$ is governed by the Reynolds averaged Navier–Stokes (RANS) equations

(2.10) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle 0=-\unicode[STIX]{x1D735}\overline{P}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\overline{\boldsymbol{U}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\boldsymbol{U}}\otimes \overline{\boldsymbol{U}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\boldsymbol{u}\otimes \boldsymbol{u}}),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\boldsymbol{U}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with inhomogeneous boundary conditions (2.3). The presence of the Reynolds stress tensor $\overline{\boldsymbol{u}\otimes \boldsymbol{u}}$ of the fluctuation means that these equations are not closed. We compute the mean flow $\overline{\boldsymbol{U}}(\boldsymbol{x})$ of a period- $T$ limit cycle by carrying out a full nonlinear simulation via (2.1), (2.3) and time averaging,

(2.11) $$\begin{eqnarray}\displaystyle \displaystyle \overline{\boldsymbol{U}}(\boldsymbol{x})=\frac{1}{T}\int _{0}^{T}\boldsymbol{U}(\boldsymbol{x},t)\,\text{d}t. & & \displaystyle\end{eqnarray}$$

The equations governing the dynamics of the fluctuation $\boldsymbol{u}(\boldsymbol{x},t)$ are obtained by substituting $\overline{\boldsymbol{U}}+\boldsymbol{u}$ for $\boldsymbol{U}$ into (2.1) and subtracting (2.10), leading to

(2.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\boldsymbol{U}}\otimes \boldsymbol{u}+\boldsymbol{u}\otimes \overline{\boldsymbol{U}})-\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{u}-\overline{\boldsymbol{u}\otimes \boldsymbol{u}})}_{\boldsymbol{f}},\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with homogeneous boundary conditions (2.6). The equations differ from the linearized Navier–Stokes equations by the presence of $\boldsymbol{f}$ . The term $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{u})$ , is the usual quadratic interaction term neglected in base flow linear stability analyses. The term $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\boldsymbol{u}\otimes \boldsymbol{u}})$ is the divergence of the Reynolds stress tensor of the fluctuation.

Studies carrying out mean flow stability analyses discard $\boldsymbol{f}$ , leading to the linearized equations

(2.13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\otimes \overline{\boldsymbol{U}}+\overline{\boldsymbol{U}}\otimes \boldsymbol{u}),\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Using once again a normal mode ansatz, this set of equations is reduced to the eigenvalue problem

(2.14) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}(\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714})\hat{\boldsymbol{u}}={\mathcal{L}}_{\overline{\boldsymbol{U}}}\hat{\boldsymbol{u}},\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}_{\overline{\boldsymbol{U}}}$ is now the Navier–Stokes operator linearized around the mean flow, with $\overline{\boldsymbol{U}}$ substituted for $\boldsymbol{U}_{b}$ in (2.8).

Although the mean flow is not an equilibrium solution of the Navier–Stokes solution, this approach has proved unexpectedly successful in characterizing the frequencies of the full nonlinear solutions of the Navier–Stokes equations. (For a counter-example, however, see Turton et al. (Reference Turton, Tuckerman and Barkley2015).) Various theoretical overlapping justifications have been proposed for dropping or modelling $\boldsymbol{f}$ in (2.12), such as:

1. (i) The quadratic interaction of the fluctuation with itself is small and its temporal mean does not appear in a linear stability analysis (Barkley Reference Barkley2006; Mantič-Lugo et al. Reference Mantič-Lugo, Arratia and Gallaire2014, Reference Mantič-Lugo, Arratia and Gallaire2015).

2. (ii) The terms can be treated via an expansion in the distance from the threshold (Sipp & Lebedev Reference Sipp and Lebedev2007).

3. (iii) The instantaneous Reynolds stress tensor $\boldsymbol{u}\otimes \boldsymbol{u}$ is approximately equal to its temporal average $\overline{\boldsymbol{u}\otimes \boldsymbol{u}}$ so that they almost cancel out (Turton et al. Reference Turton, Tuckerman and Barkley2015); see also § 4.2.

4. (iv) The resolvent operator $(\text{i}\unicode[STIX]{x1D714}-{\mathcal{L}}_{\overline{\boldsymbol{U}}})^{-1}$ is sharply peaked or of low rank. This implies that its action on any vector, including $\boldsymbol{f}$ , is approximately a projection onto the leading mode of the resolvent. The purpose of this procedure is not to eliminate $\boldsymbol{f}$ , but to approximate it as a scalar multiple of the leading resolvent mode, and then to treat $\boldsymbol{f}$ as a forcing input (McKeon & Sharma Reference McKeon and Sharma2010; Hwang & Cossu Reference Hwang and Cossu2010; Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016; Symon et al. Reference Symon, Rosenberg, Dawson and McKeon2018).

2.4 Floquet analysis

Our study of the shear-driven cavity focuses on two limit cycles, denoted by $LC_{2}$ and $LC_{3}$ , created by primary Hopf bifurcations and destabilized via secondary bifurcations. Floquet analysis will be used to characterize this destabilization. The dynamics of an infinitesimal perturbation $\boldsymbol{u}(\boldsymbol{x},t)$ evolving in the vicinity of a $T$ -periodic limit cycle $\boldsymbol{U}(\boldsymbol{x},t)$ is governed by the linearized Navier–Stokes equations

(2.15) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}={\mathcal{L}}_{\boldsymbol{U}(t)}\boldsymbol{u},\\ 0=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{ u},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with homogeneous boundary conditions (2.6). This set of equations is non-autonomous, as the operator ${\mathcal{L}}_{\boldsymbol{U}(t)}$ is $T$ -periodic. Solutions to (2.15) are of the Floquet form

(2.16) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{x},t)=\hat{\boldsymbol{u}}(\boldsymbol{x},t)\text{e}^{(\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F})t}+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

where $\hat{\boldsymbol{u}}(\boldsymbol{x},t)$ are the $T$ -periodic Floquet modes and $(\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F})$ the Floquet exponents. The stability of $\boldsymbol{U}(t)$ is determined by the Floquet multipliers

(2.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}=\text{e}^{(\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F})T}. & & \displaystyle\end{eqnarray}$$

If the moduli of the Floquet multipliers are smaller than one, perturbations will decay exponentially fast and the orbit is stable. On the other hand, if at least one of the Floquet multipliers has a modulus greater than one, then that perturbation will grow exponentially and the orbit is unstable; see, e.g. Barkley & Henderson (Reference Barkley and Henderson1996) and Gioria et al. (Reference Gioria, Jabardo, Carmo and Meneghini2009). In our study the Floquet exponents are complex and the imaginary part $\unicode[STIX]{x1D714}_{F}$ of the Floquet exponent is the argument (angle) of the Floquet multiplier. The presence of an imaginary part leads to quasi-periodic behaviour. More details and results are shown in § 3.4.

2.5 Edge state technique for computing the unstable quasi-periodic state

As will be shown in § 3, there is a range of Reynolds numbers over which limit cycles $LC_{2}$ and $LC_{3}$ co-exist. In phase space, on the boundary between the basins of attraction of these limit cycles, is an unstable quasi-periodic state $QP$ . (More specifically, $QP$ is an edge state, meaning that within the boundary, trajectories are attracted to it.) In order to compute $QP$ , we use the same technique as in Itano & Toh (Reference Itano and Toh2001) or Duguet, Willis & Kerswell (Reference Duguet, Willis and Kerswell2008) for the laminar–turbulent edge state. In such cases, whether a trajectory evolves towards a turbulent or laminar state depends on the initial condition. Some initial conditions evolve directly to turbulence, others decay directly to the laminar state. By appropriately weighting turbulent and laminar solutions, an initial condition can be constructed so that the resulting trajectory converges to and remains a long time on the edge state before diverging towards one of these two attractors.

Figure 2. Time traces of streamwise velocity at $Re=4460$ for two simulations. The initial condition of the simulation is (2.18). For $\unicode[STIX]{x1D6FC}=0.47562027$ the system evolves towards $LC_{2}$ , shown as the black curve. For $\unicode[STIX]{x1D6FC}=0.47562256$ it evolves towards $LC_{3}$ , shown as the higher-amplitude red curve.

In our problem, a quasi-periodic state separates the two stable limit cycles $LC_{2}$ and $LC_{3}$ . Therefore we construct a weighted sum of the two, seeking an initial condition $\boldsymbol{U}(\boldsymbol{x},0)$ that will evolve after some time to $QP$ and then remain as long as possible before eventually converging to either limit cycle. Using the same bisection technique as in Lopez et al. (Reference Lopez, Welfert, Wu and Yalim2017), this initial condition is given by

(2.18) $$\begin{eqnarray}\displaystyle \boldsymbol{U}(\boldsymbol{x},t)=\unicode[STIX]{x1D6FC}\boldsymbol{U}_{LC_{2}}+(1-\unicode[STIX]{x1D6FC})\boldsymbol{U}_{LC_{3}}. & & \displaystyle\end{eqnarray}$$

For $\unicode[STIX]{x1D6FC}=1$ , the initial condition is $LC_{2}$ and for $\unicode[STIX]{x1D6FC}=0$ , it is $LC_{3}$ . For each Reynolds number considered, we successively delimit an interval of $\unicode[STIX]{x1D6FC}$ by bisection to capture the quasi-periodic state. As an illustration, we plot in figure 2 the time evolution of streamwise velocity recorded by a probe located at $(x_{1},y_{1})=(1.2,0.2)$ for $\unicode[STIX]{x1D6FC}=0.47562027$ and $0.47562256$ . A slight difference in $\unicode[STIX]{x1D6FC}$ will bring the system after a long transient regime to either $LC_{2}$ or $LC_{3}$ .

2.6 Standard deviation

Figure 3. (a) Time series of streamwise velocity at $(x_{1},y_{1})=(1.2,0.2)$ and $Re=4580$ for a simulation with $\unicode[STIX]{x1D6FC}=0.829630407714844$ . The regimes corresponding to $QP$ and to $LC_{3}$ are shown in (b,c). (d) The standard deviation is computed by a sliding window.

To construct the bifurcation diagram, we seek an appropriate measure of the oscillation amplitude as a function of $Re$ . Time series from limit cycles $LC_{2}$ and $LC_{3}$ are shown in figures 2 and 3(c). Their amplitudes are easily obtained by measuring the maxima in a time series or the fundamental peak in the temporal Fourier spectrum. In contrast to these, which have maxima of constant amplitudes, the quasi-periodic state existing in the overlap region has maxima of varying heights as shown in figures 2 and 3(b). To extract a single amplitude in their study of the cubic lid-driven cavity, Lopez et al. (Reference Lopez, Welfert, Wu and Yalim2017) used the standard deviation from the mean flow, defined as

(2.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}(U)=\left[\frac{1}{N}\mathop{\sum }_{n=0}^{N}(U(x_{1},y_{1},t_{n})-\overline{U}(x_{1},y_{1}))^{2}\right]^{1/2}, & & \displaystyle\end{eqnarray}$$

where $U(x_{1},y_{1},t_{n})$ the streamwise velocity measured at $(x_{1},y_{1})=(1.2,0.2)$ and at each instant $t_{n}$ , $N$ is the number of measurements in the time series and $\overline{U}(x_{1},y_{1})$ is the temporal mean. We used the edge state technique described in § 2.5 to compute a time series in which the $QP$ is maintained for a long time. In figure 3(d), we show the standard deviation of the time series plotted in figure 3(a). The standard deviation is computed over all times of a sliding window containing fifty peaks. Once the deviation is computed, the window is shifted by ten peaks and we compute the deviation again over fifty peaks. Figure 3(d) shows two regimes of constant $\unicode[STIX]{x1D709}(U)$ corresponding to $QP$ and $LC_{2}$ , justifying the choice of $\unicode[STIX]{x1D709}(U)$ for the bifurcation diagram.

2.7 Hilbert transform

We use the Hilbert transform to obtain spatial characteristics of the flow. The Hilbert transform is a useful means of extracting the local amplitude and phase from a signal. A complex analytic signal $f_{a}(x)$ is constructed from real data $f(x)$ . In contrast to a real signal that has negative and positive frequencies, $f_{a}$ is complex and has only positive frequencies. This signal is obtained by

(2.20) $$\begin{eqnarray}\displaystyle f_{a}(x)=f(x)+\text{i}{\mathcal{H}}(f(x)). & & \displaystyle\end{eqnarray}$$

The imaginary part ${\mathcal{H}}(f(x))$ is its Hilbert transform, defined by phase shifting the positive and negative frequencies of the original real signal by $-\unicode[STIX]{x03C0}/2$ and $\unicode[STIX]{x03C0}/2$ , respectively. More details about the Hilbert transform can be found in Smith (Reference Smith2007). Equation (2.20) is written in polar form

(2.21) $$\begin{eqnarray}\displaystyle f_{a}(x)=A(x)\text{e}^{\text{i}\unicode[STIX]{x1D6F7}(x)}. & & \displaystyle\end{eqnarray}$$

Thus we can extract the envelope (amplitude $A(x)$ ) and the phase $\unicode[STIX]{x1D6F7}(x)$ from the analytic signal at each location, which is the main interest in using the Hilbert transform. We present the results in detail in § 3.3.

3.1 Overview

The flow over a shear-driven cavity at low Reynolds number consists of a free laminar shear layer and one large recirculation within the cavity. As we increase the Reynolds number, the mixing layer is fed by the shear stress and its thickness develops over the cavity. As widely presented in several studies (Rockwell & Naudascher Reference Rockwell and Naudascher1978), it is common to observe self-sustained oscillations in such configurations, in which the flow impacts on a wall. In figure 4(a) we show the bifurcation diagram over the range of Reynolds number $Re\in [4000,5000]$ . We represent the standard deviation from the mean of the streamwise velocity at a point as described in § 2.6. The standard deviation is computed over all times $t_{n}$ of the time series corresponding to $LC_{2}$ , $LC_{3}$ or $QP$ . The line $\unicode[STIX]{x1D709}(U)=0$ represents the solution of the stationary Navier–Stokes equations (base flow). We plot the stable and unstable states with bold and dashed curves, respectively.

The base flow is stable for $Re<Re_{2}$ where $Re_{2}\approx 4126$ is the critical Reynolds number of the first Hopf bifurcation. This threshold is obtained by quadratic interpolation of amplitudes and differs by only 0.34 % from that found by Sipp & Lebedev (Reference Sipp and Lebedev2007). Above this threshold, the base solution exists but is unstable. We observe a second Hopf bifurcation at $Re_{3}\approx 4348$ , also from quadratic interpolation of amplitudes, which agrees with the threshold measured by Meliga (Reference Meliga2017), differing only by 0.005 %. These successive Hopf bifurcations lead after saturation by nonlinear interactions to two limit cycles which we name $LC_{2}$ and $LC_{3}$ because they display two or three maxima of the vertical velocity fluctuations, as will be shown in the next section.

Figure 4(b) shows the schematic phase portraits corresponding to the bifurcation diagram. The stable base flow (i) loses its stability through a primary Hopf bifurcation (ii) at $Re_{2}$ producing the limit cycle $LC_{2}$ . (iii) Another primary Hopf bifurcation at $Re_{3}$ produces the limit cycle $LC_{3}$ . (iv) A secondary subcritical Hopf bifurcation from $LC_{3}$ at $Re_{3}^{\prime }$ produces the quasi-periodic state $QP$ , which moves (v) in phase space towards $LC_{2}$ until it undergoes another secondary subcritical Hopf bifurcation (vi) at $Re_{2}^{\prime }$ which destroys $QP$ and destabilizes $LC_{2}$ . Above $Re_{3}^{\prime }$ , $LC_{3}$ is stable at least until $Re=5000$ . Another Hopf bifurcation and interesting dynamics occur above $Re=5000$ but these will not be discussed in this paper.

Figure 4. (a) Bifurcation diagram of the shear-driven cavity flow for $Re\in [4000,5000]$ . On the $x$ -axis, we show the Reynolds number $Re$ and on the $y$ -axis the standard deviation from the mean of the streamwise velocity at one point. The bold dots on the curves and the thick ticks on the $x$ -axis show the critical Reynolds numbers. The integers show the number of unstable directions (counting a complex conjugate pair as a single direction). We represent stable states by bold curves and unstable ones by dashed curves. The line $\unicode[STIX]{x1D709}(U)=0$ indicates the stationary base flow. The first Hopf bifurcation occurs at $Re_{2}\approx 4126$ and the second at $Re_{3}\approx 4348$ . These two thresholds have been calculated by a quadratic interpolation from the amplitudes. The Hopf bifurcations give rise to limit cycles $LC_{2}$ and $LC_{3}$ . $LC_{2}$ is stable from its threshold until $Re_{2}^{\prime }\approx 4600$ where it loses stability to $LC_{3}$ . $LC_{3}$ is unstable from $Re_{3}$ to $Re_{3}^{\prime }\approx 4410$ and above this Reynolds number it becomes stable. Between $LC_{2}$ and $LC_{3}$ in the overlap region $Re\in [4400,4600]$ there exists a quasi-periodic state $QP$ , which has been computed by using $\unicode[STIX]{x1D6FC}LC_{2}+(1-\unicode[STIX]{x1D6FC})LC_{3}$ as an initial condition for the full nonlinear simulation where $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}(Re)$ . (b) From left to right, the schematic phase portraits corresponding to the bifurcation diagram. The ordinate and abscissa can be considered to be projections onto the eigenmodes leading to $LC_{2}$ and $LC_{3}$ , respectively. The black dots and hollow circles show the stable and unstable states. (i) Stable base flow. (ii) Limit cycle $LC_{2}$ is shown as bifurcating in the vertical direction. (iii) $LC_{3}$ bifurcates in the horizontal direction. (iv,v) The circle moving on the orbit indicates the quasi-periodic state $QP$ . (vi) $QP$ has disappeared, stabilizing $LC_{3}$ .

Figure 5. Stationary solution (base flow). (a)  $\unicode[STIX]{x1D6FA}_{b}$ over the domain for $Re=4500$ . (b) Vorticity profiles $\unicode[STIX]{x1D6FA}_{b}(x=0.5,y)$ for $Re=4200,4500,5000,6000$ . The major change in the base flow with Reynolds number is localized on the shear layer region and at the bottom of the cavity.

This bifurcation diagram is very similar to that seen in the related configuration of a lid-driven cavity (Tiesinga et al. Reference Tiesinga, Wubs and Veldman2002). These authors observed two successive primary Hopf bifurcations leading to branches with different spatial characteristics. Stability is transferred from the first branch to the second via a secondary branch, which is itself created and destroyed via subcritical secondary Hopf bifurcations. The Reynolds numbers they report for the lid-driven cavity are $Re_{2}\approx 8375$ , $Re_{3}\approx 8600$ , $Re_{3}^{\prime }=8800$ , $Re_{2}^{\prime }=9150$ , quite close to twice those we find for the shear-driven cavity. Indeed, this is a classic scenario which occurs in many hydrodynamic configurations; see Kuznetsov (Reference Kuznetsov1998), Marques, Lopez & Shen (Reference Marques, Lopez and Shen2002), Meliga, Gallaire & Chomaz (Reference Meliga, Gallaire and Chomaz2012). Representing the dynamics by two complex amplitudes $r_{2}\text{e}^{\text{i}\unicode[STIX]{x1D719}_{2}}$ , $r_{3}\text{e}^{\text{i}\unicode[STIX]{x1D719}_{3}}$ , the dynamics is governed by the normal form,

(3.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{r}}_{2}=r_{2}(\unicode[STIX]{x1D707}_{2}-a_{22}r_{2}^{2}-a_{23}r_{3}^{2})\quad \dot{\unicode[STIX]{x1D719}}_{2}=\unicode[STIX]{x1D714}_{2} & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{r}}_{3}=r_{3}(\unicode[STIX]{x1D707}_{3}-a_{32}r_{2}^{2}-a_{33}r_{3}^{2})\quad \dot{\unicode[STIX]{x1D719}}_{3}=\unicode[STIX]{x1D714}_{3}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D707}_{j}$

$a_{ij}$

$Re$

$a_{22}$

$a_{33}$

$(r_{2},r_{3})$

Figure 6. Instantaneous visualizations of $LC_{2}$ at $Re=4500$ . Panel (a–c) shows the spanwise vorticity $\unicode[STIX]{x1D6FA}_{2}$ , the horizontal velocity $u_{2}$ , and the vertical velocity $v_{2}$ . Panel (d–f) shows the fluctuations, obtained by subtracting the mean, e.g. $\unicode[STIX]{x1D6FA}_{2}^{\prime }=\unicode[STIX]{x1D6FA}_{2}-\overline{\unicode[STIX]{x1D6FA}}_{2}$ , where $\overline{\unicode[STIX]{x1D6FA}}_{2}$ is the average over one temporal period. The mean fields $\overline{\unicode[STIX]{x1D6FA}}_{2}$ , $\overline{u}_{2}$ , $\overline{v}_{2}$ are almost identical to the instantaneous fields and are therefore not shown.

Figure 7. Instantaneous vertical velocity fluctuations ( $v^{\prime }=v-\overline{v}$ ) for $LC_{2}$ over one period for $Re=4500$ . (a)  $t=0$ (b)  $T/4$ (c)  $T/2$ (d)  $3T/4$ . The structures are advected downstream, as they would be for a travelling wave. In the range of the cavity $x\in [0,1]$ we count two maxima of the vertical velocity fluctuations. The velocity contours are deformed downstream for $x>1.2$ , due to the change to a free-slip boundary condition.

Figure 8. Instantaneous vertical velocity fluctuations ( $v^{\prime }=v-\overline{v}$ ) for $LC_{3}$ over one period for $Re=4500$ . (a)  $t=0$ (b)  $T/4$ (c)  $T/2$ (d)  $3T/4$ . We observe the same dynamics as in $LC_{2}$ shown in figure 7. Over the range $x\in [0,1]$ we count three maxima of the vertical velocity fluctuations.

The base flow is shown in figure 5. Figure 5(a) shows a visualization of its spanwise vorticity $\unicode[STIX]{x1D6FA}_{b}(x,y)$ over our domain for $Re=4500$ . The change in the mixing layer with increasing $Re$ is not qualitatively visible on such a plot, so figure 5(b) presents profiles $\unicode[STIX]{x1D6FA}_{b}(x=0.5,y)$ at the midline for Reynolds numbers ranging from 4200 to 6000. These profiles show the increasing steepness of the shear layer at the top of the cavity.

In addition to base flows, we wish to show representations of the limit cycles $LC_{2}$ and $LC_{3}$ . For the flows in our study, the fluctuations are always dominated by the base or mean flow. In figure 6, we show visualizations of the vorticity $\unicode[STIX]{x1D6FA}_{2}$ , horizontal velocity $u_{2}$ and vertical velocity $v_{2}$ for an instantaneous flow from $LC_{2}$ . Panels (a–c) show the full fields and panels (d–f) show the fluctuations obtained by subtracting the temporal mean of each field. The temporal means are not shown, since they are visually indistinguishable from the full fields. Although oscillations can be seen in $\unicode[STIX]{x1D6FA}_{2}$ in figure 6(a), these are very weak, and even less visible in $u_{2}$ , $v_{2}$ in figure 6(b,c). Therefore to distinguish between $LC_{2}$ and $LC_{3}$ and the base flow, we examine the fluctuations, where the periodic structures along the top and downstream side of the cavity are quite visible. Both $\unicode[STIX]{x1D6FA}_{2}^{\prime }$ and $u_{2}^{\prime }$ in figure 6(d,e) display a two-layer structure with an abrupt vertical change of sign at the top of the cavity. Since we are interested primarily in the horizontal variation and propagation of structures along the top of the cavity, we choose to plot the quantity which most emphasizes this feature, namely the vertical velocity fluctuations $v_{2}^{\prime }$ , as in figure 6(f).

In figures 7 and 8 we show the instantaneous vertical velocity fluctuations for $LC_{2}$ and $LC_{3}$ over one period. $LC_{2}$ and $LC_{3}$ are also depicted in the supplementary movies available at https://doi.org/10.1017/jfm.2019.422. We observe four structures, i.e. two maxima and two minima of the vertical velocity fluctuations, in $LC_{2}$ and six structures, i.e. three maxima and three minima of the vertical velocity fluctuations, in $LC_{3}$ . Similar visualizations can be seen in calculations by Barbagallo et al. (Reference Barbagallo, Sipp and Schmid2009) and experiments by Basley et al. (Reference Basley, Pastur, Lusseyran, Faure and Delprat2011, Reference Basley, Pastur, Delprat and Lusseyran2013).

The behaviour resembles that of travelling waves. The structures, produced by a feedback mechanism, progress steadily to the right but the overall amplitude is not uniform. At the downstream corner, the structures split, as reported by Rockwell & Knisely (Reference Rockwell and Knisely1980): one part follows the fluid downstream, while the other is entrained by the cavity recirculation and returns to feed the flow at the upstream corner, sustaining the vortex generation. The mechanism producing the oscillations is the same for both limit cycles $LC_{2}$ and $LC_{3}$ . The temporal frequency is near 7 for $LC_{2}$ and near 10 for $LC_{3}$ , as in Sipp & Lebedev (Reference Sipp and Lebedev2007), Meliga (Reference Meliga2017). These modes are selected by the cavity length and the mean velocity of the mixing layer as described by Rossiter (Reference Rossiter1964) and as will be discussed in § 4.4.

3.2 Linearization about the base flow

Figure 9. Eigenvalues $\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}$ from linear stability analysis about the base solution with blue and red crosses, and nonlinear frequency from direct numerical simulations (DNS) with black stars. (a,c) Growth rate $\unicode[STIX]{x1D70E}$ of the most unstable eigenmode. In this range of Reynolds number, are two successive Hopf bifurcations as shown in the bifurcation diagram of figure 4. The growth rates for eigenvalues leading to $LC_{2}$ and $LC_{3}$ are shown with blue and red stars respectively. (b,d) Frequency $\unicode[STIX]{x1D714}$ of linear stability analysis about the base state and nonlinear simulation. The nonlinear and linear frequencies agree at onset, but when we increase the Reynolds number the nonlinear frequency deviates from that resulting from linear stability about the base flow.

Figure 10. Spectra of complex-conjugate eigenvalues of the base flow for (a)  $Re=4200$ , (b)  $Re=4500$ , (c)  $Re=5000$ . Blue circles designate the converged eigenvalues and red stars the eigenvalues that did not converge. In these figures we observe the evolution of modes leading to $LC_{2}$ and $LC_{3}$ . The third converged mode, which is stable in this range of Reynolds number, crosses the imaginary axis for $Re\approx 6000$ (not shown).

Figure 11. Vertical velocity of the real parts of the leading unstable eigenmodes about the base flow at $Re=4500$ . (a) $LC_{2}$ and (b) $LC_{3}$ .

In figure 9 we present the results of linear stability analysis about the base flow. We plot the growth rates $\unicode[STIX]{x1D70E}$ in figure 9(a,c) and the frequencies $\unicode[STIX]{x1D714}$ in figure 9(b,d). As previously stated, two successive Hopf bifurcations correspond to two different modes. We plot in (a,b) the eigenvalue leading to $LC_{2}$ and (c,d) that corresponding to $LC_{3}$ . The zero crossing of the growth rate marks the Hopf bifurcation at which the base flow becomes unstable to that eigenmode. As presented in the bifurcation diagram in figure 4, the base flow acquires a first unstable direction at around $Re_{2}\approx 4126$ and a second unstable direction at $Re_{3}\approx 4348$ . Figure 9(b,d) also shows the nonlinear frequencies for the two limit cycles. These agree with the eigenfrequencies at $Re_{2}$ and $Re_{3}$ , as is necessarily the case for a supercritical bifurcation, but as the Reynolds number increases, the frequencies diverge from one another. The frequency extracted at an early stage of a nonlinear simulation initialized by a small perturbation from the base flow will be equal to that given by linear analysis, but the nonlinear interactions will cause it to evolve with time to the nonlinear frequency. The frequencies leading to $LC_{2}$ and $LC_{3}$ cross at $Re=4540$ (differing by only 0.6 % from the value $Re=4567$ reported by Meliga (Reference Meliga2017)). Below this value of $Re$ , a simulation initialized with a small perturbation from the base flow will converge to $LC_{2}$ , while above, it will converge to $LC_{3}$ . In contrast, transition from $LC_{2}$ to $LC_{3}$ will take place when $Re$ is increased above $Re_{2}^{\prime }\approx 4600$ and from $LC_{3}$ to $LC_{2}$ when $Re$ is decreased below $Re_{3}^{\prime }\approx 4410$ .

Figure 10 shows a portion of the eigenvalue spectra for $Re=4200$ , $4500$ and $5000$ . We show in blue circles the eigenvalues that satisfied the convergence tolerance of $10^{-6}$ and in red stars those that did not converge. The first eigenvalue has crossed the $\unicode[STIX]{x1D70E}=0$ axis by $Re=4200$ and the second eigenvalue has crossed by $Re=4500$ . At this $Re$ there are two unstable eigenvalues with almost the same growth rate. At $Re=5000$ these two unstable eigenvalues have further increased and a third stable mode approaches the imaginary axis, becoming unstable at $Re\approx 6000$ (not shown).

Figure 11(a,b) shows the real parts of the eigenvectors leading to $LC_{2}$ and to $LC_{3}$ . We observe two vertical velocity fluctuation maxima on $LC_{2}$ and three on $LC_{3}$ , as was mentioned in the discussion of figures 7 and 8.

3.3 Spatial analysis and Hilbert transform

We have shown in the previous sections that $LC_{2}$ and $LC_{3}$ have different numbers of structures across the cavity. Figure 12(a,b) shows the vorticity fluctuations $\unicode[STIX]{x1D6FA}^{\prime }$ slightly above the cavity like those shown in figure 6(d), at $y=0.05$ and for $x\in [0,2.5]$ for these limit cycles. Curves from light to dark show the vorticity fluctuations at various phases of the temporal period. These two figures show qualitatively the behaviour of a travelling wave, but quantitatively the wavelength and amplitude are not constant. For this reason, at a fixed time, i.e. for each curve shown in figure 12(a,b), we compute an average wavelength $\overline{\unicode[STIX]{x1D706}}$ . Averaging the wavelength over only $x\in [0,1]$ is not possible because this range contains too few wavelengths. Figure 12(c) shows $\overline{\unicode[STIX]{x1D706}}$ as the dashed and solid curves for $LC_{3}$ and $LC_{2}$ . The wavelengths vary little over time and have temporal averages $\langle \overline{\unicode[STIX]{x1D706}_{2}}\rangle =0.56$ and $\langle \overline{\unicode[STIX]{x1D706}_{3}}\rangle =0.39$ .

Figure 12. Vorticity fluctuations above the cavity $\unicode[STIX]{x1D6FA}^{\prime }(x,y=0.05)$ over one period. We observe the behaviour of a travelling wave for (a) $LC_{2}$ and (b) $LC_{3}$ . Because the wavelength evolves in space, we calculate the average wavelength $\overline{\unicode[STIX]{x1D706}}$ as shown in figure (c), in solid and dashed line for $LC_{2}$ and $LC_{3}$ respectively. The mean wavelength is near $0.5$ for $LC_{2}$ and near $0.4$ for $LC_{3}$ .

Figure 13. Hilbert transform of vorticity fluctuations $\unicode[STIX]{x1D6FA}^{\prime }(x,y=0.05)$ . (a,b) vorticity fluctuation, its spline interpolation and the modulus of the Hilbert transform, which is the envelope of the signal. (c) The ratio between the real and imaginary parts of the Hilbert transform is the phase whose evolution is shown for $LC_{3}$ (red) and for $LC_{2}$ (blue). Linear fits to the phases $\unicode[STIX]{x1D6F7}_{i}$ are shown by dashed and dotted curves.

We use the Hilbert transform presented in § 2.7 to analyse in detail the final curve $\unicode[STIX]{x1D6FA}^{\prime }(x,y=0.05)$ shown in figure 12(a,b). We recall that the Hilbert transform produces from a real signal $f(x)$ a complex signal $f_{a}(x)$ which is written in polar form as $A(x)\text{e}^{\text{i}\unicode[STIX]{x1D6F7}(x)}$ . We show in figure 13(a,b) the vorticity fluctuations $\unicode[STIX]{x1D6FA}_{2}^{\prime }(x)$ and $\unicode[STIX]{x1D6FA}_{3}^{\prime }(x)$ for $LC_{2}$ and $LC_{3}$ with black curves. Because the Hilbert transform is very sensitive, we have interpolated the signals by cubic splines. The figure also shows the amplitudes $A(x)$ of the Hilbert transform of the signals. Over the range $x\in [0.5,2.3]$ , there are two influential locations: the downstream corner at $x=1$ , and the changeover of boundary condition at $x=1.75$ from no slip to free slip. The amplitudes $A(x)$ show maxima at the downstream corner. Figure 13(c) shows the phase $\unicode[STIX]{x1D6F7}(x)$ for $LC_{2}$ and $LC_{3}$ . The slope of $\unicode[STIX]{x1D6F7}(t)$ is the wavenumber $k$ . We show with dashed and dotted lines the linear regression calculated over $x\in [0.5,1.5]$ . The wavenumber for $LC_{2}$ obtained in this way is $k=12.115$ and that for $LC_{3}$ is $k=16.045$ , leading to wavelengths $\unicode[STIX]{x1D706}_{2}=0.519$ and $\unicode[STIX]{x1D706}_{3}=0.392$ . These values are fairly close to the values $\langle \overline{\unicode[STIX]{x1D706}}_{2}\rangle =0.56$ and $\langle \overline{\unicode[STIX]{x1D706}}_{3}\rangle =0.39$ obtained by measuring wavelengths, as shown in figure 12. The value $\unicode[STIX]{x1D706}_{2}\approx 0.5$ justifies our designation of $LC_{2}$ as containing two vertical velocity fluctuation maxima, since $L/\unicode[STIX]{x1D706}_{2}\approx 2$ , but the value $\unicode[STIX]{x1D706}_{3}\approx 0.4$ leads to $L/\unicode[STIX]{x1D706}_{3}\approx 2.5$ rather than 3.

3.4 Quasi-periodic state and Floquet analysis

Figure 14. Temporal Fourier spectra of the quasi-periodic state (black bold curves), computed with the edge state technique at (a)  $Re=4420$ , near $Re_{3}^{\prime }$ , (b)  $Re=4500$ and (c)  $Re=4580$ , near $Re_{2}^{\prime }$ . The vertical lines show the frequencies $\unicode[STIX]{x1D714}_{2}$ (blue dotted lines) and $\unicode[STIX]{x1D714}_{3}$ (red dashed lines) of $LC_{2}$ and $LC_{3}$ at the corresponding values of $Re$ . (a) Near $Re_{3}^{\prime }$ where $QP$ is near $LC_{3}$ , the peak $\unicode[STIX]{x1D714}_{QP_{3}}$ of the $QP$ spectrum matches its analogue $\unicode[STIX]{x1D714}_{3}$ (red dashed line) on $LC_{3}$ , while $\unicode[STIX]{x1D714}_{QP_{2}}$ is slightly to the left of $\unicode[STIX]{x1D714}_{2}$ (blue dotted line) of $LC_{2}$ . (c) Near $Re_{2}^{\prime }$ , where $QP$ is near $LC_{2}$ , the peak of $\unicode[STIX]{x1D714}_{QP_{2}}$ of the $QP$ spectrum matches its analogue $\unicode[STIX]{x1D714}_{2}$ (blue dotted line), while $\unicode[STIX]{x1D714}_{QP_{3}}$ is to the right of $\unicode[STIX]{x1D714}_{3}$ (red dashed line) from $LC_{3}$ . (b) Away from $Re_{2}^{\prime }$ and $Re_{3}^{\prime }$ , both frequencies $\unicode[STIX]{x1D714}_{QP_{2}}$ and $\unicode[STIX]{x1D714}_{QP_{3}}$ are shifted slightly from their analogues on $LC_{2}$ and $LC_{3}$ . The vertical bold and dash-dotted lines show the frequencies calculated by Floquet analysis about $LC_{2}$ by $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{2}+\unicode[STIX]{x1D714}_{F}$ (dash-dotted line) and about $LC_{3}$ by $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{3}-\unicode[STIX]{x1D714}_{F}$ (bold line). Near the thresholds at (a)  $Re=4420$ and (c)  $Re=4580$ , the frequency obtained by Floquet analysis about $LC_{3}$ and $LC_{2}$ are very close to $\unicode[STIX]{x1D714}_{QP_{2}}$ and $\unicode[STIX]{x1D714}_{QP_{3}}$ . At $Re=4500$ in figure (b) the frequencies obtained by Floquet analysis about $LC_{2}$ and $LC_{3}$ are also close to their analogues in the quasi-periodic state even though the linear analysis is only valid at the vicinity of the thresholds.

As presented in the previous sections, there is a range of $Re$ over which two limit cycles coexist, separated by a quasi-periodic state $QP$ . We mention here that this state is probably periodic rather than strictly quasi-periodic, because of the well-known nonlinear phenomenon of frequency locking, but its effective period is very long and we will continue to consider it to be quasi-periodic. Figures 2 and 3(b) show the time series corresponding to this state. Figure 14 presents temporal Fourier spectra for three values of $Re$ . The $QP$ state has two fundamental frequencies close to $\unicode[STIX]{x1D714}_{2}$ and $\unicode[STIX]{x1D714}_{3}$ as shown by the dotted and dashed lines in figure 14(a–c). $QP$ can be viewed as a nonlinear superposition of $LC_{2}$ and $LC_{3}$

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{U}_{QP}(x,y,t)=\mathop{\sum }_{n}\mathop{\sum }_{m}c_{n,m}\text{e}^{\text{i}(n\unicode[STIX]{x1D714}_{QP_{2}}+m\unicode[STIX]{x1D714}_{QP_{3}})t}, & & \displaystyle\end{eqnarray}$$

where $n,m\in \mathbb{N}$ , and $\unicode[STIX]{x1D714}_{QP_{2}}$ and $\unicode[STIX]{x1D714}_{QP_{3}}$ are the fundamental frequencies of $QP$ . The blue dotted lines show the nonlinear frequency of $LC_{2}$ and the red dashed lines that of $LC_{3}$ at the corresponding values of $Re$ .

We now interpret the spectra of figure 14 in the context of the bifurcation diagram of figure 4. Because $Re=4420$ is close to $Re_{3}^{\prime }$ , the quasi-periodic state at $Re=4420$ is close to the limit cycle $LC_{3}$ . In agreement with this, the peak at $\unicode[STIX]{x1D714}_{QP_{3}}$ matches almost exactly the nonlinear frequency $\unicode[STIX]{x1D714}_{3}=10.38$ of $LC_{3}$ indicated by the red dashed line in figure 14(a). In contrast $QP$ is not close to $LC_{2}$ at this Reynolds number and so the peak at $\unicode[STIX]{x1D714}_{QP_{2}}$ is to the left of the frequency $\unicode[STIX]{x1D714}_{2}=7.58$ of $LC_{2}$ (blue dotted line). At figure 14(c), corresponding to $Re=4580$ near $Re_{2}^{\prime }$ , the situation is naturally reversed. The peak at $\unicode[STIX]{x1D714}_{QP_{2}}$ matches almost exactly its analogue $\unicode[STIX]{x1D714}_{2}=7.634$ on $LC_{2}$ (blue dotted line) since it is close to $LC_{2}$ , while $\unicode[STIX]{x1D714}_{QP_{3}}$ is slightly to the right of $\unicode[STIX]{x1D714}_{3}=10.445$ of $LC_{3}$ . Away from $LC_{2}$ and $LC_{3}$ at $Re=4500$ , both frequencies $\unicode[STIX]{x1D714}_{2}=7.609$ and $\unicode[STIX]{x1D714}_{3}=10.412$ are slightly shifted from their analogues on the quasi-periodic state.

Figure 15. Floquet multipliers for $LC_{2}$ : (a)  $Re=4500$ , (b)  $Re=4600$ and $LC_{3}$ : (c)  $Re=4420$ , (d)  $Re=4500$ . The small boxes are enlargements of the region surrounding the dominant Floquet multiplier.

We have found from the nonlinear simulations that $LC_{2}$ loses stability towards $LC_{3}$ for $Re>Re_{2}^{\prime }\approx 4600$ and $LC_{3}$ gains stability for $Re>Re_{3}^{\prime }\approx 4420$ . To shed light on the stability of these limit cycles, we carry out a Floquet analysis. In the Floquet framework, we decompose the velocity field as

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{U}(x,y,t)=\boldsymbol{U}_{LC}(x,y,t)+\unicode[STIX]{x1D716}\text{e}^{(\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F})t}\boldsymbol{u}_{F}(x,y,t)+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

with $\boldsymbol{U}_{LC}$ the periodic solution corresponding to the limit cycle about which the Floquet analysis is performed, $\boldsymbol{u}_{F}$ the Floquet mode which is also periodic with period $T_{b}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{b}$ and $\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F}$ the Floquet exponent. We rewrite (3.3) by expressing the Floquet mode as a Fourier series, leading to

(3.4) $$\begin{eqnarray}\displaystyle \boldsymbol{U}(x,y,t)=\boldsymbol{U}_{LC}(x,y,t)+\unicode[STIX]{x1D716}\text{e}^{(\unicode[STIX]{x1D70E}_{F}+\text{i}\unicode[STIX]{x1D714}_{F})t}\mathop{\sum }_{n}\boldsymbol{u}_{F,n}(x,y)\text{e}^{\text{i}n\unicode[STIX]{x1D714}_{b}t}+\text{c.c.} & & \displaystyle\end{eqnarray}$$

For $t=T_{b}$ , we have $\text{e}^{\text{i}n\unicode[STIX]{x1D714}_{b}T_{b}}=\text{e}^{\text{i}n\unicode[STIX]{x1D714}_{b}(2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{b})}=\text{e}^{\text{i}n2\unicode[STIX]{x03C0}}=1$ and so (3.4) becomes

(3.5) $$\begin{eqnarray}\displaystyle \boldsymbol{U}(x,y,T_{b})=\boldsymbol{U}_{LC}(x,y,T_{b})+\unicode[STIX]{x1D716}\text{e}^{\unicode[STIX]{x1D70E}_{F}2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{b}}\text{e}^{\text{i}2\unicode[STIX]{x03C0}\unicode[STIX]{x1D714}_{F}/\unicode[STIX]{x1D714}_{b}}\mathop{\sum }_{n}\boldsymbol{u}_{F,n}(x,y)+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D707}\equiv \text{e}^{\unicode[STIX]{x1D70E}_{F}2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{b}}$ the modulus and $\unicode[STIX]{x1D703}\equiv 2\unicode[STIX]{x03C0}\unicode[STIX]{x1D714}_{F}/\unicode[STIX]{x1D714}_{b}$ the argument of the Floquet multipliers.

Figure 15 shows the Floquet multipliers for both limit cycles in the complex plane. All Floquet multipliers (dots) are inside the unit circle, meaning that the corresponding limit cycles are stable at these Reynolds numbers. Figure 15(a,b) shows the results for $LC_{2}$ at $Re=4500$ and $4600$ respectively. In figure 15(a) at $Re=4500$ , the dominant Floquet multiplier modulus is $|\unicode[STIX]{x1D707}|=0.981$ . On figure 15(b), by $Re=4600\lesssim Re_{2}^{\prime }$ , this multiplier has moved closer to the unit circle, with $|\unicode[STIX]{x1D707}|=0.999$ . Figure 15(c,d) shows that the moduli of the dominant Floquet multipliers for $LC_{3}$ at $Re=4420\gtrsim Re_{3}^{\prime }$ and at $Re=4500$ are $|\unicode[STIX]{x1D707}|=0.995$ and $|\unicode[STIX]{x1D707}|=0.967$ . Thus the results shown by the Floquet analysis confirm the nonlinear observations. The results of this Floquet analysis resemble those found for the lid-driven cavity by Tiesinga et al. (Reference Tiesinga, Wubs and Veldman2002).

We now turn to the argument of the Floquet multipliers. If the Floquet exponent is real $(\unicode[STIX]{x1D714}_{F}=0)$ then the Floquet multiplier is one and the bifurcating state has the same frequency as the base limit cycle. If $\unicode[STIX]{x1D714}_{F}/\unicode[STIX]{x1D714}_{b}=1/2$ then the Floquet multiplier is $-1$ , which corresponds to a subharmonic mode. In our problem the dominant Floquet multiplier is complex, and so the bifurcation is a secondary Hopf bifurcation and the solution near the threshold of $QP$ is described by (3.4). The spectrum of $QP$ near the threshold contains $\unicode[STIX]{x1D714}_{b}$ and its multiples as well as the frequencies introduced by the secondary Hopf bifurcation, namely $\pm \unicode[STIX]{x1D714}_{F}\pm n\unicode[STIX]{x1D714}_{b}$ , with a dominant contribution from $n=\pm 1$ . Indeed, near $Re_{2}^{\prime }$ , the spectrum of $QP$ contains frequencies $\unicode[STIX]{x1D714}_{2}$ and $\unicode[STIX]{x1D714}_{2}+\unicode[STIX]{x1D714}_{F}=\unicode[STIX]{x1D714}_{2}(1+\unicode[STIX]{x1D703}/2\unicode[STIX]{x03C0})=10.45$ while near $Re_{3}^{\prime }$ , $QP$ contains frequencies $\unicode[STIX]{x1D714}_{3}$ and $\unicode[STIX]{x1D714}_{3}-\unicode[STIX]{x1D714}_{F}=\unicode[STIX]{x1D714}_{3}(1-\unicode[STIX]{x1D703}/2\unicode[STIX]{x03C0})=7.37$ .

These calculations are confirmed by figure 14(a,c). For (a) $Re=4420\gtrsim Re_{3}^{\prime }$ , the Floquet analysis about $LC_{3}$ yields a frequency $\unicode[STIX]{x1D714}_{3}-\unicode[STIX]{x1D714}_{F}$ (solid vertical black line) comparable to the peak at $\unicode[STIX]{x1D714}_{QP_{2}}$ For (c) $Re=4580\lesssim Re_{2}^{\prime }$ , the Floquet analysis about $LC_{2}$ yields the frequency $\unicode[STIX]{x1D714}_{2}+\unicode[STIX]{x1D714}_{F}$ (dashed vertical black line) which is very close to  $\unicode[STIX]{x1D714}_{QP_{3}}$ .

4.1 Linearization about the mean flow

Figure 16. Distortion of the mean flow ( $\unicode[STIX]{x1D6FA}^{\ast }=\overline{\unicode[STIX]{x1D6FA}}-\unicode[STIX]{x1D6FA}_{b}$ ). For $LC_{2}$ at (a)  $Re=4200$ and (b)  $Re=4500$ . For $LC_{3}$ at (c)  $Re=4500$ and (d)  $Re=5000$ . The distortion measures the deviation of the mean flow from the base flow due to nonlinear interaction and increases with Reynolds number.

Linear stability analysis – i.e. linearizing about the base flow and solving the resulting eigenproblem – is a classic tool in hydrodynamics. Bifurcations which create new branches are determined unambiguously by linear stability analysis and, if the bifurcation is supercritical, the spatial and temporal behaviours of the new states near threshold are similar to those of the eigenvector and eigenvalue responsible for the instability. Further from the threshold, these properties evolve and may well differ substantially from those of the bifurcating eigenvector and eigenvalue. In some cases, it has been shown that linearization about the mean flow of a limit cycle can yield more accurate approximations of the nonlinear states. We have carried out a linear analysis about the temporal mean for both limit cycles $LC_{2}$ and $LC_{3}$ and compared the resulting frequencies with those obtained from linearization about the base flow and with the nonlinear frequencies of these cycles. This procedure has been carried out for $LC_{3}$ by Meliga (Reference Meliga2017); here we carry out the same procedure for $LC_{2}$ and compare the two regimes.

The mean flow greatly resembles the base flow shown in figure 5, and therefore we plot instead their difference, the distortion $\overline{\unicode[STIX]{x1D6FA}}^{\ast }\equiv \overline{\unicode[STIX]{x1D6FA}}-\unicode[STIX]{x1D6FA}_{b}$ , which is shown for $LC_{2}$ at $Re=4200$ and 4500 in figure 16(a,b) and for $LC_{3}$ at $Re=4500$ and 5000 in figure 16(c,d). Nor do we show the eigenmodes obtained by linearizing about the mean flow, since they resemble those obtained from linearizing about the base flow (figure 11) as well as the nonlinear vertical velocity fluctuations (figures 7 and 8).

Figure 17. Eigenvalues $\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}$ from linear stability analysis about the base and mean flows of $LC_{2}$ and $LC_{3}$ . The eigenvalues about the base flow are plotted with crosses. Linear stability about the mean flow is shown with circles. The nonlinear frequency from DNS is plotted with black stars. (a,c) The growth rate $\unicode[STIX]{x1D70E}$ for $LC_{2}$ for the mean flow is nearly zero but that of $LC_{3}$ is smaller than that about the base flow, but not enough to be considered to be neutrally stable. (b,d) The imaginary parts of the eigenvalues are almost exactly equal to the nonlinear frequencies, especially for $LC_{2}$ . For $LC_{3}$ , although the imaginary part of the eigenvalue and the nonlinear frequency are necessarily equal at onset, the two diverge slightly as $Re$ increases.

Figure 17 and table 1 compare the eigenvalues resulting from linearization about the base flow and the mean flow and those of the nonlinear simulation. Figure 17 plots the frequencies and growth rates over the Reynolds number range $[4000,5000]$ that we have studied, while table 1 shows numerical data extracted from these figures for three representative Reynolds numbers, 4200, 4500 and 5000. We note that unlike for the cylinder wake (Barkley Reference Barkley2006), the frequencies obtained from the usual linear stability analysis are already not very far from those of the nonlinear limit cycles. Table 1 shows a deviation of less than 0.6 % for the frequencies in $LC_{2}$ and of less than 1.4 % for $LC_{3}$ over this Reynolds number range. In contrast, for the cylinder wake (Barkley Reference Barkley2006), the difference between the nonlinear frequencies and those obtained by linear stability analysis reaches 15 % by $Re=60$ , a Reynolds number comparable to the distance above criticality studied here and reaches 100 % by $Re=180$ , a frequent upper limit of such studies. This difference will be discussed further in §§ 4.4 and 5.

Figure 17(b) for $LC_{2}$ (circles) shows that the frequency obtained by linearization about the mean flow is nonetheless much closer to the nonlinear temporal frequency (stars) than that given by linear stability about the base flow (crosses). Quantitatively, table 1 shows the relative difference at $Re=4500$ between the nonlinear frequency and the frequency obtained from the base flow to be 0.7 %; this difference is reduced to 0.04 % when the linearization is performed about the mean flow. Moreover, the growth rate obtained about the mean flow (circles) in figure 17(a) is nearly zero, as found by Barkley (Reference Barkley2006) for the cylinder wake. Table 1 shows a growth rate at $Re=4500$ of 0.073 for linearization about the base flow; this is reduced by a factor of 5 to 0.017 for the linearization about the mean.

For $LC_{3}$ , the frequency obtained by linearizing about the mean presented in figure 17(d) (circles) agrees well with the nonlinear frequency (stars). The curves begin to diverge slightly for $Re\geqslant 4600$ , and although the agreement is not as good as it is for $LC_{2}$ , the frequencies are still very close. Quantitatively, table 1 shows the relative difference at $Re=5000$ between the nonlinear frequency and the frequency obtained from the base flow to be 1.4 %; this difference is reduced to 0.02 % when the linearization is performed about the mean flow. Table 1 shows a growth rate at $Re=5000$ of 0.247 for linearization about the base flow; this is again reduced by a factor of 5 to 0.053 for the linearization about the mean.

Table 1. Linear and nonlinear frequencies for cavity modes. The first row shows frequencies from linearization about the base and from linearization about the mean and nonlinear simulations, where available. Second row for $\unicode[STIX]{x1D714}_{2}$ and $\unicode[STIX]{x1D714}_{3}$ shows deviation from frequencies from nonlinear simulations. When the RZIF property is satisfied, linearization about the nonlinear mean yields frequencies close to nonlinear frequencies. The third row for $\unicode[STIX]{x1D714}_{2}$ and $\unicode[STIX]{x1D714}_{3}$ shows linear growth rates for cavity modes. Linearization about the base flow yields growth rates which cross zero transversely as the bifurcation threshold is crossed. When the RZIF property is satisfied, linearization about the nonlinear mean yields growth rates near zero, i.e. the mean flow is nearly marginally stable. Last columns show the frequency increment between consecutive eigenvalues, which is constant to two digits, regardless of the Reynolds number or which type of frequency is used.

4.2 RZIF and SCM

We now present an argument for the validity of linearization about the mean flow. The derivation is closely related to the idea called harmonic balance in the aerodynamic literature (Hall, Thomas & Clark Reference Hall, Thomas and Clark2002; McMullen, Jameson & Alonso Reference McMullen, Jameson and Alonso2006; McMullen & Jameson Reference McMullen and Jameson2006). Turton et al. (Reference Turton, Tuckerman and Barkley2015) argued that the RZIF property holds exactly if the time dependence is monochromatic, meaning that higher harmonics are negligible compared to the fundamental frequency. Consider the evolution equation

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{U}={\mathcal{L}}\boldsymbol{U}+{\mathcal{N}}(\boldsymbol{U},\boldsymbol{U}), & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}$ is linear and ${\mathcal{N}}(\cdot ,\cdot )$ is a quadratic nonlinearity. Let

(4.2) $$\begin{eqnarray}\displaystyle \boldsymbol{U}=\overline{\boldsymbol{U}}+\mathop{\sum }_{n\neq 0}\boldsymbol{u}_{n}\text{e}^{\text{i}n\unicode[STIX]{x1D714}t}, & & \displaystyle\end{eqnarray}$$

(with $\boldsymbol{u}_{-n}=\boldsymbol{u}_{n}^{\ast }$ ) be the temporal Fourier decomposition of a periodic solution to (4.1) with mean $\overline{\boldsymbol{U}}$ and frequency $\unicode[STIX]{x1D714}$ . In (4.2) and throughout this subsection, subscripts refer to temporal harmonics of a single limit cycle, rather than serving to distinguish between $LC_{2}$ and $LC_{3}$ .

The $n=0$ (mean) component of (4.1) is

(4.3) $$\begin{eqnarray}\displaystyle 0={\mathcal{L}}\overline{\boldsymbol{U}}+{\mathcal{N}}(\boldsymbol{u}_{1},\boldsymbol{u}_{-1})+{\mathcal{N}}(\boldsymbol{u}_{-1},\boldsymbol{u}_{1})+{\mathcal{N}}(\boldsymbol{u}_{2},\boldsymbol{u}_{-2})+{\mathcal{N}}(\boldsymbol{u}_{-2},\boldsymbol{u}_{2})+\cdots & & \displaystyle\end{eqnarray}$$

while the $n=1$ component is

(4.4) $$\begin{eqnarray}\displaystyle \text{i}\unicode[STIX]{x1D714}\boldsymbol{u}_{1} & = & \displaystyle \underbrace{{\mathcal{L}}\boldsymbol{u}_{1}+{\mathcal{N}}(\overline{\boldsymbol{U}},\boldsymbol{u}_{1})+{\mathcal{N}}(\boldsymbol{u}_{1},\overline{\boldsymbol{U}})}_{{\mathcal{L}}_{\overline{U}}\boldsymbol{u}_{1}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{{\mathcal{N}}(\boldsymbol{u}_{2},\boldsymbol{u}_{-1})+{\mathcal{N}}(\boldsymbol{u}_{-1},\boldsymbol{u}_{2})+{\mathcal{N}}(\boldsymbol{u}_{3},\boldsymbol{u}_{-2})+{\mathcal{N}}(\boldsymbol{u}_{-2},\boldsymbol{u}_{3})+\cdots }_{{\mathcal{N}}_{1}}.\end{eqnarray}$$

If, as is often the case, $\Vert u_{n}\Vert \sim \unicode[STIX]{x1D716}^{|n|}$ , then ${\mathcal{N}}_{1}=O(\unicode[STIX]{x1D716}^{3})$ may be neglected and RZIF is satisfied: the linear operator ${\mathcal{L}}_{\overline{U}}$ in (4.4) has the pure imaginary eigenvalue $\text{i}\unicode[STIX]{x1D714}$ , corresponding to the frequency of the periodic solution. Hence the RZIF property is satisfied for near-monochromatic oscillations in a system with quadratic nonlinearity.

We mention here that RZIF is not predictive, since it requires a full nonlinear direct numerical simulation to be carried out in order to compute the temporal mean $\overline{U}$ . An approach which is actually predictive, i.e. which does not require a DNS, has been proposed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014, Reference Mantič-Lugo, Arratia and Gallaire2015), who called it the self-consistent model (SCM). The SCM truncates the mean flow equation (4.3) as well as the $n=1$ equation (4.4), leading to the closed system

(4.5) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}0={\mathcal{L}}\overline{\boldsymbol{U}}+{\mathcal{N}}(\boldsymbol{u}_{1},\boldsymbol{u}_{-1})+{\mathcal{N}}(\boldsymbol{u}_{-1},\boldsymbol{u}_{1}),\\ \text{i}\unicode[STIX]{x1D714}\boldsymbol{u}_{1}={\mathcal{L}}\boldsymbol{u}_{1}+{\mathcal{N}}(\overline{\boldsymbol{U}},\boldsymbol{u}_{1})+{\mathcal{N}}(\boldsymbol{u}_{1},\overline{\boldsymbol{U}}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

This system is then solved for $\overline{\boldsymbol{U}}$ and $\boldsymbol{u}_{1}$ by various iterative methods; see Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014, Reference Mantič-Lugo, Arratia and Gallaire2015). The next higher truncation, i.e. retaining $\overline{\boldsymbol{U}}$ , $\boldsymbol{u}_{1}$ , and $\boldsymbol{u}_{2}$ , has been studied by Meliga (Reference Meliga2017). It may happen, however, that RZIF is satisfied, while SCM is not, i.e. that while higher-order modes may be neglected in (4.4), they are essential to forming the correct mean flow and cannot be neglected in (4.3); see Bengana (Reference Bengana2019).

For thermosolutal convection, Turton et al. (Reference Turton, Tuckerman and Barkley2015) showed that for travelling waves, the RZIF property is satisfied and the spectrum is highly peaked, while for standing waves the spectrum is broad and the RZIF property is not satisfied. We now wish to see if the temporal spectra of $LC_{2}$ and $LC_{3}$ also explain the fact that the RZIF property is better satisfied for $LC_{2}$ than for $LC_{3}$ . We show in figure 18 the temporal spectra of streamwise velocity normalized by the fundamental frequency for $LC_{2}$ and $LC_{3}$ , for various values of Reynolds number. In figure 19 we plot the ratio of the second harmonic to the fundamental frequency. The ratio $\Vert \hat{u} _{2}\Vert /\Vert \hat{u} _{1}\Vert$ is consistently less than 0.05 for both flows over the range of our investigation, while $(\Vert \hat{u} _{2}\Vert +\Vert \hat{u} _{3}\Vert +\Vert \hat{u} _{4}\Vert )/\Vert \hat{u} _{1}\Vert$  remains below 0.07, consistent with the fact that RZIF is satisfied. We observed that RZIF is closer to being valid for $LC_{2}$ than for $LC_{3}$ , but the ratios in figure 19 follow the opposite tendency. Thus, the explanation proposed by Turton et al. (Reference Turton, Tuckerman and Barkley2015) in terms of the temporal Fourier amplitudes does not explain this difference.

Figure 18. Temporal spectra of streamwise velocity normalized by the fundamental frequency for $LC_{2}$ and $LC_{3}$ . (a) Spectra of $LC_{2}$ at $Re=4200$ and $Re=4500$ . (b) At $Re=4500$ for $LC_{2}$ and $LC_{3}$ . (c) At $Re=4500$ and $Re=5000$ for $LC_{3}$ . (d) Spectra of $LC_{2}$ and $LC_{3}$ at $Re=4200$ and $Re=5000$ , respectively.

4.3 Cross-eigenvalues

Like the base flow, the mean flow has a full spectrum of eigenvalues and eigenvectors. Thus, the mean flow of $LC_{2}$ has not only eigenvectors with two vertical velocity fluctuation maxima with corresponding eigenvalues shown in figure 17(a,b), like those which lead to $LC_{2}$ , but it also has also eigenvectors with three maxima like those which lead to $LC_{3}$ and their corresponding eigenvalues. Similarly, the mean flow of $LC_{3}$ has eigenvectors containing two maxima. We refer to these as cross-eigenvalues. Figure 20 shows the cross-eigenvalues corresponding to mode two, obtained by linearization about the mean flow of $LC_{3}$ (circles, figure 20(a,b)) and those of mode three about the mean flow of $LC_{2}$ (circles, figure 20(c,d)). The eigenvalues obtained from the base and from the mean necessarily agree at $Re_{2}$ for $LC_{2}$ and at $Re_{3}$ for $LC_{3}$ , since when the limit cycles are created, the base and mean flows are equal.

Figure 19. Amplitude of second harmonic (stars) and the sum of the amplitudes of the three lowest harmonics (crosses) normalized by the amplitude of the fundamental frequency as a function of relative Reynolds number for $LC_{2}$ and $LC_{3}$ . As in figure 18, the fundamental frequency dominates the second harmonic. These ratios are always below $10^{-1}$ and are slightly higher for $LC_{2}$ than that for $LC_{3}$ .

Figure 20. Eigenvalues $\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}$ from linearization about the base and mean flows. Here the spatial form of the eigenvector does not correspond to that of the mean flow. The eigenvalues about the base flow are plotted with crosses, blue for $LC_{2}$ and red for $LC_{3}$ . The eigenvalues of mode two about the mean flow of $LC_{3}$ are shown with red circles and that of mode three about the mean flow of $LC_{2}$ with blue circles. The nonlinear frequencies from DNS are plotted with black stars for $LC_{2}$ and $LC_{3}$ . (a) The growth rate of mode two about the mean of $LC_{3}$ (red circles) decreases from the threshold $Re_{3}\approx 4348$ of $LC_{3}$ and becomes negative at $Re\approx 4681$ . This may correspond qualitatively to the fact that $LC_{3}$ is unstable to mode 2 perturbations when it is created at $Re_{3}$ and becomes stable at $Re_{3}^{\prime }\approx 4410$ . (c) The growth rate of mode three about the mean of $LC_{2}$ (blue circles) increases from the threshold $Re_{2}\approx 4126$ of $LC_{2}$ and becomes positive at $Re\approx 4418$ . This may correspond qualitatively to the fact that $LC_{2}$ is stable when created at $Re_{2}$ and becomes unstable at $Re_{2}^{\prime }\approx 4600$ .

Focusing on figure 20(a), we recall that $LC_{3}$ is created at $Re_{3}\approx 4348$ and is unstable to eigenmodes of type 2 until $Re_{3}^{\prime }\approx 4410$ , i.e. there is a Floquet mode with positive Floquet exponent for $Re\in [4348,4410]$ . This is qualitatively the behaviour that is seen in figure 20(a), although here the cross-eigenvalue $\unicode[STIX]{x1D70E}$ is positive over a higher range, for $Re\in [4348,4681]$ . Focusing now on figure 20(c), we recall that $LC_{2}$ is created at $Re_{2}\approx 4126$ and becomes unstable to eigenmodes of type 3 at $Re_{2}^{\prime }\approx 4600$ . This is again qualitatively close to the behaviour is seen in figure 20(c), except that here the cross-eigenvalue $\unicode[STIX]{x1D70E}$ becomes positive at the lower value of $Re\approx 4418$ .

These results indicate that for a limit cycle, linearization about its mean flow may be able to convey information about the growth rate, frequency, and spatial characteristics of its stability to secondary bifurcations. Referring back to figure 4(b), these phase portraits represent the limit cycles as fixed points, whose existence and stability are governed by the dynamics of the two-dimensional normal form (3.1a ), (3.1b ). Stability or instability of $LC_{2}$ to $LC_{3}$ or vice versa is represented by the horizonal arrows emanating to or from $LC_{2}$ , located on the vertical axis, and the vertical arrows emanating to or from $LC_{3}$ , located on the horizontal axis. The temporal mean would then be interpreted as a fixed point approximation to a limit cycle, and the cross-eigenvalues would then correspond to their relative stability. Although this is plausible, we emphasize that there exists at present no mathematical framework for this interpretation. Moreover, although the trends seen in figure 20 support this interpretation, the agreement is only qualitative.

4.4 Rossiter formula

Figure 21. Spectra of the base flow for (a)  $Re=4200$ , (b)  $Re=4500$ , (c)  $Re=5000$ . Blue circles designate the converged eigenvalues and red stars the eigenvalues that did not converge. Horizontal lines emphasize the fact that the frequencies are equally spaced.

We return to table 1, the last column of which shows that the frequencies of successive modes increase by a constant interval. (For a given mode, the various versions of its frequency differ by at most 1–2 %.) We emphasize this again by reproducing the eigenspectra in figure 21, adding horizontal lines which emphasize visually the constant difference between the frequencies. Spectra similar to those in figure 21 are also seen for the lid-driven cavity in Tiesinga et al. (Reference Tiesinga, Wubs and Veldman2002).

In flows over shear-driven cavities, Rossiter (Reference Rossiter1964) observed that the temporal frequencies for self-sustained oscillations were quantized and proposed the following empirical formula:

(4.6) $$\begin{eqnarray}\displaystyle f_{n}=\frac{U_{\infty }}{L}\frac{n-\unicode[STIX]{x1D6FE}}{M+1/\unicode[STIX]{x1D705}}\quad \Longrightarrow \quad \frac{U_{\infty }}{L}\unicode[STIX]{x1D705}\hspace{2.22198pt}(n-\unicode[STIX]{x1D6FE})\quad \text{for }M=0, & & \displaystyle\end{eqnarray}$$

where $U_{\infty }$ and $L$ are the free-stream speed away from the cavity and the length of the cavity, and $M$ is the Mach number, here set to zero. The phenomenological constant $\unicode[STIX]{x1D6FE}$ is a phase lag, while the increment $\unicode[STIX]{x1D705}$ will be discussed below. The essence of (4.6) is not only that the temporal frequencies $f_{n}$ observed are quantized (which is to be expected in a finite cavity) but that they are separated by a fixed increment $\unicode[STIX]{x0394}f$ . Heuristically, if the limit cycle consists of $n$ structures advected horizontally at velocity $U_{adv}$ , then the average structure occupies a length $L/n$ and strikes the cavity corner with frequency $U_{adv}/(L/n)$ . Since the frequencies in table 1 are non-dimensionalized by $U_{\infty }/L$ , we have

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}f=\frac{U_{\infty }}{L}\unicode[STIX]{x1D705}=\frac{U_{adv}}{L}=0.45\quad \Longrightarrow \quad \unicode[STIX]{x1D705}=\frac{U_{adv}}{U_{\infty }}=0.45. & & \displaystyle\end{eqnarray}$$

Thus, the frequency is seen to be determined primarily by the geometry and the free-stream velocity.

At the same time as Rossiter, Kulikovskii (Reference Kulikovskii1966) proposed an approach, involving amplification of the perturbations along the pathway around the cavity, for determining the set of discrete frequencies that appear in a finite domain. Tuerke et al. (Reference Tuerke, Sciamarella, Pastur, Lusseyran and Artana2015) applied this approach to the open cavity in the incompressible regime and showed very good agreement with the experimental results of Basley et al. (Reference Basley, Pastur, Lusseyran, Faure and Delprat2011).