On the Stability of a Rotating Cylindrical Shear Flow: Consequences for Differentially Rotating Vortical Structures

1. Introduction

Direct numerical simulations of a variety of rotating turbulent systems, such as convection (Guervilly et al. 1, Julien et al. 2), stratified turbulence (Marino et al. 3, Oks et al. 4) and double-diffusive convection (Moll and Garaud 5, Sengupta and Garaud 6) have reported the emergence of large scale vortices (LSVs) that is indicative of a basic underlying instability that works towards sustaining these coherent structures in rotating environments. In an attempt to possibly identify this inherent mechanism that drives LSVs in rotating turbulence, we present a simplified analytical model of a vertical shear flow in a uniformly rotating system and study the growth of interfacial instabilities across the boundary of the flow. Through an asymptotic method, we investigate the linear stability of the flow and quantify the radial flux of axial angular momentum to predict regimes that could possibly spin up the core of the flow.

The stability of a column of fluid in a rotating environment has been a topic of many studies beginning with the early works of Hocking and Michael [7], Ponstein [8], Gillis [9], Alterman [10], Pedley [11] to later investigations using linear stability theory ([12,13,14]) as well as experiments reporting helical modes of instabilities in the context of rotating jets ([15]). Most of these works investigated the development of the capillary instabilities which are dominated by effects of surface tension balancing the outward centrifugal force leading to the saturation of the growing instability. The effect of viscosity has also been investigated on the growth and saturation of these instabilities by Hocking [16], Kubitschek and Weidman [17], Coppola and Semeraro [18]. The effect of a uniform circulation surrounding a cylindrical vortex sheet has been explored by Rotunno [19] with linear theory predictions for growth of unstable modes. In this work, we explore the effect of a uniform background rotation on the stability of a analogous flow as a simplified model for studying the growth of large scale vortices (LSVs) reported in recent numerical studies and also observed on various scales in nature. Recent efforts to explain the formation of such patterns in the oceanic context have concluded of its dependence to relevant parameters ([20]) albeit with semi-analytic or two-dimensional numerical simulations. This work attempts to simultaneously present DNSs of a idealized model for a vortex in form of a shear flow in a uniformly rotating triply-periodic domain.

In Section 2, we present a linear stability analysis of our toy model setup in a perturbative approach to deal with the full rotating problem in order to estimate the effect of the instability on the angular momentum budget of the shear flow. Section 3 presents numerical calculations incorporating effects of viscosity (for moderate values of Reynolds numbers) to test the predictions of our linear theory and also predicts the final non-linear outcome of this simplified picture of our vortex model. We conclude in Section 4 by discussing our findings and its implications for generation of LSVs in rotating environments and suggest improvements possible to make our proposed model approach more realistic conditions encountered for these environments

2. Linear Stability of a Rotating Cylindrical Shear Flow

We begin by investigating the linear stability of a simple cylindrical shear flow in a rotating incompressible fluid, ignoring effects of viscosity and stratification. We consider a fluid rotating at an angular velocity ${\Omega }_0$, whose direction defines the vertical (z-axis). Within that fluid, we consider a cylindrical region of radius R, aligned in the z-direction, centered at $(x,y)=(0,0)$. Within this cylinder, the fluid is assumed to rotate with an angular velocity, ${\Omega }_{+}$ (with respect to the background), and flow along the z- direction uniformly with velocity $u_{+}$. Outside the cylinder, the fluid has vertical velocity $u_{-}$, and is rotating at a rate ${\Omega }_{-}$ with respect to the background (see Figure 1). We calculate the unstable eigenmodes associated with this flow field and study their role in the evolution of the net angular momentum of the fluid in the cylinder. For simplicity, we assume the system is periodic in z, and study it in the interval $[0,L_z]$ only.

2.1. Governing Equations

For the model described above, the background velocity of the fluid (on top of the global rotation) is given by

