2. Computer experiments

We perform axisymmetric DNS of corotating Taylor–Couette flow (table 1) with periodic BC in $z$ for moderate Reynolds numbers, allowing for both large computational domains and long integration times at affordable computing costs. We integrate the incompressible Navier–Stokes equations (subject to no-slip BC in $r$ ) forward in time ( $t$ ) using our pseudo-spectral DNS code nsCouette (López et al. Reference López, Feldmann, Rampp, Vela-Martín, Shi and Avila2020). The equations are formulated in cylindrical coordinates ( $r,\theta ,z$ ) and rendered dimensionless using $d$ , ${\nu }/{d}$ and ${d^2}/{\nu }$ (i.e. the characteristic viscous time scale of the problem), as unit length, unit speed and unit time, respectively. In axisymmetric DNS, the $\theta$ -dependence is dropped, but all three velocity components are computed.

Table 1. Taylor-roll drift dynamics for $\Gamma = 24$ and $\eta = 0.99$ . Listed are control parameters ( ${R\hspace {-0.1em}a}, {R\hspace {-0.1em}e_{i}}, {R\hspace {-0.1em}e_{o}}, {R\hspace {-0.1em}e_{\hspace {-0.1em}s}}, {R_{\hspace {-0.1em}\Omega }}$ ), response parameters ( $N\hspace {-0.20em}u_{\hspace {-0.1em}s}$ , $N\hspace {-0.20em}u$ , $R\hspace {-0.1em}e_{\tau }$ ), standard deviations $V(\alpha )^{{1}/{2}}$ of the drift speed, $\alpha =v_R$ , the net axial flux, $\alpha =\overline {u}_z$ , and the effective diffusion coefficient, $D_R$ , estimated as in figure 2. NoFX stands for no axial flux ( ${\overline u}_z=0$ ) enforced in the simulation.

Motivated by an exact analogy (Eckhardt et al. Reference Eckhardt, Doering and Whitehead2020) between two-dimensional RBC and axisymmetric TCF in the narrow-gap limit ( $\eta ={r_i}/{r_o}\rightarrow {1}$ ), we set $\eta = 0.99$ and vary ${R\hspace {-0.1em}e_{\hspace {-0.1em}s}} = {Ud}/{\nu }$ , ${R_{\hspace {-0.1em}\Omega }} = {2d\Omega }/{U}$ and $\Gamma = {L_z}/{d}$ . Here, $d=r_o-r_i$ , is the gap width between the inner ( $i$ ) and the outer ( $o$ ) cylinder with radius $r_{i/o}$ , $L_z$ is their axial length, $U = u_{\theta ,i}-\eta \,u_{\theta ,o}$ and $\Omega = {u_{\theta ,o}}/{r_o}$ , where $u_{\theta ,{i}/{o}}$ denotes the azimuthal speed of the cylinders. An important response parameter is the Nusselt number ( $N\hspace {-0.20em}u_{\hspace {-0.1em}s}$ ), which quantifies the transport of angular momentum across the fluid layer (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007; Brauckmann et al. Reference Brauckmann, Salewski and Eckhardt2016; Eckhardt et al. Reference Eckhardt, Doering and Whitehead2020). According to the exact Navier–Stokes mapping of Eckhardt et al. (Reference Eckhardt, Doering and Whitehead2020), flows are mathematically identical as long as ${R\hspace {-0.1em}a}={R\hspace {-0.1em}e_{\hspace {-0.1em}s}}^2\,{R_{\hspace {-0.1em}\Omega }} (1-{R_{\hspace {-0.1em}\Omega }} )= \text {const.}$ , with ${N\hspace {-0.20em}u} = 1+{({N\hspace {-0.20em}u_{\hspace {-0.1em}s}}-1)}/{(1-{R_{\hspace {-0.1em}\Omega }})} $ and ${N\hspace {-0.20em}u}$ being the usual Nusselt number in RBC.

