Double maxima of angular momentum transport in small gap $\eta =0.91$ Taylor–Couette turbulence

Abstract We use experiments and direct numerical simulations to probe the phase space of low-curvature Taylor–Couette flow in the vicinity of the ultimate regime. The cylinder radius ratio is fixed at $\eta =r_i/r_o=0.91$, where $r_i \, (r_o)$ is the inner (outer) cylinder radius. Non-dimensional shear drivings (Taylor numbers ${\textit {Ta}}$) in the range $10^7\leq {\textit {Ta}}\leq 10^{11}$ are explored for both co- and counter-rotating configurations. In the ${\textit {Ta}}$ range $10^8\leq {\textit {Ta}}\leq 10^{10}$, we observe two local maxima of the angular momentum transport as a function of the cylinder rotation ratio, which can be described as either ‘co-’ or ‘counter-rotating’ due to their location or as ‘broad’ or ‘narrow’ due to their shape. We confirm that the broad peak is accompanied by the strengthening of the large-scale structures, and that the narrow peak appears once the driving (Ta) is strong enough. As first evidenced in numerical simulations by Brauckmann et al. (J. Fluid Mech., vol. 790, 2016, pp. 419–452), the broad peak is produced by centrifugal instabilities and that the narrow peak is a consequence of shear instabilities. We describe how the peaks change with ${\textit {Ta}}$ as the flow becomes more turbulent. Close to the transition to the ultimate regime when the boundary layers (BLs) become turbulent, the usual structure of counter-rotating Taylor vortex pairs breaks down and stable unpaired rolls appear locally. We attribute this state to changes in the underlying roll characteristics during the transition to the ultimate regime. Further changes in the flow structure around ${\textit {Ta}}\approx 10^{10}$ cause the broad peak to disappear completely and the narrow peak to move. This second transition is caused when the regions inside the BLs which are locally smooth regions disappear and the whole boundary layer becomes active.


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 1923;Chandrasekhar 1981;Andereck, Liu & Swinney 1986;Lewis & Swinney 1999;van Gils et al. 2011;Paoletti & Lathrop 2011;Fardin, Perge & Taberlet 2014;Ostilla-Mónico et al. 2014a;Grossmann, Lohse & Sun 2016). The basic TC geometry is characterized by two parameters: the first is the radius ratio η = r i /r o , where r i and r o are the inner and outer radii, respectively. The second is the aspect ratio Γ = 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 2007) where ω i,o are the inner and outer angular velocities, respectively, and ν is the kinematic viscosity of the fluid. The second control parameter is the rotation ratio 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 Ta, η and a, one could alternatively describe the parameter space in a convective reference frame as proposed by Dubrulle et al. (2005) such that the cylinders rotate with opposite velocities ±U/2 and the entire system then rotates with angular velocity Ω = Ω rf e z around the central axis. Here, U is the characteristic velocity U = 2(u i − ηu o )/(1 + η), Ω rf = (r i ω i + r o ω o )/(r i + r o ) is the mean angular velocity, e z is the unit vector in the axial direction and u i,o = r i,o ω 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 Re S for the driving strength, the curvature number R C and the rotation number R Ω 3) (1.5) We remark that Re S ∝ √ Ta, and that the rotation number R Ω is connected with the negative rotation ratio a by (1.6) While the choice of one set of parameters might seem arbitrary at first, we note that R Ω 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 2016).
In statistically stationary TC flow, the flux of angular momentum J ω = r 3 ( u r ω A,t − ν∂ r ω A,t ) is exactly conserved (Eckhardt et al. 2007); here, u r is the radial velocity, ω the angular velocity of the fluid, r is the radial coordinate and the symbol · A,t denotes a time average on a cylindrical surface coaxial with the cylinder axis. The transported quantity J ω is independent of r; any flux going through an imaginary cylinder of radius r also goes through any other imaginary cylinder, or mathematically dJ ω /dr = 0. The response of the system can then be characterized by normalizing J ω with its value for non-vortical laminar flow J ω lam = 2νr 2 i r 2 o (ω i − ω o )/(r 2 o − r 2 i ), which gives rise to the pseudo-Nusselt number in TC flow (Eckhardt et al. 2007), (1.7) The key scientific question is to accurately describe the transport throughout the parameter space, i.e. Nu ω = Nu ω (Ta, η, a, Γ ). For low Ta, the boundary layers (BLs) remain laminar and as a consequence Nu ω effectively scales roughly as Nu ω ∝ Ta 1/3 (Grossmann et al. 2016). In the ultimate regime of turbulence, in which both boundary layers and bulk are turbulent, we have Nu ω ∝ Ta 1/2 / log Ta (Grossmann & Lohse 2011). The transition to the ultimate regime happens when the boundary layers undergo a shear instability and has been observed at Ta ≈ Ta c = 3.0 × 10 8 for small and medium gaps (η ≥ 0.714) Ostilla-Mónico et al. 2014;Grossmann et al. 2016). For large gaps, where curvature dominates, the transition is postponed to large values of Ta, i.e. Ta ≈ 10 10 for η = 0.5. If one estimates the logarithmic correction, a theoretical estimate for the effective scaling Nu ω ∝ Ta 0.4 is obtained at O(Ta) ≈ 10 10 , which has been confirmed experimentally and numerically. We note that, even if the value of Ta for which the transition to the ultimate regime depends on the radius ratio η and the rotation ratio a, the Nu ω (Ta) effective scaling is not affected by these parameters after the transition (van Gils et al. 2011;Paoletti & Lathrop 2011;Merbold, Brauckmann & Egbers 2013;Ostilla-Mónico et al. 2014a).
While the rotation ratio does not affect the effective scaling Nu ω ∝ 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 Nu ω . For a fixed geometry (η), and constant driving strength (Ta), a maximum in angular momentum transport (Nu ω ) can be found for a certain rotation ratio denoted a opt (van Gils et al. 2011;Paoletti & Lathrop 2011;Grossmann et al. 2016). For the case η < 0.9, the maximum has been associated with the strengthening of the large-scale wind Brauckmann & Eckhardt 2013b;Huisman et al. 2014) and the presence of turbulent intermittent bursts originating from the BLs . 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 2013b).
The value of a opt depends on the curvature of the flow, ranging from a opt ≈ 0.2 at η = 0.5 (Merbold et al. 2013) to a opt ≈ 0.4 at η = 0.714 (Huisman et al. 2014). Ostilla-Mónico et al. (2014b) showed that, as η increases starting from 0.5, a opt becomes larger, corresponding to a system with stronger counter-rotation. However, as η increases, the peak becomes broader. To disentangle the effect of the rotation ratio a from the curvature of the flow, Brauckmann et al. (2016) numerically studied the transition from TC flow to rotating plane Couette flow (RPCF), namely the limit η → 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 Ω . When expressed in terms of R Ω , the asymptotic value (for η → 1) of R Ω,opt remains approximately constant. On the other hand, in the limit η → 1, a(η) 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. (2016) for further details). Strikingly, Brauckmann et al. (2016) find that, for η > 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. 2016). 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 Nu ω (a) relationship using R Ω as a control parameter, Brauckmann et al. (2016) were able to show that the peaks appearing in the counter-rotating regime at η = 0.5 and η = 0.714, and the broad peak for corotating cylinders for η > 0.8, were basically the same phenomenon, as they both contained strong ordered motions and fell into the same R Ω 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 (2017) divided the TC system into three sub-systems in the spirit of Malkus (1954): 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 Ω space, finding good agreement for the prediction at moderate Ta. Using the same argument, Brauckmann & Eckhardt (2017) 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 (2017) 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. (2014b). This was postulated to happen once the BLs become turbulent for the value of a in the broad peak. Brauckmann & Eckhardt (2017) predicted this to happen at Ta > 4.95 × 10 9 , close to the transition to the ultimate regime for that η.
In this study we set out to globally and locally probe the angular momentum transport in a wide range of driving strength 10 7 ≤ Ta ≤ 10 11 for the case of low curvature η = 0.91, focusing in particular in the Ta range 10 8 ≤ Ta ≤ 10 10 , where the transition to the ultimate regime happens, and where Brauckmann & Eckhardt (2017) 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 Nu ω (Γ ) dependence (Martinez-Arias et al. 2014), the vanishing of the broad peak and the changing role of vortical motions (Sacco, Verzicco & Ostilla-Mónico 2019) which all happen around this transition. We will use an experimental set-up with very large aspect ratios Γ , 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. (2014b), and this could have an effect on the two peaks.
Secondly, we will test the predictions of Brauckmann et al. (2016) and Brauckmann & Eckhardt (2017) 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 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 Ω . The final § 6 contains the conclusions and an outlook for future works.
2. Experimental set-up and measurement procedure

