2.1. Direct numerical simulations in the annular–parallelogram domain

The nonlinear boundary value problem ((2.8)–(2.9)) is discretised using a solenoidal Petrov–Galerkin scheme formerly formulated by Meseguer et al. (Reference Meseguer, Avila, Mellibovsky and Marques2007), and suitably adapted to the annular–parallelogram domain (2.15). In the transformed domain, the solenoidal velocity perturbation is approximated by means of a Fourier $\times$ Fourier $\times$ Chebyshev spectral expansion ${\boldsymbol u}_{s}$ of order $N \times L \times M$ in $\xi \times \zeta \times r$, respectively, of the form

(2.19)\begin{equation} {\boldsymbol u}_{s}(r,\xi,\zeta;t)=\sum_{\ell={-}L}^{L} \sum_{n={-}N}^{N} \sum_{m=0}^{M} a_{\ell nm}^{(1)}(t) \boldsymbol{\varPhi}_{\ell nm}^{(1)}(r,\xi,\zeta)+a_{\ell nm}^{(2)}(t) \boldsymbol{\varPhi}_{\ell nm}^{(2)}(r,\xi,\zeta). \end{equation}

Our aim here is to derive the dynamical system satisfied by the coefficients $a_{\ell nm}^{(\imath )}(t)$, as symbolically represented by the $2\times (M+1)\times (2L+1)\times (2N+1)$-dimensional state vector $\boldsymbol {a}(t)$. The binary superindex $\imath =\{1,2\}$ and the factor 2 in the count of unknowns follow from the two degrees of freedom per grid point that remain after taking into consideration that the three velocity components are not independent, but linked by the solenoidal condition. The vector fields $\boldsymbol {\varPhi }_{\ell n m}^{(\imath )}$ constitute the elements of the trial basis of solenoidal vector fields of the form

(2.20)\begin{equation} \boldsymbol{\varPhi}_{\ell n m}^{(\imath)}(r,\xi,\zeta) = \mathrm{e}^{\mathrm{i}(n\xi+\ell\zeta)} \boldsymbol{u}_{\ell nm}^{(\imath)}(r) = \mathrm{e}^{\mathrm{i}(n\xi+\ell\zeta)} (u_{\ell n m}^{(\imath)},v_{\ell n m}^{(\imath)},w_{\ell n m}^{(\imath)}), \end{equation}

where $u_{\ell n m}^{(\imath )}$, $v_{\ell n m}^{(\imath )}$ and $w_{\ell n m}^{(\imath )}$ are the radial, azimuthal and axial components of $\boldsymbol {u}_{\ell nm}^{(\imath )}(r)$, respectively. Each element of the trial basis satisfies the divergence-free condition (2.9) that, in the $(r,\xi,\zeta )$ variables, explicitly reads

(2.21)\begin{equation} \left(\partial_r + \frac{1}{r}\right) u_{\ell nm}^{(\imath)} + \frac{\mathrm{i}}{r}\left(n_1n + n_2\ell \right)v_{\ell n m}^{(\imath)} + {\rm i}\left(k_1n+k_2\ell \right) w_{\ell nm}^{(\imath)}=0. \end{equation}

Since $\boldsymbol {u}$ represents the perturbation of the velocity field, it must therefore vanish at the inner ($r=r_{i}$) and outer ($r=r_{o}$) walls of the cylinders. Therefore, $\boldsymbol {\varPhi }_{\ell n m}^{(\imath )}$ must also satisfy the homogeneous boundary conditions

(2.22)\begin{equation} \boldsymbol{\varPhi}_{\ell n m}^{(\imath)}(r_{i},\xi,\zeta;t)=\boldsymbol{\varPhi}_{\ell n m}^{(\imath)}(r_{o},\xi,\zeta;t)=\boldsymbol{0}. \end{equation}

In what follows, we define the transformed radial coordinate

(2.23)\begin{equation} x(r) = 2r-\frac{1+\eta}{1-\eta}, \end{equation}

that maps the radial domain $r\in [r_{i},r_{o}]$ to the interval $x\in [-1,1]$. In addition, we define the radial functions

(2.24a,b)\begin{equation} {h}_m(r)=(1-x^2){T}_m(x),\quad {g}_m(r)=(1-x^2)^2{T}_m(x), \end{equation}

where ${T}_m(x)$ is the Chebyshev polynomial of degree $m$. Finally, we introduce the Chebyshev weight function ${w}(x)=(1-x^2)^{-1/2}$, defined over the interval $(-1,1)$. The functions introduced in (2.24a,b) satisfy

(2.25)\begin{equation} {h}_m(r_{i})={h}_m(r_{o})={g}_m(r_{i})={g}_m(r_{o})={\rm D}{g}_m(r_{i})={\rm D}{g}_m(r_{o})=0, \end{equation}

where ${\rm D}$ stands for the radial derivative ${\rm d}/{\rm d}r$. The solenoidal spectral method consists in devising complete sets of vector fields (trial functions) satisfying (2.21) and (2.22). For $n n_1+\ell n_2 = 0$, two such vector fields are

(2.26a,b)\begin{equation} \boldsymbol{u}_{\ell n m}^{(1)}(r)=(0,{h}_m,0),\quad \boldsymbol{u}_{\ell n m}^{(2)}(r)=( -\mathrm{i}(n k_1+\ell k_2) r {g}_m,0,{\rm D}[r {g}_m]+{g}_m), \end{equation}

with the third component of $\boldsymbol {u}_{\ell n m}^{(2)}$ replaced by ${h}_m$ whenever $n k_1+\ell k_2=0$. Finally, for ${n n_1+\ell n_2 \neq 0}$, the solenoidal basis is

(2.27)$$\begin{gather} \boldsymbol{u}_{\ell nm}^{(1)}(r) =(-\mathrm{i}(n n_1+\ell n_2) {g}_m , {\rm D}[r {g}_m] , 0), \end{gather}$$

(2.28)$$\begin{gather}\boldsymbol{u}_{\ell nm}^{(2)}(r) =(0,-\mathrm{i} (n k_1+\ell k_2) r{h}_m , \mathrm{i}(n n_1+\ell n_2) {h}_m), \end{gather}$$

except that the third component of $\boldsymbol {u}_{\ell nm}^{(2)}$ is replaced by ${h}_m$ when $n k_1+\ell k_2=0$. The Petrov–Galerkin solenoidal weak formulation is completed by introducing the Hermitian product of two arbitrary solenoidal trial and dual fields $\boldsymbol {\varPhi }$ and $\boldsymbol {\varPsi }$, respectively, over the annular–parallelogram domain (2.15)

(2.29)\begin{equation} (\boldsymbol{\varPsi},\boldsymbol{\varPhi})=\int_{0}^{2{\rm \pi}}\int_{0}^{2{\rm \pi}} \int_{r_{i}}^{r_{o}} \; \boldsymbol{\varPsi}^{\dagger} \boldsymbol{\cdot} \boldsymbol{\varPhi} \;r\,\mathrm{d}r\,\mathrm{d}\xi\,\mathrm{d}\zeta. \end{equation}

Accordingly, we consider the dual basis for the projection space. In particular, the basis for the case $n n_1+\ell n_2 = 0$ is

(2.30)$$\begin{gather} \boldsymbol{\tilde{u}}_{\ell 0m}^{(1)}(r) = {w}(0, r {h}_m, 0), \end{gather}$$

(2.31)$$\begin{gather}\boldsymbol{\tilde{u}}_{\ell 0m}^{(2)}(r) = r^{{-}2}{w} (\mathrm{i}(n k_1+\ell k_2) {g}_m,0, {\rm D}_+{g}_m +2r^{{-}1}(1-x^2+rx){h}_m), \end{gather}$$