In a first set of DNS (compiled later in figure 3 a), we fix all parameters but $\Gamma$ to investigate the onset of spatio–temporal chaos with respect to the lateral domain size. In a second set (compiled in table 1), we fix $\Gamma$ and explore the effect of shear ( $\textit{Re}_{s\hspace {-0.1em}}$ ) and rotation ( $R_\Omega$ ). The initial conditions are chosen to trigger the desired number of Taylor rolls ( $N_R$ ) necessary to maintain their aspect ratio constant throughout all DNS runs ( ${N_R}/{\Gamma } = 1$ ). This is important because the dynamics is known to depend on ${N_R}/{\Gamma }$ (Ostilla-Mónico et al. Reference Ostilla-Mónico, Verzicco and Lohse2015, Reference Ostilla-Mónico, Lohse and Verzicco2016a ; Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020; Zwirner et al. Reference Zwirner, Tilgner and Shishkina2020). The highest friction Reynolds number ( ${R\hspace {-0.1em}e_{\tau }} = {u_\tau d}/{\nu }$ , where $u_\tau$ is the friction velocity at the cylinder walls) measured in all DNS is 408 (table 1). The spatial resolution in terms of wall units (i.e. based on $Re_{\tau }$ and denoted by ${}^+$ ) is at least $0.07 \leqslant \Delta r^+\leqslant 4.03$ and $\Delta z^+= 4.89$ , which is state of the art in DNS of wall-bounded turbulence (Ostilla-Mónico et al. Reference Ostilla-Mónico, Lohse and Verzicco2016a ,Reference Ostilla-Mónico, Verzicco, Grossmann and Lohse b ; Feldmann et al. Reference Feldmann, Morón and Avila2021).

3. Drift dynamics

Figure 1. Spatio–temporal dynamics of Taylor rolls for different domain sizes ( $\Gamma$ ). Shown are contours of the wall-normal velocity component ( $u_r$ ) at midgap position extracted from DNS ( ${R\hspace {-0.1em}e_{\hspace {-0.1em}s}}= 9475$ and ${R_{\hspace {-0.1em}\Omega }}=0.14$ in all cases). Arrows represent 2000 convective time units. (a) Stationary Taylor rolls in a small (subcritical) domain. (b–d) Axially drifting Taylor rolls in larger (supercritical) domains exhibiting large excursions on a slow time scale. The black line ( $z_R$ ) plotted on top of the $u_r$ space–time data in (d) represents the temporal evolution of the phase angle of the dominant Fourier mode (here $k_z = 12$ , corresponding to 12 pairs of rolls). It serves as a proxy for the collective axial displacement of the entire stack of rolls. (e) Temporal evolution of the modal kinetic energy, $\langle E_{k_{z}}\rangle _r$ , contained in mode number $k_z$ for the same case as in (d); angled brackets denote averaging in the radial direction.

Figure 2. Time series of the axial displacement of Taylor rolls ( $z_R$ ) for different $\Gamma$ at ${R\hspace {-0.1em}e_{\hspace {-0.1em}s}}= 9475$ , ${R_{\hspace {-0.1em}\Omega }}= 0.14$ . (a) Chaotic small-scale oscillations about a fixed axial location in short (subcritical) domains ( $\Gamma \in [4,8]$ ). The simulations span more than 200 viscous time units ( ${d^2}/{\nu }$ ), without reflecting any change of this behaviour; except for the initial transients in the first two viscous time units of the simulations, which we discard for all further analyses in all cases. (b) Huge erratic axial drifts in a larger (supercritical) domain close to the critical point ( $\Gamma = 10$ ). (c) Large erratic axial drifts with qualitatively similar dynamics in larger (supercritical) domains ( $\Gamma \in [16,24,48]$ ). (d) Temporal evolution of the displacement variance, $V(z_R)$ , for the data sets shown in (c). The slope of the corresponding linear fit (broken lines) serves as an estimate for the effective diffusion coefficient for the drift dynamics.

Small domains restrict the dynamics of the system, resulting in nearly stationary rolls. This is apparent from the space–time diagram of the wall-normal velocity ( $u_r$ ) for $\Gamma = 8$ (figure 1 a). If we now enlarge the domain ( $\Gamma \geqslant 10$ , ${N_R}/{\Gamma }=1$ ), the Taylor rolls undergo large, erratic, collective drifts in $z$ that evolve on a slow time scale (figure 1 b–d). In a domain with $N_R = 24$ rolls, for example, the most energetic axial mode is always $k_z = 12$ (figure 1 e), confirming that the space–time representation of $u_r$ is indeed a robust way to identify Taylor rolls and to track their dynamics. Every few viscous time units (e.g. at $t\approx 17.5$ ), the competition with neighbouring modes (here $k_z\in \{11,13\}$ ) represents rare attempts to switch to another state with 11 or 13 pairs of rolls. These attempts, however, remain unsuccessful in all our simulations.

