Statistics, plumes and azimuthally traveling waves in ultimate Taylor-Couette turbulent vortices

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 $\text{Re}_S=3\times 10^5$}. The velocity field in the gap confined by coaxial and independently rotating cylinders at a radius ratio of $\eta=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\'{e}nard convection, which emphasizes the analogies of these both systems. By performing a {complex proper orthogonal decomposition}, we remarkably detect azimuthally traveling 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 center, is similar to the well-known wavy Taylor vortex flow{,} which has comparable wave speeds, but much larger azimuthal wave numbers.


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 et al. 1986;Donnelly 1991). More recent reviews on the hydrodynamic instabilities can be found in Fardin et al. (2014), and on fully turbulent Taylor-Couette flows in Grossmann et al. (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 η = r 1 /r 2 , and the aspect ratio Γ = ℓ/d, with ℓ 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) Re S = 2r 1 r 2 d (r 1 + r 2 )ν |ω 2 − ω 1 | = u S d ν , (1.1) with the cylinder angular velocities ω 1,2 , ν 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 implying µ > 0 for co-rotation of the cylinders, µ = 0 for pure inner cylinder rotation and µ < 0 for counter-rotation. The most important response parameter of the TC system to the cylinder driving is the angular velocity transport (Eckhardt et al. 2007a) J ω = r 3 ( u r ω ϕ,z,t − ν∂ r ω ϕ,z,t ) , (1.3) where r denotes the radial coordinate, ϕ the azimuthal coordinate, t the time coordinate and · ϕ,z,t the azimuthal-axial-time average. This quantity is conserved along r and can be directly measured by the torque T acting on either the inner (IC) or the outer cylinder wall (OC). Normalizing J ω with its corresponding laminar non-vortical value J lam ω = 2νr 2 1 r 2 2 (ω 1 − ω 2 )/(r 2 2 − r 2 1 ) yields a quasi Nusselt number Nu ω = J ω /J lam ω (Eckhardt et al. 2007a), which is analogous to the Nusselt number Nu in Rayleigh Bénard flow (RB) flow, i.e. the buoyancy-driven flow which is heated from below and cooled from above. There, Nu is a measure for the amount of transported heat flux normalized by the purely conductive heat transfer. Eckhardt et al. (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 Nu ω ∼ Re α S , and on the rotation ratio µ has been widely investigated (Lathrop et al. 1992;Lewis & Swinney 1999;van Gils et al. 2011;Paoletti & Lathrop 2011;Merbold et al. 2013;Ostilla-Mónico et al. 2014c;Brauckmann et al. 2016a;Grossmann et al. 2016). For pure inner cylinder rotation (µ = 0), a change in the local scaling exponent α 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 η and is located around Re S,crit ≈ 1.6 × 10 4 for η = 0.714. In the ultimate regime, the scaling exponent becomes α ≈ 0.76 independent on η (Ostilla et al. 2013). When the driving (Re S ) is kept constant and only µ is changed, the TC flow features a maximum in the angular momentum transport (Nu ω ). For η = 0.714, the maximum is located in the counter-rotating regime at µ max ≈ −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 et al. (2017) by means of the decomposition of Nu ω 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 PIV measurements in a TC facility at η = 0.917 and Re S ∈ [1.1 × 10 4 ; 2.9 × 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 µ = 0 at η = 0.875 and for different aspect ratios. They find a regime of reappearance for 28 Re/Re C 36 and for Γ 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 η = 0.904 and Γ = 20. The azimuthal waves are identified based on Fourier analysis and proper orthogonal decomposition (POD) in the range of 23 Re/Re C 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 η = 0.733 and Γ = 34. They capture the reappearance of azimuthal waves for 20 Re/Re C 38 based on spatial correlations. Note that the three aforementioned studies are all based on pure inner cylinder rotation (µ = 0). More recently, Merbold et al. (2014) performed flow visualizations in TC flow with η = 0.5 at Re S = 5000. They find an axial oscillation of the turbulent Taylor vortices in the range of µ ∈ [−0.15, −0.3], which includes the rotation ratio for optimum transport µ 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 µ = 0 and η = 0.714 the radial profiles of the angular velocity ω at specific axial positions of the large-scale Taylor vortices; namely at the vortex inflow, vortex center 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 1a). 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 1c). In between the in-and outflow, the radial velocity component becomes zero in the middle of the gap, which denotes the location of the center of the vortex, as shown in figure 1b. 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 ω have a logarithmic shape in the corresponding boundary layer (Ostilla-Mónico et al. 2016). It is worth to mention that in the ultimate regime, the profiles become logarithmic also in the absence of dominant large-scale rolls (Huisman et al. 2013b). For η = 0.5 in the classical regime, van der Veen et al. (2016b) calculated the velocity of such plumes based on planar particle image velocimetry (PIV) measurements performed at different heights, which further confirms the connection between small-scale plumes and large-scale vortices.
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 et al. (2016b). As a reference point for the comparison, they choose Re S = 2 × 10 4 for TC and a Rayleigh number of Ra = α P gH 3 ∆T /(κν) = 10 7 for RB, where the Nusselt number is identical for both flows. Here, Ra is the dimensionless driving parameter in RB flow with α p the thermal expansion coefficient, g the gravitational acceleration, H the height of the RB cell, ∆T the temperature difference, and κ 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 volumeaveraged Nusselt number, along with a broadening of the tails with increasing Ra. Shang et al. (2004) revealed that the asymmetry of the PDFs of 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 Nu ω were analyzed for µ = 0 by Huisman et al. (2012) without any connection to specific flow structures.
Another approach investigating structures in TC flow is to analyze the kinetic energy spectra. Dong (2007) performed direct numerical simulations (DNSs) for η = 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 η = 0.909 and µ = 0 at Re S 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 et al. 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 center and vortex outflow regions for the momentum transport, how do small-scale structures affect this transport, what is the lengthscale of these structures and do Wavy-vortexlike 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 η = 0.714, in the range of Re S ∈ [9.3 × 10 3 , 3.5 × 10 5 ] and µ ∈ [0, −0.36]. We want to stress that such quasi-threedimensional experimental investigations of the local angular momentum transport statistics and flow structures in the ultimate turbulent TC flow at µ max are unique, while the measurement setup 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 setup, 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 largescale turbulent Taylor vortices onto the global flow statistics ( § 3). In § 4 the global and local angular momentum transport are analyzed 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, center, and outflow in § 5. The energy content and azimuthal lengthscale of these structures are calculated in § 6 based on azimuthal energy co-spectra, while azimuthally traveling waves superimposed to 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).