where ${\rm D}_+={\rm D}+r^{-1}$, and the third component of $\boldsymbol {\tilde {u}}_{\ell 0m}^{(2)}$ is replaced by $r {h}_m$ if $n k_1+\ell k_2=0$. Similarly, the basis for the case $n n_1+\ell n_2 \neq 0$ is

(2.32)$$\begin{gather} \boldsymbol{\tilde{u}}_{\ell n m}^{(1)}(r) = {w} ((n n_1+\ell n_2)\,r{g}_m , r{\rm D}_+[r {g}_m]+2xr^2{h}_m ,0), \end{gather}$$

(2.33)$$\begin{gather}\boldsymbol{\tilde{u}}_{\ell n m}^{(2)}(r) = {w} (0,i (n k_1+\ell k_2) r^2{h}_m,-i (n n_1+\ell n_2)\,r{h}_m). \end{gather}$$

These projection basis elements contain the Chebyshev weight function ${w}(x)$ so that the resulting radial integration involved in (2.29) can be computed exactly by suitable quadrature formulae (Moser, Moin & Leonard Reference Moser, Moin and Leonard1983; Meseguer et al. Reference Meseguer, Avila, Mellibovsky and Marques2007; Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2007, Reference Canuto, Hussaini, Quarteroni and Zang2010; Meseguer Reference Meseguer2020).

Formal substitution of the spectral expansion (2.19) into (2.8), followed by Hermitian projection onto each one of the dual basis elements leads to

(2.34)\begin{equation} (\boldsymbol{\varPsi}_{\ell n m}^{(\imath)} , \partial_t{\boldsymbol u}_{S}) = (\boldsymbol{\varPsi}_{\ell n m}^{(\imath)} , \nabla^2{\boldsymbol u}-({\boldsymbol v}_{b} \boldsymbol{\cdot} \boldsymbol{\nabla}){\boldsymbol u}-({\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}){\boldsymbol v}_{b}-({\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}){\boldsymbol u} + f\,\hat{{\boldsymbol{z}}}),\end{equation}

for $\ell = -L,\ldots,L,\, n=-N,\ldots,N$, $m = 0,\ldots, M$ and $\imath =1,2$. The pressure deviation field drops upon projection, $(\boldsymbol {\varPsi }_{\ell n m}^{(\imath )},\boldsymbol {\nabla } q)=0$, by virtue of Stokes’ theorem and has therefore been omitted. Equation (2.34), subject to the constraint (2.3), constitutes a dynamical system

(2.35)$$\begin{gather} {\mathsf{A}}_{\ell nm(\jmath)}^{pqr(\imath)} \, \dot{a}_{pqr}^{(\jmath)} = {\mathsf{B}}_{\ell n m(\jmath)}^{pqr(\imath)} \, a_{pqr}^{(\jmath)} + \left[\boldsymbol{N}(\boldsymbol{a})\right]_{\ell n m}^{(\imath)} + f{F}_{\ell nm}^{(\imath)}, \end{gather}$$

(2.36)$$\begin{gather}Q^r a_{00r}^{(2)} = 0, \end{gather}$$

for the axial forcing $f(t)$ and amplitudes $a_{\ell n m}^{(\imath )}(t)$, where repeated indices must be interpreted following the index summation convention. The zero-net-massflux constraint reduces to a mere linear equation for the $n=l=0$, $\imath =2$ coefficients upon substitution of (2.19). The quadratic form $\left[\boldsymbol{N}(\boldsymbol{a})\right]_{\ell n m}^{(\imath)}$ appearing in (2.35) corresponds to the projection of the nonlinear convective term, $(\boldsymbol {\varPsi }_{\ell n m}^{(\imath )} , ({\boldsymbol u} \boldsymbol {\cdot } \boldsymbol {\nabla }){\boldsymbol u})$, which is computed pseudospectrally using Orszag's $3/2$-dealiasing rule (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2007, Reference Canuto, Hussaini, Quarteroni and Zang2010). Overall, the resulting stiff system of ordinary differential equations is integrated in time by means of a fourth-order linearly implicit backwards differentiation scheme with explicit polynomial extrapolation of the nonlinear terms, conveniently started with a fourth-order Runge–Kutta method.

2.2. Computation and stability analysis of invariant solutions

In the $(\xi,\zeta )$ coordinate system, a travelling wave is represented by the spectral expansion

(2.37)\begin{align} &{\boldsymbol u}_{s}(r,\xi,\zeta;t)=\nonumber\\ &\qquad =\displaystyle{\sum_{\ell={-}L}^{L} \sum_{n={-}N}^{N} \sum_{m=0}^{M} [\check{a}_{\ell nm}^{(1)} \boldsymbol{\varPhi}_{\ell nm}^{(1)}(r,\xi,\zeta) + \check{a}_{\ell nm}^{(2)} \boldsymbol{\varPhi}_{\ell nm}^{(2)}(r,\xi,\zeta)] \mathrm{e}^{\mathrm{i} n(\xi-c_\xi t)}\mathrm{e}^{\mathrm{i} \ell(\zeta-c_\zeta t)}}, \end{align}

where $c_\xi$ and $c_\zeta$ are the unknown wave speed components in the $\xi$ and $\zeta$ directions, respectively, so that the Fourier–Chebyshev spectral coefficients of this particular type of solution read

(2.38)\begin{equation} a_{\ell nm}^{(\imath)}(t) = \check{a}_{\ell nm}^{(\imath)} \mathrm{e}^{-\mathrm{i} n c_\xi t}\mathrm{e}^{- \mathrm{i} \ell c_\zeta t}, \end{equation}

with the complex constant $\check {a}_{\ell nm}^{(\imath )}$ representing the wave shape unambiguously except for arbitrary rotations and shifts. In this case, formal introduction of expansion (2.37) in (2.8) followed by Hermitian projection leads to the following system of nonlinear algebraic equations for the travelling wave coefficients $\check {a}_{\ell nm}^{(\imath )}$ ($\check {\boldsymbol {a}}$ for brevity):

(2.39)$$\begin{gather} \left[{\mathsf{B}}_{\ell n m(\jmath)}^{pqr(\imath)} + \mathrm{i} (n c_\xi + \ell c_\zeta) {\mathsf{A}}_{\ell n m (\jmath)}^{pqr(\imath)}\right] \check{a}_{pqr}^{(\jmath)} + \left[\boldsymbol{N}(\boldsymbol{\check{a}})\right]_{\ell n m}^{(\imath)} + \check{f}{F}_{\ell n m}^{(\imath)} = 0, \end{gather}$$

(2.40)$$\begin{gather}Q^r \check{a}_{00r}^{(2)} = 0, \end{gather}$$

where $\check {\boldsymbol {a}}$ appearing in the nonlinear term is the state vector representing the Fourier–Chebyshev coefficients $\check {a}_{\ell nm}^{(\imath )}$ of the travelling wave, and $\check {f}$ is the axial pressure gradient required to enforce the zero-net-flux condition. The azimuthal and axial degeneracy of solutions, associated with drift speeds $c_\xi$ and $c_\zeta$, is removed using two additional phase constraints in the same way as is done for the computation of rotating–travelling waves in pipe flow (Mellibovsky & Eckhardt Reference Mellibovsky and Eckhardt2011). The system is solved using a matrix-free Newton–Krylov method (Kelley Reference Kelley2003), implicitly using the generalized minimal residual method (GMRES) as a matrix-free solver (Trefethen & Bau Reference Trefethen and Bau1997). The converged nonlinear solutions are tracked in parameter space using pseudoarclength continuation schemes (Kuznetsov Reference Kuznetsov2004).

