Statistics, plumes and azimuthally travelling waves in ultimate Taylor–Couette turbulent vortices

Abstract

In this paper, we experimentally study the influence of large-scale Taylor rolls on the small-scale statistics and the flow organization in fully turbulent Taylor–Couette flow for Reynolds numbers up to $Re_{S}=3\times 10^{5}$. The velocity field in the gap confined by coaxial and independently rotating cylinders at a radius ratio of $\unicode[STIX]{x1D702}=0.714$ is measured using planar particle image velocimetry in horizontal planes at different cylinder heights. Flow regions with and without prominent Taylor vortices are compared. We show that the local angular momentum transport (expressed in terms of a Nusselt number) mainly takes place in the regions of the vortex in- and outflow, where the radial and azimuthal velocity components are highly correlated. The efficient momentum transfer is reflected in intermittent bursts, which becomes visible in the exponential tails of the probability density functions of the local Nusselt number. In addition, by calculating azimuthal energy co-spectra, small-scale plumes are revealed to be the underlying structure of these bursts. These flow features are very similar to the one observed in Rayleigh–Bénard convection, which emphasizes the analogies of these systems. By performing a complex proper orthogonal decomposition, we remarkably detect azimuthally travelling waves superimposed on the turbulent Taylor vortices, not only in the classical but also in the ultimate regime. This very large-scale flow pattern, which is most pronounced at the axial location of the vortex centre, is similar to the well-known wavy Taylor vortex flow, which has comparable wave speeds, but much larger azimuthal wavenumbers.

1 Introduction

The flow in between two independently rotating cylinders, known as Taylor–Couette (TC) flow, is a commonly used model for general rotating shear flows. It features rich and diverse flow states, which have been explored for nearly a century (Taylor 1923; Wendt 1933; Coles 1965; Andereck, Liu & Swinney 1986; Donnelly 1991). More recent reviews on the hydrodynamic instabilities can be found in Fardin, Perge & Taberlet (2014), and on fully turbulent Taylor–Couette flows in Grossmann, Lohse & Sun (2016). The geometry of a Taylor–Couette system is defined by the gap width $d=r_{2}-r_{1}$ , where $r_{1}$ and $r_{2}$ are the inner and outer radii respectively; the radius ratio $\unicode[STIX]{x1D702}=r_{1}/r_{2}$ , and the aspect ratio $\unicode[STIX]{x1D6E4}=\ell /d$ , with $\ell$ the height of the cylinders. The external driving of the flow can be quantified by the shear Reynolds number according to (Dubrulle et al. 2005)

(1.1) $$\begin{eqnarray}Re_{S}={\displaystyle \frac{2r_{1}r_{2}d}{(r_{1}+r_{2})\unicode[STIX]{x1D708}}}|\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1}|={\displaystyle \frac{u_{S}d}{\unicode[STIX]{x1D708}}},\end{eqnarray}$$

with the cylinder angular velocities $\unicode[STIX]{x1D714}_{1,2}$ , $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the fluid and $u_{S}$ is the shear velocity. A further dimensionless control parameter is the ratio of the angular velocities

(1.2) $$\begin{eqnarray}\unicode[STIX]{x1D707}={\displaystyle \frac{\unicode[STIX]{x1D714}_{2}}{\unicode[STIX]{x1D714}_{1}}},\end{eqnarray}$$

implying $\unicode[STIX]{x1D707}>0$ for co-rotation of the cylinders, $\unicode[STIX]{x1D707}=0$ for pure inner cylinder rotation and $\unicode[STIX]{x1D707}<0$ for counter-rotation. The most important response parameter of the TC system to the cylinder driving is the angular velocity transport (Eckhardt, Grossmann & Lohse 2007a )

(1.3) $$\begin{eqnarray}J_{\unicode[STIX]{x1D714}}=r^{3}(\langle u_{r}\unicode[STIX]{x1D714}\rangle _{\unicode[STIX]{x1D711},z,t}-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{r}\langle \unicode[STIX]{x1D714}\rangle _{\unicode[STIX]{x1D711},z,t}),\end{eqnarray}$$

where $r$ denotes the radial coordinate, $\unicode[STIX]{x1D711}$ the azimuthal coordinate, $t$ the time coordinate and $\langle \cdot \rangle _{\unicode[STIX]{x1D711},z,t}$ the azimuthal–axial time average. This quantity is conserved along $r$ and can be directly measured by the torque ${\mathcal{T}}$ acting on either the inner (IC) or the outer cylinder wall (OC). Normalizing $J_{\unicode[STIX]{x1D714}}$ with its corresponding laminar non-vortical value $J_{\unicode[STIX]{x1D714}}^{lam}=2\unicode[STIX]{x1D708}r_{1}^{2}r_{2}^{2}(\unicode[STIX]{x1D714}_{1}-\unicode[STIX]{x1D714}_{2})/(r_{2}^{2}-r_{1}^{2})$ yields a quasi-Nusselt number $\mathit{Nu}_{\unicode[STIX]{x1D714}}=J_{\unicode[STIX]{x1D714}}/J_{\unicode[STIX]{x1D714}}^{lam}$ (Eckhardt et al. 2007a ), which is analogous to the Nusselt number $\mathit{Nu}$ in Rayleigh Bénard flow (RB) flow, i.e. the buoyancy-driven flow which is heated from below and cooled from above. There, $\mathit{Nu}$ is a measure for the amount of transported heat flux normalized by the purely conductive heat transfer. Eckhardt, Grossmann & Lohse (2007b ) worked out the fundamental similarities between TC and RB flow in terms of the Nusselt number and the energy dissipation rate, which we will use in this paper.

The dependence of the Nusselt number on the shear Reynolds number, commonly expressed in terms of an effective power law $\mathit{Nu}_{\unicode[STIX]{x1D714}}\sim Re_{S}^{\unicode[STIX]{x1D6FC}}$ , and on the rotation ratio $\unicode[STIX]{x1D707}$ has been widely investigated (Lathrop, Fineberg & Swinney 1992; Lewis & Swinney 1999; van Gils et al. 2011; Paoletti & Lathrop 2011; van Gils et al. 2012; Merbold, Brauckmann & Egbers 2013; Ostilla-Mónico et al. 2014c ; Brauckmann, Salewski & Eckhardt 2016a ; Grossmann et al. 2016). For pure inner cylinder rotation ( $\unicode[STIX]{x1D707}=0$ ), a change in the local scaling exponent $\unicode[STIX]{x1D6FC}$ with increasing $Re_{S}$ has been found which is caused by a transition from laminar (classical regime) to turbulent boundary layers (ultimate regime) (Ostilla-Mónico et al. 2014c ). The transition point depends on the radius ratio $\unicode[STIX]{x1D702}$ and is located around $Re_{S,crit}\approx 1.6\times 10^{4}$ for $\unicode[STIX]{x1D702}=0.714$ . In the ultimate regime, the scaling exponent becomes $\unicode[STIX]{x1D6FC}\approx 0.76$ independent of $\unicode[STIX]{x1D702}$ (Ostilla et al. 2013). When the driving ( $Re_{S}$ ) is kept constant and only $\unicode[STIX]{x1D707}$ is changed, the TC flow features a maximum in the angular momentum transport ( $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ ). For $\unicode[STIX]{x1D702}=0.714$ , the maximum is located in the counter-rotating regime at $\unicode[STIX]{x1D707}_{max}\approx -0.36$ and originates from a strengthening of the turbulent Taylor vortices (Brauckmann & Eckhardt 2013b ; Ostilla-Mónico et al. 2014b ). The contribution of these rolls to the angular momentum transport has been evaluated numerically by Brauckmann & Eckhardt (2013b ) and experimentally by Froitzheim, Merbold & Egbers (2017) by means of the decomposition of $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ into its turbulent fluctuation and large-scale circulation contributions. They find that the large-scale contribution dominates the transport in the region of the torque maximum. Besides, Tokgoz et al. (2011) could show by performing direct torque measurements and tomographic particle image velocimetry (PIV) measurements in a TC facility at $\unicode[STIX]{x1D702}=0.917$ and $Re_{S}\in [1.1\times 10^{4};2.9\times 10^{4}]$ that the torque is strongly affected by the rotation ratio, which determines whether large-scale or small-scale structures are dominant in the flow. These findings reflect that turbulent Taylor vortices can play a prominent role in the fully turbulent regime.

Hence, the morphology and physical mechanisms behind these roll structures have been in the focus of different studies throughout the literature. An interesting phenomenon regarding Taylor vortices is the reappearance of azimuthal waves in the turbulent Taylor vortex regime. Walden & Donelly (1979) measured the point-wise radial velocity component close to the OC boundary layer for $\unicode[STIX]{x1D707}=0$ at $\unicode[STIX]{x1D702}=0.875$ and for different aspect ratios. They find a regime of reappearance for $28\leqslant Re/Re_{C}\leqslant 36$ and for $\unicode[STIX]{x1D6E4}\geqslant 25$ , based on sharp peaks in the power spectrum. Here, $Re_{C}$ is the critical Reynolds number for the onset of Taylor vortex flow. Later, Takeda (1999) acquired time-resolved axial profiles of the axial velocity component using an ultrasonic measurement technique for $\unicode[STIX]{x1D702}=0.904$ and $\unicode[STIX]{x1D6E4}=20$ . The azimuthal waves are identified based on Fourier analysis and proper orthogonal decomposition (POD) in the range of $23\leqslant Re/Re_{C}\leqslant 36$ . Their results show good agreement with those of Walden & Donelly (1979). In another study, Wang et al. (2005) performed planar PIV measurements in a meridional plane for $\unicode[STIX]{x1D702}=0.733$ and $\unicode[STIX]{x1D6E4}=34$ . They capture the reappearance of azimuthal waves for $20\leqslant Re/Re_{C}\leqslant 38$ based on spatial correlations. Note that the three aforementioned studies are all based on pure inner cylinder rotation ( $\unicode[STIX]{x1D707}=0$ ). More recently, Merbold, Froitzheim & Egbers (2014) performed flow visualizations in TC flow with $\unicode[STIX]{x1D702}=0.5$ at $Re_{S}=5000$ . They find an axial oscillation of the turbulent Taylor vortices in the range of $\unicode[STIX]{x1D707}\in [-0.15,-0.3]$ , which includes the rotation ratio for optimum transport $\unicode[STIX]{x1D707}_{max}=-0.2$ . In summary, the large-scale turbulent Taylor rolls seem to feature an instability mechanism similar to the one in the laminar regime, which, however, has not yet been detected in highly turbulent TC flows. Based on numerical simulations, Ostilla et al. (2013), Ostilla-Mónico et al. (2014c ) further showed that the large-scale rolls consist and are driven by small-scale unmixed plumes, in analogy to RB flow. They calculated for $\unicode[STIX]{x1D707}=0$ and $\unicode[STIX]{x1D702}=0.714$ the radial profiles of the angular velocity $\unicode[STIX]{x1D714}$ at specific axial positions of the large-scale Taylor vortices; namely, at the vortex inflow, vortex centre and vortex outflow. The vortex inflow is characterized by the ejection of plumes from the OC in conjunction with a mean radial velocity that points away from the OC (see the sketch in figure 1 a). In contrast, the outflow features plume ejections from the IC with a mean radial velocity component directed from the IC to the OC (see the sketch in figure 1 c). In between the in- and outflow, the radial velocity component becomes zero in the middle of the gap, which denotes the location of the centre of the vortex, as shown in figure 1(b). In regions where plumes are ejected from the cylinder walls, i.e. in the vortex inflow at the outer wall and in the vortex outflow at the inner wall, the radial profiles of $\unicode[STIX]{x1D714}$ have a logarithmic shape in the corresponding boundary layer (Ostilla-Mónico et al. 2016). It is worth mentioning that in the ultimate regime, the profiles become logarithmic also in the absence of dominant large-scale rolls (Huisman et al. 2013b ). For $\unicode[STIX]{x1D702}=0.5$ in the classical regime, van der Veen et al. (2016b ) calculated the velocity of such plumes based on planar PIV measurements performed at different heights, which further confirms the connection between small-scale plumes and large-scale vortices.

