Centrifugal instability of Taylor–Couette flow in stratified and diffusive fluids

Abstract

The linear and nonlinear dynamics of centrifugal instability in Taylor–Couette flow are investigated when fluids are stably stratified and highly diffusive. One-dimensional local linear stability analysis (LSA) of cylindrical Couette flow confirms that the stabilising role of stratification in centrifugal instability is suppressed by strong thermal diffusion (i.e. low Prandtl number $Pr$). For $Pr\ll 1$, it is verified that the instability dependence on thermal diffusion and stratification with the non-dimensional Brunt–Väisälä frequency $N$ can be prescribed by a single rescaled parameter $P_{N}=N^{2}Pr$. From direct numerical simulation (DNS), various nonlinear features such as axisymmetric Taylor vortices at saturation, secondary instability leading to non-axisymmetric patterns or transition to chaotic states are investigated for various values of $Pr\leqslant 1$ and Reynolds number $Re_{i}$. Two-dimensional bi-global LSA of axisymmetric Taylor vortices, which appear as primary centrifugal instability saturates nonlinearly, is also performed to find the secondary critical Reynolds number $Re_{i,2}$ at which the Taylor vortices become unstable by non-axisymmetric perturbation. The bi-global LSA reveals that $Re_{i,2}$ increases (i.e. the onset of secondary instability is delayed) in the range $10^{-3}\lt Pr\lt 1$ at $N=1$ or as $N$ increases at $Pr=0.01$. Secondary instability leading to highly non-axisymmetric or irregular chaotic patterns is further investigated by three-dimensional DNS. The Nusselt number $Nu$ is also computed from the torque at the inner cylinder for various $Pr$ and $Re_{i}$ at $N=1$ to describe how the angular momentum transfer increases with $Re_{i}$ and how $Nu$ varies differently for saturated and chaotic states.

1. Introduction

Thermal diffusion in fluid flows is characterised by the (thermal) Prandtl number $Pr=\nu _{0}/\kappa _{0}$ , which denotes the ratio between kinematic viscosity $\nu _{0}$ and thermal diffusivity $\kappa _{0}$ . The Prandtl number $Pr$ varies from fluid to fluid; for instance, we typically consider $\textit {Pr}\sim O(10^{3})$ for glycerol, $\textit {Pr}\simeq 7$ for water, $\textit {Pr}\simeq 0.7$ for air, $\textit {Pr}\sim O(10^{-1})$ in the Earth’s liquid outer core, $\textit {Pr}\sim O(10^{-2})$ for liquid metals, $\textit {Pr}\sim O(10^{-6})$ or less in the interior of the Sun and stars (Calkins et al. 2012; Horn et al. 2013; Garaud 2020a ; Legaspi & Waite 2020). The Prandtl-number dependence has been examined for various flows coupled with heat transfer. In convection, Kerr & Herring (2000) studied a thermal convection problem with no-slip vertical boundaries and free-slip lateral boundaries and they proposed different scaling laws for the Nusselt number $Nu$ versus the Rayleigh number $Ra$ transitions as $Nu\sim Ra^{2/7}$ for $Pr\sim O(1)$ and $Nu\sim Ra^{1/4}$ for low $Pr\ll 1$ . Miquel et al. (2020) investigated the convection in a different configuration with internal sources and sinks and they proposed a different scaling law as $Nu\sim Ra^{1/2}Pr^{\chi }$ , where the exponent $\chi$ transitions from $\chi \simeq 1/2$ for $Pr\leqslant 0.04$ to $\chi \simeq 1/6$ for $Pr\gt 0.04$ . The Prandtl-number dependence has also been explored in studies of stratified flows. For instance, it is found that turbulence in stably stratified fluids permits high-wavenumber temperature fluctuation that leads to the appearance of small-scale turbulence structures as $Pr$ increases when $Pr\gt 1$ (Legaspi & Waite 2020). The effect of thermal diffusion is examined in a study of exact coherent structures in stratified plane Couette flow in the limits of either $Pr\rightarrow 0$ or $Pr\rightarrow \infty$ , the former in which the thermal diffusion is dominant over stratification and the density variation away from a linear profile of the stratification vanishes and the latter in which the density can be mixed and homogenised by advection and its stratification profile deviates from a linear profile (Langham et al. 2020).

Fluid flows with stratification and strong thermal diffusion with $Pr\ll 1$ have been the subject of interest in geophysics and astrophysics due to their relevance to the interior of the Earth and stars including the Sun (Lignières 1999; Calkins et al. 2012; Garaud 2020a ). Figure 1 illustrates a configuration in which a fluid parcel is hypothetically displaced upwards over a length scale $L_{0}$ in a stably stratified fluid with the temperature gradient along the vertical direction of gravity. The temperature of the parcel is lower than the surrounding temperature so the parcel absorbs heat at a rate with a characteristic diffusion time scale $\tau _{{diff}}=L_{0}^{2}/\kappa _{0}$ . Simultaneously, the stratification acts as a restoring force on the parcel leading to internal oscillation with a period $\tau _{{int}}=1/N_{0}$ , where $N_{0}$ is the dimensional Brunt–Väisälä frequency derived from the stratification. On the one hand, if the thermal diffusivity is low as $\kappa _{0}\ll N_{0}L_{0}^{2}$ (i.e. the diffusion time scale is much larger than the internal oscillation period as $\tau _{{diff}}\gg \tau _{{int}}$ ), the restoring force is dominant and the fluid parcel descends to the original position due to the gravitational force, a well-known mechanism for the generation of internal gravity waves. On the other hand, if the thermal diffusivity is high as $\kappa _{0}\gg N_{0}L_{0}^{2}$ (i.e. $\tau _{{diff}}\ll \tau _{{int}}$ ), the gravitational force will be weakened due to the rapidly increased temperature of the fluid parcel. In this case, the parcel descends and stops at a position higher than the original position (figure 1 a). This implies that strong thermal diffusion can suppress the internal oscillation and the effect of stratification suppressing the vertical fluid motion becomes weak as a consequence.

Such a stratification effect affected by strong thermal diffusion has been considered in prescribing shear instability-driven turbulence and associated angular momentum transport in stellar radiation zones where fluids are stably stratified and the Prandtl number $Pr$ is low. For instance, for low Péclet number $Pe=RePr$ , where $Re$ is the Reynolds number, Lignières (1999) proposed a small-Péclet-number approximation that provides an asymptotic form of the Boussinesq equations in the small- $Pe$ limit where thermal diffusion and stratification can be combined and described as a single physical process. For vertical shear instability in stratified and highly diffusive fluids, Lignières et al. (1999) revealed that instability characteristics for various $N_{0}$ and $Pr\ll 1$ can be expressed in terms of a single rescaled number $R_{P}=RiPe$ , where $Ri=N_{0}^{2}/S_{0}^{2}$ is the Richardson number with a characteristic shear $S_{0}$ . Instead of the classical Richardson number criterion for shear instability in stratified flow, $Ri\lt 1/4$ , the modified Richardson number criterion $R_{P}\lt O(1)$ is adopted in the prescription of turbulent effective viscosity in the low- $Pe$ limit for simulations of the evolution of stars (e.g. Zahn 1992; Mathis et al. 2018; Prat & Mathis 2021). For horizontal shear flows in stratified-rotating fluids, Park et al. (2020, 2021b ) reported similar dynamical behaviours such as the stratification effect suppressed by high thermal diffusivity or instability characteristics described by the rescaled parameter $R_{P}$ for shear instabilities in the low- $Pe$ limit. The small-Péclet-number approximation with the rescaled parameter $R_{P}$ is also used in turbulence simulations of horizontal and vertical shear flows in stratified and highly diffusive fluids (Prat & Lignières 2013, 2014) and more attention has been paid recently to the low- $Pr$ and low- $Pe$ regimes in stratified turbulence (Cope et al. 2020; Garaud 2020a ; Chang & Garaud 2021; Garaud et al. 2024b ). It is noteworthy that Péclet numbers of actual astrophysical systems are too high (e.g. $Pe\sim O (10^{7} )$ with $Re\sim O (10^{13} )$ and $Pr\sim O (10^{-6} )$ for the solar tachocline; see also Garaud 2020b ) and fully resolved astrophysical turbulence for such high $Re$ and $Pe$ is not achievable yet from any state-of-the-art simulations. The current study also cannot address high- $Re$ /high- $Pe$ turbulence but still aims to explain the role of strong thermal diffusion at low $Pr$ .

Figure 1. (a) Illustration of how the internal oscillation of a fluid parcel in stably stratified fluid is suppressed by a fast thermal diffusion process. (b) Schematic of Taylor–Couette flow with stable temperature stratification.

Despite the increasing interest in the effect of strong thermal diffusion on stratified shear flows, the Prandtl-number dependence has not been studied considerably in the context of Taylor–Couette flow, a canonical shear flow between two concentric cylinders rotating independently, in a stably stratified fluid (figure 1 b). Stratification in Taylor–Couette flow tends to suppress the vertical fluid motion and, as a result, the onset of centrifugal instability is delayed and the vertical length scale of axisymmetric Taylor vortices is reduced as the stratification becomes strong (Withjack & Chen 1975; Boubnov et al. 1995; Caton et al. 2000). An interesting phenomenon in stratified Taylor–Couette flow is non-axisymmetric strato-rotational instability (SRI), which occurs due to the resonance between inertia–gravity waves confined between the two cylinders (Molemaker et al. 2001; Yavneh et al. 2001). The SRI has been explored extensively by theoretical investigations (e.g. Rüdiger & Shalybkov 2009; Le Dizès & Riedinger 2010; Park & Billant 2013; Leclercq et al. 2016; Wang & Balmforth 2018; Robins et al. 2020), experiments (e.g. Le Bars & Le Gal 2007; Ibanez et al. 2016; Rüdiger et al. 2017; Park et al. 2018; Seelig et al. 2018) and numerical simulations (e.g. Lopez & Marques 2020; Meletti et al. 2021; Lopez & Marques 2022). These instability studies considered the Prandtl number of $O(1)$ for fluids like water or air or simply $Pr=\infty$ to neglect the effect of thermal diffusion in theoretical analyses, or they considered the analogous Schmidt number $Sc=\nu _{0}/D_{0}$ (where $D_{0}$ is a diffusivity of scalars like density, salinity, etc.) of $O(100)$ if stratification of water by salt is considered in experiments. Investigating the dynamics of stratified Taylor–Couette flow under the influence of thermal diffusion, in particular in the low- $Pr$ limit, is important as the results can further be used to provide insights into large-scale flows and propose turbulent viscosity models in multi-physics simulations of astrophysical and geophysical systems such as simulations of the evolution of stars (Richard & Zahn 1999; Dubrulle et al. 2005). However, Taylor–Couette flow in stratified and highly diffusive fluids is still poorly understood and this motivates our study of stratified Taylor–Couette flow and firstly its centrifugal instability under the influence of strong thermal diffusion at low Prandtl number $Pr\leqslant 1$ .

The paper consists of the following sections to investigate the centrifugal instability of Taylor–Couette flow in stratified and diffusive fluids. In § 2, the problem formulation and details of numerical methods are provided. In § 3, one-dimensional (1-D) local linear stability analysis (LSA) is performed to present LSA results on centrifugal instability of cylindrical Couette flow for various parameters. In § 4, two-dimensional (2-D) and three-dimensional (3-D) direct numerical simulations (DNS) are conducted to investigate the nonlinear dynamics of centrifugal instability such as saturation, secondary instability or transition to chaotic states. In § 5, conclusion and discussion are provided.

2. Problem formulation and methodology

2.1. Navier–Stokes equations under the Boussinesq approximation

In this study, we consider the Boussinesq approximation in which the reference density $\rho _{0}$ is assumed to be much larger than the density variation $\varrho$ (i.e. $\rho _{0}\gg \varrho$ ). The density variation $\varrho$ is assumed to satisfy a linear relation with the total temperature $\vartheta$ as $\varrho /\rho _{0}=-\alpha _{0} (\vartheta -\vartheta _{0})$ where $\alpha _{0}\gt 0$ is the thermal expansion coefficient and $\vartheta _{0}$ is the reference temperature. For velocity ${\boldsymbol{U}}=(U_{r},U_{\theta },U_{z})$ , temperature variation $\varTheta =\vartheta -\vartheta _{0}$ and associated pressure variation $P$ in cylindrical coordinates $(r,\theta ,z)$ , we consider the following continuity, momentum and energy equations:

(2.1) \begin{align} \nabla \cdot {\boldsymbol{U}}=0, \\[-24pt] \nonumber \end{align}

(2.2) \begin{align} \frac {\partial {\boldsymbol{U}}}{\partial t}+{\boldsymbol{U}}\cdot \nabla {\boldsymbol{U}}=-\frac {1}{\rho _{0}}\nabla P +\alpha _{0}g\varTheta \vec{\boldsymbol{e}}_{z}+\nu _{0}\nabla ^{2}{\boldsymbol{U}}, \\[-24pt] \nonumber \end{align}

(2.3) \begin{align} \frac {\partial \varTheta }{\partial t}+{\boldsymbol{U}}\cdot \nabla \varTheta =\kappa _{0}\nabla ^{2}\varTheta , \\[-2pt] \nonumber \end{align}

where $g$ is the gravitational acceleration in the vertical direction $z$ , $\nu _{0}$ is the kinematic viscosity, $\kappa _{0}$ is the thermal diffusivity and $\nabla ^{2}$ is the Laplacian operator. We consider cylindrical Couette flow in a stably stratified fluid as a base state with base velocity $\boldsymbol{U}_{B}= (0,V_{B}(r),0 )$ and base temperature $T_{B}(z)$ as

(2.4) \begin{align} V_{B}(r)=Ar+\frac {B}{r},\quad A=\varOmega _{i}\frac {\mu -\eta ^{2}}{1-\eta ^{2}},\quad B=\varOmega _{i}R_{i}^{2}\frac {1-\mu }{1-\eta ^{2}},\quad T_{B}(z)=\frac {\Delta T}{\Delta z}z, \end{align}

