4.1 The saddle-node bifurcation

Following the B-type steady states branch backwards in Reynolds number, the solution would seem to vanish into thin air at $Re_{\mathit{SN}}\simeq 4908$ . A thorough analysis of the transient dynamics reveals that the final decay onto the steady state, while it exists, is exponential and gets slower as $Re_{\mathit{SN}}$ is approached. This is suggestive of a real negative eigenvalue gradually advancing towards a zero crossing. Moreover, the last few solutions before the eventual zero crossing are well fitted by a square-root law following ${\sim}\sqrt{Re-Re_{\mathit{SN }}}$ . Both facts put together are indicative of a saddle-node bifurcation, whence node (state B) and saddle (henceforth referred to as state C) steady-state solution branches are issued. We choose not to provide a detailed description of the saddle-node bifurcation here, as it adds little value to the discussion, but deem it pivotal instead to produce evidence of the existence of the saddle branch.

As a saddle solution that pairs up with a stable nodal solution, state C possesses a one-dimensional unstable manifold that drives trajectories towards nodal state B on one side and pushes dynamics away from itself on the other. It is therefore not reachable by time evolution but its stable manifold bounds, if only locally, the basin of attraction of state B, and separates it from the basin of attraction of every other existing attractor, in this case state A. We have employed the edge-tracking technique (Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006) to explore the dynamics on the basin boundary, which in cases like the one at hand is known to eventually converge onto the saddle solution. Thus, starting from the A-type and B-type solutions at $Re=6000$ , we have followed a refinement process by picking intermediate initial conditions along a straight line in phase space connecting both states $(\boldsymbol{u},p)_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6FC}(\boldsymbol{u},p)_{A}+(1-\unicode[STIX]{x1D6FC})(\boldsymbol{u},p)_{B}$ . The refining bisection process on parameter $\unicode[STIX]{x1D6FC}$ leads to the determination of neighbouring initial conditions that evolve together on the basin boundary (critical threshold, also known as edge) for very long times before eventually parting, one towards state A, the other towards state B. The refinement can be taken down to machine precision, at which time a new refinement sequence can be started between the last bounding trajectories at a later time. This alternate succession of refinement procedures and adequate time leaps allows for tracking edge trajectories indefinitely.

Figure 13(a) shows the result of three consecutive refinement processes with two intervening time advances (indicated by $r_{1}$ and $r_{2}$ ). The black lines correspond to bounding trajectories obtained from each one of the three bisection refinement processes undertaken, the solid line always finally departing towards state A while the dashed line eventually leaving for state B. The locations of states A and B are indicated with grey dashed lines to guide the eye. The trajectory on the basin boundary clearly converges onto a steady state after some initial rather wild transients and a subsequent gradual oscillatory damping (see the inset) that might be indicative of a possibly oncoming Hopf bifurcation at larger $Re$ . In any case, after three refinements and at around $t\gtrsim 150$ the steady state can be considered as sufficiently converged. This steady state, which governs the dynamics on the boundary separating the basins of attraction of states A and B, is shown in figure 13(b) and indicated with an empty circle in figure 13(a). Its relative phase-space position with respect to states A and B is marked with an isolated grey cross labelled C in figure 4, in representation of the full branch of unsteady solutions that would be extremely costly to track. State C is similar but not identical to state B. Their resemblance results from their being connected at the saddle-node bifurcation not so far away at $Re_{SN}$ , but their nonlinear evolution at larger values of the parameter can be largely independent.

Figure 13. Steady state C, resulting from edge tracking at $Re=6000$ between bounding solutions A and B. (a) Time evolution of $u_{c}$ for several initial conditions tightly confining the dynamics to the basin boundary. The A- and B-labelled horizontal grey dashed lines indicate states A and B, respectively. Solid (dashed) lines indicate bounding trajectories eventually escaping towards state A (B). The vertical grey dotted lines indicate two successive edge-tracking refinements (labels $r_{1}$ and $r_{2}$ ). The inset shows an extreme magnification of the closest approach to steady-state C resulting from refinement $r_{2}$ . (b) Steady-state C, taken at the point along the basin boundary indicated with a circle in panel (a).

4.2 The Hopf bifurcation

A mere inspection of the transients that lead to steady-state B of § 3.1 (sufficiently far from the saddle-node bifurcation) reveals that final asymptotic convergence follows a spiralling trend, as evidenced by the $(u_{c}^{\prime },v_{c}^{\prime })$ phase map projection in the inset of figure 14(a). Here the primes denote deviation from steady-state values. This is an indication that the leading eigenmode is a complex-conjugate pair, and that the steady state reduces to a stable focus when the dynamics is projected on the manifold spanned by the leading eigenpair, all other eigendirections contracting faster.

Figure 14. Asymptotic evolution of the perturbation field around the stable steady state. (a) Sequences of ${u_{c}^{\prime }}^{P}$ , normalised by ${u_{c}^{\prime }}^{P}(0)$ , corresponding to the Poincaré map defined by (4.7) for $Re\in [7500,8035]$ . Least-squares fits to the exponential tails are indicated by the grey straight lines. The inset shows, as an example, the phase map that generates the sequence at $Re=7500$ . (b) Vorticity ( $\unicode[STIX]{x1D714}^{\prime }$ , top panel) and streamfunction ( $\unicode[STIX]{x1D713}^{\prime }$ , bottom panel) on a Poincaré section at $Re=8000$ .

The instantaneous flow field can be expressed in terms of the stable steady state $(\bar{\boldsymbol{u}},\bar{p})$ plus a perturbation field $(\boldsymbol{u}^{\prime },p^{\prime })$ :