Figure 1. Sketch of the ejecting regions bounded by the inner (IC) and outer (OC) cylinders within a Taylor roll: (a) vortex inflow, (b) vortex centre, (c) vortex outflow.

Within the context of both TC and RB flow, many studies in the literature focus on establishing a distinct connection between the transport of angular momentum (or heat in RB flow) and the structures inherent to these flows, both for large-scale rolls and small-scale plumes. A comparative study of the probability density functions (PDFs) of the Nusselt number in TC and RB calculated over cylindrical surfaces and horizontal planes respectively has been performed by Brauckmann, Eckhardt & Schumacher (2016b ). As a reference point for the comparison, they choose $Re_{S}=2\times 10^{4}$ for TC and a Rayleigh number of $Ra=\unicode[STIX]{x1D6FC}_{P}gH^{3}\unicode[STIX]{x0394}T/(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708})=10^{7}$ for RB, where the Nusselt number is identical for both flows. Here, $Ra$ is the dimensionless driving parameter in RB flow with $\unicode[STIX]{x1D6FC}_{p}$ the thermal expansion coefficient, $g$ the gravitational acceleration, $H$ the height of the RB cell, $\unicode[STIX]{x0394}T$ the temperature difference and $\unicode[STIX]{x1D705}$ the thermal diffusivity. They find that the PDFs of the net transport have the same asymmetric shape with differences in the width of the tails in the boundary layer regions. These differences can be attributed to different shapes and detachment frequencies of the plumes. Moreover, the PDFs of the angular momentum and temperature fluctuations depict a cusp-like shape with pronounced exponential tails as a consequence of the effect of intermittent bursting plumes. Similar analyses of heat flux PDFs in RB flow have been performed by Shishkina & Wagner (2007). They find that the instantaneous heat flux fluctuates around zero and not around the volume-averaged Nusselt number, along with a broadening of the tails with increasing $Ra$ . Shang et al. (2004) revealed that the asymmetry of the PDFs of $\mathit{Nu}$ , which mainly occurs in the tails, arises from correlated temperature and velocity signals produced by thermal plumes. These plumes lead to large but rare positive events of heat flux. However, the results of Shang et al. (2004) are only based on point-wise measurements. For TC flow, the statistics of $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ were analysed for $\unicode[STIX]{x1D707}=0$ by Huisman et al. (2012) without any connection to specific flow structures.

Another approach investigating structures in TC flow is to analyse the kinetic energy spectra. Dong (2007) performed direct numerical simulations (DNSs) for $\unicode[STIX]{x1D702}=0.5$ , $Re=8000$ and the OC at rest. Strikingly, he finds a small-scale peak in the axial spectra of the radial velocity component. According to Dong (2007) the underlying structures can be specified as herringbone streaks. The DNSs of Ostilla-Mónico et al. (2016) in the boundary layer regions reveal a peak in the azimuthal and axial spectra of the radial velocity component at large wavenumbers, which indicates the existence of small-scale plumes. Their simulations were done for $\unicode[STIX]{x1D702}=0.909$ and $\unicode[STIX]{x1D707}=0$ at $Re_{S}\geqslant 10^{5}$ . We stress that the energy spectra in TC flow do not follow the classical Kolmogorov scaling for homogeneous and isotropic turbulence (HIT) of $-5/3$ (Lewis & Swinney 1999; Dong 2007; van Hout & Katz 2011; Huisman, Lohse & Sun 2013a ; Ostilla-Mónico et al. 2016).

Based on this literature review, the following open questions are addressed within this manuscript: how do turbulence-dominated and vortex-dominated flow states differ with respect to their velocity field statistics, how important are the vortex inflow, vortex centre and vortex outflow regions for the momentum transport, how do small-scale structures affect this transport, what is the length scale of these structures and do wavy-vortex-like turbulent Taylor vortices exist in the ultimate turbulent regime? To answer these questions, we make use of PIV measurements in horizontal planes at different cylinder heights for $\unicode[STIX]{x1D702}=0.714$ , in the range of $Re_{S}\in [9.3\times 10^{3},3.5\times 10^{5}]$ and $\unicode[STIX]{x1D707}\in [0,-0.36]$ . We want to stress that such quasi-three-dimensional experimental investigations of the local angular momentum transport statistics and flow structures in the ultimate turbulent TC flow at $\unicode[STIX]{x1D707}_{max}$ are unique, while the measurement set-up consisting of a TC apparatus with a transparent top plate and a horizontal PIV configuration has already been used successfully by van der Veen et al. (2016a ) and Froitzheim et al. (2017).

The paper is organized as follows. In § 2, the experimental set-up, the measurement technique and the investigated parameter space are described in detail. Thereafter, the flow states, statistical profiles and velocity PDFs are shown and compared to the literature to prove the quality of the measurements and discuss the influence of large-scale turbulent Taylor vortices on the global flow statistics (§ 3). In § 4 the global and local angular momentum transport are analysed based on the net convective Nusselt number and the contributions of the vortex in- and outflow to the overall transport are worked out. To detect intermittent bursting small-scale structures which influence this transport, the PDFs of the net convective Nusselt number are evaluated over cylindrical surfaces as well as at the axial height of the vortex inflow, centre and outflow in § 5. The energy content and azimuthal length scale of these structures are calculated in § 6 based on azimuthal energy co-spectra, while azimuthally travelling waves superimposed on the turbulent Taylor vortices are extracted from the flow field based on a complex proper orthogonal decomposition in § 7. The paper ends with a summary and a conclusion (§ 8).

2 Experimental set-up

The PIV experiments were performed in the boiling Twente Taylor–Couette facility (BTTC) at the University of Twente. The BTTC is an ideal facility to perform PIV experiments due to its transparent outer cylinder and top plate. The inner and outer radius of the set-up are $r_{1}=75~\text{mm}$ and $r_{2}=105~\text{mm}$ , respectively, and thus the radial gap is $d=r_{2}-r_{1}=30~\text{mm}$ . The height of the cylinders is $\ell =549~\text{mm}$ , which gives an aspect ratio of $\unicode[STIX]{x1D6E4}=\ell /d=18.3$ . The radius ratio is then $\unicode[STIX]{x1D702}=r_{1}/r_{2}=0.714$ . A more detailed overview of the set-up can be found in Huisman et al. (2015). The flow consists of water whose viscosity $\unicode[STIX]{x1D708}$ and density $\unicode[STIX]{x1D70C}$ can be controlled throughout the experiments due to the temperature control of the BTTC. We fix the temperature of the experiments to $20\,^{\circ }\text{C}$ , leading to $\unicode[STIX]{x1D708}=1.002~\text{mm}^{2}~\text{s}^{-1}$ and $\unicode[STIX]{x1D70C}=0.998~\text{g}~\text{cm}^{-3}$ . The standard deviation of the temperature is 15 mK. The flow is seeded with fluorescent polyamide particles with diameters up to $20~\unicode[STIX]{x03BC}\text{m}$ with an average particle density of approximately $0.01\,\text{particles}/\text{pixel}$ . These particles are coated with Rhodamine B which has a maximum emission centred at approximately $565\,\text{nm}$ . The illumination is provided by a laser sheet from a double-pulsed cavity laser (Quantel Evergreen 145 laser, 532 nm). The thickness of the laser sheet is ${\approx}1~\text{mm}$ . The laser is mounted on a traverse system (Dantec lightweight traverse) which allows us to precisely change the location of the laser sheet along the vertical direction. The camera we use for the recordings is an Imager sCMOS ( $2560\times 2160~\text{px}$ ) 16 bit with a Carl Zeiss Milvus $2.0/100$  lens. We capture the velocity fields with a framerate of $f=15~\text{Hz}$ . Since the camera is operated in double frame mode, we can have very small interframe times, i.e. $\unicode[STIX]{x0394}t\ll 1/f$ . In order to maximize the contrast of the images, we use a long-pass filter in front of the lens (Edmund High-Performance Longpass Filter, $550~\text{nm}$ ), which collects only the emitted light from the fluorescent particles. In figure 2(a), we show a sketch of the experimental set-up.

Figure 2. (a) Sketch of the experimental apparatus (BTTC). Shown is a vertical section. A mirror set at 45 degrees is used so the camera captures the velocity field in the $r-\unicode[STIX]{x1D711}$ plane. The picture shows an exaggeration of the 23 positions of the laser sheet along $\unicode[STIX]{x0394}\ell$ . (b) Temporally averaged flow field in a horizontal plane at the axial height of the vortex centre for $Re_{S}=3.51\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ . Colours represent the azimuthal velocity component and arrows the velocity field. Only each tenth vector is shown. Black solid lines indicate the position of the inner and outer cylinder.

The velocity fields are calculated with commercial software (Davis 8.0) using a multi-pass method. The algorithm uses windows of size $64\times 64~\text{px}$ for the first pass and windows of $24\times 24\,\text{px}$ for the last iteration with a 50 % overlap of the windows. This process yields the velocity fields in Cartesian coordinates. In order to have access to the velocity fields in polar coordinates, we map the Cartesian velocity fields onto a polar grid using bilinear interpolation. The mapping is done such that the radial $\unicode[STIX]{x0394}r$ and azimuthal $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ resolution is the same as the spatial resolution in Cartesian coordinates $\unicode[STIX]{x0394}x$ , i.e. $\unicode[STIX]{x0394}r=\unicode[STIX]{x0394}x$ and $r\unicode[STIX]{x0394}\unicode[STIX]{x1D711}=\unicode[STIX]{x0394}x$ . In this way, the resultant velocity fields are of the form $\boldsymbol{u}=u_{r}(r,\unicode[STIX]{x1D711},t)\boldsymbol{e}_{\boldsymbol{r}}+u_{\unicode[STIX]{x1D711}}(r,\unicode[STIX]{x1D711},t)\boldsymbol{e}_{\unicode[STIX]{x1D711}}$ , where $u_{r}$ and $u_{\unicode[STIX]{x1D711}}$ are the radial and azimuthal velocity components which depend on the radial coordinate $r$ , the azimuthal coordinate $\unicode[STIX]{x1D711}$ and the time coordinate $t$ ; $\boldsymbol{e}_{\boldsymbol{r}}$ and $\boldsymbol{e}_{\unicode[STIX]{x1D711}}$ are the unit vectors in the radial and azimuthal direction, respectively.

In table 1, we present a summary of the measurements we performed. In total, eight cases were investigated which will be addressed in the following sections. Each case contains measurements done at 23 different heights which are separated by $\unicode[STIX]{x0394}z=4~\text{mm}$ , and each height contains 1500 different velocity fields. Thus, the height of the experiments spans a length of $\unicode[STIX]{x0394}\ell =22\unicode[STIX]{x0394}z=88~\text{mm}$ . The resolution of the velocity fields $\unicode[STIX]{x0394}x$ depends on the height but lies within $\unicode[STIX]{x0394}x\in [0.607,0.752]~\text{mm}$ , where the smallest value corresponds to the height closest to the camera at $(z-\ell /2)/d=1.5$ and the smallest to $(z-\ell /2)/d=-1.5$ .

