# Self-sustainment of coherent structures in counter-rotating Taylor-Couette
flow

B. Wang1  R. Ayats1  K. Deguchi2  F. Mellibovsky1 and A. Meseguer1
<EMAIL_ADDRESS>1 Departament de Física, Universitat Politècnica
de Catalunya, 08034, Barcelona, Spain 2 School of Mathematics, Monash
University, VIC 3800, Australia

###### Abstract

We investigate the local self-sustained process underlying spiral turbulence
in counter-rotating Taylor-Couette flow using a periodic annular domain,
shaped as a parallelogram, two of whose sides are aligned with the cylindrical
helix described by the spiral pattern. The primary focus of the study is
placed on the emergence of drifting-rotating waves (drw) that capture, in a
relatively small domain, the main features of coherent structures typically
observed in developed turbulence. The transitional dynamics of the subcritical
region, far below the first instability of the laminar circular Couette flow,
is determined by the upper and lower branches of drw solutions originated at
saddle-node bifurcations. The mechanism whereby these solutions self-sustain,
and the chaotic dynamics they induce, are conspicuously reminiscent of other
subcritical shear flows. Remarkably, the flow properties of drw persist even
as the Reynolds number is increased beyond the linear stability threshold of
the base flow. Simulations in a narrow parallelogram domain stretched in the
azimuthal direction to revolve around the apparatus a full turn confirm that
self-sustained vortices eventually concentrate into a localised pattern. The
resulting statistical steady state satisfactorily reproduces qualitatively,
and to a certain degree also quantitatively, the topology and properties of
spiral turbulence as calculated in a large periodic domain of sufficient
aspect ratio that is representative of the real system.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see Keyword PDF for the full list). Other
classifications will be added at the same time.

MSC Codes  (Optional) Please enter your MSC Codes here

##  1 Introduction

Nonlinear shear flow instabilities often lead to intermittency, a phenomenon
that involves the coexistence of laminar and turbulent flow regions. The
resulting turbulent-laminar patterns take either the form of localised
turbulent patches surrounded by otherwise laminar flow or, surprisingly more
orderly, of coherent stripes that exhibit sharp interfaces tilted with respect
to the main direction of the flow. A paradigm for coherent intermittency in
rotating shear flows is the so-called Spiral Turbulence regime (spt) that
appears in Taylor-Couette flow, i.e. the fluid flow between independently-
rotating coaxial cylinders. The spt regime consists of a turbulent helix that
forms a coil within the apparatus gap with a well-defined pitch and rotates at
a fairly constant angular speed. This peculiar flow structure, first
discovered in the 1960s (Coles & Van Atta, 1967) and declared a puzzling
phenomenon by Feynman himself (Feynman, 1964), had not been reproduced
numerically until rather recently by means of very costly direct numerical
simulation (dns) (Meseguer et al., 2009b; Dong, 2009; Dong & Zheng, 2011).

Similar oblique laminar-turbulent stripe patterns are common in other shear
flows featuring two extended space directions (Prigent et al., 2002; Duguet et
al., 2010; Tuckerman et al., 2014). For a thorough review on the rich variety
of existing intermittent shear flow phenomena, we refer the reader to the
monograph by Tuckerman et al. (2020) and references therein. Recent attempts
at elucidating the stripe formation mechanism have mostly been confined to
relatively simple parallel shear flows such as plane Couette or plane
Poiseuille flows. Both problems exhibit subcritical transition to turbulence,
namely transition in the absence of a linear instability of the base flow,
which is best tackled employing dynamical systems theory. In this framework,
simple solutions to the Navier-Stokes equations, often called exact coherent
structures (ecs), are shown to naturally play a central role in organising the
transitional and turbulent dynamics (Kerswell, 2005; Eckhardt et al., 2007;
Kawahara et al., 2012; Graham & Floryan, 2021).

The last three decades have witnessed overwhelming scientific activity in the
search for ecs in many subcritical shear flows (Nagata, 1990; Clever & Busse,
1992, 1997; Waleffe, 2003; Faisst & Eckhardt, 2003; Wedin & Kerswell, 2004),
following the discovery of a Self-Sustained Process (ssp) for coherent
structures in the absence of linear instability of the laminar flow (Waleffe,
1997). The process, which occurs at relatively short length-scales and can
therefore be observed in small periodic domains, consists in a cyclic feedback
mechanism whereby streamwise vortices generate streaks through the lift-up
mechanism, which in turn become unstable to three-dimensional waves that feed
energy back into the streamwise vortices (Boberg & Brosa, 1988; Hamilton et
al., 1995; Grossmann, 2000). The instability of the streaks responds to an
inviscid mechanism by which three-dimensional waves are strongly amplified at
a critical layer, i.e. where the streak speed coincides with the wave speed
(Wang et al., 2007; Hall & Sherwin, 2010; Deguchi & Hall, 2015).

For relatively small periodic boxes, the relationship between the exact
solution and the dynamics is understood to some extent. Of the ecs emanated
from saddle-node bifurcations, the lower (or saddle) branch typically dictates
the topology and amplitude of flow perturbations that are capable of
triggering transition, while the upper (or nodal) branch generally regulates
-or at least participates in- the formation of the turbulent set. In the
subcritical transition problem, the infinite-dimensional Navier-Stokes phase
space typically contains two stable (or metastable) invariant sets: the steady
laminar base flow and the chaotic/turbulent state. The basins of attraction of
these two sets meet along a codimension-1 manifold, usually known as edge of
chaos (Itano & Toh, 2001; Skufca et al., 2006; Schneider et al., 2008), where
the lower-branch solutions belong. The upper branch ecs are almost always
unstable and their participation in the generation of the chaotic set is not
easily dissected. Occasionally, however, upper-branch ecs are linearly stable,
if only within sufficiently small domains or tightly constrained symmetry
conditions, in a neighbourhood of the saddle-node bifurcation (Clever & Busse,
1997; Mellibovsky & Eckhardt, 2011, 2012). In these cases, further theoretical
progress can be made as the path towards chaotic dynamics admits a simpler
analysis that can draw from a parallel with low-dimensional dynamical systems
(Mellibovsky & Eckhardt, 2012; Kreilos & Eckhardt, 2012; Lustro et al., 2019).

Transition and pattern formation in Taylor-Couette flow are more involved than
for parallel shear flows due to the interplay of shear and rotation (Andereck
et al., 1986). A notable difference with respect to merely shear-driven flows
is that spt also persists in the supercritcal regime of counter-rotating
Taylor-Couette flow, beyond the linear instability of the laminar circular
Couette flow (Prigent et al., 2002; Meseguer et al., 2009b). Therefore, both
the shear and the centrifugal instabilities contribute their share to the
generation of streamwise vorticity, but this fact has gone largely unnoticed
in the literature, where the origin of the stripe may be solely explained by
the stability of both the basic flow and the autonomous vortex emerged from
the ssp.

The aim of this paper is to investigate the dynamics induced by ecs driven by
the ssp and to ascertain whether they may be held responsible for the
formation of the spt regime observed in the centrifugally-unstable region of
parameter space. The most natural place to look for nonlinear solutions is at
the linear critical point of the base flow, circular Couette flow in our case.
However, the nonlinear spiral solutions thus identified by Meseguer et al.
(2009a) are only very mildly subcritical, which implies that it is the
centrifugal instability rather than the ssp that drives them. A wealth of
solutions predicted by weakly nonlinear theory (Chossat & Iooss, 1994)
followed the discovery of the subcritical spirals Deguchi & Altmeyer (2013),
but none contributed to enlarge the known region of subcriticality. Shortly
after, however, a highly subcritical three-dimensional rotating-wave solution
was found by Deguchi et al. (2014), although its dynamical relevance was not
investigated.

As noted earlier, the computation of laminar-turbulent banded patterns is very
costly, such that using narrow orthogonal domains suitably tilted to align
with the stripes has been decisive to the study of this kind of laminar-
turbulent patterns (Barkley & Tuckerman, 2005; Shi et al., 2013; Reetz et al.,
2019; Paranjape et al., 2020). The approach is easily undertaken for parallel
shear flows, as a mere change in the direction of the base flow suffices, but
more fundamental code modifications are necessary for cylindrical and annular
geometries. The required coordinate change was generalised by Deguchi &
Altmeyer (2013) to compute ecs in parallelogram-shaped domains wrapped within
an annular geometry. Since their method of directly solving the nonlinear
algebraic equations was only applicable to travelling-wave solutions,
developing a dns code in generalised parallelogram-shaped periodic domains is
requisite to efficiently capture spt dynamics, but has hitherto not been
attempted to our best knowledge.

The outline of the paper is as follows. The problem formulation is given in
§2, alongside a description of the generalised parallelogram-shaped domain,
its application to the spectral space discretisation, and the numerical
methods employed for evolving the equations in time, and the coupling of the
time-stepper with a Poincaré-Newton-Krylov solver. Then §3 briefly summarises
the geometrical and physical parameters used for the numerical calculations,
and justifies the specific choice of the domain shape. Since the rotating-wave
solutions found by Deguchi et al. (2014) were computed in a classical
orthogonal small periodic domain, the initial task in §4 is the exploration of
bifurcation scenarios leading to solutions of the same family in small
parallelogram domains. In §5, the dynamical relevance of the solution is first
analysed in the subcritical regime to identify the ssp and the onset of
chaotic dynamics. The possible relevance of these solutions to the
supercritical spt regime is then discussed. Finally, the main findings are
summarised in §6 along with concluding remarks.

##  2 Formulation of the problem

Consider an incompressible fluid of dynamic viscosity $\mu$ and density
$\varrho$ (kinematic viscosity $\nu=\mu/\varrho$) completely filling the gap
between two concentric rotating cylinders whose inner and outer radii and
angular velocities are $r^{*}_{\rm{i}}$, $r^{*}_{\rm{o}}$ and
$\Omega_{\rm{i}}$, $\Omega_{\rm{o}}$, respectively. A full set of independent
dimensionless parameters characterising the problem are the radius ratio
$\eta=r^{*}_{\rm{i}}/r^{*}_{\rm{o}}$, which fixes the geometry of the annulus,
and the circular Couette flow inner and outer Reynolds numbers
${\mathrm{R}_{\rm i}}=dr^{*}_{\rm{i}}\Omega_{\rm{i}}/\nu$ and
${\mathrm{R}_{\rm o}}=dr^{*}_{\rm{o}}\Omega_{\rm{o}}/\nu$, where
$d=r^{*}_{\rm{o}}-r^{*}_{\rm{i}}$ is the gap between the cylinders.
Henceforth, all variables will be rendered dimensionless using $d$,
$d^{2}/\nu$, and $\nu^{2}/d^{2}$ as units for space, time and the reduced
pressure ($p=p^{*}/\varrho$), respectively. The Navier–Stokes equation, and
the incompressibility and the zero axial net massflux conditions become

$\displaystyle\partial_{t}\mbox{\bf v}+(\mbox{\bf v}\cdot\nabla)\mbox{\bf v}$
$\displaystyle=$ $\displaystyle-\nabla p+\nabla^{2}\mbox{\bf
v}+f\,\hat{{\bm{z}}},$ (1) $\displaystyle\quad\nabla\cdot\mbox{\bf v}$
$\displaystyle=$ $\displaystyle 0,$ (2) $\displaystyle Q(\mbox{\bf
v})=\int_{0}^{2\pi}{\int_{r_{\rm{i}}}^{r_{\rm{o}}}{(\mbox{\bf
v}\cdot\hat{{\bm{z}}})\,r\,dr\,dz}}$ $\displaystyle=$ $\displaystyle 0,$ (3)

