[2] $\int$[0pt][l][]$-$$\int_{#1}^{#2}$

# On high Taylor number Taylor vortices

Kengo Deguchi  School of Mathematics, Monash University, VIC 3800, Australia

###### Abstract

Axisymmetric steady solutions of Taylor-Couette flow at high Taylor numbers
are studied numerically and theoretically. As the axial period of the solution
shortens from about one gap length, the Nusselt number goes through two peaks
before returning to laminar flow. In this process, the asymptotic nature of
the solution changes in four stages, as revealed by the asymptotic analysis.
When the aspect ratio of the roll cell is about unity, the solution captures
quantitatively the characteristics of the classical turbulence regime.
Theoretically, the Nusselt number of the solution is proportional to the
quarter power of the Taylor number. The maximised Nusselt number obtained by
shortening the axial period can reach the experimental value around the onset
of the ultimate turbulence regime, although at higher Taylor numbers the
theoretical predictions eventually underestimate the experimental values. An
important consequence of the asymptotic analyses is that the mean angular
momentum should become uniform in the core region unless the axial wavelength
is too short. The theoretical scaling laws deduced for the steady solutions
can be carried over to Rayleigh-Bénard convection.

##  1 Introduction

Taylor-Couette flow (TCF) is perhaps among the most studied flows in fluid
mechanics. In the 100 years since Taylor’s monumental work (Taylor 1923), it
has provided an excellent testing ground for theoretical, experimental and
numerical studies of rotating shear flows. How shear and Coriolis forces alter
flow characteristics is important in various applications, and TCF was
designed so that they can easily be adjusted by changing the rotation speed of
the inner and outer cylinders. Researchers have long been fascinated by the
numerous metastable flow patterns observed in the relatively low Taylor number
regime (e.g. Andereck et al. 1986). On the other hand, it was only a decade
ago that the study of high Taylor number flows became active. Great efforts
were made to investigate the nature of turbulence in the parameter space by
means of high Taylor number experiments (e.g. Paoletti & Lathrop 2011; Dennis
et al. 2011; van Gils et al. 2011, 2012; Huisman et al. 2012) and direct
numerical simulations (e.g. Ostilla et al. 2013; Brauckmann & Eckhardt 2013,
2017; Ostilla-Mónico et al. 2014ab). As summarised in the review paper by
Grossmann et al. (2016), and in fact seen in the pioneering experiments by
Lathrop et al. (1992), Lewis & Swinney (1999), fully developed turbulence has
a surprisingly clean asymptotic character, while there seems to be no
definitive Navier-Stokes-based theory to explain it.

This study aims to reveal the asymptotic properties of steady axisymmetric
solutions at high Taylor numbers and to compare them with the experimental and
numerical results. The analysis of such solutions, known as Taylor vortex
solutions, goes back to the weakly nonlinear analysis, for example, by Davey
et al. (1962). With modern computational power, it is possible to calculate
solutions up to the Taylor number used in the experiments. Of course, the use
of Newton’s method is essential as the solution is unstable in the high Taylor
number regime.

The idea that an unstable solution with a relatively simple structure with
respect to time can capture some characteristics of turbulence is not as
absurd as one might think. It is well-known in dynamical systems theory that
chaotic dynamics can be approximated by a sufficiently large number of
periodic orbits (see Cvitanović et al. 2012, for example). For moderate
Reynolds number shear flows, it was repeatedly reported that a good
approximation of the turbulent dynamics can be obtained with a small number of
periodic orbits (Kawahara et al. 2012; Willis et al. 2013; Krygier et al.
2021). In phase space, these periodic orbits usually appear with stationary or
travelling wave solutions in their vicinity. The advantage of focusing on
simple solutions is that their asymptotic nature may be justified
theoretically by mathematical analyses. Over the past decade, it has been
established that the matched asymptotic expansion is a powerful tool in
understanding the behaviour of steady or travelling wave solutions in shear
flows (Hall & Sherwin 2010, Deguchi et al. 2013, Deguchi & Hall 2014ab,
Deguchi 2015, Dempsey et al. 2017, Deguchi & Walton 2018). Therefore if we are
allowed to assume that there is a simple solution that roughly captures the
scaling properties of turbulence within the vast phase space, then there is
hope for a logical explanation of the scaling from first principles.

In the high Taylor number numerical and experimental studies, the parameter
dependence of angular momentum transport was a major focus. The driving force
behind those studies was the ‘analogy’ between turbulent Rayleigh-Bénard
convection (RBC) and TCF.  This analogy was well known at least in the 1960s
and has risen and fallen throughout the history of turbulence research
(Bradshaw 1969, Dubrulle & Hersant 2002). The recent studies have been
influenced by Eckhardt et al. (2007) who argued that the phenomenology of RBC
turbulence proposed by Grossmann & Lohse (2000) can be applied to TCF as well.
Shortly later the aforementioned high Taylor number TCF experiments confirmed
that the scaling of the Nusselt number $Nu$ is similar to that observed for
RBC by He et al. (2012) ($Nu$ for TCF is defined as the torque on the cylinder
wall normalised by its laminar value).  It should be remarked that despite the
analogy that has been believed, this similarity in the ultimate scaling is
actually not at all obvious. As pointed out, for example, by Chandrasekhar
(1961), Robinson (1967), Veronis (1970), and Lezius & Johnston (1976), for the
two flows to be equivalent they must be at least axisymmetric. Moreover, the
exact equivalence of the two flows requires an infinitesimally narrow cylinder
gap, co-rotating cylinders, and a Prandtl number of unity.

It is an interesting question, then, to forget the latter three conditions and
ask whether the analogy in the sense of the $Nu$ scaling holds for an
axisymmetric Taylor vortex and a two-dimensional roll cell. The equations
governing both flows are not the same, but they certainly have a similar
structure. Many researchers have theoretically studied the large Rayleigh
number nature of roll cells in RBC over the years; see Pillow (1952), Robinson
(1967), Wesseling (1969), Chini & Cox (2009), Hepworth (2014) for constant
Prandtl number flows, Roberts (1979), Jimenez & Zufiria (1987), and Vynnycky &
Masuda (2013) for asymptotically large Prandtl number flows. Waleffe et al.
(2015) and Sondak et al. (2015) shortened the wavelength of the RBC roll cell
and found that there is a special wavelength at which the Nusselt number
reaches a maximum value. Interestingly, this optimised Nusselt number is close
to that obtained experimentally, although the Rayleigh number they used is
much lower than the one used by He et al. (2012). More recently, Wen et al.
(2022) continued the same solution branch to higher Rayleigh numbers and
claimed that the maximum Nusselt number corresponds to the so-called classical
scaling, where the Nusselt number is proportional to the Rayleigh number to
the power of one-third (Malkus et al. 1954; Priestley 1954; Grossmann & Lohse
2000; Kawano et al. 2021). The study by Kooloth et al. (2021) confirmed that
roll cells with different wavelengths play important roles in two-
dimensionally restricted RBC turbulence. This paper is motivated by all the
above RBC studies.

In the classical turbulence regime of TCF, where Taylor vortices are robustly
observed in experiments, the symmetry restriction of the flow may not be
necessary for a good agreement between the steady solution and turbulence, and
at least a better agreement than in RBC can be expected. The situation is
different in the ultimate turbulence regime, where eddies of different sizes
and wavelengths are present, making the structure far more complex than
classical turbulence. However, as long as the cylinders are not strongly
counter-rotating, the experimental observations indicate that vortices of
about the scale of the gap are still present, suggesting that Taylor vortices
may play some role in the dynamics.

It should also be noted that previous theoretical studies have shown that
analytical approximations can be obtained for vortices with extremely short
wavelengths driven by thermal or Coriolis forces (Hall & Lakin 1988; Bassom &
Hall 1989; Bassom & Blennerhassett 1992; Denier 1992; Blennerhassett & Bassom
1994). In this type of asymptotic theory, the mean flow varies by a finite
amount from the base flow and is therefore called a strongly nonlinear theory.
However, the scaling of momentum and heat transport of the nonlinear state is
not so different from that of the basic flow, and in this sense, the character
of fully developed turbulence is not well captured. Attempts have been made to
construct asymptotic solutions with larger amplitudes, but a complete
understanding is still lacking. This paper proposes a solution to this long-
standing problem.

In the next section, we begin by formulating our problem. Comparisons of the
Taylor vortex solutions with previous experiments and simulations are then
carried out in section 3. In section 4, a matched asymptotic expansion
analysis is performed for the case of Taylor vortices with an aspect ratio
about unity. We will see in section 5 that the asymptotic structure of the
solution dramatically changes when the axial period becomes asymptotically
short. In section 6, how the short-period vortices develop from the laminar
solution is investigated in detail using a matched asymptotic expansion.
Finally, in section 7, the main findings are summarised and discussed.

##  2 Formulation of the problem

Taylor-Couette flow can be described by the Navier-Stokes equations in the
cylindrical coordinates $(r,\theta,z)$. If the flow is axisymmetric, the
governing equations are written as

$\displaystyle Du-r^{-1}v^{2}=-\partial_{r}p+\triangle u-r^{-2}u,$ (1a)
$\displaystyle Dv+r^{-1}uv=\triangle v-r^{-2}v,$ (1b) $\displaystyle
Dw=-\partial_{z}p+\triangle w,$ (1c) $\displaystyle
r^{-1}\partial_{r}(ru)+\partial_{z}w=0.$ (1d)

The operators $D$ and $\triangle$ are defined as
$D=\partial_{t}+u\partial_{r}+w\partial_{z}$ and
$\triangle=\partial_{r}^{2}+r^{-1}\partial_{r}+\partial_{z}^{2}.$ The length
and velocity scales are chosen so that the cylinder gap is unity and the no-
slip conditions on the cylinder walls are described as

$\displaystyle(u,v,w)=(0,R_{o},0)\qquad\text{at}\qquad r=r_{o},$ (2)
$\displaystyle(u,v,w)=(0,R_{i},0)\qquad\text{at}\qquad r=r_{i},$ (3)

using the Reynolds numbers associated with the rotation of the inner and outer
cylinders, $R_{i}$ and $R_{o}$. Note that our non-dimensionalisation implies
that using the radius ratio $\eta=r_{i}/r_{o}<1$ the inner and outer radii are
specified as

$\displaystyle r_{i}=\frac{\eta}{1-\eta},\qquad r_{o}=\frac{1}{1-\eta},$ (4)

respectively. The circular Couette flow solution is written as

$\displaystyle(u,v,w)=(0,v_{b},0),\qquad v_{b}(r)=\frac{R_{o}-\eta
R_{i}}{1+\eta}r+\frac{\eta^{-1}R_{i}-R_{o}}{1+\eta}\frac{r_{i}^{2}}{r}.$ (5)

For other nontrivial solutions of (1) periodicity is imposed in the interval
$z\in[0,2\pi/k]$ using the axial wavenumber $k$. The momentum equations
(1a)-(1c) can be simplified as

$\displaystyle r^{-1}D(rv)=\triangle v-\frac{v}{r^{2}},\qquad
rD\left(\frac{\omega}{r}\right)=\triangle\omega-\frac{\omega}{r^{2}}+\frac{\partial_{z}v^{2}}{r},$
(6a) using the azimuthal vorticity
$\displaystyle\omega=\partial_{z}u-\partial_{r}w=-\\{\partial_{r}(r^{-1}\partial_{r}\Psi)+r^{-1}\partial_{z}^{2}\Psi\\}$
(6b) and the Stokes streamfunction $\Psi$. The roll-cell velocity can be
reconstructed as $u=-r^{-1}\partial_{z}\Psi$ and $w=r^{-1}\partial_{r}\Psi$.

Steady solutions of the above system can be computed without regarding their
stability by using the numerical code used in our previous studies (Deguchi &
Altmeyer 2013; Deguchi et al. 2014). The code is based on the Newton-Raphson
method applied to the Chebyshev-Fourier discretised system. To calculate the
Taylor vortex with an aspect ratio of about unity, we typically used up to
250th Chebyshev polynomials and 250th Fourier harmonics. This spatial
resolution is more than sufficient for $Ta=O(10^{10})$. We shall see that when
$k$ is large, we can reach higher Taylor numbers. In this case, the highest
degree of Chebyshev polynomials was increased to 450 to fully resolve the very
thin near-wall boundary layer structures.

Instead of the two Reynolds numbers, the majority of high Taylor number TCF
studies summarised in Grossmann et al. (2016) used the Taylor number $Ta$ and
the rotation rate $a$. They are easily found by the standard parameters as

$\displaystyle Ta^{1/2}=\frac{(1+\eta)^{3}}{8\eta^{2}}(R_{i}-\eta
R_{o}),\qquad a=-\eta\frac{R_{o}}{R_{i}}.$ (7)

The former parameter is similar to the shear Reynolds number used in Dubrulle
et al. (2005) and is zero for the rigid body rotation case (their second
parameter, the rotation number, is a function of $a$ and $\eta$).  The Nusselt
number is defined by

$\displaystyle
Nu=\left.\frac{\partial_{r}(r^{-1}\overline{v})}{\partial_{r}(r^{-1}v_{b})}\right|_{r=r_{i}}=\left.\frac{\partial_{r}(r^{-1}\overline{v})}{\partial_{r}(r^{-1}v_{b})}\right|_{r=r_{o}},$
(8)

where

$\displaystyle\overline{v}(r)=\frac{k}{2\pi}\int^{2\pi/k}_{0}v(r,z)dz$ (9)

is the mean azimuthal velocity.

In the limit $\eta\rightarrow 1$, the gap becomes much narrower than the
cylinder radius, and the local flow can be represented in Cartesian
coordinates. As noted by Deguchi (2016) there are several variations on the
narrow gap limit, two of which are relevant to this paper. One is of course
the rotating plane Couette flow, where the system is in perfect agreement with
RBC if the Prandtl number is unity (see Chandrasekhar 1961, Robinson 1967,
Veronis 1970, and Lezius & Johnston 1976, Brauckmann & Eckhardt 2017, for
example). The other utilises an argument similar to the derivation of the
Görtler vortex (Hall 1983), and is used, for example, by Denier (1992). For
completeness, the difference between the two limits is highlighted in Appendix
A.

##  3 Comparison of the steady solutions and the experiments

Figure 1:  The variation of Nusselt number $Nu$ with respect to Taylor number
$Ta$. The outer cylinder is fixed ($a=0$). $\eta=5/7\approx 0.714$. The red
solid curve is the Taylor-vortex solution branch with the fixed wavenumber
$k=3$. The red crosses are the same solutions, but the wavenumber is optimised
to maximise $Nu$. The blue circles are the three-dimensional direct numerical
simulations by Ostilla et al. (2013) and Ostilla-Mónico (2014a). The cyan
triangles and squares are the experiments by Lewis & Swinney (1999) and Dennis
et al. (2011), respectively. The simulation and experimental results are time-
averaged data. The best theoretical upper bound known to date has the
asymptotic form $Nu=0.0075\,Ta^{1/2}$ according to Ding & Marensi (2019).

Here we use the radius ratio $\eta=5/7$ and fix the outer cylinder ($a=0$) to
compute the Taylor vortex solutions. Significant deviations from the narrow
gap approximation can be observed at this radius ratio. That parameter choice
is frequently used in experiments and numerical simulations (see section 3 of
Grossman et al. 2016) and is hence convenient for the comparison. In
experiments, it has been confirmed that the effect of endwalls on torque/mean
flow measurements is negligible (see van Gils et al. (2012), for example).

The red curve in figure 1 shows the Taylor vortex solution calculated with the
fixed axial wavenumber $k$. According to the simulations, the Taylor cells are
relatively robust even in the numerically generated classical turbulent flows,
with a cell aspect ratio of approximately unity. More specifically, direct
numerical simulations (DNS) by Ostilla et al. (2013) imposed the axial
periodicity $2\pi$ and typically observed three vortex pairs (i.e. $k=3$
modes). This motivated the choice of $k=3$ for our Taylor vortex computation.

The solution branch bifurcates from the circular Couette flow at $Ta\approx
10^{4}$, and as the Taylor number increases, it loses stability with respect
to three-dimensional perturbations. The onset of turbulence is $Ta\approx
3\times 10^{6}$. The blue circles in the figure are the DNS by Ostilla et al.
(2013) and Ostilla-Mónico et al. (2014a). Despite the relatively short axial
periodicity imposed, the simulations are known to agree with the experimental
results (the squares and triangles in figure 1). Both the experiments and
simulations indicate the existence of a transition point $Ta\approx 10^{8}$ at
which the behaviour of $Nu$ changes. The turbulence below/above the transition
point is referred to as the classical/ultimate regime. The Nusselt number of
the ultimate turbulence has the scaling $Nu\propto Ta^{\beta}$ with the
exponent $\beta\approx 0.38$, which is also seen in the RBC experiments (He et
al. 2012). The empirical asymptotic prediction $Nu=0.009\,Ta^{0.38}$ sits well
below the theoretical upper bound of $Nu$ derived by Ding & Marensi (2019).
The exponent 1/2 is typical for upper bounds using the energy method. Similar
asymptotic behaviour of $Nu$ but with a logarithmic correction has been
proposed using the phenomenology of the log law of turbulent boundary layers
(see Grossmann et al. 2016). It is, however not known whether such scaling can
be clearly observed in experiments.

(a) | (b)
---|---
|
Figure 2:  The mean angular velocity $q$ defined in (10) for $\eta=5/7,a=0$.
(a) The classical turbulence regime $Ta=9.52\times 10^{6}$. The red solid
curve is the Taylor vortex solution. The blue dot-dashed curve is the time-
averaged DNS result by Ostilla et al. (2013). (b) The ultimate turbulence
regime. The blue dot-dashed curve is the time-averaged DNS result by Ostilla-
Mónico et al. (2014b) ($Ta=10^{10}$). The other curves are the Taylor vortex
solutions shown in figure 3 (b)-(d) ($Ta=9.75\times 10^{9}$).

(a)

(b) $k=3$ (c) $k=27.4$ (d) $k=112$

Figure 3:  Axial wavenumber dependence of the Taylor vortex solutions. The
parameters are $\eta=5/7,R_{i}=8\times 10^{4},R_{o}=0$, which correspond to
$Ta=9.75\times 10^{9}$ in figure 1. (a) The bifurcation diagram. The blue
dotted curve is the Taylor vortex solution branch. This branch bifurcates from
the linear critical point of the circular Couette flow (‘L’ in the figure).
There is another linear critical point at very small $k$, but computing the
bifurcating solution branch at this Taylor number is difficult and hence it is
omitted. (b)-(d) Azimuthal velocity $v$ at the selected points in figure 3a.
The colour bar range is [0,80000]. All the solutions posses the reflection
symmetry in $z$. (e) and (f) are enlarged views of parts of (c) and (d),
respectively. The colour bar range is changed to [10000,70000] in the enlarged
figures.

As long as the solution has an aspect ratio of about unity, it captures the
nature of classical turbulence surprisingly well. This is best illustrated by
a comparison of mean flows shown in figure 2a. Here

$\displaystyle
q=\frac{\overline{v}/r-R_{o}/r_{o}}{R_{i}/r_{i}-R_{o}/r_{o}}\in[0,1]$ (10)

is the shifted and normalised mean angular velocity (denoted by
$\langle\overline{\omega}\rangle_{z}$ in Ostilla et al. 2013). At sufficiently
high Taylor numbers, the boundary layer formation near the cylinder walls is
clearly visible. As seen in figure 1, the solution gives a reasonable estimate
of $Nu$ in the classical turbulence regime. The cross-section of the Taylor
vortex (see figure 3b) reveals that the boundary layer actually surrounds the
vortex core, where the azimuthal velocity appears to be rather uniform. The
dynamics of the boundary layers in the classical turbulence are known to be
somewhat quiescent (see Ostilla et al. 2013, for example), and their
qualitative structure is reminiscent of figure 2b. The literature summarised
in Grossman et al. (2016) often implies that the Prandtl-Blasius theory could
be applied to the boundary layer, and we shall see in section 4 that this is,
in fact, true for the asymptotic limit of the steady solution. Moreover, the
theory to be presented in section 4 yields the scaling of $Nu\propto
Ta^{0.25}$ and $Re_{w}\propto Ta^{0.5}$, which are well agreeing with the
turbulent observations.  Here $Re_{w}$ is the wind Reynolds number, i.e. the
typical roll cell circulation speed normalised by the viscous velocity scale
$\nu/d$ ($\nu$ is the kinematic viscosity of the fluid and $d$ is the cylinder
gap). The precise definition of the wind Reynolds number varies according to
the literature. Huisman et al. (2012) measured the standard deviation of the
radial velocity $u$, while Ostilla et al. (2013) computed the average of
$u^{2}+w^{2}$ and then square rooted it. Both definitions give the same
scaling, but with different prefactors.

In the ultimate turbulence, on the other hand, the situation appears to be
more intricate. The snapshots of turbulence are typically accompanied by
eddies of a vast scale, whereby the large-scale Taylor vortex is blurred in
the time-averaged field. In particular, small turbulent eddies appearing in
the boundary layer have been pointed out as a critical qualitative difference
between the two turbulent regimes separated by the transition point seen in
figure 1 (Dong 2007, Ostilla-Mónico et al. 2014a). Therefore solutions with
larger wavenumbers may play a more critical role in the turbulent dynamics.
Our Nusselt number results also support this speculation. As seen in figure 1,
the Taylor vortex with fixed $k$ underestimates experimental $Nu$ in the
ultimate regime. However, if $k$ is optimised to maximise $Nu$ (the crosses in
figure 1), it can reach the experimental values for $Ta\lesssim 10^{11}$.

Figure 3a shows how $Nu$ changes as the wavenumber $k$ of the Taylor vortex is
varied. The $Nu$ curve has two local maxima. The bimodal distribution of $Nu$
and the scaling of the two peaks are remarkably similar to those seen in the
RBC roll cell computation by Waleffe et al. (2015) and Sondak et al. (2015).
Based on their calculations at several Rayleigh numbers $Ra$, the first
(second) peak has the wavenumber scaling $k\propto Ra^{0.217}$ ($Ra^{0.256}$)
with the Nusselt number scaling exponent $\beta$ slightly larger (smaller)
than 0.31. We shall provide a theoretical rationale for the scaling in
sections 5 and 6.

For the ultimate turbulence regime, the approximation of the turbulent mean
flow by the $k=3$ state is slightly worse (figure 2b).  However, the mean flow
of course varies with $k$. As the value of $k$ is increased from 3, the mean
flow of the solution in the core region approaches a turbulent profile in the
vicinity of the second peak. Figures 3c,d show the flow field at each peak,
which is also examined in detail in the following two sections.

##  4 Taylor vortices with an $O(1)$ cell aspect ratio

In figure 2a, we saw that some kind of homogenisation occurs in the core
region of the Taylor vortex, and a thin boundary layer emerges around it. Such
a flow structure looks very much like the high Rayleigh number RBC roll cell,
which was studied intensively from the 1950’s to the 1970’s, for example, by
Pillow (1952), Robinson (1967), and Wesseling (1969). More recently, it was
reported that when the slip walls are imposed, the asymptotic solution can be
calculated semi-analytically and is in good agreement with the numerical
solutions (Chini & Cox 2009; Hepworth 2014). However, as noted by Robinson
(1967), for no-slip walls, the scaling and structure of the boundary layer
have to be modified, which makes analytical progress difficult. The situation
in the Taylor vortex is close to the latter problem, as we shall see below.

First, let us examine the core region. For RBC, the roll cell vorticity in
this region is known to become constant due to the Prandtl-Batchelor theorem,
which states that the homogenisation of the vorticity occurs in the region
where the streamlines are closed in two-dimensional inviscid flows (Prandtl
1904; Feynman & Lagerstrom 1956; Batchelor 1956). Moreover, the temperature
becomes constant as well because the temperature equation has a similar
structure to that of the vorticity equation (see Moore & Weiss 1973, for
example). To better see what precise physical quantities are homogenised in
the core of the Taylor vortex, it is desirable to choose large gaps. Figure 4
shows the results for $\eta=0.5$. The colour map shown in figure 4b clearly
indicates that it is actually angular momentum $rv$ that becomes uniform in
the region where the streamlines are closed. Also, figure 4c shows that it is
not the azimuthal vorticity $\omega$ that homogenises, but $\omega/r$.

(a) (b) (c)

Figure 4: The flow field of the Taylor vortex for $\eta=0.5,k=3$. The Reynolds
numbers used are $R_{i}=8\times 10^{4},R_{o}=0.25R_{i}$, corresponding to
$Ta=1.40\times 10^{10}$, $a=-1/8$. (a) The Stokes streamfunction $\Psi$. The
colour bar range is [-7500,7500]. (b) The angular momentum $rv$. The colour
bar range is [40000,80000]. (c) The modified azimuthal vorticity $\omega/r$.
The colour bar range is [-120000,120000]. This quantity is uniformly
distributed in the core to a value of about $\pm 37$% of the colour bar.

Building on the above observations, now we introduce the new variables
$\Gamma=Ta^{-1/2}rv$ and $\Omega=\omega/r$, whereby the governing equations
(6) are rewritten as

$\displaystyle\Psi_{r}\Gamma_{z}-\Psi_{z}\Gamma_{r}=(r^{3}(r^{-2}\Gamma)_{r})_{r}+r\Gamma_{zz},$
(11a)
$\displaystyle\Psi_{r}\Omega_{z}-\Psi_{z}\Omega_{r}=(r^{-1}(r^{2}\Omega)_{r})_{r}+{r\Omega}_{zz}+Ta\frac{2\Gamma\Gamma_{z}}{r^{3}},$
(11b) $\displaystyle\Omega=-r^{-1}\\{(r^{-1}\Psi_{r})_{r}+r^{-1}\Psi_{zz}\\}.$
(11c)

Here, the subscript $r$ or $z$ denotes a partial differentiation. The no-slip
conditions become

$\displaystyle\Psi=\Psi_{r}=0,\qquad\Gamma=\Gamma_{i}\equiv\frac{8\eta^{3}}{(1-\eta)(1+\eta)^{3}(1+a)}\qquad\text{at}\qquad
r=r_{i},$ (12a)
$\displaystyle\Psi=\Psi_{r}=0,\qquad\Gamma=\Gamma_{o}\equiv-\frac{8\eta
a}{(1-\eta)(1+\eta)^{3}(1+a)}\qquad\text{at}\qquad r=r_{o}.$ (12b)

If $\Gamma$ is considered as temperature, $\Omega$ as roll cell vorticity and
$Ta$ as Rayleigh number, one may notice that the structure of the equations is
very similar to that of the two-dimensional RBC. The last term in the right
hand side of (4.1b) corresponds to the Coriolis force when the RPCF limit is
taken, and it plays a role of buoyancy. (This term is not exactly the Coriolis
force unless the system is in the RPCF limit, but for simplicity we will call
it ‘Coriolis force’ hereafter; see Appendix A also.) This analogy in the sense
of the structure of the equations would explain why $\Gamma$ and $\Omega$
become constant in the Taylor vortex core, as seen in figure 4, at least on an
intuitive level.

Figure 5:  Sketch of the asymptotic states. In the blue shaded region
viscosity is not negligible. In the dotted region Coriolis force is at work.
(a) Taylor vortex with aspect ratio of order unity ($k=O(1)$). (b,c) The first
peak state ($k=O(Ta^{2/9})$). (b) is the close up of the near wall zone
enclosed by the red lines in (c).

In fact, following Batchelor (1956), a mathematical argument can be developed.
First, assuming the typical roll cell circulation strength $Re_{w}$ (the wind
Reynolds number) is asymptotically large, we expand
$\Psi=Re_{w}\Psi_{0}+\cdots,\Gamma=\Gamma_{0}+\cdots,\Omega=Re_{w}\Omega_{0}+\cdots$.
Then the leading order part of (11a) becomes
$\Psi_{0r}\Gamma_{0z}-\Psi_{0z}\Gamma_{0r}=0$, which suggests that the
function $\Gamma_{0}$ only depends on $\Psi_{0}$. Now, in the $r$-$z$ plane,
take a region $\mathcal{A}$ enclosed by a streamline $\mathcal{C}$ oriented
counterclockwise (thus $\Psi$ is constant along $\mathcal{C}$). Integrating
(11a) over $\mathcal{A}$, we obtain

$\displaystyle\int_{\mathcal{A}}\\{\Psi_{r}\Gamma_{z}-\Psi_{z}\Gamma_{r}\\}drdz=\int_{\mathcal{A}}\\{(r^{3}(r^{-2}\Gamma)_{r})_{r}+(r^{3}(r^{-2}\Gamma)_{z})_{z}\\}drdz.$
(13)

The left hand side vanishes because upon using the Stokes theorem it becomes

$\displaystyle\int_{\mathcal{A}}\\{(\Psi_{r}\Gamma)_{z}-(\Psi_{z}\Gamma)_{r}\\}drdz=-\oint_{\mathcal{C}}\Gamma\\{\Psi_{z}\mathbf{e}_{z}+\Psi_{r}\mathbf{e}_{r}\\}\cdot
d\mathbf{l}=0.$ (14)

The second equality comes from the fact that the gradient of $\Psi$ and the
line element on $\mathcal{C}$, $d\mathbf{l}$, are orthogonal on the
streamline. While the leading order part of the right hand side can be
transformed by again applying the Stokes theorem:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\oint_{\mathcal{C}}r^{3}\\{(r^{-2}\Gamma_{0})_{r}\mathbf{e}_{z}-(r^{-2}\Gamma_{0})_{z}\mathbf{e}_{r}\\}\cdot
d\mathbf{l}$ (15) $\displaystyle=$
$\displaystyle\frac{d\Gamma_{0}}{d\Psi_{0}}\oint_{\mathcal{C}}r^{2}\\{r^{-1}\Psi_{0r}\mathbf{e}_{z}-r^{-1}\Psi_{0z}\mathbf{e}_{r}\\}\cdot
d\mathbf{l}-2\Gamma_{0}\oint_{\mathcal{C}}\mathbf{e}_{z}\cdot d\mathbf{l}$
$\displaystyle=$
$\displaystyle\frac{d\Gamma_{0}}{d\Psi_{0}}\oint_{\mathcal{C}}r^{2}\\{u_{0}\mathbf{e}_{r}+w_{0}\mathbf{e}_{z}\\}\cdot
d\mathbf{l}.$

Here $(u_{0},w_{0})=(-r^{-1}\Psi_{0z},r^{-1}\Psi_{0r})$ is the leading order
part of the roll velocity. The integral in the last line should not vanish
because we are assuming the existence of a strong circulation due to the
swirling motion of the Taylor roll. Thus the equation implies
$\frac{d\Gamma_{0}}{d\Psi_{0}}=0$ at the value of $\Psi_{0}$ on
${\mathcal{C}}$ we choose. The above argument holds for any region enclosed by
a streamline, and hence the value of $\Gamma_{0}$ must be constant in the core
region. The argument above does not change when an arbitrary constant is added
to $\Gamma$. This implies that if power series asymptotic expansion of
$\Gamma$ is considered the homogenisation also occurs in all the higher-order
terms as long as they have a spatial scale of $O(1)$. Therefore for steady
solutions, the non-uniform component in $\Gamma$ is exponentially small.

Likewise we can show that the value of $\Omega_{0}$ is constant as well in the
core. By integrating (11b) over ${\mathcal{A}}$,

$\displaystyle\int_{\mathcal{A}}\\{\Psi_{r}\Omega_{z}-\Psi_{z}\Omega_{r}\\}drdz$
$\displaystyle=$
$\displaystyle\int_{\mathcal{A}}\\{(r^{-1}(r^{2}\Omega)_{r})_{r}+(r^{-1}(r^{2}\Omega)_{z})_{z}\\}drdz$
(16) $\displaystyle+\int_{\mathcal{A}}Ta\frac{2\Gamma\Gamma_{z}}{r^{3}}drdz,$

and of course the left hand side should vanish. Since we already know that
$\Gamma_{0}$ is a constant plus an exponentially small fluctuation, the last
term in the right hand side of (16) can be neglected. From (11b) we see that
$\Omega_{0}$ is a function of $\Psi_{0}$, and thus the leading order part of
integral (16)

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\oint_{\mathcal{C}}r^{-1}\\{(r^{2}\Omega_{0})_{r}\mathbf{e}_{z}-(r^{2}\Omega_{0})_{z}\mathbf{e}_{r})\\}\cdot
d\mathbf{l}$ (17) $\displaystyle=$
$\displaystyle\frac{d\Omega_{0}}{d\Psi_{0}}\oint_{\mathcal{C}}r^{2}\\{r^{-1}\Psi_{0r}\mathbf{e}_{z}-r^{-1}\Psi_{0z}\mathbf{e}_{r}\\}\cdot
d\mathbf{l}+2\Omega_{0}\oint_{\mathcal{C}}\mathbf{e}_{z}\cdot d\mathbf{l}$
$\displaystyle=$
$\displaystyle\frac{d\Omega_{0}}{d\Psi_{0}}\oint_{\mathcal{C}}r^{2}\\{u_{0}\mathbf{e}_{r}+w_{0}\mathbf{e}_{z}\\}\cdot
d\mathbf{l},$