(4.1a,b ) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x};t)=\bar{\boldsymbol{u}}(\boldsymbol{x})+\boldsymbol{u}^{\prime }(\boldsymbol{x};t),\quad p(\boldsymbol{x};t)=\bar{p}(\boldsymbol{x})+p^{\prime }(\boldsymbol{x};t).\end{eqnarray}$$

The equations of motion for the perturbation are obtained by substituting (4.1) in (2.1):

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

with homogeneous Dirichlet boundary conditions for $\boldsymbol{u}^{\prime }$ .

For sufficiently small perturbation fields of order $O(\unicode[STIX]{x1D700})$ with $\unicode[STIX]{x1D700}\ll 1$ , i.e. for flow fields at a vanishing distance from the steady state, the perturbation advection term $(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }\sim O(\unicode[STIX]{x1D700}^{2})$ becomes negligible relative to all other terms, which remain of order $O(\unicode[STIX]{x1D700})$ , and the equations linearise:

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

The search for solutions of the form $(\boldsymbol{u}^{\prime },p^{\prime })(t)=(\tilde{\boldsymbol{u}},\tilde{p})\text{e}^{\unicode[STIX]{x1D706}t}$ results in a generalised constrained eigenproblem

(4.4) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D706} & 0\\ 0 & 0\end{array}\right]\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{u}}\\ \tilde{p}\end{array}\right)=\left[\begin{array}{@{}cc@{}}-(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\diamond -(\diamond \boldsymbol{\cdot }\unicode[STIX]{x1D735})\bar{\boldsymbol{u}}+{\displaystyle \frac{1}{Re}}\unicode[STIX]{x1D6FB}^{2}\diamond & -\unicode[STIX]{x1D735}\diamond \\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\diamond & 0\end{array}\right]\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{u}}\\ \tilde{p}\end{array}\right),\end{eqnarray}$$

with eigenmodes $(\tilde{\boldsymbol{u}}_{i},\tilde{p}_{i})$ and associated eigenvalues $\unicode[STIX]{x1D706}_{i}$ . The dynamics of the linearised problem exactly decouples into the superposition of the evolution of individual eigenmodes:

(4.5) $$\begin{eqnarray}(\boldsymbol{u}^{\prime },p^{\prime })(t)=\mathop{\sum }_{i}a_{i}(\tilde{\boldsymbol{u}}_{i},\tilde{p}_{i})\text{e}^{\unicode[STIX]{x1D706}_{i}t}\xrightarrow[{t\rightarrow \infty }]{}a_{1}(\tilde{\boldsymbol{u}}_{1},\tilde{p}_{1})\text{e}^{\unicode[STIX]{x1D706}_{1}t}+\text{c.c.},\end{eqnarray}$$

where $a_{i}$ is the projection of the initial perturbation field $(\boldsymbol{u}^{\prime },p^{\prime })(0)$ onto mode $i$ , and eigenmodes have been sorted by decreasing value of the real part of their associated eigenvalue. For a stable steady state, all eigenvalues have negative real part ( $\text{Re}(\unicode[STIX]{x1D706}_{i})<\text{Re}(\unicode[STIX]{x1D706}_{1})<0$ ), which results in exponential decay along all eigendirections. The asymptotic behaviour of the perturbation for long times tends to align with the leading eigenmode $(\tilde{\boldsymbol{u}}_{1},\tilde{p}_{1})$ (the one associated to the eigenvalue with larger real part, $\unicode[STIX]{x1D706}_{1}$ ), as the dynamics contracts faster in all other eigendirections $(\tilde{\boldsymbol{u}}_{i},\tilde{p}_{i})$ (of smaller real part eigenvalues, $\unicode[STIX]{x1D706}_{i}$ , $i>1$ ).

The fully nonlinear evolution towards a stable steady state will experience nonlinear transients that hinder a modal analysis. Once the initial transients have been overcome, however, the dynamics will asymptotically approach the linear regime and the perturbation will align with the leading eigenmode as in the linear case.

For a focus the leading eigenmode is a complex-conjugate pair $(\tilde{\boldsymbol{u}}_{1},\tilde{p}_{1})=(\tilde{\boldsymbol{u}}_{1}^{r},\tilde{p}_{1}^{r})\pm \text{i}(\tilde{\boldsymbol{u}}_{1}^{i},\tilde{p}_{1}^{i})$ with eigenvalue $\unicode[STIX]{x1D706}_{1}=\unicode[STIX]{x1D706}_{1}^{r}\pm \text{i}\unicode[STIX]{x1D706}_{1}^{i}$ , the $r$ and $i$ superscripts denoting real and imaginary part, respectively. Hence the c.c. indication denoting complex conjugation in (4.5).

Setting the time origin sufficiently deep into the linear regime, all degrees of freedom follow an exponential oscillating decay towards the stable steady state, with shared frequency ( $f=\unicode[STIX]{x1D706}_{1}^{i}/2\unicode[STIX]{x03C0}$ ) and damping rate ( $\unicode[STIX]{x1D706}_{1}^{r}$ ). For the perturbation field, this decay is towards annihilation. In particular, the perturbation velocity components at the centre of the cavity decay following

(4.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u_{c}^{\prime }(t)=a_{1}\tilde{u} _{c}\text{e}^{(\unicode[STIX]{x1D706}_{1}^{r}+\text{i}\unicode[STIX]{x1D706}_{1}^{i})t}+\text{c.c.}=e_{u}\text{e}^{\unicode[STIX]{x1D706}_{1}^{r}t}\cos (\unicode[STIX]{x1D706}_{1}^{i}t+\unicode[STIX]{x1D713}_{u}),\\ v_{c}^{\prime }(t)=a_{1}\tilde{v}_{c}\text{e}^{(\unicode[STIX]{x1D706}_{1}^{r}+\text{i}\unicode[STIX]{x1D706}_{1}^{i})t}+\text{c.c.}=e_{v}\text{e}^{\unicode[STIX]{x1D706}_{1}^{r}t}\cos (\unicode[STIX]{x1D706}_{1}^{i}t+\unicode[STIX]{x1D713}_{v}),\end{array}\right\}\end{eqnarray}$$