To analyse the drift dynamics quantitatively, we compute axial Fourier spectra of $u_r$ (in particular the $u_r$ midgap space–time data as exemplarily shown in figure 1 d) and use the phase angle information of the dominant axial mode (here $k_z= 12$ ) to approximate the collective displacement of the Taylor rolls ( $z_R$ ), as done earlier by Sacco et al. (Reference Sacco, Verzicco and Ostilla-Mónico2019). The temporal evolution of $z_R$ (figure 1 d) aligns well with $u_r$ , thereby confirming the suitability of $z_R$ to quantify the collective axial drift of the Taylor rolls. For $\Gamma \leqslant 8$ , the rolls first undergo slow transient drifts in the beginning of the simulation and then ultimately oscillate with tiny amplitudes and high frequencies about a statistically steady state (figure 2 a). This fast dynamics of the rolls was reported earlier for three-dimensional turbulent TCF in a domain accommodating one pair of rolls (Sacco et al. Reference Sacco, Verzicco and Ostilla-Mónico2019). By contrast, for $\Gamma = 10$ , the rolls wander more than $ 100d$ before turning back for the first time, and continue moving erratically thereafter (figure 2 b). With further increasing $\Gamma$ , these excursions persist but become less extreme (figure 2 c).

We quantify the Taylor-roll motion statistically by computing the variance of the axial displacement, $V(z_R, t) = \langle z_R^{2}\rangle _{t} - \langle z_R\rangle ^{2}_t$ , where angled brackets denote temporal averaging up to time $t$ . For $\Gamma \leqslant 8$ , the fast dynamics of the rolls is centred around a fixed location and $V(z_R)$ quickly saturates to a constant, in agreement with Sacco et al. (Reference Sacco, Verzicco and Ostilla-Mónico2019), who reported Gaussian fluctuations of $z_R$ with constant variance. By contrast, for $\Gamma \geqslant 10$ , $V(z_R)$ grows approximately linearly with time, as in a Wiener process (figure 2 d). The drift of the rolls can thus be characterised with an effective diffusion coefficient, $D_R$ , as the slope of a linear fit to the variance. For all our analyses, we use time series of at least 200 viscous time units and we generally discard the first $2 {d^2}/{\nu }$ to exclude initial transients.

4. Transition to large erratic drift dynamics

Figure 3. Transition to large, erratic drift dynamics of Taylor rolls. (a) Beyond a critical domain size (grey line, $\Gamma _c= 9.99$ from a power-law fit to the data), the Taylor-roll drift can be characterised by an effective diffusion coefficient, $D_R$ , as shown in figure 2(d); here for $Ra = 10^7$, $Re_s = 9475$ and $R_\Omega = 0.14$ . Note, that $\Gamma$ must take even integer values in order to maintain a unit aspect ratio of the Taylor rolls. (b) For a fixed domain size ( $\Gamma = 24$ ), the Taylor-roll drift starts at a critical Rayleigh number, ${R\hspace {-0.1em}a}_c\approx 6 \times 10^5$ , and becomes more pronounced as Ra increases. The grey region denotes the uncertainty in determining ${R\hspace {-0.1em}a}_c$ ; it spans the interval between the last run with $D_R=0$ and the first run with $D_R\gt 0$ .

Figure 4. Transition from temporal to spatio–temporal chaos for increasing $\Gamma$ . Shown are premultiplied energy spectra, $\varepsilon (\lambda )$ , at ${R\hspace {-0.1em}e_{\hspace {-0.1em}s}}= 9475$ and ${R_{\hspace {-0.1em}\Omega }}= 0.14$ from a subcritical ( $\Gamma = 8$ ) and a supercritical ( $\Gamma = 24$ ) domain. Spectra for other $\Gamma$ look very similar and are thus not shown here. (a) Premultiplied spectra of the modal kinetic energy (as, for example, in figure 1 e) versus axial wavelengths, $\lambda _z={2\pi }/{\kappa _z}$ , where $\kappa _z$ is the axial wavenumber and angled brackets denote temporal averaging. (b) Premultiplied temporal Fourier spectra of the modal kinetic energy (as, for example, in figure 1 e) for the dominant mode (here, for example, for $k_z= 4$ in the case of $\Gamma = 8$ ) versus temporal wavelengths, $\lambda _t={2\pi }/{\omega }$ , where $\omega$ is the angular frequency. FFT means fast Fourier transform.