$$\begin{array}{cc}\hfill \overline{\mathbf{u}}& =(0,{\overline{u}}_{\varphi },{\overline{u}}_z)=\left\{\begin{array}{cc}(0,r{\Omega }_{+},u_{+})\hfill & r\le R\hfill \\ (0,r{\Omega }_{-},u_{-})\hfill & r>R\hfill \end{array}\right.\hfill \end{array}$$

(1)

where $r=\sqrt{x^2+y^2}$ is the cylindrical radius. To ensure that this flow is in a steady state, the background pressure $\overline{p}$ must satisfy

$${\displaystyle \frac{\partial {\overline{p}}_{\pm }}{\partial r}}=r{({\Omega }_0+{\Omega }_{\pm })}^2,$$

(2)

with the + and − signs corresponding to the regions inside and outside $r=R$, respectively. As a result, note that the radial pressure gradient is discontinuous at $r=R$ when ${\Omega }_{+}\ne {\Omega }_{-}$, but $\overline{p}$ can be made continuous by suitable matching i.e.

$${\overline{p}}_{+}\left(R\right)={\overline{p}}_{-}\left(R\right).$$

(3)

We now consider perturbations $\mathbf{u}$ and p to the flow and pressure respectively satisfying $\mathbf{u}\ll \overline{\mathbf{u}}$ and $p\ll \overline{p}$. Linearizing the full equations around the background state and assuming incompressibility for the flow, the perturbation equations are given by

$$\begin{array}{c}{\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·\nabla \overline{\mathbf{u}}+\overline{\mathbf{u}}·\nabla \mathbf{u}+2{\Omega }_0{\widehat{\mathbf{e}}}_z\times \mathbf{u}=-\nabla p\end{array}$$

(4a)

$$\begin{array}{c}\nabla ·\mathbf{u}=0,\end{array}$$

(4b)

where, $\mathbf{u}=\left(u_r,u_{\varphi },u_z\right)$ in cylindrical coordinates $\left(r,\varphi ,z\right)$ and as mentioned earlier, viscosity has been ignored. Using the following ansatz1 for each of the velocity components and for the pressure p

$$q(r,\varphi ,z,t)=Re\left[\tilde{q}\left(r\right)exp\left(i(kz+m\varphi )+\lambda t\right)\right],$$

(6)

where $k={\displaystyle \frac{2\pi n}{L_z}},$ and both m and n are both integers, we can reduce the set of perturbation equations within and outside the cylinder to

$${\Lambda }_{\pm }{\tilde{u}}_{r_{\pm }}-2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{{\varphi }_{\pm }}=-{\displaystyle \frac{d{\tilde{p}}_{\pm }}{dr}}$$

(7a)

$${\Lambda }_{\pm }{\tilde{u}}_{{\varphi }_{\pm }}+2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{r_{\pm }}=-{\displaystyle \frac{im{\tilde{p}}_{\pm }}{r}}$$

(7b)

$${\Lambda }_{\pm }{\tilde{u}}_{z_{\pm }}=-ik{\tilde{p}}_{\pm }$$

(7c)

$${\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left(r{\tilde{u}}_{r_{\pm }}\right)+{\displaystyle \frac{im}{r}}{\tilde{u}}_{{\varphi }_{\pm }}+ik{\tilde{u}}_{z_{\pm }}=0.$$

(7d)

with

$${\Lambda }_{\pm }=\lambda +im{\Omega }_{\pm }+ik{u}_{\pm },$$

(8)

where, as before, the + and − signs correspond to inner ($r<R$) and outer ($r>R$) regions respectively. Taking the complex conjugate of the above equations, we can readily see that if a solution with growth rate $\lambda (m,k)$ satisfies these equations, then its complex conjugate also satisfies the same equations with growth rate $\lambda (-m,-k)$. This implies that it suffices to describe the unstable modes in half of the $k-m$ parameter space. Hence, we restrict our study to $k\ge 0$ but allow m to take both positive and negative values, as is commonly done in the literature in studies of temporal flow instability (e.g. [21]).

The first three equations can be combined to express the velocity fields in terms of the pressure perturbation as

$$\begin{array}{c}{\tilde{u}}_{r_{\pm }}\left(r\right)={\displaystyle \frac{-{\Lambda }_{\pm }\frac{\partial {\tilde{p}}_{\pm }}{\partial r}-2\frac{im}{r}\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{p}}_{\pm }}{{\Lambda }_{\pm }^2+4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}},\end{array}$$

(9a)

$$\begin{array}{c}{\tilde{u}}_{{\varphi }_{\pm }}\left(r\right)={\displaystyle \frac{2\left({\Omega }_0+{\Omega }_{\pm }\right)\frac{\partial {\tilde{p}}_{\pm }}{\partial r}-\frac{im}{r}{\Lambda }_{\pm }{\tilde{p}}_{\pm }}{{\Lambda }_{\pm }^2+4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}},\end{array}$$

(9b)

$$\begin{array}{c}{\tilde{u}}_{z_{\pm }}\left(r\right)=-{\displaystyle \frac{ik}{{\Lambda }_{\pm }}}{\tilde{p}}_{\pm }.\end{array}$$

(9c)

Using these expressions in the incompressibility condition (7d) yields an equation for the pressure, namely

$$r^2{\displaystyle \frac{d^2{\tilde{p}}_{\pm }}{d{r}^2}}+r{\displaystyle \frac{d{\tilde{p}}_{\pm }}{dr}}-(m^2+{\mu }_{\pm }^2{r}^2){\tilde{p}}_{\pm }=0,$$

(10)

with

$${\mu }_{\pm }^2=k^2\left[1+{\displaystyle \frac{4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}{{\Lambda }_{\pm }^2}}\right].$$

(11)

By imposing that the pressure has to be regular at $r=0$ and tends to zero for large r, we can write solutions in terms of the modified Bessel functions given by

$$\begin{array}{cc}\hfill \tilde{p}\left(r\right)& ={\tilde{p}}_{+}\left(r\right)={\widehat{p}}_{+}{I}_m\left(r{\mu }_{+}\right),forr<R\hfill \\ & ={\tilde{p}}_{-}\left(r\right)={\widehat{p}}_{-}{K}_m\left(r{\mu }_{-}\right),forr>R\hfill \end{array}$$

(12)

where, without loss of generality, we selected the root of Eq 10 so that $Re\left({\mu }_{\pm }\right)>0$. Indeed, using the parity properties of the Bessel functions, it can be readily shown that solutions for cases with $Re\left({\mu }_{\pm }\right)<0$ are obtained simply by reversing the sign of ${\mu }_{\pm }$. Using these solutions, and the standard recursion relations for the derivatives of the Bessel functions given by

$${\displaystyle \frac{d{I}_m\left(r{\mu }_{+}\right)}{dr}}={\displaystyle \frac{m}{r}}{I}_m\left(r{\mu }_{+}\right)+{\mu }_{+}{I}_{m+1}\left(r{\mu }_{+}\right),$$

(13a)

$${\displaystyle \frac{d{K}_m\left(r{\mu }_{-}\right)}{dr}}={\displaystyle \frac{m}{r}}{K}_m\left(r{\mu }_{-}\right)-{\mu }_{-}{K}_{m+1}\left(r{\mu }_{-}\right),$$

(13b)

we define

$$\begin{array}{cc}\hfill {\mathcal{I}}_m\left(r{\mu }_{+}\right)& ={\displaystyle \frac{1}{{\mu }_{+}{I}_m\left(r{\mu }_{+}\right)}}{\displaystyle \frac{d{I}_m\left(r{\mu }_{+}\right)}{dr}}={\displaystyle \frac{I_{m+1}\left(r{\mu }_{+}\right)}{I_m\left(r{\mu }_{+}\right)}}+{\displaystyle \frac{m}{r{\mu }_{+}}}\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill {\mathcal{K}}_m\left(r{\mu }_{-}\right)& =-{\displaystyle \frac{1}{{\mu }_{-}{K}_m\left(r{\mu }_{-}\right)}}{\displaystyle \frac{d{K}_m\left(r{\mu }_{-}\right)}{dr}}={\displaystyle \frac{K_{m+1}\left(r{\mu }_{-}\right)}{K_m\left(r{\mu }_{+}\right)}}-{\displaystyle \frac{m}{r{\mu }_{-}}}.\hfill \end{array}$$

(14b)

The solutions for the velocity components can then be expressed as

$$\begin{array}{c}{\tilde{u}}_r\left(r\right)=\left\{\begin{array}{cc}{\tilde{u}}_{r_{+}}\left(r\right)={\displaystyle \frac{k^2}{{\Lambda }_{+}^2{\mu }_{+}^2}}\left[-{\Lambda }_{+}{\mu }_{+}{\mathcal{I}}_m\left(r{\mu }_{+}\right)-{\displaystyle \frac{2im\left({\Omega }_0+{\Omega }_{+}\right)}{r}}\right]{\widehat{p}}_{+}{I}_m\left(r{\mu }_{+}\right)\hfill & r<R\hfill \\ {\tilde{u}}_{r_{-}}\left(r\right)={\displaystyle \frac{k^2}{{\Lambda }_{-}^2{\mu }_{-}^2}}\left[{\Lambda }_{-}{\mu }_{-}{\mathcal{K}}_m\left(r{\mu }_{-}\right)-{\displaystyle \frac{2im\left({\Omega }_0+{\Omega }_{-}\right)}{r}}\right]{\widehat{p}}_{-}{K}_m\left(r{\mu }_{-}\right),\hfill & r>R\hfill \end{array}\right.\end{array}$$

(15a)

$$\begin{array}{c}{\tilde{u}}_{\varphi }\left(r\right)=\left\{\begin{array}{cc}{\tilde{u}}_{{\varphi }_{+}}\left(r\right)={\displaystyle \frac{k^2}{{\Lambda }_{+}^2{\mu }_{+}^2}}\left[2\left({\Omega }_0+{\Omega }_{+}\right){\mu }_{+}{\mathcal{I}}_m\left(r{\mu }_{+}\right)-{\displaystyle \frac{im}{r}}{\Lambda }_{+}\right]{\widehat{p}}_{+}{I}_m\left(r{\mu }_{+}\right)\hfill & r<R\hfill \\ {\tilde{u}}_{{\varphi }_{-}}\left(r\right)={\displaystyle \frac{k^2}{{\Lambda }_{-}^2{\mu }_{-}^2}}\left[-2\left({\Omega }_0+{\Omega }_{-}\right){\mu }_{-}{\mathcal{K}}_m\left(r{\mu }_{-}\right)-{\displaystyle \frac{im}{r}}{\Lambda }_{-}\right]{\widehat{p}}_{-}{K}_m\left(r{\mu }_{-}\right),\hfill & r>R\hfill \end{array}\right.\end{array}$$

(15b)

$$\begin{array}{c}{\tilde{u}}_z\left(r\right)=\left\{\begin{array}{cc}{\tilde{u}}_{z_{+}}\left(r\right)=-{\displaystyle \frac{ik}{{\Lambda }_{+}}}{\widehat{p}}_{+}{I}_m\left(r{\mu }_{+}\right)\hfill & r<R\hfill \\ {\tilde{u}}_{z_{-}}\left(r\right)=-{\displaystyle \frac{ik}{{\Lambda }_{-}}}{\widehat{p}}_{-}{K}_m\left(r{\mu }_{-}\right),\hfill & r>R.\hfill \end{array}\right.\end{array}$$

(15c)

The so far unknown amplitudes ${\widehat{p}}_{+}$ and ${\widehat{p}}_{-}$ are related to one another through the matching conditions at the interface between the inner and outer fluid. At any point in time, this interface is located at the position $r_i(\varphi ,z,t)=R+\eta (\varphi ,z,t),$ where $\eta$ is its radial displacement with respect to the rest position at $r=R$.

The kinematic condition relates the function $\eta$ to the radial velocity of the fluid at any point along the interface, and implies

$$u_r(R+\eta ,\varphi ,z,t)={\displaystyle \frac{D\eta }{Dt}}={\displaystyle \frac{\partial \eta }{\partial t}}+{\displaystyle \frac{u_{\varphi }}{r}}{\displaystyle \frac{\partial \eta }{\partial \varphi }}+u_z{\displaystyle \frac{\partial \eta }{\partial z}},$$

(16)

Linearizing this equation for small displacements, we get

$$u_r(R,\varphi ,z,t)=\left({\displaystyle \frac{\partial }{\partial t}}+{\displaystyle \frac{{\overline{u}}_{\varphi }}{r}}{\displaystyle \frac{\partial }{\partial \varphi }}+{\overline{u}}_z{\displaystyle \frac{\partial }{\partial z}}\right)\eta ,$$

(17)

which, using the ansatz (6) and a similar one for $\eta$, namely

$$\eta (\varphi ,z,t)=Re\left[\widehat{\eta }exp\left(i(kz+m\varphi )+\lambda t\right)\right],$$

(18)

reduces to

$${\Lambda }_{\pm }\widehat{\eta }={\tilde{u}}_{r_{\pm }}\left(R\right),$$

(19)

where, for simplicity, we define

$${\tilde{u}}_{r_{\pm }}\left(R\right)={\displaystyle \genfrac{}{}{0pt}{}{lim}{r\to R_{\pm }}}{\tilde{u}}_r\left(r\right)$$

(20)

Note that without loss of generality, we can chose $\widehat{\eta }\in {\mathbb{R}}^{+}$ so that

$$\eta =\widehat{\eta }cos\left(m\varphi +kz+{\lambda }_{I}t\right)exp\left({\lambda }_{R}t\right),$$

(21)

where

$$\lambda ={\lambda }_R+i{\lambda }_I.$$

(22)

As a result, using Eq 8, we see that

$$Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]=Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]={\lambda }_R\widehat{\eta }.$$

(23)

This expression, when combined with Eq 15, provides two equations that effectively relate $\widehat{\eta }$ to ${\widehat{p}}_{+}$ and ${\widehat{p}}_{-}$, respectively.

A third equation between $\widehat{\eta }$, ${\widehat{p}}_{+}$ and ${\widehat{p}}_{-}$ can be obtained from the only other allowable interfacial matching condition in an inviscid fluid, namely the matching of the pressure. Recall however that, for a differentially rotating fluid (i.e. ${\Omega }_{+}\ne {\Omega }_{-}$), the local centrifugal force induces a jump in the background pressure gradient at $r=R$ (see Eq 2) which has to be taken into account in the matching procedure. We therefore match the total pressure $p_{tot}=\overline{p}\left(r\right)+p(r,\varphi ,z,t)$ at the interface located at $r=R+\eta$. After linearization, this implies the continuity of $\overline{p}\left(R\right)+{\left.{\displaystyle \frac{\partial \overline{p}}{\partial r}}\right|}_{r=R}\eta +p(R,\varphi ,z,t)$ across the interface, which (using Eq 2 and the ansatz for p and $\eta$) becomes

$${\widehat{p}}_{+}{I}_m\left(R{\mu }_{+}\right)+{\left({\Omega }_0+{\Omega }_{+}\right)}^{2}R\widehat{\eta }={\widehat{p}}_{-}{K}_m\left(R{\mu }_{-}\right)+{\left({\Omega }_0+{\Omega }_{-}\right)}^{2}R\widehat{\eta }.$$

(24)

The second term on both sides is the contribution from the centrifugal force, which cancels out unless ${\Omega }_{+}\ne {\Omega }_{-}$.

Combining Eqs 15a and (24), we can now eliminate all unknown amplitudes $\widehat{\eta }$, ${\widehat{p}}_{+}$ and ${\widehat{p}}_{-}$ to obtain the dispersion relation for the growth rate $\lambda$ as

$${\displaystyle \frac{\frac{{\mathcal{I}}_m\left(R{\mu }_{+}\right)}{{\mu }_{+}{\Lambda }_{+}^2}+\frac{2im({\Omega }_0+{\Omega }_{+})}{R{\Lambda }_{+}^3{\mu }_{+}^2}+\frac{{\mathcal{K}}_m\left(R{\mu }_{-}\right)}{{\mu }_{-}{\Lambda }_{-}^2}-\frac{2im({\Omega }_0+{\Omega }_{-})}{R{\Lambda }_{-}^3{\mu }_{-}^2}}{\left[\frac{{\mathcal{I}}_m\left(R{\mu }_{+}\right)}{{\mu }_{+}{\Lambda }_{+}^2}+\frac{2im({\Omega }_0+{\Omega }_{+})}{R{\Lambda }_{+}^3{\mu }_{+}^2}\right]\left[-\frac{{\mathcal{K}}_m\left(R{\mu }_{-}\right)}{{\mu }_{-}{\Lambda }_{-}^2}+\frac{2im({\Omega }_0+{\Omega }_{-})}{R{\Lambda }_{-}^3{\mu }_{-}^2}\right]}}=R{k}^2\left[{\left({\Omega }_0+{\Omega }_{-}\right)}^2-{\left({\Omega }_0+{\Omega }_{+}\right)}^2\right].$$

(25)

This equation must usually be solved numerically for growing modes ($Re\left(\lambda \right)>0$) to reveal different regimes of instability in the flow. However, in some limiting cases, analytical or semi-analytical solutions exits which are discussed next.

2.2. Limiting Cases:

2.2.1. Non-Rotating Case (${\Omega }_{+}={\Omega }_{-}={\Omega }_0=0$) :

Without any background or differential rotation of the tube, we have ${\mu }_{\pm }=k$; so, the dispersion relation becomes

$${\displaystyle \frac{{\mathcal{I}}_m^{\left(0\right)}}{{\Lambda }_{+}^{\left(0\right)2}}}+{\displaystyle \frac{{\mathcal{K}}_m^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)2}}}=0,$$

(26)

with the superscript denoting quantities associated with the non-rotating flow, i.e.

$$\begin{array}{cc}\hfill {\mathcal{I}}_m^{\left(0\right)}& ={\mathcal{I}}_m\left(Rk\right),\hfill \\ \hfill {\mathcal{K}}_m^{\left(0\right)}& ={\mathcal{K}}_m\left(Rk\right),\hfill \\ \hfill {\Lambda }_{\pm }^{\left(0\right)}& ={\lambda }^{\left(0\right)}+ik{u}_{\pm }.\hfill \end{array}$$

Eq (26) can now be solved analytically for ${\lambda }^{\left(0\right)}$, and the solution corresponding to the growing mode of instability is

$${\lambda }^{\left(0\right)}=-ik{\displaystyle \frac{u_{-}{\mathcal{I}}_m^{\left(0\right)}+u_{+}{\mathcal{K}}_m^{\left(0\right)}}{{\mathcal{I}}_m^{\left(0\right)}+{\mathcal{K}}_m^{\left(0\right)}}}+{\displaystyle \frac{\sqrt{{\mathcal{I}}_m^{\left(0\right)}{\mathcal{K}}_m^{\left(0\right)}}}{{\mathcal{I}}_m^{\left(0\right)}+{\mathcal{K}}_m^{\left(0\right)}}}\left|k\left(u_{+}-u_{-}\right)\right|$$

(27)

With a little algebra, this result can be used to show that

$${\displaystyle \frac{{\Lambda }_{+}^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)}}}=i\sqrt{{\displaystyle \frac{{\mathcal{I}}^{\left(0\right)}}{{\mathcal{K}}^{\left(0\right)}}}}\left(k(u_{+}-u_{-})\right),$$

(28)

which will be used in what follows. Eq (27) also illustrates the basic nature of the shear instability with the shortest wavelengths being the fastest growing ones as is typically obtained for the standard inviscid planar shear instability. The information about the geometry is solely contained in the quantities ${\mathcal{I}}_m^{\left(0\right)}$ and ${\mathcal{K}}_m^{\left(0\right)}$. Note that in the limit of large m,k, ${\mathcal{I}}_m^{\left(0\right)}={\mathcal{K}}_m^{\left(0\right)}\simeq 1$ ([22]), using which we further recover the growth rates for the plane parallel shear case with

$${\lambda }_{\Vert }^{\left(0\right)}\simeq -ik{\displaystyle \frac{u_{-}+u_{+}}{2}}+{\displaystyle \frac{1}{2}}\left|k\left(u_{+}-u_{-}\right)\right|.$$

(29)

2.2.2. Rigid Rotation (${\Omega }_{+}={\Omega }_{-}=0$):

In this case, the centrifugal term of the right hand side of the dispersion relation given by (25) vanishes, and so it reduces to

$${\displaystyle \frac{{\mathcal{I}}_m\left(R{\mu }_{+}\right)}{{\mu }_{+}{\Lambda }_{+}^2}}+{\displaystyle \frac{2im{\Omega }_0}{R{\Lambda }_{+}^3{\mu }_{+}^2}}={\displaystyle \frac{2im{\Omega }_0}{R{\Lambda }_{-}^3{\mu }_{-}^2}}-{\displaystyle \frac{{\mathcal{K}}_m\left(R{\mu }_{-}\right)}{{\mu }_{-}{\Lambda }_{-}^2}}.$$

(30)

This equation must usually be solved numerically, but can be solved asymptotically in the regime of weak rotation for which

