Instabilities in Cylindrical Geometry Using the Minimalist Approach: Formalism and Rotational Instabilities

1. Introduction

Cylindrical flows play a foundational role in astrophysics. Their geometry arises naturally from rotation and provides an effective framework for describing matter and magnetic field behavior in extreme environments. Astrophysical jets offer a clear example. Jets from active galactic nuclei, microquasars, gamma-ray bursts, and protostars maintain narrow, elongated shapes over vast distances. Modeling them as cylindrical plasma beams captures their essential dynamics: rotation winds magnetic fields into helical structures that confine and guide the jet. Accretion disks form a second major domain where cylindrical geometry is indispensable. Around compact objects and young stars, matter spirals inward under gravity while maintaining nearly Keplerian rotation. Cylindrical coordinates highlight the dominant azimuthal motion, strong shear between neighboring fluid layers, which drives turbulence and magnetic field amplification.

Although jets and disks operate on much larger scales, their cylindrical geometry is mirrored in laboratory plasma devices based on pinch-type configurations, rotating plasma columns, liquid-metal experiments, and laser-driven high-energy-density plasmas. These systems provide controlled analogues that help test and validate theoretical models used in astrophysics. Additional insight comes from other areas of fluid dynamics where rotating and stratified flows have been studied extensively, as in geophysical fluid dynamics. Despite operating in very different parameter regimes, the mathematical structure of the governing equations is closely related, and techniques developed in those contexts remain directly applicable to cylindrical astrophysical flows.

Stability is central to fluid dynamics because it determines whether a flow configuration can persist and how it evolves. Jets and disks exist in environments dominated by rotation, shear, magnetic fields, and stratification, making them sensitive to small disturbances that can bend, kink, fragment, or drive the flow into turbulence. Stability analysis predicts which outcomes are possible and under what conditions. Stability also governs efficiency and lifetime. For example a stable jet remains collimated and transports energy over large distances, whereas an unstable jet disrupts and deposits energy locally. In this sense, stability analysis is not merely a mathematical exercise but a tool for interpreting the structure, energetics, and evolution of astrophysical flows.

Linear stability analysis is important because every instability begins as an infinitesimal perturbation. By examining the linearized equations of motion, one can determine whether such perturbations grow or decay, identify their characteristic spatial structure, and compute their growth rates. This provides a clear and predictive framework for understanding the onset of instability before nonlinear effects become dominant. Linear theory also reveals which mechanisms drive instability and it often yields results that are universal across a wide range of physical conditions. Moreover, the early stages of instability frequently leave observable signatures, such as oscillations or coherent distortions, that correspond directly to linear modes. Numerical simulations rely on linear theory to validate their behavior and to interpret the transition from linear growth to nonlinear saturation.

Rotational instabilities arise because rotation stores free energy that can be released when perturbations interact with shear, magnetic fields, or stratification. They fall into three main categories: centrifugal instabilities (CFI), magnetorotational instabilities (MRI), and buoyancy-driven instabilities.

CFI is purely hydrodynamic and occurs when the specific angular momentum decreases outward, violating the Rayleigh criterion [1]. Fluid elements displaced outward then lack sufficient angular momentum to remain in orbit. Rayleigh’s early analysis remains the cornerstone of rotational stability theory, and CFI can appear in strongly sheared or transitional regions of accretion disks.

MRI is the magnetohydrodynamic counterpart. A weak magnetic field couples fluid elements at different radii, allowing angular-momentum exchange and destabilizing flows that are hydrodynamically stable. The seminal work of Balbus and Hawley [2], following earlier studies [3,4], established MRI as a dominant mechanism for angular-momentum transport in accretion disks around compact objects and young stars.

A distinct line of recent work by Gourgouliatos and Komissarov [5] has expanded the concept of rotational instabilities by exploring how they operate in recollimation shocks of astrophysical jets. Their studies, which have been generalized including magnetic fields, show that when a jet undergoes recollimation – typically due to pressure mismatches with the surrounding medium – the flow can develop regions locally susceptible to CFI/MRI. Thus rotational instabilities are not limited to pure cylindrical flows but may be triggered in dynamically evolving jet structures, providing a powerful mechanism for energy dissipation and magnetic field evolution.

Buoyancy-driven rotational instabilities arise when a stratified cylindrical flow is subject to an effective gravity composed of real gravitational acceleration plus the outward centrifugal acceleration generated by rotation. Unstable stratification leads to convection, where buoyant fluid rises and denser fluid sinks. In compressible media, stability depends on how density and pressure vary together, quantified by the buoyancy (or Brunt-Väisälä) frequency. In magnetized flows, buoyancy interacts with the magnetic field, producing magnetic-buoyancy (Parker) instability [6], in which magnetic field lines rise under effective gravity while plasma slides downward along them.

Together, centrifugal, magnetorotational, and buoyancy-driven instabilities form the core framework for rotational stability in cylindrical astrophysical flows, revealing how rotation, magnetic fields, and stratification combine to shape the evolution of jets, disks, and plasma columns.

A standard way to obtain the dispersion relation is to start from the governing ideal MHD equations – momentum, induction, continuity, and energy – and linearize them around a chosen one-dimensional equilibrium. Introducing small perturbations and adopting a normal-mode form converts the linearized PDEs into a set of ODEs in the cylindrical radius for the perturbation amplitudes. From this point, several established solution strategies exist, including spectral methods as in the code Legolas [7], the Spectral Web [8], and shooting techniques. In addition, closed-form dispersion relations – rarely obtainable except under strong simplifying assumptions – and local analytic approximations such as WKBJ analysis offer physical insight into the underlying instability mechanisms.

In this paper we follow the minimalist approach [9], in which one eliminates variables systematically until the system reduces to a single first-order ODE, the principal equation. The dispersion relation then follows from the requirement that this equation admit non-trivial solutions satisfying the appropriate boundary conditions.

A particularly effective way to integrate the principal equation, especially when the eigenfunctions are real, is the Schwarzian approach [10], which naturally numbers the eigenmodes providing a clean organizational structure for the spectrum.

The plan of the paper is as follows. In Section 2 we derive the principal equation in cylindrical geometry, outline the relevant boundary conditions, and describe how the Schwarzian approach can be used to find the eigenmodes. We then apply the method to specific examples, obtaining the exact dispersion relation either through numerical integration or through closed-form expressions, and we also derive an approximate dispersion relation using the WKBJ method. Section 3 treats a wall-bounded incompressible rotating flow, while Section 4 examines an incompressible rotating flow with a tangential discontinuity as a representative interface-driven instability. In Section 5 we extend these examples to the compressible regime to assess how compressibility modifies the instability properties. Section 6 offers a concise overview of the results.

2. Linear Analysis, the Principal Equation, and Boundary Conditions

The non-relativistic ideal MHD equations in Lorentz-Heaviside units are:

$$\begin{array}{c}\hfill \left({\displaystyle \frac{𝜕}{𝜕t}}+\mathit{V}·\nabla \right)\rho =-\rho \nabla ·\mathit{V},\end{array}$$

(1)

$$\begin{array}{c}\hfill \left({\displaystyle \frac{𝜕}{𝜕t}}+\mathit{V}·\nabla \right)P=c_s^2\left({\displaystyle \frac{𝜕}{𝜕t}}+\mathit{V}·\nabla \right)\rho ,\end{array}$$

(2)

$$\begin{array}{c}\hfill \rho \left({\displaystyle \frac{𝜕}{𝜕t}}+\mathit{V}·\nabla \right)\mathit{V}=-\nabla \left(P+{\displaystyle \frac{B^2}{2}}\right)+(\mathit{B}·\nabla )\mathit{B}+\rho \mathit{g},\end{array}$$

(3)

$$\begin{array}{c}\hfill \left({\displaystyle \frac{𝜕}{𝜕t}}+\mathit{V}·\nabla \right)\mathit{B}=(\mathit{B}·\nabla )\mathit{V}-\mathit{B}(\nabla ·\mathit{V}),\end{array}$$

(4)

$$\begin{array}{c}\hfill \nabla ·\mathit{B}=0.\end{array}$$

(5)

We are interested in exploring the perturbations of a steady-state that depends only on the cylindrical radius $\varpi$. In particular, we assume that the unperturbed state has density ${\rho }_0\left(\varpi \right)$, pressure $P\left(\varpi \right)$, bulk velocity ${\mathit{V}}_0=V_{0z}\left(\varpi \right)\widehat{z}+V_{0\varphi }\left(\varpi \right)\widehat{\varphi }$, magnetic field ${\mathit{B}}_0=B_{0z}\left(\varpi \right)\widehat{z}+B_{0\varphi }\left(\varpi \right)\widehat{\varphi }$ and that the gravitational field is $\mathit{g}=-g\left(\varpi \right)\widehat{\varpi }$. Introducing the total pressure

$$\begin{array}{c}\hfill \mathrm{\Pi }=P+{\displaystyle \frac{B^2}{2}},\end{array}$$

(6)

the zeroth-order equations are satisfied provided that

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_0^{\prime }=-{\rho }_{0}g+{\displaystyle \frac{{\rho }_0{V}_{0\varphi }^2}{\varpi }}-{\displaystyle \frac{B_{0\varphi }^2}{\varpi }}.\end{array}$$

(7)

Adding perturbations of the form

$$\begin{array}{c}\hfill \mathit{V}={\mathit{V}}_0+\left[V_{1z}\left(\varpi \right)\widehat{z}+V_{1\varpi }\left(\varpi \right)\widehat{\varpi }+V_{1\varphi }\left(\varpi \right)\widehat{\varphi }\right]e^{i(m\varphi +kz-\omega t)},\end{array}$$

(8)

$$\begin{array}{c}\hfill \mathit{B}={\mathit{B}}_0+\left[B_{1z}\left(\varpi \right)\widehat{z}+B_{1\varpi }\left(\varpi \right)\widehat{\varpi }+B_{1\varphi }\left(\varpi \right)\widehat{\varphi }\right]e^{i(m\varphi +kz-\omega t)},\end{array}$$

(9)

$$\begin{array}{c}\hfill \rho ={\rho }_0\left(\varpi \right)+{\rho }_1\left(\varpi \right)e^{i(m\varphi +kz-\omega t)},\end{array}$$

(10)

$$\begin{array}{c}\hfill \mathrm{\Pi }={\mathrm{\Pi }}_0\left(\varpi \right)+{\mathrm{\Pi }}_1\left(\varpi \right)e^{i(m\varphi +kz-\omega t)},\end{array}$$

(11)

with integer m, defining the wavevector on the cylindrical surfaces

$$\begin{array}{c}\hfill {\mathit{k}}_0=k\widehat{z}+{\displaystyle \frac{m}{\varpi }}\widehat{\varphi },\end{array}$$

(12)

and the Doppler-shifted frequency

$$\begin{array}{c}\hfill {\omega }_0=\omega -k{V}_{0z}-{\displaystyle \frac{m}{\varpi }}{V}_{0\varphi },\end{array}$$

(13)

we linearize the equations in the Eulerian approach following a procedure similar to Ref. [11], shown in Appendix A. The same equations can be obtained using the Frieman-Rotenberg formulation [12] following a Lagrangian approach, as shown in Appendix B.

The linearized system reduces to two differential equations for the perturbations of the velocity in the $\widehat{\varpi }$ direction and the total pressure (there are algebraic relations for all the other perturbations, connecting them to $V_{1\varpi }$, ${\mathrm{\Pi }}_1$, and their derivatives).

Instead of $V_{1\varpi }$ it is convenient to use $y_1={\displaystyle \frac{\varpi V_{1\varpi }}{-i{\omega }_0}}$, a quantity connected to the Lagrangian displacement in the $\widehat{\varpi }$ direction.1

In addition, instead of ${\mathrm{\Pi }}_1$ it is convenient to use the perturbation of the total pressure evaluated at the perturbed location of each fluid element $y_2={\mathrm{\Pi }}_1+{\displaystyle \frac{y_1}{\varpi }}{\mathrm{\Pi }}_0^{\prime }$.

The advantage of these replacements is that the new functions $y_1$ and $y_2$ are everywhere continuous, even at locations where the unperturbed state has tangential discontinuities.

2.1. The System for $y_1$, $y_2$