yields $\frac{d\Omega_{0}}{d\Psi_{0}}=0$.

In the above argument, we have assumed that the swirling speed of the rolls is
sufficiently large, but to see exactly how large it is, we need to analyse the
boundary layer. Let us now take the region $\mathcal{A}$ as large as possible
and assume that a viscous boundary layer appears around its boundary
$\mathcal{C}$. In this core region (see figure 5a), we have the estimation
$Re_{w}=O(\Psi|_{c})=O(\Omega|_{c})$, where the subscript $c$ indicates that
the physical quantities are measured in the core. At the roll cell perimeter
$\mathcal{C}$ without loss of generality, we can set $\Psi=0$. Hence, if the
thickness of the boundary layer is written as $\epsilon$, the size of the
streamfunction there can be estimated as
$O(\Psi|_{b})=O(\epsilon\Psi|_{c})=O(\epsilon Re_{w})$ (the subscript $b$
stands for boundary layer). In order to ensure the viscous-convective balance
of (11b) in this thin layer we further require
$O(\Psi|_{b})=O(\epsilon^{-1})$, and therefore $Re_{w}=O(\epsilon^{-2})$.

As seen in figure 4c, parts of the boundary layer are attached to the walls,
where the flow has to fulfil the no-slip conditions. In order for the
streamfunction to be modified in the near-wall boundary layer, both sides of
(11c) must balance, so $O(\Omega|_{b})=O(\epsilon^{-1}Re_{w})$ using
$\partial_{r}=O(\epsilon^{-1})$ and $O(\Psi|_{b})=O(\epsilon^{-1})$. Another
essential part of the boundary layer is the plume, where the layer is detached
from the wall. When the boundary layer leaves the wall, it typically forms a
sharp corner. Similar to the RBC cases, as we round the corner, the sizes of
$\Omega$ and $\Gamma$ are unchanged. This can be justified by considering a
streamline within the boundary layer. The streamline passes through the
$O(\epsilon^{1/2})$ neighbourhood of the corner, where the flow is inviscid,
and hence $\Omega$ and $\Gamma$ are functions of the streamfunction; see also
the RBC literature introduced at the beginning of this section.