$${\mu }_{\pm }=k\left[1+{\displaystyle \frac{2{\Omega }_0^2}{{\Lambda }_{\pm }^{{\left(0\right)}^2}}}\right]=k+\mathcal{O}\left({\Omega }_0^2\right).$$

(31)

Using the following ansatz for the growth rates in the above form of the dispersion relation

$$\lambda ={\lambda }^{\left(0\right)}+i{\lambda }^{\left(1\right)},$$

(32)

which implies

$${\Lambda }_{\pm }^{\left(0\right)}={\Lambda }_{\pm }^{\left(0\right)}+i{\lambda }^{\left(1\right)},$$

(33)

where ${\lambda }^{\left(1\right)}$ is assumed to be $\mathcal{O}\left({\Omega }_0\right)$. By linearizing Eq (26) in small ${\Omega }_0$, we can compute the correction to the growth rate ${\lambda }^{\left(1\right)},$ due to rotation, as

$${\lambda }^{\left(1\right)}={\displaystyle \frac{m{\Omega }_0}{kR}}{\displaystyle \frac{\frac{1}{{\Lambda }_{+}^{\left(0\right)3}}-\frac{1}{{\Lambda }_{-}^{\left(0\right)3}}}{\frac{{\mathcal{I}}_m^{\left(0\right)}}{{\Lambda }_{+}^{\left(0\right)3}}+\frac{{\mathcal{K}}_m^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)3}}}}.$$

(34)

Using the expression for the non-rotating growth rates given by Eq 27 in the limit of large k, it can further be shown that

$${\lambda }^{\left(1\right)}=i{\displaystyle \frac{m{\Omega }_0}{kR}}\left(k(u_{+}-u_{-})\right).$$

(35)

We see that, when combined with the ansatz given by Eq (32), the first-order corrections to the non-rotating growth rates have opposite sign as m i.e. positive-m are stabilized by weak rotation whereas negative-m are further destabilized. Furthermore, this ansatz also reduces to the expression given by Eq 29 in the limit of vanishing ${\Omega }_0$. In other words, weak rotation only affects the growth rates to linear order in ${\Omega }_0$ for large values of k.

2.2.3. Differential Rotation (${\Omega }_{+}\ne {\Omega }_{-}$):

For the differentially rotating case, the full dispersion relation (Eq (25)) needs to be solved numerically. It can, however, also be solved analytically when both rotation and differential rotation are weak. As before, we assume the ansatz given by Eq 32, with ${\lambda }^{\left(1\right)}$ now assumed to be $\mathcal{O}\left({\Omega }_0,{\Omega }_{\pm }\right)$. This implies

$$\lambda ={\lambda }^{\left(0\right)}+im{\Omega }_{\pm }+i{\lambda }^{\left(1\right)}.$$

(36)

To linear order in ${\Omega }_0$ and ${\Omega }_{\pm }$, we also have ${\mu }_{\pm }\simeq k$; so, linearizing the full dispersion relation (Eq 25) in the limit of small ${\Omega }_0,{\Omega }_{\pm }$, we obtain

$${\lambda }^{\left(1\right)}={\displaystyle \frac{m}{kR}}{\displaystyle \frac{\left[\frac{{\Omega }_0+{\Omega }_{+}}{{\Lambda }_{+}^{\left(0\right)3}}-\frac{{\Omega }_0+{\Omega }_{-}}{{\Lambda }_{-}^{\left(0\right)3}}-kR\left(\frac{{\Omega }_{+}{\mathcal{I}}_m^{\left(0\right)}}{{\Lambda }_{+}^{\left(0\right)3}}+\frac{{\Omega }_{-}{\mathcal{K}}_m^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)3}}\right)\right]}{\frac{{\mathcal{I}}_m^{\left(0\right)}}{{\Lambda }_{+}^{\left(0\right)3}}+\frac{{\mathcal{K}}_m^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)3}}}}.$$

(37)

Using the expression for ${\displaystyle \frac{{\Lambda }_{+}^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)}}}$ given by Eq (33), one can write the expression for the growth rates in the limit of large m and k as

$$\lambda \simeq -{\displaystyle \frac{ik}{2}}\left(u_{+}+u_{-}\right)-{\displaystyle \frac{im}{2}}\left({\Omega }_{+}+{\Omega }_{-}\right)+{\displaystyle \frac{1}{2}}\left|k\left(u_{+}-u_{-}\right)\right|+{\displaystyle \frac{m}{2}}\left({\Omega }_{+}-{\Omega }_{-}\right)sgn\left(k(u_{+}-u_{-})\right),$$

(38)

the real part of which illustrates the contributions to the growth rate in terms of the vertical and azimuthal components of shear instabilities.

2.2.4. Axisymmetric Mode ($m=0$):

When $m=0,$ Eq (25) reduces to

$${\displaystyle \frac{{\mathcal{I}}_0\left(R{\mu }_{+}\right)}{{\mu }_{+}{\Lambda }_{+}^2}}+{\displaystyle \frac{{\mathcal{K}}_0\left(R{\mu }_{-}\right)}{{\mu }_{-}{\Lambda }_{-}^2}}=R{k}^2\left[{\left({\Omega }_0+{\Omega }_{-}\right)}^2-{\left({\Omega }_0+{\Omega }_{+}\right)}^2\right]$$

(39)

with ${\mathcal{I}}_0={\displaystyle \frac{I_1\left(R{\mu }_{+}\right)}{I_0\left(R{\mu }_{+}\right)}}$ and ${\mathcal{K}}_0={\displaystyle \frac{K_1\left(R{\mu }_{-}\right)}{K_0\left(R{\mu }_{-}\right)}}$, which must again usually be solved numerically. In the special case of pure rotation with no vertical background flow i.e. $u_{+}=u_{-}=0$, for which ${\Lambda }_{\pm }=\lambda$ and ${\mu }_{\pm }=k\sqrt{1+{\displaystyle \frac{4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}{{\lambda }^2}}}$, we recover the dispersion relation for a jet in pure rotational motion inside a rotating medium, as derived by Alterman [10]:

$${\lambda }^2\left[{\mu }_{+}{\displaystyle \frac{I_0\left(R{\mu }_{+}\right)}{I_1\left(R{\mu }_{+}\right)}}+{\mu }_{-}{\displaystyle \frac{K_0\left(R{\mu }_{-}\right)}{K_1\left(R{\mu }_{-}\right)}}\right]=R\left[{\left({\Omega }_0+{\Omega }_{+}\right)}^2-{\left({\Omega }_0+{\Omega }_{-}\right)}^2\right],$$

(40)

which is identical with their Eq 66 (which follows from their Eq 62, but for an error in the sign of the first term on the right hand side) when both the fluids inside and outside the cylinder have uniform density (of 1) and there are no capillary effects.

2.2.5. Asymptotic Limit of Large k:

In the limit $k\gg 1$, $Re\left({\mu }_{\pm }\right)\gg 1$ as well, and we can use series expansion for the modified Bessel functions ([23]) to write

$$\begin{array}{cc}\hfill {\mathcal{I}}_m\left(R{\mu }_{+}\right)& =1-{\displaystyle \frac{1}{2R{\mu }_{+}}}+\mathcal{O}\left({\displaystyle \frac{1}{R^2{\mu }_{+}^2}}\right)\hfill \\ \hfill {\mathcal{K}}_m\left(R{\mu }_{-}\right)& =1+{\displaystyle \frac{1}{2R{\mu }_{-}}}+\mathcal{O}\left({\displaystyle \frac{1}{R^2{\mu }_{-}^2}}\right),\hfill \end{array}$$

(41)

independently of m, as long as $k\gg m$. Neglecting the terms of $\mathcal{O}\left({\mu }_{\pm }^{-2}\right)$ for large values of k, and noting ${\mu }_{\pm }\simeq k$, the full dispersion relation then reduces to

$${\Lambda }_{+}^2+{\Lambda }_{-}^2=Rk\left[{\left({\Omega }_0+{\Omega }_{+}\right)}^2-{\left({\Omega }_0+{\Omega }_{-}\right)}^2\right].$$

(42)

This can be solved analytically for the growing mode of instability as

$$\begin{array}{c}\lambda =-ik{\displaystyle \frac{u_{+}+u_{-}}{2}}-{\displaystyle \frac{im}{2}}\left({\Omega }_{+}+{\Omega }_{-}\right)+\hfill \\ \hfill \sqrt{{\displaystyle \frac{1}{4}}{\left[k\left(u_{+}-u_{-}\right)+m\left({\Omega }_{+}-{\Omega }_{-}\right)\right]}^2+{\displaystyle \frac{1}{2}}kR\left[{\left({\Omega }_0+{\Omega }_{+}\right)}^2-{\left({\Omega }_0+{\Omega }_{-}\right)}^2\right]},\end{array}$$

(43)

and recovers Eq 38 in the limit of weak rotation. The real part of this general expression clearly shows the contribution from the shear (both vertical and azimuthal) and centrifugal terms. The centrifugal term is always destabilizing for an angular momentum profile that decreases outwards (i.e. ${\left({\Omega }_0+{\Omega }_{+}\right)}^2>{\left({\Omega }_0+{\Omega }_{-}\right)}^2)$ which is the condition for the occurrence of the well-known centrifugal instability. The shear term is also destabilizing but its contribution is determined by the relative signs of $k,m,\left(u_{+}-u_{-}\right)\mathrm{and}\left({\Omega }_{+}-{\Omega }_{-}\right)$.

2.3. Numerical Solutions of the Dispersion Relation

In this section, we present few representative numerical solutions of the dispersion relation for given set of values of $m,k,{\Omega }_0,{\Omega }_{+}\mathrm{and}{\Omega }_{-}$, (25) using the numerical root-finding method described by Barrodale and Wilson [24] (see [25] for a discussion on the limitations of this method) combined with the Bessel function approximation routines provided by Zhang and Jin [23]. To do so, we first non-dimensionalise our system with respect to the radius R and the vertical shear $u_{+}(=-u_{-}$) and hence define the non-dimensional quantities as

$$\check{k}=kR,$$

(44)

$$\check{\lambda }={\displaystyle \frac{R\lambda }{u_{+}}},$$

(45)

$${\check{\Omega }}_0={\displaystyle \frac{R{\Omega }_0}{u_{+}}},{\check{\Omega }}_{+}={\displaystyle \frac{R{\Omega }_{+}}{u_{+}}},{\check{\Omega }}_{-}={\displaystyle \frac{R{\Omega }_{-}}{u_{+}}}.$$

(46)

We select only the roots of (25) for the growth rate $\check{\lambda }$ that have positive real part, and when more than one exist, we report only the largest one as a function of the background rotation rate ${\check{\Omega }}_0$ to illustrate the effects of rotation and differential rotation in determining the fastest growing modes.

2.3.1. Comparison with Analytical Results Without Rotation

We first test the results of our solver in absence of rotation (${\check{\Omega }}_{+}={\check{\Omega }}_{-}={\check{\Omega }}_0=0$) and compare the numerical values of maximum growth rates with the analytical expression given by Eq (27). Figure 2 illustrates the exact agreement obtained between the two as a function of the vertical wavenumber $\check{k}$ for two choices of the azimuthal wavenumber $m=0,1$. For higher values of m, both the numerical and analytic solutions become indistinguishable from the solutions for $m=1$ with increasing values of m, but the agreement between the two remains exact. Also, the solutions for the mode $\left(-m,\check{k}\right)$ is identical to that for the mode $\left(m,\check{k}\right)$ due to the symmetry properties of the modified Bessel functions, which is why they are not shown. As expected from pure inviscid shear instabilities, the growth rate of the modes increase with increasing wavenumbers (i.e. increasing $\left|m\right|$ and $|\check{k}|$).

2.3.2. Comparison with Large $\check{k}$ Asymptotic Formula for Differentially Rotating Case

We next test our numerical solutions in the limit of large $\check{k}$ against the asymptotic approximation (Eq 43) for a differentially rotating case. Figure 3 illustrates the agreement between the numerical and asymptotic predictions for $m=0,1$ and 2 with increasing values of $\check{k}$ for an example case in which we have selected ${\check{\Omega }}_{+}=-{\check{\Omega }}_{-}={\check{\Omega }}_0=1$. As is readily seen, for large enough values of $\check{k}$, the numerical and asymptotic results agree well with each other, illustrating that the numerical solver works well in this limiting case as well. The effect of rotation and differential rotation on the growth rate of the modes is discussed in Section 2.4.

2.3.3. Obtaining Solutions for High Rotation Rates or at Low Values of k

