2.1. Governing equations

The non-dimensional velocity $\boldsymbol {v}=(v_r,v_{\theta },v_z)$ and the non-dimensional pressure ${\rm \pi}$ obey the incompressible Navier–Stokes equations (2.1–2.4). Throughout the whole paper we assume that the flow is strictly axisymmetric. We consider the Navier–Stokes equations in the cylindrical coordinate system $(\boldsymbol{e}_r,\boldsymbol{e_{\theta }},\boldsymbol{e_{z}})$

(2.1)$$\begin{gather} \frac{\partial v_r}{\partial t} + v_r \frac{\partial v_r}{\partial r} + v_z \frac{\partial v_r}{\partial z} -\frac{v^2_{\theta}}{r}={-}\frac{\partial {\rm \pi}}{\partial r} + \frac{1}{Re}\left(\frac{1}{r}\frac{\partial}{\partial r} \left(r\frac{\partial v_r}{\partial r}\right) - \frac{v_{r}}{r^2} + \frac{\partial^2 v_r}{\partial z^2}\right), \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial v_{\theta}}{\partial t} + v_r \frac{\partial v_{\theta}}{\partial r} + v_z \frac{\partial v_{\theta}}{\partial z} +\frac{v_rv_{\theta}}{r}= \frac{1}{Re}\left(\frac{1}{r}\frac{\partial}{\partial r} \left(r\frac{\partial v_{\theta}}{\partial r}\right) - \frac{v_{\theta}}{r^2} + \frac{\partial^2 v_{\theta}}{\partial z^2}\right), \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial v_z}{\partial t} + v_r \frac{\partial v_z}{\partial r} + v_z \frac{\partial v_z}{\partial z} ={-}\frac{\partial {\rm \pi}}{\partial z} + \frac{1}{Re}\left(\frac{1}{r}\frac{\partial}{\partial r} \left(r\frac{\partial v_z}{\partial r}\right) +\frac{\partial^2 v_z}{\partial z^2}\right), \end{gather}$$

(2.4)$$\begin{gather}\frac{1}{r} \frac{\partial (r v_r)}{\partial r} + \frac{\partial v_z}{\partial z}=0. \end{gather}$$

The coupled system of partial differential equations is complemented with the no-slip boundary conditions expressed as (2.5)

(2.5)\begin{equation} \begin{cases} \boldsymbol{v}=r\boldsymbol{e_{\theta}}, & \text{at the rotor}\ (z=0), \\ \boldsymbol{v}=\varGamma \boldsymbol{e_{\theta}}, & \text{at the shroud}\ (r=\varGamma), \\ \boldsymbol{v}=\boldsymbol{0}, & \text{on the stator}\ (z=1). \end{cases} \end{equation}

2.2. Discretisation

The continuous problem is discretised using a second-order finite volume method on a staggered grid. The details of the discretisation as well as the staggered arrangement can be found in Appendix A of Faugaret et al. (Reference Faugaret, Duguet, Fraigneau and Martin Witkowski2022). Two types of mesh, uniform and non-uniform, are used. The non-uniform mesh is refined in the regions near the rotor, the stator and the rotating shroud. Details on the non-uniform mesh set-up can be found in Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024).

Details on the spatial discretisation used are given in table 1. In both uniform and non-uniform cases, the number of pressure cells in the radial direction $r$ (respectively the axial direction $z$) is noted $N_r$ (respectively $N_z$). For the application of the boundary conditions a layer of ghost cells is added outside the physical domain.

Table 1. Mesh resolutions used in the study (DOF is degrees of freedom).

2.3. Numerical methods

The nonlinear system of (2.1–2.4) admits for all $Re$ a steady solution. Once the system is discretised, the solution of the large algebraic nonlinear system of equations is determined numerically using a Newton–Raphson algorithm (see e.g. Appendix B in Faugaret et al. Reference Faugaret, Duguet, Fraigneau and Martin Witkowski2022). The $O(4N_rN_z)$ unknowns are the velocity and pressure values at each discretisation point. Due to the presence of the thin boundary layers close to each disk, the set of equations resulting from the discretisation of the continuous system is in general poorly conditioned. Sparse direct solvers will be therefore preferred over the iterative solvers for solving the linear systems in each Newton–Raphson iteration.

A technical remark concerning the geometry can be made at this point: at the junction between the rotating shroud and the stationary disc, the sudden change in $u_{\theta }$ imposed in the boundary condition induces a discontinuity of the velocity field in the $(r,z)=(\varGamma,1)$ corner. This singularity can have an impact on the accuracy of the numerical solutions. In order to avoid singular boundary conditions and the well-known associated Gibbs phenomenon, several teams using pseudospectral solvers (Serre et al. Reference Serre, Tuliszka-Sznitko and Bontoux2004; Lopez et al. Reference Lopez, Marques, Rubio and Avila2009) proposed to smooth out the discontinuous boundary condition between the stator and the shroud. Whenever finite volume discretisation is used the singular corner does not degrade the second order of spatial accuracy of the scheme, as demonstrated in Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024).

The linear stability of the base flow is evaluated using the Arnoldi method based on a well-validated ARPACK package (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998). A generalised eigenvalue problem can be constructed using the Jacobian matrix used in the procedure for finding the base flow stemming from linearisation of the governing equations around a base flow (§ 4). For a generalised eigenvalue problem

(2.6)\begin{equation} \boldsymbol{A} \boldsymbol{w}=\lambda \boldsymbol{B} \boldsymbol{w}, \end{equation}

the shift-and-invert method finds a subset of eigenvalues closest to a complex shift $s$ through repeated Arnoldi iteration

(2.7)\begin{equation} \nu \boldsymbol{w}=( \boldsymbol{A}-s \boldsymbol{B})^{{-}1}\boldsymbol{B}\boldsymbol{w}, \end{equation}

where the original eigenvalues $\lambda$ can be retrieved with

(2.8)\begin{equation} \nu=\frac{1}{\lambda-s}. \end{equation}