The linear stability of a travelling-wave solution of (2.39) is formulated by adding disturbances of very small amplitude $|\varepsilon _{\ell n m}^{(\imath )}|\ll |\check {a}_{\ell nm}^{(\imath )}|$ to its Fourier–Chebyshev coefficients following

(2.41)$$\begin{gather} a_{\ell nm}^{(\imath)}(t) = (\check{a}_{\ell nm}^{(\imath)}+\varepsilon_{\ell n m}^{(\imath)} \mathrm{e}^{\sigma t})\mathrm{e}^{-\mathrm{i} n c_\xi t}\mathrm{e}^{- \mathrm{i} \ell c_\zeta t}, \end{gather}$$

(2.42)$$\begin{gather}f(t) = \check{f} + \phi \mathrm{e}^{\sigma t}, \end{gather}$$

and the forcing perturbation $|\phi |\ll |\check {f}|$ is such that coefficient perturbations comply also with the zero-massflux condition $Q^r \varepsilon _{00r}^{(2)}=0$. Formal substitution of the perturbed solution (2.41) in (2.35), and subsequent neglect of quadratic perturbation terms, leads to the constrained generalised eigenvalue problem

(2.43)$$\begin{gather} \left( \sigma -\mathrm{i} n c_\xi - \mathrm{i} \ell c_\zeta \right){\mathsf{A}}_{\ell n m(\jmath)}^{p q r(\imath)}\varepsilon_{pqr}^{(\jmath)} = {\mathsf{B}}_{\ell n m(\jmath)}^{p q r(\imath)}\varepsilon_{pqr}^{(\jmath)} + \left[\boldsymbol{D}_{\boldsymbol{a}}\boldsymbol{N}(\check{\boldsymbol{a}})\right]_{\ell n m(\jmath)}^{p q r(\imath)}\!\varepsilon_{pqr}^{(\jmath)} + \phi {F}_{\ell n m}^{(\imath)}, \end{gather}$$

(2.44)$$\begin{gather}Q^r \varepsilon_{00r}^{(2)} = 0, \end{gather}$$

where $\boldsymbol{D}_{\boldsymbol{a}}\boldsymbol{N}(\check {\boldsymbol {a}})$ is the linear action of the Jacobian of $\boldsymbol{N}$ evaluated at the travelling-wave solution $\check {\boldsymbol {a}}$, and $(p,q,r,\jmath ) \in [0,L]\times [-N,N]\times [0,M]\times \{1,2\}$. This linear action therefore allows formulation of the generalised eigenvalue problem (2.43) in a matrix-free form, so that Arnoldi methods (Trefethen & Bau Reference Trefethen and Bau1997) can be applied to compute the leading eigenvalues $\sigma _j$ and their associated eigenvectors $\boldsymbol {\varepsilon }_j=\{\varepsilon _{\ell m n}^{(\imath )}\}_j$.

To compute relative periodic orbits (i.e. modulated travelling waves) beyond their region of linear stability, a Poincaré–Newton–Krylov method is devised. The method is essentially an adaptation of the one used for the computation of modulated travelling waves in plane Poiseuille flow (Mellibovsky & Meseguer Reference Mellibovsky and Meseguer2015). In this case, the method solves the nonlinear system of equations resulting from root finding for the map defined by consecutive crossings of a Poincaré section $\boldsymbol{P}$, as follows:

(2.45)\begin{equation} \boldsymbol{a} \rightarrow \tilde{\boldsymbol{a}} = \boldsymbol{P}(\boldsymbol{a})= \boldsymbol{\mathsf{T}}(\Delta \xi,\Delta \zeta) \boldsymbol{\varphi}(\boldsymbol{a};{T}), \end{equation}

where $\boldsymbol {\varphi }(\boldsymbol {\cdot };t)$ is the action of the uniparametric group or flow generated by (2.35), ${T}$ is the modulation period of the relative periodic orbit, and

(2.46)\begin{equation} \left[\boldsymbol{\mathsf{T}}(\Delta \xi,\Delta \zeta)\boldsymbol{a}\right]_{\ell n m}^{(\imath)} = \mathrm{e}^{-\mathrm{i} n \Delta \xi }\mathrm{e}^{-\mathrm{i} \ell \Delta \zeta}a_{\ell n m}^{(\imath)} \end{equation}

is a double-shift operator, removing the drift of the relative periodic orbit in the two homogeneous parallelogram coordinates $\xi$ and $\zeta$. Finally, the above-described time-stepping formulation, as well as the travelling wave and relative periodic orbit solvers, are enforced to satisfy zero net flux condition for the perturbation field.

5.1. Flow structure of DRW

Let us begin with a detailed analysis of the flow structure and properties of the several coexisting DRW solutions in the $k_2=4.5$ domain at subcritical ${{R}_{i}}=450$. Figure 6 provides the three-dimensional characterisation of the flow structure of DRW at the representative points labelled A, B and C in figure 3. Of the solutions ensuing, saddle node SN$_1$, point C is chosen to represent the lower-torque branch. Two isosurfaces of the azimuthal perturbation velocity ($v$) shown in figure 6(ci) reveal the presence of a low-speed streak close to the inner cylinder. For the higher-torque branch (point B, figure 6bi), this slow streak induces a strong velocity distortion that reaches the middle of the gap, an effect that is all the more pronounced for the solutions originated at saddle-node SN$_3$. For the solutions related to SN$_3$, only point A (figure 6ai) on the higher-torque branch is shown, as the corresponding lower-torque solution is very similar. Another characteristic feature of the flow field is the presence of a vortex sheet, which is visualised in figure 6(ii) through a couple of azimuthal vorticity ($\omega _{\theta }$) isosurfaces. The position of the sheet does not change much in the azimuthal direction, but the sign of the vorticity fluctuates in a sinuous fashion. It is precisely the interaction of these two fields, the fluctuating vortex layer and the quasiazimuthally invariant streak, that might be held responsible for the self-sustainment of DRW in the absence of a linear instability of CCF. The physical reasons for the requirement of a feedback mechanism from the wave to the streak field can be ascertained by examining the centrifugally unstable region adjacent to the inner cylinder, bounded by $r_{i}< r< r_{n}$. Although the Rayleigh criterion predicts centrifugal inviscid instability of CCF in this region, the laminar base flow remains stable because the effective Taylor number is still slightly short of the critical threshold (Esser & Grossmann Reference Esser and Grossmann1996; Deguchi Reference Deguchi2016). The roll–streak field of DRW must therefore be supported by some other mechanism, which, as we will argue, is none other than the SSP. Following the usual definitions for parallel flows (Waleffe Reference Waleffe1997), the total velocity field (base flow plus perturbation) may be additively decomposed into a streamwise-averaged field (the roll–streak), and the rest (the wave). The roll–streak is then further split into its streamwise (streak) and cross-stream (roll) velocity components, respectively. In TCF, the streak and roll fields are the toroidal (azimuthal velocity) and poloidal (radial-axial velocity) components, respectively, of the axisymmetric part of the total velocity field. Figure 7 illustrates the roll–streak–wave decomposition of DRW solutions A, B and C. The streak and roll components are naturally represented through their angular azimuthal velocity ($\langle V \rangle _{\theta }/r$) and azimuthal vorticity ($\langle \omega _{\theta }\rangle _{\theta }$) field colourmaps, respectively, on any arbitrary $r$–$z$ cross-section (see figure 7i,ii). For the wave component, azimuthal vorticity ($\omega _{\theta }^{w}=\omega _{\theta }-\langle \omega _{\theta }\rangle _{\theta }$) has been employed in figure 7(iii,iv), corresponding to two $r$–$z$ cross-sections at $\theta =\{0,{\rm \pi} /n_1\}$, exactly half an azimuthal wavelength apart.