where $\text{c.c.}$ denotes complex conjugation and ( $e_{u}$ , $\unicode[STIX]{x1D713}_{u}$ ) and ( $e_{v}$ , $\unicode[STIX]{x1D713}_{v}$ ), the initial envelope amplitudes and phases of the $u_{c}^{\prime }$ and $v_{c}^{\prime }$ signals, respectively, are univocally determined by the choice of the time origin through the projection $a_{1}$ of the initial flow state onto the leading eigenmode.

In order to estimate the critical value of the Reynolds number at which the Hopf bifurcation takes place ( $Re_{H}$ ) and the frequency at onset of the periodic orbit branch issued from the Hopf point, the leading complex-conjugate eigenpair must be monitored as the parameter, the Reynolds number, is increased towards the bifurcation point. This can be done by comparing the actual behaviour of $(u_{c}^{\prime },v_{c}^{\prime })$ with the expected behaviour according to (4.6), using $\unicode[STIX]{x1D706}_{1}^{r}$ and $\unicode[STIX]{x1D706}_{1}^{i}$ as fitting parameters.

A convenient way of doing so consists of defining a Poincaré section and examining the resulting discrete-time dynamical system. Here we have chosen to define the section as ${\mathcal{S}}=\{(\boldsymbol{u}^{\prime },p^{\prime }):v_{c}^{\prime }=0,\dot{v}_{c}^{\prime }<0\}$ . The Poincaré map is then defined as

(4.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\mathcal{P}}:{\mathcal{S}}\longrightarrow {\mathcal{S}},\\ \boldsymbol{a}^{P}\longmapsto {\mathcal{P}}(\boldsymbol{a}^{P})=\unicode[STIX]{x1D719}_{T^{P}}(\boldsymbol{a}^{P}),\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{a}=(\boldsymbol{u}^{\prime },p^{\prime })$ defines the state of the continuous-time dynamical system (the Navier–Stokes equations), $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D70F}}(\boldsymbol{a})$ is the flow generated by evolving state $\boldsymbol{a}$ for a lapse $\unicode[STIX]{x1D70F}$ and $\boldsymbol{a}^{P}$ is a state on ${\mathcal{S}}$ . Repeated application of the map generates a sequence $\boldsymbol{a}^{P}(k)$ with associated crossing times $t^{P}(k)$ , which are in fact obtained as a post-process of the continuous-time dynamical system evolution by second-order interpolation of ${\mathcal{S}}$ -intersections. Flight times between consecutive crossings are then obtained as $T^{P}(k)=t^{P}(k)-t^{P}(k-1)$ .

In the linear regime, equation (4.6) anticipates that ${\mathcal{S}}$ -crossings will result in the sequence

(4.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}t^{P}(k)={\displaystyle \frac{1}{\unicode[STIX]{x1D706}_{1}^{i}}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}-\unicode[STIX]{x1D713}_{v}\right)+{\displaystyle \frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D706}_{1}^{i}}}k=t^{P}(0)+{\displaystyle \frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D706}_{1}^{i}}}k=t^{P}(0)+T^{P}k,\\ {u_{c}^{\prime }}^{P}(k)=A_{u}\text{e}^{(\unicode[STIX]{x03C0}/2-\unicode[STIX]{x1D713}_{v})(\unicode[STIX]{x1D706}_{1}^{r}/\unicode[STIX]{x1D706}_{1}^{i})+(2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{1}^{i})\unicode[STIX]{x1D706}_{1}^{r}k}\cos \left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}+(\unicode[STIX]{x1D713}_{u}-\unicode[STIX]{x1D713}_{v})\right)={u_{c}^{\prime }}^{P}(0)\unicode[STIX]{x1D707}^{k},\end{array}\right\}\end{eqnarray}$$

with $T^{P}(k)=T^{P}=(2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{1}^{i})$ being the lapse between consecutive crossings and $\unicode[STIX]{x1D707}=\text{e}^{T^{P}\unicode[STIX]{x1D706}_{1}^{r}}$ the dominant multiplier (the eigenvalue of the Poincaré map with largest modulus), real positive here, of the stable fixed point of the Poincaré map, coincident in this case, by construction, with the steady state of the Navier–Stokes equations.

It is therefore to be expected that as convergence onto the steady state progresses into the linear regime, the actual sequences will asymptotically approach (4.8). Then, $\unicode[STIX]{x1D706}_{1}^{i}$ could be estimated from

(4.9) $$\begin{eqnarray}T^{P}=\lim _{k\rightarrow \infty }T^{P}(k)\longrightarrow \unicode[STIX]{x1D706}_{1}^{i}={\displaystyle \frac{2\unicode[STIX]{x03C0}}{T^{P}}}\end{eqnarray}$$

and $\unicode[STIX]{x1D706}_{1}^{r}$ might be reckoned from a least-squares exponential fit, for $k\rightarrow \infty$ , of the form

(4.10) $$\begin{eqnarray}\log \left({\displaystyle \frac{{u_{c}^{\prime }}^{P}(k)}{{u_{c}^{\prime }}^{P}(0)}}\right)=k\log \unicode[STIX]{x1D707}\longrightarrow \unicode[STIX]{x1D706}_{1}^{r}={\displaystyle \frac{1}{T^{P}}}\log \unicode[STIX]{x1D707}.\end{eqnarray}$$