The resulting $2\times 2$ system is

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\varpi }}\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)+\left(\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& -f_{11}\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=0,\end{array}$$

(14)

$$\begin{array}{c}\hfill f_{11}={\displaystyle \frac{k^2{\mathrm{\Pi }}_0^{\prime }}{A}}+{\displaystyle \frac{m(2T+m{\mathrm{\Pi }}_0^{\prime })}{A{\varpi }^2}}+{\displaystyle \frac{{\rho }_0^2{\omega }_0^2(F^2\varpi g-W^2)}{AS\varpi }},f_{12}={\displaystyle \frac{{\tilde{\kappa }}^2\varpi }{A}},\\ \hfill f_{21}=-{\displaystyle \frac{A}{\varpi }}+{\displaystyle \frac{k^2{\mathrm{\Pi }}_0^{\prime 2}}{A\varpi }}+{\displaystyle \frac{{(2T+m{\mathrm{\Pi }}_0^{\prime })}^2}{A{\varpi }^3}}+{\rho }_0{\left({\displaystyle \frac{g}{\varpi }}\right)}^{\prime }-{\displaystyle \frac{{\rho }_0^2{(F^2\varpi g-W^2)}^2}{AS{\varpi }^3}},\end{array}$$

(15)

$$\begin{array}{c}\hfill F={\mathit{k}}_0·{\mathit{B}}_0,T=F{B}_{0\varphi }+{\rho }_0{\omega }_0{V}_{0\varphi },W={\omega }_0{B}_{0\varphi }+F{V}_{0\varphi },\end{array}$$

(16)

$$\begin{array}{c}\hfill A={\rho }_0{\omega }_0^2-F^2,S={\rho }_0(A{c}_s^2+{\omega }_0^2{B}_0^2),{\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-k^2-{\displaystyle \frac{m^2}{{\varpi }^2}}.\end{array}$$

(17)

The connection of the sound velocity (only the unperturbed value is needed) with the unperturbed pressure essentially enters in the expression of ${\mathrm{\Pi }}_0$. For an ideal equation of state with internal energy density ${\displaystyle \frac{P}{\mathrm{\Gamma }-1}}$, the sound velocity is $c_s=\sqrt{{\displaystyle \frac{\mathrm{\Gamma }{P}_0}{{\rho }_0}}}$.

To simplify the expressions we define the quantities F, A, S, $\tilde{\kappa }$, T and W. The first four also appear in the planar case and have the same meaning, see Ref. [13]. Namely, F is connected to magnetic tension, S to compressibility, and $\tilde{\kappa }$ to the local wavenumber in the $\widehat{\varpi }$ direction. (Defining ${\lambda }^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-k^2$, one can replace a planar wave $e^{i(\tilde{\kappa }x+k_{0y}y+k_{0z}z-\omega t)}$ with $H_m^{\left(1\right)}\left(\lambda \varpi \right)e^{i(m\varphi +kz-\omega t)}$ in cylindrical geometry, involving Hankel functions. The relation ${\tilde{\kappa }}^2={\lambda }^2-m^2/{\varpi }^2$ shows that $\tilde{\kappa }$ splits into $\lambda$, representing the wavenumber, and m, the order of the Hankel function. This replacement is exact in the homogeneous case, as shown below, and gives the essential meaning of $\tilde{\kappa }$ and $\lambda$ in all cases.)

The Alfvén velocity also enters naturally in the expressions and has an important role. It is defined as ${\mathit{v}}_{\mathrm{A}}={\displaystyle \frac{{\mathit{B}}_0}{\sqrt{{\rho }_0}}}$ and its component along ${\mathit{k}}_0$ is ${\mathit{v}}_{\mathrm{A}\Vert }={\displaystyle \frac{F/k_0}{\sqrt{{\rho }_0}}}$. Using this velocity we can write $A={\rho }_0({\omega }_0^2-k_0^2{v}_{\mathrm{A}\Vert }^2)$ and $S={\rho }_0^2\left[(c_s^2+v_{\mathrm{A}}^2){\omega }_0^2-k_0^2{c}_s^2{v}_{\mathrm{A}\Vert }^2\right]$.

2.2. The Principal Equation

As explained in Ref. [9], in order to find the dispersion relation of the wave or instability it is sufficient to work with the ratio $Y=y_1/y_2$. The system (14) implies that this ratio satisfies the principal equation

$$\begin{array}{c}\hfill Y^{\prime }=f_{21}{Y}^2-2f_{11}Y-f_{12},\end{array}$$

(18)

and we only need to integrate this single first-order differential equation, requiring Y to be everywhere continuous and to satisfy the correct boundary conditions at the extreme values of $\varpi$.

Knowing Y we can find all the other perturbations from

$$\begin{array}{c}\hfill {\displaystyle \frac{y_2^{\prime }}{y_2}}=-f_{21}Y+f_{11},y_1=Y{y}_2,\end{array}$$

(19)

together with the relations (A21)–(A27) given in Appendix A. We emphasize, however, that these are not needed to obtain the dispersion relation.

Although the principal equation (18) is the only one required, it is sometimes useful – especially for analytical work – to consider equivalent expressions given in Appendix C. These include alternative expressions of the principal equation, as well as reformulations of the system (14) as a single second-order differential equation for either $y_1$ or $y_2$. Among the many possibilities, one can write them as a Sturm-Liouville problem or as a “variable frequency oscillator” equation.

Note that in cases with real oscillating eigenfunctions, Y is real and becomes infinite at some points; the same holds for $1/Y$. This difficulty can be handled numerically by working with an angular variable, as described in Section 2 of Ref. [10]. In particular, with the substitution ${\displaystyle \frac{1}{Y}}=-{\mathcal{Y}}_{3}cot{\displaystyle \frac{\mathrm{\Phi }}{2}}+{\mathcal{Y}}_4$, where ${\mathcal{Y}}_{3,4}\left(\varpi \right)$ are chosen functions, the infinities of $1/Y$ (i.e., the zeros of Y) correspond to $\mathrm{\Phi }=2n\pi$ with integer n. The principal equation then becomes a differential equation for $\mathrm{\Phi }\left(\varpi \right)$

$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\prime }=2{\mathcal{Y}}_3{f}_{12}+\left({\displaystyle \frac{{\mathcal{Y}}_3^{\prime }}{{\mathcal{Y}}_3}}-2{\mathcal{Y}}_4{f}_{12}-2f_{11}\right)sin\mathrm{\Phi }+{\displaystyle \frac{({\mathcal{Y}}_4^2-{\mathcal{Y}}_3^2)f_{12}+2{\mathcal{Y}}_4{f}_{11}-f_{21}-{\mathcal{Y}}_4^{\prime }}{{\mathcal{Y}}_3}}(1-cos\mathrm{\Phi }).\end{array}$$

(20)

In fact working with $\mathrm{\Phi }$ instead of Y may be useful for numbering the eigenvalues and eigenfunctions, as will be shown in specific examples below. The choice of ${\mathcal{Y}}_{3,4}\left(\varpi \right)$ that eliminates $\mathrm{\Phi }$ from the right-hand side of the transformed equation is of particular interest; it corresponds to the Schwarzian $\mathrm{\Phi }$ approach presented in Ref. [10]. In that case, $\mathrm{\Phi }$ is given by the system

$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\prime }=2f_{12}{\mathcal{Y}}_3,{\mathcal{Y}}_3^{\prime }=2f_{12}{\mathcal{Y}}_3{\mathcal{Y}}_4+2f_{11}{\mathcal{Y}}_3,{\mathcal{Y}}_4^{\prime }=f_{12}({\mathcal{Y}}_4^2-{\mathcal{Y}}_3^2)+2f_{11}{\mathcal{Y}}_4-f_{21},\end{array}$$

(21)

(or equivalently by the Schwarzian $\left\{cot{\displaystyle \frac{\mathrm{\Phi }}{2}},\varpi \right\}=2{\kappa }_1^2$), and the solution is

$$\begin{array}{c}\hfill y_1\propto {\displaystyle \frac{1}{\sqrt{{\mathcal{Y}}_3}}}sin{\displaystyle \frac{\mathrm{\Phi }}{2}},y_2={\displaystyle \frac{y_1}{Y}},{\displaystyle \frac{1}{Y}}=-{\mathcal{Y}}_{3}cot{\displaystyle \frac{\mathrm{\Phi }}{2}}+{\mathcal{Y}}_4.\end{array}$$

(22)

2.3. Boundary Conditions

The functions $y_1$ and $y_2$ are everywhere continuous, and the same holds for their ratio Y. Thus, in cases where the unperturbed state has tangential discontinuities, we simply continue the integration of the principal equation when passing from one medium to the next, keeping Y continuous.

If we follow the Schwarzian approach, we integrate the system (21) while keeping $\mathrm{\Phi }$, ${\mathcal{Y}}_3$, and ${\mathcal{Y}}_4$ continuous at the interfaces, which guarantees the continuity of Y and $y_1$.

When the fluid ends at some extreme value of $\varpi$, the value of Y is known there and must be uses as a boundary condition. For example, if the plasma meets a solid wall, Y vanishes because the Lagrangian displacement $y_1$ vanishes. If the plasma is in contact with a medium of constant pressure, Y becomes infinity at the interface because $y_2$ vanishes.

Two nontrivial boundary conditions that require closer examination correspond to the symmetry axis $\varpi \to 0$ and to large radial distances $\varpi \to +\infty$. Note that the non-diverging solution of the principal equation is automatically selected when the problem is integrated numerically using the Schwarzian approach, which is a significant advantage of this method. Here, however, we obtain analytical expressions for the boundary conditions for Y by simplifying the equations.

2.3.1. Boundary Conditions on the Symmetry Axis

Near the axis, the azimuthal components of the velocity and the magnetic field scale linearly with the cylindrical distance ($V_{0\varphi }\propto \varpi$, $B_{0\varphi }\propto \varpi$). The same holds for T, W, and ${\mathrm{\Pi }}_0^{\prime }$. We are interested in examining flows near the axis in the absence of gravity, although the following analysis also applies provided that the gravitational acceleration g vanishes at least linearly with $\varpi$ near the axis. To determine the behavior of $f_{ij}$ and of the solution near the axis, we must distinguish between the cases $m\ne 0$ and $m=0$.

For $m\ne 0$ we get $f_{ij}={\displaystyle \frac{d_{ij}}{\varpi }}$ with constant ${\displaystyle d_{11}=\underset{\varpi \to 0}{lim}{\displaystyle \frac{2mT+m^2{\mathrm{\Pi }}_0^{\prime }}{A\varpi }}}$, ${\displaystyle d_{12}=\underset{\varpi \to 0}{lim}{\displaystyle \frac{-m^2}{A}}}$, ${\displaystyle d_{21}=\underset{\varpi \to 0}{lim}\left[-A+{\displaystyle \frac{{(2T+m{\mathrm{\Pi }}_0^{\prime })}^2}{A{\varpi }^2}}\right]}$. Changing variable to $ln\varpi$, we conclude that the system (14) has power-law solutions $y_{1,2}\propto {\varpi }^{\pm \sqrt{d_{11}^2+d_{12}{d}_{21}}}={\varpi }^{\pm \left|m\right|}$ with ${\displaystyle \frac{y_2}{y_1}}={\displaystyle \frac{d_{11}\pm \left|m\right|}{-d_{12}}}$. The physically acceptable solution is the non-diverging (upper sign), for which

$$\begin{array}{c}\hfill {\displaystyle {\left.{\displaystyle \frac{1}{Y}}\right|}_{\varpi =0}=\underset{\varpi \to 0}{lim}{\displaystyle \frac{\left(\right|m|-1)A+{\rho }_0{(\omega -k{V}_{0z})}^2-k^2{B}_{0z}^2}{m^2}}form\ne 0.}\end{array}$$

(23)

For $m=0$ near the axis we obtain $f_{11}=\varpi b_{11}$, $f_{12}=\varpi b_{12}$, $f_{21}={\displaystyle \frac{b_{21}}{\varpi }}$ with constant $b_{11}={\displaystyle \frac{k^2{\mathrm{\Pi }}_0^{\prime }}{A\varpi }}+{\displaystyle \frac{{\rho }_0^2{\omega }_0^2(F^2\varpi g-W^2)}{SA{\varpi }^2}}$, $b_{12}={\displaystyle \frac{{\tilde{\kappa }}^2}{A}}$, $b_{21}=-A+{\displaystyle \frac{4T^2}{A{\varpi }^2}}$. The system (14) then simplifies to ${\displaystyle \frac{d(y_1/{\varpi }^2)}{dln\varpi }}+(2+b_{11}{\varpi }^2)(y_1/{\varpi }^2)+b_{12}{y}_2=0$, ${\displaystyle \frac{d{y}_2}{dln\varpi }}+b_{21}{\varpi }^2(y_1/{\varpi }^2)-b_{11}{\varpi }^2{y}_2=0$, and its solution is $y_2=$ constant, $y_1/{\varpi }^2=-{\displaystyle \frac{b_{12}}{2}}{y}_2$. Thus the ratio $Y=y_1/y_2$ vanishes near the axis as ${\displaystyle Y\approx -{\varpi }^2\underset{\varpi \to 0}{lim}{\displaystyle \frac{{\tilde{\kappa }}^2}{2A}}}$, or, equivalently,

