1. Introduction
Taylor–Couette (TC) flow, the flow in between two coaxial, independently rotating cylinders, has successfully been used as a model for shear flows to study instabilities, flow patterns, nonlinear dynamics and transitions and turbulence (Taylor Reference Taylor1923; Chandrasekhar Reference Chandrasekhar1981; Andereck, Liu & Swinney Reference Andereck, Liu and Swinney1986; Lewis & Swinney Reference Lewis and Swinney1999; van Gils et al. Reference van Gils, Bruggert, Lathrop, Sun and Lohse2011; Paoletti & Lathrop Reference Paoletti and Lathrop2011; Fardin, Perge & Taberlet Reference Fardin, Perge and Taberlet2014; Ostilla-Mónico et al. Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016). The basic TC geometry is characterized by two parameters: the first is the radius ratio 
$\eta =r_i/r_o$, where 
$r_i$ and 
$r_o$ are the inner and outer radii, respectively. The second is the aspect ratio 
$\varGamma =L/d$, where 
$L$ is the height of the cylinders and 
$d=r_o-r_i$ is the width of the gap. The shear driving of the flow is produced by the cylinders differential rotation and, in dimensionless form, is expressed by the Taylor number (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007)
\begin{equation} {\textit{Ta}}=\frac{(1+\eta)^4}{64\eta^2}\frac{d^2(r_o+r_i)^2(\omega_i-\omega_o)^2}{\nu^2}, \end{equation}
where 
$\omega _{i,o}$ are the inner and outer angular velocities, respectively, and 
$\nu$ is the kinematic viscosity of the fluid. The second control parameter is the rotation ratio
\begin{equation} a={-}\omega_o/\omega_i, \end{equation}
where 
$a<0$ denotes corotation of the cylinders while 
$a>0$ indicates counter-rotating cylinders. The value of 
$a=0$ corresponds to the case of pure inner cylinder rotation.
We note that, instead of describing the control parameters of TC flow with 
${\textit {Ta}}$, 
$\eta$ and 
$a$, one could alternatively describe the parameter space in a convective reference frame as proposed by Dubrulle et al. (Reference Dubrulle, Dauchot, Daviaud, Longaretti, Richard and Zahn2005) such that the cylinders rotate with opposite velocities 
${\pm } U/2$ and the entire system then rotates with angular velocity 
$\boldsymbol{\varOmega} = \varOmega _{rf} \boldsymbol {e}_{\boldsymbol {z}}$ around the central axis. Here, 
$U$ is the characteristic velocity 
$U=2(u_i-\eta u_o)/(1+\eta )$, 
$\varOmega _{rf}=(r_i\omega _i +r_o\omega _o)/(r_i+r_o)$ is the mean angular velocity, 
$\boldsymbol {e}_{\boldsymbol {z}}$ is the unit vector in the axial direction and 
$u_{i,o}=r_{i,o}\omega _{i,o}$ are the inner and outer cylinder streamwise velocities. This way, any combination of differential rotations of the cylinders is parametrized as a Coriolis force. In this frame the two control parameters are the shear Reynolds number 
${\textit {Re}}_S$ for the driving strength, the curvature number 
$R_C$ and the rotation number 
$R_\varOmega$
\begin{gather} Re_S = \frac{U {d}}{\nu} = 2r_ir_o \frac{{d}|\omega_i-\omega_o|}{(r_i+r_o)\nu}, \end{gather}
\begin{gather}R_C=\frac{(1-\eta)}{\sqrt{\eta} }, \end{gather}
\begin{gather}R_\varOmega = \frac{2\varOmega_{rf} d}{U} =(1-\eta)\frac{r_i \omega_i +r_o \omega_o}{r_i \omega_i - r_i \omega_o}. \end{gather} We remark that 
${\textit {Re}}_S\propto \sqrt {{\textit {Ta}}}$, and that the rotation number 
$R_\varOmega$ is connected with the negative rotation ratio 
$a$ by
\begin{equation} R_\varOmega=(1-\eta)\frac{1-a/\eta}{1+a}. \end{equation} While the choice of one set of parameters might seem arbitrary at first, we note that 
$R_\varOmega$ is the quantity that controls the magnitude of the Coriolis force when the equations are written in the rotating reference frame, and it becomes particularly relevant to elucidate certain effects, especially in the limit of low curvature (Brauckmann, Salewski & Eckhardt Reference Brauckmann, Salewski and Eckhardt2016).
In statistically stationary TC flow, the flux of angular momentum 
$J^\omega =r^3(\langle u_r \omega \rangle _{A,t}-\nu \partial _r \langle \omega \rangle _{A,t})$ is exactly conserved (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007); here, 
$u_r$ is the radial velocity, 
$\omega$ the angular velocity of the fluid, 
$r$ is the radial coordinate and the symbol 
$\langle \cdot \rangle _{A,t}$ denotes a time average on a cylindrical surface coaxial with the cylinder axis. The transported quantity 
$J^\omega$ is independent of 
$r$; any flux going through an imaginary cylinder of radius 
$r$ also goes through any other imaginary cylinder, or mathematically 
$\textrm {d} J^\omega /\textrm {d} r = 0$. The response of the system can then be characterized by normalizing 
$J^\omega$ with its value for non-vortical laminar flow 
$J^\omega _{{lam}}=2\nu r_i^2r_o^2(\omega _i-\omega _o)/(r_o^2-r_i^2)$, which gives rise to the pseudo-Nusselt number in TC flow (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007),
\begin{equation} {\textit{Nu}}_\omega=\frac{J^\omega}{J^\omega_{{lam}}}. \end{equation} The key scientific question is to accurately describe the transport throughout the parameter space, i.e. 
${\textit {Nu}}_\omega ={\textit {Nu}}_\omega ({\textit {Ta}},\eta ,a,\varGamma )$. For low Ta, the boundary layers (BLs) remain laminar and as a consequence 
${\textit {Nu}}_\omega$ effectively scales roughly as 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^{1/3}$ (Grossmann et al. Reference Grossmann, Lohse and Sun2016). In the ultimate regime of turbulence, in which both boundary layers and bulk are turbulent, we have 
${Nu}_\omega \propto {\textit {Ta}}^{1/2} / \log {\textit {Ta}}$ (Grossmann & Lohse Reference Grossmann and Lohse2011). The transition to the ultimate regime happens when the boundary layers undergo a shear instability and has been observed at 
${\textit {Ta}}\approx {\textit {Ta}}_c = 3.0 \times 10^8$ for small and medium gaps (
$\eta \geq 0.714$) (Huisman et al. Reference Huisman, van Gils, Grossmann, Sun and Lohse2012; Ostilla-Mónico et al. Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014; Grossmann et al. Reference Grossmann, Lohse and Sun2016). For large gaps, where curvature dominates, the transition is postponed to large values of 
${\textit {Ta}}$, i.e. 
${\textit {Ta}}\approx 10^{10}$ for 
$\eta =0.5$. If one estimates the logarithmic correction, a theoretical estimate for the effective scaling 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^{0.4}$ is obtained at 
$\mathcal{O}({\textit {Ta}})\approx 10^{10}$, which has been confirmed experimentally and numerically. We note that, even if the value of 
${\textit {Ta}}$ for which the transition to the ultimate regime depends on the radius ratio 
$\eta$ and the rotation ratio 
$a$, the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ effective scaling is not affected by these parameters after the transition (van Gils et al. Reference van Gils, Bruggert, Lathrop, Sun and Lohse2011; Paoletti & Lathrop Reference Paoletti and Lathrop2011; Merbold, Brauckmann & Egbers Reference Merbold, Brauckmann and Egbers2013; Ostilla-Mónico et al. Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a).
While the rotation ratio does not affect the effective scaling 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^{0.4}$, it has a strong effect on the proportionality constant. The rotation ratio influences the organization of the flow and increases or decreases the angular momentum transport 
${\textit {Nu}}_\omega$. For a fixed geometry (
$\eta$), and constant driving strength (Ta), a maximum in angular momentum transport (
${\textit {Nu}}_\omega$) can be found for a certain rotation ratio denoted 
$a_{opt}$ (van Gils et al. Reference van Gils, Bruggert, Lathrop, Sun and Lohse2011; Paoletti & Lathrop Reference Paoletti and Lathrop2011; Grossmann et al. Reference Grossmann, Lohse and Sun2016). For the case 
$\eta <0.9$, the maximum has been associated with the strengthening of the large-scale wind (van Gils et al. Reference van Gils, Huisman, Grossmann, Sun and Lohse2012; Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013b; Huisman et al. Reference Huisman, van der Veen, Sun and Lohse2014) and the presence of turbulent intermittent bursts originating from the BLs (van Gils et al. Reference van Gils, Huisman, Grossmann, Sun and Lohse2012). Beyond the point of optimal transport 
$a>a_{{opt}}$, when the counter-rotation is strong, the stabilizing effect of the outer cylinder leads to the detachment of mean vortices from the outer layer which leads to intermittent structures in the radial directions, and decreases the overall angular momentum transport (Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013b).
The value of 
$a_{opt}$ depends on the curvature of the flow, ranging from 
$a_{opt}\approx 0.2$ at 
$\eta =0.5$ (Merbold et al. Reference Merbold, Brauckmann and Egbers2013) to 
$a_{{opt}}\approx 0.4$ at 
$\eta =0.714$ (Huisman et al. Reference Huisman, van der Veen, Sun and Lohse2014). Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b) showed that, as 
$\eta$ increases starting from 
$0.5$, 
$a_{opt}$ becomes larger, corresponding to a system with stronger counter-rotation. However, as 
$\eta$ increases, the peak becomes broader. To disentangle the effect of the rotation ratio 
$a$ from the curvature of the flow, Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) numerically studied the transition from TC flow to rotating plane Couette flow (RPCF), namely the limit 
$\eta \to 1$ in a small aspect ratio domain. In this limit it is more informative to look at the rotation of the cylinders as expressed by 
$R_\varOmega$. When expressed in terms of 
$R_\varOmega$, the asymptotic value (for 
$\eta \to 1$) of 
$R_{\varOmega ,{opt}}$ remains approximately constant. On the other hand, in the limit 
$\eta \to 1$, 
$a(\eta )$ converges to 
$a = 1$ in a rotationless system, and to 
$a = -1$ for all the other cases, showing that for this parameter the transition between TC flow and RPCF is singular (see Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) for further details). Strikingly, Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) find that, for 
$\eta >0.9$ (low curvature), not one maximum of angular momentum transport is present, but two. The first peak, located in the corotating regime, was described as the broad peak. It is associated with strong vortical motions, as evidenced by the radial velocity fluctuations which show a maximum at optimal transport (Brauckmann et al. Reference Brauckmann, Salewski and Eckhardt2016). The second peak, denoted as the narrow peak, was found for counter-rotating cylinders. It appeared only when the driving is sufficiently large, and it was speculated that it supersedes the broad peak for sufficiently large driving.
The appearance of two peaks for small gaps means that several competing mechanisms for the formation of the optimum momentum transport must exist, and that these become blurred for large gaps as the stabilizing effects due to curvature add a third factor. By analysing the 
${\textit {Nu}}_\omega (a)$ relationship using 
$R_\varOmega$ as a control parameter, Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) were able to show that the peaks appearing in the counter-rotating regime at 
$\eta =0.5$ and 
$\eta =0.714$, and the broad peak for corotating cylinders for 
$\eta >0.8$, were basically the same phenomenon, as they both contained strong ordered motions and fell into the same 
$R_\varOmega$ range. As this peak survives the limit of vanishing curvature, it becomes clear that intermittency originated from the stabilizing effect of the outer cylinder does not explain its origin. Instead, Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) divided the TC system into three sub-systems in the spirit of Malkus (Reference Malkus1954): the bulk, and the two boundary layers representing marginally stable TC systems. With this simple model, they were able to predict the location of the broad peak in 
$R_\varOmega$ space, finding good agreement for the prediction at moderate 
${\textit {Ta}}$. Using the same argument, Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) also predicted that the shear in the boundary layers, and hence their transition to turbulence, depends not only on the absolute shear driving, but also on the rotation ratio, which was corroborated by experiments. In this way, they explained the appearance of the narrow peak as an enhancement of angular momentum transport in certain regions of parameter space caused by the ‘early’ transition of the BLs to turbulence. Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) also argued that the narrow peak will dominate the broad peak once the centrifugal instabilities are superseded by shear instabilities, and only one peak would be visible as in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b). This was postulated to happen once the BLs become turbulent for the value of 
$a$ in the broad peak. Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) predicted this to happen at 
${\textit {Ta}}>4.95\times 10^{9}$, close to the transition to the ultimate regime for that 
$\eta$.
In this study we set out to globally and locally probe the angular momentum transport in a wide range of driving strength 
$10^7\leq {\textit {Ta}}\leq 10^{11}$ for the case of low curvature 
$\eta =0.91$, focusing in particular in the 
${\textit {Ta}}$ range 
$10^8\leq {\textit {Ta}}\leq 10^{10}$, where the transition to the ultimate regime happens, and where Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) observed the appearance of two peaks for angular momentum transport using numerical simulations. The main motivation of this study is to elucidate the link between the change in behaviour of the 
${\textit {Nu}}_\omega (\varGamma )$ dependence (Martinez-Arias et al. Reference Martinez-Arias, Peixinho, Crumeyrolle and Mutabazi2014), the vanishing of the broad peak and the changing role of vortical motions (Sacco, Verzicco & Ostilla-Mónico Reference Sacco, Verzicco and Ostilla-Mónico2019) which all happen around this transition. We will use an experimental set-up with very large aspect ratios 
$\varGamma$, which allows for the flow to switch between states, i.e. different roll wavelengths. By doing this, not only can we experimentally confirm the appearance of multiple angular momentum optima in TC flow, which has not been reported yet, but we can also study the transition between regimes dominated by narrow and broad peaks. We will also rule out that they are an effect of artificially constraining the flow to small periodic aspect ratios: switching between two- and three-roll states for varying driving was already reported in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b), and this could have an effect on the two peaks.
Secondly, we will test the predictions of Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) and Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) regarding the mechanism underlying the occurrence of both peaks. Is the broad peak related to vortical motions, which are strengthened by centrifugal forces? Is the narrow peak a consequence of shear? And if so, will it overtake the broad peak and at what turbulence level? By carefully examining the regime where the boundary layers transition, we can better explore the mixed dynamics arising when centrifugal effects and shear are competing side by side and further understand what is happening at the transition to the ultimate regime. To address these questions we conducted both experiments (torque and local velocity measurements) and direct numerical simulations (DNS).
The structure of the paper is as follows. In § 2, we explain the experimental methods. In § 3, we introduce the numerical details of the simulations. In § 4, we experimentally study the global response of the flow throughout a large parameter space of 
${\textit {Ta}}$ and 
$a$. In particular, we reveal transitions and local maxima of the angular momentum transport. In § 5, we complement the experimental findings with numerical simulations and discuss in detail how the size and shape of the Taylor rolls change on varying the rotation parameter 
$R_\varOmega$. The final § 6 contains the conclusions and an outlook for future works.
2. Experimental set-up and measurement procedure
2.1. Set-up
The experiments were carried out in the Twente Turbulent Taylor–Couette facility (
$\text {T}^3\text {C}$) (van Gils et al. Reference van Gils, Bruggert, Lathrop, Sun and Lohse2011). In this apparatus, the ratio 
$\eta$ and aspect ratio 
$\varGamma$ can be adjusted by installing outer cylinders of different dimensions. In this study, the radius of the inner cylinder is 
$r_i={200}\ \textrm {mm}$ and the radius of the outer cylinder (OC) is set to 
$r_o={220}\ \textrm {mm}$. As a consequence, the radius ratio is 
$\eta =r_i/r_o \approx 0.91$ and the aspect ratio results in 
$\varGamma =L/d=46.35$, with 
$d=r_o-r_i={20}\ \textrm {mm}$ and 
$L={927}\ \textrm {mm}$. Two acrylic windows located at the bottom cylinder, which cover the entire gap, allow for the capture of particle image velocimetry (PIV) fields in the 
$r\text {--}\theta$ plane. The advantage of having a second window in the bottom plate is that we can capture two velocity fields for every revolution of the outer cylinder (see figure 1).
Sketch of the experimental set-up. All elements (except for the laser) are mounted on to the frame of the 
$\text {T}^3\text {C}$. This results in no relative motion between the camera and the laser sheet due to mechanical vibrations. The velocity fields obtained with PIV are measured on the 
$r\text {--}\theta$ plane.
2.2. Global measurements: torque
We measure the torque 
$\mathcal {T}$ required to drive the cylinders at constant speed. This is done via a Honeywell model 2404 hollow reaction torque sensor which connects the driving shaft and the inner cylinder. The accuracy of the sensor is 
$0.2$ Nm. From the torque measurements, the Nusselt number can be calculated as follows:
\begin{equation} {\textit{Nu}}_\omega=\left(\frac{r_o^2-r_i^2}{2 \nu r_i^2 r_o^2 {\rm \Delta} \omega} \right) \left(\frac{\mathcal{T}}{2{\rm \pi} \ell_{{eff}} \rho } \right), \end{equation}
where 
$\ell _{{eff}}={536}\ \textrm {mm}$ is the effective length along the cylinder where the torque is measured, the difference of angular frequencies is 
${\rm \Delta} \omega =2{\rm \pi} {\rm \Delta} f=2{\rm \pi} (\,f_i-f_o)$ with 
$f_{i,o}$ the driving frequency of the inner and outer cylinders, respectively, and 
$\rho$ is the fluid density. Typically, the 
$\text {T}^3\text {C}$ facility operates in the ultimate regime of turbulence, where both boundary layers (inner and outer) are turbulent; in our case, where 
$\eta =0.91$, this corresponds to a driving of 
${\textit {Ta}}>\mathcal {O}(10^8)$. Thus, in order to capture the transitional regime (
$\mathcal {O}(10^7)<{\textit {Ta}}<\mathcal {O}(10^8)$), we use working fluids with different values of the kinematic viscosity 
$\nu$. The working fluid – depending on the desired range of Ta to be resolved – is a mixture of water and pure glycerol. The percentage of glycerol in the mixture, along with its corresponding kinematic viscosity and density, can be found in table 1. The liquid temperature is kept constant at 
$21\,^{\circ }\textrm {C}$ during all the experiments.
Properties of the different mixtures used in the experiments. The percentage of glycerol is based on volume. Both the density and kinematic viscosity ratios are calculated with respect to the density 
$\rho _w$ and the kinematic viscosity 
$\nu _w$ of water at 
$21\,^\circ \textrm {C}$. Data taken from Cheng (Reference Cheng2008).
We probe the phase space of 
${Nu}_\omega$ in two different ways. The first one is what we call an 
$a$-sweep, where the angular velocity difference 
${\rm \Delta} \omega$ and thus the driving strength 
${\textit {Ta}}$ is kept constant and the angular velocity ratio 
$a=-\omega _o/\omega _i$ is varied. In this way, we can measure different states in the co- and counter-rotation regimes while the driving (
${\textit {Ta}}$) is fixed. The second type of experiments is the opposite i.e. a 
${\textit {Ta}}$-sweep, where 
$a$ is fixed and 
${\rm \Delta} \omega$, and thus the driving strength 
${\textit {Ta}}$, is increased.
2.3. Local measurements: PIV
We seed the flow with polyamide fluorescent particles with diameters of up to 
${\approx } 20\ \mathrm {\mu }\textrm {m}$ with a seeding density of 
${\approx } 0.01 \ \text {particles} / \text {pixel}$. The emission peak of these particles is centred at 
${\approx }{565}\ \textrm {nm}$. We image the particles in the flow with an Imager SCMOS (
$2560\times 2160\ \textrm{pixel})$ 16 bit camera using a Carl Zeiss Milvus 2.0/100 objective. The illumination of the particles is provided by a Quantel Evergreen 145, 532 nm dual cavity pulsed laser. A cylindrical lens is positioned at the laser output to create a thin light sheet of 
${\approx } {1}\ \textrm {mm}$ thickness. A set of mirrors and a traverse system are installed which allow the laser sheet to move with the frame of the T
$^3$C (see figure 1). Explicitly, the laser beam (the laser is not mounted onto the frame) hits mirror 1 (see figure 1 for the labelling of the mirrors) which is tilted 
${45}^{\circ }$. Light will then be redirected upwards towards to mirror 2 (also tilted at 
${45}^{\circ }$) which redirects it finally towards the OC, perpendicular to both cylinders. A third mirror (mirror 3) is attached to the traverse system of the 
$\text {T}^3\text {C}$ which can move freely in the axial direction. All elements except for the laser head, are mounted on to the frame of the 
$\text {T}^3\text {C}$. This results in no relative motion between the camera and the laser sheet due to mechanical vibrations while the system is rotating.
The experiments require the OC to move freely; thus, a special trigger for the camera is used for the acquisition of the images. This triggering is done by magnets located on top of the OC and a Hall switch mounted onto the frame of the T
$^3$C which outputs a voltage signal every time the magnets pass by. Using this signal as a trigger, we are able to capture two fields (each one corresponding to one window in the bottom plate) per revolution of the OC. The camera is operated in double frame mode with a framerate 
$f$ that depends on the rotation rates of the outer cylinder 
$f_o$. In all cases, however, 
${\rm \Delta} t \leq 1/f$, where 
${\rm \Delta} t$ is the interframe time. In order to increase the contrast between the emission of the light from the particles and the background, we use an Edmund High-Performance Longpass filter 
$550\ \textrm {nm}$ in front of the camera lens.
A total of 7 different flow states have been investigated using particle image velocimetry. These 7 flow states have different 
$a$ and 
${\textit {Ta}}$, as reported in table 2 and – as will be shown later – correspond to the local maxima of the angular momentum transport as function of 
$a$ for a variety of 
${\textit {Ta}}$. In total, 10 different heights were explored for each state, and 500 fields were recorded for each height. The heights are uniformly spaced with a separation of 
$\delta _z = {10}\ \textrm {mm}$, and span the length 
${\rm \Delta} _z = {100}\ \textrm {mm}$ along the axial direction. When normalized with the height of the cylinders 
$L$, this corresponds to a an axial span of 
$\Delta _z/L = (10\delta _z)/L = 0.108$ within the range 
$z/L\in [0.403,0.5]$. for all the experiments. The movement of the laser sheet in the axial direction results in defocusing of the images, therefore the focus is adjusted accordingly for all the explored heights. Accordingly, every velocity field is independently calibrated depending on which window of the bottom plate of the OC is used to acquire it, the height of the laser sheet, the rotation ratio and the driving. Therefore, for a fixed window, 
$a$ and 
${\textit {Ta}}$, the resolution of the PIV fields depends on how far the field is seen by the camera, i.e. the height along the cylinder (see figure 1). In the axial range explored in this study, the resolution of the cameras is within 
${\approx }[30,35] \ \mathrm {\mu } \text {m} / \text {px}$, with the lower/upper bound corresponding to the closest/furthest height from the camera, respectively.
Experimental and numerical flow parameters used in this study. The first two columns show the driving, expressed as either 
${\textit {Ta}}$ or 
$Re_S$. The third and fourth columns show the rotation parameters expressed as either 
$a$ or 
$R_\varOmega$. The last column shows the aspect ratio 
$\varGamma$.
The velocity fields are measured in the 
$r\text {--}\theta$ plane and are computed by a multi-pass algorithm using commercial software (
$\textrm{Davis 8.0}$). The first pass is set to 
$64\times 64 \ \textrm {pixels}$ and the last one is set to 
$24\times 24 \ \textrm {pixels}$ with 
$50\,\%$ overlap of the windows. The fields obtained are then expressed in cylindrical coordinates of the form 
${\boldsymbol {u}}=u_r(r,\theta ,t){\boldsymbol {e}}_{\boldsymbol {r}}+u_\theta (r,\theta ,t) \boldsymbol {e}_{\boldsymbol {\theta }}$, where 
$u_r$ and 
$u_\theta$ are the radial and azimuthal velocities, respectively, which depend on the radius 
$r$, the angular (streamwise) direction 
$\theta$ and time 
$t$. Here, 
$\boldsymbol {e}_{\boldsymbol {r}}$ and 
$\boldsymbol {e}_{\boldsymbol {\theta }}$ are the unit vectors in the radial and azimuthal directions, respectively.
3. Set-up of the direct numerical simulations
In addition to experiments, we perform DNS using an energy-conserving second-order centred finite-difference code for the spatial discretization, while a fractional time-step advancement is adopted in combination with a low-storage third-order Runge–Kutta method. The complete description of the algorithm can be found in Verzicco & Orlandi (Reference Verzicco and Orlandi1996) and van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). This code has been extensively used and validated for TC flows (Ostilla-Mónico et al. Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a).
As mentioned in the introduction, we perform the simulations in a convective reference frame (Dubrulle et al. Reference Dubrulle, Dauchot, Daviaud, Longaretti, Richard and Zahn2005), determined by the parameters 
$Re_s$ and 
$R_\varOmega$ defined in (1.3) and (1.5). According to this scaling the non-dimensional incompressible Navier–Stokes equations read
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
\begin{gather} \displaystyle\frac{\partial \boldsymbol{{u}}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + R_\varOmega \boldsymbol{e}_z\times\boldsymbol{u} ={-} \boldsymbol{\nabla} p + Re_S^{{-}1} {\nabla}^2 \boldsymbol{u}. \end{gather} We chose the same radius ratio 
$\eta =0.91$ as in experiments, which is also the same as in the numerical simulations of Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a,Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohseb). We perform two sets of simulations with fixed Reynolds numbers, 
${\textit {Re}}_S =2.25 \times 10^4$ and 
${\textit {Re}}_S =3.4 \times 10^4$ (or 
${\textit {Ta}}=5.10\times 10^8$ and 
${\textit {Ta}}=1.17 \times 10^9$) while varying 
$R_\varOmega$ (or equivalently 
$a$). Axially periodic boundary conditions are taken with a periodicity length which is similar to the height of the cylinder 
$L$, because even if the boundary conditions are different, the resulting rolls end up being approximately the same size. Non-dimensionally this is expressed by the aspect ratio 
$\varGamma$. In the azimuthal direction, the system is naturally periodic; however, an imposed artificial rotational symmetry of order 
$n_{sym}$ is chosen in order to reduce the computational costs.
We take two computational box sizes. A small box similar to the one used by Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2013a) with 
$\varGamma =2.33$ and 
$n_{sym}=20$. This small box is used for both values of 
${\textit {Re}}_S$, and is large enough to not affect the first-order statistics of the flow (Ostilla-Mónico, Verzicco & Lohse Reference Ostilla-Mónico, Verzicco and Lohse2015). For the case of 
${\textit {Re}}_S =2.25 \times 10^4$, we also run a medium-sized box with an aspect ratio of 
$\varGamma =12.56$, and a rotational symmetry of 
$n_{sym}=3$. This allows the flow some freedom to switch between different roll states, as in Ostilla-Mónico, Lohse & Verzicco (Reference Ostilla-Mónico, Lohse and Verzicco2016). A uniform discretization is used in the azimuthal and axial directions, while a Chebyshev-type clustering near the cylinders is used in the radial direction. The spatial resolution for the small boxes at 
${\textit {Re}}_S=2.25\times 10^4$ and 
${\textit {Re}}_S = 3.4\times 10^4$ was chosen as 
$n_\theta \times n_r\times n_z = 384 \times 512 \times 768$ in the azimuthal, radial and axial directions, which in wall units for the more restrictive case of 
${\textit {Re}}_S=3.4 \times 10^4$ is a resolution of 
${\rm \Delta} z^+\approx 5$, 
${\rm \Delta} x^+ = r{\rm \Delta} \theta ^+\approx 9$ and 
$0.5 \leq {\rm \Delta} r^+\leq 5$. For the medium-size box at 
${\textit {Re}}_S =2.25 \times 10^4$, a grid of 
$n_\theta \times n_r\times n_z = 1728 \times 384 \times 1728$ was chosen, which yields a resolution of 
${\rm \Delta} z^+\approx 5$, 
${\rm \Delta} x^+ = r{\rm \Delta} \theta ^+\approx 9$ and 
$0.4 \leq {\rm \Delta} r^+\leq 2.5$. In order to achieve temporal convergence, the simulations are run until the difference between the time-averaged torque of the inner and the outer cylinders is less than 1 %. The torque is then taken as the average between these two values. The simulations are then run for at least 40 large eddy turnover times 
$tU/d$.
4. Transitions and local maxima in 
${\textit {Nu}}_\omega ({\textit{Ta}},a)$
4.1. Transitions in the ${\textit {Nu}}_\omega ({\textit {Ta}})$ scaling
Firstly, we analyse the scaling laws of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ curve for pure inner cylinder rotation in figure 2. The Nusselt number is compensated by the scaling of the classical regime, i.e. 
${\textit {Nu}}_\omega {\textit {Ta}}^{-1/3}$ and plotted as a function of the driving strength 
${\textit {Ta}}$ for pure inner cylinder rotation only (
$a=0$). We also include the DNS of Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a), and observe a good agreement between the numerics and the experiments. For values of the driving 
${\textit {Ta}}<10^7$, the flow is still in the classical regime, where both BLs are still laminar, and an effective scaling of 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^{0.3}$ can be observed. When the driving strength is increased beyond 
${\textit {Ta}}=\mathcal {O}(10^7)$, the flow enters a transitional regime, with an effective scaling exponent 
$\alpha$ in 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^\alpha$ of 
$\alpha \approx 0.2$. If the driving is further increased, a minimum value of the compensated Nusselt number is reached at a critical Taylor number 
${\textit {Ta}}_c \approx 3.0 \times 10^8$, after which a clear change in the scaling exponent to 
$\alpha =0.4$ can be seen. This indicates the onset of the ultimate regime, which coincides for experiments and numerics.
Compensated Nusselt number as a function of the driving strength 
${\textit {Ta}}$ for the case of pure inner cylinder rotation 
$a=0$ at 
$\eta =0.91$. The grey data points (
${\textit {Ta}}<10^8$) correspond to DNS from Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b). The grey data points for 
${\textit {Ta}}>10^8$ are also DNS simulations but from a different study (Ostilla-Mónico et al. Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a). In addition, each coloured marker at fixed 
${\textit {Ta}}$ corresponds to the driving variation as shown in the legend of figure 3(c). The open circle in light blue corresponds to the DNS data of the current study for 
$\varGamma =2.33$. The transition to the ultimate regime is observed at 
${\textit {Ta}} = {\textit {Ta}}_\text {c} \approx 3 \times 10^8$ (vertical dashed line). The green shaded area corresponds to the region where a local change of slope can be seen due to the disappearance of quiescent wind shearing regions. The black solid lines serve as a reference to indicate the corresponding scaling.
Figure 2 reveals also a second phenomenon, which was not previously reported in experiments. In Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b), the local scaling law was found to be 
${\textit {Nu}}_\omega \sim {\textit {Ta}}^{0.4}$ for 
${\textit {Ta}}>10^{10}$. Indeed, provided 
${\textit {Ta}} > {\textit {Ta}}_c$, the local effective scaling law appears to be the same, with one caveat: there appears to be a local change of slope in the curve around 
${\textit {Ta}}\approx 10^{10}$, where the local scaling exponent increases for a small range in 
${\textit {Ta}}$. The region where this occurs is highlighted in green in figure 2. Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a) observed a similar sudden increase in the local scaling exponent in DNS simulations at 
$a=0$, and attributed it to the sudden disappearance of quiescent wind shearing regions in the boundary layer. After this increase, the dependence of 
${\textit {Nu}}_\omega$ on the roll wavelength was completely lost (Ostilla-Mónico et al. Reference Ostilla-Mónico, Verzicco and Lohse2015). Here, we find evidence that experiments see a similar, sharp increase as observed in figure 2. This will be further investigated in § 5.1, where we will explore how this behaviour is seen across the 
$a$-range, and its effect on the local maxima of angular momentum transport.
4.2. Appearance and shifting of the local maxima
Once we allow the outer cylinder to rotate, we have a more complicated three-dimensional (3-D) parameter space. In figure 3(a), we show the Nusselt number as a function of the rotation ratio 
$a$ for different 
${\textit {Ta}}$. This figure reveals – just as the DNS from Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b) and Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) – that a very pronounced maximum of angular momentum can be found in the corotating regime when the driving is 
${\textit {Ta}}<1.33\times 10^8$. As the driving exceeds the critical value 
${\textit {Ta}}_c$, for approximately a decade of 
${\textit {Ta}}$ we can temporally identify two local angular momentum maxima: the first is located in the corotating regime at 
$a\approx -0.27$ and the second in the counter-rotating regime at 
$a\approx 0.46$. These turbulent states correspond to 
$R_\varOmega =0.16$ and 
$R_\varOmega =0.03$, which are similar to the values found by Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) for 
$\eta =0.91$. The two maxima are prominent for a very small 
${\textit {Ta}}$ range: 
$3.29\times 10^8<{\textit {Ta}}<1.17\times 10^9$.
(a) Nusselt number 
${\textit {Nu}}_\omega$ as a function of rotation ratio 
$a$ for different values of the driving 
${\textit {Ta}}$. The vertical dashed line (
$a=0$) in (a) separates the co- and counter-rotating regimes. (b) Compensated Nusselt number as a function of 
$a$ for four selected 
${\textit {Ta}}$. Here, the solid red circles represent local maxima of angular momentum, where we perform PIV measurements as described in § 4.3. (c) A 3-D representation of the compensated Nusselt number as a function of 
${\textit {Ta}}$ and 
$a$. The black solid lines represent the experiments performed for fixed 
$a$, i.e. 
${\textit {Ta}}$-sweeps. An animated version of this figure can be found in the Supplementary Material. In all figures, the solid lines represent experiments while the dashed lines represent numerics. The colours represent the variation in 
${\textit {Ta}}$ as illustrated by the legend.
This measurement reveals that the two local angular velocity transport maxima for the same driving are not an artefact of the initial conditions or the finite extent of the domain of the numerics. As we further increase the driving beyond 
${\textit {Ta}}>3.31\times 10^9$, the maximum in the corotating regime (broad peak) vanishes, while the maximum in the counter-rotating regime (narrow peak) increases its magnitude. For the largest driving we explore (
${\textit {Ta}}=4.30\times 10^{10}$), only one peak can be detected in the counter-rotating regime, although it is now less sharp. In order to highlight this trend, we show in figure 3(b) the compensated Nusselt number for four selected values of 
${\textit {Ta}}$. Again, note how the value of the driving dictates the occurrence of the maximum of angular momentum transport: if 
${\textit {Ta}}$ is too small, only one peak can be found in the corotating regime. Conversely, if 
${\textit {Ta}}$ is too large, only one peak can be observed albeit for counter-rotation. There is, however, a range of 
${\textit {Ta}}$ which lies in between these two extremes, for which two maxima can be detected. In figure 3(c), we show a 3-D representation of the compensated Nusselt number as a function of 
$a$ and 
${\textit {Ta}}$. In this figure, we included the experiments for fixed 
$a$ (shown in black), i.e. 
${\textit {Ta}}$-sweeps. Note how these experiments agree remarkably well with both the 
$a$-sweeps (shown in colour) and the numerics, mutually validating each other. An animated version of this figure can be found in the supplementary material available at https://doi.org/10.1017/jfm.2020.498. We finally note that, for 
${\textit {Ta}}=2.23\times 10^9$ and 
${\textit {Ta}}=3.31\times 10^9$, discrete jumps in the 
${\textit {Nu}}_\omega (a)$ can be observed for 
$a<0$. This observation can better be seen in figure 3(a) and will be revisited in § 5.2 with the results from the numerical simulations. We note, however, that, although these jumps are similar in magnitude to noise in other curves at different 
${\textit {Ta}}$, we are confident that they are physical and not an artefact of the measurement system. This is based on the high accuracy of the torque sensor in these cases relative to the absolute value of the measured torque in Nm.
In figure 4(a), we show the location of the observed local maxima throughout the parameter space (
$a$, 
${\textit {Ta}}$). Here, we also include the DNS data of Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b) for the same radius ratio 
$\eta =0.91$, albeit for much lower values of Ta. We note that as the driving increases from 
${\textit {Ta}}=\mathcal {O}(10^4)$ towards the critical Taylor number 
${\textit {Ta}}_c$, the peak for corotation moves around, at times towards 
$a=0$, at others away from it. Past the transition, the location of this peak remains relatively stable at 
$a\approx -0.2$ until it vanishes. Regarding the peak for counter-rotation, we see that it only appears when 
${\textit {Ta}}>10^8$ and moves as the driving increases. When 
$1.33\times 10^8 \leq {\textit {Ta}} \leq 1.17\times 10^9$, the peak moves towards higher 
$a$ values of counter-rotation. However, for 
${\textit {Ta}}>\mathcal {O}(10^{10})$, when only one peak is detectable, it seems to move back towards 
$a=0$. This side effect of the disappearance of the broad peak means that the explanation by Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) can be extended. We note, however, that at this driving, 
${\textit {Nu}}_\omega {\textit {Ta}} ^{-1/3}$ becomes less 
$a$-dependent which could over- or underestimate the precise location of the maximum. However, the shifting of the narrow peak is consistent with the asymptotic value of 
$a_{opt}$ at large 
${\textit {Ta}}$ from Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b) and happens around the same 
${\textit {Ta}}$ for which the local effective exponent changes appeared in the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ relation. The reasons for this behaviour will be revisited in § 5.1.
(a) Location in the phase space of the local maxima of angular momentum transport. The blue open circles for 
$\textit {Ta} \leq 2.5\times 10^7$ are the DNS of Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b). The blue open circles located at 
${\textit {Ta}}=1.17\times 10^9$ are from DNS of the current study. The solid circles represent experimental data for the local maxima of angular momentum. The solid orange circles represent turbulent states where we perform PIV experiments as described in § 4.3 and shown in figure 3. (b) Ratios of the magnitude of the angular momentum transport peaks as defined in (4.1). The coloured points represent the experimental data. The open circle represents the DNS simulation shown in (a) for 
${\textit {Ta}}\approx 10^9$. The blue shaded areas represent turbulence levels wherein only one peak in angular momentum can be observed. The transition to the ultimate regime is observed at 
${\textit {Ta}} = {\textit {Ta}}_{c} \approx 3 \times 10^8$. The vertical red dashed line represents the prediction of Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) for the disappearance of the (broad) peak found in corotation, namely at 
${\textit {Ta}}>4.95\times 10^9$.
Interestingly, the value of the driving (
${\textit {Ta}}=5.1\times 10^8$) for which we detect two local maxima is close to the expected value of the transition to the ultimate regime, i.e. 
${\textit {Ta}}_c=3\times 10^8$, and that at this 
${\textit {Ta}}$, the relative magnitude of the peaks is very similar. In order to quantify this observation, we define the ratio of the magnitude of the peaks as
\begin{equation} \zeta\equiv \dfrac{{\textit{Nu}}_\omega(a=a_{{counter}})}{{\textit{Nu}}_\omega(a=a_{{co}})}, \end{equation}
where 
$a_{co}$ and 
$a_{counter}$ denote the 
$a$-s that correspond to the peak for co- and counter-rotation, respectively. In figure 4(b), we report 
$\zeta$ as a function of 
${\textit {Ta}}$, showing that for 
${\textit {Ta}} < {\textit {Ta}}_c$ we have 
$\zeta <1$, whereas for 
${\textit {Ta}} > {\textit {Ta}}_c$ it holds 
$\zeta >1$. This yields an alternative representation of what was originally shown in figure 3: with sufficient driving the peak for counter-rotation will surpass the peak for corotation. As mentioned previously, the magnitude of both peaks seems to be nearly the same (
$\zeta \approx 1$) close to the transition to the ultimate regime. This indicates the link between the appearance of the second peak, and the transitions of the boundary layer as postulated by Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017). Furthermore, around the same values of 
${\textit {Ta}}$, the dependence of 
${\textit {Nu}}_\omega$ on the roll wavelength, and thus on 
$\varGamma$ changes. This was seen as a crossing of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ curves for different values of 
$\varGamma$ around the transition to the ultimate regime (Martinez-Arias et al. Reference Martinez-Arias, Peixinho, Crumeyrolle and Mutabazi2014; Ostilla-Mónico et al. Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a). That these phenomena occur all at the same time indicates the complex character of the transition to the ultimate regime.
Finally, we note that for sufficiently large driving (
${\textit {Ta}}>4.95\times 10^{9}$), the narrow peak completely dominates and the broad peak cannot be detected, as was also postulated by Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017). At these values of 
${\textit {Ta}}$, Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a) showed that the torque would no longer depend on roll wavelength, indicating again that the changes in the peak behaviour are intimately linked to changes in the 
${\textit {Nu}}_\omega (\varGamma )$ relationships.
4.3. Local flow structure and its relation to the local ${Nu}_\omega$ maxima
In the previous section, we showed that the narrow peak (counter-rotation) will surpass the broad peak (corotation) for values of driving 
${\textit {Ta}}>4.95\times 10^{9}$. To further elucidate the mechanisms behind this phenomenon, we investigate the flow locally with PIV measurements. We explore a range of 
${\textit {Ta}}$ which spans values before, close to and beyond the transition to the ultimate regime. For every driving, we investigate two flow states, namely where both the narrow and broad peaks are located as is shown in figure 3(b) and table 2.
We first investigate the strength of the radial flow by looking at the time azimuthally averaged radial velocity, i.e. 
$\langle u_r \rangle _{t,\theta }$. We remind the reader that every velocity field is measured in the 
$r\text {--}\theta$ plane, i.e. 
$u_r=u_r(r,\theta ,t)$. After an average over time and streamwise direction is performed, 
$u_r$ is only a function of the radius. We repeat this operation for all the heights explored which yields finally 
$u_r=u_r(r,z)$. In figure 5, we show 
$u_r(r,z)$ for both peaks, which are located in the corotating and counter-rotating regime, respectively, and as function of 
${\textit {Ta}}$. Here, we can clearly identify regions of negative and positive radial velocity along the axial direction which indicates the presence of Taylor rolls. Strikingly, we find that for all Ta explored, 
$|u_r|$ is much larger for corotation than for counter-rotation. This confirms what was shown numerically by Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016): the broad peak (located for corotation) is accompanied by strong and coherent rolls.
Azimuthally and time-averaged normalized radial velocity obtained from PIV experiments as described in figure 3. The legend on top of each figure represents the value of (
$a,R_\varOmega$). The upper row represents measurements of the peak in the corotating regime while the bottom row shows measurements of the peak for counter-rotation. Along a single column, the Ta is fixed for both co- and counter-rotating states. The dashed lines are added to emphasize the difference between the Ta values.
In order to give a more quantitative picture of the strength of the rolls as a function of the driving, we look at the root-mean-squared value of the radial velocity
\begin{equation} \text{RMS}(\tilde{u_r})\equiv \sqrt{\langle (\langle u_r /u_S, \rangle_{t,\theta,r_{{bulk}}})^2 \rangle_{z_\lambda}}. \end{equation}
Here, 
$r_{{bulk}}$ denotes an average over the radial domain that defines the bulk region of the flow, namely 
$0.3 \leq (r-r_i)/d \leq 0.7$. Due to the presence of the vortical structures in the flow, we average over a distance 
$z_\lambda$ that corresponds to a roll wavelength. In order to find 
$z_\lambda$, we plot 
$\langle u_r \rangle _{t,\theta ,r_{{bulk}}}$ as a function of 
$z$ (not shown). We define 
$z_\lambda$ as the distance that separates two adjacent maximum values of 
$\langle u_r \rangle _{t,\theta ,r_{{bulk}}}$. In cases where there are no visible structures present – for example 
$a=0.46$, 
$a=0.49$, 
$a=0.125$ – we simply take the average over 
$\Delta _z$. For the case of 
$a=-0.21$, we also average over 
$\Delta _z$. The reason is that the maximum points of 
$\langle u_r \rangle _{t,\theta ,r_{{bulk}}}$ in this case, are located at the minimum and maximum heights explored, respectively.
In figure 6(a), we show 
$\text {RMS}(\tilde {u_r})$ as a function of the driving, where we observe that the 
$\text {RMS}(\tilde {u_r})$ of the co-rotating peak decreases with driving. The counter-rotating peak, however, has a non-monotonic behaviour as a function of 
${\textit {Ta}}$.
(a) RMS averaged along the axial direction of the azimuthally time-bulk-averaged radial velocity, normalized with the shear velocity 
$u_S$, as a function of the driving strength 
${\textit {Ta}}$. (b) Wind Reynolds number (based on the radial velocity) as a function of 
${\textit {Ta}}$. In both figures, the blue stars represent states that correspond to the peak in the corotating regime, while the red open circles represent measurements for the peak for counter-rotation.
In addition to the strength of the rolls, we look now at the so-called ‘wind’ by looking at the magnitude of the radial velocity fluctuations defined by
\begin{equation} \sigma_{{bulk}}(u_r)\equiv \langle \sigma_{t,\theta}(u_r) \rangle_{r_{{bulk}},z_\lambda}, \end{equation}
where 
$\sigma _{t,\theta }(u_r)$ is the standard deviation profile of the azimuthal velocity and depends (for a fixed height) only on 
$r$. Here, the brackets 
$\langle \cdot \rangle$ denote the same average as the one performed to 
$\text {RMS}(\tilde {u_r})$. From this characteristic velocity we construct the wind Reynolds number, i.e. 
${\textit {Re}}_w=(d \sigma _{{bulk}}(u_r))/ \nu$. This is called the ‘wind’ Reynolds number because, in analogy to heat transport in Rayleigh–Bénard flow, it is a measure of the energy of the flow that transports the conserved quantity (Grossmann & Lohse Reference Grossmann and Lohse2011). So for the case of TC flow the strength of the wind is given by the (standard deviation of the) radial or axial velocity (Huisman et al. Reference Huisman, van Gils, Grossmann, Sun and Lohse2012).
In the classical regime of turbulence, the unifying theory of Grossmann & Lohse (Reference Grossmann and Lohse2011) predicts a scaling of the Reynolds wind 
${\textit {Re}}_w\propto {\textit {Ta}}^{3/7}$. When the driving is increased, towards the ultimate regime, the logarithmic corrections remarkably cancel out and an effective scaling of 
${\textit {Re}}_w\propto {\textit {Ta}}^{1/2}$ is observed (Grossmann & Lohse Reference Grossmann and Lohse2011; Huisman et al. Reference Huisman, van Gils, Grossmann, Sun and Lohse2012). In figure 6(b), we plot the compensated Reynolds wind with the scaling of the ultimate regime 
${\textit {Re}}_w {\textit {Ta}} ^{-1/2}$ as a function of Ta. Here, we see that, indeed, beyond the transition to the ultimate regime, 
${\textit {Re}}_w$ slowly asymptotes to a 
${\textit {Ta}}^{-1/2}$ scaling for both peaks, which is consistent with the observation of Huisman et al. (Reference Huisman, van Gils, Grossmann, Sun and Lohse2012) for 
$\eta = 0.716$ at 
$a=0$.
We now draw the attention to the case of 
${\textit {Ta}}=5.10 \times 10^8$, which is slightly above 
${\textit {Ta}}_c$, and where both angular momentum peaks have approximately the same magnitude (see figure 4b). As shown in figure 5, for corotation at 
$a=-0.21$ (second image of first row in figure 5), a peculiar pattern of the radial velocity can be appreciated along the axial direction. Due to the shape and flow direction of rolls, we expect to see a succession of positive and negative radial velocities, as is seen in the first and third images. Instead, in the second image (
$a=-0.21$), we encounter two consecutive regions of negative velocity, located in between 
$0.42 \leq z/L \leq 0.48$. This unexplained observation will be revisited in § 5.3, as we unveil more data from the numerical simulations.
5. Boundary Layer Transitions and State Switching
In the previous sections, we highlighted various observations which cannot be fully explained by the limited information we can retain from the experiments. In this section, we therefore turn to direct numerical simulations instead, to provide a more quantitative description of the observations. Namely, the mechanism responsible for the disappearance of the broad peak with sufficient shear (see figure 4a), the observation of discrete jumps in the corotating regime (
$a<0$) for 
${\textit {Ta}}=2.23\times 10^9$ and 
${\textit {Ta}}=3.31\times 10^9$ (see figure 3a) and the ‘peculiar’ pattern of radial velocities which was presented in the previous section (see figure 5). We closely examine the velocities, especially in the boundary layers, and inspect the changes in roll states.
5.1. Disappearance of the broad peak
We first focus on the changing shape of the 
${\textit {Nu}}_\omega (a)$ relationship at 
${\textit {Ta}}\approx 10^{10}$ shown in figure 3(a). In § 2.2, we mentioned how, for 
$a=0$, the sharp changes in the local scaling exponent of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ relationship appeared at 
${\textit {Ta}}\approx 10^{10}$ (see figure 2). These coincided with the disappearance of the torque on the roll wavelength, and cannot be associated with the transition to the ultimate regime, as it happened at a much higher value of 
${\textit {Ta}}$. In figure 7, we show the behaviour of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ curves for three selected values of 
$a$. We see how at 
${\textit {Ta}}\approx 10^{10}$, a reordering of the curves occurs, distinct from the one seen at 
${\textit {Ta}}={\textit {Ta}}_c$. For 
$a=0$ and 
$a=0.25$, the changes in the local scaling exponents increase the effective exponent, while for 
$a=0.5$ the opposite process seems to occur.
Compensated Nusselt number as a function of 
${\textit {Ta}}$ for three rotation ratios. A reordering of the curves around 
${\textit {Ta}}\approx 10^{10}$ can be seen. The green shaded area corresponds to the region where a significant change in the local scaling exponent can be seen due to the disappearance of quiescent wind shearing regions. The black solid line representing the scaling of 
${\textit {Nu}}_\omega \propto {\textit {Ta}}^{0.4}$ is given as reference.
By relating this to the transition observed in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, van der Poel, Verzicco, Grossmann and Lohse2014a), where the quiescent regions of the boundary layer disappeared, we can explain why the counter-rotating maximum shifts. In figure 8, we show the structure of the near-wall region for several values of 
$R_\varOmega$ (and 
$a$). As 
$R_\varOmega$ is increased (i.e. from counter-rotation towards pure inner cylinder rotation), the flow structure self-organizes: the near-wall turbulent streaks occur in a stratified manner and the turbulent Taylor roll is stabilized. This can be seen not only from the visual images of figure 8(a), but also from the root-mean-squared values of 
$u^\prime _\theta$ shown in figure 8(b). For 
$a=0$ they show a significant decrease in the regions which do not generate as many streaks, while there is a much smaller variation for 
$a=0.47$. This stratification of regions streaks can be linked to the ‘plume emission’ regions discussed in Ostilla-Mónico et al. (Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014), where the regions with streaks were classified as plume emitting, and the regions with no streaks were classified as plume impacting.
(a) Instantaneous normalized azimuthal velocity at 
$r^+=15$ for 
${\textit {Ta}}=5\times 10^8$ for 
$a=0$ (
$R_\varOmega =0.09$) and 
$a=0.47$ (
$R_\varOmega =0.03$). The DNS corresponds to the case of 
$\varGamma =2.33$. (b) Azimuthal average of the root-mean-squared azimuthal velocity 
$u_\theta ^{\prime }$ as a function of axial position at 
$r^+=15$ for the two values of rotation 
$a$. A clear break in axial homogeneity can be appreciated for 
$a=0$.
The appearance of quiescent regions is not inconsistent with the idea that the creation of the narrow peak is due to shear instabilities that arise from the BLs. For our value of 
$\eta$, both before and after the transition to the ultimate regime, only parts of the boundary layer are active in producing plumes or streaks. During the transition to the ultimate regime, the plume-emitting part of the laminar boundary layer transitions to turbulence, while the quiescent part remains quiescent. For larger 
${\textit {Ta}}$, beyond 
${\textit {Ta}}=10^{10}$, the quiescent regions disappear, and the entire boundary layer becomes active and turbulent. As this happens, the 
$a=0$ curve surpasses the 
$a=0.5$ curve. Because there is no quiescent area to eliminate for 
$a>0.5$, these branches of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ curve do not show significant change in the local scaling exponents and remain at a lower value. We note that this is unlike what occurs at 
$\eta =0.714$ in Ostilla-Mónico et al. (Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014), where the disappearance of quiescent regions and the transition to the ultimate regime happens simultaneously.
5.2. Roll state switches
We now focus on the discrete jumps in the 
${\textit {Nu}}_\omega (a)$ curve at 
${\textit {Ta}} = 2.23\times 10^9$ and 
${\textit {Ta}}=3.31\times 10^9$ for 
$a<0$ from the experiments (see figure 3). We note that in previous simulations, both by Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) and Ostilla-Mónico et al. (Reference Ostilla-Mónico, Huisman, Jannink, van Gils, Verzicco, Grossmann, Sun and Lohse2014b), a small 
$\varGamma$ domain is used, which essentially fixes the roll size. The fact that no discrete jumps are detected in the numerics with small boxes indicates that roll state switching might be responsible for this. Thus, in order to explain these jumps, we need to use both a small computational box which accommodates a single roll pair, as well as medium-sized boxes which can sustain changes in roll number. In this way, we can capture how the rolls manifest as the system goes from pure counter-rotation to corotation, and whether state switching takes place or not. We note that the number of rolls does not significantly affect the value of 
${\textit {Nu}}_\omega$ and that the only dependence comes from the roll wavelength (Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013a; Ostilla-Mónico et al. Reference Ostilla-Mónico, Lohse and Verzicco2016). However, by enforcing a periodicity length, both effects will manifest simultaneously. The roll number 
$N$ and the (non-dimensional) roll wavelength 
$\lambda _z$ are related by 
$\lambda _z=\varGamma /N$.
In figure 9, we show the compensated Nusselt number for both the small and the medium boxes for two 
${\textit {Ta}}$, and the experiments for 
${\textit {Ta}}=5.10\times 10^8$. The numerical simulations at 
${\textit {Ta}}=5.10\times 10^8$ reveal that for 
$a<-0.5$ the value of 
${\textit {Nu}}_\omega$ is at most weakly dependent on 
$\varGamma$. However, for 
$a\geq -0.5$ the two numerical studies and the experimental measurements result in different values of 
${\textit {Nu}}_\omega$. We note that 
${\textit {Nu}}_\omega$ was shown to have the largest dependence on the roll wavelength, with the first- and second-order statistics also being affected in a small manner by Ostilla-Mónico et al. (Reference Ostilla-Mónico, Lohse and Verzicco2016).
Compensated Nusselt number as a function of 
$a$ obtained from the DNS of the current study. The green solid line is the experimental data performed at 
${\textit {Ta}}=5.10\times 10^8$.
To further explore this, in figure 10 we show the azimuthally and temporally averaged radial velocity for the medium-sized domain (
$\varGamma =12.56$), as we vary 
$a$ from the corotating to the counter-rotating regime. In terms of 
$R_\varOmega$ this is equivalent to varying the Coriolis force from anti-cyclonic to zero. Here, we can see that the flow self-organizes in Taylor rolls as the Coriolis force starts to become dominant, i.e. when 
$R_\varOmega \neq 0$ (Sacco et al. Reference Sacco, Verzicco and Ostilla-Mónico2019). As 
$R_\varOmega$ increases, we first observe that the rolls become sharper and more prominent, as evidenced by an increment in 
$|u_r|$, which is also observed in the experiments (see figure 5). As we approach the ‘broad peak’, the number of roll pairs switch, first from four to five, then from five to six, which sharply reduces the roll wavelength. This reduction in wavelength clearly has an effect on the torque as can be seen in figure 9 for 
$a>-0.5$ The change in the roll wavelength changes the proportion between quiescent and active regions inside the boundary layers (van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2011; van der Poel et al. Reference van der Poel, Stevens, Sugiyama and Lohse2012). As 
${\textit {Ta}}$ increases, these sharp changes in the 
${\textit {Nu}}_\omega (a)$ curve disappear, because roll state switching does not modify the fraction of boundary layer regions which emit plumes.
Azimuthally time-averaged radial velocity normalized with the characteristic velocity 
$U$ obtained from DNS for 
${\textit {Ta}}=5.10\times 10^8$ and 
$\varGamma =12.56$. Each plot represents a different rotation state, quantified by either 
$a$ or 
$R_\varOmega$. The legend on top of each figure represents the value of (
$a,R_\varOmega$). The flow state that corresponds to either the ‘broad’ or ‘narrow’ peak are highlighted in red. The arrows indicate the direction of the increment of either 
$a$ or 
$R_\varOmega$.
The influence of the large structures on the torque can also be appreciated when we compare the 
$\varGamma =12.56$ simulation with the experimental data (
$\varGamma =46.35$) shown in figure 9. Here, we see that both the experimental and numerical data coincide within the region 
$0.3<a<1$, where the Coriolis force is small and the rolls only start to become organized. However, in the region between 
$-0.5<a<0.3$ there are significant discrepancies between experiments and numerics. These are probably a consequence of state switching. In the medium size box (
$\varGamma =12.56$), a switch in the number of roll pairs can take place, but the values that the roll wavelength can take are more restricted due to the size domain and thus, leading to large variation in roll wavelength across different states. In the experiment 
$\varGamma$ is much larger, so state switches will generally have a very small effect on the roll wavelength, and as a consequence a small effect on the torque. This can explain both the jumps in the 
${\textit {Nu}}_\omega (a)$ curve, and why they are much smaller than for numerics.
5.3. Transient roll dynamics
During the transition to the ultimate regime, the size distribution of the roll changes (see figure 10). This can be seen from the changing distance between the maximum values of 
$|u_r|$. This change in the roll dynamics suggests that there could be an interval of 
$R_\varOmega$ in which the number of rolls is allowed to change, and it is connected to the pattern depicted at the end of § 4.3, represented in the second image of figure 5. We can see there the presence of two similarly signed radial velocities that are close to each other. If one would approach this from the point of view that Taylor rolls are present, it could only be possible to think that two rolls rotating in the same way are next to each other, as the region in between them is too small to allow the presence of a well formed counter-rotating roll. Otherwise, it could also be the case that this region has no rolls and we are encountering a local or transient phenomenon. In order to explore this event, we have to analyse the instantaneous velocity fields that DNS provides. Experimental results examine consecutive meridional planes (
$r\text {--}\theta$) and assume that the flow is azimuthally homogeneous. However, if the structures change only locally, this could be wiped out by an averaging operation.
As 
$R_\varOmega$ increases (see figure 10), the rolls become sharper, until at a certain point, one of them can start to split up locally, such that transiently, outflow or inflow regions with the same sign coexist very close to each other at certain values of the azimuth. As 
$R_\varOmega$ is further increased, the roll splits up completely to form a new roll pair. However, the speed at which this roll ‘dislocation’ exists and propagates to ‘fracture’ the roll is a priori unknown.
In order to closely inspect the change in the morphology of the rolls as a function of the Coriolis force, we perform DNS by slowly changing the value of 
$R_\varOmega$ over a certain number of large eddy turnover times 
$\tau =tU/d$, where 
$U$ is the characteristic velocity as defined in § 3. The simulations are performed for 
${\textit {Ta}}=5.10\times 10^8$ and for the medium-sized box (
$\varGamma =12.56$). The idea is to capture the evolution of the radial velocity and its effect on 
${\textit {Nu}}_\omega$ as we slowly increase 
$R_\varOmega$ in discrete steps. We initiate the simulation at 
$\tau =0$ with a statistical stationary TC flow and 
$R_\varOmega =0$. Next, at 
$\tau =25$, we increase the value of the Coriolis force to 
$R_\varOmega =0.02$. When 
$\tau =50$ is reached, we set finally 
$R_\varOmega$ to 0.03. In figure 11(a), we show the instantaneous 
${\textit {Nu}}_\omega$ and in figure 11(b), we show the space–time evolution of the instantaneous radial velocity 
$u_r$ at mid-gap. We note that, globally, a correlation between transient regions of 
${\textit {Nu}}_\omega$ and the change in 
$R_\varOmega$ can be observed. This is caused by the increasing acceleration of the flow, in which the number of rolls changes as well. Locally, we clearly see two instants in which two consecutive rolls begin to approach each other until they merge. During this transient phenomenon, that lasts 
${\approx } 10\tau$, the shape of the rolls changes, and at the end also the global number of rolls pairs due to the merging event. During this time two consecutive outflow regions are allowed to coexist close to each other, while the inflow region slowly disappear. This indicates, just as was observed in the experiments shown in figure 5, that close to the transition to the ultimate regime, transient events and changes in the roll dynamics are allowed under suitable conditions.
DNS results of varying 
$R_\varOmega$ for 
${\textit {Ta}}=5.10\times 10^8$ and 
$\varGamma =12.56$. (a) Temporal evolution of the dissipation 
${\textit {Nu}}_\omega$. (b) Time–space diagram of the instantaneous 
$u_r$ at mid-gap. By looking at the regions of negative and positive radial velocities, which represent the outflow or inflow of rolls, two regions of merging rolls can clearly be seen. These mergers are highlighted with the dashed circles. The vertical dashed lines in (a) represent a stage during the simulation for a fixed value of 
$R_\varOmega$.
6. Summary and Conclusions
We probe the angular momentum transport with both experiments and direct numerical simulations for 
$\eta =0.91$ as a function of the driving which we quantify with 
${\textit {Ta}}$ and the magnitude of the Coriolis force which we quantify with either the rotation ratio 
$a$ or the rotation number 
$R_\varOmega$. The range of shear driving we explore spans several decades of 
${\textit {Ta}}$, namely 
$\mathcal {O}(10^7)\text {--}\mathcal {O}(10^{10})$, which includes the transition to the ultimate regime at 
${\textit {Ta}}\approx {\textit {Ta}}_c=3\times 10^8$. In the vicinity of 
${\textit {Ta}}_c$, our experiments reveal the presence of two local angular momentum transport maxima for fixed 
${\textit {Ta}}$, as was firstly reported by the numerical works of Brauckmann et al. (Reference Brauckmann, Salewski and Eckhardt2016) and Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017). Thus we have demonstrated that this occurrence is not an artifact of the axial domain used in the numerics. We showed that the broad peak is associated with the strengthening of the large-scale structures (Taylor rolls) due to an increment in centrifugal forcing 
$R_\varOmega$, as it is evidenced by both PIV and DNS. In addition, the numerical simulations reveal that the appearance of the narrow peak is a consequence of shear instabilities that arise from the BLs. Moreover, we show that as the driving increases, the broad peak remains roughly at the same 
$R_\varOmega$ while decreasing its magnitude. Conversely, the narrow peak increases its magnitude while experiencing a shift in 
$R_\varOmega$. We attribute this to the different way in which the near-wall structures self-organizes at 
${\textit {Ta}}\approx 10^{10}$. Explicitly, the disappearance of quiescent wind sheared regions at the cylinder walls, which suddenly transition to a state where the whole BLs can emit plumes, as evidenced by a sharp localized increase in the effective exponent of the 
${\textit {Nu}}_\omega ({\textit {Ta}})$ relationship. With sufficient driving, both the experiments and numerics confirm that the shear instabilities dominate over centrifugal instabilities and the broad peak disappears at around 
${\textit {Ta}}=4.95 \times 10^9$. Beyond this value, only the narrow peak can be detected.
In the experiments, a peculiar state of the Taylor rolls is observed at 
$a=-0.210$ (
$R_\varOmega =0.141$) and 
${\textit {Ta}}=5.10\times 10^8$, which is close to 
${\textit {Ta}}_c$, where both the narrow and broad peaks are present and have roughly the same magnitude. Here, we observe two neighbouring Taylor rolls which rotate opposite to each other. We explore the possibility that this unusual state is produced by roll splitting, which occurs when the axial extent is sufficiently large and for sufficient 
$R_\varOmega$. The numerics at 
${\textit {Ta}}=5.10\times 10^8$ reveal that, as the magnitude of the Coriolis force increases via 
$R_\varOmega$, the Taylor rolls undergo a reorganization and a change in both their number and wavelength. This effect translates into a change in the global torque, i.e. 
${\textit {Nu}}_\omega$, which is evidenced by different 
${\textit {Nu}}_\omega (a)$ curves for different values of 
$\varGamma$. We explore the reorganization of the rolls in detail by performing a last numerical simulation where we slowly vary 
$R_\varOmega$ several times. Here, we observe that for the same driving 
${\textit {Ta}}$ and curvature 
$\eta$ as in the experiments (different 
$\varGamma$), the merging of two adjacent rolls is possible. This suggests that the ‘peculiar’ pattern of roll formation discussed previously can be a consequence of such long-time transient effects. Further work in this direction is encouraged in order to thoroughly explore this new feature of secondary flows in TC flow.
Acknowledgements
This paper is dedicated to B. Eckhardt, who unexpectedly passed away in July 2019. Over the last 15 years, we had various stimulating discussions with him on TC turbulence in general and over the last few years in particular on the double-peak nature of the angular momentum transport, for which we are really thankful. We acknowledge G.-W. Bruggert for technical support. This work was funded by STW, FOM and MCEC, which are part of the Netherlands Organization for Scientific Research (NWO). C.S. acknowledges financial support the Natural Science Foundation of China under grant nos 11988102, 91852202 and 11672156. We thank the Center for Advanced Computing and Data Science (CACDS) at the University of Houston for providing computing resources and we also acknowledge PRACE for awarding us access through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC), under PRACE project number 2017174146.
Declaration of interests
The authors report no conflict of interest.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2020.498.
Open access