where the axial forcing term $f=f(t)$ in (1) is instantaneously adjusted to
fulfil the constraint imposed by (3), $\mbox{\bf
v}=(U,V,W)=U\,\hat{{\bm{r}}}+V\,\hat{{\bm{\theta}}}+W\,\hat{{\bm{z}}}$ is the
velocity of the fluid expressed in cylindrical coordinates $(r,\theta,z)$,
which satisfies no-slip boundary conditions at the cylinder walls

$\mbox{\bf v}|_{r=r_{\mathrm{i}}}=(0,{\mathrm{R}_{\rm i}},0),\quad\mbox{\bf
v}|_{r=r_{\rm o}}=(0,{\mathrm{R}_{\rm o}},0),$ (4)

with

$r_{\mathrm{i}}=\frac{r^{*}_{\mathrm{i}}}{d}=\frac{\eta}{1-\eta},\quad
r_{\mathrm{o}}=\frac{r^{*}_{\mathrm{o}}}{d}=\frac{1}{1-\eta},$ (5)

the non-dimensional inner and outer radii, respectively. The basic, laminar
and steady circular Couette flow, henceforth referred to as ccf, is

$\mbox{\bf v}_{\rm b}=U_{\rm b}\,\hat{{\bm{r}}}+V_{\rm
b}\,\hat{{\bm{\theta}}}+W_{\rm
b}\,\hat{{\bm{z}}}=\displaystyle\left(Ar+\frac{B}{r}\right)\,\hat{{\bm{\theta}}},\;\;p_{\rm
b}(r)=\int\frac{V_{\rm b}^{2}}{r}\,{\rm d}r,\;\;f_{\rm b}=0,$ (6)

with $A=({\mathrm{R}_{\rm o}}-\eta{\mathrm{R}_{\rm i}})/(1+\eta)$ and
$B=\eta({\mathrm{R}_{\rm i}}-\eta{\mathrm{R}_{\rm
o}})/\left[(1-\eta)(1-\eta^{2})\right]$. In what follows we express the
velocity and pressure fields as

$\mbox{\bf v}=\mbox{\bf v}_{\rm b}(r)+\mathbf{u}(r,\theta,z;t),\quad p=p_{\rm
b}(r)+q(r,\theta,z;t).$ (7)

The fields $q$ and
$\mathbf{u}=u\,\hat{{\bm{r}}}+v\,\hat{{\bm{\theta}}}+w\,\hat{{\bm{z}}}$ are
the deviations from the equilibrium ccf solution that, after formal
substitution of (7) into (1), satisfy

$\displaystyle\partial_{t}{\bf u}$ $\displaystyle=$ $\displaystyle-\nabla
q+\nabla^{2}{\bf u}-({\bf v}_{\rm b}\cdot\nabla){\bf u}-({\bf
u}\cdot\nabla){\bf v}_{\rm b}-({\bf u}\cdot\nabla){\bf u},+f\,\hat{{\bm{z}}}$
(8) $\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle Q(\mbox{\bf u})$ $\displaystyle=$ $\displaystyle 0.$ (10)

Although this solenoidal boundary value problem is naturally formulated in
cylindrical polar coordinates $(r,\theta,z)$, the coherent flows addressed in
this work are better captured and more efficiently represented numerically on
parallelogram domains such as the one depicted in figure 1. This type of
domains have been recently vindicated as minimal flow units to capture mixed
spiral modes and rotating-travelling waves with arbitrary axial-azimuthal
wavefront orientation in tcf (Deguchi & Altmeyer, 2013; Deguchi & Hall, 2015;
Ayats et al., 2020a). The parallelogram domain is bounded by two consecutive
($2\pi$-shifted) wavefront loci and might be naturally parametrised by
introducing the new coordinates

$\xi=n_{1}\theta+k_{1}z,\quad\zeta=n_{2}\theta+k_{2}z,$ (11)

or, conversely,

$\theta=\displaystyle\frac{k_{2}\xi-k_{1}\zeta}{n_{1}k_{2}-n_{2}k_{1}},\quad
z=\displaystyle\frac{-n_{2}\xi+n_{1}\zeta}{n_{1}k_{2}-n_{2}k_{1}}.$ (12)

Henceforth, we reformulate the boundary value problem (8-9) within the
parallelogram assuming the pressure and velocity fields $q$ and $\mathbf{u}$
are $2\pi$-periodic in the two new coordinates $\xi$ and $\zeta$, thus
satisfying

$\displaystyle q(r,\xi+2\pi,\zeta;t)$ $\displaystyle=q(r,\xi,\zeta+2\pi;t)$
$\displaystyle=q(r,\xi,\zeta;t),$ (13)
$\displaystyle\mathbf{u}(r,\xi+2\pi,\zeta;t)$
$\displaystyle=\mathbf{u}(r,\xi,\zeta+2\pi;t)$
$\displaystyle=\mathbf{u}(r,\xi,\zeta;t).$ (14)

In what follows, we numerically discretise $q$ and $\mathbf{u}$ within the
annular-parallelogram domain

$(r,\xi,\zeta)\in\left[r_{\rm{i}},r_{\rm{o}}\right]\times[0,2\pi]\times[0,2\pi],$
(15)

where the inner and outer radii of the cylinders are explicitly given in (5).

Figure 1: Sketch of the parallelogram domain introducing the new variables
$(\xi,\zeta)$ that replace the usual azimuthal and axial coordinates
$(\theta,z)$.

We characterise flows by their associated normalised torque at the inner and
outer cylinders, $\tau_{\rm{i}}$ and $\tau_{\rm{o}}$,

$\displaystyle\tau_{\rm i,\rm
o}=\left.1+\frac{\partial_{r}(r^{-1}\bar{v})}{\partial_{r}(r^{-1}v_{\rm
b})}\right|_{r=r_{\rm i},r_{\rm o}},$ (16)

where $\bar{v}$ is the averaged azimuthal velocity in the angular and axial
directions. With this definition, $\tau_{\rm{i}}=\tau_{\rm{o}}=1$ for ccf.
Similarly, we will also characterise the flows by the normalised kinetic
energy of the perturbation velocity field,

$E=\frac{1}{2E_{\rm
b}}\int_{0}^{2\pi}\mathrm{d}\zeta\int_{0}^{2\pi}\mathrm{d}\xi\int_{r_{\rm{i}}}^{r_{\rm{o}}}\mbox{\bf
u}\cdot\mbox{\bf u}\,r\,\mathrm{d}r,$ (17)

where $E_{b}$ is the kinetic energy of ccf, namely

$E_{\rm
b}=\frac{1}{2}\int_{0}^{2\pi}\mathrm{d}\zeta\int_{0}^{2\pi}\mathrm{d}\xi\int_{r_{\rm{i}}}^{r_{\rm{o}}}\mbox{\bf
v}_{\rm b}\cdot\mbox{\bf v}_{\rm
b}\,r\,\mathrm{d}r=2\pi^{2}\int_{r_{\rm{i}}}^{r_{\rm{o}}}\left(Ar+\frac{B}{r}\right)^{2}r\,\mathrm{d}r.$
(18)

###  2.1 Direct numerical simulations in the annular-parallelogram domain

The nonlinear boundary value problem (8-9) is discretised using a solenoidal
Petrov-Galerkin scheme formerly formulated by Meseguer et al. (2007), and
suitably adapted to the annular-parallelogram domain (15). In the transformed
domain, the solenoidal velocity perturbation is approximated by means of a
Fourier $\times$ Fourier $\times$ Chebyshev spectral expansion ${\bf u}_{\rm
s}$ of order $N\times L\times M$ in $\xi\times\zeta\times r$, respectively, of
the form