The plume layer is no longer parallel to the wall, and the Coriolis force
acting there drives the whole Taylor cell. This plume balance allows us to
completely fix the flow scaling in terms of $Ta$. Assuming $\partial_{r}=O(1)$
and $\partial_{z}=O(\epsilon^{-1})$ in (11b), the viscous-convective terms of
$O(\epsilon^{-2}\Omega|_{b})=O(\epsilon^{-5})$ counter-balances the Coriolis
term of $O(\epsilon^{-1}Ta)$ when $\epsilon=Ta^{-1/4}$. Here we used the fact
that within the near-wall boundary layer, the size of $\Gamma$ is $O(1)$ from
the boundary conditions, and the same size must be used in the plume.

The asymptotic structure can be summarised as follows. Within the core region
we use the expansions

$\displaystyle\Psi=\epsilon^{-2}\omega_{0}\Psi_{0}(r,z)+\cdots,\qquad\Gamma=\gamma_{0}+\cdots,\qquad\Omega=\epsilon^{-2}\omega_{0}+\cdots,$
(18)

where $\omega_{0}$ and $\gamma_{0}$ are constants. The expansion is consistent
with the Prandtl-Batchelor structure to leading order, and $\Psi_{0}$ must be
determined by

$\displaystyle 1=-r^{-1}\\{(r^{-1}\Psi_{0r})_{r}+r^{-1}\Psi_{0zz}\\}$ (19)

in the aforementioned region $\mathcal{A}$ with the boundary condition
$\Psi_{0}=0$ on $\mathcal{C}$.

Let $l$ and $n$ be the distance along $\mathcal{C}$ and its inward normal,
respectively. Then using the rescaled normal variable $N=\epsilon^{-1}n$, the
flow within the boundary layer can be expanded as

$\displaystyle\Psi=\epsilon^{-1}\Psi^{(b)}(l,N)+\cdots,\qquad\Gamma=\Gamma^{(b)}(l,N)+\cdots.$
(20)

The leading order equations in the boundary layer are

$\displaystyle(\Psi_{l}^{(b)}\partial_{N}-\Psi_{N}^{(b)}\partial_{l})\Gamma^{(b)}$
$\displaystyle=$ $\displaystyle r\Gamma_{NN}^{(b)},$ (21)
$\displaystyle(\Psi_{l}^{(b)}\partial_{N}-\Psi_{N}^{(b)}\partial_{l})\Psi_{NN}^{(b)}$
$\displaystyle=$
$\displaystyle{r\Psi}_{NNNN}^{(b)}-\frac{2\cos\varphi}{r}\Gamma^{(b)}\Gamma_{N}^{(b)},$
(22)

where $\varphi$ is the angle between the $z$-axis and the $n$-axis.

When the boundary layer is in contact with the cylinder wall, the following
boundary conditions must be imposed at $N=0$:

$\displaystyle\Psi^{(b)}=\Psi_{N}^{(b)}=0,\qquad\Gamma^{(b)}=\Gamma_{i}\qquad\text{at
the inner cylinder,}$ (23)
$\displaystyle\Psi^{(b)}=\Psi_{N}^{(b)}=0,\qquad\Gamma^{(b)}=\Gamma_{o}\qquad\text{at
the outer cylinder.}$ (24)

For the near-wall boundary layer, $\cos\varphi=0$ so (22) is none other than
Prandtl’s boundary layer equation. While for large $N$, the flow must match
the core solution. Let $U(l)$ be the value of $\Psi_{0n}$ on $\mathcal{C}$,
which can be found by the core flow problem (19). Then as
$N\rightarrow\infty$, the boundary layer solution must satisfy

$\displaystyle\Psi^{(b)}\rightarrow\omega_{0}U(l)N,\qquad\Gamma^{(b)}\rightarrow\gamma_{0}.$
(25)

When the boundary layer is detached from the wall, the similar far-field
conditions should be applied on both sides $N\rightarrow\pm\infty$.

The theoretical boundary layer thickness $\epsilon=Ta^{-1/4}$ implies the
Nusselt number scaling $Nu\propto Ta^{\beta}$ with $\beta=1/4$. We can check
if this is consistent with the angular momentum transport. By averaging (11a)
with respect to $z$ and further radially integrating from $r_{i}$ to $r$, the
transport balance can be obtained as

$\displaystyle-\overline{\Psi_{z}\widetilde{\Gamma}}=r^{3}(r^{-2}\overline{\Gamma})_{r}+\frac{16\eta
r_{i}^{2}}{(1+\eta)^{4}}Nu.$ (26)

Here we used overline to denote the average with respect to $z$ and wrote
$\Gamma=\overline{\Gamma}(r)+\widetilde{\Gamma}(r,z)$. Let us consider this
balance at a sufficient distance from both walls. In the core region,
$\widetilde{\Gamma}$ would be exponentially small because of the Prandtl-
Batchelor theorem, so the transport is extremely inefficient. The plume region
of thickness $O(\epsilon)$ is hence the main contributor to the left-hand
side, which can be estimated as
$O(\epsilon(\partial_{z}\Psi|_{b})\widetilde{\Gamma}|_{b})=O(\epsilon^{-1})$
using the scaling $\widetilde{\Gamma}|_{b}=O(1)$,
$\Psi|_{b}=O(\epsilon^{-1})$. This term indeed balances with the second term
on the right-hand side.

Figure 6:  The large $Ta$ asymptotic convergence of the Taylor vortex for
$\eta=0.5,k=3,a=-1/8$. The red solid curve is almost horizontal, implying that
the $Nu\propto Ta^{1/4}$ scaling derived for $k=O(1)$ holds. The green dashed
and blue dotted curves are computed by $v$ and $\omega$ measured at the centre
of the cell $(r,z)=(r_{i}+0.5,\pi/2k)$, respectively.

The asymptotic structure derived here is consistent with the behaviour of the
numerical Taylor vortex solution. Figure 6 summarises the numerical results
for $\eta=0.5$, $a=-1/8$. The red solid curve in the figure shows the scaling
$Nu\propto Ta^{1/4}$, which can also be seen in the data presented in figure 1
($\eta=5/7,a=0$). The green dashed and the blue dotted curves are
$Ta^{-1/2}rv$ and $Ta^{-1/2}\omega/r$ measured at the centre of the roll cell.
According to the theory, they converge towards the constants $\gamma_{0}$ and
$\omega_{0}$, respectively. The scaling of $\omega/r$ corresponds to the
scaling of $u$, $w$ in the core, and hence the scaling of $Re_{w}$.

One may have noticed that the exponent $\beta=1/4$ of the Nusselt number
differs from the exponent $\beta=1/3$ deduced in the asymptotic analysis of
Chini & Cox (2009) and Hepworth (2014). This discrepancy is not due to the
differences in the flow driving mechanism but to the boundary conditions at
the walls. If the zero stress condition is imposed, the linear extrapolation
of the core streamfunction already satisfies the boundary condition to leading
order. This means that the balance $O(\Omega|_{b})=O(\epsilon^{-1}Re_{w})$ we
assumed for the no-slip case is not necessary. For the slip wall case, the
magnitude of the vorticity does not change in the core and the boundary layer,
so the strong vortex layer seen in figure 4c does not appear. If we use the
balance $O(\Omega|_{b})=O(\Omega|_{c})=O(\epsilon^{-2})$ instead for the
scaling argument of the plume we have $\epsilon=Ta^{-1/3}$, as expected.

For the slip wall RBC, further analytical progress has been made using the
fact that the boundary layer equations become linear and the roll cells are
rectangular (Chini & Cox 2009; Hepworth 2014). In our case, however, we have
to rely on numerical calculations because the boundary layer equations are
fully nonlinear. Moreover, the shape of the core region is nontrivial due to
the small vortices appearing near the corners (see figure 4c). This means that
the core and boundary layer equations (4.8), (4.10), (4.11) need to be
iteratively solved by updating the core shapes and constants $\gamma_{0}$ and
$\omega_{0}$. Such numerical calculations are too challenging and out of the
scope of this paper. Differences in the structure of the boundary layer also
affect the $Pr$ dependence of the asymptotic solution of RBC. For the slip
wall case, asymptotic solutions for different $Pr$ can be obtained by
rescaling the $Pr=1$ solution. However, this is possible because the boundary
layer is linear, which is not the case for the no-slip case.

(a) | (b)
---|---
|
Figure 7:  Uniform mean angular momentum profiles. (a) The mean flow of the
Taylor vortex solution shown in figure 4 ($\eta=0.5,a=-1/8$). The red solid
curve is the normalised angular momentum, while the blue dotted curve is the
mean angular velocity $q$ defined in (10). (b) Time average of DNS results for
$Ta=10^{10}$, $a=-0.20,0.00,0.21,0.40,0.60,1.00$. The data are from figure 4
of Ostilla-Mónico et al. (2014b).

Finally, we show that the above analysis provides some insights into the
structure of the mean flow. As already seen in figure 4b, the angular momentum
$rv$ is almost constant in the core region. Therefore, the mean angular
momentum $r\overline{v}$ is expected to become a constant away from the wall.
This is indeed the case for the numerical Taylor vortex solution (figure 7a).
In the studies of TCF turbulence, on the other hand, the mean angular velocity
$q$ is usually plotted, and it is often noted that the profile is linear in
the core. At first glance, this appears to be the case, for example, when
looking at figure 2a, but this is because the cylinder gap ($\eta\approx
0.714$) is too narrow to clearly see the radial dependence of the profile (see
figure 7a for the $q$ profile for a wide gap case $\eta=0.5$). If the
numerical results by Ostilla-Mónico et al. (2014b) are summarised in terms of
$rv$, as shown in figure 7b, they clearly show the constant angular momentum
property in the core. Note that when $a$ becomes too large, the Taylor vortex
appears to favour the vicinity of the inner cylinder, and the homogenisation
is only observed there. In the long history of TCF studies, some researchers
have also pointed out that the turbulent mean flows might have a uniform
angular momentum zone (Wattendorf 1935, Taylor 1935, Smith & Townsend 1982,
Lewis & Swinney 1999, Dong 2007, Brauckmann & Eckhardt 2017). However, this
fact is not widely known, probably because the mathematical reasons behind it
has not been elucidated.

##  5 High wavenumber Taylor vortices

The aim of this section is to examine the asymptotic behaviour of the
solutions at the first and second peaks seen in figure 3a. To this end, in
figure 8, we performed similar calculations at two higher Taylor numbers. The
results are summarised using the theoretical scaling to be derived in this
section. Essentially, the scaling of $k$ and $Nu$ represent the width of the
roll cell in the $z$ direction and the thickness of the boundary layer
adjacent to the cylinder wall, respectively. The computations of the solutions
are not easy, especially around the first peak, where very thin boundary
layers need to be resolved. Although the convergence of the numerical
solutions to the asymptotic states is still not perfect, the theoretical
scalings are consistent with the overall features of the numerical data.

(a) | (b)
---|---
|
Figure 8:  Change in the Nusselt number of the Taylor vortex solution when the wavenumber is varied. (a) and (b) use the same numerical results. The outer cylinder is stationary ($a=0$), and the radius ratio is $\eta=5/7$. Red solid, green dashed, and blue dotted curves correspond to $Ta=2.95\times 10^{11}$, $Ta=7.37\times 10^{10}$, and $Ta=9.75\times 10^{9}$, respectively. The blue dotted curve is the same as that shown in figure 3. The three crosses in figure 1 are taken from the maxima seen in (a).  (a) | (b)
---|---
|
Figure 9:  Mean flows for the solutions at the extrema of $Nu$ seen in figure
8 ($\eta=5/7,a=0$). The red curves are the numerical results at $Ta=2.95\times
10^{11}$. (a) The first peak. The black solid line is the asymptotic result
$\overline{\Gamma}=\gamma_{0}$. The value of $\gamma_{0}=0.771$ is estimated
at the mid gap. (b): The second peak. The black curve is the asymptotic result
(39). The value of $A=335.5$ is estimated at the mid gap (see (75)).