The time integration of the governing equations uses a prediction–projection algorithm in the rotational pressure correction formulation, as described in Guermond, Minev & Shen (Reference Guermond, Minev and Shen2006). The prediction step combines a backwards differentiation formula 2 scheme for the diffusion terms and an explicit treatment of the convective terms, resulting in a Helmholtz problem for each velocity component increment. These Helmholtz problems are solved using an alternating-direction implicit method in incremental form which preserves the second-order accuracy in time. In the projection step the velocity field is projected onto the space of divergence-free fields by solving a Poisson equation for the pressure. This Poisson equation is solved using a direct sparse solver. Once the base flow solution is found, the time-stepping code can be adapted to evolve the perturbation to the base flow rather than the full velocity field itself, at the expense of a triple evaluation of the convective terms.

4.1. Optimal response theory

In this subsection the linear response of the flow to the forcing is analysed using optimal response theory. Following Cerqueira & Sipp (Reference Cerqueira and Sipp2014) we conduct the input/output analysis and show that the optimal response of the flow is for most relevant values of $Re$ in the shape of circular rolls, and that associated with high levels of optimal gain. The nonlinear system of (2.1–2.4) can be linearised around the base flow when perturbation velocities are small enough. If a forcing field $\boldsymbol{f}$ is introduced, the resulting system for the perturbation field $(\boldsymbol{u},p)$ can be rewritten

(4.1)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{U_b}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}+ (\boldsymbol{u}\boldsymbol{{\cdot}} \boldsymbol{\nabla})\boldsymbol{U_b}={-}\boldsymbol{\nabla}p + \frac{1}{Re}\boldsymbol{\nabla}^2 \boldsymbol{u} + \boldsymbol{f}, \end{gather}$$