$$\begin{array}{c}\hfill {\displaystyle {Y|}_{\varpi \approx 0}\approx -{\varpi }^2\underset{\varpi \to 0}{lim}{\displaystyle \frac{1}{2{\rho }_0}}{\displaystyle \frac{{\omega }_0^2-k^2{c}_s^2}{(v_{\mathrm{A}}^2+c_s^2){\omega }_0^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}form=0.}\end{array}$$

(24)

For practical purposes, it is sufficient to impose the boundary condition ${Y|}_{\varpi =0}=0$ for $m=0$.

2.3.2. Boundary Conditions at Infinity

In the limit $\varpi \to +\infty$ the perturbations $y_{1,2}$ must vanish. In addition, if they represent waves, they must be outgoing.

An interesting subcase corresponds to homogeneous static plasma with zero $B_{0\varphi }$ (and no gravity). In that case $f_{11}=0$, $f_{12}={\displaystyle \frac{{\lambda }^2\varpi }{A}}-{\displaystyle \frac{m^2}{A\varpi }}$, $f_{21}=-{\displaystyle \frac{A}{\varpi }}$, with constant $A={\rho }_0({\omega }^2-k^2{v}_{\mathrm{A}}^2)$, and constant ${\lambda }^2={\displaystyle \frac{({\omega }^2-k^2{v}_{\mathrm{A}}^2)({\omega }^2-k^2{c}_s^2)}{(v_{\mathrm{A}}^2+c_s^2){\omega }^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}$. (The quantity $\lambda$ and its connection to $\tilde{\kappa }$ through ${\lambda }^2={\tilde{\kappa }}^2+m^2/{\varpi }^2$ has already been discussed.) The system (14) then simplifies to $\varpi y_1^{\prime }+({\lambda }^2{\varpi }^2-m^2)y_2/A=0$, $\varpi y_2^{\prime }/A-y_1=0$. Combining these, we find that $y_2$ satisfies the Bessel equation ${\varpi }^2{y}_2^{″}+\varpi y_2^{\prime }+({\lambda }^2{\varpi }^2-m^2)y_2=0$, while $Y={\displaystyle \frac{\varpi y_2^{\prime }}{A{y}_2}}$. Using standard properties of the Bessel functions, we choose for $y_2$ the Hankel function whose amplitude decreases with distance. Selecting the root of ${\lambda }^2$ with positive imaginary part, $\lambda =i\sqrt{{\displaystyle \frac{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)(k^2{c}_s^2-{\omega }^2)}{k^2{c}_s^2{v}_{\mathrm{A}}^2-(v_{\mathrm{A}}^2+c_s^2){\omega }^2}}}$ (taking the principal value of the root), the solution is $y_2\propto H_{\left|m\right|}^{\left(1\right)}\left(\lambda \varpi \right)$, which represents an outgoing wave if $\Re \lambda$ and $\Re \omega$ have the same sign. The resulting solution of the principal equation is

$$\begin{array}{c}\hfill Y={\displaystyle \frac{\lambda \varpi {\displaystyle \frac{H_{\left|m\right|+1}^{\left(1\right)}\left(\lambda \varpi \right)}{H_{\left|m\right|}^{\left(1\right)}\left(\lambda \varpi \right)}}-\left|m\right|}{{\rho }_0(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}},\lambda =i\sqrt{{\displaystyle \frac{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)(k^2{c}_s^2-{\omega }^2)}{k^2{c}_s^2{v}_{\mathrm{A}}^2-(v_{\mathrm{A}}^2+c_s^2){\omega }^2}}}.\end{array}$$

(25)

An alternative expression, more convenient when $\lambda$ is purely imaginary (using the relation between Hankel and modified Bessel functions $H_{\alpha }^{\left(1\right)}\left(\lambda \varpi \right)={\displaystyle \frac{2}{i^{\alpha +1}\pi }}{K}_{\alpha }(-i\lambda \varpi )$) is

$$\begin{array}{c}\hfill Y={\displaystyle \frac{-i\lambda \varpi {\displaystyle \frac{K_{\left|m\right|+1}(-i\lambda \varpi )}{K_{\left|m\right|}(-i\lambda \varpi )}}-\left|m\right|}{{\rho }_0(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}},-i\lambda =\sqrt{{\displaystyle \frac{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)(k^2{c}_s^2-{\omega }^2)}{k^2{c}_s^2{v}_{\mathrm{A}}^2-(v_{\mathrm{A}}^2+c_s^2){\omega }^2}}}.\end{array}$$

(26)

2.3.3. Summary of Equations and Boundary Conditions

The equations and boundary conditions required to apply the minimalist approach and obtain the dispersion relation for a given problem in cylindrical geometry are summarized in Table 1.

3. Wall-Bounded Rotational Instabilities

Wall-bounded instabilities refer to configurations in which the flow is confined between two solid cylindrical surfaces (with the inner boundary possibly coinciding with the symmetry axis and the outer boundary extending up to a finite radius or even to $\varpi \to \infty$). As an illustrative example, consider a magnetized incompressible fluid rotating between two cylindrical walls ${\varpi }_1<\varpi <{\varpi }_2$, with angular velocity $\Omega ={\displaystyle \frac{V_{0\varphi }}{\varpi }}$ and an axial magnetic field $B_0\widehat{z}$. The radial profile of $\Omega$ is determined by the gravitational and pressure fields, and force balance in the unperturbed state requires ${\mathrm{\Pi }}_0^{\prime }=-{\rho }_0{g}_{\mathrm{eff}}$, where $g_{\mathrm{eff}}=g-{\Omega }^2\varpi$ incorporates centrifugal acceleration. We focus here on axisymmetric perturbations ($m=0$) in the incompressible limit ($S\to \infty$), for which we find $f_{11}=-{\displaystyle \frac{k^2{\rho }_0{g}_{\mathrm{eff}}}{A}}$, $f_{12}=-{\displaystyle \frac{k^2\varpi }{A}}$, $f_{21}=-{\displaystyle \frac{A}{\varpi }}+{\displaystyle \frac{k^2{\rho }_0^2{g}_{\mathrm{eff}}^2}{A\varpi }}+{\displaystyle \frac{4{\rho }_0^2{\omega }^2{\Omega }^2}{A\varpi }}+{\rho }_0{\left({\displaystyle \frac{g}{\varpi }}\right)}^{\prime }$, $A={\rho }_0({\omega }^2-k^2{v}_{\mathrm{A}}^2)$.

3.1. Approximate Dispersion Relation Using WKBJ

Before presenting exact results, we briefly discuss the approximate dispersion relation obtained from the WKBJ approximation. This provides expressions that capture the correct qualitative behavior of the instabilities and aid physical interpretation. The approximation is based on partially neglecting spatial variations of the unperturbed state, while retaining variations in the quantities that control the relevant physical mechanisms. Density variations may be unimportant, but they must be retained in terms involving the effective gravity if buoyancy effects are to be captured—this is the essence of the Boussinesq approximation. Likewise, variations of the specific angular momentum ${\varpi }^2\Omega$ must be kept when studying CFI and MRI, since these mechanisms depend sensitively on the radial gradient of angular momentum.

Based on the “variable frequency” equation (A45), the local radial wavenumber is $k_{\varpi }={\kappa }_1$, with ${\kappa }_1$ given by equation (A38). Ignoring spatial variations of $v_{\mathrm{A}}$ we obtain $k_{\varpi }^2=-k^2-{\displaystyle \frac{3}{4{\varpi }^2}}+{\displaystyle \frac{{\rho }_0^{\prime }}{2\varpi {\rho }_0}}+{\displaystyle \frac{{\rho }_0^{\prime 2}}{4{\rho }_0^2}}-{\displaystyle \frac{{\rho }_0^{″}}{2{\rho }_0}}-{\displaystyle \frac{{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}}{v_{\mathrm{A}}^2-{\omega }^2/k^2}}+{\displaystyle \frac{4{\Omega }^2{v}_{\mathrm{A}}^2}{{(v_{\mathrm{A}}^2-{\omega }^2/k^2)}^2}}$. Following the Boussinesq approximation, and assuming $\left|k\right|\varpi \gg 1$ we obtain the approximate expression $k_{\varpi }^2=-k^2-{\displaystyle \frac{{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}}{v_{\mathrm{A}}^2-{\omega }^2/k^2}}+{\displaystyle \frac{4{\Omega }^2{v}_{\mathrm{A}}^2}{{(v_{\mathrm{A}}^2-{\omega }^2/k^2)}^2}}$, giving the local $k_{\varpi }$ as a function of $\varpi$. In cases where the radial wavenumber does not change significantly over the interval ${\varpi }_1<\varpi <{\varpi }_2$ – that is, when ${\kappa }_1^{\prime }/{\kappa }_1^2$ is small and $k_{\varpi }$ may be treated as the constant wavenumber of standing waves in the radial direction, equal to an integer multiple of ${\displaystyle \frac{\pi }{{\varpi }_2-{\varpi }_1}}$ – the above relation can be interpreted as an approximate dispersion relation. Under these assumptions it can be written in the form of the following quadratic

$$\begin{array}{c}\hfill {\displaystyle \frac{k_{\varpi }^2+k^2}{k^2}}{({\omega }^2-k^2{v}_{\mathrm{A}}^2)}^2-\left[{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}\right]({\omega }^2-k^2{v}_{\mathrm{A}}^2)-4{\Omega }^2{k}^2{v}_{\mathrm{A}}^2=0.\end{array}$$

(27)

In the absence of magnetic fields and for constant density, we recover the familiar result of radial oscillations with the epicyclic frequency $\omega =\sqrt{{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}}$, or centrifugal instability with growth rate $\Im \omega =\sqrt{-{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}}$, depending on the sign of ${\left({\varpi }^4{\Omega }^2\right)}^{\prime }$, i.e., the distribution of angular momentum in the flow. Similarly, in the unmagnetized and non-rotating case we recover the standard result of radial oscillations with the buoyancy frequency $\omega =\sqrt{-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}}$, or instability when the stratification is unstable. With a nonzero magnetic field and constant density, the necessary condition for instability follows immediately from the approximate relation: ${\left({\Omega }^2\right)}^{\prime }<0$, and unstable modes occur for wavenumbers $k_{\varpi }^2+k^2<{\displaystyle \frac{-\varpi {\left({\Omega }^2\right)}^{\prime }}{v_{\mathrm{A}}^2}}$.

A thin astrophysical disk around a central mass M can be approximated as a cylindrical flow in a gravitational field $g={\displaystyle \frac{GM}{{\varpi }^2}}$ with Keplerian rotation $\Omega =\sqrt{{\displaystyle \frac{GM}{{\varpi }^3}}}$. Hydrodynamic perturbations in such a disk are stable, with oscillation frequency $\omega =\Omega$. In the presence of magnetic field, the solutions of the dispersion relation (27) take the form ${\omega }^2=k^2{v}_{\mathrm{A}}^2+{\displaystyle \frac{{\Omega }^2}{2}}\pm \sqrt{{\displaystyle \frac{{\Omega }^4}{4}}+4k^2{v}_{\mathrm{A}}^2{\Omega }^2}$ and include an unstable branch (${\omega }^2<0$) provided that $k^2<{\displaystyle \frac{3{\Omega }^2}{v_{\mathrm{A}}^2}}$. This is the magnetorotational instability, first analyzed for magnetized Couette flow in Ref. [4] and later identified as the key instability in magnetized astrophysical disks in Ref. [2].

For a general $\Omega \left(\varpi \right)$ profile, the Rayleigh criterion ${\left({\varpi }^4{\Omega }^2\right)}^{\prime }<0$, meaning that the specific angular momentum decreases outward, is associated with CFI, which appears already in the hydrodynamic case. By contrast, the condition ${\left({\Omega }^2\right)}^{\prime }<0$, even when the specific angular momentum increases, is associated with MRI, which requires the presence of a magnetic field; see the related discussion in Section 6.2 of Ref. [14]. Buoyancy instability may also occur. In the incompressible limit it is linked to unstable density stratification with respect to the effective gravity, i.e. when $-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}<0$.

In the general incompressible case, only the smaller root of the quadratic dispersion relation (27) can become negative, i.e. the

$$\begin{array}{c}\hfill {\omega }^2=k^2{v}_{\mathrm{A}}^2+{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3-{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0}{2(1+k_{\varpi }^2/k^2)}}-\sqrt{{\left[{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3-{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0}{2(1+k_{\varpi }^2/k^2)}}\right]}^2+{\displaystyle \frac{4{\Omega }^2{k}^2{v}_{\mathrm{A}}^2}{1+k_{\varpi }^2/k^2}}}.\end{array}$$

(28)

It is straightforward to show that, for nonzero magnetic field, this root becomes negative for wavenumbers satisfying $k_{\varpi }^2+k^2<{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0-\varpi {\left({\Omega }^2\right)}^{\prime }}{v_{\mathrm{A}}^2}}$, provided that the necessary condition for instability $-{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0+\varpi {\left({\Omega }^2\right)}^{\prime }<0\iff {({\rho }_0\Omega )}^{\prime }<{\rho }_0^{\prime }g/\varpi$ is satisfied. For zero magnetic field, the result reduces to

$$\begin{array}{c}\hfill {\omega }^2={\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3-{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0}{2(1+k_{\varpi }^2/k^2)}}-\left|{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3-{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0}{2(1+k_{\varpi }^2/k^2)}}\right|,\end{array}$$

corresponding to instability for all wavenumbers if ${\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}-{\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}}{{\rho }_0}}<0$.)

3.2. Exact Dispersion Relation Using the Principal Equation

The previous results were based on a local analysis. The exact dispersion relation in a rotating flow whose properties vary with radius is, however, a global eigenvalue problem.

To obtain the exact dispersion relation, we must first to specify the unperturbed state in detail. As an illustrative example, consider an incompressible fluid confined between two cylindrical walls ${\varpi }_1<\varpi <{\varpi }_2$, rotating with power-law angular velocity in the absence of gravity and possessing uniform density and magnetic field. We normalize all quantities using ${\varpi }_2$ as the unit length, ${\Omega |}_{{\varpi }_2}$ as the unit frequency, and the fluid density as the unit density. Under these normalizations the flow occupies the interval $1-\Delta \varpi <\varpi <1$, has density $\rho =1$, magnetic field $v_{\mathrm{A}}\widehat{z}$, and rotates with $\Omega ={\varpi }^b$. The three parameters $\Delta \varpi$, b and $v_{\mathrm{A}}$ uniquely define the unperturbed configuration.

To determine the growth rates of the instability, we integrate the principal equation with boundary conditions ${Y|}_{{\varpi }_{1,2}}=0$ on the cylindrical walls. We restrict attention to axisymmetric perturbations, for which the eigenvalues are real (${\omega }^2=-\Im {\omega }^2$, negative for instability) and the eigenfunctions are oscillatory. To avoid singular behavior at points where Y vanishes, we employ the Schwarzian approach (see the end of Section 2.2). Starting the integration from the inner cylinder ${\varpi }_1$ with ${\mathrm{\Phi }|}_{{\varpi }_1}=0$, the eigenvalues are those for which the solution satisfies ${\mathrm{\Phi }|}_{{\varpi }_2}=2n\pi$, with integer n, upon reaching the outer cylinder.

By performing the integration over a grid in the (k–$\Im \omega$) plane, we obtain the quantity ${\displaystyle \frac{{\mathrm{\Phi }|}_{{\varpi }_2}{-\mathrm{\Phi }|}_{{\varpi }_1}}{2\pi }}$ as function of k and $\Im \omega$. The contours where this function equals an integer n form the branches of the dispersion relation. An example is shown in Figure 1 for $\Delta \varpi =0.5$, $b=-3/2$, and $v_{\mathrm{A}}=0.1$. As seen in the top-left panel, four branches exist in this case, corresponding to $n=1,2,3,4$, as indicated by the color bar. Because $\mathrm{\Phi }\left(\varpi \right)$ is monotonically increasing function, the integer n can be used to number the eigenvalues/eigenfunctions. Indeed, as shown in the bottom-left panel of Figure 1, the number of zeros of $y_1$ (and Y) in the interior of the flow equals $n-1$.

The magnetic field, through the stabilizing effect of its tension, strongly modifies both the growth rate of the instability and the range of wavenumbers over which it occurs. The entire surface ${\displaystyle \frac{{\mathrm{\Phi }|}_{{\varpi }_2}{-\mathrm{\Phi }|}_{{\varpi }_1}}{2\pi }}$ is shifted toward lower $\Im \omega$ and lower k, as seen for example in the top-left panel of Figure 2. A further increase to values $v_{\mathrm{A}}\gtrsim 0.45$ (with the other parameters unchanged) completely stabilizes the flow. Conversely, decreasing $v_{\mathrm{A}}$ increases the growth rate and, more importantly, enlarges both the number of branches and the range of unstable wavenumbers, as shown in the top-right panel of Figure 2.

The growth rate also decreases as the width $\Delta \varpi$ becomes smaller, because the variation of $\Omega$ across the flow – which drives the instability – becomes weaker. An example is shown in the bottom-left panel of Figure 2. A further reduction to values $\Delta \varpi \lesssim 0.16$ (with the other parameters unchanged) completely stabilizes the flow.