We investigate flow states for pure inner cylinder rotation as a reference for a turbulence-dominated flow and for the rotation ratio that corresponds to the torque maximum, where pronounced large-scale Taylor rolls are present in the gap (see table 1). Further, the measurements are classified based on the discovered change in the local scaling exponent $\unicode[STIX]{x1D6FC}$ by Ostilla-Mónico et al. (2014b ) at $Re_{S}(\unicode[STIX]{x1D702}=0.714)\approx 1.6\times 10^{4}$ for $\unicode[STIX]{x1D707}=0$ , which is caused by a transition of the boundary layers (BLs). Accordingly, flow states at $Re_{S}<1.6\times 10^{4}$ are assumed to be in the so-called classical regime with laminar BLs, while those at $Re_{S}>1.6\times 10^{4}$ are assumed to be in the ultimate regime with turbulent BLs. Cases 1 and 2 are in the classical regime, where, $\unicode[STIX]{x1D707}_{max}$ changes with the shear Reynolds number, which is why $\unicode[STIX]{x1D707}_{max}=-0.15$ is different to the cases in the ultimate regime (Ostilla et al. 2013). For the flow states in the ultimate regime, the torque maximum is located around $\unicode[STIX]{x1D707}_{max}=-0.36$ .

Table 1. Overview of investigated flow states, defined by the regime, the shear Reynolds number $Re_{S}$ and the rotation ratio $\unicode[STIX]{x1D707}$ . The penultimate column depicts the abbreviations for the different flow states and the last column depicts the line styles used in the study. The transition point from the classical to the ultimate regime is located at $Re_{S}\approx 1.6\times 10^{4}$ .

3 Flow states and velocity profiles

The TC flow at $\unicode[STIX]{x1D707}_{max}$ is dominated by turbulent Taylor vortices, while at $\unicode[STIX]{x1D707}=0$ a featureless turbulent flow state develops inside the gap (van Gils et al. 2011; Brauckmann & Eckhardt 2013b ; Ostilla et al. 2013; Huisman et al. 2014). In the following, this previous finding is confirmed within our measurements and statistical profiles and velocity PDFs over cylindrical surfaces are compared to the literature for validation. Therefore, the azimuthally and time-averaged azimuthal velocity components for the 23 investigated heights are shown in figure 3. The tilde symbols denote normalized quantities. The radial coordinate is normalized as $\tilde{r}=(r-r_{1})/(r_{2}-r_{1})$ , where 0 means the location of the inner cylinder ( $r_{1}$ ) and 1 the location of the outer cylinder ( $r_{2}$ ). Further, the azimuthal velocity is normalized using the cylinder speeds as $\tilde{u} _{\unicode[STIX]{x1D711}}=(u_{\unicode[STIX]{x1D711}}-u_{\unicode[STIX]{x1D711},2})/(u_{\unicode[STIX]{x1D711},1}-u_{\unicode[STIX]{x1D711},2})$ , with $u_{\unicode[STIX]{x1D711}}=r\unicode[STIX]{x1D714}$ . In order to normalize the radial velocity component, the shear velocity is used $\tilde{u} _{r}=u_{r}/u_{S}$ .

Figure 3. Representation of the flow states in terms of temporal and azimuthally averaged velocities. The contour plots depict the azimuthal velocity component for (a)  $Re_{S}=9.30\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ), (b) $9.30\times 10^{3}$ and $\unicode[STIX]{x1D707}=-0.15$ ( $\text{C}_{2}$ ), (c) $Re_{S}=2.15\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{6}$ ) and (d) $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ). (e) Axial profile of the radial velocity component for $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ) at $\tilde{r}=0.5$ with marked locations of vortex inflow, vortex centre and vortex outflow. Grey dots represent data points, which are excluded for flow quantities calculated over the axial coordinate.

For both flows in the classical regime (see figure 3 a,b), large-scale rolls are visible in the mean field filling the whole gap. Apparently, the turbulent fluctuations depicted in figure 5 are not strong enough at low shear Reynolds numbers to suppress these large-scale rolls. Further, the rolls are more pronounced at $\unicode[STIX]{x1D707}_{max}=-0.15$ , indicating an increase in strength. In the ultimate regime at $Re_{S}=2.15\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ $\left(\text{C}_{6}\right)$ , the turbulent Taylor rolls disappear and nearly no axial dependence of the azimuthal velocity is visible in the mean field, just as also found numerically (see the phase diagram, Ostilla-Mónico et al. (2014c , figure 6)). When the rotation rate is changed to $\unicode[STIX]{x1D707}_{max}$ , Taylor vortices are formed again, which are much more pronounced than in the classical regime. The contour plot reveals a mushroom-like structure with distinct in- and outflow regions. The axial profile of the radial velocity component corresponding to figure 3(d) is shown in figure 3(e). This representation exemplifies the further analysis. The recorded 23 heights capture more than one vortex pair, which is why we exclude the grey marked data points for flow quantities calculated over the axial coordinate. The axial length of evaluation therefore starts and ends at a vortex centre. In between, we use the minimum of the azimuthally and temporally averaged axial profile of the radial velocity as the location of the vortex inflow and correspondingly, the location of the vortex outflow is obtained with its maximum. The data point which is closest to a value of zero is defined as the axial location of the vortex centre. In figure 4, we show the normalized radial profiles of the azimuthal and radial velocity components averaged over space (cylindrical surfaces) and time ( $\langle \cdot \rangle _{\unicode[STIX]{x1D711},z,t}$ ). The slopes of the profiles in the bulk nearly vanish at $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{max}$ , while for pure inner cylinder rotation they show a small negative slope in good agreement with other studies (Ostilla-Mónico et al. 2014a ; Brauckmann et al. 2016a ; Froitzheim et al. 2017). With increasing shear Reynolds number, the difference of the data points close to the wall from the wall velocity becomes larger, which is due to the steepness of the velocity gradients. The averaged radial profiles of the radial velocity component are nearly zero all over the gap with absolute values smaller than 1 % of the shear velocity $u_{S}$ , as the radial velocity ideally has to vanish when averaged over one vortex pair. The deviation results from the restricted axial resolution of the individual heights.

Figure 4. Radial profiles of the (a) azimuthal and (b) radial velocity component, averaged over space (cylindrical surfaces) and time. Tilde symbols denote normalized quantities. The radial coordinate is normalized as $\tilde{r}=(r-r_{1})/(r_{2}-r_{1})$ , where 0 means the location of the inner cylinder ( $r_{1}$ ) and 1 the location of the outer cylinder ( $r_{2}$ ). The azimuthal velocity is normalized using the cylinder speeds as $\tilde{u} _{\unicode[STIX]{x1D711}}=(u_{\unicode[STIX]{x1D711}}-u_{\unicode[STIX]{x1D711},2})/(u_{\unicode[STIX]{x1D711},1}-u_{\unicode[STIX]{x1D711},2})$ with $u_{\unicode[STIX]{x1D711}}=r\unicode[STIX]{x1D714}$ . To normalize the radial velocity component, the shear velocity $u_{S}=Re_{S}\unicode[STIX]{x1D708}/d$ is used: $\tilde{u} _{r}=u_{r}/u_{S}$ . Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

Figure 5. Radial profiles of the standard deviation of the (a) azimuthal and (b) radial velocity component, calculated over space (cylindrical surfaces) and time. The profiles are normalized by the shear velocity $u_{S}$ . Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

Next, the radial profiles of the standard deviation of the azimuthal and radial velocity component calculated over space (cylindrical surfaces) and time are plotted in figure 5. In terms of the azimuthal velocity component, the radial profiles of the standard deviation show a nearly constant low value in the centre of the gap and increasing values close to the wall. This increase is more pronounced at $\unicode[STIX]{x1D707}_{max}$ ; while for the lowest shear Reynolds number, a peak close to the outer cylinder wall is visible. In the case of the radial velocity component, the standard deviation depicts a maximum in the centre of the gap and decreases towards the cylinder walls. The maximum of $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D711},z,t}(\tilde{u} _{r})$ for both flow cases in the classical regime is noticeably higher than the one in the ultimate regime. The overall shape of the profiles agrees well with the ones for pure inner cylinder rotation of Ezeta et al. (2017), measured with PIV in the same facility at mid-height.

3.1 Probability density function of velocity components

To provide further validation of our measurements and a basis for the investigation of the local PDFs of the angular momentum transport (§ 5), we first analyse the PDFs of the azimuthal and radial velocity components. In figure 6, we show the PDFs of the azimuthal and radial velocity components at $\tilde{r}=0.5$ and $Re_{S}=2.1\times 10^{5}$ for $\unicode[STIX]{x1D707}=0(\text{C}_{6})$ and $\unicode[STIX]{x1D707}=-0.36~(\text{C}_{7})$ as representatives. In addition, the underlying PDFs at the locations of the vortex inflow, vortex centre and vortex outflow are included.

Figure 6. Probability density functions of the (a,b) azimuthal and (c,d) radial velocity component calculated over space (cylindrical surfaces) and time at $\tilde{r}=0.5$ for $Re_{S}=2.1\times 10^{5}$ . Red dashed lines correspond to $\unicode[STIX]{x1D707}=0~(\text{C}_{6})$ and red solid lines to $\unicode[STIX]{x1D707}=-0.36~(\text{C}_{7})$ . Dark grey, grey and bright grey lines indicate the local PDFs for the vortex inflow, vortex centre and vortex outflow, respectively. The black dashed line represents a Gaussian distribution with zero mean and unit variance.

The local PDFs at the specific vortex locations are close to Gaussian distributions, as already shown by Huisman et al. (2013a ) for the azimuthal velocity component measured via laser Doppler velocimetry (LDV) at mid-height and $\tilde{r}=0.5$ . When all data at different heights are included in the PDFs, they still follow a Gaussian shape at $\unicode[STIX]{x1D707}=0$ for both velocity components. At $\unicode[STIX]{x1D707}_{max}$ , however, the PDF of the azimuthal velocity depicts a cusp-like form centred at the origin with approximately exponential tails in accordance with the numerical simulations at a lower Reynolds number of $Re_{S}=2\times 10^{4}$ by Brauckmann et al. (2016a ). There, and in the study of Emran & Schumacher (2008), a nearly identical shape was reported for the temperature PDF in RB flow, which reveals yet another clear evidence of the analogies between both systems, even in the small-scale statistics. In RB flow, the specific PDF shape of the temperature is induced by a combination of bursting plumes and large-scale rolls (Castaing et al. 1989; Yakhot 1989; Procaccia et al. 1991). By using the TC–RB flow analogy, we can identify here a clear fingerprint of plumes transporting angular momentum in TC flow. In addition, the vortex-dominated TC flow in the region of the torque maximum in the fully turbulent regime shows similar behaviour to the flow organization in RB flow. In the case of the radial velocity at $\unicode[STIX]{x1D707}_{max}$ , the PDF also deviates from the ideal Gaussian shape, which can be attributed again to the intermittent bursting plumes.

Figure 7. Probability density functions of the azimuthal velocity component calculated over space (cylindrical surfaces) and time for the radial locations (a) $\tilde{r}=0.1$ , (b) $\tilde{r}=0.3$ and (c) $\tilde{r}=0.5$ . Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

As the PDFs of $u_{r}$ for flow cases at pure inner cylinder rotation depict a nearly Gaussian shape, which is also valid for different radial locations in the bulk, we omit the radial velocity component for the following analysis within this section. In figure 7, we calculate the PDFs of the azimuthal velocity component for different radial locations. When the point of evaluation approaches from the gap centre to the direction of the inner cylinder wall for $\unicode[STIX]{x1D707}_{max}$ , the right-hand tail of the initial cusp-like PDF becomes more pronounced and its width increases with the shear Reynolds number. This change is caused by the dominance of the vortex outflow in the inner gap region, where the mean azimuthal velocity exhibits a large positive value. In addition, very close to the inner wall at $\tilde{r}=0.1$ , both tails become increasingly exponential. This behaviour is even more pronounced at higher $Re_{S}$ , which is another sign of the ejection of coherent plumes from the cylinder wall. Next to these local changes, with decreasing distance to the wall, the global asymmetry of the PDFs seem to increase.

Figure 8. Radial profiles of the (a) skewness and (b) kurtosis of the azimuthal velocity component, calculated over space (cylindrical surfaces) and time. Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