(4.2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0. \end{gather}$$

The coupled system in (4.1)–(4.2) is linear in $\boldsymbol{u}$, $p$ and also $\boldsymbol{f}$, which suggests to use the resolvent formalism. The theory by Trefethen & Embree (Reference Trefethen and Embree2005) forms the ideal framework except that the original linear system needs be rewritten in a form $\partial \boldsymbol{q}/\partial t=\boldsymbol{L q} + \boldsymbol{f_q}(t)$, where the unknown field $\boldsymbol{q}$, a linear operator $\boldsymbol{L}$ and a forcing $\boldsymbol{f_q}(t)$ need to be specified. We introduce the new variable $\boldsymbol {q}=(u_r,u_{\theta },u_z,p)$ which contains the values of the fields $u_r$, $u_{\theta }$, $u_z$ and $p$ at all the points of the discretised domain. The size of $\boldsymbol {q}$ is $O(4N_rN_z)$. We also introduce the rectangular prolongation operator $\boldsymbol{P}$ of size $O(4N_rN_z \times 3N_rN_z)$ which maps $\boldsymbol{u}$ into $\boldsymbol{q}$, so that $\boldsymbol {P}^{T}\boldsymbol {q}=\boldsymbol {u}$ with the property $\boldsymbol {P}^{T}\boldsymbol {P}=\boldsymbol {I}$, see e.g. Jin, Symon & Illingworth (Reference Jin, Symon and Illingworth2021). The linear system (4.1–4.2) can then be rewritten into the new form

(4.3)\begin{equation} \boldsymbol{B}\frac{\partial \boldsymbol{q}}{\partial t}=\boldsymbol{A}\boldsymbol{q}+\boldsymbol{P}{\boldsymbol{f}}, \end{equation}

where $\boldsymbol {B}=\boldsymbol {P}\boldsymbol {P}^{T}$. After a Fourier transform in time, each Fourier component $\hat {\boldsymbol {q}}$ of $\boldsymbol {q}$ satisfies

(4.4)\begin{equation} ({\rm i} \omega \boldsymbol{B}-\boldsymbol{A}) \hat{\boldsymbol{q}}=\boldsymbol{P}\hat{\boldsymbol{f}} , \end{equation}

resulting in the Fourier components $\hat {\boldsymbol {u}}$ of the velocity field $\boldsymbol {u}$ as the action of a matrix $\boldsymbol {R}$ on the forcing $\hat {\boldsymbol {f}}$

(4.5)\begin{equation} \hat{\boldsymbol{u}}=\boldsymbol{R}\hat{\boldsymbol{f}}, \end{equation}

where

(4.6)\begin{equation} \boldsymbol{R}=\boldsymbol{P}^{T}( {\rm i} \omega \boldsymbol{B}-\boldsymbol{A})^{{-}1}\boldsymbol{P}, \end{equation}

is the resolvent operator associated with the (real) angular frequency $\omega$.

An optimal gain can be evaluated by identifying an optimal forcing for a suitable norm of the resolvent $\boldsymbol {R}$. A positive symmetric linear operator $\boldsymbol {Q} \neq \boldsymbol {I}$, associated with the cylindrical coordinate system, can be used to define the inner product

(4.7)\begin{equation} \langle\boldsymbol{u},\boldsymbol{v}\rangle_{\boldsymbol{Q}}=\boldsymbol{u}^*\boldsymbol{Q}\boldsymbol{v} =\int_0^1 \int_0^{\varGamma} (u_r^*v_r+u_{\theta}^*v_{\theta}+u_z^*v_z)r \,{\rm d}r \,{\rm d}z, \end{equation}

where the asterisk denotes complex conjugate. The associated vector norm is defined by

(4.8)\begin{equation} \| \boldsymbol{u} \|_{\boldsymbol{Q}} =\sqrt{\boldsymbol{u}^*\boldsymbol{Z}^{T}\boldsymbol{Z}\boldsymbol{u}}=\| \boldsymbol{Z}\boldsymbol{u} \|_2, \end{equation}

where $\boldsymbol {Q}=\boldsymbol {Z}^{T}\boldsymbol {Z}$. Here, $\boldsymbol {Q}$ can be also used to define the following norm for the resolvent operator:

(4.9)\begin{equation} \| \boldsymbol{R} \|_{\boldsymbol{Q}}=\| \boldsymbol{Z}\boldsymbol{R}\boldsymbol{Z}^{{-}1} \|_2=\sigma_{1} (\boldsymbol{Z}\boldsymbol{R}\boldsymbol{Z}^{{-}1}), \end{equation}

where $\sigma _{1}$ denotes the largest singular value in the SVD (singular value decomposition). Note that finding the largest singular value of $\boldsymbol {ZRZ^{-1}}$ is equivalent to finding the largest eigenvalue of the eigenvalue problem (Cerqueira & Sipp Reference Cerqueira and Sipp2014)

(4.10)\begin{equation} \boldsymbol{R}^*\boldsymbol{Q}\boldsymbol{R}\hat{\boldsymbol{f}}=\lambda\boldsymbol{Q}\hat{\boldsymbol{f}}, \end{equation}

where

(4.11)\begin{equation} G=\lambda=\sigma_{1}^2=\frac{\hat{\boldsymbol{u}}^*\boldsymbol{Q}\hat{\boldsymbol{u}}}{\hat{\boldsymbol{f}}^*\boldsymbol{Q}\hat{\boldsymbol{f}}}, \end{equation}

is interpreted as the optimal energy gain.

4.2. Optimal response: results

For a range of angular frequencies $\omega$ the eigenvalue problem (4.10) is solved numerically in MATLAB with the eigs() function. The resulting eigenvalue $G=\sigma _1^2$ from (4.11) is plotted in figure 3(a) as a function of the (real) angular frequency of the forcing. The optimal values $\sigma _{max}^2$ over all angular frequencies obtained for various values of $Re$ are compared in figure 3(b), with the case $\omega =0$ singled out. The optimal value of $G$ gain over all $\omega$ values, i.e. the optimal gain, is also listed in the table 2 together with the values obtained with different mesh resolutions.

Figure 3. Optimal response for bulk-based forcing. (a) Optimal energy gain as a function of forcing angular frequency $\omega$, given by the largest singular value of the resolvent operator (4.6). (b) Optimal gain across all $\omega$ values as a function of $Re$. Numerical resolution R0. The inset shows the optimal forcing angular frequency as listed in table 2.

Table 2. Optimal forcing gain. The last column shows the optimal angular frequency $\omega$ associated with the largest energy gain (the computations are performed with a step of 0.1 in angular frequency). Four mesh resolutions (see table 1) are used to find the optimal gain value. While for $Re\approx 250$ results can be considered as satisfactory for resolution R1, increased resolution is needed for larger $Re$. This is due to increasingly thin boundary layers required to resolve the optimal forcing mode.

For low $Re<150$ the optimally amplified frequency is always $\omega =0$. This corresponds to a steady forcing $f$ in (2.1–2.4). Forcing a given flow at $\omega =0$ can be interpreted in different complementary ways. It can first be understood as the smooth $\omega \rightarrow 0$ limit of a given harmonic forcing at frequency $\omega$. It can also be understood, in the unforced problem, as the steady streaming component associated with the nonlinear self-interaction of arbitrary oscillatory perturbations (see e.g. Mantič-Lugo, Arratia & Gallaire Reference Mantič-Lugo, Arratia and Gallaire2014). In both cases the optimal response at $\omega =0$ is interpreted as an optimal steady mean flow correction, which justifies the special focus on $\omega =0$. As seen in figure 4 the most amplified steady forcing is always localised in the region near the axis. For $Re>1800$ it inherits the characteristics of the base flow in the sense that it is composed of two boundary layers and with an invariant core in between.

Figure 4. Azimuthal perturbation velocity $u_{\theta }(r,z)$ for optimal forcing and optimal response for vanishing forcing frequency $\omega =0$. For $Re=70$ and 150 the optimal forcing at $\omega =0$ is also optimal across all $\omega$ values. The colour map spans 8 equal subintervals of $[0,0.5]$. Both optimal forcing and response are normalised such that $\| \hat {\boldsymbol {f}} \|_{\boldsymbol {Q}}=\| \hat {\boldsymbol {u}} \|_{\boldsymbol {Q}}=1$.

For $Re\approx 150$ the optimally forced structures change radically as unsteady ($\omega \neq 0$) forcing takes over steady forcing. This is clearly seen in figure 3(b), where for $Re \gtrsim 150$ non-zero forcing frequencies $\omega$ start to dominate the optimal gain curve. The structure of the associated optimal forcing lies entirely within the Bödewadt layer, see figure 5, which focuses on $\omega \approx 1.7$ close to the optimal angular frequency (see table 2). For all velocity components these structures are tilted by the mean shear into the streamwise direction. While the amplitude of optimal forcing is similar in the three spatial components, the corresponding optimal response may differ from component to component (Jovanović & Bamieh Reference Jovanović and Bamieh2005), and in the present geometry it is clearly dominated by the azimuthal component. The signature of the circular rolls can be also seen in the azimuthal response although the corresponding structures should more realistically be labelled as azimuthal streaks. Their position ($r\approx 4\unicode{x2013}6$) is perfectly consistent with the experimental observations by Schouveiler et al. (Reference Schouveiler, Le Gal and Chauve2001) and the numerical studies of Lopez et al. (Reference Lopez, Marques, Rubio and Avila2009) and Do et al. (Reference Do, Lopez and Marques2010). Most importantly, and this constitutes one of the main findings of the current work, the optimal response is in the shape of circular rolls in the Bödewadt layer.

Figure 5. Optimal forcing and response for $Re=250$. (a–f) Perturbation velocity components $u_r$, $u_{\theta }$ and $u_z$ for the optimal angular frequency $\omega =1.7$. The colour map spans 8 equal subintervals of $(-0.15\ 0.15)$ for the optimal forcing, and $(-0.6\ 0.6)$ for the optimal response plots. The optimal forcing and response are normalised such that $\| \hat {\boldsymbol {u}} \|_{\boldsymbol {Q}}=\| \hat {\boldsymbol {f}} \|_{\boldsymbol {Q}}=1$.

The evolution of the optimally forced structure is now described as $Re$ increases beyond 250. It is first noted that, as shown in figure 3(b), the optimal gain $G_{opt}(Re)$ grows exponentially with $Re$. Other examples of an exponential scaling of gain include a backward-facing step (Blackburn, Barkley & Sherwin Reference Blackburn, Barkley and Sherwin2008) and an oscillatory pipe flow (Xu, Song & Avila Reference Xu, Song and Avila2021). We note that the experimental study of Gauthier et al. (Reference Gauthier, Gondret and Rabaud1999), performed for $\varGamma \approx 20$, has suggested a supercritical bifurcation as the origin of the observed rolls. Our observations do not corroborate this hypothesis since the angular frequency content of the response features that a continuum of frequencies can be excited externally, at odds with the dominant frequency stemming from a Hopf bifurcation. Moreover, the steep exponential increase of $G(Re)$ reported above might be responsible, in experimental conditions where external forcing stays uncontrolled, for the apparent bifurcating behaviour where the amplitude of the response increases rapidly with $Re$. As $Re$ continues to increase, as shown in figure 6 the optimal forcing evolves with $Re$ from a wide support within the Bödewadt layer to thinner structures in both the Bödewadt and the shrouding wall, and even, for $Re=3000$, also in the Ekman layer. For such high $Re$ values the respective supports of the forcing and the response are almost disjoint.

Figure 6. Azimuthal perturbation velocity $u_{\theta }(r,z)$ for the optimal forcing and response for $\omega \neq 0$. For ${Re=250}$, 1800 and 3000 the optimal forcing is also optimal across all $\omega$ values. The colour map spans 8 equal subintervals of $(-0.15\ 0.15)$ for the optimal forcing and $(-0.6\ 0.6)$ for the optimal response.

The results above unambiguously point towards an unsteady response in the shape of circular rolls for all $Re \gtrsim 150$. We emphasise that the crossing of the gain curves for $\omega =0$ and $\neq 0$ in figure 3(b) does not define a threshold value for $Re$ because, as for any linear receptivity mechanism, the response depends linearly on the spectral content of the forcing history. Defining a threshold value for $Re$ is demanding because it is highly dependent on the amplitude levels that an experimentalist can detect in practice. By focusing on the value of the optimal gain ($G=4.9\times 10^2$) at $Re=250$ the following evaluation can be made. Any parasitic vibration present in the experiment projected on the orthogonal basis of the optimal modes will have a non-zero component that will be optimally forced. If the amplitude of this component is, say, of order $O(10^{-2})$ it will be amplified by the linear mechanism by $\sqrt {G}\approx 20$ to yield an $O(10^{-1})$ response, which can be detected in experiments.

4.3. Boundary forcing

While the previous section offers an elegant formal manner to explain the circular rolls from (linear) optimality arguments, we note that, at the experimentally relevant values of $Re$, the optimal forcing protocol corresponds to a force field that needs to be applied to the flow away from the solid boundaries. Set-up imperfections are expected to induce forcing at the boundaries rather than away from them. For this reason, it is unlikely that the observed linear response corresponds to truly optimal forcing, and it might be more relevant to concentrate on sub-optimal forcing. Instead, we consider a forcing protocol which without being energetically optimal affects directly the flow through unsteady motion of the boundaries (and is hence not a suboptimal from the previous bulk-based optimisation). Modulations of the instantaneous (dimensional) angular velocity are considered in the form $\varOmega (t)=\varOmega _0(1+A\ \varepsilon (t))$, where $\varepsilon (t)$ represents a normalised unsteady forcing and $A \ge 0$ is a measure of its amplitude. This is similar to Lopez et al. (Reference Lopez, Marques, Rubio and Avila2009), Poncet, Serre & Le Gal (Reference Poncet, Serre and Le Gal2009) and Do et al. (Reference Do, Lopez and Marques2010), except that $\varepsilon (t)$ is not monochromatic. The Reynolds number, now based on $\varOmega _0$ rather than $\varOmega$ (which depends on time), remains by definition unaffected by the value of $A$.

4.3.1. Time integration

The time modulation is simulated in practice by updating, at the end of each timestep after the prediction–projection step, the value of $v_{\theta }$ imposed on the shroud and the rotor as

(4.12)\begin{equation} v_{\theta}(r,t)=r(1+A\ \varepsilon(t)). \end{equation}

The modulation $\varepsilon (t)$ can be chosen in many ways. For most of this study, it is chosen as a Gaussian white noise of zero mean and standard deviation 1. The parameter $A$, which multiplies the white noise signal, is therefore the root mean square value of the forcing. All simulations are initiated with a zero perturbation field.

Figure 7 shows representative azimuthal velocity snapshots obtained from nonlinear time integration of the forced flow. As will be clear later in § 6, for the aspect ratio $\varGamma =10$ and $Re<300$, nonlinearity plays a negligible role and the observations are independent of whether the time integration is linear or nonlinear. This is why those results are discussed in the context of linearly optimal results. The response to the unsteady white noise forcing is localised, for low $Re\approx 70$, next to the rotor and the shroud. For $Re \gtrsim 150$ and beyond, wavetrains of increasingly smaller circular rolls form inside the Bödewadt layer, as is visible for $Re=250$ and $Re=400$.

Figure 7. Nonlinear time integration in the presence of white noise forcing at the $z=0$ and $r=\varGamma$ boundary with amplitude $A=10^{-2}$. Snapshots of the azimuthal perturbation velocity $u_{\theta }(r,z)$. Panels show (a–d) ${Re=70},\ 150,\ 250,\ 400$. A wavetrain of circular rolls is seen for $Re>250$.

Interestingly, the response to boundary forcing is similar to the response to the optimal forcing in the sense that the response consists also of wavetrains of circular rolls located inside the Bödewadt layer. Due to the wide frequency content this response is, however, more widespread than in the case of optimal forcing (e.g. for $Re=250$ the wavetrain is detected for $r\in (2,7)$). A pronounced response near the rotor and especially the shroud is also seen, as a signature of the imposed forcing.

4.3.2. Resolvent approach to boundary forcing

The connection between the linearised time integration and the resolvent approach from (4.4) is the Fourier transform of the linearised governing equations. Selecting the additive forcing term $\hat {\boldsymbol {f}}$ corresponding to the forcing protocol of (4.12), and then solving (4.5) for $\hat {\boldsymbol {u}}$ will therefore yield directly the Fourier transform of the response. It can be in turn compared against the Fourier transform of the probe signals extracted from linearised time integration. The same analogy has been exploited by Cerqueira & Sipp (Reference Cerqueira and Sipp2014) in order to validate their input/output analysis. In our case, beyond immediate validation, this comparison is useful as it shows the quantitative consistency between boundary forcing (although based on linear time integration rather than the nonlinear integration allowing for spectral mixing) and the inherently bulk-based resolvent approach.

We therefore select a specific, non-optimal forcing term $\boldsymbol {P}\hat {\boldsymbol {f}}$ corresponding to the Fourier transform of the boundary forcing (4.12). Solving (4.4) with prescribed $\boldsymbol {P}\hat {\boldsymbol {f}}$ yields $\boldsymbol {P}\hat {\boldsymbol {u}}$, the real part of which is plotted in figure 8. The (squared) $\boldsymbol {Q}$-norm of $\hat {\boldsymbol {u}}$ is shown in the same figure. Due to the forcing being suboptimal, its dependence on $\omega$ differs from the gain curve in figure 3(a). Again, starting with $Re=250$, a hump visible in figure 8(a) around $\omega =2$ marks the preferred response of the flow in the shape of rolls.

Figure 8. Response to the boundary forcing obtained by the resolvent approach in § 4.3.2. (a) The $\boldsymbol {Q}$-norm (squared) of the velocity response. (b–e) Real part of $\hat {u}_{\theta }(r,z)$ for $\omega =2$ for $Re=70,\ 150,\ 250,\ 400$. The $\hat {\boldsymbol {u}}$ response is normalised so that $\| \hat {\boldsymbol {u}}\|_{\boldsymbol {Q}}=1$. Colour map spans $-0.3,(0.075),0.3$.

More insight into the preferred response frequencies of the flow is possible through the analysis of the probe signals extracted from linear time integration. The raw signals, corresponding to four radial positions $r=1,3,5,7$ along the Bödewadt layer, are shown in figure 9(a). The forcing is again white in time with equal amplitude for each frequency. The Fourier transform of the probe signal confirms, as already visible from the fields in figures 7 and 8, that the strength of the response to boundary forcing depends on the radial position. The radially inwards flow in the Bödewadt layer suggests that the distance to the shrouding wall, and therefore the varying thickness of the boundary layer, are the physically meaningful variables to explain this dependency, as also suggested by Gauthier et al. (Reference Gauthier, Gondret and Rabaud1999) (we will, however, stick to the variable $r$ for commodity). The preferred response frequency evolves also with the radial position. It is close to $\omega =0$ for the probe at $r=1$, but close to $\omega =3$ for the probe $r=7$. The Fourier amplitude can be directly compared with the pointwise amplitude of $\hat {\boldsymbol {u}}$, as mentioned earlier. The convincing overlap of both data (see figure 9b) demonstrates the equivalence between the linear time integration and the direct resolvent solve based on (4.5).

Figure 9. Response to boundary forcing by Gaussian white noise with $A=10^{-2}$, $Re=200$. (a) Time series of perturbation velocity $u_{\theta }$ obtained from time integration for different probes at varying $r$ inside the Bödewadt layer $z=0.9$. (b) Comparison between the time integration and the resolvent approach to boundary forcing in frequency space. Fourier amplitude spectrum of these time series (solid lines) vs ${|\hat {u}_{\theta }|}$ associated with boundary forcing, computed from resolvent analysis (dashed lines) and evaluated at the same spatial positions. Numerical parameters $dt=4\times 10^{-3}$, $nt=20\,000$, $d\omega =0.157$, $\omega _{Nyquist}=785$ corresponding to the timestep, number of timesteps, sampling angular frequency and Nyquist angular frequency (the highest angular frequency that can be reliably measured equals ${{\rm \pi} }/{dt}$), respectively. Only the later half of the signal was used for the calculation of the Fourier transform ($40< t<80$).

An agreement between current results and the numerical study of Do et al. (Reference Do, Lopez and Marques2010) is noted when comparing figures 9 and 6 from Do et al. (Reference Do, Lopez and Marques2010). Both spectra are characterised by a broad curve centred around $\omega =2.4$. The position of the rolls around the mid-radius of the cavity also agrees while comparing figures 8 and 5 from Do et al. (Reference Do, Lopez and Marques2010).

4.4. Computation of the pseudospectrum

A fundamental tool in the analysis of non-normal systems is the pseudospectrum, which is a generalisation of the concept of the (eigen)spectrum. Recall that the spectrum of the linearised system $\dot {\boldsymbol{q}}=\boldsymbol {L}\boldsymbol{q}$ is the (complex) set of eigenvalues of the associated linear operator. If a complex number $\omega$ belongs to the spectrum $\sigma$, the resolvent $({\rm i}\omega \boldsymbol {I}-\boldsymbol {L})^{-1}$ (or expressed in the present case by (4.6)) is undefined, whereas it is continuous and smooth as a function of $\omega$ in the immediate neighbourhood of its pole at ${\rm i}\omega$. Following Trefethen & Embree (Reference Trefethen and Embree2005) the $\varepsilon$-pseudospectrum $\sigma _{\varepsilon }$ is defined as the complex set where the norm of the resolvent operator exceeds a given value $\varepsilon ^{-1}$, with $\varepsilon$ a potentially small real number

(4.13)\begin{equation} \sigma_{\varepsilon}=\{\omega \in \mathbb{C}, \| \boldsymbol{R} \|_{\boldsymbol{Q}}>\varepsilon^{{-}1}\}. \end{equation}

It includes the point spectrum and can be seen as its generalisation. In the current work, it is mainly used as an indicator of the strong non-normality of the underlying linearised operator. Owing to the definition 4.13, the cut of the pseudospectrum through the real axis directly yields the gain curves plotted in figure 3. The pseudospectrum is computed at each point, analogously to the optimal gains described in § 4.1, except that $\omega$ is allowed to be complex. Note that such computations rely on the shift-and-invert algorithm, itself dependent on a shift parameter $s$. Contours of the pseudospectrum are reported in figure 10. As expected for a strongly non-normal operator, the isocontours of $\sigma _{\varepsilon }$ do not form concentric circles around the eigenvalues yet the contours encircle more than one eigenvalue. Still, very close locally to a given eigenvalue, the isocontours found for $Re=200$ form closed loops around specific eigenvalues, as can be seen in the top right panel in figure 10. Upon increasing $Re$, the non-normality of the linearised operator increases and the background level in the pseudospectrum grows, as shown in figure 10.

Figure 10. Pseudospectrum levels $\log _{10}\| \boldsymbol {R} \|_{\boldsymbol {Q}}$ in the complex plane where $\lambda _r$ and $\lambda _i$ are the real and imaginary parts of the complex number ${\rm i}\omega$. Panels show (a,b) $Re=200$, (c,d) $Re=3000$. The right plots correspond to the red insets in the left plots, where three different values of the shift $s$ in the shift-and-invert algorithm have been used.

Another useful piece of information from the pseudospectra is the sensitivity of an eigenvalue to arbitrary perturbations of the operator $\boldsymbol {L}$ (Cerqueira & Sipp Reference Cerqueira and Sipp2014; Brynjell-Rahkola et al. Reference Brynjell-Rahkola, Shahriari, Schlatter, Hanifi and Henningson2017). As shown in chapter 28 of Trefethen & Embree (Reference Trefethen and Embree2005) for a convection–diffusion operator, high levels of non-normality typically cause iterative Arnoldi methods to converge to false eigenvalues. This is easily verified by comparing the effect of different shift values $s$ in the shift–inverse Arnoldi iteration. Here, as in Cerqueira & Sipp (Reference Cerqueira and Sipp2014), a very large level of pseudospectral contours, typically above 11, prevents ARPACK from accurately converging all the eigenvalues. This algorithmic sensibility is not only a numerical convergence problem, it also highlights a physically ambiguous situation. Individual eigenvalues, when they are not robust, do not yield a specific contribution to the dynamics. In particular, the associated eigenfrequencies are not resonant, and a large response to forcing is obtained even for forcing frequencies far from the eigenfrequencies (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). This is associated with the presence of a large hump (rather than narrow peaks) in the gain curve. In such a case, frequent in shear flows such as jets or boundary layers, one speaks of pseudo-resonance, see e.g. Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013). Non-robust eigenvalues are easily spotted as soon as their location in the complex plane depends on the shift $s$ chosen by the user, see figure 10(d). Conversely, the most unstable eigenvalue $\lambda _r=0.011,\ \lambda _i=1.95$ at $Re=3000$ in figure 10 can be trusted, and the growth rate $\lambda _r$ predicted by the computation of the eigenvalues will correspond to the growth rate of instability upon time linearised integration.