The inset of figure 14(a) indicates with circles, for the sake of illustration, a collection of consecutive ${\mathcal{S}}$ -intersections as projected on the $(u_{c}^{\prime },v_{c}^{\prime })$ plane for $Re=7500$ . The main panel represents, for a range of $Re$ below critical, the sequence ${u_{c}^{\prime }}^{P}(k)/{u_{c}^{\prime }}^{P}(0)$ as obtained from figures like the one in the inset, conveniently truncated to eliminate most of the initial nonlinear transients. A least-squares fit to the exponential tails provides the multiplier.

Figure 14(b) depicts streamfunction ( ${\unicode[STIX]{x1D713}^{\prime }}^{P}$ , top panel) and vorticity ( ${\unicode[STIX]{x1D714}^{\prime }}^{P}$ , bottom panel) fields of the perturbation ( $\boldsymbol{a}^{P}$ ) on a Poincaré crossing well into the linear regime at $Re=8000$ . This flow field can be identified with the eigenmode ( $\tilde{\boldsymbol{a}}^{P}$ ) corresponding to the (positive real) dominant eigenvalue $\unicode[STIX]{x1D707}$ of the map. It is, at the same time, the imprint on ${\mathcal{S}}$ of the dynamics generated by the complex-conjugate eigenpair ( $\tilde{\boldsymbol{a}}_{1},\tilde{\boldsymbol{a}}_{1}^{\ast }$ ) of the continuous-time dynamical system. It is clear from the figure that the incipient instability concentrates in the high-shear regions of the (still stable) steady state. In particular, upcoming unsteadiness seems to originate at the strong jet that crosses the cavity through its centre and within the boundary layers that develop on the left-hand wall at either side of the stagnation point.

Figure 15(a) shows the evolution of the real ( $\unicode[STIX]{x1D706}_{1}^{r}$ , bottom panel) and imaginary ( $\unicode[STIX]{x1D706}_{1}^{i}$ , top panel) parts of the leading eigenvalue. The latter is inferred at each $Re$ from (4.9) by estimating the asymptotic saturation of $T^{P}(k)$ as the steady state is approached, while the former results directly from (4.10).

Figure 15. Hopf bifurcation. (a) Dependence of the real ( $\unicode[STIX]{x1D706}_{1}^{r}$ , bottom panel) and imaginary ( $\unicode[STIX]{x1D706}_{1}^{i}$ , top panel) parts of the leading eigenvalue. Quadratic fits of the points at $Re\geqslant 7900$ are drawn with grey dash-dotted lines. (b) Force coefficient oscillation amplitude $A_{C_{F}}$ versus $Re$ of the bifurcated periodic orbit branch. The grey line corresponds to a square-root fit of the five points closest to the saddle node.

As the critical point is approached, the dynamics around the steady state becomes extremely slow. In cases beyond $Re>8000$ , reaching the steady state was impractical, so that an approximation was obtained by averaging instead. This poses no problem as long as the estimation is done from averaging sufficiently long time series, preferably including a natural number of Poincaré crossings, already in the linear regime. The fact that the tails remain exponential are a good sign that the eigenvalue can be estimated accurately.

A quadratic fit to the last few points ( $Re\geqslant 7900$ ) of $\unicode[STIX]{x1D706}_{1}^{r}(Re)$ provides a means of estimating the critical $Re_{H}\simeq 8039.3$ . A second fit for $\unicode[STIX]{x1D706}_{1}^{i}(Re)$ to the same points forecasts a frequency $f=\unicode[STIX]{x1D706}_{1}^{i}/2\unicode[STIX]{x03C0}=0.377$ for the branch of periodic orbits at onset.

For $Re\geqslant 8040$ , the steady state has become unstable and nonlinear evolution produces a branch of periodic orbits. The oscillation frequency at $Re=8040$ is $f=0.378$ , in good agreement with the critical value for $\unicode[STIX]{x1D706}_{1}^{i}$ . The amplitude of the oscillation, however, does not seem to follow the square-root dependency that would be expected from a supercritical Hopf bifurcation. Figure 15(b) shows how the force-coefficient amplitude of the periodic orbits, measured as $A_{C_{F}}=C_{F_{\mathit{max}}}-C_{F_{\mathit{min}}}$ , varies with $Re$ . Not only does the amplitude fail to vanish at the bifurcation point, but it remains finite for $Re<Re_{H}$ disclosing a region of coexistence of the steady and periodic states. The approach to the periodic orbits is extremely slow in this regime, such that the amplitude has not been measured from converged solutions but extrapolated from a power fit of the form

(4.11) $$\begin{eqnarray}A_{C_{F}}(k)=A_{C_{F}}+\left(A_{C_{F}}-A_{C_{F}}(0)\right)|\unicode[STIX]{x1D707}|^{k},\end{eqnarray}$$

where $A_{C_{F}}$ is the asymptotic value of the orbit amplitude and $\unicode[STIX]{x1D707}$ the multiplier ( $\unicode[STIX]{x1D707}<1$ for stable periodic orbits), to the sequence

(4.12) $$\begin{eqnarray}A_{C_{F}}(k)=\max _{t}C_{F}(t)-\min _{t}C_{F}(t),\quad \text{with }t\in [t^{P}(k-1),t^{P}(k)),\end{eqnarray}$$

once the dynamics has reached the linear regime of convergence onto the periodic orbit.

Quasi-static reduction of $Re$ along the periodic orbit branch abruptly falls onto the steady-state branch somewhere in the range $Re_{FC}\in [8030,8031)$ , which bounds the hysteresis region to $Re\in [Re_{FC},Re_{H}]$ . The disappearance of the periodic orbit branch can be explained by the presence of a fold-of-cycles bifurcation, as a square-root fit $A_{C_{F}}=A_{C_{F}}^{FC}+\sqrt{Re-Re_{FC}}$ to the five lowest $Re$ periodic orbits seems to confirm. The bifurcation occurs at a critical $Re_{FC}=8031.0$ with $A_{C_{F}}^{FC}\neq 0$ . In a fold-of-cycles, the stable (nodal) branch of orbits collides with an unstable (saddle) branch and the orbits disappear. It is this unstable branch, not reachable by mere time evolution, that presumably connects with the steady state at $Re_{H}$ . The Hopf bifurcation is therefore subcritical.