To account for such global properties of the PDFs, we show in figure 8 the radial profiles of skewness and kurtosis of the azimuthal velocity component for different Reynolds numbers. When the outer cylinder is at rest and the ultimate regime is reached, the skewness is close to zero and the kurtosis close to three, confirming the nearly Gaussian shape. The small radius dependent deviations in skewness for large $Re_{S}$ may result from remnants of turbulent Taylor vortices (Lathrop et al. 1992; Huisman et al. 2014; van der Veen et al. 2016a ). At $\unicode[STIX]{x1D707}_{max}$ the skewness increases from the inner cylinder wall to a maximum around $\tilde{r}=0.16$ , then decreases to negative values with a minimum around $\tilde{r}=0.86$ for $Re_{S}=2.14\times 10^{5}$ $\left(\text{C}_{7}\right)$ and then increases again in the direction of the outer cylinder. Furthermore, the absolute skewness increases with the shear Reynolds number. This behaviour is similar to that of the temperature fields in RB flows reported by Emran & Schumacher (2008). The kurtosis profiles at $\unicode[STIX]{x1D707}_{max}$ depict two maxima, one in the inner gap region at $\tilde{r}=0.22$ and one in the outer gap region at $\tilde{r}=0.68$ for the highest $Re_{S}$ . This reflects that the non-Gaussianity of the PDFs is most pronounced in the regions dominated by the vortex in- and outflows due to coherent plumes. In summary, the global and local velocity field statistics of our measurements agree very well with those of the mentioned literature and exceed the state of the art especially for $\unicode[STIX]{x1D707}_{max}$ to higher forcings.

4 Angular momentum transport

As the angular momentum transport, expressed in terms of the quasi-Nusselt number $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ , is strongly influenced by the local and global flow organization inside the gap, it is an appropriate parameter to statistically investigate the existence of flow structures. Therefore, within this section, we first analyse the global angular momentum transport to compare its amount with the literature, and second, we analyse the axial dependent radial profiles of $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ to reveal the most relevant axial locations for the transport. The subsequent results are fundamental for the small-scale statistical analysis of the local momentum transport in the next section and enable new insights into the vortex-dominated momentum transport. The Nusselt number is composed of a convective and a viscous part (Eckhardt et al. 2007a ),

(4.1) $$\begin{eqnarray}\mathit{Nu}_{\unicode[STIX]{x1D714}}=\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c}(r)+\mathit{Nu}_{\unicode[STIX]{x1D714}}^{\unicode[STIX]{x1D708}}(r)={\displaystyle \frac{r^{2}}{J_{lam}}}\langle u_{\unicode[STIX]{x1D711}}u_{r}\rangle _{\unicode[STIX]{x1D711},z,t}-{\displaystyle \frac{\unicode[STIX]{x1D708}r^{3}}{J_{lam}}}\unicode[STIX]{x2202}_{r}\left\langle {\displaystyle \frac{u_{\unicode[STIX]{x1D711}}}{r}}\right\rangle _{\unicode[STIX]{x1D711},z,t}.\end{eqnarray}$$

While the viscous term $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{\unicode[STIX]{x1D708}}$ dominates in the boundary layers, the convective term $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c}$ dominates in the bulk. Since the focus of our investigation is set to the bulk region, we thus neglect the viscous part. Furthermore, as is shown in figure 4(b), the radial velocity component nearly vanishes when averaged over cylindrical surfaces with an axial length of one vortex pair: $\langle u_{r}\rangle _{\unicode[STIX]{x1D711},z,t}\approx 0$ . Therefore, $r^{2}\langle u_{r}u_{\unicode[STIX]{x1D711}}\rangle _{\unicode[STIX]{x1D711},z,t}\approx r^{2}\langle u_{r}^{\prime }u_{\unicode[STIX]{x1D711}}^{\prime }\rangle _{\unicode[STIX]{x1D711},z,t}$ is valid, which means that only the fluctuations of the azimuthal velocity component around its mean profile contribute to the net momentum flux through these cylindrical surfaces (Brauckmann et al. 2016a ). Accordingly, the net convective flux in the bulk flow can be calculated as

(4.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}={\displaystyle \frac{r^{2}}{J_{lam}}}\langle u_{\unicode[STIX]{x1D711}}^{\prime }u_{r}^{\prime }\rangle _{\unicode[STIX]{x1D711},z,t},\quad \text{with}\\ u_{r}^{\prime }(r,\unicode[STIX]{x1D711},z,t)=u_{r}(r,\unicode[STIX]{x1D711},z,t)-\langle u_{r}(r,\unicode[STIX]{x1D711},z,t)\rangle _{\unicode[STIX]{x1D711},z,t},\\ u_{\unicode[STIX]{x1D711}}^{\prime }(r,\unicode[STIX]{x1D711},z,t)=u_{\unicode[STIX]{x1D711}}(r,\unicode[STIX]{x1D711},z,t)-\langle u_{\unicode[STIX]{x1D711}}(r,\unicode[STIX]{x1D711},z,t)\rangle _{\unicode[STIX]{x1D711},z,t}.\end{array}\right\}\end{eqnarray}$$

Figure 9. (a) Azimuthal- and time-averaged profiles of the total Nusselt number indicated by lines and of the net convective Nusselt number indicated by diamonds ( $\unicode[STIX]{x1D707}_{max}$ ) and circles ( $\unicode[STIX]{x1D707}=0$ ). The evaluation is restricted to the interval $0.1\leqslant \tilde{r}\leqslant 0.9$ to exclude the boundary layers and focus on the bulk flow. Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ . (b) Radially averaged total Nusselt number as function of $Re_{S}$ . Error bars represent the standard deviation of the total Nusselt number along the radial coordinate. Data are compared to torque measurements of Lewis & Swinney (1999) ( $\unicode[STIX]{x1D707}=0$ , $\unicode[STIX]{x1D702}=0.724$ ) and DNS of Ostilla-Mónico et al. (2014c ) ( $\unicode[STIX]{x1D707}=-0.4$ , $\unicode[STIX]{x1D702}=0.714$ ).

To avoid confusion, in the following, we will call the Nusselt number $\mathit{Nu}_{\unicode[STIX]{x1D714}}=\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c}+\mathit{Nu}_{\unicode[STIX]{x1D714}}^{\unicode[STIX]{x1D708}}$ shown in (4.1), the total Nusselt number. In figure 9(a), we show the radial profiles of $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ indicated by lines, and the net convective momentum flux $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ , indicated by diamonds ( $\unicode[STIX]{x1D707}_{max}$ ) and circles ( $\unicode[STIX]{x1D707}=0$ ). The analysis is restricted to the bulk region with $\tilde{r}\in [0.1,0.9]$ . A slight dependence of the $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ -profiles on the radial coordinate is observable, which differs partially for the different flow states. The deviation from the predicted conservation of the momentum transport along the radial direction according to Eckhardt et al. (2007a ) is probably due to the finite axial length of our experimental set-up in contrast to their theory, which is based on cylinders of infinite length. However, the effect of the cylinder length is reduced in our set-up by cutting the axial length of evaluation to the size of one vortex pair. Moreover, also the vortex aspect ratio influences the value of $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ as a function of $Re_{S}$ (Huisman et al. 2014; Martínez-Arias et al. 2014), which takes values for the current study in the range of $2.13d{-}2.67d$ . Considering these aspects, the shape of the $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ -profiles is satisfactory. The net convective momentum flux dominates the total Nusselt number and is valid to evaluate the momentum transport in the bulk region (Huisman et al. 2012). Furthermore, in figure 9(b), we show the radially averaged total Nusselt number as a function of $Re_{S}$ , which is also compared to both torque measurements of Lewis & Swinney (1999) ( $\unicode[STIX]{x1D707}=0,\unicode[STIX]{x1D702}=0.724$ ) and DNS of Ostilla-Mónico et al. (2014c ) ( $\unicode[STIX]{x1D707}=-0.4,\unicode[STIX]{x1D702}=0.714$ ). We find a very good agreement with these data, which enables a more detailed analysis of the momentum transport.

Figure 10. Contour plot of the dimensionless (a) convective and (b) net convective local angular momentum transport in a meridional plane for $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ). Black lines indicate the zero line of the depicted quantities. The colour codes are given in the legends. (c) Joint PDF of the radial and azimuthal velocity fluctuations for the same flow state at $\tilde{r}=0.5$ . The colours give the probability (in log scale), see legend.

In order to get deeper insight into the relation between the convective and net convective Nusselt number, we plot in figure 10(a,b) both local quantities in a meridional plane ( $r-z$ plane) for $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ $\left(\text{C}_{7}\right)$ . The local Nusselt number can be much larger than its average value as already noticed by Huisman et al. (2012) and Ostilla-Mónico et al. (2014b ). In the case of the convective Nusselt number $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c}$ (figure 10 a), the momentum transport in the area of the vortex outflow is positive and strongly concentrated in a small axial region, while in the area of the vortex inflow, the transport is negative and less focused. The magnitude of positive momentum flux is twice as large as the negative one, which yields on average a positive net transport. Alternatively, when the net transport is plotted (figure 10 b), a much clearer picture of the transport process is revealed. While in the in- and outflow region the net transport is positive, only in the sheared regions in between negative net transport can be detected. In order to explain this difference, in figure 10(c) we show the joint PDF of the radial and azimuthal velocity fluctuations at $\tilde{r}=0.5$ . In the inflow region, where the fluid is transported strongly in the negative radial direction, the azimuthal velocity depicts a high probability for negative fluctuation values. As a consequence, the net transport has to be positive. In the outflow region, both velocities are mainly positive, resulting also in a positive correlation. However, when the mean azimuthal velocity component is not subtracted, the joint PDF is shifted to the right, leading to negative correlations in the inflow region (not shown). In summary, the difference between $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c}$ and $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ is the transport of the mean azimuthal velocity by the turbulent Taylor vortices, which vanishes when averaged over cylindrical surfaces. By neglecting this fraction, a much clearer picture of the transport process is revealed. In addition, the representation of the Nusselt number in the meridional plane demonstrates the importance of the in- and outflow regions for the net convective transport.

Figure 11. Radial profiles of the contribution of the net convective momentum transport of the (a) inflow and (b) outflow to the total transport; $n_{z}$ represents the number of heights included into the average over cylindrical surfaces. Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

The contribution of the vortex in- and outflow as a function of the radial location and $Re_{S}$ is depicted in figure 11. As a global feature, the contribution of the vortex inflow to the total net convective momentum transport is especially pronounced in the outer gap region ( $\tilde{r}>0.5$ ) while the opposite is true for the outflow. For the two flow states in the classical regime, where the values of the total Nusselt numbers are comparable (see figure 9), the contributions of the in- and outflow are much more dominant at $\unicode[STIX]{x1D707}_{max}$ . This reflects the strong correlation of the enhanced momentum flux at $\unicode[STIX]{x1D707}_{max}$ and the turbulent Taylor vortices. Further, in the ultimate regime at $\unicode[STIX]{x1D707}=0$ , again the effect of remnants of these large vortices becomes visible in the slight dependence of the Nusselt number on $\tilde{r}$ . When the ultimate regime is reached at $Re_{S}=2.98\times 10^{4}$ at $\unicode[STIX]{x1D707}_{max}$ $\left(\text{C}_{3}\right)$ , the contributions of the vortex in- and outflow are small and become significant again at $Re_{S}=6.68\times 10^{4}$ $\left(\text{C}_{4}\right)$ . For even higher shear Reynolds numbers, the contribution of the inflow stays nearly constant with a fraction slightly above 30 % in the outer gap region; while the outflow contribution continuously increases up to approximately 60 % in the inner gap region. This value is strikingly high and demonstrates that at very high Reynolds numbers, the net convective transport shrinks to very small axial regions, where most of the momentum transport takes place. We would like to encourage further numerical or experimental studies, to confirm this finding.