The variation of $\Omega$ across the flow also depends on the exponent b. As expected, smaller $\left|b\right|$ stabilizes, since it reduces the difference ${\Omega |}_{{\varpi }_1}{-\Omega |}_{{\varpi }_2}$. An example is shown in the bottom-right panel of Figure 2. Increasing to values $\left|b\right|\lesssim 0.2$ (with the other parameters unchanged) again leads to complete stabilization.

All of these behaviors are characteristic of the MRI. (The chosen value of b does not correspond to CFI, and the absence of density stratification excludes buoyancy modes.)

The exact results can be qualitatively understood using the WKBJ analysis of Section 3.1. The range of unstable wavenumbers and the number of branches are controlled by the condition $k_{\varpi }^2+k^2<k_{\max }^2$ with $k_{\max }^2={\displaystyle \frac{{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0-\varpi {\left({\Omega }^2\right)}^{\prime }}{v_{\mathrm{A}}^2}}$, combined with the approximate standing-wave relation $k_{\varpi }={\displaystyle \frac{n\pi }{\Delta \varpi }}$ obtained by approximating the eigenfunction as $sin(k_{\varpi }\varpi +C)$ and enforcing vanishing at both boundaries. This instability condition gives the maximum k of each branch (for the first branch it is roughly the $k_{\max }$ above) and also the maximum mode number n (i.e., the number of branches) through $n<{\displaystyle \frac{k_{\max }\Delta \varpi }{\pi }}$. These estimates are necessarily rough, since $k_{\max }={\displaystyle \frac{\sqrt{{\rho }_0^{\prime }{g}_{\mathrm{eff}}/{\rho }_0-\varpi {\left({\Omega }^2\right)}^{\prime }}}{v_{\mathrm{A}}}}$ depends on radius. In the present case, $k_{\max }={\displaystyle \frac{\sqrt{-2b}}{v_{\mathrm{A}}{\varpi }^{-b}}}$, so it is not obvious which value should be used. One may employ mean values, or obtain upper and lower bounds using the extremal values of the functions entering the WKBJ expressions. The same applies to the dispersion relation itself. The WKBJ result is ${\omega }^2=k^2{v}_{\mathrm{A}}^2+{\displaystyle \frac{(4+2b){\varpi }^{2b}}{2(1+k_{\varpi }^2/k^2)}}-\sqrt{{\left[{\displaystyle \frac{(4+2b){\varpi }^{2b}}{2(1+k_{\varpi }^2/k^2)}}\right]}^2+{\displaystyle \frac{4{\varpi }^{2b}{k}^2{v}_{\mathrm{A}}^2}{1+k_{\varpi }^2/k^2}}}$ with $k_{\varpi }={\displaystyle \frac{n\pi }{\Delta \varpi }}$. The first branch is shown with dashed lines in the top-left panel of Figure 1, using values at the inner and outer radii. As expected, it reproduces only the qualitative behavior of the exact result, which lies between these two extreme WKBJ curves.

4. Interface-Driven Rotational Instabilities

A discontinuity in the unperturbed state may trigger an instability. Here we are interested in examining the stability of two rotating incompressible fluids that are in contact at some cylindrical radius ${\varpi }_{\mathrm{CD}}$ in the absence of gravity. In what follows, subscript 1 refers to the inner region $\varpi <{\varpi }_{\mathrm{CD}}$, which is assumed to extend to the symmetry axis, and subscript 2 to the outer $\varpi >{\varpi }_{\mathrm{CD}}$, which is assumed to extend to infinity. For simplicity, we assume uniform rotation with angular velocities ${\Omega }_1$ and ${\Omega }_2$ in the two regions, constant densities ${\rho }_1$ and ${\rho }_2$, and uniform magnetic fields corresponding to constant Alfvén velocities $v_{\mathrm{A}1}$ and $v_{\mathrm{A}2}$. We restrict attention to axisymmetric perturbations. Under these assumptions, the functions entering the principal equation take the form $f_{11}={\displaystyle \frac{k^2{\Omega }^2\varpi }{{\omega }^2-k^2{v}_{\mathrm{A}}^2}}$, $f_{12}=-{\displaystyle \frac{k^2\varpi }{\rho ({\omega }^2-k^2{v}_{\mathrm{A}}^2)}}$, $f_{21}=\rho {\displaystyle \frac{k^2{\Omega }^4{\varpi }^2+4{\omega }^2{\Omega }^2-{({\omega }^2-k^2{v}_{\mathrm{A}}^2)}^2}{\varpi ({\omega }^2-k^2{v}_{\mathrm{A}}^2)}}$, with the appropriate choice of parameters in each region.

4.1. Exact Dispersion Relation Using the Principal Equation

We can obtain the dispersion relation by numerically integrating the principal equation. A particularly effective method is to use the Schwarzian approach and determine the eigenvalues from the requirement that ${\displaystyle \frac{{\mathrm{\Phi }|}_2{-\mathrm{\Phi }|}_1}{2\pi }}$ be an integer. The subscripts 1 and 2 refer to the extreme values of $\varpi$ in the entire plasma domain; in the present configuration these are ${\varpi }_1=0$ and ${\varpi }_2=+\infty$. This approach has two advantages. First, as in the wall-bounded case, the integer ${\displaystyle \frac{{\mathrm{\Phi }|}_{\varpi =\infty }{-\mathrm{\Phi }|}_{\varpi =0}}{2\pi }}$ can be used to number the eigenvalues. Second, when the flow extends to asymptotic distances – either toward the rotation axis or toward infinity, as happens in the present two-fluid configuration – the correct boundary conditions are automatically satisfied. We start the integration at the interface ($\varpi ={\varpi }_{\mathrm{CD}}$) and integrate both inward toward the axis and outward toward infinity, using arbitrary but identical initial values of $\mathrm{\Phi }$ and ${\mathcal{Y}}_{3,4}$ at ${\varpi }_{\mathrm{CD}}$. The dispersion relation corresponds to the condition that ${\mathrm{\Phi }|}_{\varpi =\infty }{-\mathrm{\Phi }|}_{\varpi =0}$ be an integer multiple of $2\pi$. In practice, the integration need not reach exactly $\varpi =0$; it suffices to stop once $\mathrm{\Phi }$ has approached its asymptotic constant value. Similarly, the outward integration is terminated at a sufficiently large radius where $\mathrm{\Phi }$ has approached its asymptotic constant value.

4.2. Analytical Dispersion Relation

The principal equation for axisymmetric perturbations of a homogeneous incompressible fluid (with constant density $\rho$ and magnetic field corresponding to constant Alfvén velocity $v_{\mathrm{A}}\widehat{z}$) rotating uniformly (with constant angular velocity) is

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\varpi }}\left({\displaystyle \frac{1}{Y}}\right)=-{\displaystyle \frac{k^2\varpi }{\rho ({\omega }^2-k^2{v}_{\mathrm{A}}^2)}}{\left({\displaystyle \frac{1}{Y}}-\rho {\Omega }^2\right)}^2+{\displaystyle \frac{\rho ({\omega }^2-k^2{v}_{\mathrm{A}}^2)}{\varpi }}-{\displaystyle \frac{4\rho {\omega }^2{\Omega }^2/\varpi }{{\omega }^2-k^2{v}_{\mathrm{A}}^2}},\end{array}$$

(29)