An analogous study of the A-state transition from steady to periodic seems to indicate that a supercritical Hopf bifurcation occurs at $Re\simeq 13\,450$ , although a more thorough analysis would be required to pinpoint the exact location of the bifurcation and confirm its supercritical nature.

4.3 The Neimark–Sacker bifurcation

The periodic solutions branch issued from the subcritical Hopf is replaced by quasi-periodic solutions at higher $Re$ , as observed in § 3.3. The shared power spectral density peak at $f_{1}$ indicates, however, that the latter states might have evolved from the former at a Neimark–Sacker bifurcation. As a matter of fact, periodic states can be observed up to $Re\lesssim 8560$ and quasi-periodic states from $Re\gtrsim 8570$ on. Figure 16(a) shows phase map projections on the $(v_{c},C_{F})$ plane of stable solutions of the map defined with the Poincaré section ${\mathcal{S}}=\{(\boldsymbol{u},p):u_{c}=u_{c}^{P}=-0.21,\dot{u}_{c}^{\prime }>0\}$ . Periodic solutions of the Navier–Stokes equations appear as fixed points of the Poincaré map, while quasi-periodic solutions show as discrete periodic orbits. Thus, phase map trajectories are attracted by fixed points (triangles) up to $Re\leqslant 8564$ , and by periodic orbits (collection of circles) from $Re\geqslant 8565$ on. The scenario is suggestive of a direct relation between both solution branches, as the evolution of the periodic orbit centre follows smoothly from the evolution of the preexisting fixed points location. Also the amplitude of the discrete periodic orbits seems to grow from the zero amplitude of the discrete fixed points as $Re$ is increased beyond their onset value, which points to a supercritical nature of the intervening bifurcation.

Figure 16. The Neimark–Sacker bifurcation. (a) Phase map projections on the $(v_{c},C_{F})$ plane of the Poincaré maps defined by $u_{c}=u_{c}^{P}=-0.21$ , $\dot{u}_{c}>0$ . Periodic orbits are represented by triangles while circles denote quasi-periodic solutions. (b) Evolution across the bifurcation of the primary ( $f_{1},|{\hat{C}}_{F}(f_{1})|$ , circles) and modulational ( $f_{2},|{\hat{C}}_{F}(f_{2})|$ , squares) peaks in the spectrum of $C_{F}(t)$ , relative to their values at the bifurcation point. The grey dash-dotted line and circle correspond to a square-root fit of the modulational amplitude data $|{\hat{C}}_{F}(f_{2})|=k\sqrt{|Re-Re_{NS}|}$ and to the location of the bifurcation, respectively.

To establish the connection between both solution branches, the main spectral frequency peaks ( $f_{1}$ all along, $f_{2}$ for the quasi-periodic solutions) have been monitored alongside the amplitude of the solutions as measured by the spectral energy contained in the aforementioned spectral peaks ( $|{\hat{C}}_{F}(f_{1})|$ and $|{\hat{C}}_{F}(f_{2})|$ , respectively). These provide a means of quantifying the solution oscillation amplitude (for both periodic and quasi-periodic solutions) and the amplitude of the modulation of this oscillation (only for quasi-periodic). As shown in figure 16(b) (bottom panel), the oscillation frequency $f_{1}$ evolves continuously from the branch of periodic solutions to that of quasi-periodic orbits, which incorporates a modulational frequency $f_{2}$ that is incommensurate with $f_{1}$ . This places us in the scenario of a supercritical Neimark–Sacker bifurcation. Figure 16(b) (top panel) shows the evolution of the modulational amplitude $|{\hat{C}}_{F}(f_{2})|$ (squares). Supercriticality has been tested by attempting a square-root fit $|{\hat{C}}_{F}(f_{2})|=k\sqrt{|Re-Re_{NS}|}$ , with $k$ and $Re_{NS}$ as fitting parameters (grey dash-dotted line). The quality of the fit is remarkable, and provides a fair estimate of $Re_{NS}\simeq 8564.8$ at onset with zero modulational amplitude. A quadratic interpolation to the modulational frequency yields $f_{2}^{SN}\simeq 0.2106$ . Further quadratic interpolation to the last three data points before the bifurcation, corresponding to periodic solutions, results in $f_{1}^{SN}\simeq 0.3673$ and $|{\hat{C}}_{F}(f_{1})|\simeq 5.58\times 10^{-5}$ . The main oscillation amplitude remains fairly stable across the Neimark–Sacker point, while the oscillation frequency stops a decreasing trend and starts growing. Meanwhile, the modulational frequency tends to decrease from its value at the bifurcation as $Re$ is increased. The frequency ratio at the bifurcation is $R_{SN}=f_{1}^{SN}/f_{2}^{SN}\simeq 1.7441$ , not far from that of the quasi-periodic solution of § 3.3.

In order to provide further evidence of the occurrence of a Neimark–Sacker bifurcation and of its subcriticality, the modal decay of perturbations in the linear regime towards the stable periodic solutions has been analysed in a way analogous to that undertaken to elucidate the nature of the Hopf bifurcation in § 4.2. A Neimark–Sacker bifurcation is that of a fixed point of a map (or discrete-time dynamical system) whereby a complex-conjugate pair of eigenvalues crosses the unit circle in the complex plane. The map here is a Poincaré map, and the fixed point is the branch of periodic solutions.