5 Probability density function of the net convective angular momentum flux

With the knowledge of the total net convective transport and the relative contributions of the in- and outflow, we analyse now the PDFs of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ to identify statistical footprints of small-scale plume structures. To properly normalize these PDFs and provide an idealized shape for comparison, we firstly introduce the Gaussian distribution. According to figure 6, the PDFs of the radial and azimuthal velocity components are nearly Gaussian for pure inner cylinder rotation. In that case, the PDF of their product can be described according to Thoroddsen & van Atta (1992) and Chu et al. (1996) based on a Gaussian distribution. The PDF of two jointly Gaussian random variables $x$ and $y$ is given by

(5.1) $$\begin{eqnarray}P(x,y)={\displaystyle \frac{1}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}\sqrt{1-\unicode[STIX]{x1D70C}_{P}^{2}}}}\exp \left[{\displaystyle \frac{1}{2(1-\unicode[STIX]{x1D70C}_{P}^{2})}}\left({\displaystyle \frac{x^{2}}{\unicode[STIX]{x1D70E}_{x}^{2}}}-{\displaystyle \frac{2xy\unicode[STIX]{x1D70C}_{P}}{\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}}}+{\displaystyle \frac{y^{2}}{\unicode[STIX]{x1D70E}_{y}^{2}}}\right)\right],\end{eqnarray}$$

with the individual standard deviations $\unicode[STIX]{x1D70E}_{x}$ and $\unicode[STIX]{x1D70E}_{y}$ and the correlation coefficient $\unicode[STIX]{x1D70C}_{P}=\langle xy\rangle /(\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y})$ . Equation (5.1) represents an inclined elliptic shape for the joint PDF in contrast to the one shown in figure 10(c). Furthermore, the PDF of the product $z=xy$ is (Thoroddsen & van Atta 1992; Chu et al. 1996)

(5.2) $$\begin{eqnarray}P(z)={\displaystyle \frac{1}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}\sqrt{1-\unicode[STIX]{x1D70C}_{P}^{2}}}}\exp \left({\displaystyle \frac{\unicode[STIX]{x1D70C}_{P}z}{\left(1-\unicode[STIX]{x1D70C}_{P}^{2}\right)\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}}}\right)K_{0}\left({\displaystyle \frac{|z|}{\left(1-\unicode[STIX]{x1D70C}_{P}^{2}\right)\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}}}\right),\end{eqnarray}$$

with $K_{0}$ the modified Bessel function of the second kind. The prefactor of the exponential function in (5.2) is used to normalize the different PDFs of the angular momentum flux, i.e. to standardize the width of their tails. It is worth mentioning that the correlation coefficient $\unicode[STIX]{x1D70C}_{P}$ cannot be normalized in a proper way, which leads to different slopes of the exponential tails of the prediction depending on $\unicode[STIX]{x1D70C}_{P}$ . Therefore, in our calculation, we use the prediction according to (5.2) for the case $Re_{S}=3.51\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{8}$ ).

Figure 12. Probability density functions of the local net convective angular momentum transport calculated over space (cylindrical surfaces) and time for $\tilde{r}=0.5$ . (a) PDFs of all investigated flow states are depicted. (b) PDF for $Re_{S}=2.15\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{6}$ ) with the corresponding PDF at the axial height of vortex inflow, vortex centre and vortex outflow. (c) Same plot as in (b) for $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ). The dashed black line indicates the prediction according to (5.2). Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

Figure 13. (a) Probability density functions of the local net convective angular momentum transport calculated over space (cylindrical surfaces) and time for $Re_{S}=2.14\times 10^{5}$ at $\unicode[STIX]{x1D707}=-0.36$ at different radial positions ( $\text{C}_{7}$ ). (b,c) Same PDFs at $\tilde{r}=0.3$ and $\tilde{r}=0.7$ , respectively, with the corresponding PDFs related to the vortex inflow, vortex centre and vortex outflow.

In figure 12, we show the PDFs of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ , which are normalized with the factor $\unicode[STIX]{x1D70E}_{u_{r}}\unicode[STIX]{x1D70E}_{u_{\unicode[STIX]{x1D703}}}\sqrt{1-\unicode[STIX]{x1D70C}_{P}^{2}}$ in order to compare their shapes, at $\tilde{r}=0.5$ . For pure inner cylinder rotation, the PDFs of the net convective momentum transport agree very well with the proposed jointly Gaussian prediction. However, at $\unicode[STIX]{x1D707}_{max}$ the PDFs become highly skewed due to a change of shape in the positive tails. This deviation results in an increasing but still small number of positive extreme events of momentum flux, which becomes more pronounced with increasing $Re_{S}$ . According to figure 12(c), these rare and extreme events can be almost completely attributed to the vortex in- and especially the vortex outflow due to changes in the azimuthal velocity component (see also figure 6). Here, a strong correlation between the radial and azimuthal velocity component exists due to the coherence of the plumes (see also figure 10 c). In addition, all PDFs have in common that zero is the value with the highest probability instead of the average value of the Nusselt number. The overall shape of the PDFs at $\unicode[STIX]{x1D707}_{max}$ is in good agreement with the findings of Brauckmann et al. (2016b ), whose analysis was however restricted to low Reynolds numbers ( $Re_{S}=2\times 10^{4}$ ) and PDFs along cylinder surfaces. Furthermore, our PDFs are comparable to the corresponding PDFs of the heat flux in the RB flow. Shang et al. (2004) reported similar PDF shapes in the case of RB convection. They argue that such large rare events are footprints of heat flux fluctuations induced by thermal plumes. This feature is obviously shared with TC flows and supports the origin of large-scale turbulent vortices as the result of small-scale unmixed plumes. Even more, slight counter-rotating cylinders at $\unicode[STIX]{x1D707}_{max}$ instead of pure inner cylinder rotation seems to be the right kinematic boundary condition of TC flows for comparisons with the RB flow concerning $\mathit{Nu}_{\unicode[STIX]{x1D714}}$ PDFs.

In figure 13, we show the PDFs of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ for various different radii $r$ for the case of $Re_{S}=2.14\times 10^{5}$ at $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ). The negative tails of the PDFs are largely identical for all investigated radial positions, while the right tails strongly depend on  $\tilde{r}$ . From the inner to the outer cylinder wall, the width of the positive tail and therefore the asymmetry decreases coinciding with degeneration of the exponential tails in the $u_{\unicode[STIX]{x1D711}}$ PDFs in figure 7. As was shown in figure 11(b), the maximum of the outflow contribution to the overall momentum transport is close to $\tilde{r}=0.3$ . At this location, the right tail of the PDF reflects extreme positive events due to the emission of plumes and is nearly covered by the PDF of the outflow (see figure 13 b). However, in the outer gap region at $\tilde{r}=0.7$ , the asymmetric right tail of the PDF is comparably formed by both the in- and outflow region. From our point of view, our findings demonstrate that the comparative evaluation of axially global and local PDFs of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ , which has not been done in TC flow before, is crucial for deeper insights into the statistics of the momentum transport.

6 Azimuthal energy co-spectra and correlations

6.1 Azimuthal energy co-spectra

As shown in the previous section, the net convective Nusselt numbers suggest the presence of small-scale plumes concentrated especially in the in- and outflow regions of the turbulent Taylor vortices, which dominate the transport. Hence, we analyse the spatial energy co-spectra to detect the presence and length scale of small-scale structures in the gap. We assume velocity fluctuations at a constant radial $r_{c}$ and axial $z_{c}$ position in the homogeneous $\unicode[STIX]{x1D711}$ -direction with a total number of azimuthal points of $n_{\unicode[STIX]{x1D711}}=0,1,\ldots ,N-1$ equidistantly spaced by $\unicode[STIX]{x0394}s=\unicode[STIX]{x0394}\unicode[STIX]{x1D711}r_{c}$ i.e.

(6.1) $$\begin{eqnarray}\displaystyle & \displaystyle u_{r}^{\ast }(r_{c},\unicode[STIX]{x1D711},z_{c},t)=u_{r}(r_{c},\unicode[STIX]{x1D711},z_{c},t)-\langle u_{r}(r_{c},\unicode[STIX]{x1D711},z_{c},t)\rangle _{t}, & \displaystyle\end{eqnarray}$$

(6.2) $$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D711}}^{\ast }(r_{c},\unicode[STIX]{x1D711},z_{c},t)=u_{\unicode[STIX]{x1D711}}(r_{c},\unicode[STIX]{x1D711},z_{c},t)-\langle u_{\unicode[STIX]{x1D711}}(r_{c},\unicode[STIX]{x1D711},z_{c},t)\rangle _{t}. & \displaystyle\end{eqnarray}$$

The discrete spatial Fourier transform $U_{r,\unicode[STIX]{x1D711}}$ of both fluctuation components $u_{r,\unicode[STIX]{x1D711}}^{\ast }$ is given by

(6.3) $$\begin{eqnarray}U_{r,\unicode[STIX]{x1D711}}(n_{\unicode[STIX]{x1D711}})=\mathop{\sum }_{k=0}^{N-1}u_{r,\unicode[STIX]{x1D711}}^{\ast }(k)\exp \left(-{\displaystyle \frac{2\unicode[STIX]{x03C0}\text{i}kn_{\unicode[STIX]{x1D711}}}{N}}\right),\end{eqnarray}$$

where for simplicity we have not written out the dependences on $r_{c}$ , $z_{c}$ and $t$ . Thus, the spatial energy co-spectrum $E_{r\unicode[STIX]{x1D711}}$ can be calculated as

(6.4) $$\begin{eqnarray}E_{r\unicode[STIX]{x1D711}}(k_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}})=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{1}{N^{2}}}|U_{r}^{n_{\unicode[STIX]{x1D711}}}\cdot U_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}}|,\quad & \text{for }n_{\unicode[STIX]{x1D711}}=\left[0,{\displaystyle \frac{N}{2}}\right],\\[10.0pt] {\displaystyle \frac{1}{N^{2}}}(|U_{r}^{n_{\unicode[STIX]{x1D711}}}\cdot U_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}}|+|U_{r}^{N-n_{\unicode[STIX]{x1D711}}}\cdot U_{\unicode[STIX]{x1D711}}^{N-n_{\unicode[STIX]{x1D711}}}|),\quad & \text{for }n_{\unicode[STIX]{x1D711}}=\left[1,\ldots ,{\displaystyle \frac{N}{2}}-1\right],\end{array}\right.\end{eqnarray}$$

with the wavenumber vector $k_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}}=(\unicode[STIX]{x0394}s)^{-1}n_{\unicode[STIX]{x1D711}}/N$ . The co-spectra are determined for each time step $t$ and afterwards ensemble averaged over 1500 snapshots for every case. To enable a comparison of the co-spectra for different $Re_{S}$ and radial positions, we normalize all co-spectra with the area under the corresponding graph $A_{E}\approx (2\unicode[STIX]{x0394}s)^{-1}\sum _{n_{\unicode[STIX]{x1D711}}=0}^{N/2-1}[E_{r\unicode[STIX]{x1D711}}(k_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}})+E_{r\unicode[STIX]{x1D711}}(k_{\unicode[STIX]{x1D711}}^{n_{\unicode[STIX]{x1D711}}+1})]$ based on the trapezoidal integration method.

Figure 14. (a) Temporally and azimuthally averaged azimuthal kinetic energy co-spectra evaluated at $\tilde{r}=0.5$ . Spectra are normalized to cover an area of 1 under their respective curves and by the gap width $d$ . (b) Local scaling exponent $\unicode[STIX]{x1D6FE}$ of the co-spectra for $E_{r,\unicode[STIX]{x1D711}}\sim k_{\unicode[STIX]{x1D711}}^{\unicode[STIX]{x1D6FE}}$ calculated with a bin size of $\log _{10}(k_{\unicode[STIX]{x1D711}})=0.5$ . Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