Figure 6. Three-dimensional flow structure of DRW solutions at subcritical ${{R}_{i}}=450$, corresponding to points (a) A, (b) B and (c) C of figure 3. Shown are positive (yellow) and negative (blue) isosurfaces of (i) perturbation azimuthal velocity $v$, and (ii) azimuthal vorticity $\omega _\theta$. The corresponding isosurface levels are $v=\{-100,250\}$ and $\omega _\theta =\pm 1000$ (solution A), $v=\{-100,250\}$ and $\omega _\theta =\pm 600$ (B), and $v=\{-40,90\}$ and $\omega _\theta =\pm 300$ (C).

Figure 7. Comparison of (i) the streak field as signified by azimuthally averaged distribution of total angular azimuthal velocity $\langle V \rangle _{\theta }/r\in [-1200,450]$, (ii) the roll field, exposed through the azimuthally averaged azimuthal vorticity field $\langle \omega _{\theta } \rangle _{\theta }\in [-1000,1000]$ and (iii,iv) the wave field, represented by azimuthal vorticity $\omega _{\theta }^{w}=\langle \omega _{\theta } \rangle _{\theta }-\omega _{\theta }\in [-1000,1000]$, at $\theta =0$ (iii) and $\theta ={\rm \pi} /n_1=0.314$ (iv), for DRW solutions labelled (a) A, (b) B and (c) C in figure 3. Lines indicate the location of the critical layer (solid) and the nodal radius $r_{n}$ (dashed).

In the case of parallel shear flows, the rolls are energised by the Reynolds stress field of the waves, the streaks are generated by the lift-up effect from the rolls and the waves are driven by the instability of the streaks. However, in TCF the interaction between rolls and streaks is two-way, because they are coupled by the Coriolis force term. To prove that the SSP is at work, we must therefore verify that the roll–streak component triggers the wave component and is, at the same time, regenerated by it.

A numerical experiment has been devised to substantiate the essential part that the wave component plays in the sustainment of the roll–streak components. Starting a simulation from the DRW solution at B with the wave-component (all modes with $n \neq 0$) removed results in the rapid and monotonic decay of the roll–streak field, as clear from figure 8(a) (black dashed line). Using the full flow field as the initial condition, has the roll–streak components endure for a considerable lapse of time (blue line).

Figure 8. The role of the wave component in sustaining DRW at B of figures 3, 6(b) and 7(b). (a) Time evolution of kinetic energy $\kappa$ starting from either the DRW (solid blue line) or only its roll–streak components (dashed black). (b) Azimuthal vorticity colourmaps at $\theta =0$ (bi) and ${\rm \pi} /n_1=0.314$ (bii) of the only unstable eigenmode of the roll–streak system.

The origin of the wavy field in DRW can in turn be explained by a linear instability of the roll–streak component. The stability analysis method of (2.43) can be applied on any desired background flow field, not necessarily an exact solution of (2.39). In doing so, one is simply assuming that an appropriate forcing term is added to (2.39) that turns this flow field into an exact solution. As expected, the roll–streak component of DRW has a pair of unstable complex conjugate eigenvalues. At B, the imaginary part of the eigenvalue portends the phase speed of the unstable wave as $\tilde {c}_{\theta }\approx 26.1$, which is close to that of the actual solution ($c_{\theta }=28.1$). Furthermore, the flow field of the unstable eigenmode, shown in figure 8(b) through azimuthal vorticity colourmaps at $\theta =0$ (figure 8bi) and ${\rm \pi} /n_1$ (figure 8bii), closely resembles the wave component of DRW at B (figures 7biii,biv). This provides solid evidence that the SSP is in place.

The solid curves shown in figures 7(iii,iv) and 8(b) mark the critical layer, which will be defined shortly. The observation of the flow around this layer contributes yet another means of endorsing the SSP hypothesis. The critical layer is formally a singularity of the inviscid linear stability problem and the high Reynolds number asymptotic theory shows that the amplitude of the wave must surge around it. In parallel flows, the nature of this wave amplitude growth in ECS is well known (Wang et al. Reference Wang, Gibson and Waleffe2007; Hall & Sherwin Reference Hall and Sherwin2010; Deguchi & Hall Reference Deguchi and Hall2014). The critical layer is the place at which the streak velocity coincides with the phase speed of the wave. There the advection term becomes small and tends to vanish, such that the inviscid approximation breaks (Lin Reference Lin1955; Drazin & Reid Reference Drazin and Reid1981).

In defining critical layers in the Taylor–Couette problem, we note that subtle differences arise from the parallel flow cases. We first note that the velocity field of DRW solutions can be rewritten as a function of $r$ and two phase variables $\theta -c_{\theta } t$ and $z-c_z t$, where the constants $c_{\theta }$ and $c_z$ can be computed easily from $c_{\xi }$ and $c_{\zeta }$. The axial phase velocity $c_z$ is small and hence the critical layer is approximately determined by comparing the azimuthal phase velocity $c_{\theta }$ alone with the streak field. As it happens, the spanwise phase velocity of drifting waves can be shown to vanish in the asymptotic limit of increasing Reynolds number for parallel shear flows, and the same may presumably be expected when the geometry is curved. Since $c_{\theta }$ is an angular velocity, the angular velocity $\langle V\rangle _{\theta }/r$ of the streak must be used in defining the critical layer as the location where $\langle V\rangle _{\theta }/r=c_{\theta }$. In the streak representation of figure 7(i), the critical layer naturally coincides with one of the angular velocity contours. The azimuthal cross-sections of $\omega _{\theta }^{w}$ colourmaps at $\theta =\{0,{\rm \pi} /n_1\}$ (figure 7iii,iv) unequivocally show that the vortex sheet amplitude is strongest precisely around the critical layer. The same holds for the unstable eigenmode of the roll–streak component in figure 8(b), furnishing compelling evidence that the instability is of an inviscid nature.

The magnitude of $\omega _{\theta }^{w}$ increases as the critical layer shifts radially outwards (see figure 7). According to Rayleigh's stability condition, the closer to the outer cylinder one looks, the more centrifugally stable the flow is locally, and a stronger local wave amplitude is therefore required to sustain the streak field. Rayleigh's condition is based on the CCF profile, so that weak Taylor vortices may still arise due to nonlinearly driven mean-flow-field distortions in the vicinity of the inner cylinder. This is the case of solution C (see figure 7cii) and provides a physical explanation as to why the connections between DTVF and DRW occur in the way discussed in § 4.

The SSP revealed here is fundamentally different from what happens with WVF (Dessup et al. Reference Dessup, Tuckerman, Wesfreid, Barkley and Willis2018). In WVF, energy exchanges occur among the roll, streak and wave components and, as noted earlier by Jones (Reference Jones1985) and Martinand et al. (Reference Martinand, Serre and Lueptow2014), the wave indeed arises from a linear instability of the streak. However, the roll feeds essentially on the centrifugal instability of the base flow, and the feedback effect from the waves is only second order (Dessup et al. Reference Dessup, Tuckerman, Wesfreid, Barkley and Willis2018). The azimuthal average of WVF is barely different from the TVF whose instability triggers it, and which evidently persists even in the absence of waves. While, as we have shown, in the subcritical parameter regions of counter-rotating TCF, the roll–streak system of DRW solutions cannot be sustained without waves. Moreover, even in the supercritical regions, the roll–streak fields look different with and without SSP.