Let $(\bar{\boldsymbol{u}},\bar{p})^{P}$ be a stable fixed point of the map (the intersection of the periodic orbit with the Poincaré section ${\mathcal{S}}$ ) and $(\boldsymbol{u},p)^{P}(k)$ the sequence generated by successive application of the map starting from some initial state within the basin of attraction of the fixed point and sufficiently close to it to be considered within the linear regime. The perturbation field will evolve as

(4.13) $$\begin{eqnarray}(\boldsymbol{u}^{\prime },p^{\prime })^{P}(k)\equiv (\boldsymbol{u},p)^{P}(k)-(\bar{\boldsymbol{u}},\bar{p})^{P}=\mathop{\sum }_{i}b_{i}(\tilde{\boldsymbol{u}}_{i},\tilde{p}_{i})^{P}\unicode[STIX]{x1D707}_{i}^{k}\xrightarrow[{t\rightarrow \infty }]{}b_{1}(\tilde{\boldsymbol{u}}_{1},\tilde{p}_{1})^{P}\unicode[STIX]{x1D707}^{k}+\text{c.c.},\end{eqnarray}$$

where $b_{i}$ is the projection of the initial perturbation on eigenmode $(\tilde{\boldsymbol{u}}_{i},\tilde{p}_{i})^{P}$ of eigenvalue $\unicode[STIX]{x1D707}_{i}$ . For convenience, the eigenmodes are sorted by decreasing modulus $1<|\unicode[STIX]{x1D707}_{1}|<|\unicode[STIX]{x1D707}_{i}|$ . Note that stability of the fixed point implies that all eigenmodes lie within the unit circle in the complex plane. For large $k$ , the evolution will align with the dominant (largest modulus) eigenmode $\unicode[STIX]{x1D707}_{1}$ .

Individual perturbation degrees of freedom will therefore asymptotically decay as

(4.14) $$\begin{eqnarray}y(k)=b_{1}{\tilde{y}}\unicode[STIX]{x1D707}_{1}^{k}+\text{c.c.}=e_{y}r^{k}\cos (k\unicode[STIX]{x1D703}+\unicode[STIX]{x1D713}_{y}),\end{eqnarray}$$

with $r$ and $\unicode[STIX]{x1D703}$ the modulus and argument of the eigenvalue $\unicode[STIX]{x1D707}_{1}=r\text{e}^{\text{i}\unicode[STIX]{x1D703}}$ , and $e_{y}$ and $\unicode[STIX]{x1D713}_{y}$ the initial envelope amplitude and phase of the $y(k)$ sequence. In particular, figure 17(a) depicts, for several values of $Re$ approaching the Neimark–Sacker bifurcation from the stable fixed points branch, the decay of the perturbation force coefficient $C_{F}^{\prime P}$ and cavity centre vertical velocity $v_{c}^{\prime P}$ . The discrete-time origin ( $k=0$ ) has been set well into the linear regime. For ease of visualisation, both quantities have been normalised by their envelope value ( $e_{C_{F}^{\prime }P}$ and $e_{v_{c}^{\prime }P}$ ) at $k=0$ , obtained from fitting (4.14) to the numerical data. The left-hand panels clearly show how the decay of $C_{F}^{\prime P}$ gets progressively less damped as $Re$ is increased towards the bifurcation point (top to bottom). It is also apparent how the decay progresses along seven distinct lines, indicating that the map nearly closes (apart from the obvious contraction onto the fixed point) every seven map iterations. The right-hand panels correspond to $(v_{c}^{\prime },C_{F}^{\prime })$ phase map projections and both the decay and the discrete quasi-period of seven are clearly evidenced. The labels at $Re=8560$ indicate the visit sequence of the seven paths, numbered from $0$ to $6$ . The phase at every map iteration does not change by $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\sim 2\unicode[STIX]{x03C0}/7$ (or $-12\unicode[STIX]{x03C0}/7$ ), as one would tend to think, but rather by $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\sim 8\unicode[STIX]{x03C0}/7$ (or $-6\unicode[STIX]{x03C0}/7$ ). The rotation number is therefore in the vicinity of $R=2\unicode[STIX]{x03C0}/\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\sim 7/4=1.75$ (or $-7/3=-2.\overline{3}$ ). The frequency ratio of the quasi-periodic solution that results from the Neimark–Sacker bifurcation, as well as its winding sign on the invariant torus, suggests that the map iteration must be considered as winding in the anticlockwise direction and that the positive values of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ and $R$ are the ones to be confronted with continuous-time data. The sign change in the apparent spiral paths (not to be mistaken for the actual winding of the map trajectories) from clockwise at $Re=8555$ to anticlockwise at $Re=8560$ is the result of traversing a rational rotation number $R=1.75$ at $Re=8557.5$ , hence the apparent locking, the straight phase map paths and the perfectly exponential discrete-time decay of each of the individual lines.

Figure 17. Modal decay of perturbations in the vicinity of the Neimark–Sacker bifurcation of the Poincaré map defined by $u_{c}=u_{c}^{P}=-0.21$ , $\dot{u}_{c}>0$ . (a) Discrete-time evolution of $C_{F}^{\prime P}$ (left) and phase map projections on the $(v_{c}^{\prime },C_{F}^{\prime })$ plane (right) of the linear decay of a perturbation onto the fixed point for several $Re<Re_{NS}$ . Here $C_{F}^{\prime P}$ and $v_{c}^{\prime P}$ have been normalised by envelope amplitude at $k=0$ . (b) Evolution of the modulus ( $r$ , bottom) and argument ( $\unicode[STIX]{x1D703}$ , top) of the dominant eigenvalue $\unicode[STIX]{x1D707}_{1}=r\text{e}^{\text{i}\unicode[STIX]{x1D703}}$ of the stable fixed point as the bifurcation is approached. Values of $r$ and $\unicode[STIX]{x1D703}$ are estimated from fits to the numerical data. Grey circles indicate the Neimark–Sacker bifurcation as estimated from quadratic fits to the rightmost few points.