The structure of the flow field in both peaks can be roughly divided into a
boundary layer near the wall and a core region in the middle of the gap, as we
have seen in figures 2 and 3. On closer inspection, one further notices that
the asymptotic structures of the two flows are quite different. For example,
in the first peak solution, $\overline{\Gamma}$ has a flat profile in the core
region (figure 9a), but this is not the case in the second peak solution
(figure 9b). Figure 10 examines how the core flows develop from the near wall
region adjacent to the inner cylinder. In the first peak solution, the near-
wall structure is somewhat similar to the $k=O(1)$ case, with the wall
boundary layer becoming a sharp plume as rounding the corner (top panel). On
the other hand, in the second peak solution, no apparent plume can be
recognised away from the wall, and the flow only varies slowly in the
$z$-direction (bottom panel). The core flow inherits this property, as seen
from figure 11, where the axial structure of the fluctuation fields is plotted
at the midgap $r=r_{m}=(r_{o}+r_{i})/2$. For the second peak solution, only a
single Fourier mode plays a major role in the core, while for the first peak,
a number of harmonics participate in forming the plume. A theoretical
explanation of those differences will be deduced below, together with the
detailed scalings of the flows.

Figure 10: The colour map of $\omega/r$ for the first peak (top panel) and the
second peak (bottom panel) solutions seen in figure 8. $Ta=2.95\times
10^{11}$. The centre of the colour bar is zero.

(a) first peak | (b) second peak
---|---
|

(c) first peak | (d) second peak
---|---
|

Figure 11:  The flow fields at the mid gap $r=r_{m}$ for the first peak and
second peak solutions in figure 8. The red solid, green dashed, and blue
dotted curves correspond to $Ta=2.95\times 10^{11}$, $Ta=7.37\times 10^{10}$
and $Ta=9.75\times 10^{9}$, respectively. Note that the value of $k$ is
determined by the optimisation of $Nu$ and therefore varies from solution to
solution. The difference in the structure of the first and second peak
solutions can be more clearly seen in figure 10, where the axial coordinates
are not scaled.

###  5.1 The first peak asymptotic structure

Now let us find the scaling of the first peak solution. Hereafter $\delta=1/k$
denotes the length scale of the cell in the $z$ direction. The flow can be
divided into core and near-wall zones, as shown in figure 5c. Within the
regions $O(\delta)$ away from the cylinder walls, similar asymptotic arguments
as the $k=O(1)$ case would apply. Figure 5b shows an enlarged view of the near
wall zone where a thin boundary layer, of thickness $\epsilon\ll\delta$ say,
emerges around the inviscid region.

The size of the streamfunction in the inviscid region must be the same as that
of the core, $\Psi|_{c}$. Thus the size of the streamfunction within the
boundary layer can be found as $O(\Psi|_{b})=O(\epsilon\delta^{-1}\Psi|_{c})$,
noting that the boundary layer thickness relative to the near wall region is
$\epsilon/\delta$. Further using the viscous-convective balance in the
boundary layer $O(\Psi|_{b})=O(\epsilon^{-1}\delta)$, we can find
$O(\Psi|_{c})=O(\delta^{2}\epsilon^{-2})$. When $\overline{\Gamma}|_{c}$ is
almost constant, the viscous-convective balance in the core is satisfied if
$O(\Psi|_{c})=O(\delta^{-1})$. Hence, from
$O(\Psi|_{c})=O(\delta^{2}\epsilon^{-2})$ obtained in the near wall analysis,
the two small perturbation parameters are related as $\epsilon=\delta^{3/2}$.

Within the near-wall zone, the plume coming out of the wall boundary layer is
so thin that the Coriolis forces and viscous terms cannot balance, unlike the
$k=O(1)$ case. This means that the Coriolis force effect should appear as the
plume diffuses towards the core region.. The viscous-Coriolis balance in the
core can be described as
$O(\Psi|_{c}\delta^{-4})=O(Ta\,\delta^{-1}\widetilde{\Gamma}|_{c})$. To
estimate $O(\widetilde{\Gamma}|_{c})$, we can use the balance
$O(\widetilde{\Gamma}|_{c}\Psi|_{c})=O(\Gamma|_{b}\Psi|_{b})$ which will be
justified shortly. Using the argument of the inviscid corner
$O(\Gamma|_{b})=O(1)$. Therefore we have the estimation
$O(\widetilde{\Gamma}|_{c})=O(\epsilon\delta^{-1})$, which is consistent with
the momentum transport balance
$O((\partial_{z}\Psi|_{c})\widetilde{\Gamma}|_{c})=O(Nu)=O(\epsilon^{-1})$
deduced from (26). This is the final key to unlock the scaling
$\delta=O(Ta^{-2/9})$ and $\epsilon=O(Ta^{-1/3})$ with the latter motivated us
to use the $Nu\propto Ta^{1/3}$ scaling in figure 8a. If the wind Reynolds
number $Re_{w}$ is defined as the magnitude of the wall-normal velocity
component, $Re_{w}=O(\delta^{-1}\Psi|_{c})=O(Ta^{4/9})$.

The asymptotic analysis in the core is summarised as follows. First, we choose
$\delta$ as a small perturbation parameter and rescale the axial coordinate as
$Z=\delta^{-1}z$. From the scaling argument, the Taylor number has the
expansion $Ta=\delta^{-9/2}T_{0}+\cdots$. Using the core expansions

$\displaystyle\Psi=\delta^{-1}\Psi^{(c)}(r,Z)+\cdots,~{}~{}\Gamma=\gamma_{0}+\delta^{1/2}\Gamma^{(c)}(r,Z)+\cdots$
(27)

with a constant $\gamma_{0}$ to (11) yield the leading order equations

$\displaystyle(\Psi_{r}^{(c)}\partial_{Z}-\Psi_{Z}^{(c)}\partial_{r})\widetilde{\Gamma}^{(c)}=r_{i}\widetilde{\Gamma}_{ZZ}^{(c)},$
(28)
$\displaystyle(\Psi_{r}^{(c)}\partial_{Z}-\Psi_{Z}^{(c)}\partial_{r})\Psi_{ZZ}^{(c)}=r_{i}\Psi_{ZZZZ}^{(c)}-T_{0}\frac{2\gamma_{0}}{r_{i}}\widetilde{\Gamma}_{Z}^{(c)}.$
(29)

Here only the fluctuation part of the equations was extracted. The mean part
can be directly obtained from (26) as

$\displaystyle-\overline{\Psi_{Z}^{(c)}\widetilde{\Gamma}^{(c)}}=N_{0},$ (30)

using the scaled Nusselt number $N_{0}=\frac{16\eta
r_{i}^{2}}{(1+\eta)^{4}}\epsilon Nu$.

In the plume near the inner wall, we use the expansion

$\displaystyle\Psi=\epsilon^{-1}\delta\Psi^{(b)}(l,N)+\cdots,~{}~{}\Gamma=\Gamma^{(b)}(l,N)+\cdots$
(31)

together with the stretched variables $l=\frac{r-r_{i}}{\delta}$ and
$N=z/\epsilon$. The leading order part of (4.1a) takes the form

$\displaystyle\Psi_{l}^{(b)}\Gamma_{N}^{(b)}-\Psi_{N}^{(b)}\Gamma_{l}^{(b)}=r_{i}\Gamma_{NN}^{(b)}$
(32)

which is essentially identical to (4.10). If we consider $\Gamma^{(b)}$ is a
function of $\psi=\Psi^{(b)}$ and $l$, the equation becomes

$\displaystyle-\Gamma_{l}^{(b)}=r_{i}(\psi_{z}\Gamma_{\psi}^{(b)})_{\psi}$
(33)

where $\Gamma_{\psi}^{(b)}$ should be small when far enough away from the
plume. By integrating this equation across the plume we get

$\displaystyle\frac{d}{dl}\left(\int^{\infty}_{-\infty}\Gamma^{(b)}d\psi\right)=0.$
(34)

The term in the bracket is conserved during the plume diffusion so that the
balance $O(\widetilde{\Gamma}|_{c}\Psi|_{c})=O(\Gamma|_{b}\Psi|_{b})$ must be
satisfied. The effect of the plume entering the core has to be described by
boundary conditions for (5.2)-(5.3) at $r=r_{i},r_{o}$. The conditions could
be written using the Dirac delta function as done in Vynniycky & Masuda
(2013).

We do not solve the asymptotic problem, but the theory well explains the
beheviour of the Taylor vortex solution. For example, the numerical results
shown in figures 11a, c are summarised using the theoretical core scalings
$\Omega|_{c}=O(\delta^{-2}\Psi|_{c})=O(Ta^{2/3})$ and
$\widetilde{\Gamma}|_{c}=O(Ta^{-1/9})$. In figure 8a, the thin black line is
the asymptotic approximation $\overline{\Gamma}\approx\gamma_{0}$ with the
estimate $\gamma_{0}=Ta^{-1/2}r_{m}\overline{v}(r_{m})\approx 0.771$.

From the numerical data, the Nusselt number at the first peak seems to be
roughly approximated by $Nu=0.027\,Ta^{1/3}$. The asymptotic line intersect
with the $k=3$ curve shown in figure 1 at $Ta\approx 6\times 10^{7}$, which is
not a bad approximation of the transition point between the classical and
ultimate turbulence regimes.

###  5.2 The second peak asymptotic structure

The asymptotic structure of the second peak solution is more complex because
three layers appear in the vicinity of the walls (see figure 12a). We refer to
them as the bottom, middle and top boundary layers, in order from the cylinder
wall. Understanding their precise scaling requires a delicate matched
asymptotic expansion analysis, which we will leave for the next section. Here
we shall look only at how the flow scaling is intuitively determined. The
argument below contains one assumption that is not strictly fulfilled, but the
resultant scaling is nonetheless correct if some minor logarithmic factor
effects are omitted.

Figure 12:  Similar sketch as figure 5 but for the asymptotic states when
$k=Ta^{1/4}$. (a) The second peak state. BL stands for boundary layer. (b) The
transitional state. The shear layer occurs around the critical radius
$r=r_{c}$.

The key assumption is $\frac{d\overline{\Gamma}}{dr}=O(1)$ in the core, which
did not hold for the first peak case. If the governing equations (11a) and
(11b) are linearised around the mean flow, the balances
$O(\delta^{-1}\Psi|_{c})=O(\delta^{-2}\widetilde{\Gamma})$ and
$O(\delta^{-4}\Psi|_{c})=O(\delta^{-1}Ta\,\widetilde{\Gamma}|_{c})$ must be
satisfied. Those balances determine $\delta=O(Ta^{-1/4})$ and
$O(\widetilde{\Gamma}|_{c})=O(\delta\Psi|_{c})$. Further using
$O((\partial_{z}\Psi|_{c})\widetilde{\Gamma}|_{c})=O(Nu)=O(\epsilon^{-1})$
deduced from the mean equation (26), the core flow scaling can be written as
$O(\Psi|_{c})=O(\epsilon^{-1/2})$,
$O(\widetilde{\Gamma}|_{c})=O(\delta\epsilon^{-1/2})$.

Within the bottom boundary layer, the viscous-convective balance
$O(\Psi|_{b})=O(\epsilon^{-1}\delta)$ must, of course, be satisfied (the
subscript $b$ implies that $\Psi$ is measured at the bottom boundary layer).
If we further assume the viscous-Coriolis balance
$O(\Psi|_{b})=O(Ta\,\epsilon^{4}\delta^{-1})$, we obtain
$\epsilon=\delta^{6/5}=Ta^{-3/10}$, which gives $Nu\propto Ta^{3/10}$ in
figure 7b. It is this balance that is mostly met, but in fact, is a little
broken. As mentioned earlier, the effect of this slight imbalance appears as a
logarithmic factor in the scaling of $\epsilon$ in the matched asymptotic
expansion. The introduction of the logarithmic factor only has little
influence when examining the scaling of numerical data, and so it was ignored
in figure 7b. Furthermore, as shown in figures 11b, d, the core scalings
$O(\Omega|_{c})=O(\delta^{-2}\Psi|_{c})=O(Ta^{13/20})$ and
$\widetilde{\Gamma}|_{c}=O(Ta^{-1/10})$ without the logarithmic factor well
explain the numerical results. Note that the scaling result suggests
$Re_{w}=O(\delta^{-1}\Psi|_{c})=O(\delta^{-1}\epsilon^{-1/2})=O(Ta^{2/5})$.

In the middle boundary layer, the flow coming towards the wall turns back in
the region of aspect ratio about unity, so its thickness should be
$O(\delta)$.

The top boundary layer is necessary for the nonlinearity of the fluctuation
component to disappear towards the core. Its thickness
$\Delta=O(\delta\epsilon^{-1/2})=O(Ta^{-1/10})$ can be obtained by the
viscous-convective balance $O(\Delta^{-1}\delta^{-1}\Psi|_{t})=O(\delta^{-2})$
and $O(\Psi|_{t})=O(\Psi|_{c})$, where $O(\Psi|_{t})$ is the size of the
streamfunction in the top boundary layer. One of the ways to ascertain the
presence of the top boundary layer in the visualised flow field is to look at
the extreme values of $\Omega$ (see figure 9, bottom panel) or
$\widetilde{\Gamma}$. Figure 13a shows the scaled fluctuation
$\widetilde{\Gamma}$ at the plume position $z=\pi/k$. It can be seen that the
minima approach the wall as the Taylor number increases. In figure 13b, the
distance between the wall and the minima is scaled by the theoretical top
boundary layer thickness $\Delta$. As expected, all the curves almost overlap
around the minima.

##  6 Matched asymptotic expansion when $k=O(Ta^{1/4})$

(a) | (b)
---|---
|
Figure 13:  The profile of $\widetilde{\Gamma}$ at the plume position
$z=\pi/k$ for the Taylor vortex solution at the second peak. $\eta=5/7,a=0$.
The red solid, green dashed, and blue dotted curves correspond to
$Ta=2.95\times 10^{11}$, $Ta=7.37\times 10^{10}$ and $Ta=9.75\times 10^{9}$,
respectively.

The goal of this section is to derive a matched asymptotic expansion
consistent with the behaviour of the second peak solutions. This asymptotic
state can be reached by decreasing the wavenumber from the linear critical
point, as seen in figure 8b. However, for $k/Ta^{1/4}$ ranging from about 0.6
to 0.8, the appropriate scaling of the Nusselt number is $O(Ta^{0})$, which is
different from that observed for the second peak state; see figure 14a. We
shall study this transitional state in section 6.1, followed by the analysis
of the second peak in section 6.2.

The most important previous work throughout this section is due to Denier
(1992), who applied the Hall & Lakin (1988) type high wavenumber asymptotic
theory to Taylor-Couette flow in the narrow gap limit. As briefly commented in
that paper, the extension of the theory to wide-gap cases should be
straightforward. In section 6.1, we derive an analytic expression of the
Nusselt number and compare it with the numerical solution for the first time.

As the solution moves away from the linear critical point in the transitional
state, the vortex grows from the vicinity of the inner cylinder (figure 14c).
The asymptotic state corresponding to the second peak appears when the vortex
fills the entire gap. A similar scenario was suggested in Denier (1992) using
a matched asymptotic expansion, but the scaling obtained is unfortunately not
consistent with the behaviour of the numerical solutions. The main reason for
this is that the matching is not actually possible in the two-layer boundary
layer structure assumed in Denier (1992). In section 6.2, we will resolve that
difficulty by modifying the structure of the near wall-zone to three layers.

We choose $\delta=k^{-1}$ as a small perturbation parameter and introduce the
scaled axial variable $Z=\delta^{-1}z$. According to Hall (1982), the linear
critical point of curved flow problems typically has the Taylor number
expansion $Ta=\delta^{-4}T_{0}+\cdots$. This is also true for Taylor Couette
flow from the linear critical point to the second peak.