Set-up
The experiments were carried out in the Twente Turbulent Taylor-Couette facility (T 3 C) (van Gils et al. 2011). In this apparatus, the ratio η and aspect ratio Γ can be adjusted by installing outer cylinders of different dimensions. In this study, the radius of the inner cylinder is r i = 200 mm and the radius of the outer cylinder (OC) is set to r o = 220 mm. As a consequence, the radius ratio is η = r i /r o ≈ 0.91 and the aspect ratio results in Γ = L/d = 46.35, with d = r o − r i = 20 mm and L = 927 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-θ 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).

Global measurements: torque
We measure the torque 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: where eff = 536 mm is the effective length along the cylinder where the torque is measured, the difference of angular frequencies is (O(10 7 ) < Ta < O(10 8 )), we use working fluids with different values of the kinematic viscosity ν. The working fluid -depending on the desired range of Ta to be resolvedis 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 • C during all the experiments. We probe the phase space of Nu ω in two different ways. The first one is what we call an a-sweep, where the angular velocity difference Δω and thus the driving strength Ta is kept constant and the angular velocity ratio a = −ω o /ω i is varied. In this way, we can measure different states in the co-and counter-rotation regimes while the driving (Ta) is fixed. The second type of experiments is the opposite i.e. a Ta-sweep, where a is fixed and Δω, and thus the driving strength Ta, is increased.