and can be solved analytically by employing the equivalent Sturm-Liouville formulation (A43). With ${\displaystyle \frac{1}{Y}}={\displaystyle \frac{-f_{11}-y_1^{\prime }/y_1}{f_{12}}}=\rho {\Omega }^2+{\displaystyle \frac{A{y}_1^{\prime }}{k^2\varpi y_1}}$ the equation for $y_1$ is $y_1^{″}-{\displaystyle \frac{y_1^{\prime }}{\varpi }}-k^2{y}_1+{\displaystyle \frac{4{\Omega }^2/{\omega }^2}{{(1-k^2{v}_{\mathrm{A}}^2/{\omega }^2)}^2}}{k}^2{y}_1=0$, and can be transformed to a modified Bessel equation ${\varpi }^2{\left({\displaystyle \frac{y_1}{\varpi }}\right)}^{″}+\varpi {\left({\displaystyle \frac{y_1}{\varpi }}\right)}^{\prime }-\left(1+{\mathrm{\Lambda }}^2{\varpi }^2\right)\left({\displaystyle \frac{y_1}{\varpi }}\right)=0$, with $\mathrm{\Lambda }=\sqrt{k^2-{\displaystyle \frac{4{\omega }^2{\Omega }^2{k}^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}}$ (${\mathrm{\Lambda }}^2$ equals $-{\kappa }_0^2$, and ${\kappa }_0^2$ is negative in this case – to avoid using the complex ${\kappa }_0$ we define $\mathrm{\Lambda }$ which is real). The general solution is $y_1=C_{+}\varpi K_1\left(\mathrm{\Lambda }\varpi \right)+C_{-}\varpi I_1\left(\mathrm{\Lambda }\varpi \right)$, and therefore the general solution of the principal equation is ${\displaystyle \frac{1}{Y}}=\rho {\Omega }^2+{\displaystyle \frac{A\mathrm{\Lambda }}{k^2\varpi }}{\displaystyle \frac{C_{-}{I}_0\left(\mathrm{\Lambda }\varpi \right)-C_{+}{K}_0\left(\mathrm{\Lambda }\varpi \right)}{C_{-}{I}_1\left(\mathrm{\Lambda }\varpi \right)+C_{+}{K}_1\left(\mathrm{\Lambda }\varpi \right)}}$.

We can use the previous result to obtain the dispersion relation for two uniformly rotating incompressible fluids in contact. In each region we select the modified Bessel function that is regular at the corresponding asymptotic limit: $I_1\left(\mathrm{\Lambda }\varpi \right)$ for the inner region (regular at $\varpi \to 0$) and $K_1\left(\mathrm{\Lambda }\varpi \right)$ for the outer region (regular at $\varpi \to \infty$). Thus the solution is

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{Y}}=\left\{\begin{array}{ccc}{\rho }_1{\Omega }_1^2-{\displaystyle \frac{{\rho }_1(k^2{v}_{\mathrm{A}1}^2-{\omega }^2){\mathrm{\Lambda }}_1}{k^2\varpi }}{\displaystyle \frac{I_0\left({\mathrm{\Lambda }}_1\varpi \right)}{I_1\left({\mathrm{\Lambda }}_1\varpi \right)}},& \varpi <{\varpi }_{\mathrm{CD}},& {\mathrm{\Lambda }}_1=\sqrt{k^2-{\displaystyle \frac{4{\omega }^2{\Omega }_1^2{k}^2}{{(k^2{v}_{\mathrm{A}1}^2-{\omega }^2)}^2}}},\\ {\rho }_2{\Omega }_2^2+{\displaystyle \frac{{\rho }_2(k^2{v}_{\mathrm{A}2}^2-{\omega }^2){\mathrm{\Lambda }}_2}{k^2\varpi }}{\displaystyle \frac{K_0\left({\mathrm{\Lambda }}_2\varpi \right)}{K_1\left({\mathrm{\Lambda }}_2\varpi \right)}},& \varpi >{\varpi }_{\mathrm{CD}},& {\mathrm{\Lambda }}_2=\sqrt{k^2-{\displaystyle \frac{4{\omega }^2{\Omega }_2^2{k}^2}{{(k^2{v}_{\mathrm{A}2}^2-{\omega }^2)}^2}}}.\end{array}\right.\end{array}$$

(30)

The dispersion relation results from the continuity of Y

$$\begin{array}{c}\hfill {\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}-{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}={\rho }_1{\displaystyle \frac{k^2{v}_{\mathrm{A}1}^2-{\omega }^2}{k^2}}{\mathrm{\Lambda }}_1{\displaystyle \frac{I_0\left({\mathrm{\Lambda }}_1{\varpi }_{\mathrm{CD}}\right)}{I_1\left({\mathrm{\Lambda }}_1{\varpi }_{\mathrm{CD}}\right)}}+{\rho }_2{\displaystyle \frac{k^2{v}_{\mathrm{A}2}^2-{\omega }^2}{k^2}}{\mathrm{\Lambda }}_2{\displaystyle \frac{K_0\left({\mathrm{\Lambda }}_2{\varpi }_{\mathrm{CD}}\right)}{K_1\left({\mathrm{\Lambda }}_2{\varpi }_{\mathrm{CD}}\right)}}.\end{array}$$

(31)

4.3. Approximate Dispersion Relation

It is instructive to simplify the above in the limit of large $\left|k\right|{\varpi }_{\mathrm{CD}}$, using the asymptotic approximations $K_0\approx K_1$ and $I_0\approx I_1$. The resulting approximate dispersion relation is

$$\begin{array}{c}\hfill {\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}-{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}={\rho }_1{\displaystyle \frac{k^2{v}_{\mathrm{A}1}^2-{\omega }^2}{\left|k\right|}}\sqrt{1-{\displaystyle \frac{4{\omega }^2{\Omega }_1^2}{{(k^2{v}_{\mathrm{A}1}^2-{\omega }^2)}^2}}}+{\rho }_2{\displaystyle \frac{k^2{v}_{\mathrm{A}2}^2-{\omega }^2}{\left|k\right|}}\sqrt{1-{\displaystyle \frac{4{\omega }^2{\Omega }_2^2}{{(k^2{v}_{\mathrm{A}2}^2-{\omega }^2)}^2}}}.\end{array}$$

(32)

It can be written as a quartic equation in $\omega$ and solved analytically (selecting a posteriori the roots that satisfy the dispersion relation in the form derived above). Clearly, instability (${\omega }^2<0$) is possible only if ${\rho }_1{\Omega }_1^2>{\rho }_2{\Omega }_2^2$. This criterion is the discontinuous analogue of the condition ${({\rho }_0\Omega )}^{\prime }<{\rho }_0^{\prime }g/\varpi$ derived in Section 3.1. It carries the signature of decreasing ${\Omega }^2$, as in MRI, but also of decreasing ${\varpi }^4{\Omega }^2$, as in CFI. In the discontinuous case these two contributions cannot be distinguished. In addition, when the densities differ, the criterion contains a buoyancy-like term: in the absence of gravity, where $g_{\mathrm{eff}}=-{\Omega }^2\varpi$, a positive density jump ${\rho }_1-{\rho }_2$ is destabilizing.

To gain further physical insight, we examine several simple subcases in the following.

For the hydrodynamic case and same densities equation (32) gives $\sqrt{1-{\displaystyle \frac{4{\Omega }_1^2}{{\omega }^2}}}-\sqrt{1-{\displaystyle \frac{4{\Omega }_2^2}{{\omega }^2}}}={\displaystyle \frac{4}{\left|k\right|{\varpi }_{\mathrm{CD}}}}$. The requirement the left-hand side be positive means ${\displaystyle \frac{{\Omega }_2^2}{{\omega }^2}}>{\displaystyle \frac{{\Omega }_1^2}{{\omega }^2}}$ showing the influence of the angular momentum contrast in the centrifugal instability, since the sign of ${\omega }^2$ is the same with the sign of ${\Omega }_2^2-{\Omega }_1^2$.

For the hydrodynamic case with equal angular velocities, the approximate dispersion relation reduces to ${\omega }^2\sqrt{1-{\displaystyle \frac{4{\Omega }^2}{{\omega }^2}}}={\displaystyle \frac{{\rho }_2-{\rho }_1}{{\rho }_2+{\rho }_1}}{\Omega }^2{\varpi }_{\mathrm{CD}}\left|k\right|$, showing explicitly the influence of the density contrast. In this limit the perturbation behaves as a Rayleigh-Taylor mode, with the role of gravity replaced by the centrifugal acceleration ${\Omega }^2\varpi \widehat{\varpi }$.

The stabilizing effect of the magnetic field (through its tension) becomes evident in the simplified case of equal densities and equal magnetic fields. In that case we obtain $\sqrt{1-{\displaystyle \frac{4{\omega }^2{\Omega }_1^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}}-\sqrt{1-{\displaystyle \frac{4{\omega }^2{\Omega }_2^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}}={\displaystyle \frac{4{\omega }^2}{({\omega }^2-k^2{v}_{\mathrm{A}}^2)\left|k\right|{\varpi }_{\mathrm{CD}}}}$. Instability (${\omega }^2<0$) is possible only when ${\Omega }_1^2>{\Omega }_2^2$. After squaring and rearranging the above equation, one obtains a quadratic equation in ${\omega }^2$. Solving this quadratic and selecting the correct root (the one that satisfies the original dispersion relation) we obtain

$$\begin{array}{c}\hfill {\omega }^2={\displaystyle \frac{-{\displaystyle \frac{{\Omega }_1^2-{\Omega }_2^2}{2}}{k}^2{\varpi }_{\mathrm{CD}}^2}{{\displaystyle \frac{2/\mathcal{A}+k^2{\varpi }_{\mathrm{CD}}^2/M_{\mathrm{A}}^2}{1-k^2{\varpi }_{\mathrm{CD}}^2/M_{\mathrm{A}}^4}}+\sqrt{{\left({\displaystyle \frac{2/\mathcal{A}+k^2{\varpi }_{\mathrm{CD}}^2/M_{\mathrm{A}}^2}{1-k^2{\varpi }_{\mathrm{CD}}^2/M_{\mathrm{A}}^4}}\right)}^2+{\displaystyle \frac{k^2{\varpi }_{\mathrm{CD}}^2-4}{1-k^2{\varpi }_{\mathrm{CD}}^2/M_{\mathrm{A}}^4}}}}},\end{array}$$

(33)

$$\begin{array}{c}\hfill \mathcal{A}={\displaystyle \frac{{\Omega }_1^2-{\Omega }_2^2}{{\Omega }_1^2+{\Omega }_2^2}},M_{\mathrm{A}}^2={\displaystyle \frac{{\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}^2-{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}^2}{B_1^2+B_2^2}}.\end{array}$$

(34)

The above expression always gives real ${\omega }^2$: it is negative for $\left|k\right|{\varpi }_{\mathrm{CD}}<M_{\mathrm{A}}^2$, but positive for $\left|k\right|{\varpi }_{\mathrm{CD}}>M_{\mathrm{A}}^2$. Thus the magnetic field not only reduces the growth rate compared to the hydrodynamic case, but for sufficiently small wavelengths its tension becomes strong enough to completely stabilize the perturbation. The unstable regime corresponds to $\left|k\right|<\rho {\displaystyle \frac{{\Omega }_1^2-{\Omega }_2^2}{2B^2}}{\varpi }_{\mathrm{CD}}$.

It is also instructive to derive an approximate dispersion relation by assuming that the perturbation is local, i.e. in the limit of large k and ${\varpi }_{\mathrm{CD}}$. In each region the radial acceleration satisfies Newton’s law, $\ddot{\xi }-{\Omega }^2(\varpi +\xi )=-{\displaystyle \frac{{\mathrm{\Pi }}_0^{\prime }+\delta {\mathrm{\Pi }}^{\prime }}{\rho }}-k^2{v}_{\mathrm{A}}^2\xi$, where the last term represents the magnetic-tension contribution. We assume $\delta {\mathrm{\Pi }}^{\prime }=i\tilde{\kappa }\delta \mathrm{\Pi }$ with radial wavenumber $\tilde{\kappa }=\mp i\left|k\right|$, in analogy with the planar Rayleigh-Taylor instability in the incompressible limit (where $k_x^2+k_y^2+k_z^2=0$). The upper/lower sign corresponds to the inner/outer part, giving exponential decay of the amplitude away from the interface. A more accurate approximation would use the WKBJ form of equation (A44), which also yields radial wavenumber $\mp i\left|k\right|$ in the limit of slow rotation. Likewise, applying WKBJ to equation (A43) gives radial wavenumber $\mp i{\mathrm{\Lambda }}_{1,2}$ from equation (30), and in the slow-rotation limit ${\mathrm{\Lambda }}_{1,2}\approx \left|k\right|$. However, the goal here is not to repeat the full derivation, but to isolate the essential physics and obtain a compact approximate relation.

Substituting into Newton’s law gives, in both regimes, $-{\omega }^2\xi -{\Omega }^2\xi =\mp \left|k\right|{\displaystyle \frac{\delta \mathrm{\Pi }}{\rho }}-k^2{v}_{\mathrm{A}}^2\xi$. From this we obtain $\delta \mathrm{\Pi }$, and hence the Lagrangian perturbation of the total pressure at the displaced position of each fluid element, $\Delta \mathrm{\Pi }=\delta \mathrm{\Pi }+{\mathrm{\Pi }}_0^{\prime }\xi \approx \mp {\displaystyle \frac{k^2{v}_{\mathrm{A}}^2-{\omega }^2}{\left|k\right|}}\rho \xi +\rho {\Omega }^2\varpi \xi$. At the tangential discontinuity, both $\xi$ and $\Delta \mathrm{\Pi }$ must be continuous, giving $-{\displaystyle \frac{k^2{v}_{\mathrm{A}1}^2-{\omega }^2}{\left|k\right|}}{\rho }_1+{\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}={\displaystyle \frac{k^2{v}_{\mathrm{A}2}^2-{\omega }^2}{\left|k\right|}}{\rho }_2+{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}$. This reproduces the approximate dispersion relation (32), with the square-root factors replaced by unity. The resulting expression is

$${\omega }^2=-{\displaystyle \frac{{\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}-{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}}{{\rho }_1+{\rho }_2}}\left|k\right|+{\displaystyle \frac{B_1^2+B_2^2}{{\rho }_1+{\rho }_2}}{k}^2,$$

(35)

which is essentially identical to the planar Rayleigh-Taylor dispersion relation (see § 4 of Ref. [13]), with the gravitational force densities ${\rho }_{1,2}g$ replaced by the centrifugal force densities ${\rho }_{1,2}{\Omega }_{1,2}^2{\varpi }_{\mathrm{CD}}$. This approximate expression predicts that ${\omega }^2<0$ for $0<\left|k\right|<{\displaystyle \frac{M_{\mathrm{A}}^2}{{\varpi }_{\mathrm{CD}}}}$, and that the maximum growth rate is ${\displaystyle \frac{M_{\mathrm{A}}}{2}}\sqrt{{\displaystyle \frac{{\rho }_1{\Omega }_1^2-{\rho }_2{\Omega }_2^2}{{\rho }_1+{\rho }_2}}}$ at $k={\displaystyle \frac{M_{\mathrm{A}}^2}{2{\varpi }_{\mathrm{CD}}}}$.

4.4. Application of WKBJ

The simplest variant of the WKBJ approach is to treat the array of the system (14) as approximately constant. In that case the system admits exponential solutions with a constant ratio $y_1/y_2$. Equivalently, the principal equation has approximately constant solutions. By dropping the $Y^{\prime }$ term, these solutions are the roots of $f_{12}{Y}^2-2f_{11}Y-f_{12}=0$, i.e., the ${\displaystyle \frac{1}{Y}}={\displaystyle \frac{f_{11}\mp \sqrt{f_{11}^2+f_{12}{f}_{21}}}{-f_{12}}}=\rho {\Omega }^2\pm {\displaystyle \frac{\rho (k^2{v}_{\mathrm{A}}^2-{\omega }^2)\mathrm{\Lambda }}{k^2\varpi }}$ with $\mathrm{\Lambda }=\sqrt{k^2-{\displaystyle \frac{4{\omega }^2{\Omega }^2{k}^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}}$. The two signs correspond to ${\displaystyle \frac{y_1^{\prime }}{y_1}}={\displaystyle \frac{y_2^{\prime }}{y_2}}=\mp \mathrm{\Lambda }$, as follows directly from system (14). Applying this in the two regions, the upper sign must be used for $\varpi >{\varpi }_{\mathrm{CD}}$ and the lower sign for $\varpi <{\varpi }_{\mathrm{CD}}$, ensuring exponential decay away from the interface on both sides. The dispersion relation follows from continuity of Y, giving ${\rho }_1{\Omega }_1^2-{\displaystyle \frac{{\rho }_1(k^2{v}_{\mathrm{A}1}^2-{\omega }^2){\mathrm{\Lambda }}_1}{k^2{\varpi }_{\mathrm{CD}}}}={\rho }_2{\Omega }_2^2+{\displaystyle \frac{{\rho }_2(k^2{v}_{\mathrm{A}2}^2-{\omega }^2){\mathrm{\Lambda }}_2}{k^2{\varpi }_{\mathrm{CD}}}}$. This is identical to equation (32). As expected, the WKBJ approach reproduces the large $\left|k\right|{\varpi }_{\mathrm{CD}}$ limit of the exact dispersion relation.

(We note that other variants of WKBJ exist, see Section 4 of Ref. [15]. The one used here is the “zeroth-order WKBJ for the system”.)

Summarizing, the exact dispersion relation for axisymmetric perturbations in the incompressible limit can be obtained either by numerical integration of the principal equation or analytically through equation (31). In the limit of large $k{\varpi }_{\mathrm{CD}}$, the dispersion relation reduces to the approximate form (32), which – when the two sides contain identical fluids – yields the explicit solution (33). Furthermore, in the combined limit of large $k{\varpi }_{\mathrm{CD}}$ and slow rotation, the local (planar-type) approximation leads to equation (35).

Figure 3 shows results for the case of two fluids with the same density and the same magnetic field. This is purely a CFI associated with the decrease of ${\Omega }^2$ across the interface. (Equivalently, one may call it MRI, since in this discontinuous configuration the criteria for CFI and MRI reduce to the same inequality ${\Omega }_1^2>{\Omega }_2^2$.) There is no buoyancy contribution because the density is uniform in the examined case.

The exact region of instability – i.e., the maximum wavenumber $k_{\max }$ below which ${\omega }^2<0$ – is obtained by setting ${\omega }^2=0$ in equation (31). This yields

$$\begin{array}{c}\hfill M_{\mathrm{A}}^2={\displaystyle \frac{k_{\max }{\varpi }_{\mathrm{CD}}}{2}}\left[{\displaystyle \frac{I_0\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}{I_1\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}}+{\displaystyle \frac{K_0\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}{K_1\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}}\right],\end{array}$$

a relation depending only on the parameter $M_{\mathrm{A}}^2$. The resulting marrginal-stability curve is shown in the top-left panel of Figure 4. For relatively weak magnetic fields corresponding to large $M_{\mathrm{A}}^2$ and large $k_{\max }{\varpi }_{\mathrm{CD}}$, the solution agrees with the estimation $k_{\max }{\varpi }_{\mathrm{CD}}=M_{\mathrm{A}}^2$ discussed earlier. For strong magnetic fields this approximation breaks down because the assumption $k_{\max }{\varpi }_{\mathrm{CD}}\gg 1$ no longer holds. The exact result is that $k_{\max }{\varpi }_{\mathrm{CD}}$ approaches zero as $M_{\mathrm{A}}\to 1$, implying that for $M_{\mathrm{A}}<1$ the perturbation is stable at all wavelengths.

In the more general cases where the densities and magnetic fields differ between the two regions, the marginal stability curve also depends on the magnetic field contrast. The maximum wavenumber is given by

$$\begin{array}{c}\hfill {\displaystyle \frac{{\rho }_1{\Omega }_1^2{\varpi }_{\mathrm{CD}}^2-{\rho }_2{\Omega }_2^2{\varpi }_{\mathrm{CD}}^2}{{\rho }_1{v}_{\mathrm{A}1}^2+{\rho }_2{v}_{\mathrm{A}2}^2}}={\displaystyle \frac{B_1^2{k}_{\max }{\varpi }_{\mathrm{CD}}}{B_1^2+B_2^2}}{\displaystyle \frac{I_0\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}{I_1\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}}+{\displaystyle \frac{B_2^2{k}_{\max }{\varpi }_{\mathrm{CD}}}{B_1^2+B_2^2}}{\displaystyle \frac{K_0\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}{K_1\left(k_{\max }{\varpi }_{\mathrm{CD}}\right)}},\end{array}$$

(36)

which depends on two parameters: $M_{\mathrm{A}}^2$, and $B_2^2/B_1^2$. These dependencies are illustrated in the first three panels of Figure 4. At large k the result agrees with the approximate dispersion relation (35), giving $k_{\max }{\varpi }_{\mathrm{CD}}=M_{\mathrm{A}}^2$. At small k, however, the exact solution shows that $k_{\max }{\varpi }_{\mathrm{CD}}$ vanishes at a finite value of $M_{\mathrm{A}}^2$, as displayed in the bottom-right panel of Figure 4.

5. Compressibility Effects

5.1. Approximate Dispersion Relation

The WKBJ approximation gives a compact way to understand how rotational instabilities are modified when the sound velocity is finite (and therefore S is also finite). We focus on axisymmetric disturbances ($m=0$) in rotating flows where the unperturbed magnetic field is purely axial.

The local wavenumber in the radial direction can be approximated by ${\kappa }_1$ from equation (A42), which in the present case becomes ${\kappa }_1^2={\tilde{\kappa }}^2+{\displaystyle \frac{{\tilde{\kappa }}^2{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}-{\displaystyle \frac{4{\tilde{\kappa }}^2{k}^2{v}_{\mathrm{A}}^2{\Omega }^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}-{\displaystyle \frac{{\tilde{\kappa }}^2{g}_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}-{\displaystyle \frac{k^2{g}_{\mathrm{eff}}^2}{S/{\rho }_0^2}}-{\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}\right)}^2+{\displaystyle \frac{f_{12}^{\prime }}{f_{12}}}{\displaystyle \frac{{\omega }^2{g}_{\mathrm{eff}}}{S/{\rho }_0^2}}+{\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}-{\displaystyle \frac{{\omega }^2{g}_{\mathrm{eff}}}{S/{\rho }_0^2}}\right)}^{\prime }$. Following the Boussinesq approximation and assuming $\left|k\right|\varpi \gg 1$, the last three terms become negligible. Substituting $S/{\rho }_0^2=(c_s^2+v_{\mathrm{A}}^2){\omega }^2-k^2{v}_{\mathrm{A}}^2{c}_s^2$ and approximating ${\tilde{\kappa }}^2={\displaystyle \frac{{\omega }^4}{S/{\rho }_0^2}}-k^2\approx -k^2$, which remains valid for large $\left|k\right|$ even with compressibility is included, the approximate dispersion relation becomes

$$\begin{array}{c}\hfill {\displaystyle \frac{k_{\varpi }^2+k^2}{k^2}}=-{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}+{\displaystyle \frac{4k^2{v}_{\mathrm{A}}^2{\Omega }^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}+{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}+{\displaystyle \frac{g_{\mathrm{eff}}^2}{k^2{v}_{\mathrm{A}}^2{c}_s^2-(c_s^2+v_{\mathrm{A}}^2){\omega }^2}}.\end{array}$$

(37)

This is a cubic with respect to ${\omega }^2$, but can be further simplified in various limits.

In the hydrodynamic case it gives ${\omega }^2=\left[{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}\right]{\displaystyle \frac{k^2}{k_{\varpi }^2+k^2}}$, and we see that compressibility destabilizes. Inside the bracket we recognize the square of the epicyclic frequency and the square of the effective buoyancy (or Brunt-Väisälä) frequency in the compressible limit $N_{\mathrm{eff}}^2=-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}$, with the last term corresponding to the variation of density of a fluid parcel as it evolves adiabatically in comparison with the variation of the density of its surroundings. So the result is a mixture of CFI and buoyancy instability in the compressible limit. Note that, using the equilibrium condition in the hydrodynamic case $P_0^{\prime }=-{\rho }_0{g}_{\mathrm{eff}}$ we can write $N_{\mathrm{eff}}^2=-g_{\mathrm{eff}}{\left(ln{\displaystyle \frac{{\rho }_0}{P_0^{1/\mathrm{\Gamma }}}}\right)}^{\prime }$, and thus buoyancy destabilizes when ${\rho }_0/P_0^{1/\mathrm{\Gamma }}$ increases in the direction opposite to the effective gravity.

For zero rotation we recover the known result ${\omega }^2={\displaystyle \frac{k^2{N}^2}{k_{\varpi }^2+k^2}}$ with $N^2=-{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g^2}{c_s^2}}=g{\left(ln{\displaystyle \frac{P_0^{1/\mathrm{\Gamma }}}{{\rho }_0}}\right)}^{\prime }$, and the Schwarzschild instability criterion for convection ${\left(ln{\displaystyle \frac{P_0^{1/\mathrm{\Gamma }}}{{\rho }_0}}\right)}^{\prime }<0$.

With rotation but without gravity (when the effective gravity is only the centrifugal $g_{\mathrm{eff}}=-{\Omega }^2\varpi$) we get ${\omega }^2=\left[{\displaystyle \frac{{\left({\rho }_0{\varpi }^4{\Omega }^2\right)}^{\prime }}{{\rho }_0{\varpi }^3}}-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2}}\right]{\displaystyle \frac{k^2}{k_{\varpi }^2+k^2}}$, a relation that gives the instability condition ${\displaystyle \frac{dln\left({\rho }_0{\varpi }^4{\Omega }^2\right)}{dln\varpi }}<{\displaystyle \frac{{\Omega }^2{\varpi }^2}{c_s^2}}$ found in Ref. [5] using a heuristic analysis.

In the more general hydrodynamic case with gravity included the instability condition becomes $N_{\mathrm{eff}}^2+{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3<0$, or, ${\displaystyle \frac{{\left({\rho }_0{\varpi }^4{\Omega }^2\right)}^{\prime }}{{\rho }_0{\varpi }^3}}-{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}<{\displaystyle \frac{{(g-{\Omega }^2\varpi )}^2}{c_s^2}}$. Notably, it does not involve the Mach number $\Omega \varpi /c_s$, as in the case of zero gravity, but only the sum of the squares of the epicyclic frequency and of the effective buoyancy frequency. The criterion can also be written as ${\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}+g_{\mathrm{eff}}{\left(ln{\displaystyle \frac{P_0^{1/\mathrm{\Gamma }}}{{\rho }_0}}\right)}^{\prime }<0$, generalizing the Schwarzschild instability criterion for convection in the presence of rotation.

With magnetic field included we can find a simple approximate dispersion relation in the limit of relatively weak magnetic field for which $k^2{v}_{\mathrm{A}}^2{c}_s^2-(c_s^2+v_{\mathrm{A}}^2){\omega }^2\approx c_s^2(k^2{v}_{\mathrm{A}}^2-{\omega }^2)$. In that limit we get ${\displaystyle \frac{k_{\varpi }^2+k^2}{k^2}}{({\omega }^2-k^2{v}_{\mathrm{A}}^2)}^2-\left[{\displaystyle \frac{{\left({\varpi }^4{\Omega }^2\right)}^{\prime }}{{\varpi }^3}}+N_{\mathrm{eff}}^2\right]({\omega }^2-k^2{v}_{\mathrm{A}}^2)-4{\Omega }^2{k}^2{v}_{\mathrm{A}}^2=0$, where $N_{\mathrm{eff}}^2=-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}$. Comparing with equation (27) we see that only the full expression of the buoyancy frequency has replaced the one in the incompressible limit. Only the smallest root of the above quadratic can be negative, i.e., the

$$\begin{array}{c}\hfill {\omega }^2=k^2{v}_{\mathrm{A}}^2+{\displaystyle \frac{N_{\mathrm{eff}}^2+{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3}{2(1+k_{\varpi }^2/k^2)}}-\sqrt{{\left[{\displaystyle \frac{N_{\mathrm{eff}}^2+{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3}{2(1+k_{\varpi }^2/k^2)}}\right]}^2+{\displaystyle \frac{4{\Omega }^2{k}^2{v}_{\mathrm{A}}^2}{1+k_{\varpi }^2/k^2}}}.\end{array}$$

(38)

It is straightforward to show that this is indeed negative for wavenumbers $k_{\varpi }^2+k^2<{\displaystyle \frac{-N_{\mathrm{eff}}^2-\varpi {\left({\Omega }^2\right)}^{\prime }}{v_{\mathrm{A}}^2}}$, provided that the necessary condition for instability $N_{\mathrm{eff}}^2+\varpi {\left({\Omega }^2\right)}^{\prime }<0$ holds.

The fact that $N_{\mathrm{eff}}^2=-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}$ is smaller compared to its incompressible value $-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}$ and that ${\omega }^2$ is an increasing function of $N_{\mathrm{eff}}^2$ means that compressibility (i.e., to have finite $c_s$ in comparison with $c_s\to \infty$, while keeping everything else the same) always destabilizes.

The resulting dispersion relation is a mixture of MRI corresponding to the differential rotation, and magnetic buoyancy (or Parker) instability [6] corresponding to $N_{\mathrm{eff}}^2$, which indirectly depends on the magnetic field through the equilibrium condition ${\mathrm{\Pi }}_0^{\prime }=-{\rho }_0{g}_{\mathrm{eff}}$. The twin role of the magnetic field is evident in the instability condition as well as in the dispersion relation, since it appears both through $v_{\mathrm{A}}$, but also indirectly inside $N_{\mathrm{eff}}^2$. The appearance of $v_{\mathrm{A}}$ is connected to the stabilizing nature of the magnetic tension and is responsible for the fact that unstable modes exist up to some maximum wavenumber. To see the effect of the magnetic field through the effective buoyancy frequency we need to think if and how it affects the density stratification by solving the equilibrium condition. For example in the subcase of a polytropic relation in the unperturbed state characterized by a constant ${\mathrm{\Pi }}_0/{\rho }_0^{{\mathrm{\Gamma }}_{\mathrm{pol}}}$ and constant index ${\mathrm{\Gamma }}_{\mathrm{pol}}$ (as is usually assumed in studies of the Parker instability) the equilibrium condition ${\mathrm{\Pi }}_0^{\prime }=-{\rho }_0{g}_{\mathrm{eff}}$ gives ${\rho }_0^{\prime }=-{\displaystyle \frac{{\rho }_0^2{g}_{\mathrm{eff}}}{{\mathrm{\Gamma }}_{\mathrm{pol}}{\mathrm{\Pi }}_0}}$, and we get $N_{\mathrm{eff}}^2={\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}\left({\displaystyle \frac{\mathrm{\Gamma }/{\mathrm{\Gamma }}_{\mathrm{pol}}}{1+B_0^2/2P_0}}-1\right)$, showing the destabilization effect of the magnetic field (increasing $B_0^2$ decreases $N_{\mathrm{eff}}^2$). Essentially, a strong magnetic field in the case of constant ${\mathrm{\Pi }}_0/{\rho }_0^{{\mathrm{\Gamma }}_{\mathrm{pol}}}$ decreases the effect of the stable stratification that corresponds to the positive term $-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}$ in $N_{\mathrm{eff}}^2$.