where $V_{B}(r)$ is the base azimuthal velocity, $A$ and $B$ are constants as a function of the angular velocities $\varOmega _{i}$ and $\varOmega _{o}$ and radii $R_{i}$ and $R_{o}$ (where the subscripts $i$ and $o$ denote the inner and outer cylinders, respectively), $\mu =\varOmega _{o}/\varOmega _{i}$ is the angular velocity ratio, $\eta =R_{i}/R_{o}$ is the radius ratio and $\Delta T$ is the temperature difference along the vertical distance $\Delta z$ . Throughout this paper, we consider only the case in which the outer cylinder is fixed and only the inner cylinder rotates (i.e. $\varOmega _{o}=0$ and $\mu =0$ ). The temperature gradient $\Delta T/\Delta z$ is assumed to be positive and constant so the fluid is stably stratified and the base temperature $T_{B}$ increases linearly with $z$ . The base pressure $P_{B}(r,z)$ balances the velocity $V_{B}$ and temperature $T_{B}$ by satisfying the following relations:

(2.5) \begin{align} -\frac {V_{B}^{2}}{r}=-\frac {1}{\rho _{0}}\frac {\partial P_{B}}{\partial r},\quad 0=-\frac {1}{\rho _{0}}\frac {\partial P_{B}}{\partial z}+\alpha _{0}g T_{B}. \end{align}

Equations (2.1)–(2.3) can be expressed in a non-dimensional form by considering $R_{i}$ as the reference length, $\varOmega _{i}^{-1}$ as the reference time, $R_{i}\varOmega _{i}$ as the reference velocity, $\rho _{0}R_{i}^{2}\varOmega _{i}^{2}$ as the reference pressure and $\Delta T$ as the reference temperature. The distance $\Delta z$ can be chosen arbitrarily without loss of generality; thus we choose $\Delta z=R_{i}$ for simplicity. The coordinates $z$ and $r$ are considered non-dimensional hereafter, and hence the base temperature can simply be expressed as $T_{B}(z)=z$ .

We consider perturbation with its velocity $\boldsymbol{u}={\boldsymbol{U}}-\boldsymbol{U}_{B}= (u_{r},u_{\theta },u_{z} )$ , pressure $p=P-P_{B}$ and temperature $T=\varTheta -T_{B}$ . By applying the base state and perturbation to (2.1)–(2.3) and subtracting the base-state equations (2.4)–(2.5), we obtain the Navier–Stokes equations for perturbation as follows:

(2.6) \begin{align} \frac {\partial u_{r}}{\partial r}+\frac {u_{r}}{r}+\frac {1}{r}\frac {\partial u_{\theta }}{\partial \theta }+\frac {\partial u_{z}}{\partial z}=0, \\[-26pt] \nonumber \end{align}

(2.7) \begin{align} \frac {\partial u_{r}}{\partial t}+\varOmega \frac {\partial u_{r}}{\partial \theta }-2\varOmega u_{\theta }+{N}_{r}=-\frac {\partial p}{\partial r}+\frac {1}{Re}\left (\nabla ^{2}u_{r}-\frac {u_{r}}{r^{2}}-\frac {2}{r^{2}}\frac {\partial u_{\theta }}{\partial \theta }\right ), \\[-26pt] \nonumber \end{align}

(2.8) \begin{align} \frac {\partial u_{\theta }}{\partial t}+\varOmega \frac {\partial u_{\theta }}{\partial \theta }+Zu_{r}+{N}_{\theta }=-\frac {1}{r}\frac {\partial p}{\partial \theta }+\frac {1}{Re}\left (\nabla ^{2}u_{\theta }-\frac {u_{\theta }}{r^{2}}+\frac {2}{r^{2}}\frac {\partial u_{r}}{\partial \theta }\right ), \\[-26pt] \nonumber \end{align}

(2.9) \begin{align} \frac {\partial u_{z}}{\partial t}+\varOmega \frac {\partial u_{z}}{\partial \theta }+{N}_{z}=-\frac {\partial p}{\partial z}+N^{2}T+\frac {1}{Re}\nabla ^{2}u_{z}, \\[-26pt] \nonumber \end{align}

(2.10) \begin{align} \frac {\partial T}{\partial t}+\varOmega \frac {\partial T}{\partial \theta }+u_{z}+{N}_{T}=\frac {1}{RePr}\nabla ^{2}T, \\[1pt] \nonumber \end{align}

where $\varOmega (r)=V_{B}/r$ is the base angular velocity, $Z(r)=(1/r)\mathrm {d}(r^{2}\varOmega )/\mathrm {d}r$ is the base axial vorticity and $N_{r}$ , $N_{\theta }$ , $N_{z}$ and $N_{T}$ are nonlinear terms:

(2.11) \begin{align} {N}_{r}=u_{r}\frac {\partial u_{r}}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{r}}{\partial \theta }+u_{z}\frac {\partial u_{r}}{\partial z}-\frac {u_{\theta }^{2}}{r}, \\[-26pt] \nonumber \end{align}

(2.12) \begin{align} {N}_{\theta }=u_{r}\frac {\partial u_{\theta }}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{\theta }}{\partial \theta }+u_{z}\frac {\partial u_{\theta }}{\partial z}+\frac {u_{r}u_{\theta }}{r}, \\[-26pt] \nonumber \end{align}

(2.13) \begin{align} {N}_{z}=u_{r}\frac {\partial u_{z}}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{z}}{\partial \theta }+u_{z}\frac {\partial u_{z}}{\partial z}, \\[-26pt] \nonumber \end{align}

(2.14) \begin{align} {N}_{T}=u_{r}\frac {\partial T}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial T}{\partial \theta }+u_{z}\frac {\partial T}{\partial z}. \\[6pt] \nonumber \end{align}

The parameters $Re$ , $N$ and $Pr$ are the Reynolds number, non-dimensional Brunt–Väisälä frequency and Prandtl number, respectively, which are defined as

(2.15) \begin{align} Re=\frac {R_{i}^{2}\varOmega _{i}}{\nu _{0}},\quad N=\sqrt {\frac {\alpha _{0}g}{\varOmega _{i}^{2}}\frac {\Delta T}{\Delta z}},\quad Pr=\frac {\nu _{0}}{\kappa _{0}}. \end{align}

For convenience in comparison with other literature, we also use a conventional inner-cylinder-based Reynolds number $Re_{i}$ defined as

(2.16) \begin{align} Re_{i}=\frac {R_{i}\varOmega _{i}d}{\nu _{0}}=Re\left (\frac {1-\eta }{\eta }\right ), \end{align}

where $d=R_{o}-R_{i}$ is the gap size between the two cylinders.

2.2. Pseudo-spectral formulation for direct numerical simulation

In this study, spectral-based DNS are conducted to solve (2.6)–(2.10) numerically. To do so, we decompose the perturbation into the sum of modes using the following Fourier representation:

(2.17) \begin{align} \left ( \begin{array}{c} u_{r}\\ u_{\theta }\\ u_{z}\\ T\\ p \end{array} \right ) = \sum _{j=-M}^{M} \sum _{l=-K}^{K} \left ( \begin{array}{c} \tilde {u}_{\! jl}(r,t)\\ \tilde {v}_{\! jl}(r,t)\\ \tilde {w}_{\! jl}(r,t)\\ \tilde {T}_{\! jl}(r,t)\\ \tilde {p}_{\! jl}(r,t) \end{array} \right )\exp \left (\mathrm {i}m_{\! j}\theta +\mathrm {i}k_{l}z\right ), \end{align}

where $M$ and $K$ are the cut-off numbers of modes considered in the azimuthal direction $\theta$ and axial direction $z$ , respectively, $\tilde {u}_{\! jl}$ , $\tilde {v}_{\! jl}$ , $\tilde {w}_{\! jl}$ , $\tilde {p}_{\! jl}$ and $\tilde {T}_{\! jl}$ are the time-dependent mode shapes, $m_{j}=jm$ is the $j$ th multiple of the principal azimuthal wavenumber $m$ and $k_{l}=lk$ is the $l$ th multiple of the principal axial wavenumber $k$ . The above ansatz (2.17) is used for various problems such as the Taylor–Couette formulation in nsCouette code (López et al. 2020) or convection problems (Saltzman 1962; Park et al. 2021a ). The formulation is also analogous to semi-linear models, which are applied to centrifugal instability of anti-cyclonic vortices (Yim et al. 2020, 2023). The semi-linear theory allows us to investigate directly nonlinear interaction between base flow and an instability mode, and the method is generalised by considering nonlinear interaction among multiple instability modes in low-order harmonics while neglecting the triad interaction leading to harmonics of orders higher than the cut-off numbers. For each mode with indices $j$ and $l$ , we apply the ansatz (2.17) to (2.6)–(2.10) and obtain

(2.18) \begin{align} \frac {\partial \tilde {u}_{\! jl}}{\partial r}+\frac {\tilde {u}_{\! jl}}{r}+\frac {\mathrm {i}m_{\! j}\tilde {v}_{\! jl}}{r}+\mathrm {i}k_{l}\tilde {w}_{\! jl}=0, \\[-26pt] \nonumber \end{align}

