1. Introduction
Taylor–Couette flow (TCF) refers to the flow within the annulus between two independently rotating coaxial cylinders (Couette Reference Couette1888; Mallock Reference Mallock1888; Taylor Reference Taylor1923). For small values of inner-cylinder rotation speed (
$\omega _i$), the bulk flow (away from the two endwalls) is close to purely azimuthal circular Couette flow (CCF), which, on increasing 
$\omega _i$, gives birth to a non-azimuthal flow consisting of an array of stationary vortices stacked along the axial direction, called the Taylor-vortex flow (TVF). Rayleigh's inviscid stability criterion (Rayleigh Reference Rayleigh1917) explained that an inviscid rotating flow is unstable to axisymmetric perturbations if the angular momentum of the fluid decreases outwards from the axis of rotation – this implies that the CCF is unstable for any non-zero rotation rate of the inner cylinder, with the outer one being at rest. Taylor (Reference Taylor1923) presented a linear stability analysis, accounting for viscosity and backed by experimentally observed data, which predicts the primary instability at a finite rotation rate of the inner cylinder, confirming the stabilizing role of viscosity on the instability onset.
The comprehensive experimental studies of Coles (Reference Coles1965) and Andereck, Liu & Swinney (Reference Andereck, Liu and Swinney1986) established that a plethora of patterns, like stationary Taylor vortices, azimuthally propagating wavy Taylor vortices (WTVs), axially azimuthally propagating spiral vortices, azimuthally propagating ribbons (or interpenetrating spirals), corkscrews and ripples, including the non-uniqueness of certain time-dependent patterns (Coles Reference Coles1965), are possible beyond the onset of the primary instability. The experimental studies of Benjamin (Reference Benjamin1978a,Reference Benjaminb) revealed that a multiplicity of steady states and an odd number of Taylor rolls are possible in this configuration and that the conditions at the endwalls (Schaeffer Reference Schaeffer1980; Benjamin & Mullin Reference Benjamin and Mullin1981; Jones Reference Jones1982; Cliffe Reference Cliffe1983; Lücke et al. Reference Lücke, Mihelcic, Wingerath and Pfister1984; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989) play a crucial role in the development of patterns in finite-length cylinders.
In small-aspect-ratio (
$\varGamma ={\textit{O}}(1)$) Taylor–Couette (TC) cells, Benjamin & Mullin (Reference Benjamin and Mullin1981) showed that the two-roll state bifurcates into a pair of single-roll states with increasing 
$Re$. This corresponds to 
$Z_2$-symmetry-breaking pitchfork bifurcation, since the bifurcated state breaks the midplane axial symmetry. Such asymmetric two-roll states, with one large vortex coexisting with a smaller vortex, are often called ‘single-roll’ states. The latter terminology was coined by Cliffe (Reference Cliffe1983) since the asymmetric two-roll state was found to degenerate into a single-vortex state in computations of an axially periodic model. The single-roll states have been extensively probed in subsequent studies (Lücke et al. Reference Lücke, Mihelcic, Wingerath and Pfister1984; Pfister et al. Reference Pfister, Schmidt, Cliffe and Mullin1988; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989; Furukawa et al. Reference Furukawa, Watanabe, Toya and Nakamura2002; Mullin, Toya & Tavener Reference Mullin, Toya and Tavener2002). In particular, the experimental and numerical work of Mullin et al. (Reference Mullin, Toya and Tavener2002) demonstrated that the symmetry-breaking bifurcation sequence coalesce into a symmetric two-roll state with decreasing 
$\varGamma$, implying that the midplane symmetric steady states represent unique states for very short cylinders (
$\varGamma <1$) with stationary (no-slip) endwalls. The above studies belong to near-critical steady-state flows, for which the Reynolds number is of moderate value, 
$Re < {\textit{O}}(10^3)$, an order of magnitude lower than that where turbulence is expected to set in. Interestingly, a recent experimental study (Huisman et al. Reference Huisman, van Der Veen, Sun and Lohse2014) demonstrated that the multiplicity of flow structures can persist even in the turbulent regime of TCF at very high Reynolds numbers (
${\sim }10^{6}$).
While Benjamin (Reference Benjamin1978a,Reference Benjaminb) discovered that steady states with an odd number of rolls could exist with symmetric axial boundary conditions, such asymmetric Taylor vortices can be expected when the endwall conditions are asymmetric. The axial asymmetry can be imposed either by having a free surface at the top (Cole Reference Cole1976; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989) or by having a rotating endplate. For example, Mullin & Blohm (Reference Mullin and Blohm2001) studied the bifurcation phenomena in a TC set-up with one endplate rotating with the inner cylinder while the other endplate is held stationary. They probed the selection process of the preferred steady state using both experiments and numerical bifurcation analysis of the Navier–Stokes equations. They found three vortices at 
$\varGamma =O(3)$ for smaller values of 
$Re$, which collapsed into a single vortex at some critical value of 
$Re$ via a saddle-node bifurcation. For the latter case, the vortex near the stationary endwall was found to grow in size with increasing 
$Re$, eventually replacing the three-vortex state with a large single vortex. They further showed that the three-vortex state can be recovered by reducing 
$Re$ in their numerical model, but this process was interrupted by the appearance of a three-dimensional time-dependent state in experiments.
Using both experiments and numerical analysis, Tavener, Mullin & Cliffe (Reference Tavener, Mullin and Cliffe1991) investigated the TCF, with two endplates that are rotated along with the inner cylinder. The rotating endwalls produce a strong driving effect on two vortices adjacent to the walls; the rotating endwall-induced driving increases with increasing Reynolds number. This added driving (compared to the stationary endwalls) makes the two-roll configuration more stable, but the higher-roll configurations (with 
$\varGamma > 2$) tend to lose stability with increasing 
$Re$ since the two end rolls tend to grow and squeeze out their interior neighbours. They identified several codimension-two bifurcation points with interesting dynamical behaviour. Another noteworthy finding of Tavener et al. (Reference Tavener, Mullin and Cliffe1991) is that the rotating endwalls readily admit asymmetric Taylor vortices more so than the stationary endwalls.
The above works pertain to incompressible TCF. Prior studies on the effect of fluid compressibility (Kao & Chow Reference Kao and Chow1992; Welsh, Kersal'e & Jones Reference Welsh, Kersal'e and Jones2014) and rarefaction (Yoshida & Aoki Reference Yoshida and Aoki2006; Manela & Frankel Reference Manela and Frankel2007) dealt with linear stability of purely azimuthal CCF of a dilute/ideal gas, but neither the roles of axial boundary conditions nor nonlinear simulations were carried out to address (i) the bifurcation structure, (ii) the resulting nonlinear (non-CCF) states and (iii) the anomalous modes in compressible TCF. The present work deals with molecular dynamics (MD) simulations of the TCF of a compressible ‘dense’ gas, with the primary focus to understand symmetry-breaking bifurcations and hysteresis in this flow. The effectiveness of MD simulations to capture the transition from CCF to TVF was demonstrated first by Hirshfeld & Rapaport (Reference Hirshfeld and Rapaport1998); the variation of radial velocity over a range of Taylor numbers was shown to be consistent with the theoretical predictions of critical Taylor number from continuum hydrodynamics and laboratory experiments. Building upon this, the torque scaling with angular velocity and the rate of growth and decay of Taylor vortices were analysed subsequently (Hirshfeld & Rapaport Reference Hirshfeld and Rapaport2000). Travelling waves in the form of WTVs were found (Trevelyan & Zaki Reference Trevelyan and Zaki2016) in recent MD simulations of a narrow-gap (of radius ratio of 
$0.873$) TC cell.
In an altogether different context, Taylor-like vortices have also been uncovered in MD simulations of both dense (Krishnaraj & Nott Reference Krishnaraj and Nott2016) and moderately dense (Mahajan & Alam Reference Mahajan and Alam2016) granular flows – while the former employed the soft-particle discrete-element method, the latter used the traditional inelastic hard-sphere model for the simulation technique. The compressibility effects (Savage Reference Savage1988; Krishnaraj & Nott Reference Krishnaraj and Nott2016) are important in fluidized particulate flows in which the Taylor-like vortices have been realized in experiments too (Conway, Shinbrot & Glasser Reference Conway, Shinbrot and Glasser2004).
More recently, experiments on suspension TCF (Majji, Banerjee & Morris Reference Majji, Banerjee and Morris2018; Ramesh, Bharadwaj & Alam Reference Ramesh, Bharadwaj and Alam2019; Ramesh & Alam Reference Ramesh and Alam2020) have unveiled a variety of new patterns in addition to WTVs. Ramesh et al. (Reference Ramesh, Bharadwaj and Alam2019) showed that the radial segregation of particles led to shear-band-type profiles for the azimuthal velocity, which are responsible for the non-monotonic variations of the angular momentum along the radial direction. The latter finding implies that the compressibility of the particle phase in suspension TCF can lead to a scenario in which a part of the flow is centrifugally stable or unstable (Rayleigh Reference Rayleigh1917) – a theoretical understanding of this issue remains elusive due to the unavailability of universal constitutive models for suspensions. In any case, the underlying issues can be systematically addressed via a bottom-up approach in which the simpler flow configurations are first studied in detail before a fluid–particle system is taken up for study.
In this work, the compressible TCF is analysed via MD simulations with the Lennard-Jones potential with a normal damping force being used as a fluid model. Details of the simulation technique, the particle-level boundary conditions and the statistical steady state of the TCF are described in § 2. Our major objectives are to (i) quantify the primary bifurcation from CCF to TVF, (ii) identify the role of fluid compressibility on Taylor vortices and (iii) analyse the secondary bifurcation from the TVF state with increasing rotation. The secondary bifurcation is found to break the 
$Z_2$-symmetry, leading to asymmetric vortices and a hysteresis phenomenon that are analysed in detail. The rest of this paper is organized as follows. The primary bifurcation from CCF to TVF is analysed and the procedure to determine the critical rotation rate is described in § 3.1; the observed patterns are presented in terms of a phase diagram as a function of (
$\omega _i, \varGamma$) in § 3.2. The structural characteristics of compressible Taylor vortices are deciphered by analysing the hydrodynamic fields and related order parameters in § 4. The onset of secondary bifurcation from the TVF state and the emergence of multiple solutions for specified (
$\omega _i, \varGamma$) are analysed in terms of bifurcation diagrams in § 4.1.1; the patterns in small-aspect-ratio cylinders (
$\varGamma \leq 1$) are analysed in § 4.2. The role of stationary endwalls and the underlying bifurcation scenario are discussed in § 5.1; the role of endwall conditions on the single-vortex mode and the emergence of unsteady patterns are briefly discussed in § 5.2. A summary of all results with conclusions and possible future works are given in § 6. While the results in the main text are presented for a non-zero value of the particle-level damping force (see § 2), the generality of our findings is demonstrated in appendix A by setting the microscopic damping force to zero.
2. Molecular dynamics simulation of Taylor–Couette flow
The schematic of the simulation domain is shown in figure 1: it is described by an annular space of thickness 
$\delta =(R_o -R_i)$ between two concentric cylinders of equal height 
$h$, with outer and inner radii of 
$R_o$ and 
$R_i$, respectively. The outer cylinder is kept stationary (
$\omega _o=0$) and the inner cylinder rotates with a constant rotational speed (
$\omega _i \neq 0$). The soft-sphere fluid particles interact through a short-range repulsive Lennard-Jones (LJ) potential 
$V(r)$, which results in a force of the form
\begin{equation} F_e(r)= -\frac{\textrm{d}V(r)}{\textrm{d}r} = 24 \left(\frac{\epsilon}{r}\right) \left[ 2 \left( {\frac{\sigma}{r}}\right)^{12} - \left( {\frac{\sigma}{r}}\right)^{6} \right], \end{equation}
where 
$r$ is the distance between two particle centres, 
$\sigma$ is the distance at which the potential is zero (representing the effective diameter of the particles) and 
$\epsilon$ is the depth of the potential well. The simulations discussed herein have 
$\sigma =\epsilon =1$ (reduced LJ units) and a cutoff at 
$r_c=2^{1/6} \sigma$, which refers to the ‘repulsive’ part of the LJ potential. It may be noted that removing the attractive part and retaining only the repulsive part of the LJ potential provides a model for a hard-sphere gas. In addition to this, a normal damping force acts over the duration of each collision and is given by 
$F_d= - m\gamma _n v_n$ (Rapaport Reference Rapaport2004), where 
$v_n$ is the normal component of the relative velocity of two colliding particles, 
$m$ denotes the mass of a particle and 
$\gamma _n$ is the damping coefficient. This particle-level damping is akin to (i) a ‘velocity-dependent’ normal restitution coefficient in a granular gas (Goldhirsch Reference Goldhirsch2003; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004) and/or (ii) a drag force in a gas–solid suspension (Jackson Reference Jackson2000; Saha & Alam Reference Saha and Alam2017) at large Stokes numbers; the specific role of this damping mechanism on both primary and secondary bifurcations and the emerging patterns are discussed in appendix A.
Figure 1. 
$(a)$ Schematic and the domain binning of the TC cell; 
$h$ is the height of the cylinders, 
$R_i$ and 
$R_o$ are the inner and outer radii, and 
$\delta =R_o-R_i$ is the gap width; the rotational speeds of the inner and outer cylinders are 
$\omega _i$ and 
$\omega _o=0$, respectively. 
$(b)$ Schematic of wall conditions for the simulation system.
2.1. Boundary conditions
As in the work of Hirshfeld & Rapaport (Reference Hirshfeld and Rapaport1998), the inner and outer cylinders are treated as thermal walls with a wall temperature 
$T_w$. Any particle that collides with the inner or outer cylindrical wall has all memory of its velocity erased and is assigned a random velocity in addition to a mean velocity corresponding to the tangential velocity of the cylinder (
$v_\theta =\omega _iR_i$). Referring to figure 1
$(b)$ and using the transformation rule between Cartesian (
$x,y,z$) and polar (
$r, \theta , z$) coordinates, a reflected particle from the cylindrical walls is assigned the following velocity:
\begin{equation} (c'_x, c'_y, c'_z )= \sqrt{\frac{k_B T_w}{m}} \,(v_{Gx}, v_{Gy}, v_{Gz} ) + \omega_iR_i(-\sin\theta, \cos\theta, 0), \end{equation}
where 
$v_{Gi}$ is sampled from the Gaussian distribution
\begin{equation} p(v_{Gi}) = \left(\frac{m}{2{\rm \pi} k_BT_w}\right)^{1/2}\exp\left(- \frac{mv_{Gi}^2}{2k_B T_w}\right), \end{equation}
with 
$k_B$ being the Boltzmann constant. If the assigned velocity is such that the particle leaves the simulation domain, then the directional sense of 
$v_{Gx}$ and 
$v_{Gy}$ is inverted (
$v_{Gx}=-v_{Gx}$ and 
$v_{Gy}=-v_{Gy}$).
The axial boundary condition is implemented using the following rule: any particle that moves outside the top or bottom endwalls, by a distance 
$s_z$ during time integration, is put back into the box by the same distance and the 
$z$-component of the velocity of that particle is inverted such that
\begin{equation} \boldsymbol{c}' \boldsymbol{\cdot} \hat{\boldsymbol{k}} = - \boldsymbol{c} \boldsymbol{\cdot} \hat{\boldsymbol{k}}, \end{equation}
where 
$\hat {\boldsymbol {k}}$ is the unit vector along the axial direction; note that the reflecting particle's 
$x$- and 
$y$-component velocities remain unchanged.
The above thermal-wall condition (2.2) and reflective endwall condition (2.4) are similar to those employed in the previous MD simulations (Hirshfeld & Rapaport Reference Hirshfeld and Rapaport1998; Trevelyan & Zaki Reference Trevelyan and Zaki2016) of the TCF; similar boundary conditions were also used by Stefanov & Cercignani (Reference Stefanov and Cercignani1993) to study the onset of vortices in Monte Carlo simulations of the TCF of a rarefied gas. The role of stationary endwalls, resulting in no-slip velocity, will be briefly analysed for few specific cases in § 5.
2.2. Simulation methodology with CPU and GPU
The open-source LAMMPS MD code (Plimpton Reference Plimpton1995) is employed to integrate the equations of motion in the Cartesian coordinate system using a velocity–Verlet integrator. To speed up computations, LAMMPS has been built from source with the graphics processing unit (GPU) package enabled and optimized for performance on a GPU with 3072 CUDA (Compute Unified Device Architecture) cores and 12 GB GPU random-access memory (RAM). The central processing unit (CPU) computational component consists of a system total of 64 GB CPU RAM equally distributed amongst two sockets, each housing a 2.4 GHz Haswell architecture processor having eight cores, of which 14 (
$2\times 7$) cores are employed for computations. The average scale-up obtained by this ‘CPU + GPU’ combination for the family of runs corresponding to 
$\varGamma =2$ and 
$\delta =25 \sigma$ is found to be approximately 
$2.3$ times that of the CPU-only runs employing 14 cores. The scale-up performance is seen to decay with both increasing 
$\omega _i$ and increasing number of particles. The latter can be attributed to the more frequent neighbour list rebuilds, which can occur in the CPU only, owing to the nature of the potential used.
All simulations have been carried out for a mean reduced density of 
$\rho _{av}=0.5$ as in previous studies on TCF (Hirshfeld & Rapaport Reference Hirshfeld and Rapaport1998; Trevelyan & Zaki Reference Trevelyan and Zaki2016); the aspect ratio (
$\varGamma = h/\delta = 0.5, 0.75, 1.0, 1.33, 2, 4, 10$) was varied by changing the height (
$h$) of the cylinders. In most simulations, the outer and inner cylinder radii are taken as 
$R_o=75\sigma$ and 
$R_i=50\sigma$, such that the gap width is 
$\delta =(R_o-R_i)=25\sigma$, the radius ratio is 
$\eta =R_i/R_o=2/3$ and the curvature is 
$\kappa =\delta /R_i=1/2$. The latter two quantities (
$\eta , \kappa$) suggest that our simulations correspond to ‘wide-gap’ TC geometry. The total number of particles within the annular gap is calculated from
\begin{equation} N = \rho_{av}{\rm \pi} \left(\frac{R_o^2}{\sigma^2}\right) (1-\eta^2 ) \varGamma\left(\frac{\delta}{\sigma}\right), \end{equation}
which varied from 
$60\,000$ to 
$1.25\times 10^6$, depending on the value of 
$\varGamma \in (0.5,10)$. At smaller values of 
$\varGamma < 1$, a four-fold increase of 
$N$ (by keeping 
$\eta$, 
$\kappa$ and 
$\varGamma$ fixed) did not affect the observed patterns.
The simulation procedure consists of two stages: (i) an equilibration run that consist of 
$10^{5}$ time steps during which the inner cylinder is uniformly accelerated from rest to the required angular velocity and (ii) a production run consisting of another 
$2\times 10^6$ time steps during which the (instantaneous) hydrodynamic fields are calculated in each bin at every step, and then time-averaging over these bin-averaged properties is carried out over every 
$2000$ time steps as detailed in § 2.3. Most simulations have been carried out with a time step of 
${\rm \Delta} t=0.005$; a few simulations were checked with smaller time steps of 
${\rm \Delta} t=0.0025$ (at 
$\varGamma =2$ and 
$4$) and 
$0.001$ (at 
$\varGamma =4/3$), which did not alter the results.
2.3. Hydrodynamic fields and the statistical steady state of TVF
To calculate hydrodynamic fields from the instantaneous particle positions and velocities, the simulation domain is binned along the radial and axial directions with bin sizes of 
$\textrm {d}z= \textrm {d}r =1.25\sigma$. Since we will be presenting results on CCF and TVF, both of which being azimuthally invariant states, the hydrodynamic fields will be calculated by averaging particle data over all 
$\theta$. Note that the mapping from the Cartesian 
$(x,y,z)$ to the polar 
$(r, \theta , z)$ coordinate system is done via the planar rotation matrix to obtain the velocity components of a particle (
$c_r, c_\theta , c_z$) in polar coordinates.
The instantaneous ‘reduced’ density is defined as 
$\rho ={\sigma ^3 N(t)}/{V_{bin}}$, where 
$V_{bin}$ is the volume of the bin and 
$N(t)$ is the instantaneous number of particles in a bin; the reduced density is related to the particle volume fraction (
$\phi$) via 
$\phi =({\rm \pi} /6)\rho$. For a bin in the (
$r, z$) plane, the local reduced density is calculated from
\begin{equation} \rho(r,z)=\frac{{\sigma^3}}{{\rm \pi}(R_{bo}^2 - R_{bi}^2)\,\textrm{d}z} \left\langle\sum_{i=1}^{N}{\delta(r-r_i)\delta(z-z_i)}\right\rangle_\theta , \end{equation}
where 
$R_{bo}$ and 
$R_{bi}$ denote the outer and inner radii of the bin whose axial height is 
$\textrm {d}z$, and (
$r_i, z_i$) denote the coordinates of all particles residing in this particular bin. The subscript 
$\theta$ on the angular brackets in (2.6) implies a spatial averaging over all particles in the 
$\theta$-direction, in addition to temporal averaging over a few thousand snapshots of the system in the steady state.
The instantaneous velocity of a group of particles within a bin, at any instant in time, is defined as 
${\sum _{i=1}^{N}\boldsymbol {c}}/{N(t)}$, where 
$\boldsymbol {c}$ denotes the instantaneous velocity vector of the particles in that bin. Thus, the local hydrodynamic velocity is calculated from
\begin{equation} \boldsymbol{v}(r,z)=\frac{\left\langle {\sum_{i=1}^{N} \boldsymbol{c_{i}}(t) \delta(r-r_i)\delta(z-z_i)} \right\rangle_\theta}{\left\langle\sum_{i=1}^{N}{\delta(r-r_i) \delta(z-z_i)}\right\rangle_\theta} , \end{equation}
where 
$\boldsymbol {c}_i$ is the instantaneous velocity of the 
$i$th particle in the bin, and the angular brackets denote an averaging over time as well as 
$\theta$.
The local kinetic temperature in a bin, which is a measure of the fluctuation velocities of the particles within that bin, is defined as
\begin{equation} T(r,z)=\frac{\left\langle m \sum_{i=1}^{N} \left({\boldsymbol{c}_i(t)-\boldsymbol{v}(r,z)}\right) \boldsymbol{\cdot} ({\boldsymbol{c}_i(t)-\boldsymbol{v}(r,z)})\delta(r-r_i)\delta(z-z_i) \right\rangle_\theta}{3k_B \left\langle\sum_{i=1}^{N}{\delta(r-r_i)\delta(z-z_i)}\right\rangle_\theta} . \end{equation}
For the simulations considered here, 
$k_B=m=1$ and 
${\boldsymbol C}=(\boldsymbol {c}-\boldsymbol {v})$ denotes the peculiar velocity of the particles within the bin during the time frame of interest. Using the gap width (
$\delta =R_o-R_i$) and the inverse of the inner-cylinder rotation rate 
$(\omega _i^{-1})$ as the length and time scales, respectively, the following dimensionless fields are defined:
\begin{equation} z^* = \frac{z}{\delta}, \quad t^* = {\omega_i}t, \quad {\boldsymbol v}^* = \frac{\boldsymbol v}{ \omega_i\delta }, \quad T^* = \frac{T}{(\omega_i\delta)^2}. \end{equation}To ascertain whether the system reached a statistical steady state, we monitor the global temperature defined as
\begin{equation} T_s=\frac{\sum^{r}\sum^{z} T(r,z)}{n_r n_z }, \end{equation}
where 
$n_r$ and 
$n_z$ denote the number of bins along the radial and axial directions, respectively. The constancy of the global temperature 
$T_s$ with time is taken as the statistical steady state of the flow – see the main panels in figure 2. The inset of each panel in figure 2 displays a pair of counter-rotating rolls at the steady state for parameter values of 
$(a)$
$\gamma _n=5$ and 
$(b)$
$\gamma _n=0$ (i.e. with zero microscopic damping force), with 
$\omega _i=0.1$ and 
$\varGamma =h/\delta =2$. For both cases, there is an ‘outward’ jet at the mid-height of the cylinder and two ‘inward’ jets near the two endwalls – the resulting pattern is called TVF or Taylor rolls (Taylor Reference Taylor1923). All results presented below correspond to a damping coefficient 
$\gamma _n=5$ (Rapaport Reference Rapaport1998, Reference Rapaport2004) but the robustness of our reported results is verified in appendix A by performing simulations with 
$\gamma _n=0$.
Figure 2. Main panel: temporal evolution of the global kinetic temperature (dimensionless), equation (2.10), of the fluid. Inset: snapshot of steady Taylor vortices, displaying radial–axial velocity vectors. Parameter values are 
$(a)$
$\gamma _n=5$ (microscopic damping coefficient) and 
$(b)$
$\gamma _n=0$, with 
$\varGamma =2$ and 
$\omega _i=0.1$.
3. Onset of Taylor vortices, and the phase diagram of patterns
Referring to the velocity-vector plot in the inset of figure 2
$(a)$, the profiles of the radial and axial velocities can be extracted as functions of the rotation rate. The axial variations of the scaled radial velocity (
$v_r^*={v_r}/{\omega _i\delta }$) are shown in figure 3
$(a)$ for different 
$\omega _i$; the aspect ratio is set to 
$\varGamma =h/\delta =2$. Each profile represents ‘radially averaged’ data over all radial bins, which have also been time averaged (over 100 000 time steps after the production run). The flow is found to be almost purely azimuthal at 
$\omega _i\leq 0.05$, with no motion on the meridional plane. It is seen in figure 3
$(a)$ that, on increasing the angular velocity of the inner cylinder to 
$\omega _i =0.055$, two radial jets (inward flow from the outer towards the inner cylinder) are formed near the top (
$z^*=2$) and bottom (
$z^*=0$) walls and a radial outward flow at the mid-height (
$z^*=1$) of the TC cell is observed. On increasing 
$\omega _i\leq 0.25$, it is seen that the strength of the mid-height radial outward flow increases. As 
$\omega _i$ is increased further to 
$\omega _i \geq 0.3$, the radial velocity 
$v_r(z)$ exhibits two local maxima, indicating two regions of radial flow reversal, which implies the presence of three stationary Taylor vortices; see the velocity profile marked by triangles in figure 3
$(a)$ that corresponds to 
$\omega _i=0.4$. Clearly, such odd-numbered Taylor vortices break the 
$Z_2$ symmetry – this signals the outcome of a secondary bifurcation from the TVF state and will be analysed in detail in § 4.
Figure 3. 
$(a)$ Variations of scaled radial velocity (
$v_r^*=v_r/\omega _i \delta$) with axial distance for different 
$\omega _i$. 
$(b)$ Bifurcation diagram: variation of 
${\rm \Delta} v_r$, (3.1), with 
$\omega _i$; the symbols represent simulation data and the line is a least-squares fit of the form (3.2). The aspect ratio of the TC cell is 
$\varGamma =2$.
3.1. Determination of the critical rotation rate and Reynolds number
To identify the location of the primary bifurcation point 
$\omega _i=\omega ^{cr}$ for the onset of TVF, the simulations have been performed for a range of 
$\omega _i$ for specified values of the aspect ratio 
$\varGamma$. Based on the radial velocity variations such as in figure 3
$(a)$, the maximum radial velocity,
\begin{equation} {\rm \Delta} v_r =\max _{0\leq z^* \leq \Gamma} [v_r (z)], \end{equation}
is chosen as an order parameter to identify the transition from CCF to TVF. The variation of 
${\rm \Delta} v_r$ with 
$\omega _i$ is displayed in figure 3
$(b)$. The critical rotation rate 
$\omega _i^{cr}$ for transition from CCF to TVF is obtained by fitting all data for 
${\rm \Delta} v_r$ using the following function:
\begin{equation} {\rm \Delta} v_r = \alpha(\omega_i - \omega^{cr})^{\beta} , \end{equation}
with 
$\alpha \approx 2.05$ and 
$\beta \approx 0.565$. Estimating the bifurcation point (figure 3
$b$) at which 
${\rm \Delta} v_r\approx 0$ yields a value of 
$\omega ^{cr}\approx 0.054$ for the present system with 
$\varGamma =2$ and 
$\omega _o=0$; this critical rotation rate holds for all 
$\varGamma >1$ (see § 3.2).
Based on the inner-cylinder speed 
$v_{i}=\omega _{i} R_i$ as the velocity scale and the gap width 
$\delta =R_o-R_i$ as the length scale, the Reynolds number can be estimated from
\begin{equation} Re = \frac{\varrho v_{i} \delta}{\mu} , \end{equation}
where 
$\varrho =\rho _p\phi =({\rm \pi} \rho /6)\rho _p$ is the mass density of the fluid, 
$\rho _p$ is the material density of particles and 
$\mu$ is the shear viscosity. The expression of the latter for a ‘dense’ hard-sphere fluid is given by Chapman & Cowling (Reference Chapman and Cowling1970) and Saha & Alam (Reference Saha and Alam2016) as
\begin{equation} \frac{\mu}{\varrho \sigma\sqrt{(k_B/m)T}} = \frac{5\sqrt{\rm \pi}}{96\phi g_0}\left(1 + \frac{8}{5} \phi g_0\right)^2 + \frac{8}{5\sqrt{\rm \pi}}\phi g_0 =f_2(\phi), \end{equation}
where 
$g_0(\phi ) =(1-\phi /2)/(1-\phi )^3$ is the radial distribution function and 
$f_2(\phi )$ is a dimensionless function of particle volume fraction (
$\phi$). Combining (3.3) and (3.4) we arrive at the following expression for the Reynolds number:
\begin{equation} Re = \frac{\varrho \omega_{i}R_i \delta}{\varrho \sigma \sqrt{(k_B/m)T} \, f_2(\phi) } \equiv \frac{(R_i/\sigma)}{f_2(\phi)\sqrt{T^*}}, \end{equation}
with 
$T^*=T/(\omega _i\delta )$ being the dimensionless temperature, and we have set 
$k_B=1=m$ as in the present simulations. The mean/global temperature at the onset of Taylor vortices is 
$T^*=T_{av}^*\approx 0.28$ (i.e. at 
$\omega _i=\omega ^{cr}\approx 0.054$). Substituting the value of 
$T^*$, along with 
$f_2(\phi _{av})\approx 1.1029$ (evaluated at a mean volume fraction of 
$\phi =\phi _{av}={\rm \pi} \rho _{av}/6\approx 0.2618$) and 
$R_i/\sigma =50$ in (3.5), we find a critical Reynolds number of 
$Re^{cr}\approx 85.6$.
On the other hand, for the same system without microscopic damping (
$\gamma _n=0$, see appendix A), we have 
$\omega ^{cr}=0.07$ and 
$T^*=T_{av}^*\approx 0.4028$, resulting in a critical Reynolds number of 
$Re^{cr}(\gamma _n=0)\approx 71.4 < Re^{cr}(\gamma _n\neq 0)$. The above estimates do not depend on the aspect ratio (
$\varGamma =h/\delta$) of the TC cell as long as 
$\varGamma > 1$. Overall, the microscopic dissipation (
$\gamma _n>0$) increases the value of the critical Reynolds number, thereby delaying the onset of Taylor vortices, which is expected. Note that our estimate of 
$Re^{cr}(\gamma _n=0)\approx 71.4$ is close to the predicted value of 
$Re^{cr}=76.4$ (Sparrow, Munro & Jonsson Reference Sparrow, Munro and Jonsson1964) for an incompressible TCF with a radius ratio of 
$\eta =2/3$.
3.2. Phase diagram of patterns: symmetric and asymmetric vortices
Following the procedure described in § 3.1, the critical rotation rate 
$\omega ^{cr}$ for the bifurcation from CCF to TVF has been determined for different values of the aspect ratio 
$\varGamma$. The neutral/critical locus of 
$\omega _i=\omega ^{cr}$ is marked in figure 4 by the dashed and dot-dashed lines. It is seen that the critical rotation rate is nearly constant (
$\omega ^{cr}\approx 0.054$) for 
$\varGamma \geq 1$, and the data points for 
$\varGamma \in (0.2, 1)$ can be fitted via
\begin{equation} \omega^{cr} = \frac{0.051}{\varGamma}, \quad \forall \ \varGamma <1, \end{equation}
marked by the dashed line in figure 4. Equation (3.6) implies that the critical rotation rate increases for short cylinders or, in other words, the transition from CCF to TVF is delayed in short cylinders with 
$\varGamma < 1$.
Figure 4. Phase diagram of patterns in the (
$\varGamma , \omega _i$) plane. Different symbols denote different numbers of Taylor vortices as shown in the legend (bottom right corner); the empty symbols denote a sense of rotation opposite to the one observed in the case of the corresponding filled symbols. Each symbol at (
$\varGamma , \omega _i$) corresponds to a steady-state pattern that has been obtained from an individual run (such as figure 2) by uniformly accelerating the inner cylinder to a specified 
$\omega _i$ in 
$10^5$ time steps (
${\rm \Delta} t$) from the rest state, and subsequently allowing the simulations to run for 
$2\times 10^6{\rm \Delta} t$ to calculate hydrodynamic fields; see § 2.2 for details on simulation protocols. The vector visualizations corresponding to points A, B, D and F have been scaled by factors of 
$15$, 
$2$, 
$5$ and 
$2$, respectively, to improve clarity.
The main panel of figure 4 summarizes all patterns for selected aspect ratios 
$\varGamma \in (0.2,4)$ that have been obtained from long-time simulations (such as in figure 2) by varying the inner-cylinder rotation rate 
$\omega _i \in (0.05,0.4)$ that approximately spans a range of inner-cylinder Reynolds number of 
$Re=\omega _iR_i\delta /(\mu /\varrho )\in (50,\ 500)$. The Taylor vortices having different numbers of rolls are marked in figure 4 by different symbols at each 
$\varGamma$. The steady-state snapshots with increasing 
$\omega _i>\omega ^{cr}$, depicting transitions from one pattern to another, are shown on the left and right panels of figure 4 in terms of the radial–axial velocity-vector plots. For example, the snapshots I and J represent a transition from a two-roll to a three-roll TVF (marked by diamonds and triangles, respectively, in figure 4) that occurs at 
$\omega _i\approx 0.27$ at 
$\varGamma =2$. A transition from a four-roll to a five-roll state can be ascertained from the snapshots K and L displayed on the lower right panels that occurs at 
$\varGamma =4$.
The snapshots A 
$\to$ B and C 
$\to$ D represent transition from (i) an ‘unsteady’ state to a ‘steady’ single-roll state and (ii) a single-roll state to an ‘asymmetric’ two-roll state at 
$\varGamma =0.5$ and 
$0.75$, respectively. When the aspect ratio is increased to 
$\varGamma =1$, a transition from a single-roll state to a ‘symmetric’ two-roll state is found to occur; see the snapshots E 
$\to$ F. On the other hand, the snapshots G 
$\to$ H represent a transition between two possible two-roll configurations, with one having ‘outward’ flows (G) near the endwalls and the other having ‘inward’ flows (H) near the endwalls. The latter transition will be analysed in detail in § 5.1.
Among the patterns identified in figure 4, the odd-numbered or asymmetric rolls break the 
$Z_2$-symmetry (reflection about the mid-height 
$z=h/2$). Such asymmetric states, along with the midplane symmetric Taylor vortices with outward flow at the endwalls, are collectively known as ‘anomalous modes’ according to the classification of Benjamin (Reference Benjamin1978a,Reference Benjaminb). The latter-type vortices have been extensively studied in incompressible TCF with no-slip, free-slip and asymmetric boundary conditions at the endwalls (Benjamin & Mullin Reference Benjamin and Mullin1981; Cliffe Reference Cliffe1983; Lücke et al. Reference Lücke, Mihelcic, Wingerath and Pfister1984; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989; Mullin & Blohm Reference Mullin and Blohm2001). In the context of the present findings, pertinent questions are these: Are the present asymmetric roll states similar to those known in incompressible fluids? How does the compressibility of the fluid affect the primary and secondary bifurcation scenario? Does the bifurcation between different types of roll states occur via supercritical or subcritical transitions? These questions are addressed below.
4. Characterization of Taylor vortices and the bifurcation scenario
The hydrodynamic fields are used to characterize the Taylor-vortex patterns (such as those in figure 4) as discussed in § 4.1. Various order parameters are considered in § 4.1.1 to identify the bifurcation scenario among different patterns as a function of the inner-cylinder rotation. The results pertaining to short cylinders (
$\varGamma \leq 1$), leading to the emergence of single-vortex and asymmetric two-vortex states and their interrelations, are discussed separately in § 4.2.
4.1. Results for aspect ratios 
$\varGamma =2$ and $4$: stratification and asymmetric states
For the case of 
$\varGamma =2$, we display the radial–axial velocity field in figure 5
$(a)$ for four representative values of inner-cylinder rotation 
$\omega _i$. A pair of (weak) counter-rotating rolls, with an outward jet at the mid-height of the TC cell, is seen in the leftmost panel, which corresponds to 
$\omega _i=0.055>\omega ^{cr}$. With increasing 
$\omega _i$, the rolls and the outward jet become stronger, as seen in the second (
$\omega _i=0.1$) and third (
$\omega _i=0.25$) panels of figure 5
$(a)$. Increasing rotation beyond 
$\omega _i=0.25$ results in a three-vortex configuration; see the rightmost panel of figure 5
$(a)$, which corresponds to 
$\omega _i=0.4$. The corresponding variations of density and temperature in the meridional plane can be ascertained from figure 5
$(b)$ and 
$(c)$, respectively. In particular, the outward jets correspond to relatively dilute and hotter flows through which the angular momentum is transferred from the rotating inner cylinder to the stationary outer cylinder (see below).
Figure 5. 
$(a)$ Radial–axial velocity vectors, 
$(b)$ reduced density and 
$(c)$ temperature in the (
$r,z$) plane for different 
$\omega _i$. 
$(d)$ Colour maps of the Mach number 
$Ma=v_\theta /c_s$, where 
$c_s$ is the local sound speed. The aspect ratio is 
$\varGamma =2$; reflective boundary conditions are imposed at axial endwalls.
Figure 5
$(d)$ represents the colour maps of the local Mach number based on the azimuthal velocity,
\begin{equation} Ma=\frac{v_{\theta}}{c_s}, \end{equation}
where 
$c_s$ is the local sound speed (Savage Reference Savage1988):
\begin{equation} c_s = \sqrt{\left(\frac{\partial p}{\partial\varrho}\right)_s } \equiv \sqrt {T \left(1+\frac{2}{3}\,f_1(\rho)+ \rho\frac{\textrm{d} \ln f_1(\rho)}{\textrm{d}\rho} \right)f_1(\rho)}. \end{equation}
The latter expression follows from the equation of state for a dense gas, 
$p= \varrho f_1(\rho )T$, where 
$\varrho =\rho _p({\rm \pi} /6)\rho$ is the mass density, 
$\rho _p$ is the material density of particles, 
$f_1(\rho )=1+ (2{\rm \pi} /3) \rho g(\rho )$ is a dimensionless factor incorporating excluded-volume effects and 
$g(\rho )$ is the radial distribution function,
\begin{equation} g(\rho)=\frac{\rho_{max}^2(\rho_{max}- \rho/2)}{(\rho_{max}-\rho)^3}, \end{equation}
with 
$\rho _{max}\leq 6/{\rm \pi}$. For a dilute gas (
$\rho \to 0$), the dense-gas factor is 
$f_1(\rho )\to 1$ and hence we recover the well-known expression for the sound speed,
\begin{equation} c_s = \sqrt{\gamma T} = \sqrt {5T/3}, \end{equation}
with 
$\gamma =5/3$ being the ratio of specific heats for a monatomic gas. It is clear from the panels in figure 5
$(d)$ that there are considerable variations of Mach number in the (
$r,\!z$) plane, with its maximum being located near the rotating inner cylinder and across the outward jet. In fact, the azimuthal flow velocity can exceed the sound speed, and the flow being subsonic seems to hold only for rotation rates below the onset of Taylor vortices.
Figure 6 displays the time-averaged profiles of (
$a{,}d$) the reduced density, (
$b{,}e$) the kinetic temperature and (
$c{,}f$) the azimuthal velocity. Note that the radial profiles of the hydrodynamic fields represent spatial averaging over both azimuthal (
$\theta$) and axial (
$z$) directions, i.e.
\begin{equation} v_\theta(r) = \langle v_\theta(r,z,\theta) \rangle_{\theta z}, \quad \rho(r) = \langle \rho(r,z,\theta) \rangle_{\theta z} \quad \text{and} \quad T(r) = \langle T(r,z,\theta) \rangle_{\theta z}, \end{equation}while their axial profiles in the lower panel of figure 6,
\begin{equation} v_\theta(z) = \langle v_\theta(r,z,\theta) \rangle_{\theta r}, \quad \rho(z) = \langle \rho(r,z,\theta) \rangle_{\theta r} \quad \text{and} \quad T(z) = \langle T(r,z,\theta) \rangle_{\theta r}, \end{equation}have been averaged over both azimuthal and radial directions.
Figure 6. Top row: radial variations of the axially averaged 
$(a)$ reduced density 
$\rho (r)=\langle \rho (r,z)\rangle _z$, 
$(b)$ temperature 
$T(r)=\langle T(r,z)\rangle _z$ and 
$(c)$ azimuthal velocity 
$v_{\theta }(r)=\langle v_\theta (r,z)\rangle _z$. Bottom row: axial variations of the radially averaged 
$(d)$ reduced density 
$\rho (z)=\langle \rho (r,z)\rangle _r$, 
$(e)$ temperature 
$T(z)=\langle T(r,z)\rangle _r$ and 
$(\,f)$ azimuthal velocity 
$v_{\theta }(z)=\langle v_\theta (r,z)\rangle _r$. Parameter values are the same as in figure 5.
The distinguishing features of the compressible TCF can be ascertained from figure 6(
$a{,}b$) and (
$d{,}e$). There are considerable density and temperature variations along both the radial and axial directions; clearly, the fluid is compressible and non-isothermal. It is seen in figure 6
$(a)$ that the peak of 
$\rho (r)$ shifts towards the outer cylinder with increasing rotation rate of the inner cylinder, implying that the material is thrown outwards due to increasing centrifugal force. The competition between the centrifugal force and the radial pressure gradient would decide the location of the density maximum. The radial variation of the temperature profile in figure 6
$(b)$ has a minimum within the annular gap – the origin of this temperature minimum can be explained by three competing factors: (i) the thermal diffusion, (ii) the added shear work due to the rotating inner cylinder and (iii) a cooling mechanism due to the microscopic damping force (
$\gamma _n>0$); see (A 1) in appendix A. The discussion on temperature variations in the meridional plane and a comparison with the undamped case (
$\gamma _n=0$) are deferred to appendix A (see figures 23
$b$ and 24
$b$).
The axial density profiles in figure 6
$(d)$ indicate the emergence of axial density stratification, i.e. alternate regions of dilute and dense layers are formed along the axial direction. Such axial stratification of density was also observed in the context of ‘granular’ TCF where the particle collisions are inelastic (Mahajan Reference Mahajan2016; Mahajan & Alam Reference Mahajan and Alam2016). Comparing figure 6
$(d)$ with corresponding radial velocity profiles in figure 3
$(a)$, we find that the vortex centres are identified with the average density (
$\rho _{av}=0.5$), but the outward and inward jets correspond to minimum and maximum density, respectively. Similar inferences can be made about the temperature field in figure 6
$(e)$.
With increasing inner-cylinder rotation, the azimuthal velocity 
$v_\theta (r)$ profile becomes flatter (i.e. a low shear region) around the mid-gap of the annulus (see figure 6
$c$), while its axial variations 
$v_\theta (z)$ in figure 6
$(\,f)$ display similar geometrical features to its radial velocity profile (figure 3
$a$), with larger and smaller velocities representing outward and inward jets, respectively. The appearance of a plateau in 
$v_\theta (r)$ can be explained by the azimuthal momentum balance of steady, axisymmetric, azimuthal (with 
$\partial /\partial z(\boldsymbol{\cdot} ) =0$ and 
$\partial /\partial \theta (\boldsymbol{\cdot})=0$) flow:
\begin{equation} \frac{\textrm{d}}{\textrm{d}r} (r^2\Sigma_{r\theta} )= 0 , \end{equation}
where 
$\Sigma _{r\theta } = -\mu r ({\textrm {d}}/{\textrm {d}r}) (v_\theta /r)$ is the shear stress. A decreasing shear rate would necessarily imply an increased shear viscosity 
$\mu \sim f_2(\rho )\sqrt {T}$ (see (3.4)),which can happen if 
$f_2(\rho )\sim g(\rho )$ (the radial distribution function, (4.3)) increases strongly with increasing density. We have verified that the decrease of 
$\mu (\rho , T)$ due to decreasing 
$T$ with increasing 
$\omega _i$ (see figure 6
$b$) is compensated by its stronger increase (
$\propto g(\rho )$) due to increased density in the same limit. The latter finding can be ascertained from a comparison of the temperature profiles in figure 6
$(b)$ between 
$\omega _i=0.1$ and 
$0.25$.
The colour maps of the time-averaged specific angular momentum, averaged over the azimuthal direction, defined as
\begin{equation} {\mathcal{L}}(r,z) = \langle \rho(r,z,\theta) v_\theta(r,z,\theta) r \rangle_{\theta}, \end{equation}
are displayed in figure 7
$(a)$, with parameter values as in figure 5. The angular momentum is transported via the outward jets that correspond to larger values of 
${\mathcal {L}}$; the related signatures of two outward jets are seen in the last panel of figure 7
$(a)$, which represents a three-roll state. With increasing rotation rate of the inner cylinder, the radial variation of 
${\mathcal {L}}(r,z)$ along the outward jet becomes non-monotonic; see the second and third panels of figure 7
$(a)$. The above non-monotonicity is also evident in the radial profiles of axially averaged 
${\mathcal {L}}(r) = \langle {\mathcal {L}}(r,z) \rangle _{z}$; see figure 7
$(b)$. It is clear that 
${\mathcal {L}}(r)$ decreases monotonically towards the stationary outer cylinder for 
$\omega _i=0.055$ and 
$0.1$ (representing the centrifugally unstable (Rayleigh Reference Rayleigh1917) or ‘Rayleigh-unstable’ regime), but, at higher values of 
$\omega _i=0.25$ and 
$0.4$, 
${\mathcal {L}}(r)$ increases up to some radial distance beyond which it decreases towards the outer cylinder. The latter profiles can be divided into two parts: (i) a ‘Rayleigh-stable’ part located near the inner cylinder and (ii) a ‘Rayleigh-unstable’ part near the outer cylinder. Figure 7
$(c)$ shows the axial profiles of the radially averaged angular momentum 
${\mathcal {L}}(z) = \langle {\mathcal {L}}(r,z) \rangle _{r}$ – the structural features of 
${\mathcal {L}}(z)$ with increasing 
$\omega _i$ resemble that of the azimuthal velocity field (figure 6
$f$).
Figure 7. 
$(a)$ Colour maps of the specific angular momentum 
${\mathcal {L}}(r,z)=\langle \rho v_\theta r \rangle _\theta$ in the (
$r, z)$ plane for different 
$\omega _i$. Bottom row: variations of specific angular momentum 
$(b)$ axially averaged 
$\langle {\mathcal {L}}(r,z)\rangle _z (r)$ and 
$(c)$ radially averaged 
$\langle {\mathcal {L}}(r,z)\rangle _r(z)$. Parameter values are the same as in figure 5.
For a longer TC cell with an aspect ratio of 
$\varGamma =4$, the time-averaged radial–axial velocity fields in the (
$r, z$) plane are shown in figure 8
$(a)$ for different values of inner-cylinder rotation 
$\omega _i$. The colour maps of the azimuthal velocity 
$v_{\theta }(r,z)$ and the specific angular momentum field 
${\mathcal {L}}(r,z)$ are displayed in figures 8
$(b)$ and 8
$(c)$, respectively. The image in the leftmost panel of figure 8
$(a)$ indicates that the primary bifurcation from the CCF state leads to a weak four-roll state at 
$\omega _i=0.055$. The vortices become stronger with increasing rotation rate, as seen in the second panel of figure 8
$(a)$ at 
$\omega _i=0.1$, for which two outward jets are also evident in the angular momentum field in the second panel of figure 8
$(c)$. For both 
$\omega _i=0.055$ and 
$0.1$, the flow is radially inwards (towards the inner cylinder) near endwalls and they represent normal Taylor vortices.
Figure 8. Evolution of steady Taylor rolls with increasing rotation rate of the inner cylinder for 
$\varGamma =4$: 
$(a)$ axial–radial velocity vectors, 
$(b)$ azimuthal velocity and 
$(c)$ specific angular momentum field. While the four-roll state at 
$\omega _i=0.1$ represents Taylor vortices with inward flows at endwalls, the four-roll (
$\omega _i=0.15$), five-roll (
$\omega _i=0.25$) and six-roll (
$\omega _i=0.4$) states have outward flows near endwalls.
At a slightly higher rotation of 
$\omega _i=0.15$ (see the third panel in each of figure 8
$a$–
$c$), however, we find a four-roll state with outward jets near the endwalls – we shall show in § 5.1 that the present system with reflective endwall conditions admits ‘perfect’ pitchfork bifurcation (figures 10, 11 and 14) and hence the even-numbered roll state with either inward or outward jets near the endwalls are equivalent solutions. Note in figure 8(
$b{,}c$) that the outward jets are characterized by maxima in 
$v_\theta$ and 
${\mathcal {L}}$. This four-roll mode transitions to a five-roll state at 
$\omega _i=0.25$; see the fourth panels of figure 8(
$a$–
$c$). With increasing 
$\omega _i$, a four-roll state, with inward jets near endwalls, reappears at 
$\omega _i=0.3$ (not shown) and a five-roll state at 
$\omega _i=0.35$ (not shown) – these states have been marked in the phase diagram in figure 4. With further increasing rotation rate to 
$\omega _i=0.4$, a six-roll state (with outward jets near endwalls) is found whose meridional-plane structure can be ascertained from the last panels of figure 8(
$a$–
$c$).
Figure 9. Top row: radial variations of 
$(a)$ the reduced density 
$\rho (r)=\langle \rho (r,z)\rangle _z$, 
$(b)$ the azimuthal velocity 
$v_{\theta }(r)=\langle v_\theta (r,z)\rangle _z$ and 
$(c)$ the specific angular momentum 
${\mathcal {L}}(r)=\langle {\mathcal {L}}(r,z)\rangle _{z}$. Bottom row: axial variations of 
$(d)$ the reduced density 
$\rho (z)=\langle \rho (r,z)\rangle _r$, 
$(e)$ the azimuthal velocity 
$v_{\theta }(z)=\langle v_\theta (r,z)\rangle _r$ and 
$(\,f)$ the specific angular momentum 
${\mathcal {L}}(z)=\langle {\mathcal {L}}(r,z)\rangle _{r}$. The aspect ratio is 
$\varGamma =4$.
Figure 10. Bifurcation diagrams in terms of 
$(a)$ (
${\rm \Delta} v_r, \omega _i$), 
$(b)$ (
$\zeta , \omega _i$) and 
$(c)$
$({\rm \Delta} \rho , \omega _i$). The inset in each panel represents a magnified sketch showing the hysteretic region. The aspect ratio is 
$\varGamma =2$.
Figure 11. Bifurcation diagrams in terms of 
$(a)$ (
${\rm \Delta} v_r, \omega _i$) and 
$(b)$ (
$\zeta , \omega _i$). The inset zoom in panel 
$(a)$ shows a sketch of the possible hysteresis loop in the (
${\rm \Delta} v_r, \omega _i$) plane. The aspect ratio is 
$\varGamma =4$.
Figure 9 displays the radial and axial profiles of the reduced density, the azimuthal velocity and the specific angular momentum field in the upper and lower panels, respectively, for 
$\varGamma =4$. The overall radial and axial variations of these hydrodynamic fields look similar to those for 
$\varGamma =2$ (figures 6 and 7
$b$,
$c$). Note that the axial density stratification becomes stronger with increasing 
$\omega _i$ (figure 9
$d$), with the density maxima and minima coinciding across the inward and outward jets (figure 9
$e$), respectively, which correspond to channels of the (inward/outward) radial transport of angular momentum (figure 9
$f$). The latter reconfirms our previous finding that the angular momentum is transported from the inner cylinder to the outer cylinder via the relatively rarefied (dilute) outward jets.
One key finding of this section is that the compressible Taylor vortices are embedded with stratifications in density, temperature, viscosity and specific angular momentum fields along both the axial and radial directions that become stronger with increasing rotational speed of the inner cylinder. In particular, the non-monotonic radial stratification of density (figures 6
$a$ and 9
$a$), with a density maximum around the mid-gap, is a key structural signature of compressible Taylor vortices in a dense gas. This is in contrast to the TCF of a dilute gas (Welsh et al. Reference Welsh, Kersal'e and Jones2014; Aghor & Alam Reference Aghor and Alam2020), where the density maximum occurs at the outer cylinder when the Prandtl number is of order one. Related issues on the temperature field are discussed in appendix A. The other key finding is a secondary transition from a symmetric even-roll state to an asymmetric odd-roll state beyond some critical value of 
$\omega _i\gg \omega ^{cr}$. The signature of this secondary bifurcation is implicated in the non-monotonic radial variation of angular momentum (figures 7
$b$ and 9
$c$), with a Rayleigh-stable (centrifugally stable) part being located near the inner cylinder and a Rayleigh-unstable part near the outer cylinder. The nature of the latter bifurcation is characterized in terms of different order parameters in § 4.1.1.
4.1.1. Secondary bifurcation and the hysteretic transition to asymmetric vortices
Figure 10 displays bifurcation diagrams in terms of different order parameters for the case of 
$\varGamma =2$. Figure 10(
$a$) shows the variation of the maximum radial velocity 
${\rm \Delta} v_r$, (3.1), with 
$\omega _i$, whereas figure 10(
$b$) shows the variation of the net circulation as a function of 
$\omega _i$. Here, the net circulation of vortices is calculated from
\begin{equation} \zeta= \frac{1}{n_v n_l} \sum_{i=1}^{n_v}\sum_{j=1}^{n_l}{\boldsymbol{v}(r,z)\boldsymbol{\cdot}\,\textrm{d}l}, \end{equation}
where 
$n_v$ is the number of vortices and 
$n_l$ is the number of loops around each vortex centre. To evaluate (4.9), the vortex centres are first identified (i.e. the point at which 
$(v_r, v_z)\approx 0$) and two square loops of sides 
${\rm \Delta} l=0.32\delta$ and 
$0.42\delta$ around each vortex centre are defined. The net circulation 
$\zeta$, (4.9), is then calculated over each of these loops, resulting in 
$n_l=2$. This procedure is repeated over all vortices in the system. If the system has an even number of symmetric Taylor vortices, such as in figure 2, the net circulation is zero (
$\zeta =0$). Clearly, the zero and non-zero values of 
$\zeta$, (4.9), demarcate the transition from a symmetric two-roll state to a three-roll state that occurs at 
$\omega _i\approx 0.26$ (see figure 10
$b$).
Figure 10(
$a$) indicates that, while the maximum radial velocity 
${\rm \Delta} v_r$ increases with increasing 
$\omega _i$ in the two-roll state (implying that the strength of the inward/outward jets increases), there is a discontinuous drop in 
${\rm \Delta} v_r$ when the number of vortices increases from two to three. The latter is expected since the angular momentum of two vortices must be redistributed to three vortices at the ‘
$2\leftrightarrow 3$’ roll transition. The above transition is further clarified in the inset (schematic) of figure 10
$(a)$, which indicates that the ‘
$2\leftrightarrow 3$’ roll transition appears to be hysteretic over a small range of inner-cylinder rotational rates about 
$\omega _i\approx 0.26$. Note that to exactly locate the limit point and thereby confirm the hysteretic nature of the ‘
$2\leftrightarrow 3$’ roll transition would require a ‘numerical continuation’ in which the rotational speed of the inner cylinder is quasi-statically increased from rest over a large range of 
$\omega _i$ successively, which has not been checked for this case (see figure 21 and the text in appendix A on related issues).
The bifurcation diagram in figure 10
$(a)$ is replotted as figure 10
$(c)$ in terms of the following order parameter:
\begin{equation} {\rm \Delta}\rho = \max_{z} \rho(z) - \min_{z} \rho(z), \end{equation}
called the axial density contrast; refer to figure 6
$(d)$ for related axial variations of density for 
$\varGamma =2$. It is clear that 
${\rm \Delta} \rho$ increases with increasing 
$\omega _i$ in the Taylor-vortex regime, and its behaviour with 
$\omega _i$ mirrors that of 
${\rm \Delta} v_r$ in figure 10
$(a)$. The inset (schematic) of figure 10
$(c)$ reconfirms that the ‘
$2\leftrightarrow 3$’ roll transition is indeed hysteretic, with a mild increase in 
${\rm \Delta} \rho$ at 
$\omega _i\approx 0.263$, which implies that the axial density contrast increases in the three-roll state compared to its value in the two-roll state.
The analogue of figure 10(
$a$,
$b$) for the case of 
$\varGamma =4$ is displayed in figure 11(
$a$,
$b$). Here the primary bifurcation from CCF results in a four-roll state at 
$\omega _i\approx 0.054$ that persists over a range of 
$\omega _i\leq 0.3$. The inset of figure 11
$(a)$ represents its zoomed-in version, with the left and right dashed lines indicating possible hysteretic transitions, respectively, between the ‘
$4\leftrightarrow 5$’ and ‘
$5\leftrightarrow 6$’ roll states. The hysteretic nature of secondary bifurcations (the ‘
$4\leftrightarrow 5$’ and ‘
$5\leftrightarrow 6$’ rolls) is made clear in figure 11
$(b)$, which shows the variation of the net circulation 
$\zeta$, (4.9), with 
$\omega _i$. The non-zero circulation branches in figure 11
$(b)$ refer to the five-roll state that persists over the range of 
$0.18 \leq \omega _i\leq 0.37$; on the other hand, the four-roll and six-roll states persist over 
$0.054\leq \omega _i \leq 0.3$ and 
$\omega _i\geq 0.34$, respectively. It is interesting to note that the range of 
$\omega _i$ over which the system is bistable (the ‘
$4+5$’ roll and ‘
$5+6$’ roll states) is much larger at 
$\varGamma =4$ compared to the bistable regime for 
$\varGamma =2$ (see figure 10). The above hysteretic transition scenario needs to be verified via numerical continuation as discussed in appendix A.
It should be noted that the bifurcation diagrams depicted in different panels of figures 10 and 11 are, in fact, appropriate projections of a multidimensional hypersurface (
${\boldsymbol A}, \omega _i$), with the order parameter 
${\boldsymbol A}=(\zeta , {\rm \Delta} v_r, {\rm \Delta} \rho , {\rm \Delta} T, \ldots )$ being a vector-valued function (Golubitsky & Schaeffer Reference Golubitsky and Schaeffer1985). This can lead to different solution branches nearly touching each other (such as in figure 11
$a$) in a two-dimensional projection of the full bifurcation diagram. This has been further verified by carrying out simulations in a much larger TC cell with 
$\varGamma =10$ (not shown) – we found that the bifurcation sequence among different-numbered (
$10$, 
$11$, 
$12$, 
$13$ and 
$14$ rolls) states for 
$\varGamma =10$ appears to be smooth when 
${\rm \Delta} v_r$ is used as an order parameter, although the odd-numbered rolls are clearly picked up when the net circulation 
$\zeta$ is used as an order parameter.
4.2. Single-vortex and asymmetric two-vortex states: $\varGamma \leq 1$
Here we discuss results pertaining to short cylinders with aspect ratios of 
$\varGamma \leq O(1)$; some of the observed patterns are marked in the phase diagram in figure 4. The bifurcation diagram for 
$\varGamma =1$ is shown in the middle panel of figure 12 in the (
${\rm \Delta} v_r, \omega _i$) plane. The primary bifurcation from the CCF state yields a ‘single-vortex’ state for which three snapshots are displayed on the left (
$\omega _i=0.065$) and top (
$\omega _i=0.15$ and 
$0.3$) panels of figure 12. We use the terminology ‘single vortex’ to refer to the existence of one vortex spanning the axial height of the TC cell. Each snapshot represents streamline patterns (at the steady state) in the (
$r, z$) plane. This has been drawn by employing the ‘surface-LIC’ (surface line integral convolution) module of the ParaView software (Ahrens, Geveci & Law Reference Ahrens, Geveci and Law2005). Note that the single-vortex state is slightly off-centred, indicating its geometrical asymmetry, which persists over the range of 
$\omega _i$ studied. The bottom row and the right column of figure 12 display streamline patterns that belong to the lower branch of the bifurcation diagram. An asymmetric pair of vortices is found at 
$\omega _i=0.085$ (the bottom left panel) that persists up to a rotation rate of 
$\omega _i=0.1$. Both these asymmetric two-roll states consist of a stronger vortex lying below a weaker one. With increasing rotation rate, the asymmetric two-vortex state degenerates into a symmetric two-roll state; see the snapshots for 
$\omega _i=0.25$ (the bottom right panel), 
$\omega _i=0.3$ (middle right panel) and 
$\omega _i=0.4$ (top right panel) in figure 12.
Figure 12. Evolution of Taylor rolls with increasing 
$\omega _i$ along the upper and lower branches of the bifurcation diagram (middle panel) for 
$\varGamma =1$. The middle panel represents the bifurcation diagram in the (
${\rm \Delta} v_r, \omega _i$) plane, with the circles and diamonds representing single-roll and two-roll states, respectively, obtained from simulations in which the inner cylinder is started from the rest state to a specified 
$\omega _i$ as described in § 2.2. The triangles on the upper and lower branches were obtained from two separate ‘continuation’ runs in which the inner-cylinder rotation 
$\omega _i$ is (i) increased from 
$\omega _i=0.25$ (upper branch) to 
$0.35$ and (ii) decreased from 
$\omega _i=0.17$ (lower branch) to 
$0.06$ in steps of 
${\rm \Delta} \omega _i=0.01$. The ramping rate is 
$|\textrm {d}\omega _i/\textrm {d}t|= 2\times 10^{-5}$ for both cases. Refer to figure 21 for the protocol on ‘
$\omega _i$ continuation’ runs.
The above sequence of pattern evolution can be understood from the bifurcation diagram (viz. the middle panel in figure 12) on which the upper and lower branch solutions are connected via two limit points through an unstable branch (marked by the dashed line). It is clear that the underlying bifurcation is mediated via a hysteretic transition, resulting in (i) ‘1-roll 
$\rightarrow$ 2-rolls’ transition with increasing 
$\omega _i$ from the CCF, with the transition to the two-roll state occurring at the right limit point 
$\omega _i\approx 0.3$, and (ii) ‘symmetric 2-rolls 
$\to$ asymmetric 2-rolls 
$\to$ 1-roll’ transition with decreasing 
$\omega _i$ from the symmetric two-roll state, with the transition to the single-roll state occurring at the left limit point 
$\omega _i\approx 0.08$. We have verified that the above bifurcation sequence also holds at smaller values of 
$\varGamma =0.75$ for which the bifurcation diagram looks similar to that in figure 12. These findings on the overall bifurcation scenario for the single-vortex mode and the impact of stationary endwalls on them will be compared and contrasted with related works on incompressible TCF (Benjamin & Mullin Reference Benjamin and Mullin1981; Furukawa et al. Reference Furukawa, Watanabe, Toya and Nakamura2002) in § 5.2.
The hydrodynamic fields for 
$\varGamma =1$ are displayed in figure 13. Note that the radial profiles of the azimuthal velocity 
$v_\theta (r)$ look structurally similar (with the formation of a plateau around the mid-gap at 
$\omega _i\geq 0.15$) to those for 
$\varGamma =2$ (figure 6
$c$) and 
$\varGamma =4$ (figure 9
$b$). The asymmetric 
$v_\theta (z)$ profiles for 
$\omega _i=0.065$ and 
$0.15$ in figure 13(d) represent the single-vortex mode (viz. figure 12). The remaining profiles for 
$\omega _i=0.1$, 
$0.3$ and 
$0.4$ represent two-roll modes and confirm the existence of a mid-height ‘inward’ jet and two ‘outward’ jets near two walls in each case. The key signature of compressible Taylor vortices is strong density stratification along both radial and axial directions that persists also in short cylinders with 
$\varGamma \leq O(1)$ as is evident from figure 13
$(b{,}e)$. The specific angular momentum field 
${\mathcal {L}}(r)$ in figure 13
$(c)$ indicates the emergence of its non-monotonic radial variation with increasing 
$\omega _i\geq 0.25$, similar to those uncovered for 
$\varGamma =2$ (figure 7
$b$) and 
$\varGamma =4$ (figure 9
$c$).
Figure 13. Radial (top row) and axial (bottom row) variations of 
$(a{,}d)$ azimuthal velocity, 
$(b{,}e)$ reduced density and 
$(c{,}f)$ specific angular momentum for 
$\varGamma =1$. The corresponding snapshots can be found in figure 12, which refer to single-roll (
$\omega _i=0.065, 0.15$), asymmetric two-roll (
$\omega _i=0.1$) and symmetric two-roll (
$\omega _i=0.25, 0.4$) states.
About the last issue, the origin of the non-monotonic angular momentum profiles (viz. figures 7
$b$, 9
$c$ and 13
$c$) can be tied to the underlying density stratification, as we explain below. Recall that we have found strong radial variation of density (see figures 6
$a$, 9
$a$ and 13
$b$), with the density maximum shifting towards the stationary outer cylinder with increasing 
$\omega _i$. As explained in § 4.1, the constancy of the shear stress across the Couette gap immediately leads to the formation of a plateau in the azimuthal velocity (figures 6
$c$, 9
$c$ and 13
$a$) around the mid-gap with increasing 
$\omega _i$. Therefore, the coupling of the density and azimuthal velocity via the shear stress is likely to be responsible for the non-monotonic radial variation of the angular momentum (
${\mathcal {L}}(r)=\langle \rho v_\theta r \rangle$) in the compressible TCF of a dense gas. Collectively, the signature of the fluid compressibility is implicated as stratifications in all hydrodynamic fields, and the structural characteristics of the hydrodynamic profiles with increasing 
$\omega _i$ look similar at all 
$\varGamma$ (irrespective of the number of vortices along the axial direction).
5. Discussion: role of axial boundary conditions
5.1. Stationary endwalls and the imperfect bifurcation
So far all the presented results have been obtained with reflection boundary conditions (2.4) at both endwalls. The role of ‘stationary’ endwalls was checked by implementing the zero-slip velocity by using the following protocol:
\begin{equation} c'_x=0=c'_y \quad \text{and} \quad c'_z=-c_z, \end{equation}
i.e. the 
$x$- and 
$y$-component velocities of the reflecting particles are set to zero. The bifurcation diagrams for reflective and stationary endwalls are compared in figure 14
$(a)$ for the representative case of 
$\varGamma =2$, with other parameters being identical as in figure 3
$(b)$. Figure 14
$(a)$ indicates that the reflection boundary condition admits a perfect pitchfork bifurcation at CCF 
$\to$ TVF, with its canonical amplitude equation being given by
\begin{equation} \frac{\textrm{d} {\mathcal{A}}}{\textrm{d}t} = \varOmega {\mathcal{A}} - {\mathcal{A}}^3 = f({\mathcal{A}}, \varOmega) \equiv - f(-{\mathcal{A}}, \varOmega), \end{equation}
where 
${\mathcal {A}}={\rm \Delta} v_r/\alpha$ is the order parameter (amplitude) and 
$\varOmega =\omega _i-\omega ^{cr}$ is the control parameter (the distance from the bifurcation point). Both 
${\mathcal {A}}={\pm } \sqrt {\varOmega }$ (i.e. 
${\rm \Delta} v_r = {\pm }\alpha (\omega _i-\omega ^{cr})^{1/2}$) are stable solutions of (5.2) at 
$\omega >\omega ^{cr}$, indicating the symmetric/perfect nature of pitchfork bifurcation (see figure 14
$b$). In the present context, the two solutions 
${\pm }{\mathcal {A}}$ correspond to Taylor-vortex pairs with inward and outward jets, respectively, adjacent to the two endwalls.
Figure 14. 
$(a)$ Bifurcation diagrams for stationary/no-slip (open squares) and reflection (filled circles) endwalls; the solid line is a least-squares fit, 
${\rm \Delta} v_r =\alpha (\omega _i-\omega ^{cr})^{\beta }$ with 
$\omega ^{cr}=0.054$, 
$\alpha =2.05$ and 
$\beta =0.576$, through the circles. The aspect ratio of the TC cell is 
$\varGamma =2$, with other parameters as in figure 3. 
$(b)$ Schematic of imperfect bifurcation. While the reflection boundary condition (BC) admits perfect/symmetric pitchfork bifurcation at CCF 
$\to$ TVF, the no-slip axial BC makes this bifurcation ‘imperfect’, with normal (upper disconnected branch, red colour) and anomalous (lower disconnected branch, red colour) modes being separated from the bifurcation point for ideal pitchfork; see the text in § 5.1 for details.
In the case of stationary (no-slip) endwalls, the two-roll configuration persists at arbitrarily small values of 
$\omega _i$ (see the square symbols in figure 14
$a$), which is due to the existence of Ekman vortices (Coles Reference Coles1965; Andereck et al. Reference Andereck, Liu and Swinney1986). The centrifugal force is weaker near the no-slip endwalls and hence a shear/boundary layer forms near each endwall – the latter constrains the fluid to flow towards the inner cylinder, leading to Ekman vortices near the two endwalls. The square symbols in figure 14
$(a)$ indicate that the bifurcation point is not identifiable and seems to be shifted to 
$\omega _i\ll \omega ^{cr}$ for stationary endwalls. This is an example of ‘imperfect’ pitchfork bifurcation, which is described by the following normal-form equation (Golubitsky & Schaeffer Reference Golubitsky and Schaeffer1985):
\begin{equation} \frac{\textrm{d} {\mathcal{A}}}{\textrm{d}t} = \varOmega {\mathcal{A}} - {\mathcal{A}}^3 + \alpha_0 + \alpha_1{\mathcal{A}}^2 = f({\mathcal{A}}, \varOmega; \alpha_0, \alpha_1), \end{equation}
where 
$\alpha _0$ and 
$\alpha _1$ are the imperfection parameters that measure the degree of deviation from the ideal (
$\alpha _0=0=\alpha _1$) pitchfork bifurcation. The stationary solutions of (5.3) have been obtained with 
$\alpha _0=5\times 10^4$ and 
$\alpha _1=0$, marked by the red lines in figure 14
$(b)$. It is clear that the no-slip axial boundary condition makes this bifurcation imperfect, with normal (upper disconnected branch, red colour) and anomalous (lower disconnected branch, red colour) modes being separated from the bifurcation point (
$\omega _i=\omega ^{cr}$) for ideal pitchfork. The normal modes are preferred solutions with stationary endwalls at which the tangential velocity is zero and consequently the centrifugal force becomes weaker as one approaches the two endwalls, enabling the formation of inward jets near the endwalls. The lower disconnected branch of figure 14
$(b)$ is called ‘anomalous’ mode (Benjamin Reference Benjamin1978a,Reference Benjaminb) – the corresponding vortex pairs have an inward jet at mid-height and two outward jets near the two endwalls. In the language of bifurcation theory (Schaeffer Reference Schaeffer1980; Golubitsky & Schaeffer Reference Golubitsky and Schaeffer1985), the no-slip boundary condition plays the role of imperfection, thereby destroying the ideal pitchfork bifurcation that holds for reflection axial boundary condition.
That the reflection endwalls admit a pair of symmetric Taylor vortices with both inward and outward jets near the endwalls during its time evolution is demonstrated in figure 15 for the case of 
$\varGamma =1.35$ and 
$\omega _i=0.075$. While the main panel of figure 15
$(a)$ displays the evolution of system temperature over a long time, the snapshots in its inset correspond to patterns at four specific time instants. It is seen that the initial equilibrium state corresponds to a pair of Taylor rolls with an outward jet at the mid-height (see the snapshot at A) which persists up to a dimensionless time of approximately 
$t=4500$, but at later times (
$t> 5000$) this state evolves to a Taylor-vortex pair with an inward jet at the mid-height as seen in the snapshot D. The above transition is mediated via a single-roll state (snapshot B) and an asymmetric two-roll state (C) over a short time interval, and the system temperature also undergoes a sharp change during this change-over.
Figure 15. 
$(a)$ Transition from ‘inward’ to ‘outward’ jets near the endwalls of the TC cell for 
$\varGamma =1.35$ and 
$\omega _i=0.075$. 
$(b)$ Axial variations of the azimuthal velocity 
$v_{\theta }(z)$ averaged over time windows of interest, with parameter values as in panel 
$(a)$. The profiles at A and D are averaged over a time window of 
$100\,000$ time steps, while those at B and C are averaged over 
$20\,000$ time steps. The inset of panel 
$(b)$ displays the radial variations of 
$v_\theta (r)$ at four time instants (A, B, C, D).
Figure 15
$(b)$ and its inset display the radial and axial variations of the azimuthal velocity field at times that correspond to A, B, C and D. While the inset of figure 15
$(b)$ indicates that the azimuthal velocity profiles at A and D are almost indistinguishable, its main panel clarifies that the patterns A and D (see the inset of panel 
$a$) are in fact mirror images of each other and that there is a finite fluid velocity near the endwalls. Clearly, both patterns A and D are valid solutions, representing upper and lower branch solutions arising out of an ideal pitchfork bifurcation CCF 
$\to$ TVF for the present problem with reflection axial boundary conditions; these two possible solutions are marked by the black lines in figure 14
$(b)$. With parameter values as in figure 15, changing the axial boundary condition to stationary endwalls resulted in a vortex pair with an outward jet at the mid-height and two inward jets near the two endwalls; see the velocity-vector plots in the inset of figure 16
$(a)$. This solution belongs to the upper disconnected branch of the bifurcation diagram in figure 14(
$a{,}b$). Note that the axial variation of 
$v_\theta (z)$ in figure 16
$(b)$ confirms the no-slip condition at stationary endwalls for this case.
Figure 16. Same as figure 15, but with no-slip endwalls that yield ‘normal’ Taylor rolls with inward flows near the endwalls.
Based on the above results we conclude that the role of stationary (no-slip) endwalls is to break the perfect pitchfork bifurcation (figure 14) that holds with reflecting endwalls. We have not attempted to track the anomalous branch of the bifurcation diagram (figure 14
$b$) with stationary endwalls. The two sets of axial boundary conditions can be tied via the following protocol:
\begin{equation} c'_T = \epsilon c_T \quad \text{and} \quad c'_z =- c_z, \end{equation}
where 
$c_T=(c_r, c_\theta )$ is the tangential velocity of reflecting particles (
$c_z=-c'_z$), and 
$\epsilon =1$ and 
$0$ refer to reflective and stationary (no-slip) endwalls, respectively. The mixed boundary condition (5.4a,b) can be implemented in future to shed light on the emergence of anomalous modes in compressible TCF.
5.2. Role of endwalls on single-vortex mode and the emergence of twin-cell pattern
To ascertain the effect of stationary endwalls on the single-vortex mode (such as those in figure 12), we show the bifurcation diagram in terms of (
${\rm \Delta} v_r, \omega _i$) in the centre panel of figure 17; the corresponding streamline patterns for different values of 
$\omega _i$ (with corresponding Reynolds number being calculated from 
$Re=\omega _i R_i\delta /(\mu /\varrho )$ by setting the kinematic viscosity to 
$\nu =\mu /\varrho =1$) are displayed in the surrounding panels. The aspect ratio is set to 
$\varGamma =0.75$ and the analogous bifurcation diagram for 
$\varGamma =1$ looks similar (not shown). Note that these simulations were carried out using a larger gap width of 
$\delta =50\sigma$ with 
$R_i=100\sigma$ such that 
$\eta =2/3$ (the radius ratio) and 
$\kappa =1/2$ (the curvature) are fixed – this amounts to an increase in the number of particles, (2.5), by a factor of 
$4$ compared to the case of 
$\delta =25\sigma$ and 
$R_i=50\sigma$. We have verified that the overall topological features of patterns in figure 17 remain unaffected even for a smaller system with gap width of 
$\delta =25\sigma$.
Figure 17. Centre panel: bifurcation diagram with ‘stationary’ (no-slip) endwalls for 
$\varGamma =3/4$, with a gap width of 
$\delta =50\sigma$ and other parameters 
$\eta =2/3$ and 
$\kappa =1/2$. Surrounding panels: surface-LIC plots, depicting streamline patterns, for different values of (
$\omega _i, Re$): (i) (0.01625, 81.25), (ii) (0.02125, 106.25), (iii) (0.03, 150), (iv) (0.0375, 187.5), (v) (0.05, 250), (vi) (0.1, 500), (vii) (0.15, 750) and (viii) (0.25, 1250).
It is clear from figure 17 that the single-vortex mode found at 
$\varGamma \leq 1$ for reflective endwalls (see figure 12 that holds for 
$\varGamma =1$) disappears in the case of stationary endwalls. Instead, each meridional-plane pattern at 
$Re= 81, 106$ and 
$150$ consists of two symmetric large vortices near the inner cylinder coexisting with an array of smaller vortices near the stationary outer cylinder. The pattern near the outer cylinder looks irregular and changes with time – see the three snapshots in figure 18 at different times. The related velocity-vector plots (bottom row of figure 18) confirm that the smaller transient vortices near the outer cylinder are of much weaker strength compared to the two primary vortices near the inner cylinder. With increasing 
$Re$, two large vortices grow laterally in size and the smaller vortices are squeezed towards the outer cylinder – see the snapshot for 
$Re=150$ in figure 17. With further increasing 
$Re$, we find an asymmetric state of one large vortex coexisting with smaller corner vortices (
$Re=188, 250, 500$ and 
$750$). The centre of the large vortex shifts towards the outer cylinder with increasing 
$Re$, a consequence of increased centrifugal force in conjunction with density stratification along the radial direction. At large enough values of 
$Re$, only a two-roll state (with vortex centres being located closer to the outer cylinder) survives with stationary endwalls – see the snapshot for 
$Re=1250$ in figure 17.
Figure 18. Surface-LIC patterns 
$(a\text {--}c)$ and velocity-vector maps 
$(d\text {--}f)$ at different times for the case of 
$Re=106$ with stationary endwalls; other parameters as in figure 17. 
$(a) \ t = 1000; \ (b) \ t = 3000; \ (c) \ t = 5000; \ (d) \ t = 1000; \ (e) \ t = 3000; \ (\,f) \ t = 5000$.
To quantify the time dependence of patterns at smaller 
$Re$ (such as those in figure 18), we show the temporal evolutions of the radial and axial velocities, measured at the mid-height (
$z=h/2$), in figure 19
$(a)$ and 
$(d)$ at two radial locations of 
$r=R_i + \delta /4$ and 
$r=R_i+3\delta /4$, respectively; the Reynolds number is 
$Re=106$. The corresponding phase portraits in terms of (
$v_r, v_z$) in figure 19
$(b)$ and 
$(e)$ suggest the existence of an outward jet (
$v_r>0$) near the inner cylinder (at 
$r=R_i + \delta /4$) and a transient state near the outer cylinder (at 
$r=R_i + 3\delta /4$). The above phase portraits do not show any closed orbit and hence there is no characteristic frequency associated with the coexisting states shown in figure 18. The latter finding is reconfirmed from the power spectral density (PSD) plots displayed in figure 19
$(c)$ and 
$(\,f)$, which are featureless, having no distinct peak.
Figure 19. 
$(a{,}d)$ Time series of radial and axial velocities (
$v_r(t)$ and 
$v_z(t)$), 
$(b{,}e)$ their phase portraits in the (
$v_r, v_z$) plane and 
$(c{,}f)$ power spectral density. All measurements are taken at the mid-height (
$z=h/2$) at radial locations of (
$a$–
$c$) 
$r=\delta /4$ and (
$d$–
$f$) 
$r=3\delta /4$ for 
$Re=106$; other parameters as in figure 18.
By analysing the hydrodynamic fields of figure 18 we have verified that the azimuthal velocity decays faster (compared to the reflective endwall case) close to the rotating inner cylinder and the velocity levels are very small away from the mid-gap location. This is due to the added viscous dissipation near the stationary endwalls, which affects the bulk flow, especially in short cylinders (
$\varGamma \leq 1$). Consequently, the outward angular momentum transport is diminished, resulting in the formation of a twin-cell pattern, which is able to penetrate partially into the Couette gap that coexists with a transient irregular pattern near the outer cylinder. The above scenario holds in short cylinders at smaller values of 
$Re$ (see figure 17), and the angular momentum transport is enhanced with increasing 
$Re$ and the irregular pattern near the outer cylinder eventually disappears at large enough 
$Re$.
In the context of incompressible fluids, Furukawa et al. (Reference Furukawa, Watanabe, Toya and Nakamura2002) reported a twin-cell mode, i.e. two vortices spanning across the Couette gap, in short cylinders (
$\varGamma =2/3$ and 
$\eta =1/2$) at a Reynolds number of approximately 
$Re=300$; they also reported an unsteady mode, consisting of vortex splitting and merging, at a much higher Reynolds number of approximately 
$Re=1500$. These states are qualitatively different from the present finding of the coexistence of two large symmetric vortices near the inner cylinder with an array of transient (smaller) vortices near the outer cylinder that appear at near-critical 
$Re$. On the other hand, the ‘asymmetric’ two-roll states (see the right side panels in figure 17), comprising a primary (larger) vortex and a secondary (smaller) vortex, look structurally similar to those in incompressible TCF. This was first uncovered in experiments by Benjamin & Mullin (Reference Benjamin and Mullin1981) and later verified by Cliffe (Reference Cliffe1983) in direct numerical simulations. Note that the structure of the present twin-cell mode coexisting with a transient irregular pattern can be studied as a function of the homotopy parameter 
$\epsilon$ in (5.4a,b), which is likely to help in establishing the connection between the two bifurcation diagrams in figure 12 (reflective endwalls) and figure 17 (no-slip endwalls). This, along with the possible effect of asymmetric axial boundary conditions (no-slip at one end and free-slip at the other end (Cole Reference Cole1976; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989)), is left to a future work.
6. Conclusions and outlook
The pattern formation scenario in compressible Taylor–Couette flow (TCF) of a dense gas, with the inner cylinder rotating (
$\omega _i>0$) and the outer one at rest (
$\omega _0=0$), was investigated using large-scale molecular dynamics simulations for the repulsive Lennard-Jones potential embodied with a microscopic damping mechanism (representing a dissipative or athermal system). The results on both primary and secondary bifurcations were found to be qualitatively similar even in the absence of microscopic damping, as demonstrated in appendix A. Most simulations were carried out with reflection boundary conditions at two axial endwalls. This represents symmetric boundary conditions that are expected to yield midplane symmetric Taylor vortices beyond a minimum value of the inner-cylinder rotation (or Reynolds number, (3.3)). The critical Reynolds number for the onset of TVF matched well with that of incompressible TCF since the stratifications remained small in the CCF regime. However, the analogy with incompressible TCF cannot be carried over to the Taylor-vortex regime of a dense gas since the degree of stratification (in density, temperature, angular momentum and related hydrodynamic fields, along both radial and axial directions) increases markedly with increasing rotation, thus making the present system fundamentally different from its incompressible analogue. The axial locations of the maximum and minimum density (and angular momentum) were found to be along the inward and outward jets, respectively, and therefore the outward jets are relatively rarefied or dilute compared to inward jets.
All observed patterns, summarized in figure 4 in terms of a phase diagram in the (
$\omega _i, \varGamma$) plane, are believed to be new in the context of the compressible TCF of a dense athermal gas. Along with even-roll states, we uncovered a single-vortex state and an asymmetric pair of vortices at 
$\varGamma \sim O(1)$, along with three and five vortices at higher values of 
$\varGamma$. For 
$\varGamma =2$ and 
$4$, the primary bifurcation from the CCF states resulted in an array of symmetric Taylor vortices with two and four vortices, respectively. With increasing rotation rate, a secondary transition from a ‘symmetric’ even-roll state to an ‘asymmetric’ odd-roll state was found beyond some critical value of 
$\omega _i$; this transition occurred via a saddle-node bifurcation, leading to hysteresis loops. Such hysteretic transitions were more prominent in short cylinders with 
$\varGamma \leq O(1)$ that admit a single-vortex coexisting with an asymmetric two-roll or symmetric two-roll states, respectively, at smaller and larger values of 
$\omega _i$. In fact, the single-roll branch represented the primary bifurcation from the CCF state for 
$\varGamma \leq 1$, which was found to be connected to the two-roll branch via an unstable branch with two limit points. In the latter case, the transition from an asymmetric two-vortex state to a symmetric two-vortex state was smooth with increasing 
$\omega _i$.
The underlying bifurcations were quantified in terms of the net circulation 
$\zeta$, (4.9), or the maximum radial velocity 
${\rm \Delta} v_r$, (3.1), or the axial density contrast 
${\rm \Delta} \rho$, (4.10), as order parameters. The last parameter (
${\rm \Delta} \rho$) was identified as a key structural signature of compressible TVF, as it characterized the degree of axial density stratification, which became stronger with increasing rotational speed of the inner cylinder. In the Taylor-vortex regime, the radial variation of the specific angular momentum field was found to be non-monotonic, with a centrifugally/Rayleigh-stable part being located near the inner cylinder and a Rayleigh-unstable part near the outer cylinder. The neutral surface, demarcating the centrifugally stable and unstable regimes, gradually shifted towards the outer cylinder with increasing rotation. The density maximum occurring away from the outer cylinder was shown to be a salient feature of the present system that arises due to the ‘denseness’ of the gas, which is in contrast to the TCF of a dilute/ideal gas where the density maximum occurs at the outer cylinder. Since the radial density stratification is tied to the stratification in shear viscosity along the radial direction, the constancy of the shear stress led to the formation of a low-shear region (high viscosity) around the mid-gap, i.e. a plateau on the azimuthal velocity profile, which, in turn, was shown to be responsible for the genesis of non-monotonic radial variations of the angular momentum at large 
$\omega _i\gg \omega ^{cr}$.
While most results were presented for reflecting endwalls, the role of stationary (no-slip) endwalls was checked for a few cases. The noteworthy findings are: (i) the stationary endwalls preferred to pick the so-called ‘normal’ mode vortices (figure 14), having inward jets near the endwalls (figure 16), in agreement with the imperfect bifurcation scenario of Benjamin (Reference Benjamin1978a,Reference Benjaminb); and (ii) in short cylinders (
$\varGamma \leq 1$), the single-roll branch (representing primary bifurcation) disappeared and was replaced by a new state. The latter consisted of two symmetric large vortices (near the inner cylinder) coexisting with an array of transient (smaller and weaker) vortices near the stationary outer cylinder at the onset of instability. With increasing 
$\omega _i$, this state degenerated into (i) an asymmetric two-roll state with two corner eddies and finally to (ii) a symmetric two-roll state at large enough values of 
$\omega _i$. We conclude that the bistability between single-roll and two-roll states (figure 12 for reflective endwalls) is inhibited by stationary endwalls in small-aspect-ratio compressible TCF.
Comparing results between undamped (appendix A) and damped (see figure 3
$b$ and the analysis at the end of § 3.1) cases, we find that the critical Reynolds number for the onset of TVF is larger (
$Re^{cr}(\gamma _n=5)\approx 85.6 > Re^{cr}(\gamma _n=0)\approx 71.4$) in the damped case, for which the density variations are also found to be relatively larger compared to its undamped counterpart (see figures 23
$a$, 24
$a$ and 6
$a$). The latter finding indicates that the effect of compressibility is, in general, to stabilize the TCF of a dense gas. This is similar to the linear stability results (Manela & Frankel Reference Manela and Frankel2007; Welsh et al. Reference Welsh, Kersal'e and Jones2014) and direct numerical simulations (Aghor & Alam Reference Aghor and Alam2020) for a ‘dilute’ gas undergoing TCF. A systematic study to reveal compressibility effects on the genesis of Taylor vortices would require additional simulations by varying the Mach number, which can be achieved in MD simulations by changing the mean density (
$\rho =0.5$ in present simulations) – this can be taken up in a future work. A crucial effect of the microscopic damping force is to make the secondary transition to asymmetric vortices (figures 11 and 20) into a disconnected branch that coexists with the primary branch over a range of 
$\omega _i$ as discussed in appendix A.
The implications of the observed radial and axial stratifications in hydrodynamic fields (such as density and temperature) on the stability of Taylor vortices at large enough 
$Re$ have not been investigated; this can also be done from a theoretical analysis of continuum equations for dense gases and would be an interesting avenue for future work. In this regard, a complementary work on the TCF of a ‘dilute’ gas is currently under consideration (Aghor & Alam Reference Aghor and Alam2020). Apart from the fluid being a dilute molecular gas, the latter work also differs from the present work on four major points: (i) its solves the compressible Navier–Stokes–Fourier equations numerically via direct numerical simulation; (ii) it focuses solely on the effect of ‘stationary’ endwalls on flow multiplicity and bifurcation; (iii) the effect of Mach number is studied, which confirmed its stabilizing role; and (iv) the ramping protocol of increasing/decreasing the cylinder height, as in the experiments of Benjamin (Reference Benjamin1978a,Reference Benjaminb), is followed, which is difficult to implement in MD simulations with a pre-specified mean density. From the viewpoint of MD simulations, it would be worthwhile to map out the phase diagram such as figure 4 in the (
$\omega _i, \omega _o$) plane by considering both counter-rotation and co-rotation of the two cylinders (Coles Reference Coles1965; Andereck et al. Reference Andereck, Liu and Swinney1986) as functions of (
$\varGamma , \eta$) in compressible TCF.
Figure 20. Same as figure 11(
$a{,}b$), but with 
$\gamma _n=0$. The squares and pentagons denote four-roll and five-roll states, respectively, that have been obtained by slowly increasing the inner-cylinder rotation rate from rest to a specified value 
$\omega _i$; the solid lines, representing stable solutions, and the dashed line (conjectured), representing unstable solutions, are drawn to guide the eye. The filled triangles have been obtained from an up-sweep ‘
$\omega _i$ continuation’ run as described in figure 21, while the inverted triangles correspond to a ‘down-sweep’ run starting from 
$\omega _i=0.25$.
Acknowledgements
N.G. was supported, in part, by the Science and Engineering Research Board (SERB) National Post Doctoral Fellowship (Ref.: PDF/2015/001068) and would like to gratefully acknowledge the same. We thank Mr A. Mahajan for preliminary works on this topic using hard-sphere simulations. One of us (M.A.) acknowledges Professor P. R. Nott for suggesting relevant references on single-vortex mode.
Declaration of interests
The authors report no conflict of interests.
Appendix A. Role of damping force on pattern transitions and hydrodynamic fields
All the results presented in the main text (except in figure 2
$b$) correspond to a finite value of the microscopic damping coefficient, 
$\gamma _n=5$, and the corresponding damping force is given by 
$F_d=-m\gamma _n v_n$, with 
$v_n$ being the relative velocity of the colliding pair along its contact vector. As discussed in § 2, 
$\gamma _n$ can be tied to a ‘velocity-dependent’ normal restitution coefficient and hence the present microscopic model is akin to a dense ‘granular’ gas with a velocity-dependent restitution coefficient. To isolate the effect of microscopic damping, especially on the onset of asymmetric vortices and the underlying hysteretic bifurcations, the simulations at 
$\varGamma =3/4, 1, 2$ and 
$4$ (viz. figure 4) were repeated over a range of rotation rates 
$\omega _i\leq 0.4$ by setting 
$\gamma _n=0$, which represents a dense ‘molecular’ gas. In general, we found that the onset of bifurcation is delayed at 
$\gamma _n=0$ but (i) the nature of the primary bifurcation as well as (ii) the hysteretic transition to asymmetric states remain qualitatively similar to those for 
$\gamma _n=5$. A comparative discussion of these results, along with an explanation for the radial temperature profiles in figure 6
$(b)$, is presented below.
Figure 20 displays the bifurcation diagrams in the 
$(a)$ (
${\rm \Delta} v_r, \omega _i$) and 
$(b)$ (
$|\zeta |, \omega _i$) planes for 
$\varGamma =4$ and 
$\gamma _n=0$ – these should be compared with figure 11
$(a{,}b)$, which holds for the damped case (
$\gamma _n=5$). The major finding is that the critical rotation rate for the onset of Taylor vortices is 
$\omega _i^{cr}\approx 0.07$, which is larger than that in the presence of finite damping, 
$\omega ^{cr} (\gamma _n=5)\approx 0.054$ (viz. figure 4). As discussed in § 3.1, for this case, the average temperature is 
$T_{av}^*\approx 0.403$ and hence the critical Reynolds number, (3.5), is 
$Re^{cr} (\gamma _n=0)\approx 71.4 < Re^{cr} (\gamma _n=5)=85.6$. That the normal damping force increases the value of the critical Reynolds number has been verified for the range of 
$\varGamma \leq 4$ studied.
Comparing figures 20
$(b)$ and 11
$(b)$, we find that the five-roll branch (denoted by pentagons and triangles) seems to be disconnected from the four-roll branch having a lower limit point at 
$\omega _i \approx 0.08$. On the other hand, in contrast to the damped case (
$\gamma _n=5$, figure 11
$b$), we did not find six-roll states over the range of 
$\omega _i\leq 0.4$ – possibly the onset of 4-rolls 
$\to$ 6-rolls transition is delayed (
$\omega _i>0.4$) for the undamped case. It should be noted that, while the squares (representing four-roll states) and pentagons (five-roll state) in figure 20 have been obtained by running simulations at each 
$\omega _i$ from the rest state as described in the last paragraph in § 2.2, the triangles have been obtained from two independent (up-sweep and down-sweep) ‘
$\omega _i$ continuation’ runs. The ‘up-sweep’ protocol is summarized in figure 21. First, the inner cylinder is uniformly accelerated over 
$10^5{\rm \Delta} t$ (where 
${\rm \Delta} t$ is the time step) from the rest state to the prescribed rotation rate 
$\omega _i=0.25$, and the coarse-grained fields are evaluated next over a time window of 
$5\times 10^5{\rm \Delta} t$. The ‘up-sweep’ 
$\omega _i$ continuation protocol then sets in, wherein the inner cylinder is uniformly accelerated from 
$\omega _i=0.25$, over the next 
$10^4$ time steps, to 
$\omega _i=0.26$, which is maintained over another time window of 
$4\times 10^5{\rm \Delta} t$ to evaluate hydrodynamic fields before increasing further to 
$\omega _i=0.27$. This procedure is continued up to 
$\omega _i=0.35$ in steps of 
${\rm \Delta} \omega _i=0.01$; see figure 21 for details. The related ‘down-sweep’ protocol is the same but with decreasing rotation rate from 
$\omega _i=0.25$, and these data are marked by inverted triangles in figure 20
$(a{,}b)$. The latter figure confirms that the four- and five-roll states can coexist with each other over a large range of 
$\omega _i\in (0.08, 0.4)$ even in the absence of damping (
$\gamma _n=0$), with the left limit point being located at 
$\omega _i \approx 0.08$ – these overall findings are similar to those for the damped case in figure 11. Note that we have not been able to identify the right limit point (if any) in figure 20, since these simulations are extremely time-consuming, especially for 
$\gamma _n=0$ at 
$\omega _i>0.3$ that requires a much smaller time step (
${\rm \Delta} t \leq 0.0025$) for convergence – this can be studied in future.
Figure 21. Evolution of system temperature (2.10) for the ‘
$\omega _i$ continuation’ run over 
$\omega _i\in (0.25,0.35)$, with parameter values as in figure 20. The inner-cylinder rotation rate 
$\omega _i$ is increased quasi-statically from 
$\omega _i=0.25$ to 
$0.35$ at a ramping rate of 
$\textrm {d}\omega _i/\textrm {d}t= 2\times 10^{-5}$. The two insets display the zoomed part of the ramp-up region when the inner cylinder is accelerated uniformly over 
$10^4$ time steps to reach 
$\omega _i=0.26$ (left inset) and 
$\omega _i=0.32$ (right inset).
Figure 22 shows a comparison between updamped (
$\gamma _n=0$) and damped (
$\gamma _n=5$ as in the main text) cases for four-roll (panels 
$a$–
$c$ at 
$\omega _i=0.1$) and five-roll (panels 
$d$,
$e$ at 
$\omega _i=0.28$) states. The corresponding (axially averaged) radial profiles of density, temperature and azimuthal velocity are compared in figure 23(
$a$–
$c$) for both symmetric and asymmetric vortices. The axial–radial velocity vectors in figure 22(
$a{,}d$) and the density maps in figure 22(
$b{,}e$) in the meridional plane look structurally similar for 
$\gamma _n=0$ and 
$5$. The temperature fields in figure 22(
$c{,}f$) look qualitatively different: while the temperature maximum occurs within the annular gap at any 
$\omega _i$ with 
$\gamma _n=0$, the temperature is maximum near the inner cylinder and decays towards the outer cylinder at larger rotation rates for 
$\gamma _n=5$. The above characteristic features are made clearer in figure 23
$(b)$, which reconfirms that the role of microscopic dissipation (
$\gamma _n>0$) is to qualitatively change the concave shape of 
$T(r)$ with its maximum in the annular gap to a convex 
$T(r)$ profile having a minimum within the gap at 
$\omega _i=0.1$ that gradually shifts towards the outer cylinder with increasing 
$\omega _i=0.28$. Recall from § 4.1 that the latter findings are similar to those in figure 6
$(b)$ that holds for 
$\varGamma =2$ with 
$\gamma _n=5$. A closer look at figure 22
$(c,\,f)$ reveals that the fluid is in thermal equilibrium with both walls (
$T_w=1\approx T_{fluid}(r=r_i, r_o)$ at both cylinders) for the undamped case (
$\gamma _n=0$), but there is a considerable temperature slip (
$T_w\neq T_{fluid}(r=r_i, r_o)$) in the presence of finite damping force (
$\gamma _n=5$).
Figure 22. Comparison between undamped (
$\gamma _n=0$, left part of each pair) and damped (
$\gamma _n=5$, right part) cases for 
$(a{,}d)$ radial–axial velocity map, 
$(b{,}e)$ density and 
$(c{,}f)$ temperature in the (
$r,z$) plane. The upper panels (
$a$–
$c$) represent a four-roll state at 
$\omega _i=0.1$ and the lower panels (
$d$–
$f$) represent a five-roll state at 
$\omega _i=0.28$. The aspect ratio is 
$\varGamma =4$. Refer to the bifurcation diagrams in figures 11 and 20 for 
$\gamma _n=5$ and 
$\gamma _n=0$, respectively.
Figure 23. Comparison between undamped (
$\gamma _n=0$, filled symbols) and damped (
$\gamma _n=5$, open symbols) cases for the radial profiles of 
$(a)$ reduced density, 
$(b)$ temperature (
$T^*=T/(\omega _i\delta )^2$) and 
$(c)$ azimuthal velocity; other parameters as in figure 22.
Because of the thermal equilibrium of the fluid with two isothermal walls in the absence of any damping mechanism, the added shear work (
$\boldsymbol{\mathsf{P}}\,\boldsymbol {:}\,{\boldsymbol \nabla }{\boldsymbol u}$, where 
$\boldsymbol{\mathsf{P}}$ is the stress tensor) due to the rotating inner cylinder naturally leads to a rise of the fluid temperature in the annular gap, resulting in a concave 
$T(r)$ profile. This argument holds in both dilute (Welsh et al. Reference Welsh, Kersal'e and Jones2014; Aghor & Alam Reference Aghor and Alam2020) and dense gases wherein the thermal diffusion (
${\boldsymbol \nabla }\boldsymbol {\cdot }(\kappa (\rho , T){\boldsymbol \nabla }T)$, where 
$\kappa (\rho ,T)$ is thermal conductivity) balances the shear work term. In contrast, the presence of a microscopic damping force (
$\gamma _n>0$) makes the fluid ‘athermal’ (such as an inelastic or granular gas) for which the energy equation (Goldhirsch Reference Goldhirsch2003; Saha & Alam Reference Saha and Alam2016),
\begin{align} {\boldsymbol\nabla}\boldsymbol{\cdot}\left(\kappa(\rho, T){\boldsymbol\nabla}T\right) + \boldsymbol{\mathsf{P}}\,\boldsymbol{:}\,{\boldsymbol\nabla}{\boldsymbol u} - \gamma_n{\mathcal{D}}(\rho, T) =0, \end{align}
contains an additional sink term (
$\mathcal {D}$) that takes care of the cooling mechanism due to the finite microscopic damping force (
$\gamma _n>0$). Clearly, the last term in (A 1) must dominate over the shear work term for the decrease of the gas temperature in the bulk (such as in figure 6
$b$). Therefore, we conclude that the competition among the three mechanisms in (A 1) is responsible for the anomalous radial temperature profiles found in figures 23
$(b)$ and 6
$(b)$ for 
$\gamma _n=5$. A detailed theoretical analysis would require a solution of (A 1) along with the momentum equation and appropriate boundary conditions, and is beyond the scope of the present paper.
The comparison of density profiles in figure 23
$(a)$ between 
$\gamma _n=0$ and 
$5$ indicates that the density variations are relatively milder when the damping force is absent (
$\gamma _n=0$). This may be understood by drawing an analogy of the present damped system with a granular gas (or a gas–solid suspension), which is known to readily admit strong density inhomogeneities due to a dissipation-driven clustering instability. The axially averaged azimuthal velocity profiles in figure 23
$(c)$ are found to develop a plateau with increasing 
$\omega _i$ even at 
$\gamma _n=0$ – this is similar to those discussed in figures 6
$(c)$ and 9
$(b)$ for 
$\gamma _n=5$.
The above characteristic features about the radial variations of density, temperature and azimuthal velocity can be further confirmed from a comparison between figures 24(
$a$–
$c$) and 6(
$a$–
$c$) that hold for 
$\gamma _n=0$ and 
$5$, respectively, with 
$\varGamma =2$. Note, in particular, that the radial temperature profiles in figure 24
$(b)$ possess a maximum within the annular gap in contrast to a temperature minimum in figure 6
$(b)$ for the damped case. Therefore, the athermal nature (
$\gamma _n>0$) of the gas is responsible for the qualitative changes in the radial variations of the kinetic temperature.
Figure 24. Same as in figure 6(
$a$–
$c$) but with 
$\gamma _n=0$: 
$(a)$
$\rho (r)$, 
$(b)$
$T^*(r) = T(r)/(\omega _i\delta )^2$ and 
$(c)$
$v_\theta ^*(r) = v_\theta (r)/(\omega _i\delta )$.
Lastly, recall from the bifurcation diagram in figure 12 that single-roll and two-roll states can coexist with each other over a large range of 
$\omega _i$ in the presence of a normal damping force (
$\gamma _n=5$). Additional simulations at 
$\varGamma =1$ with 
$\gamma _n=0$ confirmed that the overall bifurcation diagram remains similar (not shown). Figure 25(
$a$–
$f$) displays a comparison between undamped (
$\gamma _n=0$) and damped (
$\gamma _n=5$) cases for single-roll (panels 
$a{,}c{,}e$ at 
$\omega _i=0.15$) and two-roll (panels 
$b{,}d{,}f$ at 
$\omega _i=0.4$) states. While the axial–radial velocity maps are strikingly similar, the differences show up in the density and temperature maps in the meridional plane as discussed in the preceding paragraphs.
Figure 25. Comparison between undamped (
$\gamma _n=0$, left part of each pair) and damped (
$\gamma _n=5$, right part) cases for 
$(a{,}b)$ radial–axial velocity map, 
$(c{,}d)$ density and 
$(e{,}f)$ temperature in the (
$r,z$) plane. While panels 
$(a{,}c{,}e)$ represent a one-roll state at 
$\omega _i=0.15$, panels (
$b{,}d{,}f$) denote two-roll states at 
$\omega _i=0.4$. The aspect ratio is 
$\varGamma =1$. Refer to the corresponding bifurcation diagram in figure 12 for 
$\gamma _n=5$.