${\bf u}_{\rm
s}(r,\xi,\zeta;t)=\sum_{\ell=-L}^{L}\sum_{n=-N}^{N}\sum_{m=0}^{M}a_{\ell
nm}^{(1)}(t)\bm{\Phi}_{\ell nm}^{(1)}(r,\xi,\zeta)+a_{\ell
nm}^{(2)}(t)\bm{\Phi}_{\ell nm}^{(2)}(r,\xi,\zeta).$ (19)

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
$\mathbf{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 condiseration that the three velocity components
are not independent, but linked by the solenoidal condition. The vector fields
$\bm{\Phi}_{\ell nm}^{(\imath)}$ constitute the elements of the trial basis of
solenoidal vector fields of the form

$\bm{\Phi}_{\ell
nm}^{(\imath)}(r,\xi,\zeta)=\mathrm{e}^{\mathrm{i}(n\xi+\ell\zeta)}\mathbf{u}_{\ell
nm}^{(\imath)}(r)=\mathrm{e}^{\mathrm{i}(n\xi+\ell\zeta)}\left(u_{\ell
nm}^{(\imath)},v_{\ell nm}^{(\imath)},w_{\ell nm}^{(\imath)}\right),$ (20)

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

$\left(\partial_{r}+\frac{1}{r}\right)u_{\ell
nm}^{(\imath)}+\frac{\mathrm{i}}{r}\left(n_{1}n+n_{2}\ell\right)v_{\ell
nm}^{(\imath)}+{\rm i}\left(k_{1}n+k_{2}\ell\right)w_{\ell nm}^{(\imath)}=0.$
(21)

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

$\bm{\Phi}_{\ell nm}^{(\imath)}(r_{\rm{i}},\xi,\zeta;t)=\bm{\Phi}_{\ell
nm}^{(\imath)}(r_{\rm{o}},\xi,\zeta;t)=\mathbf{0}.$ (22)

In what follows, we define the transformed radial coordinate

$x(r)=2r-\frac{1+\eta}{1-\eta},$ (23)

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

${\rm h}_{m}(r)=(1-x^{2}){\rm T}_{m}(x),\ {\rm g}_{m}(r)=(1-x^{2})^{2}{\rm
T}_{m}(x),$ (24)

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

${\rm h}_{m}(r_{\rm i})={\rm h}_{m}(r_{\rm o})={\rm g}_{m}(r_{\rm i})={\rm
g}_{m}(r_{\rm o})={\rm D}{\rm g}_{m}(r_{\rm i})={\rm D}{\rm g}_{m}(r_{\rm
o})=0,$ (25)

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 (21) and (22). For $nn_{1}+\ell n_{2}=0$, two
such vector fields are

$\mathbf{u}_{\ell nm}^{(1)}(r)=(0\,,\,{\rm h}_{m}\,,\,0),\quad\mathbf{u}_{\ell
nm}^{(2)}(r)=(-\mathrm{i}(nk_{1}+\ell k_{2})r{\rm g}_{m}\,,\,0\,,\,{\rm
D}[r{\rm g}_{m}]+{\rm g}_{m}),$ (26)

with the third component of $\mathbf{u}_{\ell nm}^{(2)}$ replaced by ${\rm
h}_{m}$ whenever $nk_{1}+\ell k_{2}=0$. Finally, for ${nn_{1}+\ell n_{2}\neq
0}$, the solenoidal basis is

$\displaystyle\mathbf{u}_{\ell nm}^{(1)}(r)$
$\displaystyle=(-\mathrm{i}(nn_{1}+\ell n_{2}){\rm g}_{m}\,,\,{\rm D}[r{\rm
g}_{m}]\,,\,0),$ (27) $\displaystyle\mathbf{u}_{\ell nm}^{(2)}(r)$
$\displaystyle=(0,-\mathrm{i}(nk_{1}+\ell k_{2})r{\rm
h}_{m}\,,\,\mathrm{i}(nn_{1}+\ell n_{2}){\rm h}_{m}),$ (28)

except that the third component of $\mathbf{u}_{\ell nm}^{(2)}$ is replaced by
${\rm h}_{m}$ when $nk_{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 $\bm{\Phi}$ and $\bm{\Psi}$, respectively,
over the annular-parallelogram domain (15)

$(\bm{\Psi},\bm{\Phi})=\int_{r_{\rm{i}}}^{r_{\rm{o}}}r\,\rm{d}r\int_{0}^{2\pi}\rm{d}\xi\int_{0}^{2\pi}\rm{d}\zeta\;\bm{\Psi}^{\dagger}\cdot\bm{\Phi}.$
(29)

Accordingly, we consider the dual basis for the projection space. In
particular, the basis for the case $nn_{1}+\ell n_{2}=0$ is

$\displaystyle\mathbf{\tilde{u}}_{\ell 0m}^{(1)}(r)$ $\displaystyle=$
$\displaystyle{\rm w}(0\,,\,r{\rm h}_{m}\,,\,0),$ (30)
$\displaystyle\mathbf{\tilde{u}}_{\ell 0m}^{(2)}(r)$ $\displaystyle=$
$\displaystyle r^{-2}{\rm w}\left(\mathrm{i}(nk_{1}+\ell k_{2}){\rm
g}_{m}\,,\,0\,,\,{\rm D}_{+}{\rm g}_{m}+2r^{-1}(1-x^{2}+rx){\rm
h}_{m}\right),$ (31)

where ${\rm D}_{+}={\rm D}+r^{-1}$, and the third component of
$\mathbf{\tilde{u}}_{\ell 0m}^{(2)}$ is replaced by $r{\rm h}_{m}$ if
$nk_{1}+\ell k_{2}=0$. Similarly, the basis for the case $nn_{1}+\ell
n_{2}\neq 0$ is

$\displaystyle\mathbf{\tilde{u}}_{\ell nm}^{(1)}(r)$ $\displaystyle=$
$\displaystyle{\rm w}\left((nn_{1}+\ell n_{2})\,r{\rm g}_{m}\,,\,r{\rm
D}_{+}[r{\rm g}_{m}]+2xr^{2}{\rm h}_{m}\,,\,0\right),$ (32)
$\displaystyle\mathbf{\tilde{u}}_{\ell nm}^{(2)}(r)$ $\displaystyle=$
$\displaystyle{\rm w}\left(0\,,\,\mathrm{i}(nk_{1}+\ell k_{2})r^{2}{\rm
h}_{m}\,,\,-\mathrm{i}\,(nn_{1}+\ell n_{2})\,r{\rm h}_{m}\right).$ (33)

These projection basis elements contain the Chebyshev weight function ${\rm
w}(x)$ so that the resulting radial integration involved in (29) can be
computed exactly by suitable quadrature formulas (Moser et al., 1983; Meseguer
et al., 2007; Canuto et al., 2007, 2010; Meseguer, 2020).

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

$\left(\bm{\Psi}_{\ell nm}^{(\imath)}\;,\;\partial_{t}{\bf u}_{\rm
S}\right)=\left(\bm{\Psi}_{\ell nm}^{(\imath)}\;,\nabla^{2}{\bf u}-({\bf
v}_{\rm b}\cdot\nabla){\bf u}-({\bf u}\cdot\nabla){\bf v}_{\rm b}-({\bf
u}\cdot\nabla){\bf u}+f\;\hat{{\bm{z}}}\right),$ (34)

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

$\displaystyle\mathbb{A}_{\ell
nm\imath}^{pqr\jmath}\;\dot{a}_{pqr}^{(\jmath)}$ $\displaystyle=$
$\displaystyle\mathbb{B}_{\ell
nm\imath}^{pqr\jmath}\;a_{pqr}^{(\jmath)}+\mathbb{N}_{\ell
nm}^{(\imath)}(\mathbf{a})+f\,\mathbb{F}_{\ell nm}^{(\imath)}$ (35)
$\displaystyle Q_{m}^{r}a_{00r}^{(2)}$ $\displaystyle=$ $\displaystyle 0$ (36)

for the axial forcing $f(t)$ and amplitudes $a_{\ell nm}^{(\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 (19). The quadratic form
$\mathbb{N}_{\ell nm}^{(\imath)}(\mathbf{a})$ appearing in (35) corresponds to
the projection of the nonlinear convective term, $\left(\bm{\Psi}_{\ell
nm}^{(\imath)}\;,({\bf u}\cdot\nabla){\bf u}\right)$, that is computed
pseudospectrally using Orszag’s $3/2$-dealiasing rule (Canuto et al., 2007,
2010). Overall, the resulting stiff system of ode 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

$\begin{array}[]{l}{\bf u}_{\rm s}(r,\xi,\zeta;t)=\\\
\displaystyle{=\sum_{\ell=-L}^{L}\sum_{n=-N}^{N}\sum_{m=0}^{M}\left[\check{a}_{\ell
nm}^{(1)}\bm{\Phi}_{\ell nm}^{(1)}(r,\xi,\zeta)+\check{a}_{\ell
nm}^{(2)}\bm{\Phi}_{\ell
nm}^{(2)}(r,\xi,\zeta)\right]\mathrm{e}^{\displaystyle\mathrm{i}n(\xi-
c_{\xi}t)}\mathrm{e}^{\displaystyle\mathrm{i}\ell(\zeta-
c_{\zeta}t)}},\end{array}$ (37)

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

$a_{\ell nm}^{(\imath)}(t)=\check{a}_{\ell
nm}^{(\imath)}\mathrm{e}^{\displaystyle-\mathrm{i}nc_{\xi}t}\mathrm{e}^{\displaystyle-\mathrm{i}\ell
c_{\zeta}t},$ (38)

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 (37) in (8) followed by Hermitian
projection leads to the system of nonlinear algebraic equations for the
travelling wave coefficients $\check{a}_{\ell nm}^{(\imath)}$
($\check{\mathbf{a}}$ for brevity):

$\displaystyle\left[\mathbb{B}_{\ell
nm\imath}^{pqr\jmath}+\mathrm{i}(nc_{\xi}+\ell c_{\zeta})\mathbb{A}_{\ell
nm\imath}^{pqr\jmath}\right]\check{a}_{pqr}^{(\imath)}+\mathbb{N}_{\ell
nm}^{(\imath)}(\check{\mathbf{a}})+\check{f}\,\mathbb{F}_{\ell nm}^{(\imath)}$
$\displaystyle=$ $\displaystyle 0,$ (39) $\displaystyle
Q_{m}^{r}\check{a}_{00r}^{(2)}$ $\displaystyle=$ $\displaystyle 0,$ (40)

where $\check{\mathbf{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 to 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, 2011). The system is solved using a matrix-free
Newton-Krylov method (Kelley, 2003), implicitly using gmres as matrix-free
solver (Trefethen & Bau, 1997). The converged nonlinear solutions are tracked
in parameter space using pseudo-arclength continuation schemes (Kuznetsov,
2004).

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

$\displaystyle a_{\ell nm}^{(\imath)}(t)$ $\displaystyle=$
$\displaystyle\left(\check{a}_{\ell nm}^{(\imath)}+\varepsilon_{\ell
nm}^{(\imath)}\mathrm{e}^{\sigma
t}\right)\mathrm{e}^{\displaystyle-\mathrm{i}nc_{\xi}t}\mathrm{e}^{\displaystyle-\mathrm{i}\ell
c_{\zeta}t},$ (41) $\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\check{f}+\phi\mathrm{e}^{\sigma t},$ (42)

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

$\displaystyle\left(\sigma-\mathrm{i}nc_{\xi}-\mathrm{i}\ell
c_{\zeta}\right)\mathbb{A}_{\ell nm\imath}^{pqr\jmath}\varepsilon_{\ell
nm}^{(\imath)}$ $\displaystyle=$ $\displaystyle\mathbb{B}_{\ell
nm\imath}^{pqr\jmath}\varepsilon_{\ell
nm}^{(\imath)}+\mathbf{D}_{\mathbf{a}}\mathbb{N}_{\ell
nm\imath}^{pqr\jmath}(\check{\mathbf{a}})\varepsilon_{\ell
nm}^{(\imath)}+\phi\,\mathbb{F}_{\ell nm}^{(\imath)},$ (43) $\displaystyle
Q_{m}^{r}\varepsilon_{00r}^{(2)}$ $\displaystyle=$ $\displaystyle 0$ (44)

where $\mathbf{D}_{\mathbf{a}}\mathbb{N}_{\ell
nm\imath}^{pqr\jmath}(\check{\mathbf{a}})$ is the linear action of the
Jacobian of $\mathbb{N}$ evaluated at the travelling wave solution
$\check{\mathbf{a}}$, and
$(p,q,r,\jmath)\in[0,L]\times[-N,N]\times[0,M]\times\\{1,2\\}$. This linear
action therefore allows to formulate the generalised eigenvalue problem (43)
in a matrix-free form, so that Arnoldi methods (Trefethen & Bau (1997)) can be
applied to compute the leading eigenvalues $\sigma_{j}$ and their associated
eigenvectors $\mathbf{\varepsilon}_{j}=\\{\varepsilon_{\ell
mn}^{(\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,
2015). 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 $\mathbb{P}$:

$\mathbf{a}\rightarrow\widetilde{\mathbf{a}}=\mathbb{P}(\mathbf{a})=\mathbb{T}(\Delta\xi,\Delta\zeta)\bm{\varphi}(\mathbf{a};\mathrm{T}),$
(45)

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

$\left[\mathbb{T}(\Delta\xi,\Delta\zeta)\mathbf{a}\right]_{\ell
nm}^{(\imath)}=\mathrm{e}^{\displaystyle-\mathrm{i}n\Delta\xi}\mathrm{e}^{\displaystyle-\mathrm{i}\ell\Delta\zeta}a_{\ell
nm}^{(\imath)}$ (46)

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.

##  3 Choice of geometrical and physical parameters

The present study is done for the same cylinder radius ratio of $\eta=0.883$
as employed by Andereck et al. (1986) in their Taylor-Couette apparatus. The
outer Reynolds number is fixed to ${\mathrm{R}_{\rm o}}=-1200$, at which
Meseguer et al. (2009b) found numerically that spt is sustained within
sufficiently large domains for values of the inner cylinder Reynolds number
ranging in ${\mathrm{R}_{\rm i}}\in[540,640]$ (see line $\Gamma_{2}$ in figure
3 of that paper). Below the lower threshold for spt, spot-like intermittency
is observed for an interval, bordering on the other side with a regime of
interpenetrating-spirals that prevails in the range ${\mathrm{R}_{\rm
i}}\in[450,480]$.

The critical point above which circular Couette flow becomes linearly unstable
with respect to non-axisymmetric perturbations occurs at ${\mathrm{R}_{\rm
i}}=447.35$. The bifurcating nonlinear spiral-wave solution branches found by
Meseguer et al. (2009a) can only be continued down to ${\mathrm{R}_{\rm
i}}=445.65$, which makes them qualify as only mildly subcritical. While none
of the mixed-mode solutions detected by Deguchi & Altmeyer (2013) reached this
value of ${\mathrm{R}_{\rm i}}$, the subcritical rotating waves computed by
Deguchi et al. (2014) do indeed exist in a much wider region of the linearly
stable regime, extending to as low as ${\mathrm{R}_{\rm i}}>377.3$.

Figure 2a shows radial vorticity, ${\omega_{r}=(\nabla\times\mathbf{v})_{r}}$,
colour maps for a snapshot of the statistically steady spt obtained at
$\rm{R}_{\rm{i}}=600$ in a sufficiently large computational domain following
Meseguer et al. (2009b).

(a) | (c)
---|---
|  |
---

(b)

Figure 2: Instantaneous fields of statistically steady states obtained with
dns for $(\eta,{\mathrm{R}_{\rm o}},{\mathrm{R}_{\rm i}})=(0.883,-1200,600)$.
All pannels show radial vorticity $\omega_{r}\in[-4000,4000]$ colourmaps at
the intermediate radius $r_{\rm m}=(r_{\rm{i}}+r_{\rm{o}})/2\approx 8.05$. (a)
spt in the full orthogonal domain of periodic aspect ratio $\Lambda=31.4$,
with spectral resolution $(L,N,M)=(322,322,42)$ (See online movie 1). The
white dashed line indicates the tilt of the stripe and of the minimal narrow
domain that can capture spt (solid red enclosure). (b) spt in a narrow long
parallelogram domain $(n_{1},k_{1},n_{2},k_{2})=(1,0.2,0,4.5)$ with resolution
$(L,N,M)=(18,322,42)$ (See online movie 2). The solid cyan parallelogram
delimits the minimal domain required for self-sustained turbulent dynamics and
used throughout the paper. (c) Turbulent (top) and laminar (bottom) snapshots
at different time instants along the same simulation within a narrow and short
parallelogram domain $(n_{1},k_{1},n_{2},k_{2})=(10,2,0,4.5)$, using
$(L,N,M)=(18,32,42)$ spectral modes (See online movie 3).

The unwrapped domain is a rectangle in the $\theta$-$z$ plane as usually
employed in spt calculations, and its height-to-gap aspect ratio is
$\Lambda\approx 31.4$ (corresponding to an axial wavenumber
$k=2\pi/\Lambda=0.2$). Within the turbulent stripe, wavy vortices of
relatively short axial and azimuthal wavelength emerge. These small-scale flow
structures seem to decorrelate rather fast along the spiral slope direction.
This statistical property of spt suggests that a narrow periodic domain
winding horizontally around the full apparatus gap but sheared to preserve the
tilt of the spiral (dashed white line in panel a) would in principle suffice
to capture minimally its main topological and dynamical features.

In fact, a localised turbulent stripe arises when dns is performed in an
azimuthally long but axially narrow parallelogram domain such as shown in
figure 2b, whose comparative size is indicated in figure 2a as the enclosure
bound by a red line. The corresponding generalised wavenumbers for this domain
(see figure 1) were chosen as follows. First we note that the pitch of spt is
enforced by the dimensions of the full domain through the ratio
$k/n=2\pi/\Lambda$. A periodic parallelogram domain that can sustain spirals
of the same pitch must have the slope of two of its sides prescribed by
azimuthal and axial wavenumbers keeping the same ratio $k_{1}/n_{1}=k/n=0.2$.
Picking $n_{1}=1$, and therefore $(n_{1},k_{1})=(1,0.2)$, the parallelogram
domain winds exactly once around the apparatus circumference, with its two
opposing sides that are aligned with the slope of the spiral falling on the
same helical surface. The other pair of sides, defined by the values
$(n_{2},k_{2})$, can be chosen arbitrarily at this stage. For simplicity we
take $n_{2}=0$, thus forcing the parallelogram to extend horizontally and
align with the azimuthal coordinate $\xi\equiv\theta$. Notice that the axial
height of the domain is then simply specified as $2\pi/k_{2}$. Inspecting the
typical size of flow structures in figure 2a, shows that $k_{2}=4.5$ is an
adequate choice for the narrow domain size to fit exactly one single streak
most of the time.

We note here that in the full domain, the symmetry of the system in the axial
direction allows for a turbulent spiral band structure with an opposite tilt.
However, in the narrow but long parallelogram domain, the symmetry is broken.
The choice of the wavenumber ratio $k_{1}/n_{1}=0.2$ makes the realisation of
the opposite spiral impossible, which could be computed using instead the
wavenumber ratio $k_{1}/n_{1}=-0.2$.

Finally, if the focus is to be placed exclusively on the small-scale flow
structures that exist within the core of spt rather than on the large-scale
interactions that result in intermittency and laminar-turbulent interfaces,
but want to keep at the same time the potential influence of the
characteristic tilt of the spiral pattern, the long domain of figure 2b can be
further reduced by shortening it along the azimuthal direction into the small
parallelogram enclosed by the solid cyan/blue line. The two snapshots of
figure 2c illustrate the dynamics in the small parallelogram domain. The flow
is temporally chaotic and alternates rather quiescent, spatially-coherent and
smoothly evolving laminar phases (bottom), with sudden bursts of violent,
spatially-disordered and rapidly fluctuating (i.e. spatiotemporally chaotic)
turbulent transients (top). The former periods are seemingly representative of
the laminar regions of spt, while the latter display characterisitic features
typical of the vortical flow structures found in the core of turbulent
spirals. Here we have chosen to shorten the coiling length of the long domain
by setting an azimuthal wavenumber $n_{1}=10$ following Deguchi et al. (2014),
such that the first pair of wave numbers is replaced by
$(n_{1},k_{1})=(10,2.0)$ and a periodical tiling of the parallelogram along
the coil fits exactly 10 times.

It is precisely this minimal domain of figure 2c, with $k_{2}=4.5$, that we
employ in studying the ecs in §4, although some continuation in the $k_{2}$
parameter has also been performed to ellucidate the complex bifurcation
scenario of drw. The dns results shown in figure 2b,c will be discussed in §5.

##  4 Bifurcation scenarios of (relative) equilibria

Figure 3 summarises the collection of ecs that we have managed to capture in
the small parallelogram domain described above (see figure 2c), including
stationary, drifting and rotating-drifting wave solutions.

---
Figure 3: Bifurcation diagram of ecs in the small box of figure 2c. Shown are
solution branches of stationary (tvf, blue) and drifting (dtvf, red) Taylor
vortex flow, as well as of drifting-rotating waves (drw, black). White bullets
denote saddle-node (sn), pitchfork (p), Hopf (h), azimuthal invariance
breaking (ibθ) and axial-modulational subharmonic (shz) bifurcation points.
The number of unstable real eignenvalues (${\rm R}^{+}$) and complex conjugate
pairs (${\rm C}^{+}$) of some solutions along the branches are given in
parentheses.

Solution branches have been tracked via Newton-Krylov arclength continuation.
The spectral resolutions used lie within the range
$(L,N,M)\in[8,16]\times[8,16]\times[24,50]$, always ensuring that no further
increase resulted in qualitative or significant quantitative variation as to
flow properties or bifurcation scenarios.

Of particular interest in figure 3 is a family of three-dimensional
travelling-wave solutions that have non-zero phase speeds in both the
azimuthal and axial directions and will henceforth be referred to as Drifting
Rotating Waves (drw, black). Continuation reveals that drw branches are
remarkably subcritical, with the corresponding saddle-node points (sn1 and
sn3) located at inner Reynolds numbers ${\mathrm{R}_{\rm i}}^{\rm
SN_{1}}=391.5$ and ${\mathrm{R}_{\rm i}}^{\rm SN_{3}}=372.9$ preceding by far
the linear instability of the base flow (ccf) at ${\mathrm{R}_{\rm
i}}=447.35$. We will see later that, by tuning the geometrical parameters that
define the domain size and shape, the drw solution can be shown to essentially
bifurcate from Taylor Vortex Flow (tvf), which in turn bifurcates from ccf, a
multi-stage bifurcation scenario that was already reported by Deguchi et al.
(2014). In our particular choice of domain, however, the bifurcation details
are more involved.

The tvf solution branch (dark blue curve) bifurcates subcritically from ccf
for $k_{2}=9$ following an axisymmetric centrifugal instability at
${\mathrm{R}_{\rm i}}=484.7$. The branch evolves for a short ${\mathrm{R}_{\rm
i}}$-range and turns in a saddle-node at ${\mathrm{R}_{\rm i}}=470.8$ before
undergoing an axial-modulational subharmonic bifurcation (shz) at
${\mathrm{R}_{\rm i}}=483.8$, whence another branch of tvf (light blue), of
fundamental wave number $k_{2}=4.5$, is issued. Systematical Arnoldi linear
stability analysis along the tvf branches shows that the $k_{2}=9$ solution
stabilises briefly across the saddle-node, and that this stability is
transfered to the $k_{2}=4.5$ solution at the supercritical shz point. The
$k_{2}=4.5$ tvf solution, which bifurcates at ${\mathrm{R}_{\rm i}}=508.4$
from ccf at the other end, exhibits a pichfork bifurcation (p) half way along
the branch (${\mathrm{R}_{\rm i}}^{\rm P}=492.1$) that breaks the Z2 axial
reflection symmetry and produces a pair of mutually-symmetric branches of
axially-Drifting Taylor Vortex Flow solutions (dtvf, red). The dtvf branch
soon undergoes an azimuthal invariance breaking bifurcation (ibθ) that
generates the highly subcritical drw (black). The ibθ bifurcation is formally
a Hopf bifurcation, as a complex conjugate pair of eigenvalues crosses into
the positive-real-half of the complex plane, but it introduces no dynamics
other than the trivial drift along the group orbit corresponding to axial
shifts.

The flow topology of the $k_{2}=9$ and $k_{2}=4.5$ tvf solutions is
illustrated in figure 4 at points shz and p, respectively, through a couple of
azimuthal vorticity ($\omega_{\theta}$) iso-surfaces.

tvf
at shz |  tvf
at p
---|---
dtvf
at ibθ |  drw
at sn1
Figure 4: Distinct types of ecs at specific points labelled in figure 3. Shown
are positive (yellow) and negative (blue) isosurfaces of azimuthal vorticity
at $\omega_{\theta}=\pm 100$ for tvf at shz, tvf at p and dtvf at ibθ, and at
$\omega_{\theta}=\pm 600$ for drw at sn1.

They both possess a reflectional symmetry in the axial direction that blocks
all possibility of axial drift, which, in combination with azimuthal
invariance, keeps them stationary. Although the $k_{2}=4.5$ tvf preserves the
vortical arrangement of the $k_{2}=9$ solution, contiguous vortex pairs are no
longer identical as a consequence of the two-fold axial modulation enacted by
the shz bifurcation. The reflectional symmetry is nonetheless preserved and
the solutions remain stationary.

The symmetry is finally broken at the pitchfork point p, such that the
resulting solution drifts axially, as expected from bifurcation theory in the
presence of symmetries (Chossat & Iooss, 1994). Azimuthal invariance is
preserved and the flow structure is still very much alike that of Taylor
vortices, which accounts for its being dubbed Drifting-tvf (dtvf). The final
breaking of the mirror symmetry is clear from the $\omega_{\theta}$
isosurfaces of figure 4 for dtvf.

All these axisymmetric flows are characterised by the presence of vortical
structures in the vicinity of the inner cylinder, the reason being that the
centrifugal instability is confined to $r<r_{\rm n}=\sqrt{-B/A}$, according to
Rayleigh’s stability condition. Here $r_{\rm n}$ is the nodal radius, where
$A$ and $B$ are the constants in the base flow (6). On the other hand, the
three-dimensional structures of the drw, also depicted in figure 4 at sn1,
exhibit remarkably different properties that will be discussed in the next
section.

A second, seemingly unrelated branch of drw is laid out in figure 3 that
bifurcates in a saddle-node sn3 at ${\mathrm{R}_{\rm i}}^{\rm SN_{3}}=372.9$.
The two apparently disconnected families of drw are, in point of fact, one and
the same, as evinced by continuation in the $k_{2}$ parameter in figure 5.

(a) |
---|---
(b) |  | (c) |
Figure 5: Bifurcation diagram of ecs for (a) $k_{2}=5$, (b) $k_{2}=5.1$, and
(c) $k_{2}=4$. Colour coding and labels as for figure 3.

Increasing the axial wavelength to $k_{2}=5$ (figure 5a) the lower-torque
branch of the sn3-related wave approaches the higher-torque branch of the
sn1-related solution. At slightly higher $k_{2}$ the two branches collide in a
straightforward codimension-2 double-zero bifurcation point, where each splits
in two separate segments that are spliced in a new arrangement. As a result,
two mutually-facing saddle-node points arise, sn2 and sn4, which are
illustrated for $k_{2}=5.1$ in figure 5b. This type of structural instability
of pairs of saddle-nodes with respect to small changes in the size of the
domain has also been identified in other shear flows (Deguchi & Nagata, 2011;
Mellibovsky & Meseguer, 2015; Ayats et al., 2020b). The new branch containing
sn4 becomes disconnected, while the other one exhibits a triplet of saddle-
nodes and connects directly to $k_{2}=5.1$ tvf at the azimuthal invariance
breaking point ibθ. The transfer of ibθ from the dtvf to the tvf branch, which
follows a codimension-2 Hopf-pitchfork bifurcation, has in fact already been
effected at $k_{2}=5$, as is clear from the inset of figure 5a. In fact, the
subharmonic tvf branch (light blue) has outdone, in terms of subcriticality,
the tvf branch from which it bifurcates (dark blue). And to complicate things
further, the branch has also bent to include a second saddle-node that moves
fast towards higher ${\mathrm{R}_{\rm i}}$ as $k_{2}$ is increased.

Another interesting phenomenon also occurs when $k_{2}$ of drw is reduced from
4.5. The branch associated with sn1 disconnects from dtvf and becomes an isola
as illustrated in figure 5c for $k_{2}=4$ (black line). In addition, the
Reynolds number at which sn3 occurs increases, making it disappear from the
range of the diagram.

We conclude from the above parametric study that although the topological
arrangement of solution branches is rather sensitive to small changes in
$k_{2}$, their amply subcritical nature is a robust feature of drw solutions.

##  5 Structure and stability of DRW, and their role in the dynamics

We shall see presently some evidence that drw solutions appear to play a
central role in organising the flow dynamics both in the subcritical and
supercritical regimes of counter-rotating Taylor-Couette flow. In §5.1 we
discuss the particular flow structure of drw and relate it to the ssp. We then
analyse in §5.2 the stability of the waves and their contribution to the
formation of a chaotic set in the subcritical regime. Finally, we inspect in
§5.3 the early supercritical regime in order to gauge the part drw solutions
play in driving supercritical turbulent dynamics.

###  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
${\mathrm{R}_{\rm i}}=450$. Figure 6 provides the three-dimensional
characterisation of the flow structure of drw at the representative points
labeled a, b, and c in figure 3.

| (a) | (b)
---|---|---
drw
at a |  |
drw
at b |  |
drw
at c |  |
Figure 6: Three-dimensional flow structure of drw solutions at subcritical
${\mathrm{R}_{\rm i}}=450$, corresponding to points a, b and c of figure 3.
Shown are positive (yellow) and negative (blue) isosurfaces of (a)
perturbation azimuthal velocity $v$, and (b) 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).

Of the solutions ensuing saddle node sn1, point c is chosen to represent the
lower-torque branch. Two isosurfaces of the azimuthal perturbation velocity
($v$) shown in the top row of figure 6a reveal the presence of a low-speed
streak close to the inner cylinder. For the higher-torque branch (point b),
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 sn3. For the solutions related to sn3, only point a
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 6b 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 quasi-azimuthally-
invariant streak, that might be held responsible for the self-sustaiment 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_{\rm i}<r<r_{\rm 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, 1996; Deguchi,
2016). The streak field of drw must therefore be supported by the Reynolds
stresses induced by the wave-like vortex layer, rather than centrifugal
effects.

The comparison of both fields, shown together in figure 7, clearly discloses
the similarity of their interaction with the ssp that is typical of parallel
shear flows.

| (a) | (b) | (c)
---|---|---|---
drw
at a |  |  |
drw
at b |  |  |
drw
at c |  |  |
Figure 7: Comparison of (a) the streak field as signified by azimuthally-
averaged distribution of total azimuthal velocity $\langle
V\rangle_{\theta}/r\in[-1200,450]$, and the wave field, represented by
azimuthal vorticity $\omega_{\theta}\in[-1000,1000]$, at (b) $\theta=0$ and
(c) $\theta=\pi/n_{1}=0.314$, for drw solutions labeled a, b and c in figure
3. Lines indicate the location of the critical layer (solid) and the nodal
radius $r_{\rm n}$ (dashed).

Here we conjecture that the generation of the fluctuating field might likely
be explained by the inviscid linear instability of the streak field, where the
critical layer holds the key to the validation of this hypothesis. The
critical layer is formally a singularity of the inviscid problem and
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., 2007; Hall & Sherwin, 2010; Deguchi & Hall, 2015). 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, 1955; Drazin & Reid,
1981). In defining streaks and critical layers in the Taylor-Couette problem,
we note that subtle differences arise from the parallel flow cases. 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 streak field must be defined as $\langle
V\rangle_{\theta}/r$, using the azimuthal average
$\langle\bullet\rangle_{\theta}$ of the total (base plus perturbation)
velocity $V=v_{b}+v$. The black line in figure 7a indicates the critical layer
thus computed on top of the streak field. The azimuthal cross-sections of
$\omega_{\theta}$ colour maps at $\theta=\\{0,\pi/n_{1}\\}$ (panels b and c)
unequivocally show that the vortex sheet amplitude is strongest precisely
around the critical layer.

The magnitude of $\omega_{\theta}$ 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 and provides a physical explanation as to why the connections
between (d)tvf and drw occur in the way discussed in §4.

###  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
sn1 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 Taylor-
Couette flow. 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., 2014) 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 sn1 is lost in a
Hopf bifurcation at ${\mathrm{R}_{\rm i}}^{\rm 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, hence the acronym p-drw. The Poincaré-Newton-Krylov method
described in § 2 has been deployed to perform natural continuation of the
p-drw solution branch. Figure 8a signifies the oscillation amplitude of the
p-drw as the area (shaded blue region) bounded by the maximum and minimum
inner torque $\tau_{\rm i}$ along a full cycle (black dots).

(a) | (b)
---|---
|
Figure 8: Onset of p-drw solutions. (a) The bifurcation diagram near the
saddle-node sn1 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_{\rm{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 (gray dashed
curve). (b) The $(\tau_{\rm{o}},\tau_{\rm{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 ${\mathrm{R}_{\rm i}}$.

The solution amplitude obeys a square-root law of the form $A_{\tau_{\rm
i}}=\tau_{\rm i}^{\rm max}-\tau_{\rm i}^{\rm min}\sim({\mathrm{R}_{\rm
i}}-{\mathrm{R}_{\rm i}}^{\rm H})^{1/2}$, as revealed by the least-squares fit
(gray 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,
p-drw solutions are stable at onset and could therefore have been computed by
mere time-stepping. The evolution of p-drw with ${\mathrm{R}_{\rm i}}$ is more
clearly illustrated by the outer-vs-inner torque
$(\tau_{\rm{o}},\tau_{\rm{i}})$ phase-map projections of figure 8b. With each
increase in ${\mathrm{R}_{\rm i}}$, the lower-torque (crosses) and higher-
torque (circles) drw solutions become further apart in phase space. A small
limit cycle (p-drw) clearly orbits the higher-toque drw at ${\mathrm{R}_{\rm
i}}=393$, and its amplitude grows as ${\mathrm{R}_{\rm i}}$ is further
increased.

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)\mathbf{a}^{\rm DRW}$, with $|\gamma|\ll 1$. Figure 9a shows the
$(\tau_{\rm{o}},\tau_{\rm{i}})$-phase map projections of a couple runs at
${\mathrm{R}_{\rm i}}=392$.

(a) | (b)
---|---
|
Figure 9: Phase-map projections of trajectories issued from the unstable
lower-torque drwLT solution on the $(\tau_{\rm o},\tau_{\rm i})$-plane at
inner Reynolds number (a) ${\mathrm{R}_{\rm i}}=392$ and (b) ${\mathrm{R}_{\rm
i}}=394$. The represented solutions are ccf (square), lower-torque (drwLT,
cross) and higher-torque (drwHT, circle) drw solutions, and p-drw (blue line).
The red/black lines depict dns trajectories starting from initial conditions
taken by scaling drwLT by $(1+\gamma)$, with $\gamma=\pm 10^{-4}$,
respectively.

As we have already noted, the lower-torque drwLT 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
drwHT solution (red line). The lower-torque drw solution is therefore acting
as an edge state (Itano & Toh, 2001; Skufca et al., 2006), separating the
basins of attraction of ccf on one side and the higher-torque drwHT solution
on the other. The dynamics remain qualitatively the same as illustrated here
within the parameter range $391.52<{\mathrm{R}_{\rm i}}<392.85$, while the
upper branch solution preserves its linear stability.

Slightly above the Hopf bifurcation (${\mathrm{R}_{\rm i}}>{\mathrm{R}_{\rm
i}}^{\rm H}=392.85$), the higher-torque drwHT solution has become unstable and
the stable p-drw has popped up into existence, as illustrated at
${\mathrm{R}_{\rm i}}=394$ by figure 9b. The dynamics are qualitatively
unaltered as regards ccf or drwLT, which remains an edge state, but the phase-
map trajectory that previously led to drwHT 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
drwHT, which remains a focus across the Hopf bifurcation, but switches from
stable to unstable. Sufficiently close to sn1 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, 2004).

At even higher Reynolds number ${\mathrm{R}_{\rm 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 p2-drw solution (black line) is shown, along with the unstable p-drw
limit cycle, at the same ${\mathrm{R}_{\rm i}}=395.5$ (dashed red) in figure
10a.

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

The period doubling that occurs just short of ${\mathrm{R}_{\rm i}}=395.5$ is
followed by a number of ensuing period-doublings upon further increasing
${\mathrm{R}_{\rm i}}$ that eventually lead to the chaotically-modulated drw
(ch-drw) of figure 10b at ${\mathrm{R}_{\rm 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, 2012; Lustro et al., 2019). The stable state at ${\mathrm{R}_{\rm
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.

The stability analysis reflected in figure 3 reveals at least two more Hopf
bifurcations of the already unstable drw as ${\mathrm{R}_{\rm 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 ${\mathrm{R}_{\rm i}}$
(Krygier et al., 2021).

###  5.3 Turbulent dynamics at ${\mathrm{R}_{\rm i}}=600$

As shown in figure 11 at ${\mathrm{R}_{\rm 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.

(a) | (b)
---|---
|
Figure 11:  Spatiotemporally chaotic dynamics at ${\mathrm{R}_{\rm i}}=600$ in
the short parallelogram domain of figure 2c
$(n_{1},k_{1},n_{2},k_{2})=(10,2,0,4.5)$. (a) Phase-map projections on the
$(\tau_{\rm o},\tau_{\rm i})$ plane. Gray bullet and error bars denote mean
and root-mean-square of the fluctuations, respectively. (b) Time series of
kinetic energy $E$. The gray line indicates mean kinetic energy. The colour
circles/squares and lines indicate drw, tvf and dtvf solutions labeled in
figure 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 ${\mathrm{R}_{\rm 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
sn1 (drwD, blue) has a perturbation kinetic energy $E$ somewhat larger than
the mean for the statistically steady turbulent state. The other two
solutions, resulting from sn3 (higher-torque solution drwE in green and lower-
torque solution drwF 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.

Interesting new phenomena occur when the computational domain is extended in
the azimuthal direction to cover the whole circumference of the annulus.
Figure 12 shows dns results in the long computational domain of figure 2b
characterised by $(n_{1},k_{1},n_{2},k_{2})=(1,2,0,4.5)$.

(a) | (b)
---|---
|  |
---
Figure 12: Transient stages of the development of spt in the long
parallelogram domain of figure 2b $(n_{1},k_{1},n_{2},k_{2})=(10,2,0,4.5)$ at
${\mathrm{R}_{\rm i}}=600$, starting from a 10-fold azimuthal replication of
drwF (red bullet in figure 11a). (a) Time series of kinetic energy $E$. (b)
Selected instantaneous snapshots of radial vorticity
$\omega_{r}\in[-4000,4000]$ fields at the mid gap $r_{m}$. The panels T1–T5
correspond to the points indicated in (a), while T6 is taken beyond all
transients for $t\gg 1$. See online movie 4, presented as supplementary
material.

The flow field has been initialised with a ten-fold azimuthal replication of
the drwF solution and thence left to evolve freely. The simulation deviates
very slowly from drwF in the beginning, such that snapshot T1, taken at
$t=0.2$, is still indistinguishable from drwF. Soon after, instability leads
to a sudden surge of kinetic energy. Snapshot T2 is taken at the peak and, as
clear from panel 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.

The flow structure of the instantaneous snapshot of the statistically steady
state T6 is examined in greater detail in figure 13.

(a)
---
| (b) | (c) | (d) | (e)
---|---|---|---
|  |  |
Figure 13: An instantaneous snapshot of the statistically steady spt state in
the long domain at ${\mathrm{R}_{\rm i}}=600$. (a) A close up of snapshot f of
figure 12b. (b)-(e) Azimuthal vorticity $\omega_{\theta}$ distribution in the
$r$–$z$ plane at azimuthal cross-sections $\theta_{1}=0.5$, $\theta_{2}=2.3$,
$\theta_{3}=3.1$, and $\theta_{4}=4.3$ (black vertical lines in panel a).
Colour ranges between $\omega_{\theta}\in[-1500,1500]$. From left to right,
the black curves in (c)-(e) correspond to selected isocontours with (c)
$V/r=-20,-80$, (d) $V/r=0,-30,-100$, (e) $V/r=0$. See online movie 5,
presented as supplementary material, showing a cross-sectional sweep of the
instantaneous flow field along the full perimeter.

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\left(c_{\theta}-\langle
c_{\theta}\rangle_{t}\right)^{2}\rangle_{t}}=15.2}$, while the local wave
propagation speed is strongly location-dependent. This behaviour of the flow
field is typical of localised solutions in shear flows (see for example Avila
et al., 2013; Brand & Gibson, 2014; Zammert & Eckhardt, 2014; Mellibovsky &
Meseguer, 2015), and may be mathematically justified by a Wentzel-Kramers-
Brillouin-type approach (see Bender & Orszag (1999), for example).

The local appearance of ssp-related flow structures in streamwise-extended
domains has been reported by Deguchi & Hall (2014) in plane Couette flow using
ecs. 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
similar phenomenon is observed in figure 13 despite the flow being turbulent.
The azimuthal vorticity colour map on a $r$–$z$ azimuthal cross-section near
the trailing edge of the turbulent stripe (figure 13e) 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 figures 13c and d. In the core of the turbulent band, several strong
vortex sheets are visible (figure 13d), while only those close to the outer
cylinder seem to persist toward the leading edge (figure 13c). The presence of
vortex sheets at various radial locations is consistent with the multiplicity
of drw solutions existing at this value of ${\mathrm{R}_{\rm i}}$, which
suggests that each one of them might participate locally in the deployment of
the ssp.

The basic flow is linearly unstable at ${\mathrm{R}_{\rm i}}=600$, which would
in principle be at odds with the observation of a laminar flow region in
figures 13a and b. 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 13b 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_{\rm
cu}=(r_{\rm i}+r_{\rm n})/2=7.705$, in the midst of the centrifugally unstable
region (figure 14a). A spiral-like flow topology clearly arises away from the
turbulent stripe.

(a) | (c)
---|---
|  |
---

(b)

Figure 14: The morphology of spt within the centrifugally unstable region of
the annulus gap. Same as figure 2, but at $r=r_{\rm cu}=(r_{\rm i}+r_{\rm
n})/2=7.705<r_{\rm n}$. Colour ranges between $\omega_{r}\in[-4000,4000]$. See
online supplementary materials: (a) movie 1, (b) movie 2 and (c) movie 3.

Despite the dislocations and irregularities, the vortical structures in the
laminar region compare favourably with exact spiral solutions (Meseguer et
al., 2009b) of the right axial-azimuthal wave-numbers (and therefore tilt).
The base ccf is unstable to a wide range of axial and azimuthal wavenumbers
for $(\eta,{\mathrm{R}_{\rm i}},{\mathrm{R}_{\rm 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 15 shows a selection of centrifugally-
driven states that fit exactly in either of the three domains considered.

(a) | (c)
---|---
|  |
---

(b)

Figure 15: Examples of centrifugally-driven exact solutions in the three
domains for $({\mathrm{R}_{\rm o}},{\mathrm{R}_{\rm i}})=(-1200,600)$,
visualised through $\omega_{r}\in[-4000,4000]$ colourmaps at $r=r_{\rm
cu}=(r_{\rm i}+r_{\rm n})/2=7.705<r_{\rm n}$. (a) $(n,k)=(5,7.4)$ spi. (b)
$(n,k)=(6,5.7)$ spi. (c) $k=9$ tvf (top) and $k=4.5$ dtvf (bottom).

In particular, figure 15a 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 wave numbers 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 14b, for which dislocated spiral patterns of wave numbers compatible
with those of the exact spi of figure 15b, 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 14c), however, quasi-axisymmetric tvf-like
structures arise during the laminar stages of the time evolution (bottom)
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 15a) and narrow but long parallelogram (figure 15b)
domains, such solutions do not exist in the short and narrow parallelogram
domain (figure 15c). 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 toward them.
Axisymmetric modes with $k=4.5$, 9 and 13.5 are the three only linear
instabilities of ccf in the small parallelogram domain at $({\mathrm{R}_{\rm
o}},{\mathrm{R}_{\rm 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 quasi-modal 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. Online
movie 3, presented as supplementary material, clearly shows how the energy-
growing quasi-axisymmetric 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 (1996) in Taylor-Couette flow 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_{\textrm{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 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 Taylor-Couette flow. This 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., 2016). 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 Taylor-
Couette flow when considering a suitably small parallelogram domain. However,
detection of bistability at ${\mathrm{R}_{\rm 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 (1996), 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 (1996) 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 spatio-temporal 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 to 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 to 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 to domain shape, size and discretisation. Here we have
used the kinetic energy ($E$), inner ($\tau_{\rm i}$) and outer ($\tau_{\rm
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, 1979, 1981) and the mean and
root-mean-square of the fluctuation component (standard deviation) computed.
Using stationary bootstrapping (Politis & Romano, 1994), 95% confidence
intervals for the two statistics estimators have been determined with
automatic block size optimisation (Politis & White, 2004; Patton et al.,
2009). 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.

case | full domain | long narrow domain
---|---|---
$(n_{1},k_{1},n_{2},k_{2})$ | $(1,0,0,0.2)$ | $(1,0.2,0,4.5)$
$(L,N,M)$ | $(322,322,42)$ | $(18,322,42)$
estimator | $\langle\bullet\rangle$ | $\sigma_{\bullet}=\sqrt{\langle\left(\bullet-\langle\bullet\rangle\right)^{2}\rangle}$ | $\langle\bullet\rangle$ | $\sigma_{\bullet}=\sqrt{\langle\left(\bullet-\langle\bullet\rangle\right)^{2}\rangle}$
$E$ | $16390\pm 40$ | $570\pm 30$ | $16960\pm 90$ | $520\pm 20$
$\tau_{\rm i}$ | $0.658\pm 0.001$ | $0.0195\pm 0.0007$ | $0.689\pm 0.003$ | $0.0182\pm 0.0005$
$\tau_{\rm o}$ | $0.658\pm 0.002$ | $0.0232\pm 0.0009$ | $0.690\pm 0.004$ | $0.0209\pm 0.0006$
$c_{\theta}$ | $-34.18\pm 0.08$ | $3.47\pm 0.07$ | $-35.4\pm 0.2$ | $3.2\pm 0.2$
Table 1: Statistics of several time signals (mean $\langle\bullet\rangle$ and
root-mean-square of the fluctuation component
$\sigma_{\bullet}=\sqrt{\langle\left(\bullet-\langle\bullet\rangle\right)^{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 2b provides reasonable quantitative agreement, wider parallelogram
domains would be required to further enhance accuracy, particularly so
regarding second-order statistics.

##  6 Conclusions

We have identified numerically three-dimensional nonlinear-wave solution
branches that spread over a wide parameter-space region of both subcritical
and supercritical counter-rotating Taylor-Couette flow. The states have been
computed in suitable minimal parallelogram-shaped domains incorporating the
oblique pseudo-invariance of the spiral turbulence (spt) regime, as observed
both experimentally (Coles & Van Atta, 1967) and numerically (Meseguer et al.,
2009b; Dong, 2009; Dong & Zheng, 2011). We have undertaken here to ascertain
their dynamical relevance as precursors of spiral turbulence beyond reasonable
doubt.

The waves, which mainly rotate in the azimuthal direction but possess also a
mild axial drift, spring from saddle-node bifurcations that considerably
anticipate the primary instability of the base laminar flow. The lower-torque-
waves branch is indirectly connected to circular Couette flow through an
intricate succession of intervening secondary solutions progressively breaking
the symmetries of the base flow in a complex bifurcation scenario. The
sequence involves the subharmonic bifurcation of slightly subcritical Taylor
vortices, an axial symmetry breaking that introduces the axial drift, and a
rupture of the azimuthal invariance that sets the wave into rotation. The
actual arrangement of solution branches and the way they connect with one
another is highly dependent on the width of the parallelogram domain, as also
are the stability properties of the solutions themselves.

The strong subcriticality of these waves implies that some mechanism, other
than centrifugal instability, must be accountable for their self-sustainment.
Axisymmetric solutions bifurcated directly from the base flow cannot be
continued far into the subcritical region, which suggests that the three-
dimensional streamwise-dependent component of the drifting-rotating waves
holds the key to self-sustainment. The origin of the three-dimensional wave
might be explained by the inviscid instability of the streak, understood as
the axisymmetric component of the azimuthal velocity field. Compelling
evidence that this is the underlying mechanism at play is provided by the
strong wave-like vortex sheet that concentrates around the critical layer,
where the inviscid problem becomes singular. It is precisely this vortex sheet
that allegedly drives the axisymmetric roll-streak field through the action of
the Reynolds stresses, thereby closing the interaction feed-back loop between
the axisymmetric and three-dimensional components. This is reminiscent of the
roll-streak-wave cycle that is characteristic of parallel shear flows
(Waleffe, 1997; Wang et al., 2007; Hall & Sherwin, 2010), except that the
streak and the roll mutually interact through both the lift-up and Coriolis
effects, and not only the former. The latter coupling term might probably
affect the formal asymptotic structure, a detailed analysis of which is beyond
the scope of this paper.

Arnoldi stability analysis and direct numerical time integration of the
Navier-Stokes equations reveal that the lower-torque branch of drifting-
rotating waves acts as an edge state separating the basin of attraction of the
base circular Couette flow from that of nontrivial nonlinear states. Moreover,
in a neighbourhood of the saddle-node and for suitably chosen streamwise width
of the parallelogram domain, the higher-torque branch happens to be initially
linearly stable. This situation is very convenient in that it can be
fruitfully exploited to one’s advantage for the detailed analysis of the onset
of chaotic motion as the governing parameter is increased further (Mellibovsky
& Eckhardt, 2012; Kreilos & Eckhardt, 2012; Lustro et al., 2019).

In the case under scrutiny, the route to chaos commences with a supercritical
Hopf bifurcation of the drifting-rotating wave that issues a branch of stable
time-periodic solutions. We have found evidence that this relative periodic
orbit undergoes a period-doubling cascade that eventually engenders a chaotic
set. As the supercritical regime is approached, the drifting-rotating wave
undergoes a succession of additional linear instabilities whence a number of
bifurcation cascades akin to the one just mentioned are expected to ensue. It
is the concurrent action of several chaotic sets thus generated that would
ultimately beget turbulent dynamics. Notably, drifting-rotating waves, albeit
unstable, can be continued all the way up to the supercritical region of
parameter space where spt is ubiquitous. A comparison of these solutions with
turbulence as computed in the small parallelogram domain for the same values
of the parameters suggests that the former are fairly descriptive of the self-
sustained process (ssp) that drives the complex dynamics of the latter.

Simulations within a narrow parallelogram domain, extended in the azimuthal
direction to wrap completely around the apparatus gap, have been deployed to
explore the relationship between the localised turbulent stripe and the ssp.
In the elongated domain, streamwise inhomogeneity naturally arises in the form
of banded localisation of turbulence. The final statistically-steady state
captures surprisingly well the features of spt as computed in the full
orthogonal domain, both qualitatively and, to a large extent, also
quantitatively. The active part of inhomogeneous turbulence within spt is
characterised by the ssp, where streaks and vortex layers, analogous to those
we report in small periodic domains, are also observed. Meanwhile, in the
apparently quiescent flow region, spiral vortices arise in the close proximity
of the inner cylinder following the centrifugal instability to which this
region is subject.

Unlike what happens for subcritical parallel shear flows featuring laminar-
turbulent patterns, spt is observed in supercritical counter-rotating Taylor-
Couette flow despite the instability of the base laminar flow. In the case of
parallel flows, the laminar-turbulent banded pattern alternates patches of the
laminar and the turbulent states, both stable in suitably chosen small
computational domains. Present results show that this seems to not be the case
for mildly supercritical counter-rotating Taylor-Couette flow. Neither the
laminar nor the turbulent states can be considered permanent, and therefore
stable, in minimal flow units. In particular, no stable state has been found
that may account for the spiral-like waves that appear in the laminar region
of spiral turbulence due to the centrifugal instability. Instead, the dynamics
in small domains are invariably chaotic, and alternate in time short
relaxation periods onto some sort of centrifugally-driven axisymmetric laminar
state resembling Taylor vortex flow, and sudden bursts into short-lived shear-
driven turbulence. To some extent, this behaviour is similar to that
identified by Coughlin & Marcus (1996) in axially narrow orthogonal domains
for a slightly wider gap, except that the laminar transients consisted of a
dislocated spirals pattern instead of axisymmetric tvf-like vortices. This
might require reconsidering the mechanisms that underly laminar-turbulent
pattern formation when centrifugal effects are at stake. To give but one
example, directed percolation theory, which has been lately employed in
explaining subcritical turbulence in parallel flows (Lemoult et al., 2016;
Sano & Tamai, 2016; Chantry et al., 2017), does not apply to the supercritical
regime. It might therefore be worthwhile examining how the theory is affected
by the centrifugal instability.

The exact spiral solutions found by Meseguer et al. (2009a) are highly
unstable within the small periodic domain, so that tvf-type axisymmetric modes
take the lead in governing the nearly-laminar phases of the chaotic dynamics.
Nevertheless, it might still be the case that spiral-like structures govern
the dynamics within the laminar regions in azimuthally-extended domains, and
that the patterns observed therein, including wave dislocations, arise from
the competition of a continuum of such spirals of varying wavelength and tilt.
The key point is that, since neither full-fledged turbulence nor stable
centrifugally-driven exact solutions are to be observed in the supercritical
regime within the small periodic domain, some large-scale effect must be
responsible for the permanent sustainment of turbulence in localised bands and
for the stabilisation of laminar spiral patterns within the laminar patches of
spt. Large-scale pattern formation such as that observed in these laminar
regions is sometimes explained by amplitude equations derived from weakly
nonlinear theory, together with some external noise term (Prigent et al.,
2002; Berghout et al., 2020). Whether this artificial noise term can somehow
be replaced by the highly nonlinear and autonomous ssp structures observed
here remains unclear.

Another question that arises naturally is whether the long but narrow domain
can accomodate exact solutions that may be capable of explaining the actual
topology of laminar-turbulent patterns. Exact localised solutions resulting
from subharmonic bifurcation of periodic wave-train solutions sharing with the
laminar-turbulent stripes their main large-scale features, have been found for
several parallel shear flows (Reetz et al., 2019; Paranjape et al., 2020).
These solutions typically develop from modulational instability of some
subcritical wave to long-wavelength perturbations (Chantry et al., 2014;
Mellibovsky & Meseguer, 2015), and their nonlinear evolution leads to the
formation of localised wave fronts that connect a patch of the nontrivial
state with the base flow at either side. Besides direct stability analysis and
branch continuation, edge-tracking has emerged as a powerful method to compute
boundary states. Unfortunately, this approach requires stable laminar and
turbulent states, something that is missing in supercritical Taylor-Couette
flow as investigated here.

Whether one employs a deterministic or a statistical approach, the
parallelogram-shaped domain we have developed here opens the path to a more
detailed yet affordable analysis of a wide variety of problems in Taylor-
Couette flow that could hitherto only be addressed for parallel shear flows
and, therefore, in the absence of centrifugal effects.

[Supplementary data]Supplementary material and movies are available at
https://doi.org/10.1017/jfm.2019…

[Funding] KD’s research was supported by Australian Research Council Discovery
Early Career Researcher Award (DE170100171). BW, RA, FM and AM research was
supported by the Spanish Ministerio de Economía y Competitivdad (grant numbers
FIS2016-77849-R and FIS2017-85794-P) and Ministerio de Ciencia e Innovación
(grant number PID2020-114043GB-I00), and the Generalitat de Catalunya (grant
2017-SGR-785). BW’s research was also supported by the Chinese Scholarship
Council (grant CSC No. 201806440152). Where no specific funding has been
provided for research, please provide the following statement: ”This research
received no specific grant from any funding agency, commercial or not-for-
profit sectors.”

[Declaration of interests] The authors report no conflict of interest.

[Data availability statement] The data that support the findings of this study
may be obtained from the corresponding author upon reasonable request.

## References

  * Andereck et al. (1986) Andereck, C. D., Liu, S. S. & Swinney, H. L. 1986 Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech. 164, 155–.
  * Avila et al. (2013) Avila, M., Mellibovsky, F., Roland, N. & Hof, B. 2013 Streamwise-localized solutions at the onset of turbulence in pipe flow. Phys. Rev. Lett. 110, 224502.
  * Ayats et al. (2020a) Ayats, R., Deguchi, K., Mellibovsky, F. & Meseguer, A. 2020a Fully nonlinear mode competition in magnetised Taylor-Couette flow. J. Fluid Mech. 897, A14.
  * Ayats et al. (2020b) Ayats, R., Meseguer, A. & Mellibovsky, F. 2020b Symmetry-breaking waves and space-time modulation mechanisms in two-dimensional plane Poiseuille flow. Phys. Rev. Fluids 5, 094401.
  * Barkley & Tuckerman (2005) Barkley, D. & Tuckerman, L. S. 2005 Computational study of turbulent laminar patterns in Couette flow. Phys. Rev. Lett. 94, 014502.
  * Bender & Orszag (1999) Bender, C. M. & Orszag, S. A. 1999 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory, 1st edn. Springer.
  * Berghout et al. (2020) Berghout, P., Dingemans, R. J., Zhu, X., Verzicco, R., Stevens, R. J. A. M., van Saarloos, W. & Lohse, D. 2020 Direct numerical simulations of spiral taylor–couette turbulence. J. Fluid Mech. 887, A18.
  * Boberg & Brosa (1988) Boberg, L. & Brosa, U. 1988 Onset of turbulence in a pipe. Z. Naturforsch. 43a, 697–726.
  * Brand & Gibson (2014) Brand, E. & Gibson, J. F. 2014 A doubly localized equilibrium solution of plane Couette flow. J. Fluid Mech. 750.
  * Canuto et al. (2007) Canuto, C. G., Hussaini, M. Y., Quarteroni, A. & Zang, T. A. 2007 Spectral methods: evolution to complex geometries and applications to fluid dynamics. Springer.
  * Canuto et al. (2010) Canuto, C. G., Hussaini, M. Y., Quarteroni, A. & Zang, T. A. 2010 Spectral methods: fundamentals in single domains, 2nd edn. Springer.
  * Chantry et al. (2017) Chantry, M., Tuckerman, L. S. & Barkley, D. 2017 Universal continuous transition to turbulence in a planar shear flow. J. Fluid Mech. 824.
  * Chantry et al. (2014) Chantry, M., Willis, A. P. & Kerswell, R. R. 2014 Genesis of streamwise-localized solutions from globally periodic traveling waves in pipe flow. Phys. Rev. Lett. 112, 164501\.
  * Chossat & Iooss (1994) Chossat, P. & Iooss, G. 1994 The Couette-Taylor Problem. Springer-Verlag.
  * Clever & Busse (1992) Clever, R. M. & Busse, F. H. 1992 Three-dimensional convection in a horizontal fluid layer subjected to a constant shear. Journal of Fluid Mechanics 234, 511.
  * Clever & Busse (1997) Clever, R. M. & Busse, F. H. 1997 Tertiary and quaternary solutions for plane couette flow. J. Fluid Mech. 344, 137–153.
  * Coles & Van Atta (1967) Coles, D. & Van Atta, C. W. 1967 Digital experiment in spiral turbulence. Phys. Fluids Suppl. 121, S120–.
  * Coughlin & Marcus (1996) Coughlin, K. & Marcus, P.S. 1996 Turbulent bursts in couette-taylor flow. Phys. Rev. Lett. 77(11), 2214–2217.
  * Deguchi (2016) Deguchi, K. 2016 The rapid-rotation limit of the neutral curve for Taylor-Couette flow. J. Fluid Mech. 808, R2.
  * Deguchi & Altmeyer (2013) Deguchi, K. & Altmeyer, S. 2013 Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow. Phys. Rev. E 87, 043017.
  * Deguchi & Hall (2014) Deguchi, K. & Hall, P. 2014 The high-reynolds-number asymptotic development of nonlinear equilibrium states in plane Couette flow. J. Fluid Mech. 750, 99–112.
  * Deguchi & Hall (2015) Deguchi, K & Hall, P. 2015 Asymptotic descriptions of oblique coherent structures in shear flows. J. Fluid Mech. 782, 356–367.
  * Deguchi et al. (2014) Deguchi, K., Meseguer, A. & Mellibovsky, F. 2014 Subcritical equilibria in Taylor-Couette flow. Phys. Rev. Lett. 112, 184502.
  * Deguchi & Nagata (2011) Deguchi, K. & Nagata, M. 2011 Bifurcations and instabilities in sliding Couette flow. J. Fluid Mech. 678, 156–178.
  * Dickey & Fuller (1979) Dickey, D. A. & Fuller, W. A. 1979 Distribution of the estimators for autoregressive time-series with a unit root. J. Am. Stat. Assoc. 74 (366), 427–431.
  * Dickey & Fuller (1981) Dickey, D. A. & Fuller, W. A. 1981 Likelihood ratio statistics for autoregressive time-series with a unit-root. Econometrica 49 (4), 1057–1072.
  * Dong (2009) Dong, S. 2009 Evidence for internal structures of spiral turbulence. Phys. Rev. E 80, 067301.
  * Dong & Zheng (2011) Dong, S. & Zheng, X. 2011 Direct numerical simulation of spiral turbulence. J. Fluid Mech. 668, 150–.
  * Drazin & Reid (1981) Drazin, P. G. & Reid, W. H. 1981 Hydrodynamic Stability. Cambridge University Press.
  * Duguet et al. (2010) Duguet, Y., Schlatter, P. & Henningson, D. S. 2010 Formation of turbulent patterns near the onset of transition in plane Couette flow. J. Fluid Mech. 650, 119.
  * Eckhardt et al. (2007) Eckhardt, B., Schneider, T. M., Hof, B. & Westerweel, J. 2007 Turbulence transition in pipe flow. Ann. Rev. Fluid Mech. 39, 447–468.
  * Esser & Grossmann (1996) Esser, A. & Grossmann, S. 1996 Analytic expression for Taylor-Couette stability boundary. Phys. Fluids 8, 1814–1819.
  * Faisst & Eckhardt (2003) Faisst, H. & Eckhardt, B. 2003 Travelling waves in pipe flow. Phys. Rev. Lett. 91, 224502.
  * Feynman (1964) Feynman, R. P. 1964 Lecture Notes in Physics, , vol. 2. Reading, MA: Addison-Wesley.
  * Graham & Floryan (2021) Graham, M. D & Floryan, D. 2021 Exact coherent states and the nonlinear dynamics of wall-bounded turbulent flows. Ann. Rev. Fluid Mech. 53, 227–253.
  * Grossmann (2000) Grossmann, S. 2000 The onset of shear flow turbulence. Rev. Mod. Phys. 72, 603–618.
  * Hall & Sherwin (2010) Hall, P. & Sherwin, S. 2010 Streamwise vortices in shear flows: harbingers of transition and the skeleton of coherent structures. J. Fluid Mech. 661, 178–205.
  * Hamilton et al. (1995) Hamilton, J. M., Kim, J. & Waleffe, F. 1995 Regeneration mechanisms of near-wall turbulence. J. Fluid Mech. 287, 317–348.
  * Itano & Toh (2001) Itano, T. & Toh, S. 2001 The dynamics of bursting process in wall turbulence. J. Phys. Soc. Japan 70 (3), 703–716.
  * Kawahara et al. (2012) Kawahara, G., Uhlmann, M. & Van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. Ann. Rev. Fluid Mech. 44, 203–225.
  * Kelley (2003) Kelley, C.T. 2003 Solving nonlinear equations with Newton’s method. Philadelphia: SIAM.
  * Kerswell (2005) Kerswell, R. R. 2005 Recent progress in understanding the transition to turbulence in a pipe. Nonlinearity 18, R17–R44.
  * Kreilos & Eckhardt (2012) Kreilos, T. & Eckhardt, B. 2012 Periodic orbits near onset of chaos in plane Couette flow. CHAOS 22, 047505\.
  * Krygier et al. (2021) Krygier, M. C., Pughe-Sanford, J. L. & Grigoriev, R. O. 2021 Exact coherent structures and shadowing in turbulent taylor–couette flow. J. Fluid Mech. 923, A7.
  * Kuznetsov (2004) Kuznetsov, Y. A. 2004 Elements of Applied Bifurcation Theory, 3rd Ed.. Springer-Verlag.
  * Lemoult et al. (2016) Lemoult, G., Shi, L., Avila, K., Jalikop, S. V., Avila, M. & Hof, B. 2016 Directed percolation phase transition to sustained turbulence in couette flow. Nature Physics 12, 254–258.
  * Lin (1955) Lin, C. 1955 The theory of hydrodynamic stability. Cambridge University Press.
  * Lustro et al. (2019) Lustro, J.R.T., Kawahara, G., van Veen, L., Shimizu, M. & Kokubu, H. 2019 The onset of transient turbulence in minimal plane Couette flow. J. Fluid Mech. 862, R2.
  * Mellibovsky & Eckhardt (2011) Mellibovsky, F. & Eckhardt, B. 2011 Takens–Bogdanov bifurcation of travelling-wave solutions in pipe flow. J. Fluid Mech. 670, 96–129.
  * Mellibovsky & Eckhardt (2012) Mellibovsky, F. & Eckhardt, B. 2012 From travelling waves to mild chaos: a supercritical bifurcation cascade in pipe flow. J. Fluid Mech. 709, 149–190.
  * Mellibovsky & Meseguer (2015) Mellibovsky, F. & Meseguer, A. 2015 A mechanism for streamwise localisation of nonlinear waves in shear flows. J. Fluid Mech. 779, R1.
  * Meseguer (2020) Meseguer, A. 2020 Fundamentals of Numerical Mathematics for Physicists and Engineers. Hoboken, N.J.: John Wiley & Sons.
  * Meseguer et al. (2007) Meseguer, A., Avila, M., Mellibovsky, F. & Marques, F. 2007 Solenoidal spectral formulations for the computation of secondary flows in cylindrical and annular geometries. Eur. Phys. J. Special Topics 146, 249–259.
  * Meseguer et al. (2009a) Meseguer, A., Mellibovsky, F., Avila, M. & Marques, F. 2009a Families of subcritical spirals in highly counter-rotating Taylor-Couette flow. Phys. Rev. E 79, 036309.
  * Meseguer et al. (2009b) Meseguer, A., Mellibovsky, F., Avila, M. & Marques, F. 2009b Instability mechanisms and transition scenarios of spiral turbulence in Taylor-Couette flow. Phys. Rev. E 80, 046315.
  * Moser et al. (1983) Moser, R. D., Moin, P. & Leonard, A. 1983 A spectral numerical method for the Navier-Stokes equations with applications to Taylor-Couette flow. J. Comp. Phys. 52, 524–544.
  * Nagata (1990) Nagata, M. 1990 Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity. J. Fluid Mech. 217, 519–527.
  * Paranjape et al. (2020) Paranjape, C. S., Duguet, Y. & Hof, B. 2020 Oblique stripe solutions of channel flow. J. Fluid Mech. 897, A7.
  * Patton et al. (2009) Patton, A., Politis, D. N. & White, H. 2009 Correction to “automatic block-length selection for the dependent bootstrap” by d. politis and h. white. Econom. Rev. 28 (4), 372–375.
  * Politis & Romano (1994) Politis, D. N. & Romano, J. P. 1994 The stationary bootstrap. J. Am. Stat. Assoc. 89 (428), 1303–1313.
  * Politis & White (2004) Politis, D. N. & White, H. 2004 Automatic block-length selection for the dependent bootstrap. Econom. Rev. 23 (1), 53–70.
  * Prigent et al. (2002) Prigent, A., Grégoire, G., Chaté, H., Dauchot, O. & van Saarloos, W. 2002 Finite-wavelength modulation within turbulent shear flows. Phys. Rev. Lett. 89, 014501.
  * Reetz et al. (2019) Reetz, F., Kreilos, T. & Schneider, T. M. 2019 Exact invariant solution reveals the origin of self-organized oblique turbulent-laminar stripes. Nat. Commun. 10, 2277.
  * Sano & Tamai (2016) Sano, M. & Tamai, K. 2016 A universal transition to turbulence in channel flow. Nat. Phys. 12(3), 249–253.
  * Schneider et al. (2008) Schneider, T. M., Gibson, J. F., Lagha, M., De Lillo, F. & Eckhardt, B. 2008 Laminar-turbulent boundary in plane Couette flow. Phys. Rev. E 78, 037301.
  * Shi et al. (2013) Shi, L., Avila, M. & Hof, B. 2013 Scale invariance at the onset of turbulence in Couette flow. Phys. Rev. Lett. 110, 204502.
  * Skufca et al. (2006) Skufca, J. D., Yorke, J. A. & Eckhardt, B. 2006 Edge of chaos in a parallel shear flow. Phys. Rev. Lett. 96, 174101.
  * Trefethen & Bau (1997) Trefethen, L. N. & Bau, D. 1997 Numerical Linear Algebra. Philadelphia: SIAM.
  * Tuckerman et al. (2020) Tuckerman, L. S., Chantry, M. & Barkley, D. 2020 Patterns in wall-bounded shear flows. Ann. Rev. Fluid Mech. 52, 010719–060221.
  * Tuckerman et al. (2014) Tuckerman, L. S., Kreilos, T., Schrobsdorff, H., Schneider, T. M. & Gibson, J. F. 2014 Turbulent-laminar patterns in plane Poiseuille flow. Phys. Fluids 26, 114103\.
  * Waleffe (1997) Waleffe, F. 1997 On a self-sustaining process in shear flows. Phys. Fluids 9, 883–900.
  * Waleffe (2003) Waleffe, F. 2003 Homotopy of exact coherent structures in plane shear flows. Phys. Fluids 15(6), 1517–1534.
  * Wang et al. (2007) Wang, J., Gibson, J. & Waleffe, F. 2007 Lower branch coherent states in shear flows: Transition and control. Phys. Rev. Lett. 98, 204501.
  * Wedin & Kerswell (2004) Wedin, H. & Kerswell, R. R. 2004 Exact coherent structures in pipe flow: travelling wave solutions. J. Fluid Mech. 508, 333–371.
  * Zammert & Eckhardt (2014) Zammert, S. & Eckhardt, B. 2014 Streamwise and doubly-localised periodic orbits in plane Poiseuille flow. J. Fluid Mech. 761, 348–359.