5.2. Stability of DRW and the onset of chaotic dynamics

Let us now turn our attention to the stability of DRW solutions and their role in engendering a chaotic set. The count and type of eigenvalues, as computed through linear stability analysis, are indicated in brackets for a number of sampled solutions along the various DRW branches in figure 3. From saddle-node SN$_1$ onwards, the higher-torque branch is stable while the lower-torque solutions have a single real unstable eigenvalue. This result, which is a must for one-dimensional dynamical systems exhibiting a saddle-node bifurcation, is nevertheless rarely encountered in high-dimensional systems such as TCF. In fact, the same (or closely related) DRW solutions present a higher number of unstable eigenvalues when considered in the usual orthogonal domain (Deguchi et al. Reference Deguchi, Meseguer and Mellibovsky2014) or in parallelograms of size $k_2$ other than 4.5. Our choice of domain is therefore expressly convenient for the simplest possible exploration of the onset of chaotic dynamics.

The stability of the higher-torque DRW branch emerged from SN$_1$ is lost in a Hopf bifurcation at ${{R}_{i}}^{H} = 392.85$ (point H in figure 3), as implied by the crossing of a pair of complex eigenvalues into the positive-real-half of the complex plane. The resulting relative periodic orbit has a structure similar to DRW but is subject to time-periodic amplitude oscillations (P-DRW). The Poincaré–Newton–Krylov method described in § 2 has been deployed to perform natural continuation of the P-DRW solution branch. Figure 9(a) signifies the oscillation amplitude of the P-DRW as the area (shaded blue region) bounded by the maximum and minimum inner torque $\tau _{i}$ along a full cycle (black dots). The solution amplitude obeys a square root law of the form $A_{\tau _{i}}=\tau _{i}^{max}-\tau _{i}^{min}\sim ({{R}_{i}}-{{R}_{i}}^{H})^{1/2}$, as revealed by the least-squares fit (grey dashed line) to a few of the closest points to the Hopf bifurcation (H). This attests to the supercritical nature of the Hopf bifurcation. As a result, solutions are stable at onset and could therefore have been computed by mere time stepping. The evolution of P-DRW with ${{R}_{i}}$ is more clearly illustrated by the outer-versus-inner torque $(\tau _{o},\tau _{i})$ phase-map projections of figure 9(b). With each increase in ${{R}_{i}}$, the lower-torque (crosses) and higher-torque (circles) DRW solutions become farther apart in phase space. A small limit cycle (P-DRW) clearly orbits the higher-toque DRW at ${{R}_{i}}=393$, and its amplitude grows as ${{R}_{i}}$ is further increased.

Figure 9. Onset of P-DRW solutions. (a) The bifurcation diagram near the saddle-node SN${_1}$ shown in figure 3. Line styles denote DRW branches with different stability properties: stable (solid), one unstable real eigenvalue (dashed) and one unstable complex pair (dash–dotted). The small black circles indicate the maximum and minimum $\tau _{i}$ attained by the relative periodic orbit P-DRW, whose oscillation amplitude is delimited by the shaded blue region. The P-DRW emerges supercritically at the Hopf bifurcation point H, such that the amplitudes obey locally a square root fit (grey dashed curve). (b) The $(\tau _{o},\tau _{i})$-phase map projections of upper and lower branches of DRW (bullets and crosses, respectively) and the P-DRW limit cycles (solid curves) for different values of ${{R}_{i}}$.

In order to establish the dynamical connections among the various solutions, DNS from small perturbations to the lower-torque DRW solution have been run. The initial condition has been set by scaling the exact solution following $(1+\gamma ) \boldsymbol {a}^{DRW}$, with $|\gamma |\ll 1$. Figure 10(a) shows the $(\tau _{o},\tau _{i})$-phase map projections of a couple runs at ${{R}_{i}}=392$. As we have already noted, the lower-torque DRW$_{LT}$ solution has only a single real unstable eigenvalue, so that the diagram reflects but a close approximation to its one-dimensional unstable manifold. For $\gamma =-10^{-4}$, the flow uneventfully departs towards CCF (black line), while for $\gamma =10^{-4}$ the flow is captured by the linearly stable higher-torque DRW$_{HT}$ solution (red line). The lower-torque DRW solution is therefore acting as an edge state (Itano & Toh Reference Itano and Toh2001; Skufca et al. Reference Skufca, Yorke and Eckhardt2006), separating the basins of attraction of CCF on one side and the higher-torque DRW$_{HT}$ solution on the other. The dynamics remain qualitatively the same as illustrated here within the parameter range $391.52<{{R}_{i}}<392.85$, while the upper branch solution preserves its linear stability.

Figure 10. Phase-map projections of trajectories issued from the unstable lower-torque DRW$_{LT}$ solution on the $(\tau _{o},\tau _{i})$-plane at inner Reynolds number (a) ${{R}_{i}}=392$ and (b) ${{R}_{i}}=394$. The represented solutions are CCF (square), lower-torque (DRW$_{LT}$, cross) and higher-torque (DRW$_{HT}$, circle) DRW solutions, and P-DRW (blue line). The red/black lines depict DNS trajectories starting from initial conditions taken by scaling DRW$_{LT}$ by $(1+\gamma )$, with $\gamma =\pm 10^{-4}$, respectively.

Slightly above the Hopf bifurcation (${{R}_{i}}>{{R}_{i}}^{H}=392.85$), the higher-torque DRW$_{HT}$ solution has become unstable and the stable P-DRW has popped up into existence, as illustrated at ${{R}_{i}}=394$ by figure 10(b). The dynamics are qualitatively unaltered as regards CCF or DRW$_{LT}$, which remains an edge state, but the phase-map trajectory that previously led to DRW$_{HT}$ is now only able to transiently approach it to some extent before being repelled towards P-DRW (blue line) in a spiralling fashion. The simulation eventually converges onto P-DRW, now the local attractor on this side of phase space. The insets clarify the nature of DRW$_{HT}$, which remains a focus across the Hopf bifurcation, but switches from stable to unstable. Sufficiently close to SN$_1$ the solution is instead a node, stable within a neighbourhood of $k_2=4.5$, but unstable beyond a certain threshold. These facts put together suggest that the saddle-node and Hopf bifurcations are in fact constituent pieces of a codimension-2 Takens–Bogdanov bifurcation (Kuznetsov Reference Kuznetsov2004).

At even higher Reynolds number ${{R}_{i}} = 395.5$, P-DRW has become unstable and the new stable limit cycle makes two similar but not identical revolutions in phase space before closing on itself. The resulting period-doubled P$_2$-DRW solution (black line) is shown, along with the unstable P-DRW limit cycle, at the same ${{R}_{i}}=395.5$ (dashed red) in figure 11(a). The period doubling that occurs just short of ${{R}_{i}}=395.5$ is followed by a number of ensuing period-doublings upon further increasing ${{R}_{i}}$ that eventually lead to the chaotically modulated DRW (CH-DRW) of figure 11(b) at ${{R}_{i}}=395.7816$. The transition route to chaos observed here is suggestive of a period-doubling cascade scenario, as has been reported for several other parallel shear flows (Kreilos & Eckhardt Reference Kreilos and Eckhardt2012; Lustro et al. Reference Lustro, Kawahara, van Veen, Shimizu and Kokubu2019). The stable state at ${{R}_{i}}=395.7816$ is only mildly chaotic and the unstable P-DRW solution provides a fair approximation of the attractor properties. Several complex global bifurcations can be identified along the period-doubling cascade that are determinant to the dynamics, but their nature will be discussed elsewhere on account of the associated intricacies.