Experimental setup
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 setup are r 1 = 75 mm and r 2 = 105 mm, respectively and thus the radial gap is d = r 2 − r 1 = 30 mm. The height of the cylinders is ℓ = 549 mm, which gives an aspect ratio of Γ = ℓ/d = 18.3. The radius ratio is then η = r 1 /r 2 = 0.714. A more detailed overview of the setup can be found in Huisman et al. (2015). The flow consists of water whose viscosity ν and density ρ can be controlled throughout the experiments due to the temperature control of the BTTC. We fix the temperature of the experiments to 20 • C, leading to ν = 1.002 mm 2 /s and ρ = 0.998 g/cm 3 . The standard deviation of the temperature is 15 mK. The flow is seeded with fluorescent polyamide particles with diameters up to 20 µm with an average particle density of approximately 0.01 particles/pixel. These particles are coated with Rhodamine B which has a maximum emission centered at around 565 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 ≈ 1 mm. The laser is mounted to 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 used for the recordings is an Imager sCMOS (2560 × 2160 px) 16 bit with a Carl Zeiss Milvus 2.0/100 lens. We capture the velocity fields with a framerate of f = 15 Hz. Since the camera is operated in double frame mode, we can have very small interframe times, i.e. ∆t ≪ 1/f . In order to maximize the contrast of the images, we use a long-pass filter in front of the lens (Edmund High- The velocity fields are calculated with commercial software (Davis 8.0) using a multipass method. The algorithm uses windows of size 64×64 px for the first pass and windows of 24 × 24 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 ∆r and azimuthal ∆ϕ resolution is the same as the spatial resolution in Cartesian coordinates ∆x, i.e. ∆r = ∆x and r∆ϕ = ∆x. In this way, the resultant velocity fields are of the form u = u r (r, ϕ, t)e r + u ϕ (r, ϕ, t)e ϕ , where u r and u ϕ are the radial and azimuthal velocity components which depend on the radial coordinate r, the azimuthal coordinate ϕ and the time coordinate t. e r and e ϕ 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, 8 cases were investigated which will be addressed in the following sections. Each case contains measurements done at 23 different heights which are separated by ∆z = 4 mm, and each height contains 1500 different velocity fields. Thus, the height of the experiments spans a length of ∆ℓ = 22∆z = 88 mm. The resolution of the velocity fields ∆x depends on the height but lies within ∆x ∈ [0.607, 0.752] mm, where the smallest value corresponds to the height closest to the camera at (z − ℓ/2)/d = 1.5 and the smallest to (z − ℓ/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 α by Ostilla-Mónico et al. (2014b) at Re S (η = 0.714) ≈ 1.6 × 10 4 for µ = 0, which is caused by a transition of the boundary layers (BLs). Accordingly, flow states at Re S < 1.6 × 10 4 are assumed to be in the so-called classical regime with

Case Regime
ReS µ Abbreviation C # | µ Re S Linestyle 1 classical 9.32 × 10 3 0 C1| 0 9.32×10 3 ---  laminar BLs, while those at Re S > 1.6 × 10 4 are assumed to be in the ultimate regime with turbulent BLs. Case 1 and 2 are in the classical regime, where, µ max changes with the shear Reynolds number, which is why µ 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 µ max = −0.36.

Flow states and velocity profiles
The TC flow at µ max is dominated by turbulent Taylor vortices, while at µ = 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 asr = (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 u ϕ = (u ϕ − u ϕ,2 )/(u ϕ,1 − u ϕ,2 ), with u ϕ = rω. In order to normalize the radial velocity component, the shear velocity is usedũ r = u r /u S .
For both flows in the classical regime (see figure 3a,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 µ max = −0.15, indicating an increase in strength. In the ultimate regime at Re S = 2.15 × 10 5 and µ = 0 (C 6 ), 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, figure  6 in Ostilla-Mónico et al. (2014c). When the rotation rate is changed to µ 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 3d is shown in   figure 3e. This representation exemplifies the further analysis. The recorded 23 heights capture more than one vortex pair, which is why we exclude the grayly marked data points for flow quantities calculated over the axial coordinate. The axial length of evaluation therefore starts and ends at a vortex center. In between, we use the minimum of the azimuthally-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 center. In figure 4, we show the normalized radial profiles of the azimuthal and radial velocity components averaged over space (cylindrical surfaces) and time ( · ϕ,z,t ). The slopes of the profiles in the bulk nearly vanish at µ = µ 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. 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 center of the gap and increasing values close to the wall. This increase is more pronounced at µ max ; while for the lowest shear Reynolds number, a peak close to the outer cylinder wall is visible. In case of the radial velocity component, the standard deviation depicts a maximum in the center of the gap and decreases towards the cylinder walls. The maximum of σ ϕ,z,t (ũ 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 midheight.

Probability density function (PDF) 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 atr = 0.5 and Re S = 2.1 × 10 5 for µ = 0 (C 6 ) and µ = −0.36 (C 7 ) as representatives. In addition, the underlying PDFs at the locations of the vortex inflow, vortex center and vortex outflow are included.
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 LDV at midheight andr = 0.5. When all data at different heights are included in the PDFs, they still follow a Gaussian shape at µ = 0 for both velocity components. At µ max however, the PDF of the azimuthal velocity depicts a cusp-like form centered at the origin and approximately exponential tails in accordance with the numerical simulations at a lower Reynolds number of Re S = 2 × 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 on 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 behavior to the flow organization in RB flow. In the case of the radial velocity at µ max , the PDF also deviates from the ideal Gaussian shape, which can be attributed again to the intermittent bursting plumes.
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 center into the direction of the inner cylinder wall for µ 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 atr = 0.1, both tails become increasingly exponential. This behavior 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.
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 µ max the skewness increases from the inner cylinder wall to a maximum aroundr = 0.16, then decreases to negative values with a minimum aroundr = 0.86 for Re S = 2.14 × 10 5 (C 7 ) and than increases again in the direction of the outer cylinder. Furthermore, the absolute skewness increases with the shear Reynolds number. This behavior is similar to the one of the temperature fields in RB flows reported by Emran & Schumacher (2008). The kurtosis profiles at µ max depict two maxima, one in the inner gap region atr = 0.22 and one in the outer gap region atr = 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 µ max to higher forcings.

Angular momentum transport
As the angular momentum transport, expressed in terms of the quasi-Nusselt number Nu ω , 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 analyze the global angular momentum transport to compare its amount with the literature, and second, we analyze the axial dependent radial profiles of Nu ω 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 vortexdominated momentum transport. The Nusselt number is composed of a convective and a viscous part (Eckhardt et al. 2007a): While the viscous term Nu ν ω dominates in the boundary layers, the convective term Nu c ω dominates in the bulk. Since the focus of our investigation is set to the bulk region, we neglect thus the viscous part. Furthermore, as it is shown in figure 4b, the radial velocity component nearly vanishes when averaged over cylindrical surfaces with an axial length of one vortex pair: u r ϕ,z,t ≈ 0. Therefore, r 2 u r u ϕ ϕ,z,t ≈ r 2 u ′ r u ′ ϕ ϕ,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 To avoid confusion, in the following, we will call the Nusselt number Nu ω = Nu c ω +Nu ν ω shown in equation 4.1, the total Nusselt number. In figure 9a, we show the radial profiles of Nu ω indicated by lines, and the net convective momentum flux Nu c,net ω , indicated by diamonds (µ max ) and circles (µ = 0). The analysis is restricted to the bulk region withr ∈ [0.1, 0.9]. A slight dependence of the Nu ω -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 setup in contrast to their theory, which is based on cylinders of infinite length. However, the effect of the cylinder length is reduced in our setup by cutting the axial length of evaluation to the size of one vortex pair. Moreover, also the vortex aspect ratio influences the value of Nu ω 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 Nu ω -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 . Furthermore, in figure 9b, 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) (µ = 0, η = 0.724) and DNS of Ostilla-Mónico et al. (2014c) (µ = −0.4, η = 0.714). We find a very good agreement with these data, which enables a more detailed analysis of the momentum transport.
In order to get deeper insight into the relation between the convective and net convective Nusselt number, we plot in figures 10a,b both local quantities in a meridional plane (r − z plane) for Re S = 2.14 × 10 4 and µ = −0.36 (C 7 ). 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 Nu c ω (figure 10a), the momentum transport in the area of the vortex outflow is positive and strongly concentrated to 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 in average a positive net transport. Alternatively, when the net transport is plotted (figure 10b), 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 10c we show the joint PDF of the radial and azimuthal velocity fluctuations atr = 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 Nu c ω and Nu 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. 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 (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 µ max . This reflects the strong correlation of the enhanced momentum flux at µ max and the turbulent Taylor vortices. Further, in the ultimate regime at µ = 0, again the effect of remnants of these large vortices becomes visible in the slight dependence of the Nusselt number onr. When the ultimate regime is reached at Re S = 2.98 × 10 4 at µ max (C 3 ), the contributions of the vortex in-and outflow are small and become significant again at Re S = 6.68 × 10 4 (C 4 ). 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.

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 analyze now the PDFs of Nu 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 with the individual standard deviations σ x and σ y and the correlation coefficient ρ P = xy /(σ x σ y ). Equation (5.1) represents an inclined elliptic shape for the joint PDF in contrast to the one shown in figure 10c. Furthermore, the PDF of the product z = xy is ( with K 0 the modified Bessel function of the second kind. The prefactor of the exponential function in equation (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 ρ P cannot be normalized in a proper way, which leads to different slopes of the exponential tails of the prediction depending on ρ P . Therefore, in our calculation, we use the prediction according to equation (5.2) for the case Re S = 3.51×10 5 and µ = 0 (C 8 ).
In figure 12, we show the PDFs of Nu c,net ω , which are normalized with the factor σ ur σ u θ 1 − ρ 2 P in order to compare their shapes, atr = 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 µ max the PDFs become highly skewed due to a change of shape in the positive tails. This deviation results in an increasing but still rare number of positive extreme events of momentum flux, which becomes more pronounced with increasing Re S . According to figure 12c, 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 10c). 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 µ 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 × 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 µ 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 Nu ω -PDFs.
In figure 13, we show the PDFs of Nu c,net ω for various different radii r for the case of Re S = 2.14 × 10 5 at µ = −0.36 (C 7 ). The negative tails of the PDFs are largely identical for all investigated radial positions, while the right tails strongly depend onr. 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 ϕ -PDFs in figure 7. As was shown in figure 11b, the maximum of the outflow contribution to the overall momentum transport is close tor = 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 13b). However, in the outer gap region at 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 Nu 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 analyze the spatial energy co-spectra to detect the presence and lengthscale 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 ϕ-direction with a total number of azimuthal points of n ϕ = 0, 1, ..., N − 1 equidistantly spaced by ∆s = ∆ϕr c i.e.
The discrete spatial Fourier transform U r,ϕ of both fluctuation components u * r,ϕ is given by 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ϕ can be calculated as with the wavenumber vector k nϕ ϕ = (∆s) −1 n ϕ /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 its graph First of all, we show the temporally and axially averaged azimuthal energy co-spectra atr = 0.5 in figure 14a 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 ϕ d ≈ 20. It is worth mentioning that we do not see any peak in the large-scale regime, Figure 15. Temporally averaged pre-multiplied azimuthal kinetic energy co-spectra in the classical regime for ReS = 9.30×10 3 and µ = −0.15 (C2) at different radial positions at the axial height of (a) vortex inflow, (b) vortex center and (c) vortex outflow. The region of kϕd ∈ [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.
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 ϕ = 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 k ϕ d 20. When Re S is increased into the ultimate regime (with µ fixed at µ = µ max ), this kink continuously diminishes and energy is redistributed from the large to the small scales. Based on the power-law ansatz E r,ϕ ∼ k γ ϕ , the local scaling exponent γ corresponding to the co-spectra is shown in figure 14b. γ is calculated with a bin size of log 10 (k ϕ ) = 0.5. It is worth mentioning that γ is sensitive to the data processing and the accompanying error propagation. γ is given here for this choice of bin size and to show the trend of γ 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 µ fixed at µ = µ max ), the exponent decreases with increasing k ϕ before it strongly drops down in the viscous regime beyond k ϕ d ≈ 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 × 10 5 and µ = 0 (C 8 ), the local exponent is nearly constant for k ϕ d 10 with a value of γ ≈ −1.52, which is slightly above −5/3.
Next, we focus on the local pre-multiplied energy co-spectra at the axial height of the vortex inflow, vortex center, and vortex outflow respectively, where the PDF analysis of Nu c,net ω suggests the occurrence of small-scale intermittent plumes. Pre-multiplied means that the co-spectra are multiplied with the wavenumber vector k ϕ , such that the area under its graph corresponds to the kinetic energy (Smits et al. 2011). In figure 15, we show the pre-multiplied energy co-spectra in the classical regime for Re S = 9.32 × 10 3 and µ = −0.15 (C 2 ) at the three vortex positions. Firstly, we address the peak that corresponds to the large scales. For the vortex inflow in figure 15a, the large-scale peak is located around k ϕ d ≈ 1.1 within the bulk for 0.2 r 0.8, while close to the outer cylinder atr = 0.9, the peak lies outside of our resolvable scales. Close to the IC wall, the peak shifts to k ϕ d ≈ 2.2 atr = 0.1. In figure 15c, at the location of the vortex outflow, the opposite behavior 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 center in figure 15b, the large-scale peak is shifted from smaller wavenumbers in the center of the gap (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 ϕ d ≈ 2.2 are formed. Since the axisymmetric flow has a wavenumber of k ϕ = 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-scales peak, we observe that it is located at wavenumbers around k ϕ d ∈ [10, 20] for all three depicted heights, independent of the radial coordinate. However, its amplitude strongly varies withr and the height z. In the region of the vortex inflow shown in figure 15a, the small-scale peak is most pronounced near the OC and its amplitude decreases monotonically with decreasingr. On the contrary, at the height of the vortex outflow in figure 15c, the small-scale peak amplitude is largest close to the IC and decreases in amplitude towards the OC. In the region of the vortex center in figure  15c, we find the highest amplitude of the small-scale peak in the center 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 center -where no predominant flow direction concerning the radial velocity component is present-the behavior is different: Plumes rise from both cylinder walls and travel towards the bulk flow which results in the highest peak amplitude atr = 0.5. Considering the fact, that the contribution of the vortex center to extreme and strong events of momentum flux is almost negligible (see figures 12, 13), these detected plumes seem to compensate there total radial momentum transport.
The pre-multiplied energy co-spectra in the ultimate regime at Re S = 2.1 × 10 5 for the three vortex locations are depicted in figure 16 for µ = 0 (C 6 ,a-c) and µ max (C 7 ,d-f). For µ = 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 3c, 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 onr. However, the co-spectra show a prominent change in the slope around k ϕ d ≈ 20: In figure 16a (vortex inflow), a small-scale peak is formed around k ϕ d ∈ [10, 20] close to the OC atr = 0.9. This is comparable to the previous case in the classical regime. Also in the region of the vortex outflow (figure 16c), we observe a peak within the same range of scales atr = 0.1. For the vortex center however (figure 16b), no peak is visible.
Also when µ is changed to µ max (see figures 16d-e), we identify a similar behavior 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 16d, a large-scale peak is present which shifts to k ϕ d ≈ 3.9 atr = 0.1. At the vortex outflow in figure 16f, this shift appears close to the OC wall atr = 0.9 for the same wavenumber. In the region of the vortex center (figure 16e), the large-scale peak is most pronounced atr = 0.5 at 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 whenr is varied from 0.1 to 0.9, as it is clearly visible in the marked regime of k ϕ d ∈ [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 lengthscale of these plumes is in the order of k ϕ d ∈ [10, 20].

Spatial correlation coefficients
As it 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 center of the gap are worked out in more detail. The spatial two-point autocorrelation coefficient of the velocity fluctuations between two points separated by ∆ϕ in the azimuthal direction is given by R rr (r, z, ∆ϕ) = u * r (r, ϕ, z, t) u * r (r, ϕ + ∆ϕ, z, t) ϕ,t u * 2 r (r, ϕ, z, t) ϕ,t , (6.5) (6.6) In figure 17a, we show the spatial two-point autocorrelation function of u r atr = 0.1 in the ultimate regime for three different flow cases at µ 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ϕ/d ≈ 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 17b), the opposite behavior is observed with a pronounced minimum of the autocorrelation function for the vortex outflow at rϕ/d ≈ 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 figures 17c,d, we show both the azimuthal and radial two-point autocorrelation function for the center of the vortex evaluated atr = 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 (C 6 ). These figures reveal that in all vortex dominated cases, a large-scale oscillating behavior 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.

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, 16 and 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 et al. (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 π/2, based on the algorithm described by Marple (1999). We further use a combined analysis of both velocity components u r and u ϕ for the POD analysis. The applied algorithm is based on three steps . At first we use the Hilbert transform H, to make the velocity fluctuations complex: Secondly, the imaginary, transformed velocity fluctuations (u * r,c , u * ϕ,c ) are arranged in a data matrix D, where columns are assigned to the spatial grid points (1 : M ) and rows to the time signals (1 : N ), while M > N : At the end, a singular value decomposition (SVD) of the data matrix is performed as D = ΦΣΨ T . The columns of the matrix Φ contain the left singular vectors of D and therefore the N complex modes CPOD(r, ϕ). The diagonal elements of Σ 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 The time-dependent coefficients a i (t) result from a projection of the complex modes onto the data matrix, i.e. a i = D T CPOD i . In the following, we will present the CPOD results for two flow states at the height of the vortex center, where the azimuthal two-pointcorrelation showed a large-scale oscillating behavior. The first case is in the classical regime at µ = 0 and the second one in the ultimate regime at µ max .

CPOD in the classical regime
The first CPOD mode for Re S = 9.32 × 10 3 and µ = 0 (C 1 ) is depicted in figure 18 for both velocity components. With the first temporal coefficient a 1 (t) and the first mode CPOD 1 , the corresponding velocity fluctuation field is reconstruct by u * 1 (r, ϕ, t) = a 1 (t)CPOD 1 (r, ϕ). In figure 18a, the real part of 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 18c, indicating an azimuthally traveling wave. The real part of the azimuthal velocity of CPOD 1 in figure 18b, 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 18d 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.
In figure 19a, we show the energy fraction captured by the CPOD modes for Re S = 9.32 × 10 3 and µ = 0 (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 19b, together with the temporal energy spectrum of the full field radial velocity component, calculated atr = 0.5 and averaged over ϕ. The dominant frequency of CPOD 1 is f /(f 1 − f 2 ) = 3.03, identical to the one of the radial velocity component. Thus we can assume that the first mode captures the temporal behavior of the full field quite accurately. The frequencies f 1 and f 2 represent the rotation frequencies of the inner and outer cylinder, respectively.
To illustrate the spatiotemporal character of the first mode, we depict in figure 20a its reconstructed radial velocity component as well as the corresponding full field at r = 0.5 as a function of the azimuthal coordinate ϕ and time t. The space-time plot of 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 η = 0.733 has already been found by Wang et al. (2005) in the Reynolds number regime of 20 Re/Re C 38. In our study in the classical regime for µ = 0 at η = 0.714, we find Re S /Re S,C ≈ 99, where the critical shear Reynolds number is located around Re S,C ≈ 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 λ w and phase speed c w of CPOD 1 . The temporal phase function φ i is defined as φ i (t) = arctan ℑ(a i (t)) ℜ(a i (t)) . (7.5) Its temporal derivative is equal to the angular frequency of the CPOD mode (Suanto et al. 1998), leading to a wave frequency of f w,1 = ∂ t φ 1 (t)/(2π) = 1.80 Hz. In addition, the spatial phase function Φ i (x) is defined as ϕ)) , (7.6) whose azimuthal derivative is a measure for the local wavenumber k ϕ,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 r = 0.5, is given by λ w,1 = 2π/(∂ ϕ Φ 1 (r = 0.5, ϕ)) = 0.68 rad, leading to an azimuthal wavenumber of around k ϕ,w,1 = 2π/λ w,1 ≈ 9. In this way, the phase speed becomes c w,1 = λ w,1 f w,1 = 1.22 rad/s = 0.35ω 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 η = 0.73, R/R c ≈ 14, and an average axial wavelength ofλ w /d = 2.4 a wave speed of c w ≈ 0.2ω 1 , while the corresponding pattern features k ϕ,w = 2. Additionally, they show for a slightly larger radius ratio of η = 0.84, that the wave speed increases when Re/Re c > 18. Thus, we can conclude that the wave speeds superimposed to the Taylor vortices in the wavy Taylor vortex regime and in the turbulent Taylor vortex regime are in 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).

CPOD in the ultimate regime
We now turn our attention to the ultimate regime at the torque maximum rotation rate. In figure 21a-d, the first CPOD mode, representing the real and imaginary part of the radial and azimuthal velocity component, for Re S = 6.68×10 4 and µ = −0.36 (C 4 ) at the axial height of the vortex center is depicted. The flow in the ultimate regime features nearly the same pattern as already shown in 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 ℜ(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 equations 7.5 and 7.6), the azimuthally traveling wave can be described by a wave frequency of f w,1 = 3.67 Hz, a wavelength of λ w,1 = 0.74 rad and a phase speed of c w,1 = 2.72 rad/s = 0.11∆ω. Remarkably, the turbulent Taylor vortices feature azimuthal waves also in the ultimate turbulent regime at µ max . To the best of our knowledge, this is η ReS µ kϕ,w,1 cw,1 [rad/s] Classical regime (C1) 0.714 9.32 × 10 3 0 ≈ 9 0.35∆ω Ultimate regime (C4) 0.714 6.68 × 10 4 −0.36 ≈ 9 0.11∆ω King et al. (1984) 0.730 1.05 × 10 3 0 2 0.21∆ω Table 2. Overview of detected azimuthally traveling waves in the classical and ultimate regime in comparison to the findings of King et al. (1984).
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 × 10 5 and µ = 0 (C 6 ), where the large Taylor rolls disappear in the mean field, does not yield any sign of traveling waves. As a comparison tool, we show in table 2 the results concerning the azimuthally traveling waves found in this study and the findings of King et al. (1984) for η = 0.73.
To what extent the detection of azimuthal traveling waves is dependent on the turbulence level (Re S ) and the geometry of the system (η, Γ ) 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.

Summary & 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 µ 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 µ = 0 and µ < 0, they fade out when the ultimate regime is reached for µ = 0, but can be formed again when sufficient counter-rotation (µ < 0) is introduced, in particular at µ 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 µ = 0 for both velocity components but changes drastically for µ max . There, the PDFs of u ϕ 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 Nu c,net ω fluctuate around zero and not the averaged value of the Nusselt number and can be described (for µ = 0) by the shape of a distribution of the product of two Gaussian variables. At µ 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 around ≈ 30% at µ 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 ≈ 60% in the inner gap region at Re S = 2.14 × 10 5 .
The energy content and azimuthal lengthscale 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 ϕ d ∈ [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 behavior instead of being axisymmetric and stimulate mid-scale azimuthal structures of sizes rϕ/d ≈ 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 center. The turbulent Taylor vortices are superimposed by large-scale azimuthally traveling waves not only in the classical regime at µ = 0, but also in the ultimate regime at µ 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∆ω in the classical, and 0.11∆ω in the ultimate regime are in the same order as the one for the wavy Taylor vortices measured by King et al. (1984), the azimuthal wave number found of k ϕ,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, center 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 (Stevens et al. 2018;Pandey 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 lengthscale of the organization of such superstructures in spanwise direction remains unclear.