(2.19) \begin{align} \frac {\partial \tilde {u}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {u}_{\! jl}-2\varOmega \tilde {v}_{\! jl}+\tilde {N}_{r,jl}=-\frac {\partial \tilde {p}_{\! jl}}{\partial r}+\frac {1}{Re}\left (\tilde {\nabla }_{jl}^{2}\tilde {u}_{\! jl}-\frac {\tilde {u}_{\! jl}}{r^{2}}-\frac {2\mathrm {i}m_{\! j}\tilde {v}_{\! jl}}{r^{2}}\right ), \\[-26pt] \nonumber \end{align}

(2.20) \begin{align} \frac {\partial \tilde {v}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {v}_{\! jl}+Z\tilde {u}_{\! jl}+\tilde {N}_{\theta ,jl}=-\frac {\mathrm {i}m_{\! j}\tilde {p}_{\! jl}}{r}+\frac {1}{Re}\left (\tilde {\nabla }_{jl}^{2}\tilde {v}_{\! jl}-\frac {\tilde {v}_{\! jl}}{r^{2}}+\frac {2\mathrm {i}m_{\! j}\tilde {u}_{\! jl}}{r^{2}}\right ), \\[-26pt] \nonumber \end{align}

(2.21) \begin{align} \frac {\partial \tilde {w}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {w}_{\! jl}+\tilde {N}_{z,jl}=-\mathrm {i}k_{l}\tilde {p}_{\! jl}+N^{2}\tilde {T}_{\! jl}+\frac {1}{Re}\tilde {\nabla }_{jl}^{2}\tilde {w}_{\! jl}, \\[-26pt] \nonumber \end{align}

(2.22) \begin{align} \frac {\partial \tilde {T}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {T}_{\! jl}+\tilde {w}_{\! jl}+\tilde {N}_{T,jl}=\frac {1}{RePr}\tilde {\nabla }_{jl}^{2}\tilde {T}_{\! jl}, \end{align}

where $\tilde {N}_{r,jl}$ , $\tilde {N}_{\theta ,jl}$ , $\tilde {N}_{z,jl}$ and $\tilde {N}_{T,jl}$ are the terms convoluted from the nonlinear terms (2.11)–(2.14) using the Fourier transform and $\tilde {\nabla }_{jl}^{2}=\partial ^{2}/\partial r^{2}+(1/r)\partial /\partial r-k_{l}^{2}-m_{j}^{2}/r^{2}$ is the modal Laplacian operator. For the boundary conditions, we consider

(2.23) \begin{align} \tilde {u}_{\! jl}=\tilde {v}_{\! jl}=\tilde {w}_{\! jl}=\frac {\partial \tilde {T}_{\! jl}}{\partial r}=0, \end{align}

at both cylinders $r=1$ and $r=1/\eta$ .

An advantage of using the equations in a modal form is that (2.18)–(2.22) can further be simplified by eliminating the pressure $\tilde {p}_{\! jl}$ as follows:

(2.24) \begin{align} \mathcal {A}_{jl}\frac {\partial \tilde {\boldsymbol{q}}_{\! jl}}{\partial t}=\mathcal {B}_{jl}\tilde {\boldsymbol{q}}_{\! jl}+\tilde {\boldsymbol{N}}_{\! jl}, \end{align}

where $\tilde {\boldsymbol{q}}_{\! jl}$ is the vector of three variables, $\mathcal {A}_{jl}$ and $\mathcal {B}_{jl}$ are the differential operator matrices and $\tilde {\boldsymbol{N}}_{\! jl}$ is the vector comprised of the nonlinear terms, all of which depend on indices $j$ and $l$ such as whether $l=0$ or not and whether $j=0$ or not. We refer the reader to Appendix A for more details of these vectors and matrices. For spatial discretisation in the radial direction $r$ , a spectral method using the Chebyshev collocation points is considered (Antkowiak & Brancher 2004; Park 2012). The vector $\tilde {\boldsymbol{N}}_{\! jl}$ is computed at each step using a pseudo-spectral method, similar to that in Deloncle et al. (2008). Equation (2.24) is written in a spectral form in the azimuthal and axial directions and we use the backward Fourier transform to obtain the vector $\boldsymbol{q}_{\! jl}$ in the physical space $(r,\theta ,z)$ . Then $\boldsymbol{q}_{\! jl}$ is used to compute the nonlinear term $\boldsymbol{N}_{\! jl}$ in the physical space and apply the forward Fourier transform to $\boldsymbol{N}_{\! jl}$ to compute $\tilde {\boldsymbol{N}}_{\! jl}$ in the spectral space $(r,m_{j},k_{l})$ . In this pseudo-spectral approach, we consider the numbers of collocation points $N_{\theta }=2M+1$ and $N_{z}=2K+1$ in the azimuthal and axial directions, respectively. The pseudo-spectral approach is similar to that in Guseva et al. (2015) who also consider the axial periodicity with the Fourier method but utilise a fractional time-stepping method by carefully computing an intermediate pressure using the influence-matrix method, which is especially effective in solving equations for the magnetic field. On the contrary, we avoid the use of the pressure $\tilde {p}_{\! jl}$ by considering different operator matrices $\mathcal {A}_{jl}$ and $\mathcal {B}_{jl}$ that depend on the indices $j$ and $l$ . Unlike Deloncle et al. (2008), the de-aliasing technique, which is required in turbulence simulations, is not implemented as the Reynolds number $Re_{i}$ considered in this work is not large enough to observe small-scale turbulence. For time advancement, we use a conventional semi-implicit scheme (see e.g. Kim et al. 1987) with the Crank–Nicolson method for the linear term and the Adam–Bashforth method for the nonlinear term. For instance, at each step $n$ , we find the next step solution $\tilde {\boldsymbol{q}}_{\! jl}^{(n+1)}$ by solving numerically the discretised version of (2.24):

(2.25) \begin{align} \left (\mathcal {A}_{jl}-\frac {\Delta t}{2}\mathcal {B}_{jl}\right )\tilde {\boldsymbol{q}}_{\! jl}^{(n+1)}=\left (\mathcal {A}_{jl}+\frac {\Delta t}{2}\mathcal {B}_{jl}\right )\tilde {\boldsymbol{q}}_{\! jl}^{(n)}+\frac {\Delta t}{2}\left (3\tilde {\boldsymbol{N}}_{\! jl}^{(n)}-\tilde {\boldsymbol{N}}_{\! jl}^{(n-1)}\right ), \end{align}

where $\Delta t$ is the time step. In this study, the time step $\Delta t$ is fixed to $\Delta t=0.01$ , which is found to be sufficiently small for parameters considered in this study. It is verified that the Courant–Friedrichs–Lewy condition for perturbation velocity is satisfied and every DNS demonstrates no numerical divergence with this $\Delta t$ .

2.3. One-dimensional local linear stability analysis

By assuming that the perturbation is infinitesimally small, we can neglect nonlinear terms and perform 1-D local LSA using the normal mode:

(2.26) \begin{align} \left ( \begin{array}{c} u_{r}\\ u_{\theta }\\ u_{z}\\ T\\ p \end{array} \right ) = \left ( \begin{array}{c} \hat {u}(r)\\ \hat {v}(r)\\ \hat {w}(r)\\ \hat {T}(r)\\ \hat {p}(r) \end{array} \right )\exp \left (\mathrm {i}m\theta +\mathrm {i}kz-\mathrm {i}\omega t\right )+\text{c.c.}, \end{align}

where $\text{c.c.}$ denotes the complex conjugate, $\hat {u}$ , $\hat {v}$ , $\hat {w}$ , $\hat {p}$ and $\hat {T}$ are the mode shapes, $m$ is the azimuthal wavenumber, $k$ is the axial wavenumber and $\omega$ is the complex frequency $\omega =\omega _{r}+\mathrm {i}\omega _{i}$ , where $\omega _{r}={\textrm {Re}}(\omega )$ is the temporal frequency and $\omega _{i}={\textrm {Im}}(\omega )$ is the temporal growth rate. Throughout the paper, the non-dimensional wavenumber $k_{d}=kd/R_{i} = k(1-\eta )/\eta$ rescaled by the gap size $d$ is also used for convenience in comparison with other literature. The normal mode (2.26) is applied to (2.6)–(2.10) with the nonlinear terms neglected and, after manipulations, we can eliminate the pressure mode shape and obtain the following simplified eigenvalue problem:

(2.27) \begin{align} -\mathrm {i}\omega \mathcal {A}\hat {\boldsymbol{q}}=\mathcal {B}\hat {\boldsymbol{q}}, \end{align}

where $\hat {\boldsymbol{q}}= (\hat {u},\hat {v},\hat {T} )^{\mathrm {T}}$ and $\mathcal {A}$ and $\mathcal {B}$ are the operator matrices, which are essentially the same as $\mathcal {A}_{11}$ and $\mathcal {B}_{11}$ , respectively, the operator matrices in (2.24) with $j=l=1$ . To solve the eigenvalue problem (2.27), the MATLAB routine eig is used with the following boundary conditions imposed at both cylinders $r=1$ and $r=1/\eta$ :

(2.28) \begin{align} \hat {u}=\hat {v}=\hat {w}=\frac {\mathrm {d}\hat {T}}{\mathrm {d}r}=0. \end{align}

The Chebyshev spectral method is used for discretisation in the radial direction $r$ and a number of collocation points $N_{r}$ between 60 and 120 is considered in this study. The choice of $N_{r}$ depends on how large are the parameters such as the Reynolds number $Re$ or the Péclet number $Pe=RePr$ , how the eigenmodes are confined near the boundaries, etc. For more details of the LSA and numerical methods, we refer the reader to our previous work (e.g. Park 2012; Park & Billant 2013) that used the same code.

2.4. Two-dimensional bi-global linear stability analysis

For parameters considered in this study, 1-D LSA reveals that the axisymmetric mode with $m=0$ is the most unstable one for cylindrical Couette flow. As demonstrated in the next sections, centrifugal instability develops nonlinearly and saturates leading to axisymmetric Taylor vortices as a new base state. The new 2-D base state $\bar {\boldsymbol{Q}}(r,z)=(\bar {\boldsymbol{U}},\bar {T})$ is comprised of the new base velocity $\bar {\boldsymbol{U}}(r,z)= (\bar {U}(r,z),\bar {V}(r,z),\bar {W}(r,z) )$ and the new base temperature $\bar {T}(r,z)$ . For this Taylor vortex flow, we can analyse its secondary instability through 2-D bi-global LSA by considering a non-axisymmetric perturbation $\bar {\textbf {q}}$ with $m\neq 0$ expressed in the following ansatz:

(2.29) \begin{align} \bar {\textbf {q}}= \left ( \begin{array}{c} \bar {u}_{r}\\ \bar {u}_{\theta }\\ \bar {u}_{z}\\ \bar {T}\\ \bar {p} \end{array} \right ) = \left ( \begin{array}{c} \hat {u}_{m}(r,z)\\ \hat {v}_{m}(r,z)\\ \hat {w}_{m}(r,z)\\ \hat {T}_{m}(r,z)\\ \hat {p}_{m}(r,z) \end{array} \right )\exp \left (\mathrm {i}m\theta -\mathrm {i}\omega _{m} t\right )+\text{c.c.}, \end{align}

where $\hat {u}_{m}$ , $\hat {v}_{m}$ , $\hat {w}_{m}$ , $\hat {T}_{m}$ and $\hat {p}_{m}$ are the 2-D global mode shapes and $\omega _{m}$ is the complex frequency of the global mode. Similar to the 1-D local LSA, an eigenvalue problem can be formulated as

(2.30) \begin{align} -\mathrm {i}\omega _{m}\mathcal {A}_{m}\hat {\textbf {q}}_{m}=\mathcal {B}_{m}\hat {\textbf {q}}_{m}, \end{align}

where $\mathcal {A}_{m}$ and $\mathcal {B}_{m}$ are the operator matrices detailed in Appendix A.2 and $\hat {\textbf {q}}_{m}= (\hat {u}_{m},\hat {v}_{m},\hat {T}_{m} )^{\mathrm {T}}$ is the mode shape. The matrices $\mathcal {A}_{m}$ and $\mathcal {B}_{m}$ in the 2-D bi-global LSA have a much larger size than the size of the matrices $\mathcal {A}$ and $\mathcal {B}$ in (2.27). Thus we use the MATLAB routine eigs that computes only a few eigenvalues near a specified value using the Krylov–Schur algorithm (Stewart 2002). To solve the eigenvalue problem (2.30), we consider the following boundary conditions:

(2.31) \begin{align} \hat {u}_{m}=\hat {v}_{m}=\hat {w}_{m}=\frac {\partial \hat {T}_{m}}{\partial r}=0. \end{align}

The Chebyshev and Fourier spectral methods are used for discretisation in the radial and axial directions, respectively. The Fourier method imposes the periodic boundary condition in the axial direction $z$ given that the axisymmetric Taylor vortices as a new base state are also periodic in $z$ with wavelength $\lambda _{z}=2\pi /k$ .

3. Linear centrifugal instability in stratified and diffusive fluids

Linear instability of stratified flows in thermally diffusive fluids, especially with low $Pr$ , has been investigated for various planar shear flows in a linear profile (Barker et al. 2019, 2020; Dymott et al. 2023), a periodic sinusoidal profile (Chang & Garaud 2021; Garaud et al. 2024b ) and a hyperbolic tangent profile (Park et al. 2020, 2021b ). They consider flows with vertical shear, horizontal shear or mixed shear (where the ‘vertical’ direction in these studies implies the direction of gravity and stratification) and investigate linear and nonlinear properties of shear instabilities. In this section, we similarly investigate linear centrifugal instability of Taylor–Couette flow, a horizontally sheared rotating flow, in stratified and diffusive fluids. The linear analysis results are followed by nonlinear simulation results presented in § 4.

3.1. Neutral stability curves

Figure 2(a) shows the growth rate $\omega _{i}$ of the axisymmetric mode ( $m=0$ ), which is found to be most unstable, versus the wavenumber $k_{d}$ for different parameter sets of $(N,Pr)$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ . For every case presented in figure 2(a), the corresponding frequency $\omega _{r}$ is zero. It is known that centrifugal instability reaches its maximum $\omega _{i,max }$ as $k\rightarrow \infty$ in the inviscid limit $Re\rightarrow \infty$ as the growth rate scales as $\omega _{i,{inviscid}}=\omega _{i,max } -A_{0}/k^{2}$ (Billant & Gallaire 2005; Park et al. 2017) while the viscous growth rate scales as $\omega _{i,{viscous}}=\omega _{i,{inviscid}}-A_{1}k^{2}/Re$ (Yim et al. 2016), where $A_{0}$ and $A_{1}$ are positive constants that depend on stratification and other parameters. Due to this characteristic varying with $k$ , the viscous growth rates are positive in a finite range of $k_{d}$ and each curve reaches its peak at a certain wavenumber $k_{d,max }$ . As the stratification increases from $N=0$ (grey) to $N=1$ (black) for $Pr=1$ , the growth-rate curves descend while the wavenumber $k_{d,max }$ at the peak of $\omega _{i}$ increases. For fixed $N=1$ , as the Prandtl number $Pr$ decreases from $Pr=1$ , the growth rate increases and, remarkably, the growth-rate curves for $Pr=10^{-2}$ and $10^{-4}$ overlap with the grey curve, the unstratified case with $N=0$ . This implies that the centrifugal instability in stratified and highly diffusive fluids with $Pr\ll 1$ behaves as the instability in unstratified fluids as the effect of stratification is suppressed by strong thermal diffusion. Another remarkable result in figure 2(a) is that the growth-rate curves overlap for $(N,Pr)=(10,10^{-2})$ and $(10^{2},10^{-4})$ , the cases with the same $P_{N}=N^{2} Pr=1$ . The growth rate for other values of $N$ and $Pr\ll 1$ is invariant if the rescaled parameter $P_{N}$ is the same. Such an invariant feature at the same $P_{N}$ but different $N$ and $Pr$ is similarly reported for both vertical and horizontal shear instabilities in vertically stratified fluids (Lignières et al. 1999; Park et al. 2020, 2021b ).

Figure 2. (a) Growth-rate curves for various sets of $(N,Pr)$ at $\mu =0$ , $\eta =0.9$ , $Re_{i}=200$ and $m=0$ . (b) Neutral stability curves in the parameter space $(N,Re_{i})$ for different $Pr$ at $\mu =0$ , $\eta =0.9$ and $m=0$ . (c) The curves for different $Pr$ same as (b) but over a wider range of $N$ . (d) The curves same as (c) overlapped due to the rescaled parameter $P_{N}=N^{2}Pr$ on the abscissa. An additional thick grey line is a neutral stability curve obtained from the small- $Pr$ approximation.

Figure 2 $(b)$ display neutral stability curves, which denote the critical Reynolds number $Re_{i,c}$ at which the growth rate $\omega _{i}$ of the most unstable mode is zero, in the parameter space $(N,Re_{i})$ for different values of $Pr$ at $\mu =0$ , $\eta =0.9$ and $m=0$ . For $Pr\geqslant 1$ , the neutral stability curves increase rapidly with $N$ , a feature similarly found in Park et al. (2017), and the $Pr=1$ case is more stable than other cases with $Pr\gt 1$ when $N\gt 3$ . This is expected as thermal dissipation is proportional to the diffusivity $\kappa _{0}$ ; thus lower diffusivity $\kappa _{0}$ (i.e. higher $Pr$ ) implies less thermal dissipation and more instability. In the range $0\lt N\lt 3$ , the situation is more complicated since the $Pr=10$ case is more stable than other $Pr$ cases for the axisymmetric perturbation. Although such high- $Pr$ dynamics at a moderate $N$ should be further investigated due to its great importance in geophysical and other contexts (e.g. oceanic flows where the analogous Schmidt number $Sc$ is around 700), the current study will only focus on the highly diffusive regime with low $Pr\leqslant 1$ . For $Pr\lt 1$ , the curves increase very slowly as $N$ increases and it is difficult to see in figure 2(b) the increase of the curves at $Pr=10^{-4}$ and $10^{-6}$ over the range $0\leqslant N\leqslant 20$ . Figure 2(c) displays the same neutral stability curves for $Pr\leqslant 1$ over a wider range of $N$ to see how slowly the neutral stability curves increase with $N$ as $Pr$ decreases. This increasing trend confirms that, while stratification enhances the stability of stratified Taylor–Couette flow, strong thermal diffusion destabilises by preventing the stabilising role of stratification and making the flow behave like an unstratified flow. In figure 2(d), we plot again the same neutral stability curves but over the rescaled parameter $P_{N}$ on the abscissa. All the curves now overlap each other, even for the case with $Pr=1$ , and can be described by a linear relation as $Re_{i,c}= 131.6+8.965P_{N}$ , where $Re_{i,c}=131.6$ is the critical Reynolds number for the unstratified case $N=0$ at $m=0$ and $\eta =0.9$ , the number agreeing with DiPrima et al. (1984).

Figure 3. (a,b) Real (solid) and imaginary (dashed) parts of the mode shape $\hat {v}(r)$ and rescaled mode $\hat {T}(r)/Pr$ for $Pr=1$ (red), $Pr=10^{-2}$ (blue) and $Pr=10^{-4}$ (grey) at $\mu =0$ , $\eta =0.9$ , $Re_{i}=200$ , $N=1$ , $m=0$ and $k_{d}=3.91$ . (c) Perturbation temperature $T(r,z)$ reconstructed from $\hat {T}(r)$ for $Pr=10^{-2}$ .

Figure 3 shows the real and imaginary parts of eigenmode shapes $\hat {v}(r)$ and $\hat {T}(r)$ for different $Pr$ for the wavenumber set $(m,k_{d})=(0,3.91)$ . The mode shapes are normalised by the maximum value of $\hat {v}$ and it is found that the mode shape $\hat {T}$ has a smaller amplitude than $\hat {v}$ for small $Pr$ . As the Prandtl number $Pr$ decreases, the amplitude of $\hat {T}$ decreases and scales as $O(Pr)$ . For a better comparison, figure 3(b) displays the mode shape $\hat {T}(r)$ divided by $Pr$ and we see that the rescaled $\hat {T}/Pr$ has the same mode shape for the cases with $Pr=10^{-2}$ and $Pr=10^{-4}$ . In figure 3(c), the perturbation temperature $T(r,z)=\hat {T}\exp (\mathrm {i}kz)+\hat {T}^{*}\exp (-\mathrm {i}kz)$ for $Pr=10^{-2}$ is plotted. Two in-phase waves, one near the inner cylinder and the other near the outer one, are clearly shown. The shape of this axisymmetric mode is different from that of a non-axisymmetric centrifugal instability mode, which is weakly sheared, or that of a SRI mode that is out of phase (Park et al. 2017, 2018). The perturbation temperature has one node (i.e. the zero crossing) around $r\simeq 1.053$ as it is the first mode which becomes the most unstable. Higher-order modes with a greater number of nodes are, on the other hand, found to be stable for the parameter set in figure 3 with $Pr=10^{-2}$ .

3.2. Small- $Pr$ approximation

The single dependence on $P_{N}$ for different $N$ and $Pr$ in the limit $Pr\ll 1$ can be understood by taking the small- $Pr$ approximation. Consider the Taylor expansions in terms of small $Pr$ as $\hat {\textbf {u}}= (\hat {u},\hat {v},\hat {w} )^{\mathrm {T}}=\hat {\textbf {u}}^{(0)}+Pr\hat {\textbf {u}}^{(1)}+O (Pr^{2} )$ and similarly for $\hat {T}$ and $\hat {p}$ as $\hat {T}=\hat {T}^{(0)}+Pr\hat {T}^{(1)}+\cdots$ and $\hat {p}=\hat {p}^{(0)}+Pr\hat {p}^{(1)}+\cdots$ . For the Péclet number $Pe$ being small, $Pe=RePr\ll 1$ (i.e. $Re\sim O(1)$ and $Pr\ll 1$ ), we obtain the leading-order equation for $\hat {T}^{(0)}$ as

(3.1) \begin{align} \frac {1}{RePr}\hat {\nabla }^{2}\hat {T}^{(0)}=0, \\[-6pt] \nonumber \end{align}

which leads to the analytic solution $\hat {T}^{(0)}(r)=a_{1}\mathrm {I}_{m}(kr)+a_{2}\mathrm {K}_{m}(kr)$ , where $a_{1}$ and $a_{2}$ are constants and $\mathrm {I}_{m}$ and $\mathrm {K}_{m}$ are modified Bessel functions of the first and second kinds, respectively (Abramowitz & Stegun 1965). The constants $a_{1}$ and $a_{2}$ become zero if the no-flux conditions $\mathrm {d}\hat {T}^{(0)}/\mathrm {d}r=0$ are imposed at both cylinders; thus the leading-order solution simply becomes zero (i.e. $\hat {T}^{(0)}(r)=0$ ). At the next order, we obtain the following equations:

(3.2) \begin{align} \frac {\mathrm {d} \hat {u}^{(0)}}{\mathrm {d} r}+\frac {\hat {u}^{(0)}}{r}+\frac {\mathrm {i}m\hat {v}^{(0)}}{r}+\mathrm {i}k\hat {w}^{(0)}=0, \\[-24pt] \nonumber \end{align}

(3.3) \begin{align} -\mathrm {i}\omega \hat {u}^{(0)}+\mathrm {i}m\varOmega \hat {u}^{(0)}-2\varOmega \hat {v}^{(0)}=-\frac {\mathrm {d} \hat {p}^{(0)}}{\mathrm {d} r}+\frac {1}{Re}\left (\hat {\nabla }^{2}\hat {u}^{(0)}-\frac {\hat {u}^{(0)}}{r^{2}}-\frac {2\mathrm {i}m\hat {v}^{(0)}}{r^{2}}\right ), \\[-24pt] \nonumber \end{align}

(3.4) \begin{align} -\mathrm {i}\omega \hat {v}^{(0)}+\mathrm {i}m\varOmega \hat {v}^{(0)}+Z\hat {u}^{(0)}=-\frac {\mathrm {i}m\hat {p}^{(0)}}{r}+\frac {1}{Re}\left (\hat {\nabla }^{2}\hat {v}^{(0)}-\frac {\hat {v}^{(0)}}{r^{2}}+\frac {2\mathrm {i}m\hat {u}^{(0)}}{r^{2}}\right ), \\[-24pt] \nonumber \end{align}

(3.5) \begin{align} -\mathrm {i}\omega \hat {w}^{(0)}+\mathrm {i}m\varOmega \hat {w}^{(0)}=-\mathrm {i}k\hat {p}^{(0)}+P_{N}\hat {T}^{(1)}+\frac {1}{Re}\hat {\nabla }^{2}\hat {w}^{(0)}, \\[-24pt] \nonumber \end{align}

(3.6) \begin{align} \hat {w}^{(0)}=\frac {1}{Re}\hat {\nabla }^{2}\hat {T}^{(1)}, \\[-6pt] \nonumber \end{align}

where $P_{N}=N^{2}Pr$ . We see that the two parameters $N$ and $Pr$ representing thermal stratification and diffusion are simplified into a single parameter $P_{N}$ , as observed in other studies under the small-Péclet-number approximation (e.g. Lignières 1999; Park et al. 2020). By considering the continuity equation and eliminating the pressure, (3.2)–(3.6) can further be simplified into the following eigenvalue problem:

(3.7) \begin{align} -\mathrm {i}\omega \mathcal {A}^{(0)}\hat {\boldsymbol{q}}^{(0)}=\mathcal {B}^{(0)}\hat {\boldsymbol{q}}^{(0)}, \end{align}

where $\hat {\boldsymbol{q}}^{(0)}$ is the mode shape vector and $\mathcal {A}^{(0)}$ and $\mathcal {B}^{(0)}$ are the operator matrices detailed in Appendix A. In figure 2(d), it is clearly shown that a neutral stability curve computed from (3.7) overlaps with other neutral stability curves when they are plotted over the rescaled Prandtl number $P_{N}=N^{2}Pr$ as the abscissa.

3.3. Perturbation energy analysis

The instability characteristics can also be understood by examining the evolution of perturbation energy. We first define the total perturbation energy $E(t)$ as

(3.8) \begin{align} E(t)=\frac {1}{2}\int _{0}^{L_{z}}\int _{0}^{2\pi }\int _{1}^{1/\eta }\left (u_{r}^{2}+u_{\theta }^{2}+u_{z}^{2}+N^{2}T^{2}\right )r\,\mathrm {d}r\,\mathrm {d}\theta \,\mathrm {d}z=\frac {1}{2}\left \lt \textbf {q}_{E}\,\textbf {:}\,\textbf {q}_{E}\right \gt , \end{align}

where $L_{z}$ is a vertical length of the domain assumed to be periodic (e.g. $L_{z}=2\pi /k$ if we consider one periodic length of a normal mode with the axial wavenumber $k$ ), the angle brackets denote the volume integral defined as $\left \lt \textbf {X}\right \gt =\int _{0}^{L_{z}}\int _{0}^{2\pi }\int _{1}^{1/\eta }\textbf {X} \,r\mathrm {d}r\mathrm {d}\theta \mathrm {d}z$ , $\textbf {q}_{E}= (u_{r},u_{\theta },u_{z},NT )^{\mathrm {T}}$ and $\textbf {:}$ denotes the Frobenius product (for more details, see also Park et al. 2017). For the momentum and energy equations (2.7)–(2.10), we multiply both sides by $ (u_{r}, u_{\theta }, u_{z}, N^{2}T )^{\mathrm {T}}$ and, after some manipulations, we obtain the following equation for perturbation energy:

(3.9) \begin{align} \frac {\partial E}{\partial t}=\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\left (ru_{r}u_{\theta }\right )-\frac {1}{Re}\left [\nabla \textbf {u}\,\textbf {:}\,\nabla \textbf {u}+\frac {N^{2}}{Pr}\left (\nabla T\,\textbf {:}\,\nabla T\right )\right ]\right \gt , \end{align}

where the first term on the right-hand side within the angle brackets denotes the contribution from the mean angular shear $\mathrm {d}\varOmega /\mathrm {d}r$ and the second term corresponds to the kinetic and potential energy dissipation. The contribution from the mean angular shear is a main source of energy production if the momentum transfer term $u_{r}u_{\theta }$ is anti-correlated with the mean angular shear $\mathrm {d}\varOmega /\mathrm {d}r$ , a well-known mechanism for instability called the Orr mechanism (Orr 1907). The momentum and thermal dissipation terms are always negative and thus stabilise the perturbation energy. We note that this mechanism is valid for both linear and nonlinear cases as the evolution equation (3.9) is derived from the nonlinear perturbation equations (2.7)–(2.10) and the nonlinear terms are cancelled out in the derivation process in which continuity is taken into account. If we apply the normal mode (2.26) to the evolution equation (3.9), we obtain the following expression for the growth rate:

(3.10) \begin{align} \omega _{i}=\frac {1}{\hat {E}}\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\frac {r\left (\hat {u}^{*}\hat {v}+\hat {u}\hat {v}^{*}\right )}{2}-\frac {1}{Re}\left [\hat {\nabla }\hat {\textbf {u}}\,\textbf {:}\,\hat {\nabla }\hat {\textbf {u}}+\frac {N^{2}}{Pr}\left (\hat {\nabla }\hat {T}\,\textbf {:}\,\hat {\nabla }\hat {T}\right )\right ]\right \gt _{r}, \end{align}

where $*$ denotes the complex conjugate, $\lt \gt _{r}$ denotes the line integral over the radial coordinate $r$ as $\left \lt \hat {\textbf {X}}\right \gt _{r}=\int _{1}^{1/\eta }\hat {\textbf {X}} \,r\mathrm {d}r$ and $\hat {E}=\left \lt |\hat {u}|^{2}+|\hat {v}|^{2}+|\hat {w}|^{2}+|N\hat {T}|^{2}\right \gt _{r}$ .

Table 1 presents examples of the maximum growth rate $\omega _{i,max }$ and corresponding wavenumber $k_{d,max }$ for various parameter sets of $(N,Pr,m)$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ . The table also details the production and dissipation terms $\mathcal {P}_{\Omega }$ , $\epsilon _{k}$ and $\epsilon _{p}$ , which are defined as

(3.11) \begin{align} \mathcal {P}_{\Omega }=\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\frac {r\left (\hat {u}^{*}\hat {v}+\hat {u}\hat {v}^{*}\right )}{2}\right \gt _{r},\quad \epsilon _{k}=-\frac {1}{Re}\left \lt \hat {\nabla }\hat {\textbf {u}}:\hat {\nabla }\hat {\textbf {u}}\right \gt _{r},\quad \epsilon _{p}=-\frac {N^{2}}{RePr}\left \lt \hat {\nabla }\hat {T}\,\textbf {:}\,\hat {\nabla }\hat {T}\right \gt _{r}. \end{align}

The production and dissipation terms are calculated by implementing the eigenfunction into the expressions (3.11) and their sum shows a good agreement with the eigenvalue $\omega _{i,max }$ . For every case, the contribution from the thermal dissipation $\epsilon _{p}$ to the growth rate is small as $Pr\leqslant 1$ and the growth rate $\omega _{i}$ is mainly determined by a difference between the production $\mathcal {P}_{\Omega }$ and the viscous dissipation $\epsilon _{k}$ . For the cases with $Pr\ll 1$ , if we apply the small- $Pr$ approximation, we can re-express the thermal dissipation $\epsilon _{p}$ as

(3.12) \begin{align} \epsilon _{p}\simeq -\frac {P_{N}}{Re}\left \lt \hat {\nabla }\hat {T}_{1}\,\textbf {:}\,\hat {\nabla }\hat {T}_{1}\right \gt _{r}\sim O\left (\frac {P_{N}}{Re}\right ), \end{align}

where the dependence on the two parameters $N$ and $Pr$ is now expressed by the single parameter $P_{N}$ under the small- $Pr$ approximation. The scaling for thermal dissipation $\epsilon _{p}\sim O(P_{N}/Re)$ is based on the assumption $\hat {T}_{1}\sim O(1)$ under the small- $Pr$ approximation.

Table 1. Values of the maximum growth rates $\omega _{i,max }$ with the corresponding wavenumber $k_{d,max }$ , production and dissipation terms for various $N$ , $Pr$ and $m$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ .

3.4. Parametric investigations

Figure 4. Neutral stability curves for different $m$ for (a) $Pr=1$ and (b) $Pr=10^{-4}$ at $\mu =0$ and $\eta =0.9$ . (c) Neutral stability curves for different $Pr$ over the rescaled parameter $P_{N}$ at $\mu =0$ , $\eta =0.9$ and $m=2$ . A thick grey line denotes the neutral stability curve from the small- $Pr$ approximation (SPA).

Figures 4(a) and 4(b) show neutral stability curves for axisymmetric and non-axisymmetric cases in the parameter space $(N,Re_{i})$ at $\mu =0$ and $\eta =0.9$ for $Pr=1$ and $10^{-4}$ , respectively. For both $Pr$ values, all the neutral stability curves increasing with $N$ ascend as $m$ increases; thus the lowest curves correspond to the axisymmetric case with $m=0$ . A difference between the two Pr values is that, for $Pr=1$ , neutral stability curves for non-axisymmetric cases ( $m\gt 0$ ) stay closer to the curve for the axisymmetric case ( $m=0$ ) as $N$ increases, while for $Pr=10^{-4}$ , the curves for $m\gt 0$ stay further from the curve for $m=0$ as $N$ increases. This implies that for strongly stratified fluids with $N\gg 0$ for $Pr=1$ , the axisymmetric mode will appear at the instability onset $Re_{i}=Re_{i,c}$ but then non-axisymmetric perturbations will also become unstable immediately after the onset; thus competition between the axisymmetric and non-axisymmetric modes will occur right above $Re_{i}\gt Re_{i,c}$ . On the other hand, for low $Pr$ and above the instability onset $Re_{i}\gt Re_{i,c}$ , the axisymmetric mode will become more dominant over non-axisymmetric modes as $N$ increases. Figure 4(c) displays neutral stability curves over the rescaled $P_{N}$ for different $Pr$ at $m=2$ . As similarly observed for the axisymmetric case $m=0$ in figure 2(c), the curves for $Pr\leqslant 10^{-2}$ overlap and agree with the prediction from the small- $Pr$ approximation. It is found that the line from the small- $Pr$ approximation for $m=2$ scales as $Re_{i,c}= 134.4+9.971P_{N}$ , where $Re_{c}=134.4$ as the $Re_{i}$ -axis intercept at $P_{N}=0$ corresponds to the critical Reynolds number for unstratified case $N=0$ . The slope 9.971 for $m=2$ is higher than the slope 8.965 of the $m=0$ case. It is verified that the critical Reynolds number and the slope increase with $m$ and thus they are at the lowest for $m=0$ . This implies that for low $Pr$ , the axisymmetric mode with $m=0$ is expected to be the most unstable one for $P_{N}=N^{2}Pr\geqslant 0$ .

Figure 5. Neutral stability curves obtained from the small- $Pr$ approximation at $\mu =0$ .

Figure 5 shows neutral stability curves obtained from the small- $Pr$ approximation in the wide parameter space $(\eta ,Re_{i})$ for different values of $P_{N}$ at $\mu =0$ . The curve for $P_{N}=0$ agrees with the result in DiPrima et al. (1984). The azimuthal wavenumber of the curves corresponds to $m=0$ as the axisymmetric perturbation is found to be more unstable than non-axisymmetric perturbations. An advantage of using the small- $Pr$ approximation is that, for $Pr\ll 1$ , the stability curves depend solely on a single parameter $P_{N}=N^{2}Pr$ instead of two parameters $N$ and $Pr$ ; thus parametric investigations become simplified. As $P_{N}$ increases, the curves ascend and, in particular, the wide-gap Taylor–Couette flow with small $\eta$ is strongly stabilised. It is also shown that the curves are less sensitive to the change in $P_{N}$ when $\eta$ is close to 1 (i.e. Taylor–Couette flow with a small gap).

This section has described how linear centrifugal instability of stratified Taylor–Couette flow is affected by strong thermal diffusion. It is shown that the stratification having a stabilising effect on the instability is suppressed by strong thermal diffusion. For the cases with $Pr\ll 1$ , the dependence on $N$ and $Pr$ is simplified further by a single rescaled parameter $P_{N}=N^{2}Pr$ as derived by the small- $Pr$ approximation. In § 4, we study how the instability modes develop nonlinearly in stratified and diffusive fluids. Nonlinear dynamical behaviours such as nonlinear saturation, secondary instability of the saturated state or transition to chaotic states are examined.

4. Nonlinear development of centrifugal instability

In this section, we investigate via DNS how the centrifugal instability develops nonlinearly. Table 2 provides details of physical and numerical parameters used for the main 3-D DNS cases that are thoroughly analysed. There are more 3-D and 2-D DNS results with similar parameters in the paper; however, we focus on the analysis of these three main cases featuring nonlinear saturation, secondary instability and transition to chaotic states. These nonlinear features are thoroughly discussed in the following subsections. In the table, parameters to be noted are the axial wavenumber $k$ and the domain length $L_{z}=2\pi /k$ (i.e. one periodic length). The wavenumber $k$ in table 2 is chosen as $k=k_{c}$ , which is the critical wavenumber at the critical Reynolds number $Re_{i}=Re_{i,c}$ at the onset of primary instability. This wavenumber choice is coherent with previous studies that investigated axisymmetric Taylor vortices for $Re_{i}\gt Re_{i,c}$ by considering the characteristic wavenumber $k$ close to $k_{c}$ as $k\simeq k_{c}$ , although a better agreement with experiments can be met for the torque of Taylor-vortex flow and wavy-vortex flow at high $Re_{i}\gg R_{i,c}$ if a suitable variation of the wavenumber is taken into account (see e.g. Meyer 1966; Davey et al. 1968; DiPrima et al. 1984). Our study considers the Reynolds number $Re_{i}$ not too far from $Re_{i,c}$ (i.e. the Reynolds-number ratio $\mathcal {R}_{c}=Re_{i}/Re_{i,c}\lt 2$ ); thus the flow is not turbulent. In this case, the choices of the wavenumber $k=k_{c}$ and the corresponding domain length $L_{z}=2\pi /k$ are seemingly appropriate for comparison with other studies. We support this presumption by providing in Appendix B the validation against experiments and numerical verification of the domain length dependence. The principal azimuthal wavenumber $m$ in DNS is $m=1$ so that the entire angle $\theta =[0,2\pi ]$ with its azimuthal periodicity is considered as the azimuthal domain.

Table 2. Physical and numerical parameters for representative 3-D DNS cases.

If the Fourier representation (2.17) is taken into account, the perturbation energy $E(t)$ can be expressed as the sum: $E(t)=\sum _{j=-M}^{M}\sum _{l=-K}^{K}\tilde {E}_{jl}(t)$ , where $\tilde {E}_{jl}$ is the modal energy defined as

(4.1) \begin{align} \tilde {E}_{jl}=\pi L_{z}\int _{1}^{1/\eta }\left (|\tilde {u}_{\! jl}|^{2}+|\tilde {v}_{\! jl}|^{2}+|\tilde {w}_{\! jl}|^{2}+N^{2}|\tilde {T}_{\! jl}|^{2}\right )r\,\mathrm {d}r. \end{align}

For 3-D DNS conducted in this study, we consider a controlled initial condition as the combination of axisymmetric modes with $(j,l)=(0,\pm 1)$ and $\tilde {E}_{jl}=5\times 10^{-7}$ for each mode, which is most unstable, and other unstable non-axisymmetric modes with $(j,l)=(\pm 1,\pm 1)$ at a smaller amplitude with $\tilde {E}_{jl}=5\times 10^{-9}$ for each mode. These modes are computed from the 1-D local LSA and the non-axisymmetric modes are normalised to be out of phase, the case in which secondary instability can be promoted at the equilibrium state (i.e. axisymmetric Taylor vortices) reached by nonlinear saturation of the axisymmetric mode (see also Davey et al. 1968; Eagles 1974). Consideration of the non-axisymmetric modes with $j=\pm 1$ (i.e. $m=\pm 1$ ) is essential to allow the nonlinear energy transfer to modes with higher azimuthal wavenumber $|m|\gt 1$ and axial wavenumber $jk$ with $|j|\geqslant 1$ . For the analysis of secondary instability with the 2-D bi-global LSA presented in the following subsections, 2-D DNS are also conducted by considering the axisymmetric modes only (i.e. $M=0$ and $N_{\theta }=1$ ).

4.1. Nonlinear saturation to axisymmetric Taylor vortices

Figure 6. (a) Time evolution of the total energy $E(t)$ (black) and representative modal energy components: $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green) and $\tilde {E}_{00}$ (red) for case 1. The dashed lines denote 1-D LSA predictions on the growth of the modes with ( $m=0$ ) (blue) and $m=1$ (green). The green dotted line denotes a 2-D bi-global LSA prediction on the decay of the non-axisymmetric mode ( $m=1$ ). (b) Velocity field on the plane $(r,z)$ at $\theta =0$ and $t=500$ with contours denoting the total azimuthal velocity $U_{\theta }$ and vector plot denoting the transverse velocity field $(u_{r},u_{z})$ . (c) The total temperature profile $\varTheta (r,z)$ at $t=500$ .

Figure 6(a) shows an example of the time evolution of the total energy $E(t)$ for case 1 and a few representative examples of the modal energy $\hat {E}_{jl}$ for the axisymmetric mode with $(j,l)=(0,1)$ , the non-axisymmetric mode with $(j,l)=(1,1)$ and the mean-flow distortion with $(j,l)=(0,0)$ . It is recalled that the indices $j$ and $l$ are from the azimuthal wavenumber $m_{j}=jm$ and axial wavenumber $k_{l}=lk$ . At the initial stage, the perturbation energy is small enough and the total energy grows exponentially as $E(t)\sim \exp (2{\textrm {Im}}(\omega _{01})t)$ , where $\omega _{01}$ is the growth rate of the most unstable axisymmetric mode (i.e. $j=0$ and $l=1$ ). A good agreement between the LSA prediction and DNS is clearly shown in figure 6 $(a)$ for the growth of the $m=0$ and $m=1$ modes at early stage. Once the energy of the axisymmetric mode increases and saturates, the flow reaches an equilibrium state featuring axisymmetric Taylor vortices as shown in figure 6( $b$ ). In this saturation process, the mean-flow distortion also grows and its energy $\hat {E}_{00}$ saturates. While the energies of the axisymmetric mode and mean-flow distortion remain constant at large $t$ , the energy of the non-axisymmetric mode with $(j,l)=(1,1)$ decays as the Taylor vortices appear. It is not shown here but other modes with higher $m$ or $k$ also decay as the Taylor vortices appear. Using the Taylor vortices at saturation as a new base state, the bi-global LSA can predict the decay of non-axisymmetric modes with the index $l=1$ as we see a good agreement on the decay rate between the DNS and the bi-global LSA predictions. In figure 6( $c$ ), the total temperature $\varTheta (r,z)=T_{B}+T$ at the equilibrium state is displayed. For case 1, thermal diffusion is moderate with $Pr=1$ and thus the temperature perturbation $T$ can grow and have an amplitude comparable to $T_{B}$ at equilibrium. In this case, the total temperature $\varTheta$ no longer varies linearly but is affected and mixed by the Taylor vortices; i.e. the temperature $\varTheta$ increases (decreases) in the region where the direction of the velocity is upwards (downwards). The Taylor vortices in stratified fluids lead to baroclinicity with a radial temperature gradient; however, for case 1 with moderate $N=1$ and $Re_{i}=145$ , the overturning of the temperature does not occur and it does not become secondarily unstable. Supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.261 demonstrates the nonlinear development of the Taylor vortices and perturbation through the saturation process. It is not shown here but with the same physical parameters of case 1, different sets of numerical parameters are tested (i.e. a higher-resolution case with $(N_{r},N_{\theta },N_{z})=(120,65,65)$ and another case with a longer domain $L_{z}=8\pi /k$ with resolution $(N_{r},N_{\theta },N_{z})=(60,33,129)$ ). The results are found to be the same qualitatively with the appearance of axisymmetric Taylor vortices without secondary instability and quantitatively with insignificant differences in the energy growth curves.

4.2. Bi-global linear stability analysis and secondary instability

Figure 7. Growth-rate curves from 1-D LSA (dashed lines) and 2-D bi-global LSA (solid lines) for $(a)$ unstratified case $N=0$ with numbers denoting $m$ , $(b)$ $(N,Pr)=(1,1)$ , $(c)$ $(N,Pr)=(1,0.01)$ and $(d)$ $(N,Pr)=(1.5,0.01)$ . Black dashed lines denote the growth rate of the axisymmetric mode $m=0$ and other coloured curves in descending order (for $\mathcal {R}_{c}\lt 1$ ) denote the growth rate for $m=1$ and higher $m$ .

For the Taylor vortex flow, which is axisymmetric, we can explore its secondary instability by performing a 2-D bi-global LSA. To obtain the axisymmetric base state $\bar {\boldsymbol{Q}}(r,z)=(\bar {\boldsymbol{U}}(r,z),\bar {T}(r,z))$ for the bi-global LSA, 2-D DNS with $N_{\theta }=1$ (i.e. $M=0$ ) are conducted instead of 3-D DNS for computational efficiency. Another reason for using 2-D DNS is that the base state $\bar {\boldsymbol{Q}}$ does not become secondarily unstable but remains saturated, and thus we can use this steady and axially periodic base state in the bi-global LSA (see also Park et al. 2011). Figure 7 shows examples of the growth rates of the most unstable modes versus the Reynolds-number ratio $\mathcal {R}_{c}=Re_{i}/Re_{i,c}$ . Results are computed from the 1-D local LSA (dashed lines) and 2-D bi-global LSA (solid lines). Highly non-axisymmetric modes not appearing in figure 7 are more stable than the modes shown in figure 7. For every 1-D case in figure 7, the axisymmetric mode with $m=0$ becomes primarily unstable, as shown by black dashed lines. In both 1-D and 2-D LSA, the axial wavenumber $k=k_{c}$ of the axisymmetric mode at the onset of instability $\mathcal {R}_{c}=1$ . In the bi-global LSA, we use the Taylor vortices with axial wavenumber $k=k_{c}$ , the value at the onset of instability $Re_{i}=Re_{i,c}$ , and compute growth rates of non-axisymmetric modes with $m\gt 0$ by increasing $\mathcal {R}_{c}$ . For $\mathcal {R}_{c}\lt 1$ , growth rates of non-axisymmetric modes computed from the bi-global LSA are the same as those from the 1-D LSA since the axisymmetric Taylor vortices are not developed and the base state obtained from 2-D DNS is essentially cylindrical Couette flow (2.4). For the unstratified case with $N=0$ in figure 7(a), the axisymmetric mode becomes unstable at $Re_{i,c}\simeq 131.6$ . As the axisymmetric Taylor vortices appear for $Re_{i}\gt Re_{i,c}$ , the growth rates of non-axisymmetric modes for $m\geqslant 1$ are attenuated by this new base state as shown by coloured solid lines in figure 7(a). In the presence of the Taylor vortices, growth rates of weakly non-axisymmetric modes for $1\leqslant m\leqslant 4$ increase slowly with $\mathcal {R}_{c}$ and the non-axisymmetric mode with $m=1$ becomes unstable at $\mathcal {R}_{c}=1.075$ (i.e. $Re_{i}=Re_{i,2}\simeq 141.5$ , where $Re_{i,2}$ denotes the secondary critical Reynolds number at which a non-axisymmetric mode becomes secondarily unstable). When $\mathcal {R}_{c}=1.075$ , the growth rates of highly non-axisymmetric modes increase faster with $\mathcal {R}_{c}$ for the base flow case with Taylor vortices than the case with cylindrical Couette flow. This implies that highly non-axisymmetric modes are strongly promoted by the axisymmetric Taylor vortices for large $\mathcal {R}_{c}$ . However, the $m=1$ mode is the second unstable mode and it is difficult to predict in advance which non-axisymmetric mode becomes the next dominant mode for $Re_{i}\gt Re_{i,2}$ as nonlinear interactions involving the growth of the $m=1$ mode will lead to a new non-axisymmetric base state.

Characteristics of secondary instability are similar for a stratified case with $(N,Pr)=(1,1)$ in figure 7(b), the case where the axisymmetric mode becomes primarily unstable at a higher $Re_{i,c}\simeq 140.10$ and the non-axisymmetric mode with $m=1$ becomes secondarily unstable at a higher Reynolds-number ratio $\mathcal {R}_{c}=1.085$ (i.e. $Re_{i,2}\simeq 152.1$ ) than the unstratified case. For a highly diffusive case with $Pr=0.01$ in figure 7(c), characteristics of secondary instability change as weakly non-axisymmetric modes with $m=1$ and $2$ are stabilised by the Taylor vortices while a highly non-axisymmetric mode with $m=7$ becomes secondarily unstable at $\mathcal {R}_{c}=1.3$ (i.e. $Re_{i,2}=171.31$ ). This implies that the onset of secondary instability is delayed by strong thermal diffusion at $Pr=0.01$ and highly non-axisymmetric modes become dominant while weakly non-axisymmetric modes are suppressed. The onset of secondary instability is further delayed as stratification becomes stronger with $N=1.5$ as shown in figure 7(d), the case where the Taylor vortices become unstable by the $m=9$ mode at a higher ratio with $\mathcal {R}_{c}=1.466$ (i.e. $Re_{i,2}=193.25$ ).

The prominence of highly non-axisymmetric modes, which become dominant over weakly non-axisymmetric modes with low $m$ , can be understood by conducting an energetics analysis. Similar to (3.9), one can obtain an equation for the temporal evolution of perturbation energy from the linearised perturbation equations (A11)–(A15) as

(4.2) \begin{align} \frac {\partial \bar {E}}{\partial t}=\mathcal {P}_{\bar {U}}+\mathcal {P}_{\bar {V}}+\mathcal {P}_{\bar {W}}+\mathcal {P}_{\bar {\mathcal {T}}}-\mathcal {D}_{k}-\mathcal {D}_{p}, \end{align}

where $\mathcal {P}_{\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}}}$ are the production terms induced by axisymmetric Taylor vortices with $(\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}})$ and $\mathcal {D}_{K,P}$ are the viscous and thermal dissipation terms, respectively, all of which are defined as follows:

(4.3) \begin{eqnarray} &&\mathcal {P}_{\bar {U}}=-\left \lt \frac {\partial \bar {U}}{\partial r}\bar {u}^{2}+\frac {\bar {U}}{r}\bar {v}^{2}+\frac {\partial \bar {U}}{\partial z}\bar {u}\bar {w}\right \gt ,\quad \mathcal {P}_{\bar {V}}=-\left \lt \left (\frac {\partial \bar {V}}{\partial r}-\frac {\bar {V}}{r}\right )\bar {u}\bar {v}+\frac {\partial \bar {V}}{\partial z}\bar {v}\bar {w}\right \gt ,\nonumber \\ && \mathcal {P}_{\bar {W}}=-\left \lt \frac {\partial \bar {W}}{\partial r}\bar {u}\bar {w}+\frac {\partial \bar {W}}{\partial z}\bar {w}^{2}\right \gt ,\quad \mathcal {P}_{\bar {\mathcal {T}}}=-N^{2}\left \lt \frac {\partial \bar {\mathcal {T}}}{\partial r}\bar {u}\bar {T}+\frac {\partial \bar {\mathcal {T}}}{\partial z}\bar {w}\bar {T}\right \gt ,\nonumber \\ &&\mathcal {D}_{K}=\frac {1}{Re}\left \lt {\nabla }\bar {\textbf {u}}\,\textbf {:}\,{\nabla }\bar {\textbf {u}}\right \gt ,\quad \mathcal {D}_{P}=\frac {N^{2}}{RePr}\left \lt {\nabla }\bar {T}\,\textbf {:}\,{\nabla }\bar {T}\right \gt . \end{eqnarray}