Figure 11. Two stages of the period-doubling cascade of P-DRW as illustrated by phase-map projections on the $(\tau _{o},\tau _{i})$ plane at (a) ${{R}_{i}}=395.5$ and (b) ${{R}_{i}}=395.7816$. Shown are the stable attractor (black line) along with the unstable P-DRW solution (red dashed line).

The stability analysis reflected in figure 3 reveals at least two more Hopf bifurcations of the already unstable DRW as ${{R}_{i}}$ is increased beyond the value for the onset of CH-DRW. Solution branches issued from these Hopf points and subsequent bifurcation cascades like the one we report may probably contribute, along with CH-DRW, through global bifurcations involving crises and mergers, to the formation of the strange set that sustains spatiotemporally chaotic dynamics at higher ${{R}_{i}}$ (Krygier, Pughe-Sanford & Grigoriev Reference Krygier, Pughe-Sanford and Grigoriev2021).

5.3. Turbulent dynamics at ${{R}_{i}}=600$

As shown in figure 12 at ${{R}_{i}}=600$, the dynamics in the short parallelogram domain exhibits wild fluctuations alternating rather low-torque/low-kinetic-energy, ostensibly dormant, laminar phases with actively turbulent stages (see also Supplementary Movie 3). If these turbulent transients are governed by the coalescence of a number of temporally chaotic sets such as CH-DRW, it is to be expected that the periodic orbits at their origin should contribute their part to the dynamics. Unfortunately, current computational resources have not allowed continuation of the P-DRW branch this far up in Reynolds number. Nonetheless, existing DRW solutions at ${{R}_{i}}=600$ (points D, E and F indicated in figure 3), all of them unstable, seem to still play a role in scaffolding the transient turbulent state. The only solution that remains of those issued from SN$_1$ (DRW$_{D}$, blue) has a perturbation kinetic energy $\kappa$ somewhat larger than the mean for the statistically steady turbulent state. The other two solutions, resulting from SN$_{3}$ (higher-torque solution DRW$_{E}$ in green and lower-torque solution DRW$_{F}$ in red), have smaller perturbation energy and torque. We will argue that, in spite of their misleadingly high torque and kinetic energy values, the laminar stages of the dynamics are indirectly linked to the TVF (magenta) and DTVF (cyan) solutions.

Figure 12. Chaotic dynamics at ${{R}_{i}}=600$ in the short parallelogram domain of figure 2(c) $(n_1,k_1,n_2,k_2)=(10,2,0,4.5)$. (a) Phase-map projections on the $(\tau _{o},\tau _{i})$ plane. Grey bullet and error bars denote mean and root-mean-square of the fluctuations, respectively. (b) Time series of normalised kinetic energy $\kappa$. The grey line indicates mean kinetic energy. The colour circles/squares and lines indicate DRW, TVF and DTVF solutions labelled in figure 3.

Interesting new phenomena occur when the computational domain is extended in the azimuthal direction to cover the whole circumference of the annulus. Figure 13 shows DNS results in the long computational domain of figure 2(b) characterised by $(n_1,k_1,n_2,k_2)=(1,2,0,4.5)$. The flow field has been initialised with a 10-fold azimuthal replication of the DRW$_{F}$ solution and thence left to evolve freely. The simulation deviates very slowly from DRW$_{F}$ in the beginning, such that snapshot T1, taken at $t=0.2$, is still indistinguishable from DRW$_{F}$. Soon after, instability leads to a sudden surge of kinetic energy. Snapshot T2 is taken at the peak and, as clear from figure 13(b), still markedly preserves the wavelength of the original pattern. The instability driving the burst is therefore superharmonic. However, as the vortex collapses and the energy of the disturbance drops, the wavelength becomes strongly modulated as illustrated by snapshot T3 at the valley of the decay. From this point on, the azimuthal inhomogeneity of the flow field becomes increasingly pronounced as time evolves. Snapshots T4 and T5 show how this inhomogeneity gradually develops into strong localisation. After all remnants of transient dynamics, the statistically turbulent state exhibits azimuthal localisation as exemplified by the instantaneous snapshot T6.

Figure 13. Transient stages of the development of SPT in the long parallelogram domain of figure 2(b) $(n_1,k_1,n_2,k_2)=(10,2,0,4.5)$ at ${{R}_{i}}=600$, starting from a 10-fold azimuthal replication of DRW$_{F}$ (red bullet in figure 12a). (a) Time series of normalised kinetic energy $\kappa$. (b) Selected snapshots of radial vorticity $\omega _r\in [-4000,4000]$ fields at the mid gap $r_{m}$. Subpanels T1–T5 correspond to the points indicated in (a), while T6 is taken beyond all transients for $t\gg 1$. See Supplementary Movie 4.

The flow structure of the instantaneous snapshot of the statistically steady state T6 is examined in greater detail in figure 14. The basic structure within the turbulent region locally displays sinuous wavy vortices resembling DRW but the wavelength exhibits azimuthal modulation. The turbulent band moves from right to left with mean azimuthal phase speed $\langle c_{\theta }\rangle _t=-35.4$ and root-mean-square of the fluctuations ${\sqrt {\langle (c_{\theta }-\langle c_{\theta }\rangle _t)^2\rangle _t}=15.2}$, while the local wave propagation speed is strongly location-dependent, as the travelling-wave solutions that are intermittently and locally approached have each their own characteristic phase speed. This behaviour of the flow field is typical of localised solutions in shear flows (see for example Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013; Brand & Gibson Reference Brand and Gibson2014; Zammert & Eckhardt Reference Zammert and Eckhardt2014; Mellibovsky & Meseguer Reference Mellibovsky and Meseguer2015), and may be mathematically justified by a Wentzel–Kramers–Brillouin-type approach (see Bender & Orszag Reference Bender and Orszag1999, for example). The relation between the stability of small-domain solutions and large-domain flow structure have been the object of recent analysis (Barnett, Gurevich & Grigoriev Reference Barnett, Gurevich and Grigoriev2017; Ritter et al. Reference Ritter, Zammert, Song, Eckhardt and Avila2018).

Figure 14. An instantaneous snapshot of the statistically steady SPT state in the long domain at ${{R}_{i}}=600$. (a) A close up of snapshot T6 of figure 13(b) at mid gap $r_{m}$. (b) Azimuthal vorticity $\omega _{\theta }$ distribution in the $r$–$z$ plane at azimuthal cross-sections $\theta _1 = 0.5$ (bi), $\theta _2 = 2.3$ (bii), $\theta _3 = 3.1$ (biii) and $\theta _4 = 4.3$ (biv) (black vertical lines in panel a). Colour ranges between $\omega _{\theta }\in [-1500,1500]$. From left to right, the black curves in (b) correspond to selected isocontours with (bii) $V/r = \{-20, -80\}$, (biii) $V/r = \{0,-30,-100\}$, (biv) $V/r = 0$. See Supplementary Movie 5, showing a cross-sectional sweep of the instantaneous flow field along the full perimeter.

The local appearance of SSP-related flow structures in streamwise-extended domains has been reported by Deguchi & Hall (Reference Deguchi and Hall2014) in plane Couette flow when analysing the high-Reynolds number asymptotic development of nonlinear equilibria. The observation of the local SSP must in this case be based on the strong correlation existing between the vortex sheet and a critical layer defined by the local (rather than the streamwise-averaged) streamwise velocity field. A similarly high correlation is observed in figure 14 despite the flow being turbulent. The azimuthal vorticity colourmap on a $r$–$z$ azimuthal cross-section near the trailing edge of the turbulent stripe (figure 14biv) reveals a vortex sheet in the vicinity of the inner cylinder. Its spatial arrangement is closely outlined by an appropriate choice of $V/r$ contour (black line), which is also true for figure 14(bii,biii). In the core of the turbulent band, several strong vortex sheets are visible (figure 14biii), while only those close to the outer cylinder seem to persist towards the leading edge (figure 14bii). This supports the shear-driven, rather than centrifugally driven, nature of SSP. The presence of vortex sheets at various radial locations is consistent with the multiplicity of DRW solutions existing at this value of ${{R}_{i}}$, which suggests that each one of them might participate locally in the deployment of the SSP.