First of all, we show the temporally and axially averaged azimuthal energy co-spectra at $\tilde{r}=0.5$ in figure 14(a) to illustrate the scaling of our spectra and compare it with other studies. The kinetic energy co-spectra show that most of the energy lies within the large scales, corresponding to small wavenumbers. The co-spectra depict a noticeable drop for $k_{\unicode[STIX]{x1D711}}d\approx 20$ . It is worth mentioning that we do not see any peak in the large-scale regime, because of the limited range of the azimuthal coordinate in the experiments. Moreover, we cannot resolve the energy content of axisymmetric Taylor rolls, as they correspond to an azimuthal wavenumber of $k_{\unicode[STIX]{x1D711}}=0$ .

In the case of $Re_{S}$ being in the classical regime, the spectra show, compared to the other flow states, a stronger decrease of the spectral energy at mid scales and a kink in the region of $10\leqslant k_{\unicode[STIX]{x1D711}}d\leqslant 20$ . When $Re_{S}$ is increased into the ultimate regime (with $\unicode[STIX]{x1D707}$ fixed at $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{max}$ ), this kink continuously diminishes and energy is redistributed from the large to the small scales. Based on the power-law ansatz $E_{r,\unicode[STIX]{x1D711}}\sim k_{\unicode[STIX]{x1D711}}^{\unicode[STIX]{x1D6FE}}$ , the local scaling exponent $\unicode[STIX]{x1D6FE}$ corresponding to the co-spectra is shown in figure 14(b); $\unicode[STIX]{x1D6FE}$ is calculated with a bin size of $\log _{10}(k_{\unicode[STIX]{x1D711}})=0.5$ . It is worth mentioning that $\unicode[STIX]{x1D6FE}$ is sensitive to the data processing and the accompanying error propagation; $\unicode[STIX]{x1D6FE}$ is given here for this choice of bin size and to show the trend of $\unicode[STIX]{x1D6FE}$ with $d$ . We observe that the spectra neither show $-1$ nor $-5/3$ scaling, which is consistent with the results of Lewis & Swinney (1999), Ostilla-Mónico et al. (2016) and Huisman et al. (2013a ). In the ultimate regime (with $\unicode[STIX]{x1D707}$ fixed at $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{max}$ ), the exponent decreases with increasing $k_{\unicode[STIX]{x1D711}}$ before it strongly drops down in the viscous regime beyond $k_{\unicode[STIX]{x1D711}}d\approx 20$ . This decrease becomes smaller with increasing $Re_{S}$ .

However, in the case of the highest investigated shear Reynolds number at $Re_{S}=3.51\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{8}$ ), the local exponent is nearly constant for $k_{\unicode[STIX]{x1D711}}d\leqslant 10$ with a value of $\unicode[STIX]{x1D6FE}\approx -1.52$ , which is slightly above $-5/3$ .

Figure 15. Temporally averaged pre-multiplied azimuthal kinetic energy co-spectra in the classical regime for $Re_{S}=9.30\times 10^{3}$ and $\unicode[STIX]{x1D707}=-0.15$ ( $\text{C}_{2}$ ) at different radial positions at the axial height of (a) vortex inflow, (b) vortex centre and (c) vortex outflow. The region of $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ with enhanced plume emission is marked in light blue. Sketches above the co-spectra for the different vortex regions are added for clarity.

Next, we focus on the local pre-multiplied energy co-spectra at the axial height of the vortex inflow, vortex centre and vortex outflow, respectively, where the PDF analysis of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ suggests the occurrence of small-scale intermittent plumes. Pre-multiplied means that the co-spectra are multiplied with the wavenumber vector $k_{\unicode[STIX]{x1D711}}$ , such that the area under its graph corresponds to the kinetic energy (Smits, McKeon & Marusic 2011). In figure 15, we show the pre-multiplied energy co-spectra in the classical regime for $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=-0.15$ ( $\text{C}_{2}$ ) at the three vortex positions. Firstly, we address the peak that corresponds to the large scales. For the vortex inflow in figure 15(a), the large-scale peak is located around $k_{\unicode[STIX]{x1D711}}d\approx 1.1$ within the bulk for $0.2\leqslant \tilde{r}\leqslant 0.8$ , while close to the outer cylinder at $\tilde{r}=0.9$ , the peak lies outside of our resolvable scales. Close to the IC wall, the peak shifts to $k_{\unicode[STIX]{x1D711}}d\approx 2.2$ at $\tilde{r}=0.1$ . In figure 15(c), at the location of the vortex outflow, the opposite behaviour is observed. Here, the large-scale peak shifts to smaller wavenumbers when the radial position increases from the IC to the OC. At the vortex centre in figure 15(b), the large-scale peak is shifted from smaller wavenumbers in the centre of the gap ( $\tilde{r}=0.5$ ) to larger ones at both cylinder walls, in an almost symmetric manner. This suggests that the formation of this large-scale peak is connected to the mean radial velocity field, which is in itself caused by the large-scale turbulent Taylor rolls. When fluid impacts on the cylinder walls due to these rolls, structures of the size $k_{\unicode[STIX]{x1D711}}d\approx 2.2$ are formed. Since the axisymmetric flow has a wavenumber of $k_{\unicode[STIX]{x1D711}}=0$ , the existence of a large-scale peak may be an indication of modified turbulent Taylor vortices. This point will be discussed in more detail in § 7.

With respect to the small-scale peak, we observe that it is located at wavenumbers around $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ for all three depicted heights, independent of the radial coordinate. However, its amplitude strongly varies with $\tilde{r}$ and the height $z$ . In the region of the vortex inflow shown in figure 15(a), the small-scale peak is most pronounced near the OC and its amplitude decreases monotonically with decreasing  $\tilde{r}$ . On the contrary, at the height of the vortex outflow in figure 15(c), the small-scale peak amplitude is largest close to the IC and decreases in amplitude towards the OC. In the region of the vortex centre in figure 15 c), we find the highest amplitude of the small-scale peak in the centre of the gap. This is due to the emission of coherent plumes from the cylinder walls which give rise to the formation of Taylor rolls, as was already mentioned before. Thus, this peak should indeed be most pronounced in the ejecting regions, i.e. in the outflow region at the IC and in the inflow region at the OC, respectively. At the height of the vortex centre – where no predominant flow direction concerning the radial velocity component is present – the behaviour is different: plumes rise from both cylinder walls and travel towards the bulk flow, which results in the highest peak amplitude at $\tilde{r}=0.5$ . Considering the fact that the contribution of the vortex centre to extreme and strong events of momentum flux is almost negligible (see figures 12 and 13), these detected plumes seem to compensate their total radial momentum transport.

Figure 16. (a–c) Temporally averaged pre-multiplied azimuthal kinetic energy co-spectra for $Re_{S}=2.15\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{6}$ ) at different radial positions at the axial height of (a) vortex inflow, (b) vortex centre and (c) vortex outflow. Region of $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ is marked in light blue. (d–f) Same spectra for $Re_{S}=2.14\times 10^{5}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{7}$ ) at the axial height of (d) vortex inflow, (e) vortex centre and (f) vortex outflow.

The pre-multiplied energy co-spectra in the ultimate regime at $Re_{S}=2.1\times 10^{5}$ for the three vortex locations are depicted in figure 16 for $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{6}$ , a–c) and $\unicode[STIX]{x1D707}_{max}$ ( $\text{C}_{7}$ , d–f). For $\unicode[STIX]{x1D707}=0$ , no large-scale peak can be seen in any of the co-spectra at the investigated heights. This is consistent with the finding shown in figure 3(c), where we observed that the Taylor rolls have faded away. Accordingly, the shape of the co-spectra is less dependent on both the axial coordinate and on $\tilde{r}$ . However, the co-spectra show a prominent change in the slope around $k_{\unicode[STIX]{x1D711}}d\approx 20$ : In figure 16(a) (vortex inflow), a small-scale peak is formed around $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ close to the OC at $\tilde{r}=0.9$ . This is comparable to the previous case in the classical regime. Also in the region of the vortex outflow (figure 16 c), we observe a peak within the same range of scales at $\tilde{r}=0.1$ . For the vortex centre however (figure 16 b), no peak is visible.

Also when $\unicode[STIX]{x1D707}$ is changed to $\unicode[STIX]{x1D707}_{max}$ (see figure 16 d,e), we identify a similar behaviour as in the classical regime, although the energy is more homogeneously distributed over all scales. At the height of the vortex inflow as seen in figure 16(d), a large-scale peak is present which shifts to $k_{\unicode[STIX]{x1D711}}d\approx 3.9$ at $\tilde{r}=0.1$ . At the vortex outflow in figure 16(f), this shift appears close to the OC wall at $\tilde{r}=0.9$ for the same wavenumber. In the region of the vortex centre (figure 16 e), the large-scale peak is most pronounced at $\tilde{r}=0.5$ , a smaller wavenumber. For all the heights explored, the co-spectra do not reveal a peak in the small-scale regime. However, at the vortex inflow the energy is redistributed continuously from large to small scales when $\tilde{r}$ is varied from 0.1 to 0.9, as it is clearly visible in the marked regime of $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ . We also observe that the energy contained in the small scales increases in the vortex outflow region with decreasing radial position.

In summary, should prominent turbulent Taylor rolls exist in the flow, a peak at large scales exists and small-scale structures are present throughout the gap, but most prominently in the vortex ejecting regions. However, when the Taylor rolls have faded, the peak at large scales vanishes and small-scale structures are only detectable close to the cylinder walls. These findings reveal yet another clear evidence of the existence of turbulent non-axisymmetric Taylor rolls and support the idea that the large-scale rolls consist of small-scale unmixed plumes. Furthermore, the azimuthal length scale of these plumes is of the order of $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ .

6.2 Spatial correlation coefficients

As was shown in the previous section, the pre-multiplied energy co-spectra demonstrate the existence of small-scale structures (coherent plumes) and large-scale azimuthal structures, which are both connected thanks to the presence of turbulent Taylor vortices. For the sake of clarity, the investigation was confined to three specific cases. Within this section, however, we extend the analysis of the large-scale peak based on the azimuthal two-point autocorrelation coefficient. More specifically, the shift of the large-scale peak close to the cylinder walls as well as its characteristic in the centre of the gap are worked out in more detail. The spatial two-point autocorrelation coefficient of the velocity fluctuations between two points separated by $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ in the azimuthal direction is given by

(6.5) $$\begin{eqnarray}\displaystyle & \displaystyle R_{rr}\left(r,z,\unicode[STIX]{x0394}\unicode[STIX]{x1D711}\right)={\displaystyle \frac{\langle u_{r}^{\ast }\left(r,\unicode[STIX]{x1D711},z,t\right)u_{r}^{\ast }\left(r,\unicode[STIX]{x1D711}+\unicode[STIX]{x0394}\unicode[STIX]{x1D711},z,t\right)\rangle _{\unicode[STIX]{x1D711},t}}{\langle u_{r}^{\ast 2}\left(r,\unicode[STIX]{x1D711},z,t\right)\rangle _{\unicode[STIX]{x1D711},t}}}, & \displaystyle\end{eqnarray}$$

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle R_{\unicode[STIX]{x1D711}\unicode[STIX]{x1D711}}\left(r,z,\unicode[STIX]{x0394}\unicode[STIX]{x1D711}\right)={\displaystyle \frac{\langle u_{\unicode[STIX]{x1D711}}^{\ast }\left(r,\unicode[STIX]{x1D711},z,t\right)u_{\unicode[STIX]{x1D711}}^{\ast }\left(r,\unicode[STIX]{x1D711}+\unicode[STIX]{x0394}\unicode[STIX]{x1D711},z,t\right)\rangle _{\unicode[STIX]{x1D711},t}}{\langle u_{\unicode[STIX]{x1D711}}^{\ast 2}(r,\unicode[STIX]{x1D711},z,t)\rangle _{\unicode[STIX]{x1D711},t}}}. & \displaystyle\end{eqnarray}$$