Finally, the pseudospectra can also be post-processed to extract a lower bound on the maximal transient growth associated with given initial perturbations. Following chapter 14 in Trefethen & Embree (Reference Trefethen and Embree2005), the Kreiss constant is the lower bound of the maximal growth of the initial perturbation and is defined by

(4.14)\begin{equation} \mathcal{K}=\sup_{\varepsilon>0}\frac{\alpha_{\varepsilon}}{\varepsilon}, \end{equation}

where $\alpha _{\varepsilon }$ is the pseudospectral abscissa (the maximal real part of the $\varepsilon$ contour in the complex plane). Computation of $\mathcal {K}$ for the pseudospectra in figure 10 gives $\mathcal {K}_{Re=200}=1.4\times 10^3$ and $\mathcal {K}_{Re=3000}=2.2\times 10^8$, suggesting very large growth potential, as demonstrated by Daube & Le Quéré (Reference Daube and Le Quéré2002) and Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024) for the same set-up.

6.1. Effect of the nonlinearity

We begin by describing the effect of increasing the forcing amplitude on the gain curves, such as displayed in figure 9(b). Since nonlinear effects are expected to become more important with increasing $Re$, we first illustrate their effect by choosing two representative values of $Re=600$ and $2100$. Rather than choosing two different values of the forcing amplitude in (4.12), we choose only one such value $A=10^{-2}$ and compare respectively linear and nonlinear temporal simulations started from the same initial condition $\boldsymbol{u}(t=0)=0$. The case of $Re=600$ is investigated in figure 12 while $Re=2100$ is analysed in figure 13. For $Re=600$, time series taken from probes located in the Bödewadt layer reveal visually no difference between linear and nonlinear simulations. The Fourier transforms of these signals, interpretable as local gain curves, reveal, however, that the fastest scales (large $\omega \gtrsim 4$) are affected by nonlinear effects at least for the radii $r=3$ and $5$. At other radii neither the fast nor the slow scales are affected by nonlinearity. For $Re=2100$, the situation is clearer: the velocity time series look very different and their Fourier transforms, unsurprisingly, differ at all radial locations. In all cases nonlinearity manifests itself by a slower, exponential-like decrease of the Fourier amplitudes with increasing $\omega$, interpretable as a cascade in frequency space. It is also nonlinear interactions between the frequencies that are responsible for the irregular shape of the spectrum. This is because the timestep in the simulation cannot be commensurate with all frequencies present in the forcing.