Needless to say, there are other interesting cases that warrant further exploration . For example in the absence of rotation equation (37) gives ${\displaystyle \frac{k_{\varpi }^2+k^2}{k^2}}(c_s^2+v_{\mathrm{A}}^2){({\omega }^2-k^2{v}_{\mathrm{A}}^2)}^2+\left[(k_{\varpi }^2+k^2)v_{\mathrm{A}}^4+{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}(c_s^2+v_{\mathrm{A}}^2)+g^2\right]({\omega }^2-k^2{v}_{\mathrm{A}}^2)+{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}{k}^2{v}_{\mathrm{A}}^4=0$, corresponding to the pure magnetic buoyancy (or Parker) instability, for which we find that only the smaller root of the quadratic can become negative, ${\omega }^2={\displaystyle \frac{k^2{v}_{\mathrm{A}}^2{c}_s^2}{c_s^2+v_{\mathrm{A}}^2}}+{\displaystyle \frac{(k_{\varpi }^2+k^2)v_{\mathrm{A}}^4-{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}(c_s^2+v_{\mathrm{A}}^2)-g^2}{2(c_s^2+v_{\mathrm{A}}^2)(1+k_{\varpi }^2/k^2)}}-\sqrt{{\left[{\displaystyle \frac{(k_{\varpi }^2+k^2)v_{\mathrm{A}}^4-{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}(c_s^2+v_{\mathrm{A}}^2)-g^2}{2(c_s^2+v_{\mathrm{A}}^2)(1+k_{\varpi }^2/k^2)}}\right]}^2+{\displaystyle \frac{g^2{k}^4{v}_{\mathrm{A}}^4}{{(c_s^2+v_{\mathrm{A}}^2)}^2(k_{\varpi }^2+k^2)}}}$, provided that the instability condition $c_s^2\left[(k_{\varpi }^2+k^2)v_{\mathrm{A}}^2-{\displaystyle \frac{g{\rho }_0^{\prime }}{{\rho }_0}}\right]<g^2$ is satisfied. This instability condition in the subcase of constant ${\mathrm{\Pi }}_0/{\rho }_0^{{\mathrm{\Gamma }}_{\mathrm{pol}}}$ (with constant ${\mathrm{\Gamma }}_{\mathrm{pol}}$), when the equilibrium ${\mathrm{\Pi }}_0^{\prime }=-{\rho }_{0}g$ gives ${\rho }_0^{\prime }=-{\displaystyle \frac{{\rho }_0^{2}g}{{\mathrm{\Gamma }}_{\mathrm{pol}}{\mathrm{\Pi }}_0}}$, becomes ${\displaystyle \frac{(k_{\varpi }^2+k^2)v_{\mathrm{A}}^2{c}_s^2}{g^2}}<1-{\displaystyle \frac{\mathrm{\Gamma }{P}_0}{{\mathrm{\Gamma }}_{\mathrm{pol}}{\mathrm{\Pi }}_0}}$, and requires ${\displaystyle \frac{\mathrm{\Gamma }}{{\mathrm{\Gamma }}_{\mathrm{pol}}}}<1+{\displaystyle \frac{B_0^2}{2P_0}}$, i.e., for given positive ${\mathrm{\Gamma }}_{\mathrm{pol}}$ it requires sufficiently strong magnetic field corresponding to sufficiently small plasma beta (or, alternatively, thermodynamics corresponding to a sufficiently small index $\mathrm{\Gamma }$), a known characteristic of the Parker instability.

The limit of negligible thermal pressure is also an interesting case to examine, similarly to Refs. [16,17].

It is also possible to study the general properties of the cubic dispersion relation (37) without any further approximation (to alleviate the assumption of relatively weak magnetic field that we consider above in the rotating cases, but also include finite thermal pressure). Even more generally, one could include azimuthal magnetic field as in Ref. [18], and also consider nonaxisymmetric perturbations ($m\ne 0)$, as in Ref. [19]. These generalizations are beyond the scope of the present paper.

5.2. Approximate Dispersion Relation for Interface-Driven Instability