Figure 17. (a,b) Azimuthal two-point autocorrelation function of the fluctuation radial velocity component at the axial height of vortex in- and outflow at (a) $\tilde{r}=0.1$ and (b)  $\tilde{r}=0.9$ . (c,d) Azimuthal two-point autocorrelation function of the fluctuation (c) azimuthal and (d) radial velocity component at the axial height of vortex centre at $\tilde{r}=0.5$ . Legend abbreviations represent $\text{C}_{\#}|_{Re_{S}}^{\unicode[STIX]{x1D707}}$ .

In figure 17(a), we show the spatial two-point autocorrelation function of $u_{r}$ at $\tilde{r}=0.1$ in the ultimate regime for three different flow cases at $\unicode[STIX]{x1D707}_{max}$ , and for both the vortex inflow and outflow. We observe that, independently of the flow case, the autocorrelation function at the vortex inflow shows a minimum at $r\unicode[STIX]{x1D711}/d\approx 0.13$ , which indicates the existence of azimuthal structures of that size. In contrast, for the outflow, we observe a monotonic decrease of the autocorrelation function. Close to the outer cylinder wall (see figure 17 b), the opposite behaviour is observed with a pronounced minimum of the autocorrelation function for the vortex outflow at $r\unicode[STIX]{x1D711}/d\approx 0.15$ . This suggests that when the large-scale Taylor rolls transport fluid against the cylinder walls, azimuthal structures seem to be stimulated close to these walls in good agreement with the findings of the spectral analysis shown in § 6.1.

Next, in figure 17(c,d), we show both the azimuthal and radial two-point autocorrelation function for the centre of the vortex evaluated at $\tilde{r}=0.5$ . At this position, the large-scale peak in the co-spectra is most pronounced as it was shown in figures 15 and 16. Here, we compare three flow states in the classical and ultimate regime where turbulent Taylor vortices are present and a case where they are absent ( $\text{C}_{6}$ ). These figures reveal that in all vortex dominated cases, a large-scale oscillating behaviour is developed, while in the absence of rolls, the autocorrelation simply decreases towards zero over the whole azimuthal measurement length. We thus conclude, that the large-scale peak in the spectra apparently results from a wavy azimuthal pattern connected to the Taylor rolls.

7 Complex proper orthogonal decomposition (CPOD)

7.1 CPOD method

In order to reveal the flow structure connected to the large-scale oscillation within the flow found in figures 15–17, we perform a proper orthogonal decomposition (POD) of the velocity fluctuations field. The POD, also known as empirical orthogonal function analysis (EOF), is a technique to extract modes like coherent structures from a flow field that contribute most to the energy of the flow (Taira et al. 2017). It can be further used to identify the most relevant degrees of freedom in a dynamical system. For additional information, we refer the reader to the review of Berkooz, Holmes & Lumley (1993). As a classical POD mode cannot capture propagating structures, we perform a complex POD using the Hilbert transform (see Pfeffer et al. 1990), where a mode is split into two patterns namely its real and its imaginary part with a phase difference of $\unicode[STIX]{x03C0}/2$ , based on the algorithm described by Marple (1999). We further use a combined analysis of both velocity components $u_{r}$ and $u_{\unicode[STIX]{x1D711}}$ for the POD analysis. The applied algorithm is based on three steps (Harlander et al. 2011). At first we use the Hilbert transform ${\mathcal{H}}$ , to make the velocity fluctuations complex

(7.1) $$\begin{eqnarray}\displaystyle & \displaystyle u_{r,c}^{\ast }=u_{r}^{\ast }+i{\mathcal{H}}(u_{r}^{\ast }) & \displaystyle\end{eqnarray}$$

(7.2) $$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D711},c}^{\ast }=u_{\unicode[STIX]{x1D711}}^{\ast }+i{\mathcal{H}}(u_{\unicode[STIX]{x1D711}}^{\ast }). & \displaystyle\end{eqnarray}$$

Secondly, the imaginary, transformed velocity fluctuations ( $u_{r,c}^{\ast },u_{\unicode[STIX]{x1D711},c}^{\ast }$ ) are arranged in a data matrix $\unicode[STIX]{x1D63F}$ , where columns are assigned to the spatial grid points ( $1:M$ ) and rows to the time signals ( $1:N$ ), while $M>N$

(7.3) $$\begin{eqnarray}\unicode[STIX]{x1D63F}=\left[\begin{array}{@{}ccc@{}}u_{r,c}^{\ast }(x_{1},t_{1}) & \cdots \, & u_{r,c}^{\ast }(x_{1},t_{N})\\ \vdots & \ddots & \vdots \\ u_{r,c}^{\ast }(x_{M},t_{1}) & \cdots \, & u_{r,c}^{\ast }(x_{M},t_{N})\\ u_{\unicode[STIX]{x1D711},c}^{\ast }(x_{1},t_{1}) & \cdots \, & u_{\unicode[STIX]{x1D711},c}^{\ast }(x_{1},t_{N})\\ \vdots & \ddots & \vdots \\ u_{\unicode[STIX]{x1D711},c}^{\ast }(x_{M},t_{1}) & \cdots \, & u_{\unicode[STIX]{x1D711},c}^{\ast }(x_{M},t_{N})\end{array}\right]\!.\end{eqnarray}$$

At the end, a singular value decomposition (SVD) of the data matrix is performed as $\unicode[STIX]{x1D63F}=\unicode[STIX]{x1D731}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D733}^{\text{T}}$ . The columns of the matrix $\unicode[STIX]{x1D731}$ contain the left singular vectors of $\unicode[STIX]{x1D63F}$ and therefore the $N$ complex modes $\mathbf{CPOD}(r,\unicode[STIX]{x1D711})$ . The diagonal elements of $\unicode[STIX]{x1D72E}$ hold the singular values, whose square is a measure of the kinetic energy captured by the individual modes. Note that the singular values are sorted in descending order together with the complex modes, meaning that the first mode represents a larger fraction of kinetic fluctuation energy of the full field than the second one and so forth. Next, the velocity fluctuations at a specific axial location $z$ can be reconstructed by

(7.4) $$\begin{eqnarray}\boldsymbol{u}^{\ast }(r,\unicode[STIX]{x1D711},t)=\mathop{\sum }_{i=1}^{N}a_{i}(t)\mathbf{CPOD}_{i}(r,\unicode[STIX]{x1D711}).\end{eqnarray}$$

The time-dependent coefficients $a_{i}(t)$ result from a projection of the complex modes onto the data matrix, i.e. $a_{i}=\unicode[STIX]{x1D63F}^{\text{T}}\mathbf{CPOD}_{i}$ . In the following, we will present the CPOD results for two flow states at the height of the vortex centre, where the azimuthal two-point correlation showed a large-scale oscillating behaviour. The first case is in the classical regime at $\unicode[STIX]{x1D707}=0$ and the second one in the ultimate regime at $\unicode[STIX]{x1D707}_{max}$ .

7.2 CPOD in the classical regime

Figure 18. Value of $\mathbf{CPOD}_{1}$ for $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ) at the axial height of the vortex centre for both velocity components. (a) Real and (c) imaginary parts of $\mathbf{CPOD}_{1}$ for the radial velocity component. (b) Real and (d) imaginary parts of $\mathbf{CPOD}_{1}$ for the azimuthal velocity component. Black arrows represent the resulting velocity field of $\mathbf{CPOD}_{1}$ .

The first CPOD mode for $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ) is depicted in figure 18 for both velocity components. With the first temporal coefficient $a_{1}(t)$ and the first mode $\mathbf{CPOD}_{1}$ , the corresponding velocity fluctuation field is reconstructed by $\boldsymbol{u}_{1}^{\ast }(r,\unicode[STIX]{x1D711},t)=a_{1}(t)\mathbf{CPOD}_{1}(r,\unicode[STIX]{x1D711})$ . In figure 18(a), the real part of $\mathbf{CPOD}_{1}$ , representing the radial velocity component, shows nearly circular regions of positive and negative velocity, alternating in the azimuthal coordinate direction. This azimuthal wave pattern becomes azimuthally shifted in the corresponding imaginary part in figure 18(c), indicating an azimuthally travelling wave. The real part of the azimuthal velocity of $\mathbf{CPOD}_{1}$ in figure 18(b), depicts diagonal bands with pointy edges close to the cylinder walls. These regions again show alternating positive and negative velocities. The corresponding imaginary part in figure 18(d) is also azimuthally shifted, showing the same pattern. As a result, it is revealed that the velocity field consists of counter-rotating vortices in the horizontal plane, whose vorticity axes are co-axial to the rotation axis of the system. The pattern appears similar to that of the well-known wavy Taylor vortices (WTV) at lower Reynolds number TC flow. Thus, here we show that in the classical turbulent regime, turbulent Taylor vortices can also feature azimuthal waves.

Figure 19. (a) Turbulent energy fraction captured by the CPOD modes for the case $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ) at the axial height of the vortex centre. Only the first 20 of the total 1500 modes are plotted. (b) Temporal power spectrum of the real part of the temporal coefficient $a_{1}(t)$ of $\mathbf{CPOD}_{1}$ and of the radial velocity component of the full field at $\tilde{r}=0.5$ , averaged over $\unicode[STIX]{x1D711}$ . The frequency $f$ is normalized by the difference of the cylinder frequencies $f_{1}-f_{2}$ .

In figure 19(a), we show the energy fraction captured by the CPOD modes for $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ). The first mode captures approximately 17 % of the total energy fluctuation and we observe a strong drop to the second mode down to 6.7 %. Accordingly, the first CPOD mode is dominant and is, therefore, the only one discussed. The temporal energy spectrum of the real part of the first mode’s temporal coefficient $a_{1}(t)$ is pictured in figure 19(b), together with the temporal energy spectrum of the full field radial velocity component, calculated at $\tilde{r}=0.5$ and averaged over $\unicode[STIX]{x1D711}$ . The dominant frequency of $\mathbf{CPOD}_{1}$ is $f/(f_{1}-f_{2})=3.03$ , identical to that of the radial velocity component. Thus we can assume that the first mode captures the temporal behaviour of the full field quite accurately. The frequencies $f_{1}$ and $f_{2}$ represent the rotation frequencies of the inner and outer cylinder, respectively.

Figure 20. Space–time diagram of (a) the reconstructed radial velocity component based on the first CPOD mode and (b) the full field radial velocity component for $Re_{S}=9.32\times 10^{3}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{1}$ ) at the axial height of the vortex centre as a function of $\unicode[STIX]{x1D711}$ and $t$ . The radial coordinate is fixed at $\tilde{r}=0.5$ .

To illustrate the spatio-temporal character of the first mode, we depict in figure 20(a) its reconstructed radial velocity component as well as the corresponding full field at $\tilde{r}=0.5$ as a function of the azimuthal coordinate $\unicode[STIX]{x1D711}$ and time $t$ . The space–time plot of $\mathbf{CPOD}_{1}$ depicts diagonal bands of alternating positive and negative velocity, representing azimuthally propagating waves into the direction of the mean flow. Furthermore, the same diagonal bands – superimposed by turbulent fluctuations – are observed in the full field. A reappearance of azimuthal waves for pure inner cylinder rotation at a similar radius ratio of $\unicode[STIX]{x1D702}=0.733$ has already been found by Wang et al. (2005) in the Reynolds number regime of $20\leqslant Re/Re_{C}\leqslant 38$ . In our study in the classical regime for $\unicode[STIX]{x1D707}=0$ at $\unicode[STIX]{x1D702}=0.714$ , we find $Re_{S}/Re_{S,C}\approx 99$ , where the critical shear Reynolds number is located around $Re_{S,C}\approx 94.5$ (Ostilla-Mónico et al. 2014c ). This provides evidence of the reappearance of an azimuthal wave in the classical turbulent regime for a much larger turbulence intensity. Moreover, we can extract the corresponding wave frequency $f_{w}$ , wavelength $\unicode[STIX]{x1D706}_{w}$ and phase speed $c_{w}$ of $\mathbf{CPOD}_{1}$ . The temporal phase function $\unicode[STIX]{x1D719}_{i}$ is defined as