In general, finding numerical solutions of the full dispersion relation given by Eq 25 using the root-finding method employed here (which is based on a Muller algorithm) requires the user to supply a sufficiently accurate guess to initialize the iterative solver. To do this, we use one of two possible alternatives (depending on which one is more successful). In some cases we fix m and $\check{k}$ to the desired values, and use the exact non-rotating expression given by Eq 27 as a starting guess and gradually increase the rotation rate and differential rotation towards the target parameter values, each time using the numerical solution obtained for one set of parameter as a guess for the next set of parameters. Alternatively, we can also fix the value of m and of the rotation rates ${\check{\Omega }}_0,{\check{\Omega }}_{+},{\check{\Omega }}_{-}$ and use the asymptotic approximation given by Eq 43 for a sufficiently high value of $\check{k}$ as a starting guess instead, and gradually decrease $\check{k}$ to obtain numerical results at the desired wavenumber. With this continuation technique, we have been able to find solutions of Eq 25 for almost all sets of parameters. We discuss these results in the following section illustrating the nature of the associated (centrifugal and shear) instabilities.

2.4. Effect of Rotation

In this section, we present results on the linear stability of the rotating sheared tube model described in Section 2.1 and clarify the effects of shear, rotation and differential rotation.

2.4.1. Shear with Uniform Rotation

We begin by considering the uniformly rotating case i.e. when ${\check{\Omega }}_{+}={\check{\Omega }}_{-}=0$. Figure 4 shows the growth rate, $Re\left({\check{\lambda }}_{max}\right)$ for two different vertical wavenumbers $\check{k}=1$ and 5, for $m=0,\pm 1,\pm 2$, as a function of the rotation rate ${\Omega }_0$. We also compare the asymptotic predictions for the growth rates for low rotation rates (as given by Eq (32)) with our numerical solutions which shows agreement with the numerical values for ${\Omega }_0\lesssim 0.2$. It clearly shows that the negative m- modes are destabilized while the positive m-modes are stabilized by rotation. This trend, which is predicted by the low ${\check{\Omega }}_0$ asymptotic theory, persists and intensifies at higher rotation rate. We therefore predict that for rapidly rotating systems, negative m-modes could prevail, while for slowly rotating systems, both negative and positive m-modes could be present during the linear growth phase of the instability. We also see, as expected, that higher k modes are more unstable, at least in this inviscid limit.

2.4.2. Shear with Background Rotation and Differential Rotation

We now present illustrative results showing how the growth rate of the $m=\pm 1$,$\check{k}=5$ modes vary with both rotation ${\check{\Omega }}_0$ and differential rotation ${\check{\Omega }}_{+}-{\check{\Omega }}_{-}$ . For simplicity, we assume ${\check{\Omega }}_{+}=-{\check{\Omega }}_{-}=▵\check{\Omega }$ , so positive values of $▵\Omega$ correspond to centrifugally destabilized systems, while negative values of $▵\check{\Omega }$ correspond to centrifugally stabilized systems. As in Figure 4, we also compare these numerical results with the asymptotic predictions for weak rotation (see Eqs. 32-34). Figure 5 demonstrates that the numerical results agree well with the asymptotic predictions for growth rates for ${\check{\Omega }}_0\lesssim 0.1$ with different choices of $▵\check{\Omega }\le {\check{\Omega }}_0$. We also see that the difference between the behavior of the $m=1$ and $m=-1$ modes with increasing $▵\check{\Omega }$ (at a fixed value of ${\check{\Omega }}_0$), that makes the $m=1$ modes more unstable but stabilizes the $m=-1$ mode. This effect can be understood from the asymptotic expression for the growth rates at large $\check{k}$ (given by Eq 43) with the contribution from the centrifugal instabilities given by the term $m({\Omega }_{+}-{\Omega }_{-}$), being of different signs depending on the sign of m.

2.5. Transport of Angular Momentum

Having studied the stability of the system to small perturbations, we now investigate how these perturbations influence the rotation of the cylinder itself. In this section, we go back to dimensional quantities. The specific angular momentum expressed in cylindrical coordinates is given by

$$\mathcal{L}\left(r\right)=r^2\Omega \left(r\right)+r{u}_{\varphi },$$

(47)

where

$$\begin{array}{c}\Omega \left(r\right)=\left\{\begin{array}{cc}{\Omega }_0+{\Omega }_{+},\hfill & r\le R\hfill \\ {\Omega }_0+{\Omega }_{-},\hfill & r>R.\hfill \end{array}\right.\end{array}$$

(48)

In a non-dissipative system, the conservation of angular momentum requires

$${\displaystyle \frac{\partial \mathcal{L}}{\partial t}}+\nabla ·\left(\mathbf{u}\mathcal{L}\right)=0,$$

(49)

which, upon integration over the volume of the cylinder between 0 and $L_z$, gives the net rate of change of angular momentum within as

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={{\displaystyle \frac{\partial }{\partial t}}}_{\varsigma }\mathcal{L}dV={-}_{\varsigma }\nabla ·\left(\mathbf{u}\mathcal{L}\right)dV={-}_{\partial \varsigma }\mathcal{L}\mathbf{u}·\mathbf{n}dS,$$

(50)

where $\zeta$ denotes the volume of the cylinder contained in $r\le R$, $\mathbf{n}$ is the unit vector directed normally outward on the surface element $\partial \zeta$, so ${ʃ}_{\varsigma }f(r,\varphi ,z)dV={\int }_0^{L_z}dz{\int }_0^{2\pi }d\varphi {\int }_0^{R}rf(r,\varphi ,z)dr$ and, given the periodicity in z, ${\int }_{\partial \varsigma }f(R,\varphi ,z)dS=R{\int }_0^{L_z}dz{\int }_0^{2\pi }f(R,\varphi ,z)d\varphi$ (i.e. only the integral over the $r=R$ surface matters).

Using Eq (48), at $r=R$, we then have

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=-\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}R{u}_r(R^2{\Omega }_{+}+R{u}_{\varphi })d\varphi dz=-R^2\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{u}_r{u}_{\varphi }d\varphi dz,$$

(51)

where the first integral vanishes due to mass conservation (Eq (5)). The remaining integral over $\varphi$ and z needs to be taken with care since the quantities $u_r$ and $u_{\varphi }$ are evaluated differently depending on the sign of the displacement . The calculation is presented in Appendix B, and provides the net expression for the angular momentum flux into the tube:

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=-{\displaystyle \frac{\pi L_z{R}^2}{2}}{e}^{2{\lambda }_{R}t}{\left[<{\tilde{u}}_{r_{+}}{\tilde{u}}_{{\varphi }_{+}}>+<{\tilde{u}}_{r_{-}}{\tilde{u}}_{{\varphi }_{-}}>\right]}_{r=R},$$

(52)

where we have defined

$$\begin{array}{c}<{\tilde{u}}_r{\tilde{u}}_{\varphi }>=Re\left[{\tilde{u}}_r\right]Re\left[{\tilde{u}}_{\varphi }\right]+Im\left[{\tilde{u}}_r\right]Im\left[{\tilde{u}}_{\varphi }\right]\\ ={\displaystyle \frac{1}{2}}\left({\tilde{u}}_r^{*}{\tilde{u}}_{\varphi }+{\tilde{u}}_r{\tilde{u}}_{\varphi }^{*}\right),\end{array}$$

(53)

and ${\left[\right]}_{r=R}$ indicates that the quantity is evaluated at $r=R$. Using the solutions ${\tilde{u}}_{r_{\pm }}$ and ${\tilde{u}}_{{\varphi }_{\pm }}$ given by (15), we can write

$$\begin{array}{c}{\left[<{\tilde{u}}_{r_{\pm }}{\tilde{u}}_{{\varphi }_{\pm }}>\right]}_{r=R}={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{4\left({\Omega }_0+{\Omega }_{\pm }\right){\lambda }_R}{d_{\pm }^2}}\left({\displaystyle \frac{m^2}{R^2}}{\left|{\tilde{p}}_{\pm }\right|}^2-{\left|{\displaystyle \frac{\partial {\tilde{p}}_{\pm }}{\partial r}}\right|}^2\right)+\right.\hfill \\ \hfill {\left.{\displaystyle \frac{im}{R}}\left({\tilde{p}}_{\pm }{\displaystyle \frac{\partial {\tilde{p}}_{\pm }^{*}}{\partial r}}-{\tilde{p}}_{\pm }^{*}{\displaystyle \frac{\partial {\tilde{p}}_{\pm }}{\partial r}}\right){\displaystyle \frac{{∥{\Lambda }_{\pm }∥}^2-4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}{d_{\pm }^2}}\right]}_{r=R},\end{array}$$

(54)

where we have defined

$$d_{\pm }^2={∥{\Lambda }_{\pm }^2+4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2∥}^2.$$

(55)

Using the solutions for the pressure given by Eq 12, and the kinematic boundary condition expressed by Eq (19), we can finally express the rate of change of angular momentum in the cylinder explicitly as a function of $\widehat{\eta }$ and other properties of the unstable modes (see Appendix B for detail), as

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=\left\{\begin{array}{c}\begin{array}{cc}{\displaystyle \frac{2{\pi }^2{R}^2}{k}}{e}^{2{\lambda }_{R}t}{\lambda }_R{\widehat{\eta }}^2\left(2{\Omega }_0+{\Omega }_{+}+{\Omega }_{-}\right)& \mathrm{if}m=0,\mathrm{and}\\ {\displaystyle \frac{{\pi }^2{R}^2}{2k}}{e}^{2{\lambda }_{R}t}{\widehat{\eta }}^2\left[{\displaystyle \frac{f_{+}}{{\mathcal{D}}_{+}}}+{\displaystyle \frac{f_{-}}{{\mathcal{D}}_{-}}}\right]& \mathrm{otherwise},\end{array}\hfill \end{array}\right.$$

(56)

where the functions $f_{\pm }$ and ${\mathcal{D}}_{\pm }$ are given by

$$\begin{array}{c}{f}_{\pm }={∥{\Lambda }_{\pm }∥}^2\left[4\left({\Omega }_0+{\Omega }_{\pm }\right){\lambda }_R\left({∥{\mu }_{\pm }∥}^2{∥{\zeta }_{\pm }∥}^2-{\displaystyle \frac{m^2}{R^2}}\right)\pm {\displaystyle \frac{im}{R}}({\mu }_{\pm }{\zeta }_{\pm }-{\mu }_{\pm }^{*}{\zeta }_{\pm }^{*})\left({∥{\Lambda }_{\pm }∥}^2-4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2\right)\right]\\ {\mathcal{D}}_{\pm }={∥{\Lambda }_{\pm }{\mu }_{\pm }{\zeta }_{\pm }\pm {\displaystyle \frac{2im}{R}}\left({\Omega }_0+{\Omega }_{\pm }\right)∥}^2,\end{array}$$

(57)

with ${\zeta }_{+}={\mathcal{I}}_m\left(R{\mu }_{+}\right)$ and ${\zeta }_{-}={\mathcal{K}}_m\left(R{\mu }_{-}\right)$. Using this expression, we can obtain the angular momentum change of the tube for the limiting cases considered already:

1.

Non-rotating case (${\Omega }_{+}={\Omega }_{-}={\Omega }_0=0$) : In the non-rotating limit, ${\mu }_{\pm }=k$ and thus,

$${\displaystyle \frac{f_{\pm }}{{\mathcal{D}}_{\pm }}}={\displaystyle \frac{im}{Rk}}{\displaystyle \frac{\left({\zeta }_{\pm }-{\zeta }_{\pm }^{*}\right)}{{∥{\zeta }_{\pm }∥}^2}}=0,$$

(58)

since ${\zeta }_{\pm }$ are real according to Eq 14. Hence, it is readily seen that there is no net change of angular momentum of the tube in absence of rotation, as expected.

2.

Rigid rotation (${\Omega }_{+}={\Omega }_{-}=0$) : for this scenario,

$$f_{\pm }={∥{\Lambda }_{\pm }∥}^2\left[4{\Omega }_0{\lambda }_R\left({∥{\mu }_{\pm }∥}^2{∥{\zeta }_{\pm }∥}^2-{\displaystyle \frac{m^2}{R^2}}\right)\pm {\displaystyle \frac{im}{R}}\left({\mu }_{\pm }{\zeta }_{\pm }-{\mu }_{\pm }^{*}{\zeta }_{\pm }^{*}\right)\left({∥{\Lambda }_{\pm }∥}^2-4{\Omega }_0^2\right)\right].$$

(59)

If we further assume that ${\Omega }_0$ is small and since ${\mu }_{\pm }\simeq k+\mathcal{O}\left({\Omega }_0^2\right)$, then

$${\displaystyle \frac{f_{\pm }}{{\mathcal{D}}_{\pm }}}\simeq 4{\Omega }_0{\lambda }_R\left(1-{\displaystyle \frac{m^2}{R^2{k}^2}}{\displaystyle \frac{1}{{∥{\zeta }_{\pm }∥}^2}}\right).$$

(60)

Hence, the rate of change of angular momentum can further be simplified to

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\pi }^2{R}^2}{\left|k\right|}}{e}^{2{\lambda }_{R}t}{\widehat{\eta }}^2{\Omega }_0{\lambda }_R\left[2-{\displaystyle \frac{m^2}{R^2{k}^2}}\left({\displaystyle \frac{1}{{\mathcal{I}}_m^{\left(0\right)2}}}+{\displaystyle \frac{1}{{\mathcal{K}}_m^{\left(0\right)2}}}\right)\right]+\mathcal{O}\left({\Omega }_0^2\right),$$