(a) | (b)
---|---
|

(c)

Figure 14: (a): The close-up of figure 8b around the bifurcation point. The
black solid curve is the asymptotic result (53), which has the property that
$Nu=1$ at $k/Ta^{1/4}=K_{1}$, and that $Nu\rightarrow\infty$ as
$k/Ta^{1/4}\rightarrow K_{\infty}$. (b) The mean flow profile at
$k/Ta^{1/4}=0.65,0.7$. The black solid curves are the asymptotic result (51).
The bullets indicate $r=r_{c}$ given in (47). The red points are the numerical
Taylor vortex result for $Ta=9.75\times 10^{9}$. (c) The colourmap of
$\omega/r$ for the numerical Taylor vortex. The top and bottom panels
correspond to $k/Ta^{1/4}=0.7$ and 0.65, respectively. Tics are placed at
$r=r_{c}$.

###  6.1 The transitional states around the linear critical point

For simplicity, we fix the outer cylinder (i.e. $a=\Gamma_{o}=0$) as Denier
(1992). The analysis of the general cases ($a\neq 0$) is relegated to Appendix
B as it is somewhat complex. The asymptotic structure of the transitional
state is depicted in figure 12b. The vortices are concentrated in the core
region $r_{i}<r<r_{c}$, where $r_{c}\in[r_{i},r_{o}]$ is the critical radius
to be determined analytically. We can show that the vortex amplitude decays
exponentially within a shear layer of thickness $O(\Delta)=O(Ta^{-1/6})$
appearing around the critical radius. The structure of the shear layer is
identical to that described in Denier (1992) and is not discussed in detail
here. The only important fact that will be used later is that the mean flow
and its derivative are continuous across this layer.

In the region where $r$ is greater than $r_{c}$ it is sufficient to analyse
the mean flow. To leading order, the left-hand side of (26) can be ignored and
hence the mean flow is a linear combination of a constant and $r^{2}$. Writing

$\displaystyle N_{0}=\frac{16\eta r_{i}^{2}}{(1+\eta)^{4}}Nu$ (35)

for simplicity, we get the asymptotic approximation

$\displaystyle\overline{\Gamma}=B(r_{o}^{2}-r^{2})+\cdots,\qquad
B=\frac{N_{0}}{2r_{o}^{2}}.$ (36)

Note that we used the fact that $\overline{\Gamma}$ should vanish at
$r=r_{o}$.

In the core region, the appropriate expansions are

$\displaystyle\Psi=\widehat{\Psi}^{(c)}(r)\sin
Z+\cdots,~{}~{}\widetilde{\Gamma}=\delta\widehat{\Gamma}^{(c)}(r)\cos
Z+\cdots,~{}~{}\overline{\Gamma}=\overline{\Gamma}^{(c)}(r)+\cdots.$ (37)

Substituting them into (11), the leading order parts yield

$\displaystyle\overline{\Gamma}_{r}^{(c)}\widehat{\Psi}^{(c)}=r\widehat{\Gamma}^{(c)},\qquad
0=\widehat{\Psi}^{(c)}+T_{0}\frac{2\overline{\Gamma}^{(c)}}{r^{2}}\widehat{\Gamma}^{(c)}.$
(38)

For those equations to have a nontrivial solution
$0=r^{3}+2T_{0}\overline{\Gamma}^{(c)}\overline{\Gamma}^{(c)}_{r}$ must be
satisfied. Therefore the mean flow is obtained as

$\displaystyle\overline{\Gamma}^{(c)}(r)=\sqrt{\frac{A-r^{4}}{4T_{0}}}.$ (39)

The constant $A$ must be

$\displaystyle A=r_{i}^{4}+4T_{0}\Gamma_{i}^{2}$ (40)

to satisfy $\overline{\Gamma}^{(c)}(r_{i})=\Gamma_{i}$.

A boundary layer of thickness $O(\delta)$ is needed near the inner cylinder to
satisfy the no-slip boundary conditions. In this layer we introduce the
stretched variable $\xi=\frac{r-r_{i}}{\delta}$ and assume the expansion

$\displaystyle\Psi=\Psi^{(b)}(\xi,Z)+\cdots,~{}~{}\Gamma=\Gamma_{i}+\delta\Gamma^{(b)}(\xi,Z)+\cdots.$
(41)

The leading order equations are obtained as

$\displaystyle(\Psi_{\xi}^{(b)}\partial_{Z}-\Psi_{Z}^{(b)}\partial_{\xi})\Gamma^{(b)}=r_{i}(\partial_{\xi}^{2}+\partial_{Z}^{2})\Gamma^{(b)},$
(42)
$\displaystyle(\Psi_{\xi}^{(b)}\partial_{Z}-\Psi_{Z}^{(b)}\partial_{\xi})(\partial_{\xi}^{2}+\partial_{Z}^{2})\Psi^{(b)}=r_{i}(\partial_{\xi}^{2}+\partial_{Z}^{2})^{2}\Psi^{(b)}-T_{0}\frac{2\Gamma_{i}\Gamma_{Z}^{(b)}}{r_{i}}.$
(43)

The no-slip boundary conditions are

$\displaystyle\Psi^{(b)}=\Psi_{\xi}^{(b)}=0,\qquad\Gamma^{(b)}=0\qquad\text{at}\qquad\xi=0,$
(44)

while for the far-field $\xi\rightarrow\infty$ we use the matching conditions

$\displaystyle\Psi^{(b)}\rightarrow\widehat{\Psi}^{(c)}(r_{i})\sin
Z,\qquad\widetilde{\Gamma}^{(b)}\rightarrow\widehat{\Gamma}^{(c)}(r_{i})\cos
Z,\qquad\overline{\Gamma}^{(b)}\rightarrow\overline{\Gamma}^{(c)}_{r}(r_{i})\xi.$
(45)

Finally, we determine the unknown constants $B$ and $r_{c}$ from the
continuity of $\overline{\Gamma}$ and $r^{3}(r^{-2}\overline{\Gamma})_{r}$ at
$r=r_{c}$. Using the mean flows (39) and (36), the continuity conditions can
be written as

$\displaystyle\sqrt{\frac{A-r_{c}^{4}}{4T_{0}}}=B(r_{o}^{2}-r_{c}^{2}),\qquad\frac{A}{\sqrt{T_{0}(A-r_{c}^{4})}}=2Br_{o}^{2},$
(46)

where $A$ is known by (40). From those equations it is easy to find

$\displaystyle r_{c}=\sqrt{\frac{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}{r_{o}^{2}}},$
(47)
$\displaystyle\frac{1}{B}=\sqrt{4T_{0}\left(\frac{r_{o}^{4}}{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}-1\right)}.$
(48)

Thus if $T_{0}$ is given the leading order mean flow is completely determined
as