Similarly to the incompressible case (Section 4.4) we can use the same variant of WKBJ to find the dispersion relation of an interface-driven instability. The constant solutions of the principal equation are ${\displaystyle \frac{1}{Y}}={\displaystyle \frac{f_{11}\mp \sqrt{-{\kappa }_0^2}}{-f_{12}}}$ with ${\kappa }_0^2=-f_{11}^2-f_{12}{f}_{21}$ given by equation (A41), and correspond to ${\displaystyle \frac{y_1^{\prime }}{y_1}}={\displaystyle \frac{y_2^{\prime }}{y_2}}=\mp \sqrt{-{\kappa }_0^2}$. Applying the above in the two sides of the interface, the upper sign should be used in the region $\varpi >{\varpi }_{\mathrm{CD}}$ and the lower sign in the region $\varpi <{\varpi }_{\mathrm{CD}}$, such that we get exponential decrease as we move away from the interface in both sides. The dispersion relation comes from the continuity of Y.

Using the expression of ${\kappa }_0^2$ for axisymmetric disturbances ($m=0$) of rotating flows in which the unperturbed magnetic field has only $\widehat{z}$ component and the gravity is zero, i.e.,

$$\begin{array}{c}\hfill {\kappa }_0^2=-{\displaystyle \frac{(v_{\mathrm{A}}^2-{\omega }^2/k^2)(c_s^2-{\omega }^2/k^2){\mathrm{\Lambda }}^2-{\Omega }^4{\varpi }^2}{v_{\mathrm{A}}^2{c}_s^2-(c_s^2+v_{\mathrm{A}}^2){\omega }^2/k^2}},\mathrm{\Lambda }=\sqrt{k^2-{\displaystyle \frac{4{\omega }^2{\Omega }^2{k}^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}},\end{array}$$

(39)

we get

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{Y}}={\displaystyle \frac{{\rho }_0{\Omega }^2{c}_s^2}{c_s^2-{\omega }^2/k^2}}\pm {\displaystyle \frac{{\rho }_0}{\varpi }}\sqrt{\left[(v_{\mathrm{A}}^2-{\omega }^2/k^2){\mathrm{\Lambda }}^2-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2-{\omega }^2/k^2}}\right]\left[v_{\mathrm{A}}^2-{\displaystyle \frac{c_s^2{\omega }^2/k^2}{c_s^2-{\omega }^2/k^2}}\right]}.\end{array}$$

(40)

Thus the dispersion relation is

$$\begin{array}{c}\hfill {\left[{\displaystyle \frac{{\rho }_0{\Omega }^2\varpi c_s^2}{c_s^2-{\omega }^2/k^2}}-{\rho }_0\sqrt{\left[(v_{\mathrm{A}}^2-{\omega }^2/k^2){\mathrm{\Lambda }}^2-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2-{\omega }^2/k^2}}\right]\left[v_{\mathrm{A}}^2-{\displaystyle \frac{c_s^2{\omega }^2/k^2}{c_s^2-{\omega }^2/k^2}}\right]}\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{-}}=\\ \hfill {\left[{\displaystyle \frac{{\rho }_0{\Omega }^2\varpi c_s^2}{c_s^2-{\omega }^2/k^2}}+{\rho }_0\sqrt{\left[(v_{\mathrm{A}}^2-{\omega }^2/k^2){\mathrm{\Lambda }}^2-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2-{\omega }^2/k^2}}\right]\left[v_{\mathrm{A}}^2-{\displaystyle \frac{c_s^2{\omega }^2/k^2}{c_s^2-{\omega }^2/k^2}}\right]}\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{+}}.\end{array}$$

(41)

In most cases of interest in the limit of large $\left|k\right|{\varpi }_{\mathrm{CD}}$ we may approximate $c_s^2-{\omega }^2/k^2\approx c_s^2$, and $v_{\mathrm{A}}^2{c}_s^2-(c_s^2+v_{\mathrm{A}}^2){\omega }^2/k^2\approx c_s^2(v_{\mathrm{A}}^2-{\omega }^2/k^2)$, and the expression simplifies to