Figure 12. Comparison of linear and nonlinear time integrations for $Re=600$, numerical resolution R1. (a,b) Perturbation azimuthal velocity $u_{\theta }$ at location ($r=5$, $z=0.9$) from linear (a) and nonlinear time integration (b). (c) Corresponding spectrum of nonlinear (solid lines) and linear (dashed lines) velocity signal at $z=0.9$ for various $r$.

Figure 13. Comparison of linear and nonlinear time integrations for $Re=2100$, numerical resolution R1. (a,b) Perturbation azimuthal velocity $u_{\theta }$ at location ($r=5$, $z=0.9$) from linear (a) and nonlinear time integration (b). (c) Corresponding spectrum of nonlinear (solid lines) and linear (dashed lines) velocity signal at $z=0.9$ for various $r$.

Our goal is now to map in a $(A,Re)$ diagram the different behaviours encountered in the simulations, whether linear or nonlinear. The followed strategy is illustrated in figure 14 by focusing on $Re=2100$. The time series corresponding to the observable $\| \omega _{pert} \|$ defined by

(6.1)\begin{equation} \| \omega_{pert} \|=\sqrt{\iint \vert\omega -\omega_{b}\vert^2 r \,{\rm d}r \,{\rm d}z}, \end{equation}

where $\omega$ is the azimuthal vorticity and $\omega _b$ that of the base flow, are shown in figure 14(a) in logarithmic scale, and in figure 14(b) in linear scale, for an amplitude $A$ spanning 6 decades from $A=10^{-12}$ to $A=10^{-6}$ (a) and 6 decades from $A=10^{-7}$ to $A=10^{-1}$ (b). For $A$ from $A=10^{-12}$ to $A=10^{-8}$, the time series in figure 14(a) differ only by their amplitude, which turns out to be proportional to the value of $A$. For larger $A=10^{-7}$ and especially $10^{-6}$, the time series change aspect and their amplitude no longer follows $A$ linearly. Note that the logarithmic scale is essential for this assessment, by comparison figure 14(b) is much less easy to interpret than figure 14(a).