The dependence of the diffusion coefficient $D_R$ on the domain size $\Gamma$ is shown in figure 3 for a fixed Rayleigh number ( ${R\hspace {-0.1em}a}=10^7$ ). Our data suggest a divergence of $D_R$ near a critical point ( $\Gamma _c= 9.99$ ) followed by a monotonic decrease as $\Gamma$ increases. To examine the nature of this transition, we compare spatial and temporal Fourier spectra from sub- and supercritical domains (figure 4). For $\Gamma = 8$ , the axial spectrum of $u_r$ presents discrete peaks at wavelength $\lambda _z = 2d$ and its harmonics only (figure 4 a). This implies that the flow state consists of four perfectly synchronised copies of one pair of Taylor rolls. In fact, when comparing this state to those obtained for $\Gamma \in \{2,4\}$ , the same Nusselt numbers, $\langle {N\hspace {-0.20em}u_{\hspace {-0.1em}s}}\rangle = 12.82474\pm 0.00002$ , and spectra (not shown here) are recovered. By contrast, the states obtained for $\Gamma \geqslant 10$ exhibit continuous spatial spectra (e.g. for $\Gamma = 24$ in figure 4 a), indicating spatial defects in the roll structure, i.e. there are no identical rolls in the entire stack. A similar transition to spatio–temporal chaos was reported before for axially oscillated (Avila et al. Reference Avila, Marques, Lopez and Meseguer2007) and hydromagnetic (Guseva et al. Reference Guseva, Willis, Hollerbach and Avila2015) Taylor–Couette flow. For both, however, no slow large-scale drift of the roll patterns was reported, possibly due to much shorter simulation times.

The transition to spatio–temporal chaos also alters the temporal spectra (figure 4 b). For $\Gamma \lt \Gamma _c$ , the temporal spectrum is continuous, indicating temporal chaos, and exhibits a peak at $2\times 10^{-3}{d^2}/{\nu }$ (i.e. approximately 20 convective time units) before falling sharply. This peak is associated with the fast, small-displacement dynamics with $D_R=0$ and Gaussian statistics reported by Sacco et al. (Reference Sacco, Verzicco and Ostilla-Mónico2019). For $\Gamma \gt \Gamma _c$ , the temporal spectrum features an additional broad peak at approximately $ 0.1 {d^2}/{\nu }$ , corresponding to the slow drift dynamics characterised by a Wiener process (i.e. $D_R \gt 0$ and linearly increasing variance). The transition to spatio–temporal chaos is also reflected in the mean Nusselt number, but only in the third digit; $\langle {N\hspace {-0.20em}u_{\hspace {-0.1em}s}}\rangle = 12.76\pm 0.04$ for all $\Gamma \geqslant 10$ . We note that while much effort has been dedicated to remove drifts in the analysis of turbulent dynamics of wall-bounded flows (Willis et al. Reference Willis, Cvitanović and Avila2013; Budanur et al. Reference Budanur, Cvitanović, Davidchack and Siminos2015), here the onset of spatio–temporal chaos appears intrinsically linked to the slow, erratic drift dynamics.

5. Dependence of the drift dynamics on the flow configuration