Local measurements: PIV
We seed the flow with polyamide fluorescent particles with diameters of up to ≈20 μm with a seeding density of ≈0.01 particles/pixel. The emission peak of these particles is centred at ≈565 nm. We image the particles in the flow with an Imager SCMOS (2560 × 2160 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 ≈1 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 • . Light will then be redirected upwards towards to mirror 2 (also tilted at 45 • ) which redirects it finally towards the OC, perpendicular to both cylinders. A third mirror (mirror 3) is attached to the traverse system of the T 3 C which can move freely in the axial direction. All elements except for the laser head, are mounted on to the frame of the T 3 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, Δt ≤ 1/f , where Δ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 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 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 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 δ z = 10 mm, and span the length Δ z = 100 mm along the axial direction. When normalized with the height of the cylinders L, this corresponds to a an axial span of z /L = (10δ z )/L = 0.108 within the range z/L ∈ [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 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 ≈[30, 35] μm/px, with the lower/upper bound corresponding to the closest/furthest height from the camera, respectively.
The velocity fields are measured in the r-θ plane and are computed by a multi-pass algorithm using commercial software (Davis 8.0). The first pass is set to 64 × 64 pixels and the last one is set to 24 × 24 pixels with 50 % overlap of the windows. The fields obtained are then expressed in cylindrical coordinates of the form u = u r (r, θ, t)e r + u θ (r, θ, t)e θ , where u r and u θ are the radial and azimuthal velocities, respectively, which depend on the radius r, the angular (streamwise) direction θ and time t. Here, e r and e θ are the unit vectors in the radial and azimuthal directions, respectively.

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 (1996) and van der Poel et al. (2015). This code has been extensively used and validated for TC flows (Ostilla-Mónico et al. 2014a).
As mentioned in the introduction, we perform the simulations in a convective reference frame (Dubrulle et al. 2005), determined by the parameters Re s and R Ω defined in (1.3) and (1.5). According to this scaling the non-dimensional incompressible Navier-Stokes equations read We chose the same radius ratio η = 0.91 as in experiments, which is also the same as in the numerical simulations of Ostilla-Mónico et al. (2014a,b). We perform two sets of simulations with fixed Reynolds numbers, Re S = 2.25 × 10 4 and Re S = 3.4 × 10 4 (or Ta = 5.10 × 10 8 and Ta = 1.17 × 10 9 ) while varying R Ω (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 Γ . 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 (2013a) with Γ = 2.33 and n sym = 20. This small box is used for both values of Re S , and is large enough to not affect the first-order statistics of the flow (Ostilla-Mónico, Verzicco & Lohse 2015). For the case of Re S = 2.25 × 10 4 , we also run a medium-sized box with an aspect ratio of Γ = 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 (2016). 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 Re S = 2.25 × 10 4 and Re S = 3.4 × 10 4 was chosen as n θ × n r × n z = 384 × 512 × 768 in the azimuthal, radial and axial directions, which in wall units for the more restrictive case of Re S = 3.4 × 10 4 is a resolution of Δz + ≈ 5, Δx + = rΔθ + ≈ 9 and 0.5 ≤ Δr + ≤ 5. For the medium-size box at Re S = 2.25 × 10 4 , a grid of n θ × n r × n z = 1728 × 384 × 1728 was chosen, which yields a resolution of Δz + ≈ 5, Δx + = rΔθ + ≈ 9 and 0.4 ≤ Δr + ≤ 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.

Transitions and local maxima in Nu ω (Ta, a)
4.1. Transitions in the Nu ω (Ta) scaling Firstly, we analyse the scaling laws of the Nu ω (Ta) curve for pure inner cylinder rotation in figure 2. The Nusselt number is compensated by the scaling of the classical regime, i.e. Nu ω Ta −1/3 and plotted as a function of the driving strength Ta for pure inner cylinder rotation only (a = 0). We also include the DNS of Ostilla-Mónico et al. (2014a), and observe a good agreement between the numerics and the experiments. For values of the driving Ta < 10 7 , the flow is still in the classical regime, where both BLs are still laminar, and an effective scaling of Nu ω ∝ Ta 0.3 can be observed. When the driving strength is increased beyond Ta = O(10 7 ), the flow enters a transitional regime, with an effective scaling exponent α in Nu ω ∝ Ta α of α ≈ 0.2. If the driving is further increased, a minimum value of the compensated Nusselt number is reached at a critical Taylor number Ta c ≈ 3.0 × 10 8 , after which a clear change in the scaling exponent to α = 0.4 can be seen. This indicates the onset of the ultimate regime, which coincides for experiments and numerics. Figure 2 reveals also a second phenomenon, which was not previously reported in experiments. In Ostilla-Mónico et al. (2014b), the local scaling law was found to be Nu ω ∼ Ta 0.4 for Ta > 10 10 . Indeed, provided Ta > 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 Ta ≈ 10 10 , where the local scaling exponent increases for a small range in Ta 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 Ta. This figure reveals -just as the DNS from Ostilla-Mónico et al. (2014b) and Brauckmann et al. (2016) -that a very pronounced maximum of angular momentum can be found in the corotating regime when the driving is Ta < 1.33 × 10 8 . As the driving exceeds the critical value Ta c , for approximately a decade of Ta we can temporally identify two local angular momentum maxima: the first is located in the corotating regime at a ≈ −0.27 and the second in the counter-rotating regime at a ≈ 0.46. These turbulent states correspond to R Ω = 0.16 and R Ω = 0.03, which are similar to the values found by Brauckmann et al. (2016) for η = 0.91. The two maxima are prominent for a very small Ta range: 3.29 × 10 8 < Ta < 1.17 × 10 9 .
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 Ta > 3.31 × 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 (Ta = 4.30 × 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  dictates the occurrence of the maximum of angular momentum transport: if Ta is too small, only one peak can be found in the corotating regime. Conversely, if Ta is too large, only one peak can be observed albeit for counter-rotation. There is, however, a range of 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 Ta. In this figure, we included the experiments for fixed a (shown in black), i.e. 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 Ta = 2.23 × 10 9 and Ta = 3.31 × 10 9 , discrete jumps in the Nu ω (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 -0.5 10 4 10 7 10 8 10 9 10 10 10 11 Only one peak different 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, Ta). Here, we also include the DNS data of Ostilla-Mónico et al. (2014b) for the same radius ratio η = 0.91, albeit for much lower values of Ta. We note that as the driving increases from Ta = O(10 4 ) towards the critical Taylor number 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 ≈ −0.2 until it vanishes. Regarding the peak for counter-rotation, we see that it only appears when Ta > 10 8 and moves as the driving increases. When 1.33 × 10 8 ≤ Ta ≤ 1.17 × 10 9 , the peak moves towards higher a values of counter-rotation. However, for Ta > 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 (2017) can be extended. We note, however, that at this driving, Nu ω 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 Ta from Ostilla-Mónico et al. (2014b) and happens around the same Ta for which the local effective exponent changes appeared in the Nu ω (Ta) relation. The reasons for this behaviour will be revisited in § 5.1.
Interestingly, the value of the driving (Ta = 5.1 × 10 8 ) for which we detect two local maxima is close to the expected value of the transition to the ultimate regime, i.e. Ta c = 3 × 10 8 , and that at this 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 ζ ≡ Nu ω (a = a counter ) Nu ω (a = a co ) , (4.1) 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 ζ as a function of Ta, showing that for Ta < Ta c we have ζ < 1, whereas for Ta > Ta c it holds ζ > 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 (ζ ≈ 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 (2017). Furthermore, around the same values of Ta, the dependence of Nu ω on the roll wavelength, and thus on Γ changes. This was seen as a crossing of the Nu ω (Ta) curves for different values of Γ around the transition to the ultimate regime (Martinez-Arias et al. 2014;Ostilla-Mónico et al. 2014a). 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 (Ta > 4.95 × 10 9 ), the narrow peak completely dominates and the broad peak cannot be detected, as was also postulated by Brauckmann & Eckhardt (2017). At these values of Ta, Ostilla-Mónico et al. (2014a) 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 Nu ω (Γ ) relationships.

Local flow structure and its relation to the local Nu ω maxima
In the previous section, we showed that the narrow peak (counter-rotation) will surpass the broad peak (corotation) for values of driving Ta > 4.95 × 10 9 . To further elucidate the mechanisms behind this phenomenon, we investigate the flow locally with PIV measurements. We explore a range of 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. u r t,θ . We remind the reader that every velocity field is measured in the r-θ plane, i.e. u r = u r (r, θ, 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 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. (2016): the broad peak (located for corotation) is accompanied by strong and coherent rolls.
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 Here, r bulk denotes an average over the radial domain that defines the bulk region of the flow, namely 0.3 ≤ (r − r i )/d ≤ 0.7. Due to the presence of the vortical structures in the flow, we average over a distance z λ that corresponds to a roll wavelength. In order to find z λ , we plot u r t,θ,r bulk as a function of z (not shown). We define z λ as the distance that separates two adjacent maximum values of u r t,θ,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 z . For the case of a = −0.21, we also average over z . The reason is that the maximum points of u r t,θ,r bulk in this case, are located at the minimum and maximum heights explored, respectively. In figure 6(a), we show RMS(ũ r ) as a function of the driving, where we observe that the RMS(ũ r ) of the co-rotating peak decreases with driving. The counter-rotating peak, however, has a non-monotonic behaviour as a function of Ta. 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 where σ t,θ (u r ) is the standard deviation profile of the azimuthal velocity and depends (for a fixed height) only on r. Here, the brackets · denote the same average as the one performed to RMS(ũ r ). From this characteristic velocity we construct the wind Reynolds number, i.e. Re w = (dσ bulk (u r ))/ν. 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 2011). So for the case of TC flow the strength of the wind is given by the (standard deviation of the) radial or axial velocity ).
In the classical regime of turbulence, the unifying theory of Grossmann & Lohse (2011) predicts a scaling of the Reynolds wind Re w ∝ Ta 3/7 . When the driving is increased, towards the ultimate regime, the logarithmic corrections remarkably cancel out and an effective scaling of Re w ∝ Ta 1/2 is observed (Grossmann & Lohse 2011;Huisman et al. 2012). In figure 6(b), we plot the compensated Reynolds wind with the scaling of the ultimate regime Re w Ta −1/2 as a function of Ta. Here, we see that, indeed, beyond the transition to the ultimate regime, Re w slowly asymptotes to a Ta −1/2 scaling for both peaks, which is consistent with the observation of Huisman et al. (2012) for η = 0.716 at a = 0.
We now draw the attention to the case of Ta = 5.10 × 10 8 , which is slightly above 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 ≤ z/L ≤ 0.48. This unexplained observation will be revisited in § 5.3, as we unveil more data from the numerical simulations. Nu ω Ta -1/3 × 10 3 FIGURE 7. Compensated Nusselt number as a function of Ta for three rotation ratios. A reordering of the curves around Ta ≈ 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 Nu ω ∝ Ta 0.4 is given as reference.

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 Ta = 2.23 × 10 9 and Ta = 3.31 × 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.

Disappearance of the broad peak
We first focus on the changing shape of the Nu ω (a) relationship at Ta ≈ 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 Nu ω (Ta) relationship appeared at Ta ≈ 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 Ta. In figure 7, we show the behaviour of the Nu ω (Ta) curves for three selected values of a. We see how at Ta ≈ 10 10 , a reordering of the curves occurs, distinct from the one seen at Ta = 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.
By relating this to the transition observed in Ostilla-Mónico et al. (2014a), 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 Ω (and a). As R Ω 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 θ 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. (2014), where the regions with streaks were classified as plume emitting, and the regions with no streaks were classified as plume impacting.
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 η, 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 Ta, beyond 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 Nu ω (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 η = 0.714 in Ostilla-Mónico et al. (2014), where the disappearance of quiescent regions and the transition to the ultimate regime happens simultaneously.

Roll state switches
We now focus on the discrete jumps in the Nu ω (a) curve at Ta = 2.23 × 10 9 and Ta = 3.31 × 10 9 for a < 0 from the experiments (see figure 3). We note that in previous simulations, both by Brauckmann & Eckhardt (2017) and Ostilla-Mónico et al. (2014b), a small Γ 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 Nu ω and that the only dependence comes from the roll wavelength periodicity length, both effects will manifest simultaneously. The roll number N and the (non-dimensional) roll wavelength λ z are related by λ z = Γ /N.
In figure 9, we show the compensated Nusselt number for both the small and the medium boxes for two Ta, and the experiments for Ta = 5.10 × 10 8 . The numerical simulations at Ta = 5.10 × 10 8 reveal that for a < −0.5 the value of Nu ω is at most weakly dependent on Γ . However, for a ≥ −0.5 the two numerical studies and the experimental measurements result in different values of Nu ω . We note that Nu ω 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. (2016).
To further explore this, in figure 10 we show the azimuthally and temporally averaged radial velocity for the medium-sized domain (Γ = 12.56), as we vary a from the corotating to the counter-rotating regime. In terms of R Ω 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 Ω / = 0 (Sacco et al. 2019). As R Ω 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 2011; van der Poel et al. 2012). As Ta increases, these sharp changes in the Nu ω (a) curve disappear, because roll state switching does not modify the fraction of boundary layer regions which emit plumes.
The influence of the large structures on the torque can also be appreciated when we compare the Γ = 12.56 simulation with the experimental data (Γ = 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 (Γ = 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 Γ is much larger, so state switches will generally have a very   FIGURE 10. Azimuthally time-averaged radial velocity normalized with the characteristic velocity U obtained from DNS for Ta = 5.10 × 10 8 and Γ = 12.56. Each plot represents a different rotation state, quantified by either a or R Ω . The legend on top of each figure represents the value of (a, R Ω ). 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 Ω . small effect on the roll wavelength, and as a consequence a small effect on the torque. This can explain both the jumps in the Nu ω (a) curve, and why they are much smaller than for numerics.

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 Ω 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-θ ) 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 Ω 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 Ω 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 Ω over a certain number of large eddy turnover times τ = tU/d, where U is the characteristic velocity as defined in § 3. The simulations are performed for Ta = 5.10 × 10 8 and for the medium-sized box (Γ = 12.56). The idea is to capture the evolution of the radial velocity and its effect on Nu ω as we slowly increase R Ω in discrete steps. We initiate the simulation at τ = 0 with a statistical stationary TC flow and R Ω = 0. Next, at τ = 25, we increase the value of the Coriolis force to R Ω = 0.02. When τ = 50 is reached, we set finally R Ω to 0.03. In figure 11(a), we show the instantaneous Nu ω 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 Nu ω and the change in R Ω 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 ≈10τ , 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.

Summary and Conclusions
We probe the angular momentum transport with both experiments and direct numerical simulations for η = 0.91 as a function of the driving which we quantify with Ta and the magnitude of the Coriolis force which we quantify with either the rotation ratio a or the rotation number R Ω . The range of shear driving we explore spans several decades of Ta, namely O(10 7 )-O(10 10 ), which includes the transition to the ultimate regime at Ta ≈ Ta c = 3 × 10 8 . In the vicinity of Ta c , our experiments reveal the presence of two local angular momentum transport maxima for fixed Ta, as was firstly reported by the numerical works of Brauckmann et al. (2016) and Brauckmann & Eckhardt (2017). 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 Ω , 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 Ω while decreasing its magnitude. Conversely, the narrow peak increases its magnitude while experiencing a shift in R Ω . We attribute this to the different way in which the near-wall structures self-organizes at Ta ≈ 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 Nu ω (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 Ta = 4.95 × 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 Ω = 0.141) and Ta = 5.10 × 10 8 , which is close to 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 Ω . The numerics at Ta = 5.10 × 10 8 reveal that, as the magnitude of the Coriolis force increases via R Ω , 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. Nu ω , which is evidenced by different Nu ω (a) curves for different values of Γ . We explore the reorganization of the rolls in detail by performing a last numerical simulation where we slowly vary R Ω several times. Here, we observe that for the same driving Ta and curvature η as in the experiments (different Γ ), 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.