(61)

which can be shown to be positive (see Figure ) for growing modes (${\lambda }_R>0$), implying that there is a net flux of angular momentum into the tube. Note that this expression can also be recovered independently by working with the low ${\Omega }_0,{\Omega }_{\pm }$ asymptotic expressions for $u_{r_{\pm }}$ and $u_{{\varphi }_{\pm }}$ (see Appendix A).

3.

Axisymmetric mode ($m=0$): In this case, $f_{\pm }=4\left({\Omega }_0+{\Omega }_{\pm }\right){\lambda }_R{∥{\mu }_{\pm }∥}^2{∥{\Lambda }_{\pm }∥}^2{∥{\zeta }_{\pm }∥}^2$ and ${\mathcal{D}}_{\pm }={∥{\Lambda }_{\pm }{\mu }_{\pm }{\zeta }_{\pm }∥}^2$, hence

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\pi }^2{R}^2}{\left|k\right|}}{e}^{2{\lambda }_{R}t}{\widehat{\eta }}^2{\lambda }_R\left(2{\Omega }_0+{\Omega }_{+}+{\Omega }_{-}\right),$$

(62)

which could be either positive or negative depending on the signs and relative values of ${\Omega }_0$ and ${\Omega }_{\pm }$. However, the only way to have ${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}<0$ is to have the outer cylinder counter-rotating (i.e. ${\Omega }_0+{\Omega }_{-}$ is sufficiently negative).

4.

Asymptotic limit of large k: in this limit, ${\mu }_{\pm }\simeq k$, ${\mathcal{I}}_m={\mathcal{K}}_m\simeq 1$ which implies

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\pi }^2{R}^2}{\left|k\right|}}{e}^{2{\lambda }_{R}t}{\widehat{\eta }}^2{\lambda }_R\left(k^2-{\displaystyle \frac{m^2}{R^2}}\right)\left[{\displaystyle \frac{{∥{\Lambda }_{+}∥}^2\left({\Omega }_0+{\Omega }_{+}\right)}{{\mathcal{D}}_{+}}}+{\displaystyle \frac{{∥{\Lambda }_{-}∥}^2\left({\Omega }_0+{\Omega }_{-}\right)}{{\mathcal{D}}_{-}}}\right],$$

(63)

which can only be negative if m itself is larger than k or if ${\Omega }_0+{\Omega }_{-}$ is sufficiently negative.

To find the direction of angular momentum transport in general (i.e. not in one of these above limits), we directly compute the expression given by Eq 56 using the numerical solutions for the fastest growing modes obtained in Section 2.3. Figure 6 shows illustrative results obtained for $m=\pm 1,$$\check{k}=1$ and $m=\pm 4,\check{k}=1$ modes for a range of values of the rotation rate ${\check{\Omega }}_0$ and for a range of values of ${\check{\Omega }}_{+}=-{\check{\Omega }}_{-}=\Delta \check{\Omega }$. Red regions in the figure imply net influx of angular momentum into the cylinder and blue regions imply net outflux of angular momentum from of the cylinder. For the $m=1$ case, the numerical results show a distinct region in the $\left({\Omega }_0,\Delta \Omega \right)$ parameter space with positive flux of angular momentum into the cylinder which becomes negative for sufficiently large rotation rate or sufficiently large differential rotation. This predicts a spin-up mechanism for the shear flow for low to moderate rotation rates if $m=1$ is the dominant mode of instability in the flow. For $m=-1$, we predict a similar spin-up through the linear instability across almost the entire parameter space in $\left({\Omega }_0,\Delta \Omega \right)$, except in the case of very weak rotation for a narrow range of values of $\Delta \Omega$ or $\Delta \Omega$ much greater than ${\Omega }_0$. For modes with $\check{k}\gtrsim 2$, the numerical solution always predicts an inward flux of angular momentum into the cylinder for all values of ${\Omega }_0,\Delta \Omega \le 1$, which are hence not plotted in Figure 6. This is also consistent with the asymptotic expression for large k (give by Eq 63) for $\left|m\right|<k$. In the following section, we test these predictions for the spin-up mechanism of our idealized shear flow setup through direct numerical simulations.

3. Numerical Simulations

In order to test our findings regarding the nonlinear evolution of the rotating cylindrical shear flow, we now run some direct numerical simulations (DNSs) of a model that is similar (though not exactly identical) to the one described in Section 2. For this purpose, we use the pseudo-spectral triply periodic PADDI code (Stellmach et al. [26], Traxler et al. [27]) that has been modified to include rotation (Moll and Garaud [5], Sengupta and Garaud [6]) as well as the imposed forcing term. In all that follows, we non-dimensionalize the flow velocity with ${\displaystyle \frac{u_{+}-u_{-}}{2}}$ and the distances with the radius of the cylinder so that $R=1$. With this setup, we define the background flow to be the following top-hat form for the z-component of velocity:

$${\overline{u}}_z\left(r\right)=tanh\left({\displaystyle \frac{1-r}{\delta }}\right)$$

(64)

with $\delta =0.08$ (though we also present a test case with $\delta =0.04$) and $r=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}$ being the distance from the center $(x_c,y_c)=\left({\displaystyle \frac{L_x}{2}},{\displaystyle \frac{L_y}{2}}\right)$ of the tube. The computational domain is of dimensions $L_x\times L_y\times L_z$.

3.1. Non-Dimensional Governing Equation

The non-dimensional Navier-Stokes equations for our problem setup are given by

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·\nabla \mathbf{u}+{\displaystyle \frac{1}{\mathrm{Ro}}}{\mathbf{e}}_{\mathbf{z}}\times \mathbf{u}& =-\nabla p+{\displaystyle \frac{1}{\mathrm{Re}}}{\nabla }^2\mathbf{u}+\mathbf{f},\hfill \\ \hfill \nabla ·\mathbf{u}& =0.\hfill \end{array}$$

(65)

Note that for the DNSs, the effect of viscosity must be retained to ensure stability. We see the usual non-dimensional parameters appear, such as the Reynolds number $\mathrm{Re}={\displaystyle \frac{R\left(u_{+}-u_{-}\right)}{2\nu }}$ (where $\nu$ is the viscosity) and the Rossby number, $\mathrm{Ro}={\displaystyle \frac{u_{+}-u_{-}}{4R{\Omega }_0}}$. A forcing term $\mathbf{f}=f\left(r\right){\mathbf{e}}_{\mathbf{z}}$ is imposed on the system, which drives the desired cylindrical shear flow ${\overline{u}}_z$ (given by Eq 64) provided $f\left(r\right)$ is defined as:

$$f\left(r\right)=-{\displaystyle \frac{1}{\mathrm{Re}}}\left({\displaystyle \frac{d^2{\overline{u}}_z}{d{r}^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d{\overline{u}}_z}{dr}}\right).$$

(66)

To test the predictions of our model, we investigate the sign of the angular momentum transport during the early development of the instability by computing the flux of angular momentum through the cylinder of radius r centered at $(x_c,y_c)$, at each time t, as:

$$F(r,t)=\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{r}^2{u}_r{u}_{\varphi }d\varphi dz.$$

(67)

To do so, we first calculate $u_{r}and{u}_{\varphi }$ from the Cartesian components of the velocity vector $\mathbf{u}=\left(u,v,w\right)$ returned by PADDI using the following transformations:

$$\begin{array}{cc}\hfill u_r& =ucos\varphi +vsin\varphi \hfill \\ \hfill u_{\varphi }& =-usin\varphi +vcos\varphi \hfill \end{array}$$

(68)

and, where at each position $(x,y,z)$,

$$cos\varphi ={\displaystyle \frac{x-x_c}{r}},sin\varphi ={\displaystyle \frac{y-y_c}{r}}.$$

(69)

Using the azimuthal velocity $u_{\varphi }$, we also compute the average angular velocity as a function of the radial distance from the center of the domain at a given time t:

$$\overline{\Omega }(r,t)={\displaystyle \frac{1}{2\pi L_z}}\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{\displaystyle \frac{u_{\varphi }(r,\varphi ,z,t)}{r}}d\varphi dz.$$

(70)

3.2. Typical DNS

We first present results for a single DNS at $\mathrm{Ro}=8$, and fairly large Reynolds numbers, to ensure that the results are not dominated by viscous effects (which were ignored in the analytical computations presented in Section 2).

For our DNSs at $\mathrm{Ro}=8$ (that corresponds to a ${\check{\Omega }}_0=0.0625$), Figure 7 shows isocontours of vertical components of velocities for two runs A2 and A3 (as in Table 1) during the linear growth phase (left panels). Both these runs clearly show the emergence of a preferred mode with low values of the azimuthal wavenumber ($m\sim 1$) and fairly high vertical wavenumber ($k\simeq 5$), irrespective of the domain aspect ratio (the run A1 with the same aspect ratio as A3 also emerged with a similar mode). We now look at the results of our simulations more quantitatively.

3.2.1. Early Phase:

We compare the evolution of the run A3 with the predictions from our linear theory presented in Section 2, with regards to the direction of angular momentum transport. The left panel shows the profiles of the total angular momentum flux $F(r,t)$ as a function of r, at selected times during the growth of the linear instability, showing a clear negative spike at $r=1$, which should result in a positive flux of angular momentum into the cylinder. The right panel of Figure 8 shows the evolution of $\overline{\Omega }$ during the exponential growth phase of the linear instability, clearly demonstrating that the angular velocity within the region $r<1$ grows with time, and thus validates the feedback mechanism of angular momentum predicted from our linear theory (given by Eq 56).

We also calculate the average rate of rotation for the cylinder as

$${\Omega }_{av}\left(t\right)={\displaystyle \frac{2}{R^2}}\underset{0}{\overset{R}{ʃ}}r\overline{\Omega }(r,t)dr,$$

(71)

Figure 9 follows the evolution of ${\Omega }_{av}$ (calculated according to Eq 71) for our DNSs A1 and A3. All these simulations illustrate the initial positive flux of angular momentum into the tube that significantly spins up the region inside $r=R$.

3.2.2. Growth of Centrifugal Instabilities

Figure 8 shows that $F(r,t)$ becomes positive near $r=1$ around $t=9.4$, which halts the growth of the angular velocity for $r<1$. To understand why this happens, we plot the vertically and azimuthally-averaged angular momentum profile, defined by

$$\varphi (r,t)=r^2\left(\overline{\Omega }(r,t)+{\displaystyle \frac{1}{2\mathrm{Ro}}}\right).$$

(72)

Figure 10 shows a radial decrease in the angular momentum profiles around $r=1$, starting at the earliest time snapshot presented ($t\simeq 6.8$). When $\varphi$ is decreasing with increasing r, this angular momentum distribution is well-known to be unstable according to the generalized Rayleigh criteria for centrifugal instabilities ([28]) which requires

$${\displaystyle \frac{1}{r^3}}{\displaystyle \frac{\partial {\varphi }^2}{\partial r}}>0$$

(73)

for any flow to be stable to perturbations. Any centrifugal instability is known to work against its causing mechanism and thus would act to stabilize this negative angular momentum gradient. However, from Figure 9, we see that the tube still continues to spin up (manifested in an increase of ${\Omega }_{av}$ ) until the associated instabilities reach finite amplitudes (at $t\simeq 8.6$) for them to arrest further spin-up of the tube. The competing effects of an unstable angular momentum gradient and associated instabilities control the ultimate outcome of the nonlinear evolution of the spin of the tube, as narrated in the following section.

3.2.3. Nonlinear Evolution: Emergence of Negative Helical Structure