We note that the dissipation terms $\mathcal {D}_{K,P}$ are always positive and play a stabilising role while the production terms $\mathcal {P}_{\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}}}$ are not necessarily positive as they depend on the correlation between the base state $(\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}})$ and perturbation variables $(\bar {u},\bar {v},\bar {w},\bar {T})$ . Applying the normal mode (2.29) to (4.2) and considering the normalisation $\hat {E}_{m}=1$ , where $\hat {E}_{m}$ is the modal energy defined as

(4.4) \begin{align} \hat {E}_{m}=\left \lt |\hat {u}_{m}|^{2}+|\hat {b}_{m}|^{2}+|\hat {w}_{m}|^{2}+N^{2}|\hat {T}_{m}|^{2}\right \gt _{rz},\,\left \lt \hat {\textbf {X}}_{m}\right \gt _{rz}=\int _{0}^{L_{z}}\int _{1}^{1/\eta }\hat {\textbf {X}}_{m} \,r\,\mathrm {d}r\,\mathrm {d}z, \end{align}

we obtain the following relation between the growth rate $\omega _{m,i}$ and the modal contribution terms:

(4.5) \begin{align} \omega _{m,i}=\hat {\mathcal {P}}_{m,\bar {U}}+\hat {\mathcal {P}}_{m,\bar {V}}+\hat {\mathcal {P}}_{m,\bar {W}}+\hat {\mathcal {P}}_{m,\bar {\mathcal {T}}}-\hat {\mathcal {D}}_{m,K}-\hat {\mathcal {D}}_{m,P}, \end{align}