$\displaystyle\overline{\Gamma}\approx\left\\{\begin{array}[]{c}(r_{o}^{2}-r^{2})\sqrt{\frac{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}{4T_{0}\left(r_{o}^{4}-r_{i}^{4}-4T_{0}\Gamma_{i}^{2}\right)}}\qquad\text{if}\qquad
r>r_{c},\\\
\sqrt{\frac{r_{i}^{4}-r^{4}+4T_{0}\Gamma_{i}^{2}}{4T_{0}}}\qquad\text{if}\qquad
r\leq r_{c}.\end{array}\right.$ (51)

This is the black curve in figure 14b, which predicts the numerical results
very well. Note that the analytic expression of $\widehat{\Psi}^{(c)}$ and
$\widehat{\Gamma}^{(c)}$ can also be found by (38) and the leading order part
of the mean flow equation (26)

$\displaystyle-\frac{\overline{\Gamma}^{(c)}_{r}}{2r}(\widehat{\Psi}^{(c)})^{2}=r^{3}(r^{-2}\overline{\Gamma}^{(c)})_{r}+N_{0}.$
(52)

The Nusselt number is found by (35), (36), (48) as

$\displaystyle
Nu(T_{0})=\frac{(1+\eta)^{4}}{16\eta^{3}\sqrt{T_{0}}}\sqrt{\frac{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}{r_{o}^{4}-r_{i}^{4}-4T_{0}\Gamma_{i}^{2}}}.$
(53)

This function is plotted by the black curve in figure 14a using
$k/Ta^{1/4}=T_{0}^{-1/4}$ for the horizontal axis. The other curves in the
figure, corresponding to the Taylor vortex solutions at finite $Ta$, clearly
approach the asymptotic result as $Ta\rightarrow\infty$.

As the scaled wavenumber decreases from the linear critical point
$K_{1}=T_{1}^{-1/4}$, the Nusselt number increases from the laminar value
unity. The asymptotic result
$T_{1}=\frac{r_{i}^{2}(r_{o}^{2}-r_{i}^{2})}{4\Gamma_{i}^{2}}$ can be found
from (47) by setting $r_{c}=r_{i}$ and $T_{0}=T_{1}$.

Similar to the narrow gap case studied by Denier et al. (1992), $Nu$ blows up
at finite $T_{0}$. Writing $K_{\infty}=T_{\infty}^{-1/4}$, we can deduce
$T_{\infty}=\frac{r_{o}^{4}-r_{i}^{4}}{4\Gamma_{i}^{2}}$ from (47) by using
$r_{c}=r_{o}$ and $T_{0}=T_{\infty}$. That is, in figure 14a, $K_{\infty}$ is
the point at which the flow switches from the transitional state to the second
peak state.

###  6.2 The second peak states

We now turn our attention to the second peak asymptotic states. Again $\delta$
is the small asymptotic parameter, and $Ta=\delta^{-4}T_{0}+\cdots$. For the
bottom boundary layer thickness $\epsilon$, we impose the condition

$\displaystyle\frac{\delta^{6}}{\epsilon^{5}}=\ln\frac{\delta}{\epsilon},$
(54)

which is needed for successful matching. The boundary layer thickness is
estimated as $O(\epsilon)=O(Ta^{-3/10}(\ln Ta)^{-1/5})$. Therefore the Nusselt
number and the wind Reynolds number scalings now involve the logarithmic
correction as $Nu=O(\epsilon^{-1})=O(Ta^{3/10}(\ln Ta)^{1/5})$,
$Re_{w}=O(\delta^{-1}\Psi|_{c})=O(\delta^{-1}\epsilon^{-1/2})=O(Ta^{2/5}(\ln
Ta)^{1/10})$.

Let us start with the core analysis by writing

$\displaystyle\Psi=\epsilon^{-1/2}\widehat{\Psi}^{(c)}(r)\sin Z+\cdots,~{}~{}$
(55)
$\displaystyle\widetilde{\Gamma}=\delta\epsilon^{-1/2}\widehat{\Gamma}^{(c)}(r)\cos
Z+\cdots,~{}~{}\overline{\Gamma}=\overline{\Gamma}^{(c)}(r)+\cdots.$ (56)

Those expansions, based on the argument in section 5, yield the identical
leading order equations (38) as the transitional state. As a result, the mean
flow is given by (39). However, for the second peak state, the multi-layered
near-wall structure must also be present near the outer cylinder, and thus
there is no way to determine the constant $A$ a priori. The leading order part
of the mean flow equation (26) is obtained as

$\displaystyle-\frac{\overline{\Gamma}^{(c)}_{r}}{2r}(\widehat{\Psi}^{(c)})^{2}=N_{0},$
(57)

where $N_{0}=\frac{16\eta r_{i}^{2}}{(1+\eta)^{4}}\epsilon Nu$. The unknowns
appeared in the core solutions $\overline{\Gamma}^{(c)},\widehat{\Psi}^{(c)}$,
and $\widehat{\Gamma}^{(c)}$ are $A$ and $N_{0}$.

The two near-wall regions have identical asymptotic structure, so only the
layers near the inner cylinder are examined below. The top boundary layer
expansions are

$\displaystyle\Psi=\epsilon^{-1/2}\Psi^{(t)}(\zeta,Z)+\cdots,~{}~{}\Gamma=\gamma_{i}+\delta\epsilon^{-1/2}\Gamma^{(t)}(\zeta,Z)+\cdots,$
(58)

where the constant $\gamma_{i}=\overline{\Gamma}^{(c)}(r_{i})$ comes from the
leading order core mean flow. The stretched variable
$\zeta=\frac{r-r_{i}}{\Delta}$ is defined by using the top boundary layer
thickness $\Delta=\delta\epsilon^{-1/2}=O(Ta^{-1/10}(\ln Ta)^{1/10})$. The
leading order governing equations for the top boundary layer are similar to
those for the core region in the first peak states:

$\displaystyle(\Psi_{\zeta}^{(t)}\partial_{Z}-\Psi_{Z}^{(t)}\partial_{\zeta})\widetilde{\Gamma}_{t}=r_{i}\widetilde{\Gamma}_{ZZ}^{(t)},$
(59)
$\displaystyle(\Psi_{\zeta}^{(t)}\partial_{Z}-\Psi_{Z}^{(t)}\partial_{\zeta})\Psi_{ZZ}^{(t)}=r_{i}\Psi_{ZZZZ}^{(t)}-T_{0}\frac{2\gamma_{i}}{r_{i}}\widetilde{\Gamma}_{Z}^{(t)},$
(60) $\displaystyle-\overline{\Psi_{Z}^{(t)}\widetilde{\Gamma}^{(t)}}=N_{0}.$
(61)

As $\zeta\rightarrow\infty$, the flow matches to the core solution as

$\displaystyle\Psi^{(t)}\rightarrow\widehat{\Psi}^{(c)}(r_{i})\sin
Z,\qquad\widetilde{\Gamma}^{(t)}\rightarrow\widehat{\Gamma}^{(c)}(r_{i})\cos
Z,\qquad\overline{\Gamma}^{(t)}\rightarrow\overline{\Gamma}^{(c)}_{r}(r_{i})\zeta.$
(62)

While the asymptotic behaviour of the flow towards the wall should be

$\displaystyle\Psi^{(t)}\rightarrow\zeta^{1/3}\widehat{\Psi}^{(t)}(Z),\qquad\Gamma^{(t)}\rightarrow\zeta^{-1/3}\widehat{\Gamma}^{(t)}(Z)\qquad\text{as}\qquad\zeta\rightarrow
0,$ (63)

which is similar to that used in Denier (1992). However, to match this with
the bottom boundary layer, we must insert the middle boundary layer.

Recall that the middle boundary layer is the special place where the radial
and axial derivatives of the flow have the same size. Thus the expansions
there must be written in terms of $\xi=\frac{r-r_{i}}{\delta}$ as

$\displaystyle\Psi=\epsilon^{-1/3}\Psi^{(m)}(\xi,Z)+\cdots,~{}~{}\Gamma=\gamma_{i}+\delta\epsilon^{-2/3}\Gamma^{(m)}(\xi,Z)+\cdots.$
(64)

Given the expansions it is a straightforward task to find that the leading
order equations

$\displaystyle(\Psi_{\xi}^{(m)}\partial_{Z}-\Psi_{Z}^{(m)}\partial_{\xi})\widetilde{\Gamma}^{(m)}=0,$
(65)
$\displaystyle(\Psi_{\xi}^{(m)}\partial_{Z}-\Psi_{Z}^{(m)}\partial_{\xi})(\partial_{\xi}^{2}+\partial_{Z}^{2})\Psi^{(m)}=-T_{0}\frac{2\gamma_{i}}{r_{i}}\widetilde{\Gamma}_{Z}^{(m)},$
(66) $\displaystyle-\overline{\Psi_{Z}^{(m)}\widetilde{\Gamma}^{(m)}}=N_{0},$
(67)

are inviscid. Here the Coriolis force term participating in the second
equation is critical for successful matching to the top boundary layer:

$\displaystyle\Psi^{(m)}\rightarrow\xi^{1/3}\widehat{\Psi}^{(t)}(Z),\qquad\Gamma^{(m)}\rightarrow\xi^{-1/3}\widehat{\Gamma}^{(t)}(Z)\qquad\text{as}\qquad\xi\rightarrow\infty.$
(68)

On the other hand, the appropriate behaviour of the solution towards the wall
$\xi\rightarrow 0$ can be found as

$\displaystyle\Psi^{(m)}\rightarrow\xi(-\ln\xi)^{1/3}\widehat{\Psi}^{(b)}(Z),\qquad\Gamma^{(m)}\rightarrow\xi^{-1}(-\ln\xi)^{-1/3}\widehat{\Gamma}^{(b)}(Z)$
(69)

by seeking the consistent limiting behaviour of the solution. The key
observation here is that for both matching conditions, the rate of increase in
$\Psi$ and the rate of decrease in $\Gamma$ are the same when moving away from
the wall. This is a requirement to satisfy the angular momentum transport
(67).

The leading order part of the bottom boundary layer expansions

$\displaystyle\Psi=\frac{\delta}{\epsilon}\Psi^{(b)}(Y,Z)+\cdots,~{}~{}\Gamma=\Gamma^{(b)}(Y,Z)+\cdots$
(70)

matches to the middle boundary layer if (54) holds and

$\displaystyle\Psi^{(b)}\rightarrow
Y\widehat{\Psi}^{(b)}(Z),\qquad\Gamma^{(b)}\rightarrow
Y^{-1}\widehat{\Gamma}^{(b)}(Z)$ (71)

as $Y=\frac{r-r_{i}}{\epsilon}\rightarrow\infty$. The leading order equations
in the bottom boundary layer are

$\displaystyle(\Psi_{Y}^{(b)}\partial_{Z}-\Psi_{Z}^{(b)}\partial_{Y})\Gamma^{(b)}=r_{i}\Gamma_{YY}^{(b)},$
(72)
$\displaystyle(\Psi_{Y}^{(b)}\partial_{Z}-\Psi_{Z}^{(b)}\partial_{Y})\Psi_{YY}^{(b)}=r_{i}\Psi_{YYYY}^{(b)}.$
(73)

Thanks to the radial diffusivity the no-slip boundary conditions

$\displaystyle\Psi^{(b)}=\Psi_{Y}^{(b)}=0,\qquad\Gamma^{(b)}=\Gamma_{i}\qquad\text{at}\qquad
Y=0$ (74)

can be imposed.

Although the combined boundary layers appear rather complex, their role is
merely to define the relationship between $A$ and $N_{0}$. Once $T_{0}$ and
$A$ are fixed, the flow near the inner cylinder can be calculated to find
$N_{0}=f_{i}(A,T_{0})$. A similar calculation for the near outer cylinder
region yields another functional relationship $N_{0}=f_{o}(A,T_{0})$. Hence,
in principle, $A(T_{0})$ and $N_{0}(T_{0})$ can be determined by solving the
two near-wall regions numerically. We do not perform such challenging
calculations since we have already seen in section 5 that the finite $Ta$
numerical results are consistent with the asymptotic theory.

The value of $A$ can be predicted from the finite $Ta$ numerical solutions as

$\displaystyle A=(2k^{-2}r_{m}\overline{v}(r_{m}))^{2}+r_{m}^{4}$ (75)

by using the mean flow at the mid gap. The black curve in figure 9b is the
result of using this $A$ in (39), and it can be seen that this asymptotic
prediction matches the numerical calculation in the entire core region.
Similarly, the fluctuation component in the core can be calculated explicitly.
The black curve in figure 15 is the asymptotic approximation

$\displaystyle\widetilde{\Gamma}\approx Ta^{-1/4}r\left(\frac{32\eta
r_{i}^{2}}{(1+\eta)^{4}}\frac{Nu}{\sqrt{(A-r^{4})}}\right)^{1/2}\cos Z,$ (76)

which is in good agreement with the numerical result.

Figure 15: The red dotted curve is the profile of $\widetilde{\Gamma}$ at the plume position $z=\pi/k$ for the Taylor vortex solution at the second peak. The same data as figure 13 ($\eta=5/7,a=0,Ta=2.95\times 10^{11}$). The black curve is the asymptotic result (76). (a) | (b)
---|---
|
Figure 16:  (a) The Nusselt number at $Ta=10^{10}$ for $\eta=5/7$. The red
crosses are the Taylor vortex solution with the optimised wavenumbers shown in
(b). The cyan squares are experimental results by Dennis et al. (2011). (b)
The red crosses are the local maximum of $Nu$ closest to the bifurcation
point. According to the asymptotic theory, the Taylor vortex solution
bifurcates from circular Couette flow at $K_{1}$ (see (87)). The Nusselt
number $Nu$ is $O(1)$ when $k/Ta^{1/4}\in(K_{\infty},K_{1})$, where
$K_{\infty}$ is given in (91). For $k/Ta^{1/4}<K_{\infty}$, the $Nu$ scaling
is like the second peak state.

We shall also comment briefly on the case where $a$ is non-zero. For any $a$,
it is not so difficult to decrease $k$ of the Taylor vortex solution branch
from the linear critical point until the Nusselt number takes a peak. The red
crosses in figure 16a summarise the $Nu$ at the peak for various $a$. The
numerical result shows that the optimised $Nu$ is maximised when the outer
cylinder slightly counter-rotates ($a\approx 0.08$). This is qualitatively in
line with the experimental ultimate turbulence results shown by cyan squares.
However, for the experimental results, the maximum of $Nu$ is attained at
$a\approx 0.3$ and, as noted in section 3, the Nusselt number scaling is
$Nu\propto Ta^{0.38}$. The Taylor vortex results naturally have a much smaller
$Nu$, because it is scaled like $Ta^{0.3}$ as the second peak seen in the
$a=0$ case. In figure 16b, we also plotted $K_{1}$ and $K_{\infty}$ obtained
in the transitional state analysis (see (87) and (91) in Appendix B). Roughly
speaking, the peaks occur when $k/Ta^{1/4}$ is about half of $K_{\infty}$. The
linear and nonlinear instabilities disappear at the Rayleigh line
$a=-\eta^{2}\approx-0.51$, which can be found from Rayleigh’s stability
condition.

Brauckmann et al. (2016) and Brauckmann & Eckhardt (2017) reported that when
the gap is narrow, there are two values of
$R_{\Omega}=(1-\eta)\frac{R_{i}+R_{o}}{R_{i}-\eta R_{o}}$ at which the local
maximum of $Nu$ occurs. Whether similar results can be obtained with the
Taylor vortex solution is left for future work.

##  7 Conclusions and discussion

In this study, the Taylor vortex solution branch is continued to very high
Taylor numbers, and a theoretical rationale is given for its asymptotic
properties. The limiting behaviour of the solution depends on how the axial
period of the solution $2\pi/k$ is scaled by $Ta$. The range of $k$ from
$O(1)$ to the value where it is so large that the nonlinear solution ceases to
exist is covered in the analysis. The Taylor vortex solution with a wavenumber
$k$ of $O(1)$ reproduces surprisingly well the properties of the large-scale
coherent structures in the classical turbulence regime. The numerical results
also suggest that there may be some connection between the onset of the
ultimate turbulence regime and the high wavenumber solutions, although the
precise role of the solutions in the dynamics needs to be clarified in the
future.

Our results provide evidence that the strategy of taking a simple solution
from the dynamics and applying the matched asymptotic expansion for it can be
used even for fully developed high-Reynolds-number turbulence to some extent.
This finding offers great hope for the first-principles explanation for
turbulent momentum and heat transport that is still poorly understood.

###  7.1 Summary of the asymptotic properties of the solutions

When $k=O(1)$, the asymptotic solution consists of an inviscid core and a
viscous boundary layer enveloping it. In the core region, the Prandtl-Bachelor
theorem imposes strong restrictions on the structure of the flow field. The
asymptotic scaling of the wind Reynolds number $Re_{w}=O(Ta^{1/2})$ is
consistent with the turbulent observations (Huisman et al. 2012, Ostilla et
al. 2013). On the other hand, the viscous boundary layer can be divided into a
wall boundary layer and a plume. According to the asymptotic analysis, it is
the Coriolis force acting on the plume that drives the overall flow. The
theoretical Nusselt number for the $k=O(1)$ solutions is $Nu=O(Ta^{1/4})$, and
the DNS data suggests that the same $Nu$ scaling can be applied for the
classical turbulence.

The Nusselt number of the solution reaches its maximum when $k=O(Ta^{2/9})$.
This asymptotic state, which we call the first peak state, has a structure
similar to the $k=O(1)$ flow in a neighbourhood of the walls. The near-wall
boundary layer becomes thinner than the $k=O(1)$ case, so naturally, the
Nusselt number is increased. At the same time, viscous effects in the plume
dominate Coriolis forces. This means that the Coriolis force must act in the
core to drive the vortex, and this condition determines the scaling
$Nu=O(Ta^{1/3})$ and $Re_{w}=O(Ta^{4/9})$. For large wavenumber solutions, the
radial velocity is much larger than the axial velocity, so the scaling of
$Re_{w}$ is the same for the definitions used by Huisman et al. (2012) and
Ostilla et al. (2013).

There is another local maximum for the Nusselt number, which occurs at
$k=O(Ta^{1/4})$. In this second peak state, the nonlinear interaction of the
axially fluctuating component disappears in the core. As a result, only a
single Fourier mode and mean flow play a major role. To connect this core flow
towards the cylinder walls, three near-wall layers are necessary. The bottom
boundary layer adjacent to the cylinder wall is needed in order to satisfy the
no-slip boundary conditions, and the top boundary layer must appear in order
to attenuate the Fourier harmonics towards the core region. These two boundary
layers of different natures can be matched through the middle boundary layer.
The Nusselt number and wind Reynolds number scalings are $O(Ta^{3/10})$ and
$O(Ta^{2/5})$, respectively, ignoring a minor logarithmic correction. The
Coriolis force is a leading-order effect in all regions except the bottom
boundary layer.

The second peak state only appears when the scaled wavenumber $kTa^{-1/4}$ is
smaller than the critical value $K_{\infty}$, which can be determined
analytically. As the wavenumber $k$ is further increased, the Taylor vortex
undergoes a transitional state similar to that analysed by Denier (1992) and
then becomes a laminar flow. The Nusselt number of the transitional state is
$O(Ta^{0})$ and can be found analytically, which agrees well with the
numerical results.

In summary, as $k$ is reduced from the linear critical point, the Taylor
vortex solution develops according to the following scenario. First, vortices
grow from near the inner cylinder in the transitional state. When the vortices
reach the outer cylinder, the flow field becomes the second peak state. A
further decrease in $k$ increases the thickness of the outer boundary layer,
where Fourier harmonics are seen. The first peak state appears when the outer
boundary layers on both sides come into contact. The axial scale of this state
is sufficiently large for the separation of the inviscid region and the plume
to be noticeable in the near-wall regions. When $k=O(1)$, the inviscid region
extends across the gap, and the scaling for the largest-scale vortices is
established.

How the above steady solutions relate to the turbulence dynamics is a natural
question; Kooloth et al. (2021) investigated this problem for two-dimensional
RBC and their results may also be applicable to our case. Their key
observation is that local structures in the turbulent dynamics can be
approximated by the relatively short-wavelength solutions if the wavelength is
suitably chosen. The detailed comparison of the dynamics and the solutions
revealed that local structures with a wide range of wavelengths are embedded
in turbulence. It is still unknown how often the solution of each wavelength
appear in the dynamics. But if all solutions are equally likely to appear,
then the solutions with the largest $Nu$ scaling will have a greatest impact
on the average value of $Nu$ in the dynamics.

An interesting aspect of the above asymptotic theories is that it shows what
the mean flow $\overline{v}$ in the core will look like. When $k$ is not too
large, the mean angular momentum $r\overline{v}$ becomes a constant in the
core region. In the $k=O(1)$ case, this is a consequence of the Prandtl-
Batchelor theorem. The mean angular momentum profile becomes nonuniform when
$k=O(Ta^{1/4})$, but instead, an analytical solution can be derived up to an
unknown constant. The time-averaged turbulence data support the former uniform
angular momentum law (Wattendorf 1935; Taylor 1935; Smith & Townsend 1982;
Lewis & Swinney 1999; Ostilla-Mónico et al. 2014b). The fact that
$\overline{v}$ is proportional to $1/r$ implies that the mean flow is
irrotational. By taking the narrow gap limit, it can be seen that the argument
here explains why zero absolute vorticity states are commonly found in
rotating parallel shear flows (see Johnston et al. 1972, Tanaka et al. 2000,
Suryadi et al. 2014, Kawata et al. 2016).

We further remark that the uniform angular momentum states are a natural
generalisation of the uniform momentum states often observed in non-rotating
shear flows. In particular, the staircase-like uniform momentum zones
ubiquitously seen in the near-wall turbulent boundary layer have long
attracted much attention (Meinhart & Adrian 1995; de Silva, Hutchins & Marusic
2016). To explain this phenomenon, Montemuro et al. (2020) and Blackburn et
al. (2021) attempted to construct a theory for homogeneous shear flows. At the
heart of these theories is the Prandtl-Batchelor theorem, which was first
proposed to apply for the uniform momentum state by Deguchi & Hall (2014a) for
a steady solution of plane Couette flow. It is however not yet clear how such
a core structure interacts with the near wall flow and deduces the scaling of
the momentum transport in Reynolds number.

###  7.2 Link to the RBC studies

The asymptotic analyses of the first and second peaks presented in this paper
can also be applied to the roll-cell solutions of RBC (regardless of whether
the boundary is slip or no-slip) and therefore provide a theoretical
explanation for the scaling obtained by the numerical computation in Waleffe
et al. (2015) and Sondak et al. (2015). Note, however that for RBC, the
transitional state only exists in the vicinity of the bifurcation point in the
parameter space (Blennerhassett & Bassom 1994). This is due to the fact that
in RBC, the instability occurs uniformly in the flow, whereas in TCF, the
centrifugal instability is usually most pronounced near the inner cylinder.

The computation of the roll-cell solution branch by Wen et al. (2022) reaches
$Ra=10^{14}$, which is close to the experimentally feasible limit. In their
numerical result the exponent for $k$ at the first peak is closer to 1/5 than
2/9. However, their $Nu$ exponent is also slightly off from 1/3 so a
calculation with a larger $Ra$ would be necessary to obtain the exponents
accurately. Of course, the possibility that an unknown asymptotic state exists
cannot be ruled out, but it appears from the author’s experience that it is
difficult to construct new sensible asymptotic theories.

It is worth mentioning that Iyer et al. (2020) obtained the Nusselt number
exponents 0.29 and 0.331 for moderate and high Rayleigh number regimes,
respectively, for RBC in a vertically elongated tank. These exponents are
close to those of the second and first peaks. The scaling of the first peak
matches with the so-called classical scaling, but our derivation significantly
differs from that of Priestley (1954), Malkus (1954), Grossman & Lohse (2000).
Recently, Kawano et al. (2021) derived the scaling laws by considering the
naive balance of each term in the governing equations within the boundary
layer. Although they do not assume the two-dimensionality of the flow, the
scaling is consistent with our first peak asymptotic state. In particular, the
wind Reynolds number scaling $Re_{w}=O(Ta^{4/9})$ we found is in agreement
with the results by Grossman & Lohse (2000) and Kawano et al. (2021). As
pointed out by the latter paper, the exponent $4/9\approx 0.444$ is close to
0.458 observed in the high Rayleigh number DNS by Iyer et al. (2020). Their
wind Reynolds number is defined by the root-mean square of the three velocity
components normalised by viscous velocity scale $\nu/H$, where $H$ is the
vertical length of the tank, and thus comparable with our results. Iyer et al.
(2020) used a computational domain where the horizontal scale is 1/10 of the
vertical scale, and it therefore makes sense that the result is related to our
large $k$ solutions.

The exponent of the Nusselt number obtained in Iyer et al. (2020) is smaller
than the 0.38 obtained by He et al. (2012). The tank aspect ratio might be one
of the causes of this difference. However, there are only a few reports of
$Nu$ exponents exceeding 1/3 in RBC, and the true nature of ultimate RBC
scaling is still a matter of active debate. Wen et al. (2022) summarised the
$Nu$ data of turbulent RBC and found that they were all smaller than those of
the optimumised steady roll cell solution corresponding to the first peak
state. Zhu et al. (2018) restricted the DNS of RBC to two dimensions and
increased the Rayleigh number to $10^{14}$ hoping to see the ultimate
turbulence scaling. However, their numerical data are scaled by the exponent
$\beta=1/3$, as pointed out by Doering et al. (2019). Even without restricting
to two-dimensional steady states, the existence of simple solutions with an
exponent greater than $1/3$ is still unknown. Motoki et al. (2018) computed a
three-dimensional steady solution of RBC, but the $Nu$ scaling was similar to
that seen in Waleffe et al. (2015) and Sondak et al. (2015).

###  7.3 What is missing in the theoretical results?

Unlike RBC, the emergence of the Nusselt number exponent $\beta=0.38$ is well-
established in the ultimate turbulence regime of TCF. Therefore, the
predictions by the first peak state should eventually deviate from the
turbulent measurement as the Taylor number increases. In the ultimate
turbulence, a mix of large-scale and small-scale structures was observed, so
in a sense, it is natural that this difference should appear. Experiments have
shown that large-scale structures with scaling of $Re_{w}=Ta^{1/2}$ can be
observed even in the ultimate turbulence regime (see Huisman et al. 2012).
This suggests that the $k=O(1)$ state studied in section 4 still plays some
role in the dynamics, but as we have seen in our analysis, it cannot
efficiently transport angular momentum. The reason for the inefficiency is the
Prandtl-Batchelor homogenisation of $\Gamma$ as remarked just below (4.15).
Thus it would be a natural idea to introduce small-scale vortices in the core
to enhance transport. Note however that vortices of about Kolmogorov
microscale $Re_{w}^{-1/2}$ are unlikely to show any improvement. The typical
velocity scale $Re_{w}^{1/2}$ (see Deguchi 2015, for example) yields
$O(\widetilde{\Gamma}|_{c})=O(Ta^{-1/2}Re_{w}^{1/2})$, so from the transport
balance $Nu=O(Re_{w}\widetilde{\Gamma}|_{c})$ derived by (4.15) the scaling
remains $Nu=O(Ta^{1/4})$. On the other hand, the azimuthal velocity
perturbations generated by the first and second peak states are larger and may
improve transport. If instead $O(\widetilde{\Gamma}|_{c})=O(Ta^{-1/10})$
obtained in the second peak analysis is used, the resultant scaling
$Nu=O(Ta^{0.4})$ is close to the experimental observation, though of course it
is not known whether such a flow structure is possible in the asymptotic
expansion framework.

Another possible reason for the deviation of the $Nu$ scaling would be the
axisymmetric restriction imposed for the Taylor vortex. When the three-
dimensionality of TCF is allowed, feedback through Reynolds stresses appears
from the non-axisymmetric component to the axisymmetric component. This
feedback effect is a key process in self-sustaining the coherent structures in
non-rotating shear flows (Waleffe 1997), and recently Dessup et al. (2018),
Sacco et al. (2019) attempted to extend this idea to TCF. Our results suggest
that the feedback is unimportant for classical turbulence, though it may not
be so for ultimate turbulence.

Minor modifications to the theories obtained in this paper using the existing
asymptotic results would not fully explain the scaling of the ultimate
turbulence. For example, azimuthal and time-dependent waves can be
incorporated into our axisymmetric asymptotic theories, using the method used
in Hall & Sherwin (2010), but this does not change the $Nu$ scaling. The crux
of this problem would be that all existing nonlinear asymptotic theories for
shear flows do not much change the scaling of the wall shear rate from that
for laminar flows.

Experiments and numerical simulations have shown that there is universality in
the friction factor of various high Reynolds number shear flows (see Orlandi
et al. (2015) for example). No rational asymptotic theory has been found to
explain this friction factor. If a new simple solution that reproduces the
universal friction factor scaling can be found, it might serve an important
hint for understanding the ultimate turbulence of TCF as well. For example, if
one assumes the empirical Blasius friction law (friction factor $\propto
Re^{-1/4}$), it implies the emergence of near wall boundary layer of thickness
$O(Re^{-3/4})$. In terms of the Taylor number the boundary layer thickness is
$O(Ta^{-3/8})$, so the Nusselt number exponent is $3/8=0.375$, which is not
too far from 0.38. The similarities between the friction factor of TCF and
other shear flows are also noted by Lathrop et al. (1992). If the hypothesis
that the ultimate turbulence of TCF shares a common mechanism with that of
non-rotating shear flows is true, then this could settle the debate on whether
the TCF-RBC analogy is valid in the extreme parameters.

Finally, we remark that the asymptotic structure discussed in this paper is
valid when $R_{i}$ is smaller than or equal to $O(|R_{o}|)$, and hence it does
not cover the whole parameter space. When the cylinders are counter-rotating,
on the neutral curve the scaling $R_{i}=O(|R_{o}|^{5/3})$ is established, as
noted by Donnelly & Fultz (1960) using a dimensional argument, by Esser &
Grossman (1996) using a heuristic extension of Rayleigh’s stability condition,
and by Deguchi (2016) using an asymptotic analysis. It is well-known that non-
axisymmetric disturbances are important near the neutral curve and nonlinear
patterns such as spirals and ribbons emerge. Furthermore, some solution
branches can be continued in the subcritical parameter regime, and in this
case non-axisymmetry of the flow is critically important to sustain the
nonlinear flow; see Deguchi et al. (2014), Wang et al. (2022). Developing a
nonlinear asymptotic theory capable of describing the above non-axisymmetric
phenomena is another interesting piece of future work.

The author thanks Dr Ostilla-Mónico for sharing his DNS data, and Dr David
Goluskin for inspiring discussion. This work was supported by Australian
Research Council Discovery Project DP230102188.

The author reports no conflict of interest.

##  Appendix A The narrow-gap limits

First, we shall derive the narrow gap limit of the Görtler type. Here we focus
on the stationary outer cylinder case studied by Denier (1992); see Deguchi
(2016) for general cases. Before taking the limit, we rewrite the system (6)
using $V=v/R_{i}$ and $y=r-r_{i}$. The boundary conditions become

$\displaystyle(u,V,w)=(0,0,0)\qquad\text{at}\qquad y=1,$ (77)
$\displaystyle(u,V,w)=(0,1,0)\qquad\text{at}\qquad y=0.$ (78)

While taking the limit $\eta\rightarrow 1$, or equivalently
$r_{i}\rightarrow\infty$, we assume that the transformed velocity components
and its derivatives are all $O(r_{i}^{0})$. Keeping $T=2R_{i}^{2}/r_{i}$ as
$O(r_{i}^{0})$ to balance the Coriolis term, the limiting flow is governed by

$\displaystyle DV=\triangle V,\qquad D\omega=\triangle\omega+TV\partial_{z}V,$
(79)

and $\partial_{y}u+\partial_{z}w=0$. Here
$\omega=\partial_{z}u-\partial_{y}w,D=\partial_{t}+u\partial_{y}+w\partial_{z},\triangle=\partial_{y}^{2}+\partial_{z}^{2}$.

For the rotating plane Couette flow (RPCF) limit, on the other hand, the
radial coordinate is shifted as $y=r-r_{m}$ so that the origin is at the mid
gap. Furthermore, we take the new azimuthal velocity $V$ to satisfy
$(v-v_{b})=G(V+y)$, which means that a constant $G$ scales the disturbance.
Here we expect that the base flow of the transformed system becomes almost a
linear profile when the gap is narrow. The boundary conditions become

$\displaystyle(u,V,w)=(0,-1/2,0)\qquad\text{at}\qquad y=1/2,$ (80)
$\displaystyle(u,V,w)=(0,1/2,0)\qquad\text{at}\qquad y=-1/2.$ (81)

The constant $G$ is chosen as $G=-(rv_{b})^{\prime}/r$, and the key assumption
to take the system to RPCF is $O(G)\ll O(v_{b})$. Keeping
$T=2Gv_{b}(r_{m})/r_{m}$ as $O(r_{m}^{0})$, it is straightforward to show that
in the narrow gap limit, the system (6) reduces to

$\displaystyle DV=\triangle V,\qquad D\omega=\triangle\omega+T\partial_{z}V.$
(82)

and $\partial_{y}u+\partial_{z}w=0$. The derivation of RPCF here is
mathematically simpler than the common one (see Drazin & Reid 1981), though
the physical motivation would be a bit harder to see.

The assumption $O(G)\ll O(v_{b})$ used in the RPCF limit is equivalent to
$O(R_{o}-R_{i})\ll O(R_{i})$, meaning that the inner and outer cylinders
rotate at approximately the same speed. Thus when considering a coordinate
system rotating at the average speed of the cylinders, the Coriolis force
acting on the fluid is uniform. The uniform force produces an equivalent
effect to buoyancy, allowing us to use a perfect correspondence with RBC. Of
course, in the general case this force is not uniform, which is why in the
Denier’s case the instability grew from near the inner cylinder (see figure 14
also).  Note that the forcing term in the general case corresponds to the
third term on the right-hand side of (4.1b), and strictly speaking this is not
a Coriolis force anymore. In fact, the definition of the Coriolis force
changes depending on which rotational coordinate system is chosen, and hence
it cannot be well defined when the inner and outer cylinders differentially
rotate. The forcing terms are produced by the change in the direction of the
base flow (i.e. the linear terms in the $r^{-1}v^{2}$ and $r^{-1}uv$ terms in
(2.1)). They are not centrifugal force either - centrifugal force can be
included into the pressure gradient term and does not affect the dynamics.

The difference between the two limits also affects the azimuthal development
of the flow, if any. For the Görtlier type limit, the scaling factor of the
azimuthal velocity is large, and therefore the flow must have a larger length
scale than the gap in the azimuthal direction. In contrast, for the RPCF type
limit, we can make the scaling factor $G$ to be $O(1)$ when $R_{i}=O(r_{m})$.
In this case, the azimuthal length scale of the flow is comparable to the gap.

##  Appendix B The transitional states with outer cylinder rotation

This section discusses what modifications are required if $\Gamma_{o}\neq 0$
in section 6.1. The effect of the outer cylinder rotation appears in the mean
flow layer solution so that (36) must be replaced by

$\displaystyle\overline{\Gamma}=B(r_{o}^{2}-r^{2})+\Gamma_{o}+\cdots,\qquad
N_{0}=2(r_{o}^{2}B+\Gamma_{o}).$ (83)

Noting that in the core we can still use (39) and (40), the continuity of
$\overline{\Gamma}$ and $r^{3}(r^{-2}\overline{\Gamma})_{r}$ at $r=r_{c}$ is
satisfied if

$\displaystyle\sqrt{\frac{A-r_{c}^{4}}{4T_{0}}}=B(r_{o}^{2}-r_{c}^{2})+\Gamma_{o},$
(84)
$\displaystyle\frac{A}{\sqrt{T_{0}(A-r_{c}^{4})}}=2(Br_{o}^{2}+\Gamma_{o}).$
(85)

Eliminating $B$ from those equations we get

$\displaystyle\frac{A}{r_{o}^{2}}-r_{c}^{2}=\frac{2\sqrt{T_{0}}}{r_{o}^{2}}\sqrt{A-r_{c}^{4}}\Gamma_{o}$
(86)

which links $r_{c}^{2}$ and $T_{0}$.

Now let us suppose $r_{c}=r_{i}$ when $T_{0}=T_{1}$ in the above equation.
Then the scaled wavenumber at the linear critical point can be easily found as

$\displaystyle K_{1}=T_{1}^{-1/4},\qquad
T_{1}=\frac{r_{i}^{2}(r_{o}^{2}-r_{i}^{2})}{4(\Gamma_{i}^{2}-\Gamma_{i}\Gamma_{o})}.$
(87)

Likewise setting $r_{c}=r_{o}$ and $T_{0}=T_{\infty}$ in (86),

$\displaystyle\sqrt{4T_{\infty}\Gamma_{i}^{2}+r_{i}^{4}-r_{o}^{4}}\left(\sqrt{4T_{\infty}\Gamma_{i}^{2}+r_{i}^{4}-r_{o}^{4}}-\sqrt{4T_{\infty}}\Gamma_{o}\right)=0.$
(88)

From this equation the critical scaled wavenumber where the $Nu$ blows up can
be found as

$\displaystyle K_{\infty}=T_{\infty}^{-1/4},\qquad
T_{\infty}=\left\\{\begin{array}[]{c}\frac{r_{o}^{4}-r_{i}^{4}}{4(\Gamma_{i}^{2}-\Gamma_{o}^{2})}~{}~{}~{}\text{if}~{}~{}~{}\Gamma_{i}>\Gamma_{o}>0,\\\
\frac{r_{o}^{4}-r_{i}^{4}}{4\Gamma_{i}^{2}}~{}~{}~{}\text{if}~{}~{}~{}\Gamma_{o}\leq
0.\end{array}\right.$ (91)

Note that equation (88) may have two roots, but by assuming that
$K_{\infty}(a)=T_{\infty}^{-1/4}$ is a continuous function and that
$K_{\infty}\rightarrow 0$ as the Rayleigh line ($\Gamma_{i}=\Gamma_{o}$) is
approached we can uniquely determine the solution (91).

To find $N_{0}$, we use the fact that the quadratic equation

$\displaystyle
0=\left(\frac{r_{o}^{4}}{A}-1\right)N_{0}^{2}+4\Gamma_{o}N_{0}-\left(4\Gamma_{o}^{2}+\frac{r_{o}^{4}}{T_{0}}\right)$
(92)

can be obtained by combining (84), (85), (86). The solution

$\displaystyle
N_{0}=2\frac{\sqrt{\Gamma_{o}^{2}+(A^{-1}r_{o}^{4}-1)(\Gamma_{o}^{2}+\frac{r_{o}^{4}}{4T_{0}})}-\Gamma_{o}}{A^{-1}r_{o}^{4}-1}$
(93)

and (35) yields

$\displaystyle
Nu(T_{0})=\frac{(1+\eta)^{4}}{16\eta^{3}\sqrt{T_{0}}}\frac{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}{r_{o}^{4}-r_{i}^{4}-4T_{0}\Gamma_{i}^{2}}\left(\sqrt{\frac{r_{o}^{4}+4T_{0}\Gamma_{o}^{2}}{r_{i}^{4}+4T_{0}\Gamma_{i}^{2}}-1}-\Gamma_{o}\right),$
(94)

which is the generalised version of (53).

## References

  * [] 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–183.
  * [] Bassom, A. P. & Blennerhassett, P. J. 1992 The structure of highly nonlinear vortices in curved channel flow. _Proc. R Soc. Lond._ A 439 , 317–336.
  * [] Bassom, A. P. & Hall, P. 1989 On the generation of mean flows by the interaction of Görtler vortices and Tollmien-Schlichting waves in curved channel flows. _Stud. Appl. Math._ 81 , 185–219.
  * [] Batchelor, G. K. 1956 On steady laminar flow with closed streamlines at large Reynolds number. _J. Fluid Mech._ 1(2), 177–190.
  * [] Blackburn, H. M., Deguchi, K. & Hall, P. 2021 Distributed vortex-wave interactions: the relation of self-similarity to the attached eddy hypothesis. _J. Fluid Mech._ 924, A8.
  * [] Blennerhassett, P. J. & Bassom, A. P. 1994 Nonlinear high-wavenumber Bénard convection. _IMA J. Appl. Math._ 52 , 51–77.
  * [] Bradshaw, P. 1969 The analogy between streamline curvature and buoyancy in turbulent shear flow. _J. Fluid Mech._ 36, 177–191.
  * [] Brauckmann, H. J. & Eckhardt, B. 2013 Direct numerical simulations of local and global torque in Taylor-Couette flow up to $Re=30000$. _J. Fluid Mech._ 718, 398–427.
  * [] Brauckmann, H. J. & Eckhardt, B. 2017 Marginally stable and turbulent boundary layers in low-curvature Taylor–Couette flow. _J. Fluid Mech._ 815, 149–168.
  * [] Brauckmann, H. J., Salewski, M. & Eckhardt, B. 2016 Momentum transport in Taylor-Couette flow with vanishing curvature. _J. Fluid Mech._ 790, 419–452.
  * [] Chandrasekhar, S. 1961 _Hydrodynamic and hydromagnetic stability_. Oxford University Press.
  * [] Chini, G. P. & Cox, S. M. 2009 Large Rayleigh number thermal convection: heat flux predictions and strongly nonlinear solutions. _Phys. Fluids_ 21, 083603.
  * [] Cvitanović, P., Artuso, R., Mainieri, R., Tanner, G. & Vattay, G. 2012 Chaos: Classical and Quantum. Niels Bohr Inst., www.chaosbook.org.
  * [] Davey, A. 1962 The growth of Taylor vortices in flow between rotating cylinders. _J. Fluid Mech._ 14, 336–368.
  * [] Deguchi, K. 2015 Self-sustained states at Kolmogorov microscale. _J. Fluid Mech._ 781, R6.
  * [] Deguchi, K. 2016 The rapid-rotation limit of the neutral curve for Taylor-Couette flow. _J. Fluid Mech._ 808, R2.
  * [] Deguchi, K. & Altmeyer 2013 Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow. _Phys. Rev. E_ 87(4), 043017.
  * [] Deguchi, K. & Hall, P. 2014a The high Reynolds number asymptotic development of nonlinear equilibrium states in plane Couette flow. _J. Fluid Mech._ 750, 99–112.
  * [] Deguchi, K. & Hall, P. 2014b Free-stream coherent structures in parallel boundary-layer flows. _J. Fluid Mech._ 752, 602–625.
  * [] Deguchi, K., Hall, P. & Walton, A. G. 2013 The emergence of localized vortex-wave interaction states in plane Couette flow. _J. Fluid Mech._ 721, 58–85.
  * [] Deguchi, K., Meseguer, A. & Mellibovsky, F. 2014 Subcritical equilibria in Taylor-Couette flow. _Phys. Rev. Lett._ 112(18), 184502.
  * [] Deguchi, K. & Walton, A. G. 2018 Bifurcation of nonlinear Tollmien-Schlichting waves in a high-speed channel flow. _J. Fluid Mech._ 843, 53–97.
  * [] Dempsey, L. J., Deguchi, K., Hall, P. & Walton, A. G. 2016 Localized vortex/Tollmien-Schlichting wave interaction states in plane Poiseuille flow. _J. Fluid Mech._ 791, 97–121.
  * [] Denier, J. P. 1992 The structure of fully nonlinear Taylor vortices. _IMA J. Appl. Math._ 49, 15–33.
  * [] Dennis, P. M., Huisman, S. G., Bruggert, G., Sun, C. & Lohse, D. 2011 Torque scaling in turbulent Taylor-Couette flow with co- and counterrotating cylinders. _Phys. Rev. Lett._ 106, 024502.
  * [] Dessup, T., Tuckerman, L. S., Wesfreid, J. E., Barkley, D. & Willis, A. P. 2018 The self-sustaining process in Taylor–Couette flow. _Phys. Rev. Fluid_ 3(12), 123902.
  * [] Ding, Z. & Marensi, E. 2019 Upper bound on angular momentum transport in Taylor-Couette flow. _Phys. Rev. E_ 100, 063109.
  * [] Doering, C. R., Toppaladoddi, S., & Wettlaufer, J. S. 2019 Absence of evidence for the ultimate regime in two-dimensional Rayleigh-Bénard convection. _Phys. Rev. Lett._ 123, 259401.
  * [] Dong, S. 2007 Direct numerical simulation of turbulent Taylor-Couette flow. _J. Fluid Mech._ 587, 373–393.
  * [] Donnelly, R. J. & Fultz, D. 1960 Experiments on the stability of viscous flow between rotating cylinders. II. Visual observations. _Proc. Roy. Soc._ A 258, 101–123.
  * [] Drazin, P. G. & Reid W. H. 1981 _Hydrodynamic Stability._ Cambridge University Press.
  * [] Dubrulle, B. & Hersant, F. 2002 Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection. _Eur. Phys. J._ B 26, 379–386.
  * [] Dubrulle, B., Dauchot, O., Daviaud, F., Longgaretti, P. Y., Richard, D. & Zahn, J. P. 2005 Stability and turbulent transport in Taylor-Couette flow from analysis of experimental data. _Phys. Fluids_ 17, 095103.
  * [] Eckhardt, B., Grossmann, S. & Lohse, D. 2007 Torque scaling in turbulent Taylor-Couette flow between independently rotating cylinders. _J. Fluid Mech._ 581, 221–250.
  * [] Esser, A. & Grossmann, S. 1996 Analytic expression for Taylor-Couette stability boundary. _Phys. Fluids_ 8, 1814–1819.
  * [] Feynman, R. P., & Lagerstrom, P. A. 1956 Remarks on high Reynolds number flows in finite domains. _Proc. IX International Congress on Applied Mechanics_ 3 342–343.
  * [] Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. _J. Fluid Mech._ 407, 27–56.
  * [] Grossmann, S., Lohse, D. & Sun, C. 2016 High-Reynolds number Taylor–Couette turbulence. _Annu. Rev. Fluid Mech._ 48, 53–80.
  * [] Hall, P. 1982 Taylor–Görtler vortices in fully developed or boundary layer flows: Linear theory. _J. Fluid Mech._ 124, 475–494.
  * [] Hall, P. 1983 The linear development of Görtler vortices in growing boundary layers. _J. Fluid Mech._ 130, 41–58.
  * [] Hall, P. & Lakin, W. D. 1988 The fully nonlinear development of Görtler vortices in growing boundary layers. _Proc. R. Soc. Lond._ A 415 , 421–444.
  * [] 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.
  * [] He, X., Funfschilling, D., Nobach, H., Bodenschatz, E. & Guenter, A. 2012 Transition to the ultimate state of turbulent Rayleigh-Bénard convection. _Phys. Rev. Lett._ 108, 024502.
  * [] Hepworth, B. J. 2014 Nonlinear two-dimensional Rayleigh Bénard convection. PhD thesis, Department of Applied Mathematics, University of Leeds, UK.
  * [] Huisman, S. G., van Gils, D. P. M., Grossmann, S., Sun, C. & Lohse, D. 2012 Ultimate turbulent Taylor-Couette flow. _Phys. Rev. Lett._ 108, 024501.
  * [] Iyer, K. P., Scheel, J. D., Schumacher, J. & , Sreenivasan, K. R. 2020 Classical 1/3 scaling of convection holds up to $Ra=10^{15}$. _PNAS_ 117, 7594-7598.
  * [] Jimenez, J. & Zufiria, J. A. 1987 A boundary-layer analysis of Rayleigh-Benard convection at large Rayleigh number. _J. Fluid Mech._ 178, 53–71.
  * [] Johnston, J. P., Halleen, M. & Lezius, D. K. 1972 Effects of spanwise rotation on the structure of two-dimensional fully developed turbulent channel flow. _J. Fluid Mech._ 56, 533–557.
  * [] Kawahara, G., Uhlmann, M. & van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. _Ann. Rev. Fluid Mech._ 44, 203-225.
  * [] Kawata, T. & Alfredsson, P. H. 2016 Experiments in rotating plane Couette flow – momentum transport by coherent roll-cell structure and zero-absolute vorticity state. _J. Fluid Mech._ 791, 191–213.
  * [] Kawano, K., Motoki, M., Shimizu, M. & Kawahara, G. 2021 Ultimate heat transfer in ‘wall-bounded’ convective turbulence _J. Fluid Mech._ 914, A13.
  * [] Kooloth, P., Sondak, D. & Smith, L. M. 2021 Coherent solutions and transition to turbulence in two-dimensional Rayleigh-Bénard convection. _Phys. Rev. Fluids_ 6, 013501.
  * [] 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.
  * [] Lathrop, D. P., Fineberg, J. & Swinney, H. L. 1992 Turbulent flow between concentric rotating cylinders at large Reynolds number. _Phys. Rev. Lett._ 68(10), 1515–1518.
  * [] Lewis, G. S., & Swinney, H. L. 1999 Velocity structure functions, scaling, and transitions in high-Reynolds-number Couette-Taylor flow. _Phys. Rev. E_ 59, 5457–5467.
  * [] Lezius, D., & Johnston, J. P. 1976 Roll-cell instabilities in rotating laminar and turbulent channel flows. _J. Fluid Mech._ 77, 153–175.
  * [] Malkus, W. V. R. 1954 The heat transport and spectrum of thermal turbulence. _Proc. R. Soc. Lond._ A 225 (1161), 196–212.
  * [] Meinhart, C. D. & Adrian, R. J. 1995 On the existence of uniform momentum zones in a turbulent boundary layer. _Phys. Fluids_ 7(7), 694–696.
  * [] Montemuro, B., White, C. M., Klewicki, J. & Chini, G. P. 2020 A self-sustaining process theory for uniform momentum zones and inertial shear layers in high Reynolds number shear flows. _J. Fluid Mech._ 901, A28.
  * [] Moore, D. R. & Weiss, N. O. 1973 Two-dimensional Rayleigh-Bénard convection. _J. Fluid Mech._ 58, 289–312.
  * [] Motoki, S., Kawahara, G. & Shimizu, M. 2021 Multi-scale steady solution for Rayleigh-Bénard convection. _J. Fluid Mech._ 914, A14.
  * [] Orlandi, P., Bernardini, M. & Pirozzoli, S. 2015 Poiseuille and Couette flows in the transitional and fully turbulent regime. _J. Fluid Mech._ 770, 424–441.
  * [] Ostilla, R., Stevens, R. J. A. M., Grossmann, S., Verzicco, R. & Lohse, D. 2013 Optimal Taylor-Couette flow: direct numerical simulations. _J. Fluid Mech._ 719, 14–46.
  * [] Ostilla-Mónico, R., van der Poel, E. P., Verzicco, R., Grossmann, S. & Lohse, D. 2014a Boundary layer dynamics at the transition between the classical and the ultimate regime of Taylor-Couette flow. _Phys. Fluids_ 26, 015114.
  * [] Ostilla-Mónico, R., van der Poel, E. P., Verzicco, R., Grossmann, S. & Lohse, D. 2014b Exploring the phase diagram of fully turbulent Taylor-Couette flow. _J. Fluid Mech._ 761, 1–26.
  * [] Paoletti, M. S., & Lathrop, D. P. 2011 Angular momentum transport in turbulent flow between independently rotating cylinders. _Phys. Rev. Lett._ 106, 024501.
  * [] Priestley, C. H. B. 1954 Convection from a large horizontal surface. _Austral. J. Phys._ 7(1), 176–201.
  * [] Pillow, A. F. 1952 _Aust. Dep. Supply Aero. Res. Lab. Rept._ A79.
  * [] Prandtl, L. 1904 Uber flussigkeits bewegung bei sehr kleiner reibung. _Verhaldlg III Int. Math. Kong_ 484–491.
  * [] Roberts, G. O. 1979 Fast viscous Bénard convection. _Geophys. Astrophys. Fluid Dynamics_ 12, 235–272.
  * [] Robinson, J. L. 1967 Finite amplitude convection cells. _J. Fluid Mech._ 30(3), 577–600.
  * [] Sacco, F., Verzicco, R. & Ostilla-Mónico, R. 2019 Dynamics and evolution of turbulent Taylor rolls. _J. Fluid Mech._ 870, 970–987.
  * [] Smith, G. P. & Townsend, A. A. 1982 Turbulent Couette flow between concentric cylinders at large Taylor numbers. _J. Fluid Mech._ 123, 187–217.
  * [] de Silva, C. M., Hutchins, N. & Marusic, I. 2016 Uniform momentum zones in turbulent boundary layers. _J. Fluid Mech._ 786, 309–321.
  * [] Sondak, D., Smith, M. L. & Waleffe, F. 2015 Optimal heat transport solutions for Rayleigh-Bénard convection. _J. Fluid Mech._ 784, 565–595.
  * [] Suryadi, A., Segalini, A. & Alfredsson, P. H. 2014 Zero absolute vorticity: insight from experiments in rotating laminar plane Couette flow. _Phys. Rev. E_ 89, 033003.
  * [] Tanaka, M., Kida, S., Yanase, S. & Kawahara, G. 2000 Zero-absolute -vorticity state in a rotating turbulent shear flow. _Phys. Fluids_ 12, 1979–1985.
  * [] Taylor, G. I. 1923 Stability of a viscous liquid contained between two rotating cylinders. _Phil. Trans. R. Soc. Lond. A_ 223, 289–343.
  * [] Taylor, G. I. 1935 Distribution of velocity and temperature between concentric rotating cylinders. _Proc. R. Soc. Lond. A_ 151, 494–512.
  * [] van Gils, D. P. M., Huisman, S. G., Bruggert, G. W., Sun, C. & Lohse, D. 2011 Torque scaling in turbulent Taylor-Couette flow with co- and counter-rotating cylinders. _Phys. Rev. Lett._ 106, 024502.
  * [] van Gils, D. P. M., Huisman, S. G., Grossmann, S., Sun, C. & Lohse, D. 2012 Optimal Taylor-Couette turbulence. _J. Fluid Mech._ 706, 118–149.
  * [] Veronis, G. 1970 The analogy between rotating and stratified fluids. _Annu. Rev. Fluid Mech._ 2, 37–66.
  * [] Vynnycky, M. & Masuda, Y. 2013 Rayleigh-Benard convection at high Rayleigh number and infinite Prandtl number: asymptotics and numerics. _Phys. Fluids_ 25, 113602.
  * [] Waleffe, F. 1997 On a self-sustaining process in shear flows. _Phys. Fluids_ 9, 883–900.
  * [] Waleffe, F., Boonkasame, A. & Smith, L. 2015 Heat transport by coherent Rayleigh-Bénard convection. _Phys. Fluids_ 27(5), 051702.
  * [] Wang, B., Ayats, R., Deguchi, K., Mellibovsky, F. & Meseguer, A. 2022 Self-sustainment of coherent structures in counter-rotating Taylor-Couette flow. _J. Fluid Mech._ 951, A21.
  * [] Wattendorf, F. L. 1935 A study of the effect of curvature on fully developed turbulent flow. _Proc. R. Soc. London A_ 148(5), 565–598.
  * [] Wen, B., Goluskin, D. & Doering, C. R. 2022 Steady Rayleigh-Bénard convection between no-slip boundaries. _J. Fluid Mech._ 933, R4.
  * [] Wesseling, P. 1969 Laminar convection cells at high Rayleigh number. _J. Fluid Mech._ 36(4), 625–637.
  * [] Willis, A. P., Cvitanović, P.& M. Avila 2013 Revealing the state space of turbulent pipe flow by symmetry reduction. _J. Fluid Mech._ 721, 514–540.
  * [] Zhu, X., Mathai, V., Stevens, R. J. A. M., Verzicco, R. & Lohse, D. 2018 Transition to the ultimate state of turbulent Rayleigh-Bénard convection. _Phys. Rev. Lett._ 128, 144502.