$$\begin{array}{ccc}& {\left[{\rho }_0{\Omega }^2\varpi -{\rho }_0{\displaystyle \frac{k^2{v}_{\mathrm{A}}^2-{\omega }^2}{\left|k\right|}}\sqrt{1-{\displaystyle \frac{4{\Omega }^2{\omega }^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}}\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{-}}& =\hfill \\ & {\left[{\rho }_0{\Omega }^2\varpi +{\rho }_0{\displaystyle \frac{k^2{v}_{\mathrm{A}}^2-{\omega }^2}{\left|k\right|}}\sqrt{1-{\displaystyle \frac{4{\Omega }^2{\omega }^2}{{(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}^2}}-{\displaystyle \frac{{\Omega }^4{\varpi }^2}{c_s^2(k^2{v}_{\mathrm{A}}^2-{\omega }^2)}}}\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{+}}& .\hfill \end{array}$$

(42)

A necessary condition for instability is clearly ${\left[{\rho }_0{\Omega }^2\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{-}}>{\left[{\rho }_0{\Omega }^2\right]}_{\varpi ={\varpi }_{\mathrm{CD}}^{+}}$.

5.3. Exact Results

To examine how compressibility affects the growth rate of the rotational instability, we consider axisymmetric perturbations of a rotating flow, assuming magnetic field $B_0\left(\varpi \right)\widehat{z}$ and zero gravity. For this case $f_{11}={\displaystyle \frac{k^2{c}_s^2{\Omega }^2\varpi }{\left(c_s^2+v_{\mathrm{A}}^2\right){\omega }^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}$, $f_{12}={\displaystyle \frac{({\omega }^2-k^2{c}_s^2)\varpi /{\rho }_0}{(c_s^2+v_{\mathrm{A}}^2){\omega }^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}$, $f_{21}=-{\displaystyle \frac{{\rho }_0({\omega }^2-k^2{v}_{\mathrm{A}}^2)}{\varpi }}+{\displaystyle \frac{4{\rho }_0{\omega }^2{\Omega }^2/\varpi }{{\omega }^2-k^2{v}_{\mathrm{A}}^2}}+{\displaystyle \frac{k^2{\rho }_0(c_s^2+v_{\mathrm{A}}^2){\Omega }^4\varpi }{\left(c_s^2+v_{\mathrm{A}}^2\right){\omega }^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}$, and the principal equation is

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\varpi }}\left({\displaystyle \frac{1}{Y}}\right)={\displaystyle \frac{({\omega }^2-k^2{c}_s^2)\varpi /{\rho }_0}{(c_s^2+v_{\mathrm{A}}^2){\omega }^2-k^2{c}_s^2{v}_{\mathrm{A}}^2}}{\left({\displaystyle \frac{1}{Y}}+{\displaystyle \frac{{\rho }_0{k}^2{c}_s^2{\Omega }^2}{{\omega }^2-k^2{c}_s^2}}\right)}^2+{\displaystyle \frac{{\rho }_0({\omega }^2-k^2{v}_{\mathrm{A}}^2)}{\varpi }}-{\displaystyle \frac{4{\rho }_0{\omega }^2{\Omega }^2/\varpi }{{\omega }^2-k^2{v}_{\mathrm{A}}^2}}-{\displaystyle \frac{{\rho }_0{k}^2{\Omega }^4\varpi }{{\omega }^2-k^2{c}_s^2}}.\end{array}$$

(43)

Compared to the incompressible case there is an additional difficulty connected to the requirement the total pressure ${\mathrm{\Pi }}_0={\displaystyle \frac{{\rho }_0{c}_s^2}{\mathrm{\Gamma }}}+{\displaystyle \frac{{\rho }_0{v}_{\mathrm{A}}^2}{2}}$ of the unperturbed state to satisfy the equilibrium condition ${\mathrm{\Pi }}_0^{\prime }={\rho }_0{\Omega }^2\varpi$. Still, it is possible in some cases to find analytical expressions for the solution of the principal equation and the dispersion relation under certain assumptions, see, e.g. Appendix D for the case of uniform angular velocity, and Appendix E for the case of uniform rotational velocity (in both cases constant $c_s$ and $v_{\mathrm{A}}$ are assumed). The numerical integration is straightforward though, and covers all cases.

Here we consider two cases, extending previous results by taking compressibility into account. In the first, we assume a rotating shell between two solid cylindrical boundaries, ${\varpi }_1<\varpi <{\varpi }_2$, in order to explore how compressibility modifies the results of Section 3.2.

In the second case, we assume a uniformly rotating cylindrical flow in $\varpi <{\varpi }_{\mathrm{CD}}$ surrounded by a static homogeneous plasma in $\varpi >{\varpi }_{\mathrm{CD}}$, generalizing the results of Section 4.2 (for a non-rotating outer part, $\mathcal{A}=1$). The static environment is characterized by its constant density ${\rho }_2$, sound velocity $c_{s2}$ and Alfvén velocity $v_{\mathrm{A}2}$ corresponding to uniform magnetic field $v_{\mathrm{A}2}\sqrt{{\rho }_2}\widehat{z}$. These need only to satisfy the pressure balance at the interface in the unperturbed state, i.e., continuity of ${\mathrm{\Pi }}_0={\displaystyle \frac{{\rho }_2{c}_{s2}^2}{{\mathrm{\Gamma }}_2}}+{\displaystyle \frac{{\rho }_2{v}_{\mathrm{A}}^2}{2}}$. Here we choose same density and magnetic field in the two regions, same polytropic index $\mathrm{\Gamma }=5/3$, and continuous $c_s$ at the interface.

In both cases we keep the unperturbed state the same as in the incompressible analogue, with the only difference being that the sound velocity (and the thermal pressure) is now finite, and has the proper spatial dependence required fo the equilibrium ${\mathrm{\Pi }}_0^{\prime }={\rho }_0{\Omega }^2\varpi$ with ${\mathrm{\Pi }}_0={\displaystyle \frac{{\rho }_0{c}_s^2}{\mathrm{\Gamma }}}+{\displaystyle \frac{{\rho }_0{v}_{\mathrm{A}}^2}{2}}$ to be satisfied.

In particular, for the rotating shell we assume uniform density ${\rho }_0$, uniform magnetic field $v_{\mathrm{A}}\sqrt{{\rho }_0}\widehat{z}$, and power-law rotation $\Omega ={\Omega }_2{\left({\displaystyle \frac{\varpi }{{\varpi }_2}}\right)}^b$, in which case integration of the equilibrium condition gives $c_s=\sqrt{c_{s1}^2+{\displaystyle \frac{\mathrm{\Gamma }{\Omega }_2^2}{{\varpi }_2^{2b}}}{\displaystyle \frac{{\varpi }^{2b+2}-{\varpi }_1^{2b+2}}{2b+2}}}$, where $c_{s_1}$ is the sound velocity at the inner radius ${\varpi }_1={\varpi }_2-\Delta \varpi$.

Similarly, for the rotating cylinder we assume uniform density ${\rho }_0$, uniform magnetic field $v_{\mathrm{A}}\sqrt{{\rho }_0}\widehat{z}$, and constant $\Omega$, when the integration of the equilibrium condition gives $c_s=\sqrt{c_{s1}^2+{\displaystyle \frac{\mathrm{\Gamma }{\Omega }^2{\varpi }^2}{2}}}$, where $c_{s_1}$ is the sound velocity at the symmetry axis. The thermal pressure resulting from the equilibrium is $P_0={\displaystyle \frac{{\rho }_0{\Omega }^2{\varpi }^2}{2}}+$ constant. The additive constant is determined by the value of the sound velocity $c_{s_1}$ at the symmetry axis, and thus the sound velocity as function of radius can be written as $c_s=\sqrt{c_{s1}^2+{\displaystyle \frac{\mathrm{\Gamma }{\Omega }^2{\varpi }^2}{2}}}$.

The principal equation can be integrated using the Schwarzian approach, and the eigenvalues are obtained from the requirement that ${\displaystyle \frac{{\mathrm{\Phi }|}_2{-\mathrm{\Phi }|}_1}{2\pi }}$ be an integer (with the subscripts 1 and 2 refer to the extreme values of $\varpi$ in the entire region occupied by plasma; ${\varpi }_1$ may be finite or zero, and ${\varpi }_2$ may be finite or infinity).

Figure 5 shows how the growth rate changes with decreasing $c_s$ in the case of a wall-bounded rotational instability of a flow between two cylindrical walls. In the top left panel the value of $c_{s1}$ is large enough and the growth rate is practically that of the incompressible limit (indeed, the top-left panel of Figure 5 is essentially identical to the top-left panel of Figure 1). The top-right panel of Figure 5 corresponds to the opposite extreme $c_{s1}=0$.

As expected from the WKBJ analysis of Section 5.1, compressibility is destabilizing, since it reduces the sum $N_{\mathrm{eff}}^2+{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3$ appearing in the approximate dispersion relation (38) by a factor ${\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}$, corresponding to a relative difference of the order ${\displaystyle \frac{g_{\mathrm{eff}}^2/c_s^2}{{\Omega }^2}}={\displaystyle \frac{{\varpi }^2{\Omega }^2}{c_s^2}}$. Although the unperturbed density is uniform, buoyancy effects are present once compressibility is included, and they modify the dispersion relation accordingly.

The bottom panels of Figure 5 show the stabilization effect of the magnetic field. In Section 5.1 we discussed the twin role of the magnetic field, which stabilizes through its tension, but may destabilize if it affects the density stratification (as a consequence of the equilibrium condition). In the present case, however, the destabilizing effect is absent because the unperturbed density was chosen to be a constant.

Figure 6 shows how the growth rate changes with decreasing $c_s$ in the case of a uniformly rotating cylinder surrounded by static plasma with same density and magnetic field. In the left panel the value of $c_{s1}$ is large enough that the growth rate is essentially that of the incompressible limit (indeed, the top-left panel of Figure 3 – the curve for $\mathcal{A}=1$ – is practically identical to the top-left panel of Figure 6). The right panels of Figure 6 correspond to the opposite extreme $c_{s1}=0$. Similarly to the rotating shell, we observe a higher growth rate and a larger range of unstable wavenumbers due to the influence of buoyancy, although in this case the effect is less pronounced.

For small enough $c_s$, additional modes are present, with distinctively different eigenfunctions, as seen in the bottom panel of Figure 6. The CFI corresponds to the first mode with the larger growth rate, whose eigenfunction is localized around the tangential discontinuity (and becomes increasingly localized as $c_s$ decreases).

As already discussed, compressibility is destabilizing, but in this case the relative decrease of the sum $N_{\mathrm{eff}}^2+{\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3$ is of order ${\displaystyle \frac{g_{\mathrm{eff}}^2/c_s^2}{4{\Omega }^2}}={\displaystyle \frac{{\varpi }_{\mathrm{CD}}^2{\Omega }^2}{4(c_{s1}^2+\mathrm{\Gamma }{\Omega }^2{\varpi }_{\mathrm{CD}}^2/2)}}<{\displaystyle \frac{1}{2\mathrm{\Gamma }}}$.

The newly appearing modes ($n\ge 2$) are body modes, attributable to the magnetic buoyancy/Parker-type instability that arises once compressibility is included, even though the unperturbed density is uniform. They can be qualitatively understood using the WKBJ results of Section 5.1, in particular the dispersion relation (38) which becomes ${\omega }^2=k^2{v}_{\mathrm{A}}^2+{\Omega }^2{\displaystyle \frac{4-{\Omega }^2{\varpi }^2/c_s^2}{2(1+k_{\varpi }^2/k^2)}}-\sqrt{{\left[{\Omega }^2{\displaystyle \frac{4-{\Omega }^2{\varpi }^2/c_s^2}{2(1+k_{\varpi }^2/k^2)}}\right]}^2+{\displaystyle \frac{4k^2{v}_{\mathrm{A}}^2{\Omega }^2}{1+k_{\varpi }^2/k^2}}}$, and the instability condition $k_{\varpi }^2+k^2<{\displaystyle \frac{{\Omega }^4{\varpi }^2}{v_{\mathrm{A}}^2{c}_s^2}}$. The latter explains why these modes are absent in the incompressible limit, while in the opposite extreme $c_{s1}=0$ they are present whenever ${(n-1)}^2{\pi }^2+k^2<{\displaystyle \frac{2{\Omega }^2}{\mathrm{\Gamma }{v}_{\mathrm{A}}^2}}$.

6. Summary

The main goal of this paper is to present the form of the principal equation and the associated boundary conditions for hydrodynamic and ideal magnetohydrodynamic flows in cylindrical geometry. These results, which enable the application of the minimalist approach to any cylindrical equilibrium – with radially varying density and pressure, axial and azimuthal components of both the velocity and magnetic field, and a radially directed gravitational field – are summarized in Table 1. Within this framework, integrating the principal equation alone is sufficient to obtain the dispersion relation of linear perturbations.

We applied this methodology to rotating magnetized flows, primarily for illustration rather than to exhaust all possible equilibrium configurations. Even so, the exact dispersion relations of the selected examples clearly reveal the underlying physical mechanisms, and the WKBJ approximation provides valuable insight. In particular, we derived approximate dispersion relations for axisymmetric perturbations in flows with axial magnetic fields, finding excellent agreement with exact solutions. These relations delineate the regimes in which centrifugal (CFI), magnetorotational (MRI), and buoyancy instabilities act individually or in combination.

CFI is governed by the epicyclic frequency, whose square is ${\left({\varpi }^4{\Omega }^2\right)}^{\prime }/{\varpi }^3$, and is active whenever this quantity is negative. In the presence of magnetic fields, CFI transitions to MRI, which is governed by the differential rotation frequency, $\varpi {\left({\Omega }^2\right)}^{\prime }$, and likewise operates when this quantity is negative. Buoyancy instability is controlled by the square of the effective buoyancy (or Brunt-Väisälä) frequency $N_{\mathrm{eff}}^2=-{\displaystyle \frac{g_{\mathrm{eff}}{\rho }_0^{\prime }}{{\rho }_0}}-{\displaystyle \frac{g_{\mathrm{eff}}^2}{c_s^2}}$, where rotation modifies gravity through $g_{\mathrm{eff}}=g-{\Omega }^2\varpi$. Instability occurs when $N_{\mathrm{eff}}^2<0$. In the incompressible limit this reduces to the familiar condition of unstable density stratification, while compressibility introduces an additional negative term associated with adiabatic density changes of fluid parcels, further enhancing buoyancy-driven instability.

In general, the stability of the system is determined by the sum of the squares of these characteristic frequencies; instability arises when the sum becomes negative.

Magnetic fields play a central role in the stability of rotating flows. First, flux freezing couples neighboring fluid parcels along field lines, forcing them to rotate together and effectively replacing the epicyclic frequency by the differential rotation frequency in the stability criteria. This transforms CFI into MRI and typically has a destabilizing effect – for example, even a weak magnetic field destabilizes an otherwise stable hydrodynamic Keplerian disk. Second, magnetic tension stabilizes sufficiently short-wavelength perturbations, implying that any instability operates only up to a finite wavenumber $k_{\max }$. Third, magnetic pressure can modify the equilibrium density stratification, altering the first term in $N_{\mathrm{eff}}^2$. This mechanism underlies magnetic-buoyancy (Parker) instability, where magnetic pressure weakens stabilizing stratification and promotes buoyant rise of magnetized regions.

We examined two broad classes of flows: wall-bounded configurations and flows containing tangential discontinuities. Both are generally susceptible to instabilities – global modes in the former and localized modes in the latter – provided the relevant criteria are satisfied. In the presence of discontinuities, essentially the same instability criteria remain valid when spatial derivatives are interpreted as sharp gradients (i.e. jumps) across the interface. For example, a discontinuity in rotational velocity is associated with CFI/MRI and behaves like a local Rayleigh-Taylor instability for axisymmetric perturbations when the wavenumber is normal to the velocity shear, but may exhibit Kelvin-Helmholtz character when $m\ne 0$. Current-driven and ballooning instabilities may also arise when azimuthal magnetic fields and appropriate pressure profiles are included.

The minimalist approach provides a unified and efficient way to obtain dispersion relations in all these situations, just as demonstrated in the examples of this paper. Furthermore, the alternative forms of the equations presented here may facilitate closed-form dispersion relations in additional cases.

Needless to say, although we selected rotational flows as the first applications of the minimalist approach in cylindrical geometry – focusing on axisymmetric modes and axial magnetic fields – the formalism summarized in Table 1 can be applied to far more general cylindrical equilibria of disks and jets. These include equilibria with axial and azimuthal magnetic fields, rotational and axial velocity, and studies of nonaxisymmetric modes as well. We plan to explore these extensions in future works.

Appendix C.1. Alternative Forms of the Principal Equation

The following alternative forms of the principal equation may be useful, particularly when seeking analytical solutions2

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{1}{Y}}\right)}^{\prime }=f_{12}{\left({\displaystyle \frac{1}{Y}}\right)}^2+2f_{11}\left({\displaystyle \frac{1}{Y}}\right)-f_{21}\iff \\ \hfill {\left({\displaystyle \frac{f_{12}}{Y}}+f_{11}+{\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}\right)}^{\prime }={\left({\displaystyle \frac{f_{12}}{Y}}+f_{11}+{\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}\right)}^2+{\kappa }_1^2,\end{array}$$

(A36)

$$\begin{array}{c}\hfill {\left(f_{21}Y-f_{11}+{\displaystyle \frac{f_{21}^{\prime }}{2f_{21}}}\right)}^{\prime }={\left(f_{21}Y-f_{11}+{\displaystyle \frac{f_{21}^{\prime }}{2f_{21}}}\right)}^2+{\kappa }_2^2,\end{array}$$

(A37)

with

$$\begin{array}{c}\hfill {\kappa }_1^2={\kappa }_0^2+f_{12}{\left({\displaystyle \frac{f_{11}}{f_{12}}}\right)}^{\prime }-{\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}\right)}^2+{\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}\right)}^{\prime },\end{array}$$

(A38)

$$\begin{array}{c}\hfill {\kappa }_2^2={\kappa }_0^2-f_{21}{\left({\displaystyle \frac{f_{11}}{f_{21}}}\right)}^{\prime }-{\left({\displaystyle \frac{f_{21}^{\prime }}{2f_{21}}}\right)}^2+{\left({\displaystyle \frac{f_{21}^{\prime }}{2f_{21}}}\right)}^{\prime },\end{array}$$

(A39)

$$\begin{array}{c}\hfill {\kappa }_0^2=-f_{11}^2-f_{12}{f}_{21}.\end{array}$$

(A40)

After substituting the expressions of $f_{ij}$ we can write ${\kappa }_0^2$ in the form

$$\begin{array}{c}\hfill {\kappa }_0^2={\tilde{\kappa }}^2-{\displaystyle \frac{{\tilde{\kappa }}^2{\rho }_0\varpi }{A}}{\left({\displaystyle \frac{g}{\varpi }}\right)}^{\prime }+{\displaystyle \frac{4k^2{T}^2}{{\varpi }^2{A}^2}}\\ \hfill -{\displaystyle \frac{k^2}{S/{\rho }_0^2}}{\left(g_{\mathrm{eff}}+{\displaystyle \frac{2{\omega }_0{B}_{0\varphi }W}{\varpi A}}\right)}^2-{\displaystyle \frac{1}{S/{\rho }_0^2}}{\left[{\displaystyle \frac{2T{\omega }_0^2}{\varpi A}}-{\displaystyle \frac{m}{\varpi }}\left(g_{\mathrm{eff}}+{\displaystyle \frac{2{\omega }_0{B}_{0\varphi }W}{\varpi A}}\right)\right]}^2,\end{array}$$

(A41)

and the expression of ${\kappa }_1^2$ as

$$\begin{array}{c}\hfill {\kappa }_1^2={\tilde{\kappa }}^2+{\displaystyle \frac{{\tilde{\kappa }}^2\varpi }{A}}{\left({\displaystyle \frac{B_{0\varphi }^2}{{\varpi }^2}}\right)}^{\prime }-{\displaystyle \frac{4{\tilde{\kappa }}^{2}F{B}_{0\varphi }T}{{\varpi }^2{A}^2}}+{\displaystyle \frac{{\tilde{\kappa }}^2{g}_{\mathrm{eff}}{\rho }_0^{\prime }}{A}}-{\displaystyle \frac{k^2+m^2/{\varpi }^2}{S/{\rho }_0^2}}{\left(g_{\mathrm{eff}}+{\displaystyle \frac{2{\omega }_0{B}_{0\varphi }W}{\varpi A}}\right)}^2\\ \hfill -{\displaystyle \frac{{\tilde{\kappa }}^2{\rho }_0{\left({\varpi }^2{V}_{0\varphi }^2\right)}^{\prime }}{A{\varpi }^3}}-{\displaystyle \frac{4{\tilde{\kappa }}^2{\rho }_{0}F{V}_{0\varphi }W}{{\varpi }^2{A}^2}}+{\displaystyle \frac{2{\omega }_0^2}{S/{\rho }_0^2}}\left(g_{\mathrm{eff}}+{\displaystyle \frac{2{\omega }_0{B}_{0\varphi }W}{\varpi A}}\right)\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}+{\displaystyle \frac{2mT}{{\varpi }^{2}A}}\right)\\ \hfill -{\left({\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}+{\displaystyle \frac{2mT}{{\varpi }^{2}A}}\right)}^2+{\left[{\displaystyle \frac{f_{12}^{\prime }}{2f_{12}}}+{\displaystyle \frac{2mT}{{\varpi }^{2}A}}-{\displaystyle \frac{{\omega }_0^2}{S/{\rho }_0^2}}\left(g_{\mathrm{eff}}+{\displaystyle \frac{2{\omega }_0{B}_{0\varphi }W}{\varpi A}}\right)\right]}^{\prime }.\end{array}$$

(A42)

Appendix C.2. Sturm-Liouville Formulation

We can reformulate the system (14) of the two first-order differential equations for $y_1$ and $y_2$, into a single second-order differential equation for either $y_1$ or $y_2$, and obtain Y afterwards:

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{y_1^{\prime }}{f_{12}}}\right)}^{\prime }+\left[-f_{11}^2-f_{12}{f}_{21}+f_{12}{\left({\displaystyle \frac{f_{11}}{f_{12}}}\right)}^{\prime }\right]{\displaystyle \frac{y_1}{f_{12}}}=0,{\displaystyle \frac{1}{Y}}={\displaystyle \frac{-f_{11}-y_1^{\prime }/y_1}{f_{12}}},\end{array}$$

(A43)

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{y_2^{\prime }}{f_{21}}}\right)}^{\prime }+\left[-f_{11}^2-f_{12}{f}_{21}-f_{21}{\left({\displaystyle \frac{f_{11}}{f_{21}}}\right)}^{\prime }\right]{\displaystyle \frac{y_2}{f_{21}}}=0,Y={\displaystyle \frac{f_{11}-y_2^{\prime }/y_2}{f_{21}}}.\end{array}$$

(A44)

Appendix C.3. “Variable Frequency Oscillator” Formulation

As explained in Ref. [10] we can also reformulate the above as “variable frequency oscillator” equations

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{y_1}{\sqrt{f_{12}}}}\right)}^{″}+{\kappa }_1^2{\displaystyle \frac{y_1}{\sqrt{f_{12}}}}=0,{\displaystyle \frac{1}{Y}}={\displaystyle \frac{-f_{11}-y_1^{\prime }/y_1}{f_{12}}},\end{array}$$

(A45)

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{y_2}{\sqrt{f_{21}}}}\right)}^{″}+{\kappa }_2^2{\displaystyle \frac{y_2}{\sqrt{f_{21}}}}=0,Y={\displaystyle \frac{f_{11}-y_2^{\prime }/y_2}{f_{21}}}.\end{array}$$

(A46)

The above expressions give directly the local radial wavenumber $k_{\varpi }$ in the WKBJ approximation. The corresponding approximate dispersion relation can therefore be obtained either from equation (A38) by setting ${\kappa }_1=k_{\varpi }$, or from equation (A39) by setting ${\kappa }_2=k_{\varpi }$. Notably, these two choices lead to distinct WKBJ approximations of the dispersion relation; see the related discussion in §4 of Ref. [15]. As explained there, a third option is to apply the WKBJ method to the system associated with ${\kappa }_0=k_{\varpi }$. The appropriate choice among ${\kappa }_1$, ${\kappa }_2$, and ${\kappa }_0$ depends on the relative smallness of the quantities ${\kappa }_1^{\prime }/{\kappa }_1^2$, ${\kappa }_2^{\prime }/{\kappa }_2^2$, and ${\kappa }_0^{\prime }/{\kappa }_0^2$, respectively.