Figure 14. Global observable $\| \omega _{pert} \|$ in nonlinear time integration at $Re=2100$ for varying forcing amplitude $A$. The values of $A$ are indicated both in the legend and in the plot for simplicity. Starting from $A=10^{-7}$ nonlinear effects are observed. For $A\geq 10^{-6}$ the observable jumps to a different level corresponding to the presence of a non-trivial top branch solution.

This analysis can be repeated for different values of $Re$. In figure 15(a), the time-averaged $\| \omega _{pert} \|$ is plotted directly vs the forcing amplitude $A$. This representation allows for a direct assessment of the parameter values for which $\| \omega _{pert} \|$ scales linearly with $A$, which suggests linear receptivity. All the points lying away from the related straight lines are labelled as nonlinear in the $(A,Re)$ map of figure 15(b). Beyond the simple linear/nonlinear labelling, we note that, for the three largest values of $Re=2100$, $2400$ and $2700$ in figure 15(a), the observable $\| \omega _{pert} \|$ saturates with increasing $A$ to an almost constant level. This suggests that the final state of the simulation is the same whatever $A$, in other words that the forced system is attracted towards a different region of the state space for at least $Re \gtrsim 2100$. A non-trivial chaotic attractor was reported, for the same parameters, for the unforced problem, by Daube & Le Quéré (Reference Daube and Le Quéré2002) and Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024), for $Re \gtrsim 1800$. This attractor is different from the laminar one and it has a different attraction basin in terms of initial condition. There is little doubt that the same attracting chaotic state is identified here in the presence of strong amplitude forcing (the initial condition being fixed). A closer zoom on the figure, especially in linear scale, would show that the saturated value of the observable depends weakly on $A$, which suggests that the associated attractor is gently sensitive to the unsteady forcing. When such an attractor is reached, the corresponding point $(A,Re)$ in the map of figure 15(b) changes from red to blue. The points that are red correspond to the values of $(A,Re)$ where nonlinear effects do bend the gain curves (panel a), but no finite-amplitude state is present and thus no saturation to any non-trivial state can occur. This leads to a 3-colour cartography of the parameter space $(A,Re)$ for the chosen forcing protocol parameterised by the amplitude parameter $A$. The notion of double threshold, popular in the shear flow transition community (Grossmann Reference Grossmann2000), appears here in both $A$ and $Re$ (rather than in initial perturbation amplitude and $Re$), for a forced rather than unforced problem. The boundary between linear and nonlinear behaviour in 15(b) is consistent with an approximate fit of the form $A=O(Re^{-\alpha })$, with $\alpha \approx 10$. This large value of $\alpha$ is the signature of a very steep boundary. As for the boundary between red (nonlinear) and blue (non-trivial saturation), the present data do not suggest any simple fit. Besides the dataset is known to end at $Re \approx 3000$, beyond which all points are blue because the base flow is linearly unstable and the turbulent state is the only attractor left (Gesla et al. Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024). Rather than a cartography in the $(A,Re)$ plane, where $A$ is specific to a given forcing protocol, we plot in figure 16 the time-averaged observable $\overline {\| \omega _{pert} \|}$ vs $Re$, the value of $A$ being treated implicitly. The data are represented using the same colour coding as the preceding figure, both in linear (a) and in logarithmic scales (b). In both plots the data corresponding to the non-trivial top and lower branches (i.e. edge states) obtained in Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024) are also superimposed for clarity (the same observable was used, see their figure 19). The two upper and lower branches from the unforced problem were reported to merge in a saddle-node bifurcation at $Re=Re_{SN} \approx 1800$. Interestingly, none of the red or blue data points, which correspond to the unforced problem, fall in the gap between the lower and upper branches. This suggests that the data from the deterministic problem can be used to delineate between the two types of response, even in the presence of a finite-amplitude forcing. The panel (b) of the figure shows the same data in logarithmic scale. There, for the choice of this observable, the boundary between linear and nonlinear behaviour obeys an approximate power-law scaling $\overline {\| \omega _{pert} \|}=O(Re^{-\beta })$, with $\beta \approx 3.88$.