We followed the evolution of one of our low (A1) and high (A3) Reynolds number runs until a quasi-steady state was achieved to check for the final nonlinear outcome for the direction of rotation of the cylinder. At late times, as shown in the Figure 11 for the run A3, the flow through the cylinder (inside $r=1$) does not remain perfectly coherent anymore but gradually evolves into a weakly helical structure.

12 shows the time-evolution of the average rotation rate of the tube, ${\Omega }_{av}\left(t\right)$ (as defined by Eq 71) for both runs. Beyond the early phase of positive growth, the onset of centrifugal instabilities over time appears to make the tube marginally counter-rotating. In order to determine whether the negative rotation rate is intrinsic to the tube, or may be contaminated by the presence of its helical twist, we also propose an alternative way of measuring the average angular velocity of the tube by re-defining the center of the tube at a particular height z as

$$x_c(z,t)=\frac{{ʃ}_0^{L_y}{ʃ}_0^{L_x}x{u}_z(x,y,z,t)\mathcal{H}\left(u_z\right)dxdy}{{ʃ}_0^{L_y}{ʃ}_0^{L_x}{u}_z(x,y,z,t)\mathcal{H}\left(u_z\right)dxdy}$$

(74)

$$y_c(z,t)=\frac{{ʃ}_0^{L_y}{ʃ}_0^{L_x}y{u}_z(x,y,z,t)\mathcal{H}\left(u_z\right)dxdy}{{ʃ}_0^{L_y}{ʃ}_0^{L_x}{u}_z(x,y,z,t)\mathcal{H}\left(u_z\right)dxdy}$$

(75)

where $\mathcal{H}$ is the Heaviside function. We over-plot the results (shown as symbols) in Figure 12 , which demonstrates the helical shape of the flow only influences our results beyond the initial phase of linear growth, when the center (as defined by Eq 74) does not shift much from the center of our domain. In the nonlinear regime, our alternative method gives a more continuous pattern in the evolution of the averaged rotation rate with time although the general conclusion with regards a weakly counter-rotating (i.e. with marginally negative ${\Omega }_{av}$) flow seems to hold till the tube reaches a quasi-steady state.

3.3. DNSs at lower $\mathrm{Ro}$ values

While our DNSs A1-A3 affirm the prediction from our model (described in Section ) for positive feedback of angular momentum into the tube at early times, our preliminary runs (B1-B3, C1-C2) for lower values of $\mathrm{Ro}$ did not show any signs of spin up of the core of the tube (for $r<R$) even for the initial few snaps shots. This is demonstrated in Figure 13 for two set of runs (B1 and B2) $\mathrm{Ro}=2.4$ which both show a decrease in the averaged rotation rate (as defined in Eq 71) early on in the simulations.

We also checked for the radial profiles of the Reynolds stress function $F(r,t)$ (given by Eq 67) and the average angular velocity (given by Eq 70) at early times for one of these runs (B2) as shown in Figure 14,

which indeed demonstrates a negative flux of angular momentum from the region $r<R$ at all times. These runs still emerged with similar modes (low m, fairly high k) as the higher $\mathrm{Ro}$ cases, but clearly do not conform to the expectation from the linear theory that predicts positive angular momentum feedback into the tube for such linearly unstable modes. In following section, we discuss a plausible cause of this discrepancy, to with our periodic domains that seems to problematic at low $\mathrm{Ro}$ values when the results seem to be dependent on the horizontal extent of the domain.

4. Discussions

We have studied an idealized setup of an inviscid vertical shear flow through a cylindrical tube using linear stability analysis and predicted a positive feedback mechanism transporting angular momentum into the tube. Direct numerical simulations for an analogous setup at $\mathrm{Ro}=8$ support our theoretical predictions for the positive feedback of angular momentum at early times, before the onset of centrifugal instabilities that arrest the growth of the spin-up of the tube. However, our preliminary DNSs with moderate ($\mathrm{Ro}=2.4$) to strong ($\mathrm{Ro}=0.08$) rotation seemingly do not agree with our model predictions. In context of studies of vortex dynamics, it has been pointed out in the literature by Pradeep and Hussain [29] that the use of triply periodic domains can often give incorrect results due to artificial zero circulation constraint inherent in periodic simulations. Following the recent work of Sutyrin and Radko [30], we performed additional runs (B3, B4 in Table 1) with half the vertical extent (for both B3 and B4) but double the horizontal extent (B4).

The corresponding time evolution of the average rotation rate of these two runs are compared in Figure 15 which clearly shows that the wider domain run spins up during early times (consistent with our linear theory predictions) unlike the other ones (B1-B3), all of whom showed no such signs of positive feedback of angular momentum with narrower horizontal extent of domains. These computations are substantially more expensive due to the larger resolutions (${256}^2$ in the horizontal direction for the B4,B5,C2 runs) needed to resolve the shear layers at the edge of the tube that occupies only a very small region at the center of the domain. However, from these numerical experiments, it is apparent that particular caution needs to be exercised in the choice of the domain size (compared to the radial extent of the tube) in order to ensure the DNS predictions are not influenced by the periodicity assumptions. The increasing computational cost for lower $\mathrm{Ro}$ and/or higher $\mathrm{Re}$ simulations inhibits us from performing such robustness tests for our most rapidly rotating runs (C1, C2) at $\mathrm{Ro}=0.04$, but they would be important to establish the theory at very low $\mathrm{Ro}$ values, particularly in context of the regime of emergence of LSVs as discovered in Sengupta and Garaud [6] for astrophysical flows.

Appendix A.1. Asymptotic Approximations in the Limit of Weak Rotation

In this appendix, we present approximate solutions to the full rotating problem in the limit of small rotation rates (${\Omega }_0+{\Omega }_{\pm }\ll {\displaystyle \frac{u_{\pm }}{R}}$). We use the general ansatz for all perturbations given by

$$q=q^{\left(0\right)}+i{q}^{\left(1\right)}$$

(A1)

where the (0) superscript represent non-rotating values of the perturbation $q=\{{\tilde{u}}_{r_{\pm }},{\tilde{u}}_{{\varphi }_{\pm }},{\tilde{u}}_{z_{\pm }},{\tilde{p}}_{\pm }\}$, and assume $q^{\left(1\right)}\ll q^{\left(0\right)}$. Thus,

$${\Lambda }_{\pm }={\lambda }^{\left(0\right)}+i{\lambda }^{\left(1\right)}+im{\Omega }_{\pm }+ik{u}_{\pm }={\Lambda }_{\pm }^{\left(0\right)}+i{\Lambda }_{\pm }^{\left(1\right)}$$

(A2)

with ${\Lambda }_{\pm }^{\left(0\right)}={\lambda }^{\left(0\right)}+ik{u}_{\pm }$, ${\Lambda }_{\pm }^{\left(1\right)}={\lambda }^{\left(1\right)}+m{\Omega }_{\pm }$ with ${\Lambda }_{\pm }^{\left(1\right)}\ll {\Lambda }_{\pm }^{\left(0\right)}$.

For the non-rotating problem (i.e. for ${\Omega }_0={\Omega }_{+}=-{\Omega }_{-}=0$), the velocity components (given by Eq 15) are given by

$${\tilde{u}}_{r_{+}}^{\left(0\right)}=-{\displaystyle \frac{{\mathcal{I}}_m^{\left(0\right)}k}{{\Lambda }_{+}^{\left(0\right)}}}{\widehat{p}}_{+}^{\left(0\right)}{I}_m^{\left(0\right)}\left(kr\right)$$

(A3a)

$${\tilde{u}}_{r_{-}}^{\left(0\right)}={\displaystyle \frac{{\mathcal{K}}_m^{\left(0\right)}k}{{\Lambda }_{-}^{\left(0\right)}}}{\widehat{p}}_{-}^{\left(0\right)}{K}_m^{\left(0\right)}\left(kr\right)$$

(A3b)

$${\tilde{u}}_{{\varphi }_{+}}^{\left(0\right)}=-{\displaystyle \frac{im}{r{\Lambda }_{+}^{\left(0\right)}}}{\widehat{p}}_{+}^{\left(0\right)}{I}_m^{\left(0\right)}\left(kr\right)$$

(A3c)

$${\tilde{u}}_{{\varphi }_{-}}^{\left(0\right)}=-{\displaystyle \frac{im}{r{\Lambda }_{-}^{\left(0\right)}}}{\widehat{p}}_{-}^{\left(0\right)}{K}_m^{\left(0\right)}\left(kr\right)$$

(A3d)

where ${\mathcal{I}}_m^{\left(0\right)}$ and ${\mathcal{K}}_m^{\left(0\right)}$ were defined in Eq 14. The pressure continuity condition at $r=R$ implies

$${\widehat{p}}_{+}^{\left(0\right)}{I}_m\left(kR\right)={\widehat{p}}_{-}^{\left(0\right)}{K}_m\left(kR\right).$$

(A4)

Further, imposing the kinematic condition (given by Eq 16) gives us the dispersion relation (Eq 26) and the non-rotating values of growth rates (Eq 27) from which we obtain

$${\displaystyle \frac{{\Lambda }_{+}^{\left(0\right)}}{{\Lambda }_{-}^{\left(0\right)}}}=i\sqrt{{\displaystyle \frac{{\mathcal{I}}_m^{\left(0\right)}}{{\mathcal{K}}_m^{\left(0\right)}}}}sgn\left(k(u_{+}-u_{-})\right).$$

(A5)

Using the ansatz A1 in the set of governing equations (Eq 4b), we obtain the set of governing equation for the corrections to the perturbations in presence of rotation as

$$\begin{array}{cc}\hfill i{\Lambda }_{\pm }^{\left(1\right)}{\tilde{u}}_{r_{\pm }}^{\left(0\right)}+i{\Lambda }_{\pm }^{\left(0\right)}{\tilde{u}}_{r_{\pm }}^{\left(1\right)}-2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{{\varphi }_{\pm }}^{\left(0\right)}& =-{\displaystyle \frac{d{\tilde{p}}_{\pm }^{\left(1\right)}}{dr}}\hfill \end{array}$$

(A6a)

$$\begin{array}{cc}\hfill i{\Lambda }_{\pm }^{\left(1\right)}{\tilde{u}}_{{\varphi }_{\pm }}^{\left(0\right)}+i{\Lambda }_{\pm }^{\left(0\right)}{\tilde{u}}_{{\varphi }_{\pm }}^{\left(0\right)}+2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{r_{\pm }}^{\left(0\right)}& =-{\displaystyle \frac{im{\tilde{p}}^{\left(1\right)}}{r}}\hfill \end{array}$$

(A6b)

$$\begin{array}{cc}\hfill i{\Lambda }_{\pm }^{\left(1\right)}{\tilde{u}}_{z_{\pm }}^{\left(0\right)}+i{\Lambda }_{\pm }^{\left(0\right)}{\tilde{u}}_{z_{\pm }}^{\left(1\right)}& =-ik{\tilde{p}}_{\pm }^{\left(1\right)}\hfill \end{array}$$