Fitting (4.14) to the $C_{F}^{\prime P}$ data at several values of $Re$ approaching the Neimark–Sacker bifurcation from the stable fixed points side provides an accurate estimation of the modulus $r$ and argument $\unicode[STIX]{x1D703}$ of the dominant eigenvalue of the Poincaré map. Their $Re$ dependence is shown in figure 17(b). The intersection with the ( $r=1$ )-line of a quadratic fit to the modulus of the last few points before the bifurcation (bottom panel, dash-dotted grey line) yields an estimate of $Re_{NS}\simeq 8565.0$ for the critical point, in good agreement with the value inferred from the modulational amplitude of the quasi-periodic solution resulting from the bifurcation. A second quadratic fit, this time to the argument of the same last few points that precede the Neimark–Sacker bifurcation (top panel, dash-dotted grey line), results in $\unicode[STIX]{x1D703}_{NS}\simeq 3.6023$ . The corresponding rotation ratio is $R_{NS}\simeq 1.7442$ , again in very close agreement with that obtained from estimating the oscillation and modulational frequencies at the bifurcation point.

4.4 Nonlinear evolution of the two-torus

Figure 18 depicts phase map projections of quasi-periodic (and occasionally phase-locked) solutions winding on the two-torus created at the Neimark–Sacker bifurcation. The hole in the torus is clearly perceptible in the ( $u_{c}$ , $C_{F}$ ) projection of the solution at $Re=8600$ (black), where nonlinearity is still mild and trajectories retain the most salient features of the underlying periodic orbit that has become unstable in the Neimark–Sacker bifurcation. As $Re$ is increased (increasingly lighter shades of grey), the torus evolves nonlinearly, the lid shear force decreases and the dynamics associated with the modulation gains prominence as the imprint on the solution behaviour becomes comparable to that of the oscillation dynamics. Despite the deep transformation of the solution (as evident from phase map projections) from its origin at $Re=8600$ to $Re=10\,500$ , at the verge of becoming chaotic, the evolution is gradual and uneventful, with the sole exception of a collection of frequency-locking episodes.

Figure 18. Nonlinear evolution of the two-torus. Phase map projections on the ( $u_{c}$ , $C_{F}$ ), ( $v_{c}$ , $C_{F}$ ) and ( $v_{c}$ , $u_{c}$ ) planes of quasi-periodic solutions winding on the two-torus at $Re=8600$ , $9000$ , $9500$ , $10\,000$ and $10\,500$ (from black toward increasingly lighter shades of grey).

We have detected several such episodes, characterised by the locking of trajectories onto periodic solutions winding on the torus, as the $Re$ path in parameter space crosses a series of Arnold tongues. When entering these regions, a fold of cycles occurs on the torus. The periodic solutions observed correspond to the nodal orbit, which pairs up with a saddle orbit that is not detectable through mere time evolution. Both periodic orbits collide and are destroyed in a second fold of cycles upon leaving the Arnold tongue. Phase-locking with frequency ratios 5:3, 21:13, 50:31 or 97:60 has been found at $Re=10\,200$ , $10\,500$ (lightest grey in figure 18; see supplementary movie 4 for time evolution of vorticity field), $10\,540$ and $10\,590$ , respectively. In particular, the phase-locked orbit in the tongue at $Re\simeq 10\,540$ undergoes a period-doubling cascade that triggers chaotic dynamics.

4.5 The onset of chaotic dynamics

The onset of chaotic motion appears to follow the breakdown of the two-torus created in the Neimark–Sacker bifurcation much in the way advanced by Ruelle & Takens (Reference Ruelle and Takens1971) and later extended by Newhouse et al. (Reference Newhouse, Ruelle and Takens1978). The Afraimovich–Shilnikov theorem postulates three different ways whereby this may occur (Afraimovich & Shilnikov Reference Afraimovich and Shilnikov1983), namely (1) breakdown following some ordinary bifurcation of phase-locked limit cycles; (2) sudden transition in the presence of a homoclinic (or heteroclinic) structure; and (3) soft transition due to a gradual loss of smoothness. All three alternative paths have since been confirmed as naturally occurring in autonomous physical systems (Anischenko, Safonova & Chua Reference Anischenko, Safonova and Chua1993). While the latter path entails direct transition, the former two scenarios involve torus breakdown prior to the advent of chaos, either by the appearance of a period-doubled limit cycle that no longer lies on the invariant torus (first case) or by manifold tangles (second case).

Figure 19 shows phase map projections on the ( $v_{c}$ , $C_{F}$ ) plane of the Poincaré map defined by ${\mathcal{S}}=\{(\boldsymbol{u},p):u_{c}=u_{c}^{P}=-0.36,\dot{u}_{c}^{\prime }>0\}$ . The phase portrait at $Re=10\,520$ (figure 19 a) is clearly indicative of quasi-periodic motion on the two-torus, as the crossings fill the discrete-time limit cycle densely. Its $C_{F}$ spectrum is that of a quasi-periodic signal, with clear peaks at the two incommensurate frequencies and a bunch of secondary peaks due to nonlinear effects. At $Re=10\,530$ , the system has crossed into a phase-locked region with a 51:30 resonance, meaning that a phase-locked limit cycle of frequency $f_{l}=f_{1}/51=f_{2}/30$ is now winding on the invariant torus (figure 19 b). This periodic orbit, visible as a discrete number of crossings on the Poincaré section and with a spectrum made of peaks at all integer multiples of $f_{l}$ , corresponds to the stable solution. Pairing up with it, a saddle limit cycle, also winding on the torus, exists. Its unstable manifold, which drives the flow towards the stable limit cycle, fulfils the closure of the invariant set. The torus is still intact. The invariant torus breaks down some place between $Re=10\,530$ and $Re=10\,540$ , precisely at the period-doubling bifurcation that produces the period-doubled limit cycle we observe at $Re=10\,540$ (figure 19 c). The Poincaré map now clearly depicts a period-two discrete limit cycle that, unlike the discrete cycle at $Re=10\,530$ , winds twice before closing. The spectrum is conspicuously similar, except for the fact that new peaks have appeared between every two consecutive peaks of the original cycle, the basic frequency having been halved to $f_{l_{2}}=f_{1}/102=f_{2}/60$ . A bifurcation cascade, whose detailed analysis is beyond the scope of this study, follows, such that at $Re=10\,550$ the solution is no longer periodic but chaotic (figure 19 d). The crossings on the Poincaré section do not fall on precisely located points or along well-defined lines, but now fill a narrow band around the phase-space region where the initial period-doubled cycle appeared. The invariant set has now turned into a strange attractor. The principal peaks and their nonlinearly generated secondary counterparts are still discernible in the spectrum, but they are now accompanied by broadband noise that is well above machine precision and can by no means be ascribed to numerical truncation error.