The basic flow is linearly unstable at ${{R}_{i}}=600$, which would in principle be at odds with the observation of a laminar flow region in figure 14(a,bi). Note, however, that the $\theta{-}z$ cross-section corresponds to mid gap, while the centrifugal instability is known to develop in the close neighbourhood of the inner cylinder. In fact, examination of the cross-section of figure 14(bi) reveals the presence of weak vortices there. The structure of these nonlinear vortices is best understood by examining an unwrapped $\theta{-}z$ section of the full domain simulation at $r=r_{cu}=(r_{i}+r_{n})/2=7.705$, in the midst of the centrifugally unstable region (figure 15a). A spiral-like flow topology clearly arises away from the turbulent stripe. Despite the dislocations and irregularities, the vortical structures in the laminar region compare favourably with exact spiral (SPI) solutions (Meseguer et al. Reference Meseguer, Mellibovsky, Avila and Marques2009b) of the right axial-azimuthal wavenumbers (and therefore tilt). The base CCF is unstable to a wide range of axial and azimuthal wavenumbers for $(\eta,{{R}_{i}},{{R}_{o}})=(0.883,600,-1200)$, including axisymmetric perturbations with $k\in [2.44,14.26]$ and spirals modes with ${n\leq 8.84}$ ($<9$ if the pattern is to fit exactly an integer number of times around the apparatus). Figure 16 shows a selection of centrifugally driven states that fit exactly in either of the three domains considered.

Figure 15. The morphology of SPT within the centrifugally unstable region of the annulus gap. Same as figure 2, but at $r=r_{cu}=(r_{i}+r_{n})/2=7.705< r_{n}$. Colour ranges between $\omega _{r}\in [-4000,4000]$. See Supplementary Movies: (a) Movie 1; (b) Movie 2; and (c) Movie 3.

Figure 16. Examples of centrifugally driven exact solutions in the three domains for $({{R}_{o}},{{R}_{i}})=(-1200,600)$, visualised through $\omega _{r}\in [-4000,4000]$ colourmaps at $r=r_{cu}=(r_{i}+r_{n})/2=7.705< r_{n}$. Here (a) $(n,k)=(5,7.4)$ SPI ; (b) $(n,k)=(6,5.7)$ SPI; (c) $k=9$ TVF (ci) and $k=4.5$ DTVF (cii).

In particular, figure 16(a) shows one such solution for $(n,k)=(5,7.4)$, evincing that the tilted vortical structures in the laminar region of SPT could well be the result of the interaction of a continuum of spiral solutions of varying wavenumbers around those of the state shown here. The same holds upon inspection of the same $\theta{-}z$ section for the narrow long domain in figure 15(b), for which dislocated spiral patterns of wavenumbers compatible with those of the exact SPI of figure 16(b), corresponding to $(n,k)=(6,5.7)$, also show up away from the turbulent spiral. The azimuthally extended nature of both domains, combined with the spatially adjusting role of the turbulent spiral and, to a lesser degree, of defects and dislocations, allows for the formation of spirals not necessarily restricted to $n\in \mathbb {Z}$. As a result, one must not expect that the spirals observed in the laminar region of SPT conform to a superposition of exact solutions, each one strictly fitting in the domain.

In the small domain (figure 15c), however, quasiaxisymmetric TVF-like structures arise during the laminar stages of the time evolution (figure 15cii) instead of SPI patterns. While SPI solutions of tilt and wavenumbers compatible with those observed in the time evolution exist in the full orthogonal (figure 16a) and narrow but long parallelogram (figure 16b) domains, such solutions do not exist in the short and narrow parallelogram domain (figure 16c). Unstable modes of CCF are exclusively axisymmetric in the small parallelogram, which explains why the laminar stages of DNS display TVF- and DTVF-like structures in the centrifugally unstable region and not spirals. The laminar phases of the time evolution, however, reach well below the torque and kinetic energy levels that are characteristic of TVF or even DTVF. This suggests that it is not the exact states themselves that are being approached but the unstable manifolds of CCF that lead towards them. Axisymmetric modes with $k=4.5$, 9 and 13.5 are the only three unstable modes of CCF in the small parallelogram domain at $({{R}_{o}},{{R}_{i}})=(-1200,600)$. The one fitting twice in the domain ($k=9$) is by far the most unstable, followed by that fitting just once ($k=4.5$). Mode $k=13.5$, which fits three times in the domain, is only marginally unstable and therefore difficult to observe. As the flow decays from turbulent towards CCF, it is repelled along the direction of one of these modes and the phase-map trajectory aligns with the corresponding manifold. As this occurs, the flow fields adopt the shape of one of the unstable modes of CCF and undergo a quasimodal growth that eventually deforms into the nonlinear manifolds connecting to TVF or DTVF. This is what is actually observed during the laminar period of the time evolution. If the centrifugally driven solutions were stable (or only marginally unstable), the approach would be consummate (or at close quarters), but since they are strongly unstable, the trajectory is again repelled well before the actual TVF or DTVF solution can be reached in any of their unstable directions. Supplementary Movie 3, clearly shows how the energy-growing quasiaxisymmetric flow fields of the laminar stages of time evolution indeed undergo a wavy destabilisation before eventually triggering the turbulent burst.

Cyclic turbulent bursts similar to those observed in our small parallelogram have been reported by Coughlin & Marcus (Reference Coughlin and Marcus1996) in TCF at similar values of the parameters and slightly smaller radius ratio using an orthogonal computational domain that, though extended in the azimuthal direction, was incapable of sustaining SPT due to an insufficient axial height. In their computations, a laminar pattern of dislocated/defective spirals forms in the centrifugally unstable region following the linear instability of CCF. These spirals are in turn unstable to a modulational instability that reaches beyond $r_{n}$ into the centrifugally stable region, which, at these values of the parameters, happens to behave as a subcritically unstable shear flow. When the modulational instability exceeds a certain threshold, a gap-filling turbulent burst is triggered. During the laminar part of the cycle a fair amount of energy has been stored in the flow's differential rotation but, as the turbulent burst develops, the energy is quickly transferred to the smaller scales and dissipated by viscosity. Once the mean flow energy has been depleted, the turbulent state has no energy source to rely upon and decays. As the flow relaminarises, the system tries to restore the basic CCF and clutches onto the unstable manifold that leads towards SPI, thereby recommencing the cycle. The same process seems to apply to the dynamics in the narrow and short parallelogram domain presented here, but the spiral instability of CCF is in our case replaced by an axisymmetric instability connecting to TVF due to our restricted azimuthal size. The interplay between the inner (centrifugally unstable) and outer (stable) regions have been shown to also be at play in a recent follow-up paper by Crowley et al. (Reference Crowley, Krygier, Borrero-Echeverry, Grigoriev and Schatz2020) that provided experimental and numerical evidence of a subcritical transition to turbulence in low aspect-ratio TCF. In their set-up, transition is mediated by the interim appearance of yet another type of laminar solutions, namely interpenetrating spirals.