(A6c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left(r{\tilde{u}}_{r_{\pm }}^{\left(1\right)}\right)+{\displaystyle \frac{im}{r}}{\tilde{u}}_{{\varphi }_{\pm }}^{\left(1\right)}+ik{\tilde{u}}_{z_{\pm }}^{\left(1\right)}& =0\hfill \end{array}$$

(A6d)

These equations can be combined to obtain a new equation for the pressure correction term as

$$-{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left(r{\displaystyle \frac{d{\tilde{p}}_{\pm }^{\left(1\right)}}{dr}}\right)+\left({\displaystyle \frac{m^2}{r^2}}+k^2\right){\tilde{p}}_{\pm }^{\left(1\right)}+2\left({\Omega }_0+{\Omega }_{\pm }\right)\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left(r{\tilde{u}}_{\varphi }^{\left(0\right)}\right)-{\displaystyle \frac{im}{r}}{\tilde{u}}_r^{\left(0\right)}\right)=0.$$

(A7)

The last term can be shown to vanish using the non-rotating solutions given above, so Eq A7 reduces to

$$r^2{\displaystyle \frac{d^2{\tilde{p}}_{\pm }^{\left(1\right)}}{d{r}^2}}+r{\displaystyle \frac{d{\tilde{p}}_{\pm }^{\left(1\right)}}{dr}}-(m^2+k^2{r}^2){\tilde{p}}_{\pm }^{\left(1\right)}=0$$

(A8)

which is the same equation as for the non-rotating pressure ${\tilde{p}}_{\pm }^{\left(0\right)}$, and hence, without any loss of generality, we can therefore assume ${\tilde{p}}_{\pm }^{\left(1\right)}=0$. By solving Eq A6a, we then find that the corrections to the velocity perturbations in presence of rotation are

$$\begin{array}{cc}\hfill {\tilde{u}}_{r_{\pm }}^{\left(1\right)}& ={\displaystyle \frac{2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{\varphi \pm }^{\left(0\right)}-i{\Lambda }_{\pm }^{\left(1\right)}{\tilde{u}}_{r\pm }^{\left(0\right)}}{i{\Lambda }_{\pm }^{\left(0\right)}}}\hfill \end{array}$$

(A9a)

$$\begin{array}{cc}\hfill {\tilde{u}}_{{\varphi }_{\pm }}^{\left(1\right)}& ={\displaystyle \frac{-2\left({\Omega }_0+{\Omega }_{\pm }\right){\tilde{u}}_{r\pm }^{\left(0\right)}-i{\Lambda }_{\pm }^{\left(1\right)}{\tilde{u}}_{\varphi \pm }^{\left(0\right)}}{i{\Lambda }_{\pm }^{\left(0\right)}}}\hfill \end{array}$$

(A9b)

$$\begin{array}{cc}\hfill {\tilde{u}}_{z_{\pm }}^{\left(1\right)}& =-{\displaystyle \frac{{\Lambda }_{\pm }^{\left(1\right)}}{{\Lambda }_{\pm }^{\left(0\right)}}}{\tilde{u}}_{z_{\pm }}^{\left(0\right)}.\hfill \end{array}$$

(A9c)

Finally, using the kinematic condition at $r=R$, we can obtain an expression for the correction term to the growth rate as

$${\lambda }^{\left(1\right)}={\displaystyle \frac{m}{kR}}{\displaystyle \frac{\left[\frac{{\Omega }_0+{\Omega }_{+}}{{\Lambda }_{+}^{\left(0\right)3}}-\frac{{\Omega }_0+{\Omega }_{-}}{{\Lambda }_{-}^{\left(0\right)3}}-kR\left(\frac{{\Omega }_{+}{\mathcal{I}}_m^{\left(0\right)}\left(kR\right)}{{\Lambda }_{+}^{\left(0\right)3}}+\frac{{\Omega }_{-}{\mathcal{K}}_m^{\left(0\right)}\left(kR\right)}{{\Lambda }_{-}^{\left(0\right)3}}\right)\right]}{\frac{{\mathcal{I}}_m^{\left(0\right)}\left(kR\right)}{{\Lambda }_{+}^{\left(0\right)3}}+\frac{{\mathcal{K}}_m^{\left(0\right)}\left(kR\right)}{{\Lambda }_{-}^{\left(0\right)3}}}},$$

(A10)

which recovers the earlier result (Eq 34) obtained directly from linearizing Eq 25 for weak rotation rates.

Appendix A.1.1. Angular Momentum Transport:

Using the analytical expressions derived above in the limit of small rotation rates, we can arrive at an approximate expression for the angular momentum transport. Using the following ansatz for ${\tilde{u}}_{r_{\pm }}^{\left(1\right)}$ and ${\tilde{u}}_{{\varphi }_{\pm }}^{\left(1\right)}$:

$${\tilde{u}}_{{\varphi }_{\pm }}^{\left(1\right)}=(a_{\pm }+i{b}_{\pm }){\tilde{u}}_{r_{\pm }}^{\left(0\right)},{\tilde{u}}_{r_{\pm }}^{\left(1\right)}=(c_{\pm }+i{d}_{\pm }){\tilde{u}}_{{\varphi }_{\pm }}^{\left(0\right)},$$

(A11)

it can be shown that the rate of change of angular momentum into the cylinder, as derived in Section 2.5, is then given by

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{{\pi }^2}{k}}{R}^2{e}^{2{\lambda }_{R}t}\left[b_{+}{∥{\tilde{u}}_{r_{+}}^{\left(0\right)}∥}^2+b_{-}{∥{\tilde{u}}_{r_{-}}^{\left(0\right)}∥}^2+d_{+}{∥{\tilde{u}}_{{\varphi }_{+}}^{\left(0\right)}∥}^2+d_{-}{∥{\tilde{u}}_{{\varphi }_{-}}^{\left(0\right)}∥}^2\right].$$

(A12)

With some algebra, the required coefficients in the above ansatz are determined as

$$\begin{array}{cc}\hfill b_{\pm }& ={\displaystyle \frac{\left|k\right(u_{+}-u_{-}\left)\right|}{{∥{\Lambda }_{\pm }^{\left(0\right)}∥}^2}}{\displaystyle \frac{\sqrt{{\mathcal{I}}_m^{\left(0\right)}\left(kR\right){\mathcal{K}}_m^{\left(0\right)}\left(kR\right)}}{\left({\mathcal{I}}_m^{\left(0\right)}\left(kR\right)+{\mathcal{K}}_m^{\left(0\right)}\left(kR\right)\right)}}\left({\Omega }_0+{\Omega }_{\pm }\right)\left(2-{\displaystyle \frac{m^2}{R^2{k}^2}}{\displaystyle \frac{1}{f_{\pm }^2}}\right),\mathrm{with}{f}_{+}={\mathcal{I}}_m^{\left(0\right)}\left(kR\right),f_{-}={\mathcal{K}}_m^{\left(0\right)}\left(kR\right)\hfill \\ \hfill d_{\pm }& =-{\displaystyle \frac{\left|k\right(u_{+}-u_{-}\left)\right|}{{∥{\Lambda }_{\pm }^{\left(0\right)}∥}^2}}{\displaystyle \frac{\sqrt{{\mathcal{I}}_m^{\left(0\right)}\left(kR\right){\mathcal{K}}_m^{\left(0\right)}\left(kR\right)}}{\left({\mathcal{I}}_m^{\left(0\right)}\left(kR\right)+{\mathcal{K}}_m^{\left(0\right)}\left(kR\right)\right)}}\left({\Omega }_0+{\Omega }_{\pm }\right).\hfill \end{array}$$

(A13)

Hence, we can express the angular momentum flux expression as

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\pi }^2{R}^2}{k}}{\widehat{\eta }}^2{e}^{2{\lambda }_{R}t}{\lambda }_R^{\left(0\right)}\left[2{\Omega }_0+{\Omega }_{+}+{\Omega }_{-}-{\displaystyle \frac{m^2}{R^2{k}^2}}\left({\displaystyle \frac{{\Omega }_0+{\Omega }_{+}}{{\mathcal{I}}_m^{\left(0\right)}{\left(kR\right)}^2}}+{\displaystyle \frac{{\Omega }_0+{\Omega }_{-}}{{\mathcal{K}}_m^{\left(0\right)}{\left(kR\right)}^2}}\right)\right]$$

(A14)

where we denote the real part of the non-rotating growth rate as ${\lambda }_{\mathrm{R}}^{\left(0\right)}=\left|k(u_{+}-u_{-})\right|{\displaystyle \frac{\sqrt{{\mathcal{I}}_m^{\left(0\right)}{K}_m^{\left(0\right)}}}{\left({\mathcal{I}}_m^{\left(0\right)}\left(kR\right)+{\mathcal{K}}_m^{\left(0\right)}\left(kR\right)\right)}}$. This also recovers the rigidly-rotating expression given by Eq 61 for ${\Omega }_{+}={\Omega }_{-}=0$.

Appendix A.2. Derivation of Expression for Angular Momentum Transport

Axisymmetric case (m = 0)

For the axisymmetric mode, $u_r$ and $u_{\varphi }$ at $r=R$ depend only on z and t, so

$$\underset{0}{\overset{2\pi }{ʃ}}{u}_r{u}_{\varphi }d\varphi =2\pi u_r{u}_{\varphi },$$

(A15)

and hence, from Eq 51, we get

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=-2\pi R^2\underset{0}{\overset{L_z}{ʃ}}{u}_r(R,z,t)u_{\varphi }(R,z,t)dz.$$

(A16)

The difficulty in evaluating this integral becomes apparent when we recall that

$$u_r(R,z,t)=Re\left[{\tilde{u}}_r\left(R\right)exp\left(ikz+\lambda t\right)\right]$$

(A17)

$$u_{\varphi }(R,z,t)=Re\left[{\tilde{u}}_{\varphi }\left(R\right)exp\left(ikz+\lambda t\right)\right],$$

(A18)

where

$${\tilde{u}}_r\left(R\right)=\left\{\begin{array}{c}\begin{array}{cc}{\tilde{u}}_{r_{+}}\left(R\right)& \mathrm{if}\eta >0\\ {\tilde{u}}_{r_{-}}\left(R\right)& \mathrm{if}\eta <0\end{array}\hfill \end{array}\right.$$

(A19)

and similarly for ${\tilde{u}}_{\varphi }$. So, the integral needs to be split into intervals of positive and negative values of $\eta$, respectively (illustrated in the left panel of Figure B.1). Recalling the ansatz for $\eta$ given by Eq 21, we define $z_0,z_1\mathrm{and}{z}_2$ to be three successive nodes of $cos\left(kz+{\lambda }_{I}t\right)$, such as, e.g.

$$k{z}_0+{\lambda }_{I}t=-{\displaystyle \frac{\pi }{2}},k{z}_1+{\lambda }_{I}t={\displaystyle \frac{\pi }{2}},k{z}_2+{\lambda }_{I}t={\displaystyle \frac{3\pi }{2}}.$$

(A20)

Then, we let

$$\begin{array}{cc}\hfill I_{+}& {=}_0^{z_1}{u}_r(R,z,t)u_{\varphi }(R,z,t)dz\hfill \\ & =e_{z_0}^{2{\lambda }_{R}t}{z}_1\left[\mathcal{R}e\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]\mathcal{R}e\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]{cos}^2(kz+{\lambda }_{I}t)+\mathcal{I}m\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]\mathcal{I}m\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]{sin}^2(kz+{\lambda }_{I}t)\right]dz\hfill \\ & -e_{z_0}^{2{\lambda }_{R}t}{z}_1{\displaystyle \frac{1}{2}}\left[Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]Im\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]+Im\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]Re\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]\right]sin\left(2\left[kz+{\lambda }_{I}t\right]\right)dz\hfill \\ & =e^{2{\lambda }_{R}t}{\displaystyle \frac{\pi }{2k}}\left[Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]Re\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]+Im\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]Im\left[{\tilde{u}}_{{\varphi }_{+}}\left(R\right)\right]\right],\hfill \end{array}$$

where the integral containing $sin\left(2kz+2{\lambda }_{I}t\right)$ vanishes. Likewise, we can evaluate

$$\begin{array}{cc}\hfill I_{-}& {=}_{z_1}^{z_2}{u}_r(R,z,t)u_{\varphi }(R,z,t)dz\hfill \\ & =e^{2{\lambda }_{R}t}{\displaystyle \frac{\pi }{2k}}\left[Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]Re\left[{\tilde{u}}_{{\varphi }_{-}}\left(R\right)\right]+Im\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]Im\left[{\tilde{u}}_{{\varphi }_{-}}\left(R\right)\right]\right].\hfill \end{array}$$

Finally, we find that the rate of change of the total angular momentum of the cylinder is

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}& =-2\pi R^2{e}^{2{\lambda }_{R}t}n\left[I_{+}+I_{-}\right]\hfill \\ & =-{\displaystyle \frac{\pi L_z{R}^2}{2}}{e}^{2{\lambda }_{R}t}{\left[Re\left[{\tilde{u}}_{r_{+}}\right]Re\left[{\tilde{u}}_{{\varphi }_{+}}\right]+Im\left[{\tilde{u}}_{r_{+}}\right]Im\left[{\tilde{u}}_{{\varphi }_{+}}\right]+Re\left[{\tilde{u}}_{r_{-}}\right]Re\left[{\tilde{u}}_{{\varphi }_{-}}\right]+Im\left[{\tilde{u}}_{r_{-}}\right]Im\left[{\tilde{u}}_{{\varphi }_{-}}\right]\right]}_{r=R}\hfill \end{array}$$

(A21)

$$\begin{array}{cc}& =-{\displaystyle \frac{\pi L_z{R}^2}{4}}{e}^{2{\lambda }_{R}t}{\left[{\tilde{u}}_{r_{+}}^{*}{\tilde{u}}_{{\varphi }_{+}}+{\tilde{u}}_{r_{+}}{\tilde{u}}_{{\varphi }_{+}}^{*}+{\tilde{u}}_{r_{-}}^{*}{\tilde{u}}_{{\varphi }_{-}}+{\tilde{u}}_{r_{-}}{\tilde{u}}_{{\varphi }_{-}}^{*}\right]}_{r=R},\hfill \end{array}$$