Figure 15. (a) Mean observable value $\overline {\| \omega _{pert} \|}$ as a function of forcing amplitude $A$. The linear regime is indicated using open symbols whenever the slope of the line is at most 1 % different from 1. For $Re\geq 2100$ the mean observable value level jumps to the top-branch level for strong enough forcing amplitudes. (b) Amplitudes corresponding to linear and nonlinear regime indicated respectively with green and red symbols. The top branch is reached for $Re>1800$ whenever $\overline {\| \omega _{pert} \|}>5$.

Figure 16. (a) Mean observable value $\overline {\| \omega _{pert} \|}$ as a function of $Re$. The data corresponding to the chaotic solutions from Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024) are superimposed in respectively blue (top branch) and pink (for the edge branch). Same colour coding as in the previous figure. (b) Same data plotted on logarithmic scales.

6.2. Leaky chaotic attractors and finite lifetime dynamics

From the previous subsection, we know that forcing the nonlinear problem with a sufficiently large amplitude leads, for $Re$ above $Re_{SN} \approx 1800$, to a response different from the response triggered at lower forcing amplitudes. The difference is attested for global observables such as in figure 14 and it is also clear from local velocity probes, as shown in figure 17 for $Re=2100>Re_{SN}$. In panel (a), low amplitude values $A \le 10^{-7}$ (corresponding to green points in figure 15) lead to a disordered response whose diameter increases (linearly) with the value of $A$. In panel (b), only values of $A \ge 10^{-6}$ are displayed and it is apparent again, especially from the numbers along the axes, that the response depends weakly on $A$. This last point suggests the existence of an attracting state independent of $A$, in other words, a deterministic attractor, solution of the forced problem, identical to the top-branch solutions identified by Gesla et al. (Reference Gesla, Duguet, Le Quéré and Martin Witkowski2024). Similar conclusions can be drawn for all $Re$ larger than 2100, including the supercritical values larger than $Re_c\approx 3000$ when only the non-trivial attractor is stable.