(7.5) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{i}(t)=\arctan \left({\displaystyle \frac{\text{Im}(a_{i}(t))}{Re(a_{i}(t))}}\right).\end{eqnarray}$$

Its temporal derivative is equal to the angular frequency of the CPOD mode (Suanto, Zheng & Yan 1998), leading to a wave frequency of $f_{w,1}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}_{1}(t)/(2\unicode[STIX]{x03C0})=1.80~\text{Hz}$ . In addition, the spatial phase function $\unicode[STIX]{x1D6F7}_{i}(x)$ is defined as

(7.6) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{i}(r,\unicode[STIX]{x1D711})=\arctan \left({\displaystyle \frac{\text{Im}(\text{CPOD}_{i}(r,\unicode[STIX]{x1D711}))}{Re(\text{CPOD}_{i}(r,\unicode[STIX]{x1D711}))}}\right),\end{eqnarray}$$

whose azimuthal derivative is a measure for the local wavenumber $k_{\unicode[STIX]{x1D711},w}$ . Note that the CPOD analysis yields only one temporal coefficient $a_{i}(t)$ , but two spatial modes concerning the radial and azimuthal velocity components. Therefore we can extract two wavelengths, which are however nearly identical. The averaged wavelength, evaluated for $\tilde{r}=0.5$ , is given by $\unicode[STIX]{x1D706}_{w,1}=2\unicode[STIX]{x03C0}/(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D711}}\unicode[STIX]{x1D6F7}_{1}(\tilde{r}=0.5,\unicode[STIX]{x1D711}))=0.68\,\text{rad}$ , leading to an azimuthal wavenumber of approximately $k_{\unicode[STIX]{x1D711},w,1}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{w,1}\approx 9$ . In this way, the phase speed becomes $c_{w,1}=\unicode[STIX]{x1D706}_{w,1}f_{w,1}=1.22\,\text{rad}~\text{s}^{-1}=0.35\unicode[STIX]{x1D714}_{1}$ . As the detected wave pattern reminds us of wavy Taylor vortices, we compare our result with the wave speeds for WTV measured by King et al. (1984), although one should keep in mind the large Reynolds number difference. King et al. (1984) find for the case of pure inner cylinder rotation at $\unicode[STIX]{x1D702}=0.73$ , $R/R_{c}\approx 14$ , and an average axial wavelength of $\bar{\unicode[STIX]{x1D706}}_{w}/d=2.4$ a wave speed of $c_{w}\approx 0.2\unicode[STIX]{x1D714}_{1}$ , while the corresponding pattern features $k_{\unicode[STIX]{x1D711},w}=2$ . Additionally, they show for a slightly larger radius ratio of $\unicode[STIX]{x1D702}=0.84$ , that the wave speed increases when $Re/Re_{c}>18$ . Thus, we can conclude that the wave speeds superimposed on the Taylor vortices in the wavy Taylor vortex regime and in the turbulent Taylor vortex regime are of the same order. Moreover, the wave speed seems to increase with the forcing $Re_{S}$ , as the wave speed in our case is noticeably higher than the one reported by King et al. (1984).

7.3 CPOD in the ultimate regime

Figure 21. Value of $\mathbf{CPOD}_{1}$ for $Re_{S}=6.68\times 10^{4}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{4}$ ) at the axial height of the vortex centre for both velocity components. (a) Real and (c) imaginary parts of $\mathbf{CPOD}_{1}$ for the radial velocity component. (b) Real and (d) imaginary parts of $\mathbf{CPOD}_{1}$ for the azimuthal velocity component. Black arrows represent the resulting velocity field of $\mathbf{CPOD}_{1}$ .

We now turn our attention to the ultimate regime at the torque maximum rotation rate. In figure 21(a–d), the first CPOD mode, representing the real and imaginary parts of the radial and azimuthal velocity component, for $Re_{S}=6.68\times 10^{4}$ and $\unicode[STIX]{x1D707}=-0.36$ ( $\text{C}_{4}$ ) at the axial height of the vortex centre is depicted. The flow in the ultimate regime features nearly the same pattern as already shown in figure 18 with alternating regions of positive and negative velocity. In case of the radial and azimuthal velocity component, these regions form circular patches and diagonal bands with pointy edges, respectively. In addition, the patterns of the imaginary part are azimuthally shifted relative to the real parts.

Similar to the previous case in the classical regime, the turbulent energy fraction decreases from the first to the second mode from approximately 12 % to 5 %. Further, the temporal power spectrum of the temporal coefficient $Re(a_{1}(t))$ depicts a prominent peak at $f/(f_{1}-f_{2})=0.92$ (figures are not shown). By use of the temporal and spatial phase functions (see (7.5) and (7.6)), the azimuthally travelling wave can be described by a wave frequency of $f_{w,1}=3.67\,\text{Hz}$ , a wavelength of $\unicode[STIX]{x1D706}_{w,1}=0.74\,\text{rad}$ and a phase speed of $c_{w,1}=2.72~\text{rad}~\text{s}^{-1}=0.11\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . Remarkably, the turbulent Taylor vortices feature azimuthal waves also in the ultimate turbulent regime at $\unicode[STIX]{x1D707}_{max}$ . To the best of our knowledge, this is the first study reporting such a wave pattern in that highly turbulent regime. It is further worth to mention that a CPOD analysis in the ultimate regime for $Re_{S}=2.15\times 10^{5}$ and $\unicode[STIX]{x1D707}=0$ ( $\text{C}_{6}$ ), where the large Taylor rolls disappear in the mean field, does not yield any sign of travelling waves. As a comparison tool, we show in table 2 the results concerning the azimuthally travelling waves found in this study and the findings of King et al. (1984) for $\unicode[STIX]{x1D702}=0.73$ .

Table 2. Overview of detected azimuthally travelling waves in the classical and ultimate regime in comparison to the findings of King et al. (1984).

To what extent the detection of azimuthal travelling waves is dependent on the turbulence level ( $Re_{S}$ ) and the geometry of the system ( $\unicode[STIX]{x1D702},\unicode[STIX]{x1D6E4}$ ) is not known. We note that further experimental as well as numerical investigations would be required to address this issue and reveal whether it is an intrinsic property of the flow or not.

8 Summary and conclusions

By means of planar PIV measurements performed in horizontal planes at different axial heights, we showed the dependence of the small-scale statistics and flow organization on the presence of large-scale Taylor rolls in high Reynolds number TC flow. The ratio of angular velocities $\unicode[STIX]{x1D707}$ is the appropriate parameter to control whether or not the Taylor rolls within the gap are prominent and stable in a statistical sense. While in the classical turbulent regime Taylor rolls are observed for both $\unicode[STIX]{x1D707}=0$ and $\unicode[STIX]{x1D707}<0$ , they fade out when the ultimate regime is reached for $\unicode[STIX]{x1D707}=0$ , but can be formed again when sufficient counter-rotation ( $\unicode[STIX]{x1D707}<0$ ) is introduced, in particular at $\unicode[STIX]{x1D707}_{max}$ where the flux of angular momentum is maximum. This turbulent state can then be used to compare the flow – for the same $Re_{S}$ – with and without the presence of these large-scale rolls.

To uncover the interplay between large-scale rolls and small-scales plumes in a statistical manner, PDFs of both velocity components and the net convective Nusselt number have been evaluated over cylindrical surfaces and at specific axial positions. The PDFs of the radial and azimuthal velocity component are close to Gaussian, when evaluated at a fixed height, in accordance with the results of Huisman et al. (2013a ). However, when all heights covering one vortex pair are included, the Gaussian shape is only preserved for $\unicode[STIX]{x1D707}=0$ for both velocity components but changes drastically for $\unicode[STIX]{x1D707}_{max}$ . There, the PDFs of $u_{\unicode[STIX]{x1D711}}$ feature a cusp-like shape with exponential tails, which is a fingerprint of intermittent small-scale plumes and can also be found in RB flow concerning the temperature (Emran & Schumacher 2008; Brauckmann et al. 2016b ). The PDFs of $\mathit{Nu}_{\unicode[STIX]{x1D714}}^{c,net}$ fluctuate around zero and not the averaged value of the Nusselt number and can be described (for $\unicode[STIX]{x1D707}=0$ ) by the shape of a distribution of the product of two Gaussian variables. At $\unicode[STIX]{x1D707}_{max}$ however, the right-hand tail broadens with increasing $Re_{S}$ , which can be attributed to an increasing number of rare and strong events of momentum flux at the axial location of the vortex in- and outflow. These axially dependent events originate from the ejection of coherent plumes from the cylinder walls into the gap and lead to a dominant contribution of the vortex in- and outflow to the overall momentum transport. While the contribution of the vortex inflow reaches a nearly fixed value in the ultimate regime of ${\approx}30\%$ at $\unicode[STIX]{x1D707}_{max}$ in the outer gap region, the contribution of the outflow increases monotonically with $Re_{S}$ . Here, we find a strikingly large value of the angular momentum transport of ${\approx}60\%$ in the inner gap region at $Re_{S}=2.14\times 10^{5}$ .

The energy content and azimuthal length scale of these small-scale flow structures in the presence of large-scale rolls are calculated based on azimuthal energy co-spectra. The small-scale plumes show an azimuthal extent of $k_{\unicode[STIX]{x1D711}}d\in [10,20]$ and contribute most to the overall fluctuation energy in regions where the large-scale rolls transport fluid away from the wall, i.e. the so-called ejection regions (Ostilla-Mónico et al. 2014c ). In addition, the energy spectra combined with a correlation analysis revealed that the large-scale Taylor rolls feature an azimuthally oscillating behaviour instead of being axisymmetric and stimulate mid-scale azimuthal structures of sizes $r\unicode[STIX]{x1D711}/d\approx 0.13-0.15$ , where fluid transported by these rolls impact on the cylinder walls.

In order to finally capture the underlying flow structure of non-axisymmetric Taylor rolls, we performed a complex POD analysis at the height of the vortex centre. The large-scale azimuthally travelling waves are superimposed on turbulent Taylor vortices, not only in the classical regime at $\unicode[STIX]{x1D707}=0$ , but also in the ultimate regime at $\unicode[STIX]{x1D707}_{max}$ . These waves propagate into the direction of the mean flow and are similar to the well-known wavy Taylor vortices. While the wave speeds of $0.35\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ in the classical, and $0.11\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ in the ultimate regime are of the same order as the one for the wavy Taylor vortices measured by King et al. (1984), the azimuthal wavenumber found, of $k_{\unicode[STIX]{x1D711},w,1}=9$ , is much higher than in the laminar regime.

Our findings reveal the intrinsic statistical relation between structures of different scales: large-scale Taylor rolls and small-scale plumes in high Reynolds number TC flow. Their interplay strongly relies on the specific locations along the vortex (inflow, centre or outflow) and since they play a prominent role for the flow organization and the momentum transport, our study underlines the importance of an axial exploration when studying the statistics of turbulent TC flows.

We finally note that we believe that our results are much more general than only holding for turbulent TC flow. Clearly, they will generalize to turbulent RB flow with its organization on very large scales (Pandey, Scheel & Schumacher 2018; Stevens et al. 2018), but also to pipe and channel flow with highly organized structures in spanwise direction(Smits et al. 2011; Jimenez 2012; Marusic & Monty 2019). What mechanism, however, sets the length scale of the organization of such superstructures in the spanwise direction remains unclear.