(A22)

where we recall that from Eq 6, $n={\displaystyle \frac{k{L}_z}{2\pi }}$. Note that, in the last step, we have made use of the complex identity

$$Re\left[{\tilde{u}}_r\right]Re\left[{\tilde{u}}_{\varphi }\right]+Im\left[{\tilde{u}}_r\right]Im\left[{\tilde{u}}_{\varphi }\right]={\displaystyle \frac{1}{2}}\left[{\tilde{u}}_r^{*}{\tilde{u}}_{\varphi }+{\tilde{u}}_r{\tilde{u}}_{\varphi }^{*}\right].$$

(A23)

From Eq Section 2.1, the axisymmetric solutions for ${\tilde{u}}_{r_{\pm }}$ and ${\tilde{u}}_{{\varphi }_{\pm }}$ are given by

$$\begin{array}{c}{\tilde{u}}_{r_{\pm }}={\displaystyle \frac{-{\Lambda }_{\pm }\frac{\partial {\tilde{p}}_{\pm }}{\partial r}}{{\Lambda }_{\pm }^2+4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}}\end{array}$$

(A24)

$$\begin{array}{c}{\tilde{u}}_{{\varphi }_{\pm }}={\displaystyle \frac{2\left({\Omega }_0+{\Omega }_{\pm }\right)\frac{\partial {\tilde{p}}_{\pm }}{\partial r}}{{\Lambda }_{\pm }^2+4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2}},\end{array}$$

(A25)

which further reduces the expression for angular momentum flux into the cylinder to

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\lambda }_R{\pi }^2{R}^2}{k}}{e}^{2{\lambda }_{R}t}\left[{\displaystyle \frac{{\left|\frac{\partial {\tilde{p}}_{+}}{\partial r}\right|}_{r=R}^2\left({\Omega }_0+{\Omega }_{+}\right)}{{∥{\Lambda }_{+}^2+4{\left({\Omega }_0+{\Omega }_{+}\right)}^2∥}^2}}+{\displaystyle \frac{{\left|\frac{\partial {\tilde{p}}_{-}}{\partial r}\right|}_{r=R}^2\left({\Omega }_0+{\Omega }_{-}\right)}{{∥{\Lambda }_{-}^2+4\left({\Omega }_0+{\Omega }_{-}\right)∥}^2}}\right].$$

(A26)

Thus, if ${\Omega }_0+{\Omega }_{\pm }>0$, this expression is always positive implying positive influx of angular momentum into the cylinder. However, for sufficiently negative ${\Omega }_{-}$, ${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}$ could become negative resulting in the drainage of angular momentum out of the cylinder.

We can also rewrite the above expression in terms of $\widehat{\eta }$ by using Eqs 23 and A24, as

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}={\displaystyle \frac{2{\lambda }_R{\pi }^2{R}^2}{k}}{\widehat{\eta }}^2{e}^{2{\lambda }_{R}t}\left[2{\Omega }_0+{\Omega }_{+}+{\Omega }_{-}\right].$$

(A27)

Figure B.1. Illustration of the shape of the interface perturbed from $r=R$ (blue line) with an amplitude $\eta$ for axisymmetric (left) and non-axisymmetric (right) modes.

Figure B.1. Illustration of the shape of the interface perturbed from $r=R$ (blue line) with an amplitude $\eta$ for axisymmetric (left) and non-axisymmetric (right) modes.

Appendix A.2.1. Non-Axisymmetric Case (m≠0):

Using a method similar to the one described above for the axisymmetric case, we can rewrite the integral of the Reynolds stress term over $\varphi$ as

$${ʃ}_0^{2\pi }{u}_r{u}_{\varphi }d\varphi =\left|m\right|\left[\underset{{\varphi }_0}{\overset{{\varphi }_1}{ʃ}}{u}_r(R,z,t)u_{\varphi }(R,z,t)d\varphi +\underset{{\varphi }_1}{\overset{{\varphi }_2}{ʃ}}{u}_r(R,z,t)u_{\varphi }(R,z,t)d\varphi \right]$$

(A28)

where ${\varphi }_0$, ${\varphi }_1$ and ${\varphi }_2$ are successive zeros (${\varphi }_0<{\varphi }_1<{\varphi }_2$) of $\eta$ (see Eq 21), as illustrated in the right panel of Figure B.1. For $m>0$ , they are for instance given by

$${\varphi }_0=-{\displaystyle \frac{\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}},{\varphi }_1={\displaystyle \frac{\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}},{\varphi }_2={\displaystyle \frac{3\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}},$$

(A29)

while for $m<0$, they are

$${\varphi }_0={\displaystyle \frac{3\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}},{\varphi }_1={\displaystyle \frac{\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}},{\varphi }_2=-{\displaystyle \frac{\pi }{2m}}-{\displaystyle \frac{kz+{\lambda }_{I}t}{m}}.$$

(A30)

For a given mode with wavenumbers $m,k$, we can we can evaluate the integral given by Eq A28 in a manner very similar to $I_{+}$ and $I_{-}$ in the axisymmetric case above, to obtain

$$\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{u}_r{u}_{\varphi }d\varphi dz=-{\displaystyle \frac{\pi L_z}{2}}{e}^{2{\lambda }_{R}t}{\left[<{\tilde{u}}_{r_{+}}{\tilde{u}}_{{\varphi }_{+}}>+<{\tilde{u}}_{r_{-}}{\tilde{u}}_{{\varphi }_{-}}>\right]}_{r=R}.$$

(A31)

This finally yields

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=-R^2\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{u}_r{u}_{\varphi }d\varphi =e^{2{\lambda }_{R}t}{\displaystyle \frac{\pi R^2{L}_z}{4}}{\left[{\tilde{u}}_{r_{+}}^{*}{\tilde{u}}_{{\varphi }_{+}}+{\tilde{u}}_{r_{+}}{\tilde{u}}_{{\varphi }_{+}}^{*}+{\tilde{u}}_{r_{-}}^{*}{\tilde{u}}_{{\varphi }_{-}}+{\tilde{u}}_{r_{-}}{\tilde{u}}_{{\varphi }_{-}}^{*}\right]}_{r=R}.$$

(A32)

Now, combining Eq (15) and (19), we can write

$$\begin{array}{cc}\hfill {\widehat{p}}_{+}{I}_m\left(r{\mu }_{+}\right)& ={\displaystyle \frac{d_{+}{\Lambda }_{+}\widehat{\eta }}{-{\Lambda }_{+}{\mu }_{+}\mathcal{I}\left(r{\mu }_{+}\right)-\frac{2im}{r}\left({\Omega }_0+{\Omega }_{+}\right)}},\hfill \\ \hfill {\widehat{p}}_{-}{K}_m\left(r{\mu }_{-}\right)& ={\displaystyle \frac{d_{-}{\Lambda }_{-}\widehat{\eta }}{-{\Lambda }_{+}{\mu }_{+}\mathcal{K}\left(r{\mu }_{-}\right)-\frac{2im}{r}\left({\Omega }_0+{\Omega }_{-}\right)}}.\hfill \end{array}$$

(A33)

with $d_{\pm }$ given by Eq 55. Inserting these expressions for the pressures in Eq (15), with some algebra, the net flux of angular momentum for the cylinder can be expressed in the form

$${\displaystyle \frac{\partial {\mathcal{L}}_{cyl}}{\partial t}}=e^{2{\lambda }_{R}t}{\displaystyle \frac{\pi R^2{L}_z}{4}}{\widehat{\eta }}^2{\left[{\displaystyle \frac{f_{+}}{{\mathcal{D}}_{+}}}+{\displaystyle \frac{f_{-}}{{\mathcal{D}}_{-}}}\right]}_{r=R},$$

(A34)

where the functions $f_{\pm }$ and ${\mathcal{D}}_{\pm }$ are given by

$$\begin{array}{c}{f}_{\pm }={∥{\Lambda }_{\pm }∥}^2\left[4\left({\Omega }_0+{\Omega }_{\pm }\right){\lambda }_R\left({∥{\mu }_{\pm }∥}^2{∥{\zeta }_{\pm }∥}^2-{\displaystyle \frac{m^2}{R^2}}\right)\pm {\displaystyle \frac{im}{R}}({\mu }_{\pm }{\zeta }_{\pm }-{\mu }_{\pm }^{*}{\zeta }_{\pm }^{*})\left({∥{\Lambda }_{\pm }∥}^2-4{\left({\Omega }_0+{\Omega }_{\pm }\right)}^2\right)\right]\\ {\mathcal{D}}_{\pm }={∥{\Lambda }_{\pm }{\mu }_{\pm }{\zeta }_{\pm }\pm {\displaystyle \frac{2im}{R}}\left({\Omega }_0+{\Omega }_{\pm }\right)∥}^2,\end{array}$$

(A35)

with ${\zeta }_{+}={\mathcal{I}}_m\left(R{\mu }_{+}\right)$ and ${\zeta }_{-}={\mathcal{K}}_m\left(R{\mu }_{-}\right)$.

Appendix A.3. Derivation of Mass Conservation

To show that mass is exactly conserved in the formulation of the model setup described in Section (2.1), we compute the net flux of the flow across the cylinder by evaluating the integral

$$\Phi =R\underset{0}{\overset{L_z}{ʃ}}\underset{0}{\overset{2\pi }{ʃ}}{u}_r(R,\varphi ,z,t)d\varphi dz$$

(A36)

which has to be done separately for axisymmetric and non-axisymmetric cases, in a similar manner as Appendix B.

Axisymmetric case (m = 0)

$$\begin{array}{cc}\hfill \Phi & =2\pi Rn{e}^{{\lambda }_{\mathrm{R}}t}\left(\underset{z_0}{\overset{z_1}{ʃ}}{u}_{r_{+}}(R,z,t)dz+{ʃ}_{z_1}^{z_2}{u}_{r_{-}}(R,z,t)dz\right)\hfill \\ & =k{L}_z{e}^{{\lambda }_{\mathrm{R}}t}{ʃ}_{z_0}^{z_1}\left(Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]cos\left(kz+{\lambda }_{I}t\right)-Im\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]sin\left(kz+{\lambda }_{I}t\right)\right)dz\hfill \\ & +k{L}_z{e}^{{\lambda }_{\mathrm{R}}t}{ʃ}_{z_1}^{z_2}\left(Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]cos\left(kz+{\lambda }_{I}t\right)-Im\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]sin\left(kz+{\lambda }_{I}t\right)\right)dz\hfill \\ & =\pi L_z{e}^{{\lambda }_{\mathrm{R}}t}\left(Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]-Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]\right)\hfill \\ & =0,\hfill \end{array}$$

(A37)

by using the identity given by Eq (23), where $n={\displaystyle \frac{k{L}_z}{2\pi }}$.

Appendix Non-axisymmetric case (m≠0):

As in Appendix B, we can rewrite the integral over $\varphi$ in Eq (A36) as

$$\begin{array}{cc}\hfill {\int }_0^{2\pi }{u}_{r}d\varphi & =\left|m\right|\left[{\int }_{{\varphi }_0}^{{\varphi }_1}{u}_{r_{+}}(R,z,t)d\varphi +{\int }_{{\varphi }_1}^{{\varphi }_2}{u}_{r_{-}}(R,z,t)d\varphi \right]\hfill \\ & =\left|m\right|{ʃ}_{{\varphi }_0}^{{\varphi }_1}\left(Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]cos\left(kz+{\lambda }_{I}t\right)-Im\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]sin\left(kz+{\lambda }_{I}t\right)\right)dz\hfill \\ & +\left|m\right|{ʃ}_{{\varphi }_1}^{{\varphi }_2}\left(Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]cos\left(kz+{\lambda }_{I}t\right)-Im\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]sin\left(kz+{\lambda }_{I}t\right)\right)\hfill \\ & =2sgn\left(m\right)\left(Re\left[{\tilde{u}}_{r_{+}}\left(R\right)\right]-Re\left[{\tilde{u}}_{r_{-}}\left(R\right)\right]\right)\hfill \\ & =0,\hfill \end{array}$$

(A38)

due to Eq (23), thereby implying

$$\Phi {=}_0^{L_z}{0}^{2\pi }{u}_{r}d\varphi =0.$$

(A39)

Thus, both for axisymmetric and non-axisymmetric modes, the net mass flux across the cylinder vanishes identically, guaranteeing mass conservation.