where

(4.6) \begin{eqnarray} &&\hat {\mathcal {P}}_{m,\bar {U}}=-\left \lt \frac {\partial \bar {U}}{\partial r}|\hat {u}_{m}|^{2}+\frac {\bar {U}}{r}|\hat {v}_{m}|^{2}+\frac {\partial \bar {U}}{\partial z}\left (\frac {\hat {u}_{m}^{*}\hat {w}_{m}+\hat {u}_{m}\hat {w}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ &&\hat {\mathcal {P}}_{m,\bar {V}}=-\left \lt \left (\frac {\partial \bar {V}}{\partial r}-\frac {\bar {V}}{r}\right )\left (\frac {\hat {u}_{m}^{*}\hat {v}_{m}+\hat {u}_{m}\hat {v}^{*}_{m}}{2}\right )+\frac {\partial \bar {V}}{\partial z}\left (\frac {\hat {v}_{m}^{*}\hat {w}_{m}+\hat {v}_{m}\hat {w}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ && \hat {\mathcal {P}}_{m,\bar {W}}=-\left \lt \frac {\partial \bar {W}}{\partial r}\left (\frac {\hat {u}_{m}^{*}\hat {w}_{m}+\hat {u}_{m}\hat {w}^{*}_{m}}{2}\right )+\frac {\partial \bar {W}}{\partial z}|\hat {w}_{m}|^{2}\right \gt _{rz},\nonumber \\ && \hat {\mathcal {P}}_{m,\bar {\mathcal {T}}}=-N^{2}\left \lt \frac {\partial \bar {\mathcal {T}}}{\partial r}\left (\frac {\hat {u}_{m}^{*}\hat {T}_{m}+\hat {u}_{m}\hat {T}^{*}_{m}}{2}\right )+\frac {\partial \bar {\mathcal {T}}}{\partial z}\left (\frac {\hat {w}_{m}^{*}\hat {T}_{m}+\hat {w}_{m}\hat {T}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ &&\hat {\mathcal {D}}_{m,K}=\frac {1}{Re}\left \lt \hat {\nabla }_{m}\hat {\textbf {u}}_{m}\,\textbf {:}\,\hat {\nabla }_{m}\hat {\textbf {u}}_{m}\right \gt _{rz},\quad \hat {\mathcal {D}}_{m,P}=\frac {N^{2}}{RePr}\left \lt \hat {\nabla }_{m}\hat {T}_{m}\,\textbf {:}\,\hat {\nabla }_{m}\hat {T}_{m}\right \gt _{rz}. \end{eqnarray}

Examples of the growth rate $\omega _{i}$ and instability contribution terms (4.6) computed from the eigenfunction $\hat {\textbf {q}}_{m}$ are presented in table 3. In the table, the growth rate $\omega _{m,i}$ computed from the eigenvalue problems (2.27) and (2.30) and the sum on the right-hand side of (4.5), which is obtained by integrating the eigenfunction computed from either (2.27) for 1-D LSA cases or (2.30) for 2-D LSA cases, are in good agreement with a very small difference of $O(10^{-4})$ . For 1-D LSA cases in which the base flow is cylindrical Couette flow, the axisymmetric mode with $m=0$ is most unstable due to the largest production $\hat {\mathcal {P}}_{m,\bar {V}}$ and the least total dissipation $\hat {\mathcal {D}}_{m}=\hat {\mathcal {D}}_{m,K}+\hat {\mathcal {D}}_{m,P}$ . This applies to both cases with $Pr=0.01$ and $Pr=1$ although the thermal dissipation $\hat {\mathcal {D}}_{m,P}$ is small when $Pr$ is small. As $m$ increases, the production $\hat {\mathcal {P}}_{m,\bar {V}}$ decreases while the dissipation $\hat {\mathcal {D}}_{m}$ increases, which leads to the decrease of the growth rate. For $m=1$ , we compare the 1-D LSA cases against the 2-D LSA cases where the base flow is 2-D axisymmetric Taylor vortices. It is clearly shown that the growth rate decreases as the flow becomes 2-D with decrease in the total production term $\hat {\mathcal {P}}_{m}=\hat {\mathcal {P}}_{m,\bar {U}}+\hat {\mathcal {P}}_{m,\bar {V}}+\hat {\mathcal {P}}_{m,\bar {W}}+\hat {\mathcal {P}}_{m,\bar {T}}$ and increase in the total dissipation $\hat {\mathcal {D}}_{m}$ . While other production terms such as $\hat {\mathcal {P}}_{m,\bar {U}}$ , $\hat {\mathcal {P}}_{m,\bar {W}}$ and $\hat {\mathcal {P}}_{m,\bar {T}}$ appear for 2-D cases due to the velocity and temperature gradients of the Taylor vortices in the radial and axial directions (see also figure 6 b,c), their contribution to the growth rate is smaller than the contribution from the azimuthal velocity $\bar {V}$ (i.e. $\hat {\mathcal {P}}_{m,\bar {V}}$ ). What is noteworthy in secondary instability of axisymmetric Taylor vortices is that the total dissipation $\hat {\mathcal {D}}_{m}$ increases monotonically with $m$ while the total production $\hat {\mathcal {P}}_{m}$ increases and then decreases as $m$ increases.

Table 3. Growth rates $\omega _{m,i}$ and contribution terms computed from 1-D local LSA and 2-D bi-global LSA for $(Re_{i},N,Pr)=(200,1,0.01)$ and $(Re_{i},N,Pr)=(200,1,1)$ .

Figure 8. (a,b) Growth rate $\omega _{m,i}$ (filled circles in a), production $\hat {\mathcal {P}}_{m}$ (open circles in b) and dissipation $\hat {\mathcal {D}}_{m}$ (crosses in b) versus azimuthal wavenumber $m$ for different $Pr$ : $Pr=10^{-6}$ (black), $10^{-4}$ (blue), $0.01$ (red) and $1$ (green) at $(Re_{i},N)=(200,1)$ . (c) Production and dissipation terms versus $m$ for $(Re_{i},N,Pr)=(200,1,0.01)$ .

Figure 8(a,b) demonstrates behaviours of the growth rate $\omega _{m,i}$ , production $\hat {\mathcal {P}}_{m}$ and dissipation $\hat {\mathcal {D}}_{m}$ with varying $m$ for different $Pr$ at $(Re_{i},N)=(200,1)$ . Except for $Pr=0.01$ , the growth rate $\omega _{m,i}$ increases first and then decreases as $m$ increases (figure 8 a). This is due to the fact that the production $\hat {\mathcal {P}}_{m}$ increases only for small $m$ and decreases as $m$ increases further while the dissipation $\hat {\mathcal {D}}_{m}$ increases exponentially with $m$ (figure 8 b). This results in the growth rate $\omega _{m,i}$ being positive only in a finite range of $m$ . The growth rates for $Pr=10^{-4}$ and $Pr=10^{-6}$ are very similar to each other and are higher than $\omega _{m,i}$ for $Pr=1$ . The growth rate behaviour for $Pr=0.01$ is more complicated as $\omega _{m,i}$ is negative for $m\leqslant 5$ before it shows a similar increasing/decreasing trend. Only highly non-axisymmetric modes in the range $6\leqslant m\leqslant 12$ are secondarily unstable for $Pr=0.01$ while non-axisymmetric modes for other $Pr$ cases are unstable in broader ranges of $m$ such as $1\leqslant m\leqslant 10$ for $Pr=1$ and $1\leqslant m\leqslant 13$ for $Pr=10^{-4}$ and $10^{-6}$ . This implies that a highly non-axisymmetric flow pattern, as further discussed in the next subsection, can appear for $Pr=0.01$ by secondary instability of Taylor vortices. Figure 8(c) describes variations of all the production and dissipation terms with $m$ at $(Re_{i},N,Pr)=(200,1,0.01)$ . The dominant production and dissipation terms are $\hat {\mathcal {P}}_{\bar {V}}$ and $\hat {\mathcal {D}}_{K}$ while other terms are small. The production term $\hat {\mathcal {P}}_{\bar {W}}$ becomes negative as $m\gt 1$ increases implying the stabilising role of the axial velocity $\bar {W}$ .

Figure 9 details how the growth rate, production and dissipation vary with $Pr$ for various non-axisymmetric modes at $Re_{i}=200$ and $N=1$ . In figure 9 $(a)$ , we see that, for low $m\leqslant 7$ , the growth rate $\omega _{m,i}$ decreases as $Pr$ decreases from $Pr=1$ and then it increases as $Pr$ further decreases from $Pr= 10^{-1}$ to $10^{-2}$ . For high $m\geqslant 8$ , the growth rate increases monotonically as $Pr$ decreases. Such behaviours are similar for the total production term $\hat {P}_{m}$ while the dissipation term $\hat {D}_{m}$ decreases overall monotonically as $Pr$ decreases, as shown in figure 9 $(b)$ . As revealed in figure 8 $(c)$ , the dominant contributions to the production and dissipation come from $\hat {P}_{m,\bar {V}}$ and $\hat {D}_{m,K}$ , respectively, which have the magnitude of $O(1)$ (figure 9 $c$ ). The $Pr$ behaviour of $\hat {P}_{m,\bar {V}}$ is similar to that of the total production $\hat {P}_{m}$ , which drives the growth rate. Although it is not straightforward to delineate the $Pr$ tendency on which non-axisymmetric mode becomes most unstable, one can infer from figure 9 $(c)$ that highly non-axisymmetric modes overall contribute to the secondary instability via its interaction with the azimuthal velocity $\bar {V}$ when $Pr$ is sufficiently low as $Pr\lt 10^{-2}$ , while the behaviour is more complicated in the range $10^{-2}\lt Pr\lt 1$ . In figure 9( $d$ ), we describe other contribution terms $\hat {P}_{m,\bar {U}}$ , $\hat {P}_{m,\bar {W}}$ and $\hat {D}_{m,P}$ , which are minor with a magnitude of $O(0.1)$ or less. The production $\hat {P}_{m,\bar {U}}$ is overall positive except for $m=1$ around $Pr=10^{-2}$ and behaves similar to $\hat {P}_{m,\bar {V}}$ while the contribution $\hat {P}_{m,\bar {W}}$ is overall negative except for the case with $m=1$ , which is positive, and overall decreases as $Pr$ decreases. As can be inferred from the expression in (4.6), thermal dissipation $\hat {D}_{m,P}$ overall decreases to zero as $Pr$ decreases.

Figure 9. Variations with the Prandtl number $Pr$ for $(a)$ the growth rate $\omega _{m,i}$ , $(b)$ the total production $\hat {P}_{m}$ and dissipation $\hat {D}_{m}$ , $(c)$ the dominant production term $\hat {P}_{m,\bar {V}}$ and dissipation term $\hat {D}_{m,K}$ and $(d)$ other production terms $\hat {P}_{m,\bar {U}}$ , $\hat {P}_{m,\bar {W}}$ and dissipation term $\hat {D}_{m,P}$ for different $m=1,4,7,11$ and 14 at $Re_{i}=200$ and $N=1$ .

Figure 10. Neutral stability curves from the 1-D LSA of the axisymmetric mode $m=0$ (solid lines denoting $Re_{i,c}$ ) and bi-global LSA of the non-axisymmetric modes (dashed lines denoting $Re_{i,2}$ ) for (a) $N=1$ and (b) $Pr=0.01$ . The numbers above the dashed lines indicate the azimuthal wavenumbers of the non-axisymmetric mode which becomes secondarily unstable.

Figure 10 displays different neutral stability curves obtained from the 1-D local LSA (solid lines) and 2-D bi-global LSA (dashed lines) to show how the critical Reynolds numbers $Re_{i,c}$ and $Re_{i,2}$ vary with the Prandtl number $Pr$ or the Brunt–Väisälä frequency $N$ . For $(N,Pr)=(1,1)$ in figure 10( $a$ ), primary centrifugal instability by an axisymmetric mode occurs at $Re_{i,c}=140.2$ while secondary instability occurs at $Re_{i,2}=152.0$ by a non-axisymmetric mode with $m=1$ . At $N=1$ , the critical Reynolds number $Re_{i,c}$ decreases monotonically as $Pr$ decreases while the secondary critical Reynolds number $Re_{i,2}$ increases as $Pr$ decreases from $Pr=1$ to $Pr=0.06$ and then $Re_{i,2}$ decreases monotonically from the peak at $Pr=0.06$ as $Pr$ further decreases. It is noteworthy that, in the range $10^{-3}\lt Pr\lt 1$ , the secondary instability is triggered by highly non-axisymmetric modes with $m\geqslant 4$ ; thus for $Re_{i}\gt Re_{i,2}$ , we expect to observe highly non-axisymmetric flow patterns in this range of $Pr$ . This highly non-axisymmetric pattern is further discussed in the following subsection. In the range $Pr\leqslant 10^{-3}$ , the critical Reynolds numbers $Re_{i,c}$ and $Re_{i,2}$ approach those of the unstratified case at $N=0$ (i.e. $Re_{i,c}=131.6$ and $Re_{i,2}=141.5$ ), the case with secondary instability triggered by a non-axisymmetric mode with $m=1$ . At $Pr=0.01$ , it is shown in figure 10 $(b)$ that the critical Reynolds number $Re_{i,c}$ does not change significantly in the range $0\leqslant N\leqslant 2$ (e.g. $Re_{i,c}=131.6$ at $N=0$ and $Re_{i,c}=132.0$ at $N=2$ ) while the secondary critical Reynolds number $Re_{i,2}$ increases monotonically with $N$ . The corresponding azimuthal wavenumber of the non-axisymmetric mode, which triggers secondary instability, also increases with $N$ .

Neutral stability curves in figure 10 delineate the regimes of axisymmetric and non-axisymmetric flow patterns in the parameter space $(Re_{i},Pr)$ or $(Re_{i},N)$ . The results are obtained by 1-D and 2-D LSA and we discuss in the next subsection about nonlinear simulation results how flow patterns change from axisymmetric Taylor vortices to non-axisymmetric wavy vortices via secondary instability or how the flow transitions to a chaotic and irregular state, a precursor stage prior to turbulence.

4.3. Transition of flow patterns

Figure 11. (a) Temporal evolution of the total energy $E(t)$ (black) and modal energy components $\tilde {E}_{jl}$ for case 2: $\tilde {E}_{00}$ (red), $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green), $\tilde {E}_{91}$ (yellow), and other energy components denoted by grey lines. (b) The corresponding instantaneous velocity profiles at different times on the plane $(r,z)$ at $\theta =0$ . The contours denote the azimuthal velocity $U_{\theta }(r,z)$ and the vector plot denotes the transverse velocity field $(U_{r},U_{z})$ .

Figure 12. Time evolution of the vertical velocity $u_{z}$ in the $(x,y)$ plane at $z=0$ for case 2.

We now consider case 2 at $(Re_{i},N,Pr)=(200,1,0.01)$ , a case where the axisymmetric Taylor vortices become secondarily unstable by highly non-axisymmetric modes (see also figures 7 $c$ and 10 $a$ ). The perturbation energy plot in figure 11 $(a)$ shows that the axisymmetric mode becomes unstable and saturates quickly to the state of axisymmetric Taylor vortices (e.g. figure 11 $b$ at $t=100$ ). At this first saturation state, the energies of weakly non-axisymmetric modes like the $m=1$ mode (green line in figure 11 $a$ ) decay exponentially. In contrast, highly non-axisymmetric modes like the $m=9$ mode (yellow line in figure 11 $a$ ) become unstable as the axisymmetric Taylor vortices sustain and their amplitudes become comparable to that of the axisymmetric mode as $t\gt 200$ . The flow reaches a second saturation state in which the Taylor vortices start to oscillate vertically (figure 11 $b$ and supplementary movie 2) and become non-axisymmetric. The transition at which non-axisymmetric modes become dominant is clearly illustrated by the vertical velocity $u_{z}(x,y)$ plot for different time $t$ (figure 12 and supplementary movie 3). The most energetic non-axisymmetric mode with $m=1$ loses its energy at the first transition stage (e.g. $u_{z}$ at $t=100$ ) while other non-axisymmetric modes with higher $m$ gain their energies and become comparable in the transition stage $150\lt t\lt 200$ . Weakly and highly non-axisymmetric modes compete during the transition (e.g. $u_{z}$ at $t=150$ and $190$ ) and after $t\gt 200$ , the non-axisymmetric mode with $m=9$ becomes the second dominant mode after the axisymmetric mode $m=0$ forming the wavy Taylor vortices. This non-axisymmetric pattern in $u_{z}$ has a maximum amplitude around the centreline between the two cylinders and is different from the feature of SRI, which is triggered by two out-of-phase inertia–gravity waves confined near the two cylinders and has maxima near the cylinders (Park & Billant 2013; Park et al. 2017). Contours of the azimuthal vorticity $\omega _{\theta }=\partial u_{r}/\partial z-\partial u_{z}/\partial r$ in figure 13 $(a)$ clearly show the axisymmetric Taylor vortices at $t=100$ as the first saturation state and non-axisymmetric wavy Taylor vortices as the second saturation state. Various modal energies of the perturbation are expected to either saturate or decay after $t\gt 500$ and thus non-axisymmetric Taylor vortices sustain for large $t$ .

Figure 13. $(a)$ Contours of the azimuthal vorticity $\omega _{\theta }$ at the boundary surfaces for case 2. $(b)$ Profiles of the averaged velocity $\bar {U}_{\theta }(r)$ and cylindrical Couette flow $V_{B}(r)$ for case 2. $(c)$ Nusselt number $Nu$ as a function of time $t$ for case 2 (black) and case 3 (blue).

Figure 13 $(b)$ shows profiles of the total averaged azimuthal velocity $\bar {U}_{\theta }(r,t)=V_{B}(r)+\tilde {v}_{00}(r,t)$ , the latter which is the azimuthal velocity of the mean-flow distortion denoting the azimuthal perturbation velocity averaged in the directions $\theta$ and $z$ . Compared with the profile of $\bar {U}_{\theta }$ at $t=100$ for the axisymmetric Taylor vortices at the first saturation, the distortion in the azimuthal velocity is reduced for non-axisymmetric wavy Taylor vortices at the second saturation. The reduced distortion in mean flow implies the reduction in the velocity gradient at the inner cylinder and the corresponding torque applied at the inner cylinder. Following Martínez Arias (2015), we define an axially averaged non-dimensional torque $G$ at the inner cylinder as

(4.7) \begin{align} G=\left (\frac {\eta Re_{i}}{1-\eta }\right )\frac {1}{2\pi L_{z}}\int _{0}^{2\pi }\int _{0}^{L_{z}}\left .\left (\frac {U_{\theta }}{r}-\frac {\partial U_{\theta }}{\partial r}\right )\right |_{r=1}\mathrm {d}\theta \,\mathrm {d}z,\, \end{align}

which gives the laminar torque $G_{\mathrm {lam}}$ for cylindrical Couette flow $U_{\theta }=V_{B}(r)$ as

(4.8) \begin{align} G_{\mathrm {lam}}=\frac {2\eta (1-\mu )Re_{i}}{(1+\eta )(1-\eta )^{2}}. \end{align}

We also define the Nusselt number $Nu$ as $Nu=G/G_{\mathrm {lam}}$ , which is the ratio between the transverse convective transport of angular velocity and the molecular transport of angular velocity as a measure of the non-dimensional angular momentum transfer (Dubrulle & Hersant 2002; Eckhardt et al. 2007; Martínez Arias 2015). Figure 13 $(c)$ shows how the Nusselt number $Nu$ changes over time for case 2 (black line). Compared with the axisymmetric Taylor vortices in the range $50\leqslant t\leqslant 150$ having a constant $Nu=1.77$ during the first saturation process, the wavy Taylor vortices at the second saturation in the range $t\gt 200$ have a lower $Nu=1.68$ . This implies that the non-axisymmetric Taylor vortices after secondary instability lead to reduced angular momentum transfer compared with the axisymmetric Taylor vortices after primary centrifugal instability.

Figure 14. (a) Temporal evolution of the total energy $E(t)$ (black) and modal energy components $\tilde {E}_{jl}$ for case 3: $\tilde {E}_{00}$ (red), $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green), $\tilde {E}_{41}$ (brown), and other energy components denoted by grey lines. Corresponding instantaneous velocity field $u_{z}(x,y)$ at $z=0$ at (b) $t=200$ and (c) $t=500$ .

Figure 15. (a) Evolution of the total temperature $\varTheta (r,z)$ at $\theta =0$ for case 3. (b) Corresponding spatio-temporal diagram of $\varTheta (t,z)$ measured at $r=1.05$ and $\theta =0$ .

Case 3 with $(N,Pr)=(1,1)$ considers the Reynolds number $Re_{i}=200$ at the ratio $\mathcal {R}_{2}=Re_{i}/Re_{i,2}=1.316$ , which is larger than the ratio $\mathcal {R}_{2}=1.168$ in case 2. In this case, it is observed that the secondary instability bypasses the saturation but leads to an irregular fluctuation as shown by the temporal variation of the Nusselt number $Nu$ in figure 13 $(c)$ and modal energies of perturbation in figure 14 $(a)$ . Primary centrifugal instability breaks down as $t\gt 100$ where other non-axisymmetric modes including the mode with $m=1$ become dominant with strong nonlinear modal interaction. For case 3, the most energetic mode varies over time and the flow exhibits a rather chaotic temporal fluctuation as shown by the energy evolution in figure 14 $(a)$ and by the velocity and temperature fields in the $(r,z)$ plane in supplementary movies 4 and 5. Although this irregular flow state does not show a small length-scale structure (see also figure 14 $b$ and supplementary movie 6), the levels of the modal energies of large-wavenumber modes are not negligible and turbulence is expected to occur for $(N,Pr)=(1,1)$ at high Reynolds number $Re_{i}\gt 200$ , the range to be further explored in a future study. The Prandtl number $Pr=1$ is not small and the temperature perturbation $T$ is comparable with the base temperature $T_{B}$ . Therefore, the fluctuation in the total temperature $\varTheta =T_{B}+T$ induced by the velocity fluctuation is visible as shown in figure 15 $(a)$ although the overturning of temperature, which is a potential source of baroclinic instability, does not occur. Chaotic fluctuation of the temperature $\varTheta$ after $t\gt 100$ is also clearly shown in the spatio-temporal diagram measured near the mid-point at $r=1.05$ (figure 15 $b$ ).

Figure 16. $(a)$ Nusselt number $Nu$ versus Reynolds number $Re_{i}$ for $N=0$ (black), $(N,Pr)=(1,1)$ (blue), $(N,Pr)=(1,0.01)$ (red) and $(N,Pr)=(1,10^{-4})$ (green). Filled circles and squares indicate $Nu$ for the axisymmetric Taylor vortices and non-axisymmetric wavy vortices at saturation, respectively. Open squares with error bars indicate statistically averaged $Nu$ and the minimum and maximum of $Nu$ in the time interval considered (see also figure 13 $c$ ). $(b)$ Nusselt number $Nu$ same as $(a)$ but over the rescaled Reynolds-number ratio $\mathcal {R}_{c}$ . The grey line indicates the Nusselt number from (4.9). $(c)$ Temporal evolution of the total energy $E(t)$ for different $Pr$ at $Re_{i}=200$ and $N=1$ . Black line denotes the unstratified case with $N=0$ .

Figure 16 $(a)$ shows how the Nusselt number $Nu$ changes with the Reynolds number $Re_{i}$ for different sets of $(N,Pr)$ : $(1,1)$ (blue), $(1,0.01)$ (red), $(1,10^{-4})$ (green) and the unstratified case with $N=0$ (black). The Nusselt number $Nu$ is computed at the saturation state of the axisymmetric Taylor vortices (filled circles), which are computed through 3-D DNS with a numerical resolution the same as that in case 1, or for non-axisymmetric wavy vortices at the second saturation (filled squares) computed by 3-D DNS with a numerical resolution the same as that in case 2. For the flow in which its instability does not saturate and $Nu$ fluctuates like case 3, open squares with error bars are used to indicate the mean, minimum and maximum $Nu$ averaged in the time interval of fluctuation (see e.g. figure 13 $c$ ). In these cases with fluctuations, the same numerical resolution as in case 3 is used. It is clearly shown that $Nu$ for cases with the axisymmetric Taylor vortices increases fast with $Re_{i}$ and it increases slowly down as secondary instability occurs. The increasing trend of the Nusselt number $Nu$ is similar for unstratified cases with $N=0$ and stratified cases with $(N,Pr)=(1,10^{-4})$ . For both sets, secondary instability leading to a fluctuating flow state appears for $Re_{i}\gt 160$ and there is a jump in $Nu$ around $Re_{i}=180$ . The Nusselt number $Nu$ continues to increase similarly for both as $Re_{i}$ increases further for $Re_{i}\gt 180$ . For the cases with $(N,Pr)=(1,0.01)$ , the Nusselt number $Nu$ is higher than for other cases as secondary instability occurs at a later stage as $Re_{i,2}=171.31$ for $Pr=0.01$ . Up to $Re_{i}=210$ , the breakdown of secondary instability is not observed for $Pr=0.01$ . Unlike other cases with critical Reynolds numbers $Re_{i,c}=131.6{-}131.7$ , the centrifugal instability for $(N,Pr)=(1,1)$ occurs at a higher $Re_{i,c}=140.2$ and thus $Nu$ in this case is smaller than for other cases at the same $Re_{i}$ . Although the onset of primary centrifugal instability is delayed for $(N,Pr)=(1,1)$ , it is found that the breakdown of instability occurs early and the fluctuation in $Nu$ is larger than in other cases.

For a better comparison among the cases, we display in figure 16 $(b)$ the same $Nu$ against the Reynolds-number ratio $\mathcal {R}_{c}=Re_{i}/Re_{i,c}$ . For the axisymmetric Taylor vortices, DiPrima et al. (1984) proposed the following scaling law for $Nu$ as a function of the ratio $\mathcal {R}_{c}$ :

(4.9) \begin{align} Nu=1+A\left (1-\frac {1}{\mathcal {R}_{c}^{2}}\right )+B\left (1-\frac {1}{\mathcal {R}_{c}^{2}}\right )^{2}, \end{align}

where $A=1.246$ and $B=-0.37$ are the constants for $\eta =0.9$ according to DiPrima et al. (1984). It is clearly shown in figure 16 $(b)$ that $Nu$ for every case agrees well with the scaling law (4.9) in the range $\mathcal {R}_{c}\leqslant 1.1$ where the axisymmetric Taylor vortices are present. This implies that $Nu$ for the Taylor vortices, which are induced by primary centrifugal instability, is not strongly influenced by $Pr$ . When the flow is secondarily unstable or chaotic, the cases with $(N,Pr)=(1,0.01)$ for which secondary instability appears late as $\mathcal {R}_{c}\gt 1.3$ have higher $Nu$ than the scaling law (4.9) while other cases for which secondary instability appears early as $\mathcal {R}_{c}\gt 1.1$ have lower $Nu$ than the proposed scaling (4.9). This implies that the angular momentum transport represented by $Nu$ depends on $Pr$ once the flow becomes secondarily unstable or chaotic. However, as figure 16(a,b) only considers $N=1$ , more detailed investigations are required to confirm the relation between $Nu$ and $Pr$ in a wider parameter space of $(N,Re_{i})$ , an important topic to be further studied in the future to unravel the physics of turbulent angular momentum transport that depends on stratification and thermal diffusion.

Figure 17. Instantaneous velocity contours for $U_{\theta }$ and vector plots for the transverse velocity field $(U_{r},U_{z})$ on the plane $(r,z)$ for different $Pr$ at $t=500$ , $\theta =0$ , $Re_{i}=200$ and $N=1$ . The rightmost panel corresponds to the unstratified case with $N=0$ .

As a summary of the nonlinear dynamics with varying $Pr$ , we display in figures 16 $(c)$ and 17 the temporal evolution of the total energy $E(t)$ and instantaneous velocity fields for different $Pr$ at $Re_{i}=200$ and $N=1$ . We see in figure 16 $(c)$ that, as $Pr$ is small, $Pr\lt 10^{-4}$ , the energy curves collapse and become identical to the curve of the unstratified case with $N=0$ . Instantaneous velocity fields obtained at $t=500$ from various DNS, which started from the same configuration of the initial condition described in the early part of § 4, are also identical between the stratified cases with $Pr=10^{-5}$ and $10^{-6}$ and unstratified case with $N=0$ . These results support the explanation in the introduction and figure 1 that the stratification effect is suppressed by strong thermal diffusion and the flow behaves as unstratified flow in both linear and nonlinear regimes. Nonetheless, there still remains a question as to whether our conclusion will remain valid for turbulent cases and if so, up to which Reynolds number $Re$ . Although DNS investigations like those by Prat & Lignières (2013) verified this question for $Re$ of $O(100)$ , further investigations should be conducted at higher Reynolds numbers.

5. Conclusion and discussion

In this paper, linear and nonlinear dynamics of centrifugal instability of stratified Taylor–Couette flow are studied for diffusive fluids in which the Prandtl number $Pr$ is low, $Pr\leqslant 1$ . The LSA reveals that the stratification effect, which suppresses the centrifugal instability (Boubnov et al. 1995; Caton et al. 2000), is suppressed by strong thermal diffusion. This implies that the thermal diffusion promotes centrifugal instability of stratified Taylor–Couette flow, a situation which can be found similarly in other contexts like shear instabilities in stratified flows in the low- $Pr$ limit (Lignières 1999; Garaud 2020a ). When the thermal diffusion is sufficiently strong, we apply the small- $Pr$ approximation and find that the instability characteristics are self-similar with the new scaled parameter $P_{N}=N^{2}Pr$ (see also Lignières 1999; Park et al. 2020). The thermal diffusion effect as promoting centrifugal instability is found to be the same for both axisymmetric and non-axisymmetric perturbations. For the parameters considered here (i.e. $\mu =0$ , $\eta =0.9$ and $Pr\leqslant 1$ ), it is found that the axisymmetric perturbation with $m=0$ is most unstable compared with non-axisymmetric perturbations. By conducting DNS for the Reynolds number $Re_{i}$ above the critical one, $Re_{i}\gt Re_{i,c}$ , we also study the nonlinear evolution of centrifugal instability. In the DNS, we use controlled initial conditions with both the axisymmetric perturbation and a smaller-amplitude non-axisymmetric perturbation to see how they grow and interact nonlinearly and to understand their nonlinear modal interaction leading to different states. At the initial stage, the axisymmetric perturbation grows fastest and saturates along the nonlinear interaction with base flow. Depending on $Re_{i}$ , a new mean flow in the shape of axisymmetric Taylor vortices can become secondarily unstable by the growth of non-axisymmetric perturbation. Unlike the case of primary instability, which is promoted by strong thermal diffusion, the 2-D bi-global LSA of the Taylor vortices reveals that thermal diffusion in the range $10^{-3}\lt Pr\lt 1$ at $N=1$ delays the onset of secondary instability and potentially laminar–turbulent transition triggered by highly non-axisymmetric perturbations. We note that the appearance of such a non-axisymmetric flow pattern by the onset of secondary instability in highly diffusive and stratified flows has not been explored properly in previous studies. At $Pr=1$ , we also examine irregular flow patterns leading to chaotic mixing of temperature, the state as a precursor to turbulence. Furthermore, we analyse how the Nusselt number $Nu$ as a measure of angular momentum transfer varies with the Reynolds number $Re_{i}$ for various sets of $(N,Pr)$ . It is verified that the secondary instability slows down the increase of $Nu$ , implying that the momentum transfer is slowly enhanced after the onset of secondary instability.

In recent years, stratified turbulence in the low- $Pr$ regime has been studied increasingly not only in fluid dynamics but also in astrophysics due to its relevance to the interior of stars and the Sun where the Prandtl number $Pr$ is low as $10^{-6}$ or below (Garaud 2020a ; Dymott et al. 2023; Garaud et al. 2024b ). Although the Reynolds numbers of astrophysical flows cannot be reached due to their immensely large length scales, simulations at lower Reynolds numbers and theories unveil new regimes of turbulence and give us a hint as to how stratified turbulence behaves with strong thermal diffusion (Cope et al. 2020; Skoutnev 2023; Garaud et al. 2024a ). What has not been explored yet though in these turbulence studies is how the turbulence is generated via laminar–turbulent transition processes. In this regard, the current study contributes to our knowledge of how the centrifugal instability develops nonlinearly towards turbulence in stratified and highly diffusive fluids. Beyond this paper, we will aim to study first the characteristics of fully developed stratified-rotating turbulence triggered by centrifugal instability in the low- $Pr$ regime. At low $Pr$ , an unexpected highly non-axisymmetric pattern can generate turbulent flow with characteristics different from those of unstratified turbulence. Another topic to be investigated in the future is how the SRI will behave in the presence of strong thermal diffusion. From the perspective of linear instability, the thermal diffusion suppressing the stratification effect is expected to suppress the SRI as well. Linear and nonlinear SRIs need to be further explored in the low- $Pr$ regime. The nonlinear development of the SRI is also of great interest in geophysical and astrophysical fluid dynamics, in particular for flows in the Keplerian or super-rotation regimes, the latter in which the outer part rotates faster than the inner region (Dubrulle et al. 2005; Le Dizès & Riedinger 2010; Park & Billant 2013).