The flow structure associated with the centrifugally unstable modes is manifestly different from that of a vortex sheet resulting from the SSP, and both together constitute the fundamental pieces of the SPT pattern in supercritical counter-rotating TCF. The large-scale laminar–turbulent stripe patterns decisively depend on the high-wavelength azimuthal Fourier modes engendered by the SSP, as suggested by various DNS simulations started from the same SPT state of figure 15(b), but gradually reducing the azimuthal resolution. Truncating to $N=15$ azimuthal Fourier modes, quickly smears the laminar–turbulent interfaces after some initial transients (see Supplementary Movie 6), while the banded pattern subsists permanently when relaxing to $N=30$ (Supplementary Movie 7).

The formation of the large-scale pattern may be understood, to some extent, in the light of a heteroclinic loop that connects the laminar with the turbulent state (and back) spatially in the circumferential direction, rather than temporally. Exact localised solutions in shear flows can be interpreted as homoclinic orbits of the base laminar solution if the problem is tackled as a dynamical system with the extended coordinate playing the role of the (time-like) independent variable (Kirchgässner Reference Kirchgässner1982; Iooss & Pérouème Reference Iooss and Pérouème1993). This is justified by the fact that a linearised analysis can be applied to the localised solution tails (Barnett et al. Reference Barnett, Gurevich and Grigoriev2017; Ritter et al. Reference Ritter, Zammert, Song, Eckhardt and Avila2018). However, numerical and asymptotic analyses at high Reynolds numbers also back the notion that such spatial orbits locally approach the SSP in short streamwise-periodic domains (Deguchi & Hall Reference Deguchi and Hall2014). The linear tails are therefore passive and actually driven by the nonlinear mechanism that sustains the core region of localised turbulence, which can therefore be safely studied even in small periodic domains.

This scenario is of course analogous to what happens for subcritical parallel shear flows, where alternating bands of the laminar and turbulent states conform the stripe pattern. However, there is one essential difference between the two types of banded structures. For parallel shear flows, both the laminar and turbulent flows are stable solutions of the equations when considered in minimal flow units (small periodic domains). The realisation of the banded pattern can therefore be explained by the local bistability of the system, a fact that has been recently used as the basis for a statistical approach employing directed percolation theory (Lemoult et al. Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016). One might thus expect, by analogy, that both a permanent (or at least long-lived) shear-driven turbulent regime and a stable spiral-like state of centrifugal origin must exist in TCF when considering a suitably small parallelogram domain. However, detection of bistability at ${{R}_{i}}=600$ turns out to be utterly elusive. Instead, the laminar and the turbulent states, the former in the shape of axisymmetric TVF-like structures instead of SPI due to modal selection in the minimal flow unit, seem to be embedded in some sort of heteroclinic cyclic loop that drives the dynamics in a chaotic cycle approaching, yet never fully relaxing onto, either state. It would therefore appear that the formation of SPT in supercritical flows might be the result of some sort of large-scale interaction, something that would then be at odds with the current understanding of pattern formation in parallel shear flows. One might speculate that the heteroclinic cycle the dynamics shadow in time within the small domain, is instead deployed in space for azimuthally long domains. As observed by Coughlin & Marcus (Reference Coughlin and Marcus1996), the centrifugally driven laminar state feeds a subcritical shear instability that triggers a turbulent burst. The resulting turbulent state relies on the pre-existence of a laminar flow field from which energy can be extracted for its long-term sustainment while, at the same time, gap-filling turbulence suppresses the centrifugal instability that drives the laminar state. As a result, the turbulent state is short-lived and cannot self-sustain. In the small domain this can only happen by alternating the laminar and turbulent states in time, while SPT patterns can arise in azimuthally long domains as turbulence continuously decays and forms at either one of the laminar–turbulent interfaces. The azimuthally long domain of Coughlin & Marcus (Reference Coughlin and Marcus1996) should then, in principle, be capable of sustaining laminar–turbulent patterns. It might be the case that the tilt of the laminar–turbulent interfaces, which is not compatible with their narrow axially periodic orthogonal domain, is essential to equilibrating the rates of turbulence production and decay. If the characteristic lifetime of the turbulent state is long in comparison with the growth rate of the laminar-state instability that triggers the bursts, no spatiotemporal coexistence can be expected.

In view of the reported observations so far, the narrow long domain is capable of qualitatively reproducing the SPT regime of the full domain, but nothing has been claimed so far as to quantitative accuracy. The large length scale associated with the spiral tilt has been straightforwardly addressed by the use of the parallelogram-shaped domain, and the short length scales of the underlying DRW solutions and the even shorter associated with turbulence are captured by a right choice of $k_2$ and sufficient resolution, respectively. There might, however, be intermediate length scales related to a modulational instability of DRW along the helical direction. If such intermediate scales play an important role in SPT, the elongated domain must be wide enough for the associated coherent structures to decorrelate, as would happen in an infinitely (or sufficiently) long apparatus.

The degree to which the statistics of turbulent signals converge as the domain and/or resolution are modified provides a means of quantifying the inaccuracies associated with domain shape, size and discretisation. Here we have used the kinetic energy ($\kappa$), inner ($\tau _{i}$) and outer ($\tau _{o}$) torque, and azimuthal propagation speed ($c_{\theta }$) signals in the comparison of full domain and narrow-long domain results. The simulations were run for 20 viscous time units past all foreseeable transients to reach the statistically steady turbulent state and then run for an additional 30 and 100 viscous time units, for the full and narrow-long domain, respectively, to collect statistics. A similar resolution density was kept from one domain to the other. All signals were positively checked for stationarity via the augmented Dickey–Fuller test (Dickey & Fuller Reference Dickey and Fuller1979, Reference Dickey and Fuller1981) and the mean and root-mean-square of the fluctuation component (standard deviation) computed. Using stationary bootstrapping (Politis & Romano Reference Politis and Romano1994), 95 % confidence intervals for the two statistics estimators have been determined with automatic block size optimisation (Politis & White Reference Politis and White2004; Patton, Politis & White Reference Patton, Politis and White2009). Since fluctuation amplitude depends on domain size (note that a system of infinite aspect ratio would present no fluctuations on account of the averaging properties of normalised aggregate quantities as we are monitoring here), the signal variance needs to be scaled accordingly. This is the same as assuming that the signal fully decorrelates over the height of the domain, and that the variance is therefore the composition of the variances of as many independent, randomly distributed, identical signals as times the narrow domain fits in the full domain. Fluctuation amplitudes are therefore dependent on domain size, and must be scaled. Results are summarised in table 1.

Table 1. Statistics of several time signals (mean $\langle \bullet \rangle$ and root-mean-square of the fluctuation component $\sigma _{\bullet }=\sqrt {\langle (\bullet -\langle \bullet \rangle )^2\rangle }$, this latter corrected for domain size) of the statistically steady turbulent state as computed in the full orthogonal domain and in the azimuthally long but axially narrow parallelogram domain. The uncertainty margins correspond to 95 % confidence intervals obtained through stationary bootstrapping.

The runs are long enough for the mean to be converged to within $1\,\%$ for all time series. The uncertainty in the root-mean-square is larger, ranging from $3\,\%$ to $9\,\%$. In any case, it is sufficient to claim, with $95\,\%$ certainty, that the narrow and long parallelogram domain produces results that are quantitatively different from the full domain. Taking the mean of the statistic (mean or root-mean-square) as the best estimate, however, bounds the inaccuracy to within $5\,\%$ for the mean and $10\,\%$ for the root-mean-square, which is fair but not extremely precise. While the narrow long domain of figure 2(b) provides reasonable quantitative agreement, wider parallelogram domains would be required to further enhance accuracy, particularly so regarding second-order statistics.