Figure 17. Two-dimensional state portrait $(u_r,u_z)$ from a velocity probe located $(r=5, z=0.9)$ at $Re=2100$. The width of the state portraits increases linearly with the forcing amplitude for $A\leq 10^{-7}$. For $A\geq 10^{-6}$ state portraits lose their elliptic shape. The centre of mass of the state portrait also shifts away from 0, indicating strong nonlinear interaction.

Departures from this picture can, however, be noted at lower $Re$ close to (but above) $Re_{SN}$. Figure 18 shows the outcome of several temporal simulations of the unforced problem for $Re=1840$, focusing on the time series of $\| \omega _{pert} \| (t)$. These time series demonstrate that the chaotic behaviour lasts for a finite time only, after which the observable rapidly reaches zero, indicating a return to the base flow. Some of the measured lifetimes are much longer than the duration of the return to the base flow, which justifies the label supertransient (Lai & Tél Reference Lai and Tél2011). The statistics of the lifetimes have been gathered over many such realisations of the deterministic problem, all initialised by different velocity fields drawn randomly. The cumulative distributions of the lifetimes $P(T>t)$, where $T$ is the lifetime, are shown in figure 19 for several values of $Re$ between 1800 and 1860. The distributions all look exponential for up to two decades, which suggests a memoryless process (Bottin & Chaté Reference Bottin and Chaté1998).

Figure 18. Time series of $\| \omega _{pert} \|(t)$ for $Re=1840$, unforced problem. Random initial perturbations of amplitude $8\times 10^{-3}$ added to the $u_{\theta }$ component. Four different initial conditions lead to four different lifetime values.

Figure 19. Cumulative distribution function $P(T>t)$ for the lifetime $T$. Here, $P$ is estimated as the fraction of runs that did not relaminarise in the time interval $(0:t)$. Relaminarisation is decided whenever $\| \omega _{pert} \|<10^{-3}$. For each value of $Re$, simulations with different random initial perturbations were conducted.

As noted by Avila, Willis & Hof (Reference Avila, Willis and Hof2010), the perturbations need a certain time $t_0$, before they reach a leaky attractor. The mean lifetime $\tau$ of only these perturbations that have reached the leaky attractor can be computed by discarding all of the simulations that relaminarised before the time $t_0$

(6.2)\begin{equation} \tau(t_0)=\langle t-t_0 \,|\, t>t_0\rangle, \end{equation}

where $\langle {\cdot }\rangle$ stands for the mean. Here, $\tau (t_0)$, which is a constant function of $t_0$, is a characteristic feature of a memoryless process.

In case of relaminarisation of every simulation out of $r$ total simulations, the maximum likelihood estimator (MLE) of $\tau$ can be computed using

(6.3)\begin{equation} \tau=\frac{1}{r}\sum_{i=1}^rt_i, \end{equation}

where $t_i$ is the relaminarisation time of individual simulations. However, because of the finite numerical integration time, not all simulations relaminarise. To account for the perturbations whose lifetime exceeds a censoring time $t_r$, the MLE is modified as (Lawless Reference Lawless2003)

(6.4)\begin{equation} \tau=\frac{1}{r}\left( \sum_{i=1}^rt_i+(n-r)t_r\right), \end{equation}

with $r$ the number of simulations that relaminarised and $n$ the total number of simulations. A 95 % confidence interval of the estimator (6.4) is

(6.5)\begin{equation} \tau \times \left[ \frac{2r}{\chi^2_{2r,0.975}},\frac{2r}{\chi^2_{2r,0.025}}\right], \end{equation}

where $\chi ^2_{m,p}$ is the $p$th quantile of the chi-squared distribution with $m$ degrees of freedom. The total number of simulations performed for each $Re$, together with $t_0$, the escape rate from the saddle $\kappa$ (Tél & Lai Reference Tél and Lai2008), its confidence interval and a corresponding $\tau$, are all reported in table 3. The dependence of $\kappa$ on $t_0$ and $Re$ is shown in figure 20. Linear dependence $\log (\kappa (Re))$ would suggest an exponential scaling of the mean lifetime as a function of $Re$, at least over the interval in $Re$ where the data were gathered (which lies well below the linear instability threshold).

Table 3. Statistics of lifetimes. For each value of $Re$, $r$ is the number of simulations out of total $n$ to relaminarise before the censoring time $t_r=7800$. The distribution of lifetimes is approximately exponential for $t>t_0$. The MLE of a mean lifetimes is computed using (6.4) and its confidence interval using (6.5). The escape rate from the saddle $\kappa =1/\tau$ and its confidence interval are plotted in figure 20.

Figure 20. (a) Escape rate $\kappa (t_0)$ for $Re=1800$ and $Re=1810$. The exponential distribution is reached when $\kappa$ is independent of $t_0$, which can be assumed for $t_0>800$ when $Re=1800$, and for $t_0>850$ when $Re=1810$. (b) Dependence of $\kappa$ on $Re$. The fits in the figure are $\exp (-6.97\times 10^{-2}Re+ 1.2\times 10^{2})$ (red) and $\exp (-\exp (9.13\times 10^{-3}Re -1.5\times 10^{1}))$ (blue), suggesting super-exponential dependence on $Re$.

The whole situation is reminiscent from the finite probabilities to relaminarise in 3-D turbulence in wall-bounded shear flows, notably the case pipe flow, where lifetimes were reported to increase super-exponentially with $Re$ (Hof et al. Reference Hof, De Lozar, Kuik and Westerweel2008). For the case of pipe flow, the issue of whether lifetimes truly diverge at a finite value of $Re$ or stay finite was solved by remarking that another competing phenomenon, namely turbulence proliferation, would outweigh turbulence decay above some threshold in $Re$ (Avila et al. Reference Avila, Moxey, De Lozar, Avila, Barkley and Hof2011) and keep turbulence alive for all times. There is no such equivalent here in axisymmetric rotor-stator flow, for which a scenario where lifetimes are always finite is a possibility below the threshold of linear stability of the base flow $Re_c\approx 3000$ (Daube & Le Quéré Reference Daube and Le Quéré2002). Nevertheless, the values of $T$ found in practice are so large that they usually exceed the observation times allowed for. In other words, the hypothesis that the invariant set underlying the existence of non-trivial dynamics in the forced problem is a deterministic attractor is only here for commodity: mathematically speaking, it could possibly be a leaky attractor i.e. a chaotic saddle (Lai & Tél Reference Lai and Tél2011).