Figure 19. Onset of chaotic motion. Phase map projections on the ( $v_{c}$ , $C_{F}$ ) of the Poincaré map defined by $u_{c}^{P}=-0.36$ (left) and spectrum $|{\hat{C}}_{F}|$ (right) at (a) $Re=10\,520$ , (b) $Re=10\,530$ , (c) $Re=10\,540$ and (d) $Re=10\,550$ .

As it happens, the chaotic transition is reversed as $Re$ is increased and a new phase-locked limit cycle corresponding to a different Arnold tongue is obtained at $Re=10\,590$ . Increasing $Re$ further, the region is left behind following a new period-doubling cascade similar to the one already described. At $Re=10\,600$ the flow is again chaotic and no evidence of reverting to quasi-periodicity has been found all the way up to $Re=13\,000$ , although this is difficult to ascertain without undertaking a thorough analysis of the evolution of the strange attractor. The evolution from the well-structured mildly chaotic attractor at $Re=10\,600$ towards the seemingly unstructured wild attractor at $Re=11\,000$ and beyond is gradual, as shown in figure 20. The vestiges of the original torus are progressively lost in an apparently featureless cloud of Poincaré crossings (figure 20 a) and the broadband noise in the spectrum steadily grows in energy (figure 20 b), gradually burying the characteristic frequency peaks.

Figure 20. Evolution of the chaotic set. (a) Phase map projections on the ( $v_{c}$ , $C_{F}$ ) of the Poincaré map defined by $u_{c}^{P}=-0.36$ and (b) spectrum $|{\hat{C}}_{F}|$ (right) at $Re=10\,600$ (black), $Re=10\,800$ (dark grey) and $Re=11\,000$ .

4.6 Bursting events away from the chaotic set

As the chaotic dynamics of the strange attractor becomes progressively wild, bursting episodes start occurring for $Re\gtrsim 12\,000$ (see supplementary movie 6) that take phase-space trajectories into a previously unexplored region where some kind of unstable state seems to dwell. These visits become recurrent and protracted as $Re$ approaches $13\,000$ . The left-hand panel of figure 21 compares the dynamics at $Re=13\,000$ (black) and $Re=11\,000$ (light grey) by displaying part of the respective $u_{c}$ time series. Mean value and fluctuations share similar characteristics part of the time, but regular excursions towards lower $|u_{c}|$ values that were absent at $Re=11\,000$ occur frequently at $Re=13\,000$ . One such excursion is indicated in the time series (dark grey) and zoomed in in the inset. The dynamics in this region seems to be governed by a periodic orbit that captures phase map trajectories along its stable manifold for a while before the solution escapes back following its unstable manifold. The pseudo-periodic nature of the signal is evident in the inset within the inset, which corresponds to a further zoom to the smoothest stage of the excursion. Phase map projection on the ( $u_{c},v_{c}$ ) plane (right-hand panel of figure 21) shows how the chaotic set remains pretty much in the same location where it was at $Re=11\,000$ (shown in light grey). However, phase map trajectories clearly incorporate frequent visits to phase-space regions not previously explored, following bursting events (one such event is highlighted in dark grey) that take the dynamics away from the location of the original chaotic set and then back. Especially significant is the occasional wandering at the top-right corner of the phase map, magnified in the inset, where trajectories seem to shadow the unstable manifold of some sort of mildly unstable state, possibly the continuation to this higher $Re$ of the state C presented in § 4.1, which would have undergone a Hopf bifurcation at some point earlier in Reynolds number. For a while, the flow in the centre of the cavity stays nearly quiescent, an indication that the diagonal jet that is characteristic of B-type states is momentarily dismantled. The most quiet stage of the approach is shown in the second inset, where two full pseudo-periods have been indicated (black) to convey the dynamic properties of the underlying state that trajectories appear to orbit for a while. We conjecture that the aforementioned state-C saddle solution might be responsible for piercing the chaotic attractor at a slightly higher $Re$ in a boundary crisis. We shall not explore the issue further, as the resolution and the computational resources required to fully clarify the situation are well beyond the scope of this study, but leave it for future investigation.

Figure 21. Chaotic flow dynamics at $Re=13\,000$ (black), as compared to $Re=11\,000$ (light grey). The left-hand panel contains $u_{c}$ time series, while the right-hand panel displays phase map projections on the ( $u_{c},v_{c}$ ) plane. A prototypical excursion (dark grey) is highlighted in both panels and the dashed boxes are zoomed in in the respective insets. A further zoom, shown in the second set of insets, reveals the pseudo-periodic nature of the signal in the region of phase space visited during burst events.