1. Introduction
Accretion disks are thin reservoirs of dispersed material (e.g. gas and dust) orbiting massive central objects, such as growing galaxies, planets, stars and black holes (Ji & Balbus Reference Ji and Balbus2013). Their dynamics has been traditionally modelled by a type of Taylor–Couette (TC) flow, where the inner cylinder rotates faster than the outer one, but the angular momentum increases radially outward, i.e. 
$\omega _i^*>\omega _o^*$ and 
$\omega _i^* r_i^{*2}<\omega _o^* r_o^{*2}$ (Dubrulle et al. Reference Dubrulle, Dauchot, Daviaud, Longaretti, Richard and Zahn2005; Avila Reference Avila2012; Ji & Balbus Reference Ji and Balbus2013). This type of TC flow is known as quasi-Keplerian flow, as depicted in figure 1. Based on the TC flow model, the current study aims to better elucidate the angular momentum transport mechanism in accretion disks.
(a) Sketch of the coaxial co-rotating TC flow system investigated in the present work. The gravitational acceleration points radially inwards. (b) Flow regimes in the parameter space. The shaded area corresponds to quasi-Keplerian flows that include the Keplerian case satisfying 
$\omega _i^* r_i^{*{3}/{2}}=\omega _o^* r_o^{*{3}/{2}}$. It should be noted that the Keplerian flow profile can only be approximated but not exactly satisfied by the radial structure of a TC flow even in the Keplerian case. Thus, the flow at 
$\omega _i^* r_i^{*{3}/{2}}=\omega _o^* r_o^{*{3}/{2}}$ is also quasi-Keplerian.
For the central object within a disk to continuously gain mass from the surrounding material, the ‘fluid particles’ must lose angular momentum to fall inward under gravity. Observed accretion rates cannot be fully explained by angular momentum transport due to the molecular viscosity alone (Pringle Reference Pringle1981; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016), prompting researchers to explore alternative mechanisms. Shakura & Sunyaev (Reference Shakura and Sunyaev1973) first proposed that turbulence could play a crucial role in this process. Since then, considerable effort has been devoted to searching for flow instability mechanisms that could lead to turbulence in quasi-Keplerian flows. Nevertheless, hydrodynamic considerations, excluding other magnetic or thermal effects, indicate that quasi-Keplerian flows are linearly stable to axisymmetric disturbances, according to Rayleigh's criterion (Rayleigh Reference Rayleigh1917). Recent calculations for non-axisymmetric disturbances also suggest linear stability in astrophysically interesting regimes (Deguchi Reference Deguchi2017). This hydrodynamic linear stability raises the question of how turbulence arises to facilitate the radial transport of angular momentum in accretion disks.
1.1. Magnetorotational instability and hydrodynamic subcritical nonlinear instability
A significant milestone in addressing the angular momentum transport problem was reached by Balbus & Hawley (Reference Balbus and Hawley1991), who introduced the magnetorotational instability (MRI) concept, originally proposed by Velikhov (Reference Velikhov1959), into the astrophysics community. They discovered that a quasi-Keplerian flow could be linearly unstable even in the presence of very weak magnetic fields. Over 30 years later, MRI was successfully detected in a laboratory TC flow experiment (Wang et al. Reference Wang, Gilson, Ebrahimi, Goodman, Caspary, Winarto and Ji2022b,Reference Wang, Gilson, Ebrahimi, Goodman and Jic). This success supports using quasi-Keplerian flow at laboratory scales to study the physical processes in accretion disks on astronomical scales. However, the MRI diagram struggles to explain the angular momentum transport in protoplanetary or protostellar disks. This is because in these large, cool and dusty disks, the ionisation level in certain zones (dubbed ‘dead zones’) is so low that magnetising the fluid particles is difficult or impossible, rendering this diagram inapplicable (Gammie Reference Gammie1996; Fleming & Stone Reference Fleming and Stone2003; Armitage Reference Armitage2011; Held & Latter Reference Held and Latter2018).
With the magnetic field being inoperative, researchers have explored hydrodynamic instabilities with other effects to account for radial angular momentum transport in dead zones. One line of research has focused on subcritical nonlinear instability, as linear stability does not necessarily ensure nonlinear stability (Grossmann Reference Grossmann2000; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). Unfortunately, numerical simulations and laboratory experiments have shown no signs of subcritical transition in quasi-Keplerian flows, even at Reynolds numbers as high as millions (Ji et al. Reference Ji, Burin, Schartman and Goodman2006; Ji & Balbus Reference Ji and Balbus2013; Edlund & Ji Reference Edlund and Ji2014; Ostilla-Mónico et al. Reference Ostilla-Mónico, Verzicco, Grossmann and Lohse2014; Shi et al. Reference Shi, Hof, Rampp and Avila2017). Given that typical Reynolds numbers in accretion disks can reach 
$10^{12}\sim 10^{14}$ (Ji et al. Reference Ji, Burin, Schartman and Goodman2006; Fromang & Lesur Reference Fromang and Lesur2019), much higher than currently explored values, we cannot entirely dismiss the possibility of subcritical transition therein. Another line of research has examined the thermal effects on flow stability/instability in accretion disks, which is the direction the present work follows, to be reviewed below.
1.2. Hydrodynamic instability driven by thermal effects
The thermal structure of accretion disks is complex and non-uniform, leading to intricate thermodynamics that shapes the disks’ structure (Lesur et al. Reference Lesur2022). For instance, in protoplanetary disks, numerical modelling and calculations based on astrophysical observations (Chiang & Goldreich Reference Chiang and Goldreich1997; Fromang & Lesur Reference Fromang and Lesur2019) showed that the mid-plane is relatively cool (due to radiation to space) and sandwiched between two hot surface layers (due to irradiation by the central star). These studies also demonstrated that the mid-plane temperature decreases with increasing distance from the central star, indicating radial thermal stratification.
Considering the non-uniform thermal structure, various instabilities due to the thermal effect have been identified in the modelling and analyses of accretion disks. These include the two-dimensional (2-D) subcritical baroclinic instability in the 
$r$–
$\theta$ plane due to non-axisymmetric disturbances (Lesur & Papaloizou Reference Lesur and Papaloizou2010), the three-dimensional (3-D) linear convective overstability due to axisymmetric disturbances (Klahr & Hubbard Reference Klahr and Hubbard2014) and the axisymmetric 2-D linear vertical shear instability in the 
$r$–
$z$ plane (Nelson, Gressel & Umurhan Reference Nelson, Gressel and Umurhan2013). For a comprehensive review of thermally related instability mechanisms and others, including gravitational instability, zombie vortex instability and stratorotational instability, the reader is referred to Fromang & Lesur (Reference Fromang and Lesur2019). The richness of these instabilities demonstrates that many physical effects, additional to the magnetic effects, can destabilise the canonical Keplerian flow. Nevertheless, their effectiveness depends on the disk structures and the thermal time scales of the gases (Lesur et al. Reference Lesur2022). Among these instabilities, the baroclinic instability, caused by radial buoyancy forces, is a nonlinear mechanism. The linear convective overstability has been analysed for inviscid flows within a local planar shearing box. Following the shearing box approximation adopted by Klahr & Hubbard (Reference Klahr and Hubbard2014), Lyra (Reference Lyra2014) confirmed the linear convective overstability and highlighted the need to investigate its operation in ‘global models’. To our knowledge, 3-D ‘global’ linear instability concerning the thermodynamics of accretion disks, especially in the radial direction, has rarely been reported, which is the focus of the present work. By ‘global’, we specifically mean that the boundary conditions at the two rotating cylinders are directly incorporated in the problem formulation based on the quasi-Keplerian flow model in a TC geometry. It should be mentioned that the abovementioned stratorotational instability (Yavneh, Mcwilliams & Molemaker Reference Yavneh, Mcwilliams and Molemaker2001) is a linear instability identified in TC flow, but it requires stable vertical stratification and the separation between cylinder walls must be sufficiently close for the instability to be effective.
1.3. The position and structure of the current work
In light of the efforts to uncover possible hydrodynamic instabilities without magnetic fields, this study incorporates radial thermal stratification along with stellar gravity in the modelling of accretion disks using TC flow with hot inner and cold outer walls in the quasi-Keplerian regime. The objective is to identify possible 3-D thermo-hydrodynamic linear instability mechanisms that may operate in the thermal TC flow. In our flow system, the thermal gradient aligns with the direction of stellar gravity.
While it is well known that thermal buoyancy is destabilising and quasi-Keplerian shear is strongly stabilising (Ji et al. Reference Ji, Burin, Schartman and Goodman2006), to the best of our knowledge, the quantitative aspects of their interplay in the considered flow have not been thoroughly detailed in the literature. Similar flow configurations exist but differ from our system, such as rotating Rayleigh–Bénard convection in a closed cylindrical container (Ecke & Shishkina Reference Ecke and Shishkina2023) and centrifugal TC flows with hot outer and cold inner cylinders (Becker & Kaye Reference Becker and Kaye1962; Ali & Weidman Reference Ali and Weidman1990; Meyer, Yoshikawa & Mutabazi Reference Meyer, Yoshikawa and Mutabazi2015; Kang et al. Reference Kang, Meyer, Mutabazi and Yoshikawa2017; Jiang et al. Reference Jiang, Zhu, Wang, Huisman and Sun2020; Meyer, Mutabazi & Yoshikawa Reference Meyer, Mutabazi and Yoshikawa2021; Jiang et al. Reference Jiang, Wang, Liu and Sun2022; Wang et al. Reference Wang, Jiang, Liu, Zhu and Sun2022a), where gravitational effects are ignored. To illustrate these differences, we compare our flow configuration with a representative work (Meyer et al. Reference Meyer, Mutabazi and Yoshikawa2021) in table 2 in Appendix A. Additionally, our 3-D ‘global’ thermally driven linear instability also contrasts with the local shearing box approximations and short wave approximations used in prior works such as Klahr & Hubbard (Reference Klahr and Hubbard2014), Lyra (Reference Lyra2014), Latter (Reference Latter2016) and Held & Latter (Reference Held and Latter2018). A comparison of our problem setting with Klahr & Hubbard (Reference Klahr and Hubbard2014) is provided in table 3 in Appendix A.
The rest of the paper is organised as follows. Section 2 describes the flow configuration, governing equations and linear stability analysis as an eigenvalue problem. The numerical method is introduced in § 3. Section 4 presents and discusses the results, characterising the linear instability at moderate values of the parameters (including Taylor number, Richardson number, Prandtl number and radius ratio) to examine their effects on the linear instability and the scaling relation. We also demonstrate that the linear instability persists within extreme parameter ranges relevant to accretion disks, even with weak thermal buoyancy effects, highlighting the effectiveness of this thermal instability. Comparisons to previous works are made where possible. Section 5 concludes the paper. Appendices summarise the comparison to previous works and the influence of the gravitational acceleration profile on the linear instability.
2. Problem formulation
2.1. Flow configuration and governing equations
The flow investigated is an incompressible TC flow confined between two coaxial co-rotating cylinders, differentially heated at the walls, as depicted in figure 1(a). The fluid is characterised by density 
$\rho ^*$, dynamic viscosity coefficient 
$\mu ^*$ (or kinematic viscosity 
$\nu ^*=\mu ^*/\rho ^*$) and thermal diffusivity coefficient 
$\kappa ^*$. Throughout this paper, dimensional quantities are indicated with the superscript 
$*$. For the mathematical modelling of the TC structure, a cylindrical coordinate system 
$(r,\theta,z)$ is adopted, where 
$r$, 
$\theta$ and 
$z$ represent the radial, azimuthal and axial directions, respectively. The radius of the inner cylinder is 
$r_i^*$ and that of the outer cylinder is 
$r_o^*$. The cylinders are rotating about the 
$z$ axis in the same direction with angular frequencies 
$\omega _i^*$ and 
$\omega _o^*$, respectively. The flow in the 
$z$ direction is assumed to be homogeneous and infinitely long, allowing the use of periodic boundary conditions in our linear stability analysis.
The wall temperature of the two cylinders is maintained at fixed values: 
$T_i^*$ on the inner wall and 
$T_o^*$ on the outer wall, with 
$T_i^* > T_o^*$, to mimic the temperature distribution in accretion disks (Chiang & Goldreich Reference Chiang and Goldreich1997). To account for the thermally driven buoyancy effect due to thermal stratification, stellar gravity in the radial direction is incorporated into the modelling process. Following Balbus & Hawley (Reference Balbus and Hawley1991) and Latter (Reference Latter2016), we adopt the Boussinesq approximation, assuming that the fluid density varies linearly with temperature as 
$\rho ^* = \rho _o^* - \rho _o^* \alpha ^* (T^* - T_o^*)$ and this variation is significant only in the gravitational buoyancy term. Here, 
$\alpha ^*$ is the coefficient of thermal expansion and 
$\rho _o^*$ is the reference density at temperature 
$T_o^*$. In contrast to Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021), centrifugal buoyancy is ignored in our work; the effect would be stabilising in our temperature setting.
The governing equations consist of the Navier–Stokes equation under the Boussinesq approximation, the continuity equation and the energy equation expressed in terms of temperature. For non-dimensionalisation, the gap between the two cylinders 
$d^*=r_o^* - r_i^*$ is chosen as the characteristic length scale, the characteristic velocity scale is the tangential speed of the inner wall 
$U^*=|\omega _i^* - \omega _o^*|r_i^*$ observed in the rotating frame of reference, the time scale is 
$d^*/U^*$, and the reference density, pressure and gravitational acceleration magnitude are 
$\rho _o^*$, 
$\rho _o^* U^{*2}$ and 
$g_o^*$, respectively.
The resultant non-dimensional equations, expressed in a reference frame rotating with the outer cylinder, are given as
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} ={-} \boldsymbol{\nabla}p + \frac{f(\eta)}{{Ta}^{1/2}} \nabla^2 \boldsymbol{u} + \frac{Ra f^2(\eta)}{Pr Ta}T g(r) \boldsymbol{e}_r - \frac{1}{Ro} \boldsymbol{e}_z \times \boldsymbol{u}, \end{gather}
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather} \frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} T = \frac{f(\eta)}{Pr {Ta}^{1/2}} \nabla^2 T, \end{gather}
where 
$f(\eta )={(1+\eta )^3}/{8\eta ^2}$. The form of 
$f(\eta )$ and the introduction of 
$Ta$ and 
$Ro$ (defined below) follow the derivations for isothermal flows in TC geometry in Grossmann et al. (Reference Grossmann, Lohse and Sun2016). Here, 
$\boldsymbol {u}=(u_r,u_\theta,u_z)$ is the non-dimensional velocity vector, 
$p$ is pressure and 
$T=(T^* - T_o^*)/(T_i^* - T_o^*)$ is the temperature difference. The associated boundary conditions include the no-slip and no-penetration boundary conditions for the velocity and the Dirichlet boundary conditions for the temperature, i.e.
\begin{gather} \boldsymbol{u}(r_i) = \boldsymbol{e}_{\theta} \ \ \text{when }\omega_i^*>\omega_o^*, \quad \boldsymbol{u}(r_i) ={-}\boldsymbol{e}_{\theta} \ \ \text{when }\omega_i^*<\omega_o^*, \quad \boldsymbol{u}(r_o) = \boldsymbol{0}, \end{gather}
\begin{gather} T(r_i) = 1, \quad T(r_o) = 0 \quad \text{with} \ r_i=\frac{\eta}{1-\eta}, r_o=r_i+1. \end{gather}
Here, 
$\boldsymbol {e}_{r},\boldsymbol {e}_{\theta },\boldsymbol {e}_{z}$ are the unit vectors in the cylindrical coordinate.
Regarding the gravitational acceleration profile in (2.1a), it should be noted that from Newton's Law of universal gravitation, for uniform mass distribution along the 
$z$ axis in cylindrical coordinates, the gravitational acceleration magnitude is 
$g^*(r^*)=2G^* \lambda ^* /r^*$, where 
$G^*$ is the universal gravitational constant and 
$\lambda ^*$ is the mass per unit length of the 
$z$ axis. After non-dimensionalisation with 
$g^*(r_o^*)$, it gives the following expression in the non-dimensional form
\begin{equation} g(r) = \frac{r_o}{r}. \end{equation}
To compare, the gravitational potential due to a point mass at the centre of the central object in an accretion disk is 
$\phi ^*(r^*)=-G^* M^*/r^{*}$, with the corresponding gravitational acceleration magnitude being 
$g^*(r^*)=G^* M^*/r^{*2}$ when described in spherical coordinates. Non-dimensionalisation of this 
$g^*(r^*)$ with 
$g^*(r_o^*)$ results in a profile of the form 
$g(r)=r_o^2/r^2$. The third profile that can be considered is an artificial constant profile 
$g(r)=1$, which has been shown to be experimentally realisable on thermo-electric convection in a cylindrical annulus during a sounding rocket flight (Antoine et al. Reference Antoine, Martin, Vasyl and Christoph2023). In § 4, the profile given in (2.3) is studied. As suggested by one of the reviewers, a brief comparison of the effects of the other two profiles with (2.3) is conducted; see Appendix B.
In this equation system, there are five governing parameters defined as
\begin{equation} \left. \begin{gathered}
\eta = \frac{r_i^*}{r_o^*},\quad Ta = \frac{(1+\eta)^4}{64
\eta^2}\frac{{d}^{*2}(r_i^* + r_o^*)^2 (\omega_i^* -
\omega_o^*)^2}{\nu^{*2}},\quad Pr =
\frac{\nu^*}{\kappa^*},\\ Ra = \frac{\alpha^*
g_o^*(T_i^* - T_o^*)d^{*3}}{\nu^* \kappa^*},\\
Ro=\frac{|\omega_i^* - \omega_o^*|r_i^*}{2\omega_o^* d^*} =
\frac{\eta |1-a|}{2(1-\eta)a} \quad \text{with the rotation
ratio} \ a=\frac{\omega_o^*}{\omega_i^*}=\eta^{q},
\end{gathered} \right\}
\end{equation}
where the rotation exponent parameter 
$q$ is real valued for co-rotating cylinders (
$a>0$), the case considered in the present study. Specifically, when 
$\omega _i^*>\omega _o^*>0$, we have 
$0< a<1$ and, thus, 
$q>0$. The physical significance of these governing parameters is as follows:
(i) The radius ratio
$\eta$ determines the cylinders’ gap width as well as inner and outer radii $r_i$ and $r_o$ since we require $r_o-r_i=1$ according to the non-dimensionalisation step. Small values of $\eta$ correspond to large accretion disks with small central stars.(ii) The Taylor number
$Ta$ defines the relative importance of fluid inertia relative to viscosity. For highly rarefied gases in fast rotating motions, $Ta$ can be extremely high, which will be considered in our work.(iii) The Prandtl number
$Pr$ for most gases is of order one, meaning that momentum and heat diffuse at comparable rates. For rarefied gases, $Pr$ would be smaller.(iv) The Rayleigh number
$Ra$ characterises the importance of natural thermal convection relative to thermal diffusion in the quasi-Keplerian flow. Weak thermal stratification dictates small values of $Ra$. In addition to $Ra$, we will also adopt the Richardson number $Ri\equiv Ra/(PrTa)$ that characterises the relative importance between the convection due to thermal buoyancy and the flow shearing.(v) Finally, the Rossby number
$Ro$ quantifies the relative rotation rate of the two cylinders. It should be emphasised that, for a given $\eta$, $Ro$ depends only on $q$. The Keplerian regime corresponds to a specific rotation scenario as discussed below.
The rotation exponent 
$q$ is important in delimiting the various flow regimes illustrated in figure 1(b): 
$q=0$ indicates solid-body rotation; 
$q>2$ corresponds to the conventional Rayleigh unstable regime of isothermal flows; 
$0< q<2$, resulting in 
$\omega _i^* r_i^{*2} < \omega _o^* r_o^{*2}$, represents the Rayleigh stable regime where the angular momentum increases monotonically with 
$r$. The special case at 
$q=3/2$, where the rotation rates of the two cylinders obey Kepler's third law, is referred to as the Keplerian regime. Conventionally, the entire range of 
$0< q<2$ is referred to as the quasi-Keplerian regime. The present study focuses on the case at 
$q=3/2$.
2.2. Linear stability analysis: inhomogeneous in the radial direction
The equation system (2.1)–(2.4) described above admits the following one-dimensional laminar steady solution, which can be theoretically derived by noting that the only non-zero velocity component is the azimuthal velocity and the solution is a function of 
$r$. For the investigated case of 
$\omega _i^*>\omega _o^*>0$, the laminar solution reads
\begin{align} U_{b,r}=0,\, U_{b,\theta}(r)=\frac{1}{\eta(1+\eta)}\left(\frac{r_i^2}{r} - \eta^2 r\right),\, U_{b,z}=0,\, T_b(r) = \frac{1}{\ln \eta}\ln r - \frac{\ln r_o}{\ln \eta}.\end{align}
Here the pressure profile of the laminar solution 
$P_b(r)$ is not provided since it is not needed in the linear stability analysis.
To probe the stability or instability of the above laminar base flow subject to disturbances, we decompose the total flow variables in (2.1) into their base states plus perturbations:
\begin{equation} \boldsymbol{u}=\boldsymbol{U}_b + \boldsymbol{u}', \quad p = P_b + p', \quad T = T_b + T'. \end{equation}Then linearisation around the base flow leads to the linear perturbation equation
\begin{gather} \frac{\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{u}' \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{U}_b + \boldsymbol{U}_b \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}'={-} \boldsymbol{\nabla}p' + \frac{f(\eta)}{{Ta}^{1/2}} \nabla^2 \boldsymbol{u}' + \frac{Ra f^2(\eta)}{Pr Ta}T' g(r) \boldsymbol{e}_r - \frac{1}{Ro} \boldsymbol{e}_z \times \boldsymbol{u}', \end{gather}
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}' = 0, \end{gather}
\begin{gather} \frac{\partial T'}{\partial t} + \boldsymbol{U}_b \boldsymbol{\cdot}\boldsymbol{\nabla} T' + \boldsymbol{u}' \boldsymbol{\cdot}\boldsymbol{\nabla} T_b = \frac{f(\eta)}{Pr {Ta}^{1/2}} \nabla^2 T', \end{gather}along with homogeneous Dirichlet boundary conditions for both the perturbed velocity and temperature. The solution to (2.7) can be assumed to take the normal-mode form, i.e.
\begin{equation} (\boldsymbol{u}',p',T') = \left[\tilde{\boldsymbol{u}}'(r),\tilde{p}'(r),\tilde{T}'(r)\right]\exp(-{\rm i}\omega t + {\rm i}m\theta+{\rm i}k z) + \text{c.c.}, \end{equation}
where integer valued 
$m$ and real valued 
$k$ are the azimuthal and axial wavenumbers, respectively. Variables with the symbol tilde 
$\tilde {\ }$ are shape functions in the radial direction and c.c. represents the complex conjugate of the preceding term. Additionally, 
$\omega =\omega _r + \textrm {i}\omega _i$ includes 
$\omega _r$ as the frequency and 
$\omega _i$ as the linear growth rate. At a given parameter setting, the flow is said to be linearly unstable if 
$\omega _i>0$, linearly stable if 
$\omega _i<0$ and neutral if 
$\omega _i=0$.
Inserting (2.8) into the perturbation equations (2.7) results in the linearised equations in the spectral space, which are given as
\begin{align}
-{\rm i}\omega
\tilde{u}'_r &= \left[-\frac{{\rm i}mU_{b,\theta}}{r} +
\frac{f(\eta)}{{Ta}^{1/2}} \left(\tilde{\nabla}^2 -
\frac{1}{r^2}\right)\right]\tilde{u}'_r \nonumber\\ &\quad
+ \left[\frac{2U_{b,\theta}}{r} +
\frac{f(\eta)}{{Ta}^{1/2}} \left(-\frac{2{\rm
i}m}{r^2}\right) + \frac{1}{Ro} \right]
\tilde{u}'_{\theta}- \frac{{\rm d} \tilde{p}'}{{\rm d} r} +
\frac{Ra f^2(\eta)}{Pr Ta}\tilde{T}'
g(r),
\end{align}
\begin{align}
-{\rm i}\omega
\tilde{u}'_{\theta} &= \left[-\frac{{\rm d}
U_{b,\theta}}{{\rm d} r}-\frac{U_{b,\theta}}{r} +
\frac{f(\eta)}{{Ta}^{1/2}} \left(\frac{2{\rm
i}m}{r^2}\right)- \frac{1}{Ro} \right]\tilde{u}'_r
\nonumber\\ &\quad + \left[-\frac{{\rm i}mU_{b,\theta}}{r}
+ \frac{f(\eta)}{{Ta}^{1/2}} \left(\tilde{\nabla}^2 -
\frac{1}{r^2}\right) \right] \tilde{u}'_{\theta}-
\frac{{\rm i}m
\tilde{p}'}{r},
\end{align}
\begin{align} -{\rm i}\omega
\tilde{u}'_{z} &= \left(-\frac{{\rm i}mU_{b,\theta}}{r} +
\frac{f(\eta)}{{Ta}^{1/2}} \tilde{\nabla}^2 \right)
\tilde{u}'_{z} - {\rm i}k \tilde{p}',
\end{align}
\begin{align} 0 &= \left(\frac{{\rm d}}{{\rm d} r} + \frac{1}{r} \right) \tilde{u}'_{r} +\frac{{\rm i}m }{ r}\tilde{u}'_{\theta} + {\rm i}k \tilde{u}'_{z}, \end{align}
\begin{align} -{\rm i}\omega \tilde{T}' &={-}\frac{{\rm d} T_b}{{\rm d} r} \tilde{u}'_{r} + \left(-\frac{{\rm i}mU_{b,\theta}}{r} + \frac{f(\eta)}{Pr {Ta}^{1/2}} \tilde{\nabla}^2 \right) \tilde{T}', \end{align}
where the Laplacian 
$\tilde {\nabla }^2 = \textrm {d}^2 /\textrm {d}r^2 + 1/r (\textrm {d}/\textrm {d}r) - m^2/r^2 - k^2$. The above equation system can be written in a generalised eigenvalue problem in the form of
\begin{equation} -{\rm i}\omega
\tilde{\boldsymbol{M}}\,\tilde{\boldsymbol{\gamma}}=\tilde{\boldsymbol{L}}\,
\tilde{\boldsymbol{\gamma}},
\end{equation}
or, more specifically,
\begin{equation} -{\rm i}\omega
\begin{pmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 \\ 0
& 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1
\end{pmatrix} \begin{pmatrix} \tilde{u}'_{r} \\
\tilde{u}'_{\theta} \\ \tilde{u}'_{z} \\ \tilde{p}' \\
\tilde{T}' \end{pmatrix} = \begin{pmatrix} L_{r,r} &
L_{r,\theta} & 0 & L_{r,p} & L_{r,T}\\ L_{\theta,r} &
L_{\theta,\theta} & 0 & L_{\theta,p} & 0 \\ 0 & 0 & L_{z,z}
& L_{z,p} & 0 \\ L_{p,r} & L_{p,\theta} & L_{p,z} & 0 & 0
\\ L_{T,r} & 0 & 0 & 0 & L_{T,T} \end{pmatrix}
\begin{pmatrix} \tilde{u}'_{r} \\ \tilde{u}'_{\theta} \\
\tilde{u}'_{z} \\ \tilde{p}' \\ \tilde{T}' \end{pmatrix}.
\end{equation}
Here the linear operators 
$L_{*,*}$ can be derived by matching the equations term by term. The perturbation variable vector 
$\tilde {\boldsymbol {\gamma }}=(\tilde {u}'_{r},\tilde {u}'_{\theta },\tilde {u}'_{z},\tilde {p}',\tilde {T}')^\textrm {T}$ is identified as the eigenvector. Since the eigenvector can be arbitrarily scaled in the linear analysis, to eliminate ambiguity, it is normalised to have a unit amplitude and zero phase angle in the azimuthal velocity component at the middle point of the gap 
$\tilde {u}'_{\theta }(r={(r_i+r_o)}/{2})$, unless specified otherwise.
2.3. Linear stability analysis: local (or homogeneous) in the radial direction
Local stability analysis is a widely used method for investigating flow instability in modelled accretion-disk problems (Klahr & Hubbard Reference Klahr and Hubbard2014; Lyra Reference Lyra2014; Latter Reference Latter2016; Held & Latter Reference Held and Latter2018). Among these studies, Klahr & Hubbard (Reference Klahr and Hubbard2014) and Lyra (Reference Lyra2014) employed a short wave approximation in their local analysis to explore the convective overstability in radially stratified disks under thermal relaxation. In contrast, Latter (Reference Latter2016) and Held & Latter (Reference Held and Latter2018) utilised the shearing box approximation to analyse vertically stratified disk models. Following the suggestion from one of the reviewers, we conduct a local stability analysis of the problem formulated in § 2.1 and compare the results with those from our global analysis (global in the radial direction, yet local in the axial direction).
Given that the flow under investigation is radially stratified, the shearing box approximation in a local rotating Cartesian coordinate system is unsuitable for formulating an eigenvalue problem. The main challenge is implementing the shearing periodic boundary condition in the radial direction along which thermal stratification and gravitational acceleration exist. While the implementation is feasible for numerical simulations, it is not appropriate for the present analysis. Therefore, we adopt the short wave approximation, following the approach of Klahr & Hubbard (Reference Klahr and Hubbard2014) and Lyra (Reference Lyra2014). In this approximation, perturbations are expressed in the Fourier mode form as follows:
\begin{equation} (\boldsymbol{u}',p',T') = [\hat{\boldsymbol{u}}',\hat{p}',\hat{T}']\exp(-{\rm i}\omega t + {\rm i} k_r r + {\rm i}m\theta+{\rm i}k_z z) + \text{c.c.} \end{equation}
Note the distinction from (2.8): wavelike solutions are also assumed in the radial direction here, rendering this analysis fully local in all three directions. In this context, 
$k_r$ and 
$k_z$ represent the radial wavenumber and axial wavenumber, respectively. Noting that the linear instability to be reported in § 4 below is predominantly axisymmetric (
$m=0$), in this local analysis we focus on the case where 
$m=0$, meeting the assumptions 
$m\ll k_r r$ and 
$m\ll k_z z$ for the short wave approximation.
Inserting (2.12) into the perturbation equation (2.7) leads to the following equation system for the local quasi-Keplerian flow at a radial location 
$r_0$:
\begin{gather} -{\rm i}\omega \hat{u}'_r = \frac{f(\eta)}{{Ta}^{1/2}} \hat{\nabla}^2\hat{u}'_r + \left(\frac{2U_{b,\theta}(r_0)}{r_0} + \frac{1}{Ro} \right) \hat{u}'_{\theta}- {\rm i} k_r \hat{p}' + \frac{Ra f^2(\eta)}{Pr Ta}\hat{T}' g(r_0), \end{gather}
\begin{gather} -{\rm i}\omega \hat{u}'_{\theta} = \left[-\left(\frac{{\rm d} U_{b,\theta}}{{\rm d} r}\right)_{r_0}-\frac{U_{b,\theta}(r_0)}{r_0} - \frac{1}{Ro} \right]\hat{u}'_r + \frac{f(\eta)}{{Ta}^{1/2}} \hat{\nabla}^2 \hat{u}'_{\theta}, \end{gather}
\begin{gather} -{\rm i}\omega \hat{u}'_{z} =\frac{f(\eta)}{{Ta}^{1/2}} \hat{\nabla}^2 \hat{u}'_{z} - {\rm i}k_z \hat{p}', \end{gather}
\begin{gather} 0 = {\rm i} k_r \hat{u}'_{r} + {\rm i}k_z \hat{u}'_{z}, \end{gather}
\begin{gather} -{\rm i}\omega \hat{T}' ={-}\left(\frac{{\rm d} T_b}{{\rm d} r}\right)_{r_0} \hat{u}'_{r} + \frac{f(\eta)}{Pr {Ta}^{1/2}} \hat{\nabla}^2 \hat{T}'. \end{gather}
Here the Laplacian 
$\hat {\nabla }^2 = - k_r^2 - k_z^2$. In this analysis, some less significant terms have been omitted, following Lyra (Reference Lyra2014). For instance, for an infinitesimal disturbance to the local flow at 
$r_0$, it is assumed that 
$u'_r/r\ll \partial _r u'_r$. Consequently, the term 
$u'_r/r$ in the continuity equation is neglected. Similarly, we neglect the term 
$1/r\partial _r u'_r$ compared with 
$\partial _{rr} u'_r$ in the Laplacian. The resulting equations (2.13) form a 
$5\times 5$ eigenvalue problem, which is computationally much more efficient to solve than (2.11). Note that no discretisation is needed for solving (2.13) in the local analysis. In contrast, (2.11), after discretisation based on the spectral collocation method, results in a 
$5N\times 5N$ eigenvalue problem for the global analysis presented in § 4, where the node number 
$N$ ranges from 51 for normal parameters to 601 for extreme parameters in our calculations.
3. Numerical method and validation
To solve the generalised linear eigenvalue problem expressed in (2.10), we developed an in-house MATLAB code using the spectral collocation method (Trefethen Reference Trefethen2000). This approach employs Chebyshev–Lobatto nodes, clustered near both the inner and outer cylinders while being sparse at the gap centre. For Chebyshev polynomials of 
$(N-2)$th order, there are in total 
$N$ Chebyshev–Lobatto nodes. In our calculations, a resolution of 
$N=51$ typically proves sufficient, although for extreme parameter regimes, we extend up to 
$N=601$. Central to our methodology is the construction and utilisation of Chebyshev differentiation matrices, facilitating the formation of 
$\tilde {\boldsymbol {M}}$ and 
$\tilde {\boldsymbol {L}}$, as described in Trefethen (Reference Trefethen2000). Implementing homogeneous boundary conditions for 
$\tilde {\boldsymbol {u}}'$ and 
$\tilde {T}'$ at the inner cylinder 
$r=r_i$ and outer cylinder 
$r=r_o$ is straightforward, accomplished by removing the corresponding rows and columns from 
$\tilde {\boldsymbol {M}}$ and 
$\tilde {\boldsymbol {L}}$. Subsequently, we employ MATLAB's built-in functions ‘eig()’ and ‘eigs()’ to solve the generalised linear eigenvalue problem (2.10), yielding the eigenvalue 
$\omega$ and the eigenvector 
$\tilde {\boldsymbol {\gamma }}$. A validation against the data in Deguchi (Reference Deguchi2017) and a convergence test of our code can be found in § 1 of the supplementary material; favourable agreement and good convergence have been achieved. For the local analysis, the small eigenvalue problem (2.13) is also solved using the two built-in functions.
4. Results and discussion
To clearly present the numerical results below, a summary of the following subsections is provided at the outset.
(i) Section 4.1 characterises the thermo-hydrodynamic linear instability at a moderate Taylor number of
$Ta=10^6$ in the TC flow. While this value may not reach the extreme parameters relevant to accretion-disk dynamics, its selection facilitates a clear elucidation of the primary destabilising mechanism. We demonstrate the robustness of the flow pattern in the corresponding unstable mode, noting that even at significantly higher $Ta$, the eigenvector pattern bears resemblance to that observed at $Ta=10^6$.(ii) Section 4.2 delves into an exploration of the effects of the principal governing parameters, including
$Ta$, $Ri$, $Pr$ and $\eta$. In a linear stability analysis, a key point of interest is to determine the lowest Taylor number, referred to as the critical Taylor number $Ta_c$, beyond which the linear instability sets in. Notably, we illustrate that the instability exhibits resilience even at very low yet finite Richardson numbers $Ri\sim 0.01$, indicative of weak thermal effects. The critical $Ta_c$ appears to scale with the deviation of $Ri$ from the limiting $Ri$ corresponding to infinite $Ta_c$. Additionally, another scaling law emerges at asymptotically small Prandtl numbers $Pr\sim 10^{-4}$. By extrapolating these findings to high-$Ta$ regimes, we infer the pertinence of our observations to the dynamics of accretion disks. Particularly, § 4.2.5 is the only part dedicated to the local stability analysis for comparison with the global one. Qualitatively consistent results are obtained; differences are underscored.(iii) The final section, § 4.3, is dedicated to an examination of an extreme case scenario characterised by
$\eta =0.05$, $Pr=10^{-3}$, $Ri=10^{-3}$ and $10^{9}\lessapprox Ta \lessapprox 10^{16}$. Within the linear framework, this scenario mimics the hydrodynamics in accretion disks based on the TC flow geometry.
4.1. Characteristics of the thermo-hydrodynamic linear instability
Within the Keplerian regime, characterised by 
$q=1.5$, we observe a linear instability within the thermal TC flow. The dynamics of this instability exhibit a remarkably rich complexity. To illustrate this, we maintain fixed values for the governing parameters: 
$\eta =0.3$, 
$Pr=0.7$, 
$Ri=0.1$ and 
$Ta={10}^6$. The selection of these parameters is guided by several considerations. Firstly, the chosen radius ratio, though relatively small, remains achievable in both experimental set-ups and numerical simulations. This paves the way for investigating flow instability in future works, particularly those exploring its nonlinear evolution and potential turbulence simulations. Secondly, the specific 
$Pr$ chosen reflects a typical value for most gases under atmospheric conditions. However, it is important to note that while hydrogen and helium primarily constitute the gas components in accretion disks (Pringle Reference Pringle1981; Mineshige Reference Mineshige1993; Ikoma & Hori Reference Ikoma and Hori2012; Klahr & Hubbard Reference Klahr and Hubbard2014), the corresponding 
$Pr$ therein can be significantly smaller (Latter Reference Latter2016). We investigate such extreme cases in §§ 4.2 and 4.3. Lastly, the selected 
$Ri$ aligns with commonly used values in the literature (Lesur & Ogilvie Reference Lesur and Ogilvie2010; Held & Latter Reference Held and Latter2018). The variation of these parameters will be explored in the next subsection.
Figure 2(a) displays the converged results with 
$N=(51,151)$ for the eigenspectra of the axisymmetric (
$m=0$) modes, where the two modes above 
$\omega _i=0$ (the dashed line) are linearly unstable. These are stationary eigenmodes with zero real parts, meaning that they grow over time but do not travel in space. Non-axisymmetric modes (
$m\ne 0$) can also become linearly unstable, further classified into two types. One type is axial independent (
$k=0$) with a single unstable mode, as shown in panel 2(b), and the other is axial dependent (
$k\neq 0$) with two unstable modes, also called helical modes, as displayed in panel 2(c). The frequencies of these unstable non-axisymmetric modes are non-zero, meaning they oscillate in time and travel in the azimuthal direction and also in the axial direction if 
$k\neq 0$.
Eigenspectra of the flow in the Keplerian regime at 
$\eta =0.3$, 
$Pr=0.7$, 
$Ri=0.1$ and 
$Ta={10}^6$ for (a) 
$(k,m)=(5,0)$, (b) 
$(k,m)=(0,1)$, (c) 
$(k,m)=(5,1)$. The modes above the dashed line are linearly unstable. (d) The corresponding linear growth rate 
$\omega _i$ in the 
$(k,m)$ plane where the black curve traces 
$\omega _i=0$. Note that 
$m$ is a non-negative integer; at each combination of 
$(k,m)$, 
$\omega _i$ of the most unstable/least stable mode in an eigenspectrum is recorded for the plot. The red and blue stars in panel (d) respectively mark the locations of the most unstable mode for all 
$(k,m)$ and the most unstable non-axisymmetric mode for all 
$k$ and 
$m\neq 0$. The corresponding mode patterns are visualised in figure 3.
In the above, we have reported two categories of unstable modes: stationary axisymmetric modes and oscillatory non-axisymmetric modes. However, axisymmetric modes can also become oscillatory at specific parameters, as will be shown in figure 15(b). In a different TC flow system with a geophysical application background, Jenny & Nsom (Reference Jenny and Nsom2007) reported the coexistence of these three kinds of mode. The difference between their study and ours is that they fixed the outer cylinder and examined the radial buoyancy effects due to temperature and/or salinity stratification in the Earth's equatorial region. The difference between the two settings, along with the similarity in the results, demonstrates the robustness of the instability in thermal TC flows with different configurations. Distinguishing these different unstable modes is not only of interest within the linear framework but also crucial from a (weakly) nonlinear perspective when examining flow bifurcations. For example, Kang et al. (Reference Kang, Meyer, Mutabazi and Yoshikawa2017) studied the effect of centrifugal buoyancy on circular Couette flow and found that the flow bifurcation to stationary axisymmetric modes was always supercritical, while that to oscillatory axisymmetric modes was subcritical. Therefore, it would also be interesting to examine the weak nonlinearity of the present thermally stratified quasi-Keplerian flow in future studies.
To gain a comprehensive view of the most unstable/least stable modes in the axial- and azimuthal-wavenumber space, we plot the linear growth rate 
$\omega _i$ in the 
$k$–
$m$ plane in panel 2(d). The largest growth rate 
$\omega _i \approx 0.5332736016$ is observed at 
$m=0$ and 
$k\approx 5.685$ (red star in panel d). This indicates that the most unstable mode in the quasi-Keplerian flow at 
$(q,\eta,Pr,Ri,Ta)=(1.5,0.3,0.7,0.1,10^6)$ is axisymmetric. To examine the pattern of this mode, we visualise its temperature contours superposed with the velocity field in the 
$r$–
$z$ plane in figure 3(a). From this figure, we observe that the axial scale of the mode is comparable to the gap width, with the velocity field representing a pair of counter-rotating vortices aligned in the periodic 
$z$ direction. This convection feature is reminiscent of the canonical Rayleigh–Bénard convection (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000) where gravitational buoyancy drives fluid from the hot wall to the cold wall, resulting in two counter-rotating vortices per wavelength. This suggests that the thermal effects dominate in rendering the instability in quasi-Keplerian flows, as the Keplerian shearing is known to be strongly stabilising.
Contours of disturbance temperature in the quasi-Keplerian flow at 
$\eta =0.3$, 
$Pr=0.7$, 
$Ri=0.1$ and 
$Ta={10}^6$. Panel (a) is for the most unstable mode attained at 
$k\approx 5.685$ and 
$m=0$, corresponding to the red star marked in figure 2(d). Panel (b) is for the leading non-axisymmetric mode at 
$k\approx 5.049$ and 
$m=1$, corresponding to the blue star marked in figure 2(d). In panel (a) one wavelength 
$2{\rm \pi} /k\approx 1.11$, comparable to the cylinder gap 
$d=1$, is shown in the 
$z$ direction. The green circle in panel (b) marks the location of the critical layer.
In addition to the above most unstable axisymmetric mode, the helical mode with the largest linear growth rate 
$\omega _i=0.3658437831$, attained at 
$m=1$ and 
$k\approx 5.049$ (blue star in figure 2d), is visualised in figure 3(b) in the 
$r$–
$\theta$ plane at 
$z=0$. This mode manifests a spiral structure rotating counterclockwise along with the laminar base flow. The two ends of the two spiral arms are located where hot fluids decelerate upon encountering the outer wall and then turn back towards the inner wall. Scrutinising the velocity vector field reveals that the mode structure actually consists of two counter-rotating vortices severely distorted by the Keplerian shearing in the bulk flow. The frequency of the above most unstable helical mode is 
$\omega _r=0.5662622568$, corresponding to a phase speed 
$c_r = \omega _r / m\approx 0.566$. This indicates that the disturbance travels in the same direction as the bulk flow, with a speed about 
$24\,\%$ of the maximum angular frequency 
$2.333$ at the inner cylinder wall. Since the angular frequency profile of the laminar base flow can be calculated as 
$\varOmega =U_{b,\theta }/r$, with 
$\eta =0.3$, 
$r_o=1/(1-\eta )\approx 1.429$ and 
$r_i=r_o - 1\approx 0.429$, the above phase speed 
$c_r$ is equal to the local value of 
$\varOmega$ at a location 
$r\approx 0.769$. This position, where the azimuthal phase speed of a mode is locally equal to the angular velocity of the base flow, is identified as the critical layer in curvilinear shear flows, a concept that has been employed by Leclercq et al. (Reference Leclercq, Partridge, Augier, Caulfield, Dalziel and Linden2016) in their study on stratified TC flow. In the present study, this position, marked by a green circle in figure 3(b), is located 
$0.66$ length scale from the outer wall and 
$0.34$ from the inner wall. From the perspective of the critical-layer concept developed in the context of parallel shear flows, this location is identified as the most receptive position for energy amplification (Maslowe Reference Maslowe1986; McKeon Reference McKeon2017). It can be seen that the critical layer in figure 3(b) approximately coincides with the radial location where the disturbance magnitude is the largest, suggesting its relevance in the present curvilinear context. Nevertheless, its importance remains to be examined in the future.
Although the above results pertain to a moderate Taylor number (
$Ta=10^6$), which is not high enough to be directly linked to accretion disks, the mode patterns and underlying instability mechanisms can be applied to high-
$Ta$ regimes. Numerical evidence at 
$Ta=(10^8, 10^{10})$ is provided in § 3 of the supplementary material. We have checked that the patterns persist even at 
$Ta=10^{16}$, but the axial length scale is too small and, thus, not shown. As 
$Ta$ increases, both the most unstable axisymmetric mode and the helical mode become increasingly localised around the inner rotating cylinder, suggesting that instability space is diminished by the enhanced Keplerian shearing from fast rotation. This localisation can be relevant for the accretion-disk dynamics for two reasons. First, accretion disks are typically flat and thin structures (Abramowicz & Straub Reference Abramowicz and Straub2014), making it difficult for disturbances with large axial wavelengths to survive, leaving short-wavelength disturbances as primary candidates for triggering flow transitions. For example, Klahr & Hubbard (Reference Klahr and Hubbard2014) adopted a quasi-hydrostatic approximation in their study of convective overstability in radially stratified accretion disks, requiring perturbations to be vertically thin. Additionally, Latter (Reference Latter2016) estimated that the dominant vertical wavelength of flow structures in a protoplanetary disk is about 
$1\,\%$ of the disk thickness at 1 AU. Second, Latter (Reference Latter2016) found that the convective overstability discovered by Klahr & Hubbard (Reference Klahr and Hubbard2014) was not prevalent in the outer disk regions but seemed only possible in the inner regions, which may be corresponding to the narrow region around the inner rotating cylinder at high 
$Ta$ in our study. Comparatively, it is interesting to note that, for vertical thermal stratification under vertical stellar gravity, Held & Latter (Reference Held and Latter2018) found that the fastest growing mode was at the shortest radial length scales, manifesting as radially thin elongated structures.
From these perspectives, our findings are generally consistent with previous research on thermally driven linear instability in accretion disks. However, previous studies often used 2-D local shearing rectangles (Klahr & Hubbard Reference Klahr and Hubbard2014; Latter Reference Latter2016) or 3-D local shearing boxes (Lyra Reference Lyra2014; Held & Latter Reference Held and Latter2018), assuming homogeneity in all directions of a cylindrical coordinate system. These studies bear local significance, whereas the linear instability in our work is of a 3-D global nature in the TC configuration. Clearly, the mode structures observed in our study cannot be fully contained within their local shearing boxes. As suggested by one of the reviewers, to evaluate the similarities and differences between the present global analysis and the local analysis of the same quasi-Keplerian flow, we have also conducted a local analysis based on the short wave approximations employed in Klahr & Hubbard (Reference Klahr and Hubbard2014) and Lyra (Reference Lyra2014). The results, as will be presented in § 4.2.5, are generally consistent with those from the global analysis. However, quantitatively, some scaling laws (to be discussed in the next subsection) differ. In addition, none of the unstable mode structures observed in the global analysis can be captured in the local analysis. Besides, Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) did not arrive at the same conclusion presented above as they focused on small 
$Ta$ in their stability analysis.
To sum up, the analysis in this subsection indicates that the destabilising thermally driven buoyancy can overcome the stabilising Keplerian shear to render the linear instability in the quasi-Keplerian flow. As will be further demonstrated in § 4.3, this linear instability persists even when the thermal buoyancy effect is very weak but finite, despite the strong Keplerian shear. It can arise at exceptionally large axial wavenumbers (small axial wavelengths), which are especially relevant to the flat geometry of accretion disks. This demonstrates the potential prevalence of thermo-hydrodynamic linear instability in accretion disks, particularly in the inner radii regions, which is one of the most significant findings of this study.
4.2. Effects of governing parameters on the flow instability
4.2.1. Effects of Taylor number $Ta$
In the previous subsection the Taylor number is primarily fixed at 
$Ta={10}^6$. Next, we explore the effect of varying 
$Ta$ while keeping other control parameters fixed in the Keplerian regime at 
$(q,\eta,Pr,Ri)=(1.5,0.3,0.7,0.1)$. First, figure 4(a) shows that the linear growth rate 
$\omega _{i,{max}}$ of the most unstable disturbance, optimised over all axial wavenumbers 
$k$, increases monotonically with 
$Ta$ for any azimuthal wavenumber 
$m$. This indicates that faster rotation (corresponding to larger 
$Ta$) enhances the linear instability. Second, the associated axial wavenumber 
$k$, plotted in figure 4(b), also increases with 
$Ta$. This trend features an increasingly thinner mode structure (see figure 1 of the supplementary material), which aligns more closely with the thin nature of the accretion disks. Third, the most unstable mode is almost always axisymmetric (except 
$Ta=10^5$), as shown in panel 4(a).
(a) Variations of the largest linear growth rate 
$\omega _{i,{max}}$ with the azimuthal wavenumber 
$m$ for various Taylor numbers 
$Ta$ for the quasi-Keplerian flow at 
$\eta =0.3$, 
$Pr=0.7$ and 
$Ri=0.1$. (b) Variations of the corresponding axial wavenumber 
$k$ at which 
$\omega _{i,{max}}$ is attained. The discontinuity in the variation of 
$k(\omega _{i,{max}})$ is due to a mode shift from one branch to another branch in the corresponding eigenspectrum.
Figure 4 also shows that as 
$Ta$ increases, more non-axisymmetric modes become linearly unstable, although they do not dominate. Specifically, for high 
$m$, the most unstable mode is always axial independent (
$k=0$). Typical non-axisymmetric modes, which exhibit an increasingly localised structure around the outer rotating cylinder, are briefly discussed in § 4 of the supplementary material; see figure 3 therein for the visualisation. Here, we discuss a potential connection between the observed complex thermal instability diagram and MRI. In this study, the magnetic field is omitted to isolate and focus on thermal effects. However, in real accretion disks, both magnetic and thermal fields are typically present, even if one or both are relatively weak (see § 4.2.3 of Lesur et al. (Reference Lesur2022) for a discussion on the coupling of magnetohydrodynamics and thermodynamics in protoplanetary disks). The coexistence of these fields suggests that the thermally driven linear instability observed in this study may interact – either linearly or nonlinearly – with various manifestations of the MRI, e.g. those in Hollerbach & Rüdiger (Reference Hollerbach and Rüdiger2005). Recently, a synergy between thermal effects and magnetic effects on fluid convection has indeed been observed in a TC flow experiment by Seilmayer, Ogbonna & Stefani (Reference Seilmayer, Ogbonna and Stefani2020) and theoretically confirmed in a linear stability analysis by Mishra et al. (Reference Mishra, Mamatsashvili, Seilmayer and Stefani2024). Note that the gravity considered in their analysis is along the axial direction, unlike our study that models gravity in the radial direction, specifically for accretion disks.
To determine the critical 
$Ta_c$, neutral curves along which 
$\omega _i=0$ are depicted in figure 5 for the linearised flow at 
$(q,\eta,Pr,Ri)=(1.5,0.3,0.7,0.1)$. Panel (a) illustrates contours of the linear growth rate in a 
$Ta$–
$k$ plane, zoomed in around the linear critical condition (
$Ta_c\approx 3.7\times 10^4, k_c\approx 0.96, m_c=2$) marked by the red star and the black neutral curve. For other 
$m$ values, the Taylor number is always higher than the identified 
$Ta_c$, as evident in panel (b). Furthermore, it is notable that at 
$Ta$ significantly above its critical value, there exists a broad spectrum of axial wavelengths of disturbances capable of triggering both axisymmetric and non-axisymmetric linear instabilities. The wider this range and the higher the axial wavenumber, the more relevant the instability becomes to accretion disks.
(a) Contours of the linear growth rate 
$\omega _{i}$ in the 
$Ta$–
$k$ plane at 
$m=2$ for the quasi-Keplerian flow at 
$(\eta,Pr,Ri)=(0.3,0.7,0.1)$. The black curve is the neutral curve for 
$m=2$. The flow is linearly unstable on the right-hand side of the curves. (b) Neutral curves at various 
$m$. Note that, for 
$m=(1,3)$, each of the neutral curves consists of two segments, the area confined by which is linearly unstable; the small area at the bottom right corner is linearly stable for 
$m=(1,3)$ and unstable for 
$m=(4,5,6)$. Note the log scale for the horizontal axis in panel (b). The red star marks the linear critical condition (
$Ta_c\approx 3.7\times 10^4, k_c\approx 0.96, m_c=2$).
4.2.2. Effects of Richardson number $Ri$
The 
$Ri$ characterises the thermal effect relative to the shearing in our thermally driven quasi-Keplerian flow. We find that increasing 
$Ri$ promotes the linear instability. This can be seen from figure 6(a), where neutral curves in the 
$Ta$–
$k$ plane for various 
$Ri$ are plotted, with parameters 
$(\eta,Pr)=(0.3,0.7)$ fixed for axisymmetric flows. As our focus is on the linear critical 
$Ta_c$, only the nose segments of the neutral curves are shown. With increasing 
$Ri$ from 
$0.045$ to 
$1$, the neutral curve shifts from right to left, with the onset 
$Ta$ decreasing from 
$10^7$ to 
$10^3$, four orders of magnitude smaller. This suggests that when the destabilising thermal buoyancy effect is strong, the flow is more prone to instability and, thus, requires a weaker inertial effect, corresponding to a slower rotation, to overcome the stabilising Keplerian shear to trigger the instability.
Neutral curves for various Richardson numbers for the quasi-Keplerian flow at 
$\eta =0.3$ and 
$Pr=0.7$: (a) azimuthal wavenumber 
$m=0$, (b) 
$m=1$, (c) 
$m=2$, (d) 
$m=3$. In all the calculations, only the neutral curve segments close to the left noses are shown and the stars mark the linear critical conditions at each 
$Ri$ minimised over all 
$m$ and 
$k$ values.
In addition to the axisymmetric mode, similar destabilising effects of increasing 
$Ri$ are observed for non-axisymmetric modes; see the neutral curves for 
$m=(1,2,3)$ in panels 
$(b,c,d)$, respectively. Certain non-axisymmetric modes become linearly stable, resulting in the absence of neutral curves in the investigated range of 
$Ta$. Some neutral curves consist of disconnected parts or exhibit kinks, such as 
$(m,Ri)=(1,0.06)$ (blue dashed lines in panel b) and 
$(m,Ri)=(2,0.1)$ (black curve in panel 
$c$). Our result that 
$Ta_c$ increases with smaller 
$Ri$ is consistent with the numerical observations in Held & Latter (Reference Held and Latter2018), where a simulation study of thermally stratified quasi-Keplerian flows in a shearing box revealed that weakening the thermal buoyancy results in a linear instability at a larger Rayleigh number. Again, we would like to emphasise that our ‘global’ flow configuration is different than their shearing box.
The linear critical 
$Ta_c$ can be determined by minimising 
$Ta$ along a neutral curve for all combinations of 
$(m,k)$ at a given parameter setting. To achieve this, we extract the eight leftmost data points, marked by stars in figure 6, and display 
$(Ta_c,k_c,m_c)$ as functions of 
$Ri$ in figure 7. Here, we observe that 
$Ta_c$ increases monotonically with reducing 
$Ri$, with the variation being more pronounced at small 
$Ri$. Regarding the linear critical wavenumbers, for relatively high 
$Ri\geq 0.1$, 
$Ta_c$ are obtained for helical modes with 
$m_c=2$ and 
$k_c\approx 2$. In contrast, for relatively small 
$Ri=(0.045, 0.05)$, the linear instability only exists for 
$m=0$ within the investigated range of Taylor numbers, and 
$k_c$ dramatically increases upon decreasing 
$Ri$. These results suggest that when the thermally driven buoyancy effect is strong, as manifested at relatively large 
$Ri$, both the axisymmetric and non-axisymmetric disturbances can trigger instabilities, with the non-axisymmetric mode being more dominant at the instability onset. Conversely, for a weak buoyancy effect exhibited at relatively small 
$Ri$, the instability of non-axisymmetric modes seems to be severely suppressed, while the axisymmetric mode instability can still be operative.
Variations of the linear critical parameters with Richardson number 
$Ri$ for the quasi-Keplerian flow at 
$\eta =0.3$ and 
$Pr=0.7$: (a) linear critical Taylor number 
$Ta_c$ (see also the stars in figure 6), (b) linear critical axial wavenumber 
$k_c$, (c) linear critical azimuthal wavenumber 
$m_c$. These data are extracted from the critical points of the neutral curves in figure 6(a–c); see the eight stars therein.
The rapid increase of 
$(Ta_c,k_c)$ at small 
$Ri$ implies that by further reducing thermal buoyancy through decreasing 
$Ri$, both the Taylor number and axial wavenumber would increase asymptotically. To explore how these results vary with 
$Ri$, we look for scaling results in 
$(Ta_c,k_c)$ with 
$Ri$ decreasing from 
$0.05$, as shown in figure 8. Our computational results indicate that when 
$Ri\rightarrow Ri_{limit}=0.034945$, the 
$Ta_c$ approaches infinity under the current numerical resolution at 
$\eta =0.3, Pr=0.7$. Theoretically, an infinitely large 
$Ta_c$ means linear stability. As shown in figure 8(a), a clear scaling of 
$Ta_c$ is present with the deviation of 
$Ri$ from the 
$Ri_{limit}$, as well as the scaling in 
$k_c$ as a function of 
$Ri-Ri_{limit}$ in panel (b). To further corroborate the scaling result, the case of 
$Pr=0.01$ is additionally calculated, at which 
$Ri_{limit}\approx 0.00087$, showing the same scaling exponents. These scaling results may serve as a guide for future high-fidelity numerical investigations of the studied flow. It is noted that similar scalings in the critical parameters have also been observed in other complex fluids as a function of governing parameters, such as viscoelastic fluids; see Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018).
Scaling laws of (a) linear critical Taylor number 
$Ta_c \propto (Ri-Ri_{{limit}})^{-5}$ and (b) linear critical axial wavenumber 
$k_c \propto (Ri-Ri_{{limit}})^{-6/5}$ at the low-Richardson-number 
$Ri$ limit for axisymmetric (
$m=0$) quasi-Keplerian flow at 
$\eta =0.3$. The limit values of 
$Ri_{{limit}}$ are computationally estimated to be 
$Ri_{{limit}}\approx 0.034945$ for 
$Pr=0.7$ and 
$Ri_{{limit}}\approx 0.00087$ for 
$Pr=0.01$.
4.2.3. Effects of Prandtl number $Pr$
To explore the effects of the Prandtl number 
$Pr$, the linear critical Taylor number 
$Ta_c$ is calculated for various 
$Pr$ covering five orders of magnitude, with other parameters fixed at 
$(q,\eta,Ri)=(1.5,0.3,0.1)$. Several observations emerge from this analysis.
Firstly, the results in figure 9(a) reveal that, at 
$Ri=0.1$, high momentum diffusivity or low thermal diffusivity (meaning greater 
$Pr$) has destabilising effects on the linear instability in the quasi-Keplerian flow because 
$Ta_c$ monotonically decreases with increasing 
$Pr$ from 
$0.0001$ to 
$7$. However, at 
$Ri=0.01$, the effect of 
$Pr$ is initially destabilising as it increases from 0.0001 to 0.02, but it becomes strongly stabilising with further increases. The rapid increase of 
$Ta_c$ suggests a maximal 
$Pr$ beyond which the viscous dissipation is sufficiently strong to render the flow linearly stable, where the 
$Ta_c\rightarrow \infty$. In both cases of 
$Ri=0.1$ and 
$0.01$, the asymptotic slopes of the function curves approach 
$-\tfrac {6}{5}$ in the low-
$Pr$ limit, indicating a scaling law of 
$Ta_c \propto Pr^{-6/5}$. As pointed out by one of the reviewers, caution should be taken when one attempts to make any link between this trend to the typical scaling of the linear critical Reynolds number 
$Re_c\propto Pm^{-1}$, where 
$Pm$ is the magnetic Prandtl number, for the onset of MRI with a pure axial magnetic field in a TC geometry (Rüdiger, Schultz & Shalybkov Reference Rüdiger, Schultz and Shalybkov2003; Hollerbach & Rüdiger Reference Hollerbach and Rüdiger2005). Note that 
$Ta$ is a parameter corresponding to the square of the typical definition of 
$Re$.
Variations of the linear critical parameters as functions of Prandtl number 
$Pr$ for the quasi-Keplerian flow at 
$\eta =0.3$: (a) linear critical Taylor number 
$Ta_c$, (b) linear critical axial wavenumber 
$k_c$, (c) linear critical azimuthal wavenumber 
$m_c$. Here, 
$Ta_c \propto Pr^{-6/5}$ when 
$Pr$ is asymptotically small, but no scaling law in 
$k_c$ is observed. The inset of panel (b) shows the large variation of 
$k_c$ when 
$Ri=0.01$.
Secondly, figure 9(b) shows that the linear critical axial wavenumber 
$k_c$ mostly remains of 
$O(1)$ for all 
$Pr$ values at 
$Ri=0.1$. This result is consistent with that in figure 8(b), where 
$k_c$ can be extrapolated to 
$O(1)$ at 
$Ri=0.1$. Specifically, for intermediate Prandtl numbers 
$0.5\lessapprox Pr\lessapprox 1$, figure 9(b) indicates 
$k_c \sim 1$, meaning that the typical length scale of the these critical modes is comparable to the cylinder gap width. However, for relatively small or large values of 
$Pr\lessapprox 0.1$ or 
$Pr\gtrapprox 5$, 
$k_c$ is approximately equal to 
$3$. At 
$Ri=0.01$, the variation of 
$k_c$ as a function of 
$Pr$ presents a radically different behaviour. At the low-
$Pr$ limit, 
$k_c$ still remains at about 
$3$, but when the maximal 
$Pr$ is approached, 
$k_c$ increases dramatically; see the red crosses in the inset of panel (b). Accompanying the scaling law of 
$Ta_c$ at the low-
$Pr$ limit, the variation of 
$k_c$ is small, increasing by less than 
$10\,\%$ when 
$Pr$ is reduced by two orders of magnitude. No observable scaling law of 
$k_c$ is evident within the investigated range of Prandtl number.
Note that in both 
$Ri$ cases, 
$k_c\sim 3$ at the small-
$Pr$ limit. A similar value of 
$k_c\approx 3.13$ has been reported by Jenny & Nsom (Reference Jenny and Nsom2007) in a geophysical study based on the TC flow modelling. Likewise, Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) documented an asymptotic wavenumber of about 
$k_c\approx 3.1$ at large 
$Pr$ or 
$Ra$ for a TC flow with centrifugal buoyancy effects. The close value of 
$k_c$ in these different TC configurations suggests the generic features shared by the various flow instabilities in TC flows subjected to radial buoyancy due to stratification. It should be emphasised that the above (relatively) long-wavelength feature pertains to the linear critical conditions, whereas the short-wavelength modes, which are more astrophysically relevant, will dominate in the quasi-Keplerian flow for 
$Ta\gg Ta_c$, as illustrated in figure 4.
Lastly, it is found that axisymmetric modes dominate in the linear instability at high and low 
$Pr$. As illustrated in figure 9(c), for the 
$Pr$ values studied at 
$Ri=0.1$, 
$Ta_c$ is consistently attained at 
$m=0$, except for the cases 
$Pr=(0.5,0.7)$, for which 
$Ta_c$ is obtained for helical modes with 
$m=2$. Thus, non-axisymmetric modes are only significant at the instability onset in a narrow range of 
$Pr$. At 
$Ri=0.01$, the axisymmetric mode dominates the instability onset at all the 
$Pr$ investigated.
4.2.4. Effects of radius ratio $\eta$
In the preceding analysis the radius ratio of the co-rotating cylinders is held fixed at 
$\eta =0.3$. However, the relevant radius ratio in the TC flow to model the actual accretion disks may be much smaller than 
$0.3$ due to their inherently flat and thin structures around small central objects (Abramowicz & Straub Reference Abramowicz and Straub2014). In this section we investigate how varying 
$\eta$ affects the linear instability in the thermal TC flow.
Figure 10(a) reveals that reducing the radius ratio destabilises the quasi-Keplerian flow for the four selected combinations of 
$(Pr,Ri)$. Specifically, for the case at 
$(Pr,Ri)=(0.7,0.01)$, no linear instability can be found when 
$\eta \gtrapprox 0.2$. This demonstrates the significant influence of the curvature effect on the quasi-Keplerian flow stability/instability. In contrast, Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) observed that reducing the radius ratio stabilised quasi-Keplerian flows subjected to centrifugal buoyancy in the TC geometry with a higher temperature at the outer cylinder – a configuration distinct from ours.
Variations of the linear critical parameters as functions of radius ratio 
$\eta$ for the quasi-Keplerian flow at different 
$(Pr,Ri)$ settings: (a) linear critical Taylor number 
$Ta_c$, (b) linear critical axial wavenumber 
$k_c$, (c) linear critical azimuthal wavenumber 
$m_c$.
Panel 10(b) depicts the variation of the corresponding linear critical 
$k_c$, showing that the onset of instability is consistently triggered by disturbances with finite axial wavelength. This axial dependence suggests that the linear instability sets in with axial variation, at least when 
$Ta$ is slightly beyond 
$Ta_c$. Regarding the linear critical 
$m_c$, panel 10(c) illustrates that at small 
$\eta \lessapprox 0.1$ the instability onset is consistently dominated by the first helical mode at 
$m=1$, except for the case at 
$(Pr,Ri)=(0.01,0.01)$, which is dominated by the axisymmetric mode. At relatively large 
$\eta$, the case 
$(Pr,Ri)=(0.7,0.1)$ favours 
$m_c=2$, whereas all the other cases pick the axisymmetric mode.
4.2.5. A comparison to the local stability analysis of the flow
As suggested by one of the reviewers, we next compare the results of our stability analysis to those of a local analysis (Klahr & Hubbard Reference Klahr and Hubbard2014; Lyra Reference Lyra2014). In a local analysis, flow homogeneity has been assumed in all the directions. To distinguish, we refer to our above analysis as a global analysis, as the flow in the radial direction is inhomogeneous due to the presence of cylinder walls.
In the following, we select control parameters in the local analysis as close as possible to those used in §§ 4.2.1–4.2.4. The obtained local results also reveal a linear instability in the thermally stratified quasi-Keplerian flow. The observed effects of various governing parameters on this instability are also qualitatively consistent to those in the global analysis. Differences lie in two aspects. First, the local analysis based on the short wave approximation only provides limited information on the instability growth rate due to spatially periodic disturbances of short wavelength; it gives no information on the unstable mode pattern. In contrast, the global analysis uncovers rich mode patterns of the instability, exhibiting various length scales in all three directions; see figure 3, figure 14 and more in the supplementary material. Second, the instability onset Taylor number and axial wavenumber from the two analyses show notable differences; and so do the corresponding scaling laws. All the comparisons, either qualitative or quantitative, are summarised in table 1. Below, we describe them in detail.
Comparisons of the results from the global analysis and the local stability analysis of the same quasi-Keplerian flow in the TC geometry. For all the comparisons, 
$m=0$. For numbers 2 to 8 in the comparisons, the local results are obtained at a fixed radial location of 
$r_0=(r_i+r_o)/2$.
In all the calculations for the local analysis, solving the 
$5\times 5$ eigenvalue problem (2.13) results in five eigenvalues, from which the leading eigenvalue is always selected for extracting the linear growth rate. Figure 11(a) displays contours of the linear growth rate 
$\omega _i$ in the wavenumber plane at three different radial locations for the quasi-Keplerian flow at 
$(\eta,Pr,Ri,Ta,m)=(0.3,0.7,0.1,10^{10},0)$. The area below the neutral curve indicates linear instability obtained from the local analysis. The corresponding instability in the global analysis is reported in figure 4, where the largest growth rate is 
$\omega _i\approx 1.35$ at 
$(k,m)=(47.62,0)$. In figure 11(a) the largest growth rate of 
$\omega _i\approx 1.05$ is attained at 
$(k_r,k_z)\approx (44,67)$ for 
$r_0=0.5$. Quantitatively, both the growth rate and the axial wavelength are of the same order of magnitude. The radial location 
$r_0=0.5$ also aligns with the localised mode pattern depicted in figure 1(b) of the supplementary material. Another important observation is that for a given 
$k_z$, the growth rate increases as 
$k_r$ decreases, indicating that the local instability favours relatively large radial structures. Comparatively, the global analysis is better suited for resolving large-scale unstable modes.
Linear instability observed in the local analysis of the quasi-Keplerian flow at 
$(\eta,Pr,Ri,m)=(0.3,0.7,0.1,0)$. (a) Contours of the linear growth rate in the wavenumber plane for the flow at three different radial locations 
$r_0$ and 
$Ta=10^{10}$. The solid black curve represents the neutral curve. (b) Variation of the neutral curve in the wavenumber plane with 
$Ta$ for the local flow at the centre of the cylinder gap 
$r_0=(r_i+r_o)/2$. (c) Variation of the neutral curve in the 
$Ta$–
$k_z$ plane with varying radial locations 
$r_0$. (d) Variation of the neutral curve with Richardson number 
$Ri$ for the local flow at the centre of the cylinder gap 
$r_0=(r_i+r_o)/2$.
Figure 11(b) shows how the neutral curve varies with 
$Ta$ in the wavenumber plane for the local flow at the centre of the cylinder gap 
$r_0=(r_i+r_o)/2$. The expansion of the neutral curve with increasing 
$Ta$ indicates a destabilising effect, consistent with the global analysis reported in § 4.2.1. However, this variation is specific to the radial location 
$r_0=(r_i+r_o)/2$. To determine the critical 
$Ta$ that depends on 
$r_0$, we plot neutral curves in the 
$Ta$–
$k_z$ plane for various values of 
$r_0$ in figure 11(c). In this plot, the radial wavenumber 
$k_r$ is set to its lowest allowed value at each radial location 
$r_0$, i.e. 
$k_r=\max [{\rm \pi} /(r_0-r_i),{\rm \pi} /(r_o-r_0)]$ (note that for a given radial wavenumber 
$k_r$, the disturbance wavelength is 
$L_r=2{\rm \pi} /k_r$ and the whole wave must be confined within the radial domain, requiring that 
$r_0-L_r/2\ge r_i$ and 
$r_0+L_r/2\le r_o$). This choice of 
$k_r$ ensures the highest linear growth rate maximised over 
$k_r$. Figure 11(c) shows that 
$r_0=(r_i+r_o)/2$ corresponds to the smallest critical 
$Ta$ above which the local instability occurs, with an approximate value of 
$Ta_c\approx 5\times 10^{5}$. This value is of a similar order of magnitude to that observed in figure 5(b) from the global analysis.
In light of the above effect of 
$r_0$ on the local instability, we fix 
$r_0=(r_i+r_o)/2$ to study the influence of various governing parameters on the instability onset. Figure 11(d) illustrates the destabilising effect of increasing the Richardson number 
$Ri$ on the neutral curve, consistent with the global analysis results presented in figure 6(a). However, a notable difference is that the linear critical axial wavenumber 
$k_{z,c}$ remains constant across varying 
$Ri$; see also figure 12(d,e,f), where it maintains a value of approximately 
$4.44$ despite changes in 
$(\eta,Pr,Ri)$. The reason for this constant value remains unclear. Figure 12(a) shows 
$Ta_c$ as a function of 
$Ri$. It can be seen that the scaling law of 
$Ta_c$ has significantly changed, with an exponent of 
$-1$ compared with 
$-5$ in the global analysis. Increasing the Prandtl number 
$Pr$ is also destabilising, as evidenced by figure 12(b), where the scaling law is 
$Ta_c\propto Pr^{-1}$. This is in close agreement with the scaling 
$Ta_c\propto Pr^{-6/5}$ from the global analysis shown in figure 9(a). Additionally, figure 12(c) demonstrates that reducing the radius ratio 
$\eta$ generally has a destabilising effect, consistent with the trend observed in figure 10(a) from the global analysis.
Variations of the linear critical Taylor number 
$Ta_c$ and the linear critical axial wavenumber 
$k_{z,c}$ with the control parameters 
$(Ri,Pr,\eta )$ for the axisymmetric (
$m=0$) quasi-Keplerian local flow at the cylinder gap centre 
$r_0=(r_i+r_o)/2$. (a,d) At the low-Richardson-number 
$Ri$ limit and 
$(\eta,Pr)=(0.3,0.7)$. The limit values of 
$Ri_{{limit}}$ are computationally estimated to be 
$Ri_{{limit}}\approx 0.0453873$. (b,e) At the small-Prandtl-number 
$Pr$ limit and 
$(\eta,Ri)=(0.3,0.1)$. (c,f) At varying viscosity ratio 
$\eta$ and 
$(Pr,Ri)=(0.01,0.01)$. In panels (d–f), 
$k_c\approx 4.44$.
In a nutshell, the results from the local stability analysis of the quasi-Keplerian flow are qualitatively consistent with those from the global analysis, indicating that similar underlying physics is at play. However, quantitatively, the scaling laws of the linear critical conditions differ between the local and global analyses, highlighting the distinctiveness of these two approaches. Comparatively, the global analysis provides more solid results than the local analysis. The local analysis involves numerous assumptions and the large-scale mode structures resolved in the global analysis can never be reproduced in the local analysis, underscoring the limitations of the local approach.
4.3. Linear instability at extreme parameters
The preceding subsection demonstrates that, generally, decreasing 
$Pr$ or 
$Ri$ has stabilising effects, while decreasing 
$\eta$ is destabilising. To discuss the values of these parameters relevant to accretion disks in the literature, we mention that Held & Latter (Reference Held and Latter2018) used 
$Pr=2.5$ and 
$Ri=0.05$ in their numerical simulations of thermal quasi-Keplerian flows in a 3-D local shearing box. Latter (Reference Latter2016) noted that 
$Pr$ could be as low as 
$10^{-7}$ in the inner radii and the author adopted a small 
$Ri$ value of 
$0.01$. For a TC geometry, appropriate values of 
$\eta$ are not directly extractable from the literature due to the absence of an ‘outer rotating cylinder’ in real accretion disks. For modelling purposes, one may consider small values since disks are typically thin and flat with relatively small central objects (Abramowicz & Straub Reference Abramowicz and Straub2014). This subsection illustrates that the above reported thermally driven linear instability persists with the control parameters becoming as small as 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$, which are more astrophysically relevant. The linear instability at much more severe parameters is further demonstrated in § 5 of the supplementary material. In addition, the influence of the gravitational acceleration profile is examined at the parameter setting of 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$ and discussed in Appendix B, showing that a stronger gravity field enhances the linear instability and vice versa.
4.3.1. The critical condition for the flow at $(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$
The neutral curves in the case 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$ for different values of 
$m$ are shown in figure 13, from which one can see that the linear instability sets in at approximately 
$(Ta_c,k_c,m_c)=(2.28 \times 10^{9},2.46,1)$; see the red star in panel (b). The eigenvalue of this neutral mode at the critical condition is 
$\omega \approx -0.15607945375-0.00000000002\textrm {i}$. The frequency corresponds to a negative phase speed of 
$c_r=\omega _r/m\approx -0.156$, indicating that this mode rotates in the opposite direction of the base flow when viewed in the rotating frame of reference. After transforming back to the fixed frame of reference (see § 2 of the supplementary material for a description of the transformation method), its phase speed becomes 
$c_r=[\omega _r + m/(2Ro)]/m=\omega _r/m + 1/(2Ro)\approx 0.0587$, much smaller than the lowest angular frequency 
$\omega _o =1/(2Ro)\approx 0.214$ of the base flow, which is attained at the outer rotating cylinder. The pattern of this mode at 
$m=1$ is visualised in figure 14(a iii,b iii), showing that the disturbance fills the entire cylinder gap, with the largest variations occurring in the thin vicinity near the inner rotating cylinder. In the 
$r$–
$z$ plane the vortices attached to the inner cylinder are strongly tilted, indicating that the mode travels in the negative 
$z$ direction when viewed in the rotating frame of reference. This is consistent with its phase speed in the 
$z$ direction 
$c_r=\omega _r/k_c\approx -0.0634$.
Contours of the linear growth rate 
$\omega _{i}$ in the 
$Ta$–
$k$ plane for the quasi-Keplerian flow at 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$. The black curves are the neutral curves (
$\omega _i=0$) and the red star marks the linear critical condition at about 
$(Ta_c,k_c,m_c)=(2.28 \times 10^{9},2.46,1)$.
Contours of the perturbation velocity magnitude at the four critical points identified in figure 13. Panel (a) shows the cross-sections in the 
$r$–
$z$ plane where only one wavelength is plotted at the azimuthal angle 
$\theta =0$; panel (b) shows the cross-sections in the 
$r$–
$\theta$ plane at the axial position 
$z=0$. From left to right, (a i,a ii) and (b i,b ii) correspond to a pair of complex conjugate modes at 
$m=0$, (a iii–a v) and (b iii–b v) are for 
$m=1,2,3$, respectively.
In addition to the helical mode at 
$m=1$, the axisymmetric mode (
$m=0$) warrants examination due to the proximity of its neutral curve's lowest 
$Ta$ value (
$2.41 \times 10^{9}$) to 
$Ta_c\approx 2.28 \times 10^{9}$. This closeness suggests potential interaction between these modes in nonlinear simulations of flow transition at the instability onset. The neutral mode has a pair of complex conjugate eigenvalues 
$\omega \approx \pm 0.307448512+0.000000003\textrm {i}$. The structures of this pair of modes are visualised in figure 14(a i,a ii,b i,b ii). In the 
$r$–
$z$ plane, localised vortices are seen tilting up and down near the inner rotating cylinder, indicating that these two modes propagate at the same speed but in opposite directions. In the 
$r$–
$\theta$ plane the structure is identical when visualised at the axial location 
$z=0$, as shown in figure 14(b i,b ii). In addition, figure 14(a iv,a v,b iv,b v) depicts the flow structures of the non-axisymmetric modes 
$m=(2,3)$ at the identified critical points. As 
$m$ increases, the disturbance becomes more localised and moves closer to the outer rotating cylinder. Similar behaviour can be observed in figure 3 of the supplementary material. It should be noted that the parameters used here and those for figure 3 therein are significantly different, suggesting that the shifted localisation of the mode with increasing 
$m$ may be a generic feature across the parameter space.
4.3.2. Beyond the critical condition: linear instability at $Ta={10}^{16}$
With the instability onset identified at 
$(Ta_c,k_c,m_c)=(2.28 \times 10^{9},2.46,1)$, the following question arises: Is this critical condition relevant to accretion disks? A key factor in this consideration may be the axial wavelength of the mode. In this case, the axial wavelength is on the order of the cylinder gap, making it incompatible with accretion disks because the characteristic axial length scale in real disks is usually much smaller than the radial length scale, such as the disk radius (Latter Reference Latter2016). To obtain more relevant results, we have searched for linear instabilities with very small axial wavelengths. Indeed, figure 4 demonstrates that linear instability persists at Taylor numbers significantly higher than the onset value (
$Ta\gg Ta_c$). At such high 
$Ta$, the linear instability spans a wide range of high axial wavenumbers, corresponding to small axial length scales, which are more similar to the accretion-disk conditions.
For the case at 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$, we investigate the unstable mode at 
$Ta=10^{16}$, seven orders of magnitude larger than the critical value of about 
$10^9$. The highest 
$Ta$ previously investigated in TC flows is about 
$10^{13}$ (Grossmann et al. Reference Grossmann, Lohse and Sun2016). Although 
$Ta=10^{16}$ may be challenging to achieve in numerical simulations and physical experiments, it is feasible in our linear stability analysis provided sufficient spatial resolution in the spectral method. For the following calculations, we use 
$401$ Chebyshev–Lobatto nodes. Figure 15(a) shows the dispersion relation 
$\omega _i(k)$ for the case abovementioned. At this 
$Ta$, the most unstable mode becomes the axisymmetric mode (
$m=0$), rather than the helical mode at 
$m=1$ that dominates at the linear critical condition, as discussed in figure 13(b). Notably, the flow is linearly unstable over a wide spectrum of axial wavenumbers 
$k$. For 
$m=0$, the instability range is 
$0< k<800$; for 
$m=1$, the range narrows but remains broad at 
$0\le k<300$. Higher azimuthal wavenumbers 
$m>2$ also exhibit linear instability, but with much smaller linear growth rates, so they are not shown. The non-smoothness of the blue curve for 
$m=0$ in figure 15(a) at around 
$k\approx 250$ is due to the transition of the most unstable mode from stationary to oscillatory, as seen in the frequency variations 
$\omega _r(k)$ in panel (b). Lastly, it should be mentioned that in addition to the above most unstable mode at 
$m=0$ and the leading unstable mode at 
$m=1$, there are a magnitude of axisymmetric and non-axisymmetric unstable modes in the corresponding eigenspectra for this case, signifying the high complexity of the instability diagram. As these modes are less unstable, we will not characterise them further.
(a) Linear growth rate 
$\omega _i$ as function of 
$k$ for the quasi-Keplerian flow at 
$Ta=10^{16}, \eta =0.05$, 
$Pr=10^{-3}$ and 
$Ri=10^{-3}$. (b) Frequency 
$\omega _r$ as function of 
$k$ for the same flow.
5. Conclusions
Motivated by research efforts searching for hydrodynamic instabilities without magnetic effects that may be operative in protoplanetary disks, we conducted a 3-D ‘global’ analysis of thermally stratified quasi-Keplerian flows between a hot inner rotating cylinder and a cold outer rotating cylinder.
Unlike previous work by Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021), our model incorporates stellar gravity in the radial direction and considers the temperature distribution within accretion disks, which decreases with increasing distance from the central star. The flow instability in the TC geometry should also be distinguished from that found based on short wave approximations (Klahr & Hubbard Reference Klahr and Hubbard2014; Lyra Reference Lyra2014), as the latter bears only local significance.
We identified a robust flow instability in the quasi-Keplerian flow regime driven purely by radial thermal stratification along with radial gravity, without the need for a magnetic field as required by MRI. This instability manifests in both infinitesimal axisymmetric and non-axisymmetric modes across two or three dimensions, exhibiting typical unstable eigenmode structures resembling those in canonical Rayleigh–Bénard convection but distorted by the laminar circular Couette flow, especially in non-axisymmetric cases. Our calculation of the linear critical conditions (
$Ta_c,k_c,m_c$) revealed that decreasing the Prandtl number 
$Pr$ or Richardson number 
$Ri$ generally stabilises the quasi-Keplerian flow. Conversely, reducing the radius ratio 
$\eta$ is destabilising. At small values of 
$Pr=0.001$ and 
$Ri=0.001$, the linear instability persists at sufficiently large 
$Ta\gtrsim 10^9$ with a small axial wavelength in the weakly thermally stratified quasi-Keplerian flow. This suggests that even very weak thermal buoyancy is sufficient to trigger the linear instabilities. Likewise, Balbus & Hawley (Reference Balbus and Hawley1991) showed that even a weak magnetic field suffices to cause MRI. Additionally, we revealed two scaling laws. The first one shows that 
$Ta_c\propto Pr^{-6/5}$ at the small-
$Pr$ limit, implying again the relevance to accretion disks where 
$Pr$ can be as low as 
$10^{-7}$ (Latter Reference Latter2016). The second one reveals that 
$Ta_c\propto (Ri-Ri_{{limit}})^{-5}$, indicating the stabilising effect of decreasing 
$Ri$. These scaling laws may help to guide future high-fidelity numerical investigations or experiments of the studied flow.
This thermally driven hydrodynamic linear instability and subsequent nonlinear dynamics may provide a mechanism for explaining radially outward angular momentum transport in unmagnetised accretion disks or dead zones where MRI is inoperative. We note that Lyra (Reference Lyra2014) has successfully produced self-sustained 3-D vortex flows from the saturation of a linear overstability in the modelling of protoplanetary disks; however, their linear instability is localised since they assumed homogenisation in three directions and the evolution may need to be further investigated in a global sense. In addition, Held & Latter (Reference Held and Latter2018) have demonstrated a series of accretion-disk dynamics from the linear instability (of the axisymmetric mode only) to turbulence based on a shearing box approximation of the quasi-Keplerian flow also with thermal stratification but in the axial direction. Unlike these previous localised linear instability assumptions, our 3-D global analysis complements existing discussions, furthering our understanding of accretion-disk dynamics by showing that even very weak radial thermal stratification can cause sufficiently strong gravitational buoyancy to drive the quasi-Keplerian flow linearly unstable. Meanwhile, we point out that the TC flow is a limited model for the accretion disk. One of the drawbacks is the configuration of two no-slip no-penetration cylinder walls that do not exist in real accretion disks. This work focuses on the momentum transport mechanism in accretion disks. Nevertheless, we hope that the present work on the radial thermal stratification could inspire more studies along this vein of research and contribute to a deeper understanding of the complex accretion-disk dynamics. Future research avenues may include investigating flow bifurcation beyond the linear instability, conducting direct numerical simulations of nonlinear development, exploring turbulent flows in thermally driven TC flow and assessing the effects of fluid compressibility on linear instability.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2024.976.
Acknowledgements
We thank all the reviewers for their constructive suggestions.
Funding
D.W. was supported by a PhD scholarship (no. 201906220200) from the China Scholarship Council and an NUS research scholarship. M.Z. acknowledges the financial support of NUS (Suzhou) Research Institute and National Natural Science Foundation of China (grant no. 12202300) and the MOE Tier 1 grant A-8001172-00-00. We also acknowledge the financial support from the Max-Planck Society, the Deutsche Forschungsgemeinschaft through grants 521319293 and 540422505, and the Daimler and Benz Foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Illustration of key differences between the current work and two previous studies on the same topic
In the introduction section we briefly mentioned the links of our study to some relevant prior works. Below, we describe in detail the differences between our study and Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) and Klahr & Hubbard (Reference Klahr and Hubbard2014) to demonstrate the position and novelty of our research.
Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) similarly conducted the linear stability analysis of the quasi-Keplerian flow in a TC geometry. However, their flow configuration, with a cold inner cylinder and a hot outer cylinder, does not align with the temperature distribution of accretion disks, where the temperature decreases in the radial direction (Fromang & Lesur Reference Fromang and Lesur2019). Additionally, their focus was on the effects of centrifugal buoyancy, neglecting gravitational acceleration, which may not be directly applicable for modelling accretion disks under stellar gravity. Moreover, the parameters most relevant to disk dynamics comprise low 
$Pr$ (Latter Reference Latter2016) and high 
$Ta$ (Ji et al. Reference Ji, Burin, Schartman and Goodman2006; Grossmann et al. Reference Grossmann, Lohse and Sun2016), which is at variance with the high Prandtl number and low Taylor number considered in Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021). The present study, however, takes these perspectives into account. More detailed comparisons are shown in table 2. Despite the differences mentioned, we acknowledge that the general destabilising mechanism in our study and that in Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021) are similar: density stratification due to a thermal gradient can be unstable under acceleration from external forces when the acceleration is in the opposite direction to the density stratification. Nevertheless, we believe that our study adds value to the explanation of angular momentum transport in accretion disks with more detailed results.
Differences and similarities between the present study and Meyer et al. (Reference Meyer, Mutabazi and Yoshikawa2021).
Regarding the comparison with Klahr & Hubbard (Reference Klahr and Hubbard2014), there are three main differences. First, they modelled the thermal relaxation process, where the heating and cooling time scales play a crucial role, in their discovery of the convective overstability. This is more realistic concerning the thermodynamics in accretion disks than our thermal equilibrium assumption, which was adopted for simplicity. Second, their analysis was local, performed in a small volume around a given spatial location 
$(r,\theta )$ (note that they ignored the axial velocity 
$u_z$ equation). This local assumption excludes the possibility of any linear instability triggered by large-scale structures. In contrast, our model is ‘global’, with flow boundaries in the radial direction and the entire 
$2{\rm \pi}$ radian in the azimuthal direction. Third, they made various assumptions and simplifications in deriving the dispersion relation due to the local model selection. In contrast, we numerically solve a generalised eigenvalue problem for the complete dispersion relation. Table 3 shows more detailed comparisons.
Differences and similarities between this study and Klahr & Hubbard (Reference Klahr and Hubbard2014).
In summary, both of the two abovementioned works are relevant to accretion disks, though the modelling levels differ. Subsequent works on the convective overstability by Lyra (Reference Lyra2014), Latter (Reference Latter2016) and Held & Latter (Reference Held and Latter2018) followed the same local approximation in a shearing box, so comparisons will not be repeated here.
Appendix B. Influence of the gravitational acceleration profile on the linear instability
In the present study the gravitational acceleration profile 
$g(r)=r_o/r$ is suited for the TC flow geometry described in cylindrical coordinates. However, for realistic accretion disks described in spherical coordinates, the profile due to the central object should follow 
$g(r)=r_o^2/r^2$. To the best of our knowledge, neither of these profiles has been realisable in laboratory experiments to date. This implies that numerical simulations are likely the most viable approach to further investigate the complex nonlinear dynamics of the observed linear instability. Fortunately, a constant profile 
$g(r)=1$ (non-dimensional) has been realised in experiments on thermo-electric flow convection in a cylindrical annulus (Antoine et al. Reference Antoine, Martin, Vasyl and Christoph2023). This ‘electric gravity’ can be generated due to dielectrophoretic forces when an electric field is applied radially. To minimise the impact of Earth's gravity on the experimental set-up, these experiments can be conducted during sounding rocket flights, which provide a microgravity environment (Antoine et al. Reference Antoine, Martin, Vasyl and Christoph2023). Although the ‘electric gravity’ technique is expensive and challenging, it offers a potential avenue for investigating the present thermally driven linear instability in quasi-Keplerian flows. Future work in this direction could benefit from linear stability analysis using the constant profile 
$g(r)=1$. As suggested by one of the reviewers, in this appendix we compare the effects of the above three different gravitational acceleration profiles on the linear instability.
These profiles can be expressed collectively as
\begin{equation} g(r;p) = \frac{r_o^p}{r^p} \quad \text{with} \ p=0,1,2. \end{equation}
At the outer rotating cylinder 
$r_o$, the non-dimensional gravitational acceleration is always equal to one, as 
$g_o^*$ is chosen as the reference value for non-dimensionalisation. For the flow field at 
$r< r_o$, the gravitational acceleration 
$g(r)\ge 1$. Its maximum value is attained at the inner rotating cylinder 
$r_i$, with 
$g(r_i;p=2)>g(r_i;p=1)>g(r_i;p=0)=1$.
Our calculations indicate that, with all the other parameters fixed, compared with the case with 
$p=1$, a stronger gravity field (
$p=2$) is destabilising, while a weaker gravity field (
$p=0$) is stabilising. This behaviour is illustrated by the eigenspectra shown in figure 16. Panel 16(a) displays the axisymmetric linear instability for the flow at 
$(\eta,Pr,Ri,Ta)=(0.05,10^{-3},10^{-3},10^{16})$. The leading unstable modes for 
$p=0,1,2$ are labelled as ‘LMp0a’, ‘LMp1a’ and ‘LMp2a’, respectively. Compared with mode ‘LMp1a’, mode ‘LMp2a’ exhibits a higher linear growth rate 
$\omega _i$, while mode ‘LMp0a’ shows a lower growth rate. A similar trend is observed for non-axisymmetric modes in panel 16(b). Neutral curves provide another perspective on the destabilising or stabilising effects. As shown in figure 17(a,b), the neutral curve for 
$p=2$ shifts to the left, indicating a smaller critical Taylor number and, thus, a destabilising effect due to the stronger gravity field. Conversely, for the constant profile (
$p=0$), the neutral curve shifts to the right and results in a larger critical Taylor number, demonstrating the stabilising effect of a weaker gravity field. It is worth noting that the linear critical Taylor numbers 
$Ta_c$ for the three gravitational acceleration profiles are of the same order of magnitude. This suggests that the constant gravitational acceleration profile with 
$p=0$, which is the only realisable profile in experiments to date, could be used to predict approximately the instability onset of the flow with 
$p=1$ in future experiments.
Eigenspectra for different gravitational acceleration profiles (
$g(r)=r_o^p/r^p$ with 
$p=0,1,2$, respectively) for the quasi-Keplerian flow at 
$(\eta,Pr,Ri,Ta)=(0.05,10^{-3},10^{-3},10^{16})$. Panel (a) is for the case at 
$(k,m)=(75.3178,0)$; panel (b) is for the case at 
$(k,m)=(80.4170,1)$; these two panels correspond to the two maxima respectively in figure 15(a). Labels point to the leading modes for each 
$p$.
Neutral curves for different gravitational acceleration profiles (
$g(r)=r_o^p/r^p$ with 
$p=0,1,2$, respectively) for the quasi-Keplerian flow at 
$(\eta,Pr,Ri)=(0.05,10^{-3},10^{-3})$ and (a) 
$m=0$, (b) 
$m=1$. The neutral curves for 
$p=1$ are exactly those in figure 13(a,b).
In the end, we mention in passing that the influence of the gravitational acceleration profile on the linear instability can be understood in terms of buoyancy. Increased gravity enhances the buoyancy force, which promotes flow instability, as evidenced by the gravitational buoyancy term in the governing equation (2.1a).
Open access