We exploit the analogy between TCF and RBC (Bradshaw Reference Bradshaw1969; Veronis Reference Veronis1970; Prigent et al. Reference Prigent, Dubrulle, Dauchot, Mutabazi, Mutabazi, Wesfreid and Guyon2006; Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007, Reference Eckhardt, Doering and Whitehead2020) to demonstrate that the drift dynamics is found throughout the centrifugally unstable corotating regime. According to the exact Navier–Stokes mapping of Eckhardt et al. (Reference Eckhardt, Doering and Whitehead2020), axisymmetric TCF systems in the narrow-gap limit ( $\eta \rightarrow 1$ ) are exactly identical if ${R\hspace {-0.1em}a} = {R\hspace {-0.1em}e_{\hspace {-0.1em}s}}^{2}\,{R_{\hspace {-0.1em}\Omega }} (1-{R_{\hspace {-0.1em}\Omega }} ) = \text {const.}$ , i.e. the large-scale drift dynamics is identical as well. Indeed, for moderate outer cylinder rotation ( ${R_{\hspace {-0.1em}\Omega }}\in \{0.14, 0.30\}$ ), the drift statistics are similar (table 1, compare rows 3 and 4) and the same is true for the corrected Nusselt number, ${N\hspace {-0.20em}u} = 1 + {({N\hspace {-0.20em}u_{\hspace {-0.1em}s}}-1)}/{(1-{R_{\hspace {-0.1em}\Omega }})}$ . We attribute the small deviations to small, yet finite, curvature effects ( $\eta = 0.99 \lt 1$ ), which are not included in the analogy. For very slow outer cylinder rotation ( ${R_{\hspace {-0.1em}\Omega }} = 0.05$ ), the drift statistics deviate noticeably (table 1, compare rows 3, 4 and 6). This is as expected because the exact analogy breaks down in the limit of a stationary outer cylinder ( ${R_{\hspace {-0.1em}\Omega }}= 0$ ).

Next, we fix ${R_{\hspace {-0.1em}\Omega }}= 0.14$ and $\Gamma = 24$ , and vary $R\hspace {-0.1em}e_{\hspace {-0.1em}s}$ , thereby varying the Rayleigh number (figure 3 b). At low Ra, the rolls are stationary, as observed for low $\Gamma$ in § 3. Similarly, at a critical point ( ${R\hspace {-0.1em}a}_c\approx 6\times 10^5$ ), the rolls begin to drift in the axial direction. As Ra is increased beyond this critical point, the diffusion coefficient of the drift increases. This is contrary to the effect of increasing $\Gamma$ , but consistent with RBC experiments (Xi et al. 2006), where the rate of erratic rotations of the LSC increases tenfold as Ra increases from $10^9$ to $10^{10}$ . A two-dimensional parametric study of the combined effect of $\Gamma$ and $R\hspace {-0.1em}a$ would be interesting, but expensive and thus out of the scope of this work.

The axial drift of the Taylor rolls is associated with a net mass flux in $z$ with mean speed $\overline {u}_z$ . This flux is strongly correlated to the drift speed ( $v_R = \dot{z}_R = {\partial z_R}/{\partial t}$ ) of the rolls (figure 5 a), and raises the question of whether the roll displacement causes the net axial flux or vice versa. The fact that $v_R$ is approximately $500{d}/{U}$ ahead of $\overline {u}_z$ suggests the former (figure 5 b). We probe this hypothesis by enforcing $\overline {u}_z= 0$ , as in laboratory experiments of TCF with end walls. In our DNS with axially periodic BC, we enforce $\overline {u}_z =0$ by imposing an appropriate adverse pressure gradient at each time step. This technique was previously applied to successfully compare axially periodic simulations to lab experiments for Taylor–Couette flow in the counter-rotating regime (Edwards et al. Reference Edwards, Tagg, Dornblaser, Swinney and Tuckerman1991), with radial heating (Ali & Weidman Reference Ali and Weidman1990) and with axially oscillating inner cylinder (Marques & Lopez Reference Marques and Lopez1997; Avila et al. Reference Avila, Marques, Lopez and Meseguer2007).

Figure 5. Taylor-roll dynamics with and without axial flux constraint for otherwise identical parameters ( ${R\hspace {-0.1em}e_{\hspace {-0.1em}s}} = 9475$ , ${R_{\hspace {-0.1em}\Omega }} = 0.14$ , $\Gamma = 24$ ). (a) Time series of the drift speed of the Taylor rolls ( $v_R = \dot{z}_R = {\partial z_R}/{\partial t}$ ) and the net axial flux ( $\overline {u}_z$ ) for the case in figure 1(d). (b) Close-up to the data in (a). As a visual reference, the black dash represents 500 convective time units. (c) Time series from a simulation with no axial flux (NoFX, i. e. $\overline {u}_z =0$ ) enforced.

As a result of suppressing the axial mass flux, $v_R$ , $V(v_R)$ and $D_R$ are substantially reduced (figure 5 c, table 1, compare rows 4 and 5), but when rescaled, the drift dynamics remains qualitatively unaltered (figure 6 a). Specifically, $V(z_R)$ still increases linearly with time (figure 6 e), although at a slower pace.