# Rossby wave instability in weakly ionized protoplanetary disks. I. azimuthal
or vertical B-fields

Can Cui1,2, Ashutosh Tripathi2, Cong Yu3,4, Min-Kai Lin5,6 and Andrew
Youdin7,8

1Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON
M5S 3H4, Canada
2DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
3School of Physics and Astronomy, Sun Yat-Sen University, Zhuhai 519082, China
4CSST Science Center for the Guangdong-Hong Kong-Macau Greater Bay Area,
Zhuhai 519082, China
5Institute of Astronomy and Astrophysics, Academia Sinica,Taipei 10617,
Taiwan, R.O.C.
6Physics Division, National Center for Theoretical Sciences, Taipei 10617,
Taiwan, R.O.C.
7Department of Astronomy and Steward Observatory, University of Arizona,
Tucson, AZ 85721, USA
8The Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721,
USA
<EMAIL_ADDRESS>
###### Abstract

Rossby wave instability (RWI) is considered the underlying mechanism to
crescent-shaped azimuthal asymmetries, discovered in (sub-)millimeter dust
continuum of many protoplanetary disks. Previous works on linear theory were
conducted in the hydrodynamic limit. Nevertheless, protoplanetary disks are
likely magnetized and weakly ionized. We examine the influence of magnetic
fields and non-ideal magnetohydrodynamic (MHD) effects - namely, Ohmic
resistivity, Hall drift, and ambipolar diffusion - on the RWI unstable modes.
We perform radially global linear analyses, employing constant azimuthal
($B_{\phi}$) or vertical ($B_{z}$) background magnetic fields. It is found
that, in the ideal MHD regime, magnetism can either enhance or diminish RWI
growth. Strong non-ideal MHD effects cause RWI growth rates to recover
hydrodynamic results. The sign of Hall Elsässer number subtly complicates the
results, and vertical wavenumbers generically diminish growth rates.

###### keywords:

hydrodynamics – MHD – methods: analytical – protoplanetary disks

††pubyear: 2023††pagerange: Rossby wave instability in weakly ionized
protoplanetary disks. I. azimuthal or vertical B-fields–B

##  1 Introduction

Protoplanetary disks are composed of gas and dust orbiting pre-main-sequence
stars (Armitage, 2011). They are the birth place of planets, where micron-
sized dust grains coalesce and evolve into km-sized planetesimals, eventually
giving rise to terrestrial planets or gas giants cores. Nevertheless, the
growth of dust grains faces several barriers, including bouncing,
fragmentation, and fast radial drift (Weidenschilling, 1977; Güttler et al.,
2010; Zsom et al., 2010). Rossby wave instability (RWI) is perhaps a promising
mechanism to circumvent these barriers. Its non-linear state generates large,
lopsided crescent-shaped vortices, which concentrate grains towards pressure
maximum, fostering streaming instability and subsequent gravitational collapse
(Goldreich & Ward, 1973; Youdin & Goodman, 2005). These vortices manifest as
azimuthal asymmetries observed in (sub-)millimeter dust continuum as well as
CO rotational transition lines (e.g., Huang et al., 2018; van der Marel et
al., 2021). Notable examples include IRS 48 (van der Marel et al., 2013),
HD142527 (Muto et al., 2015), and AB Aur (Fuente et al., 2017).

The RWI is triggered by local extrema in the radial profile of vortensity,
$(\Sigma\Omega/\kappa^{2})S^{2/\Gamma}$ (Lovelace et al., 1999; Li et al.,
2000). It gives exponential growth of non-axisymmetric modes
($\propto\exp[im\phi]$, where $m=1,2,...$) on each side of the corotation
radius. These unstable Rossby modes are confined between the inner and outer
Lindblad resonances, where density waves are launched and propogate away from
the resonances (Lin & Shu, 1964; Tsang & Lai, 2008). The eigenfunctions of
Rossby waves manifest as anticyclonic vortices, characterized by vorticity
($\nabla\times\delta\mathbf{u}$) directed oppositely to the rotation of the
disk, yielding a maximum perturbed pressure at the vortex core (Lovelace &
Romanova, 2014). Extensions of the classic RWI incorporating self-gravity (Lin
& Papaloizou, 2011; Lin, 2012), magnetic fields (Yu & Li, 2009; Yu & Lai,
2013), and dust-gas two-fluid drag (Liu & Bai, 2023) have been explored.

In the context of protoplanetary disks, numerical simulations of the RWI have
been commonly performed at the gap edges carved by a planet (Zhu et al., 2014;
Zhu & Stone, 2014; Hammer et al., 2017; Huang et al., 2018; Cimerman &
Rafikov, 2023), or at the dead zone edges of the magneto-rotational
instability (MRI; Lyra & Mac Low, 2012; Miranda et al., 2016, 2017), for which
local vortensity extrema take place. These simulations elucidate that the non-
linear saturation of the Rossby vortices is primarily governed by $m=1$ modes
(Godon & Livio, 1999; Li et al., 2001; Meheut et al., 2012). Furthermore, the
long-term survival of Rossby vortices are examined under secondary
instabilities (e.g., Lesur & Papaloizou, 2009), thermal relaxation (e.g., Fung
& Ono, 2021; Rometsch et al., 2021), and dust-gas two-fluid framework (e.g.,
Surville et al., 2016; Surville & Mayer, 2019; Li et al., 2020).

While the RWI is inherently a hydrodynamic phenomenon, it is crucial to
acknowledge the magnetized nature of protoplanetary disks, where magnetic
fields likely originate from primordial molecular clouds (Galli & Shu, 1993;
Girart et al., 2006; Girart et al., 2009). Owing to the weak thermal
ionization of passively heated protoplanetary disks and the modest non-thermal
ionization by, for instance, the stellar FUV, EUV, X-rays, and cosmic rays
(Wardle, 2007; Bai, 2011), the gas and magnetic fields are not perfectly
coupled. The disk material may be weakly ionized, and three non-ideal MHD
effects come into play: Ohmic resistivity, the Hall drift, and ambipolar
diffusion (Weiss et al., 2021; Lesur, 2021). At the same ionization fraction,
Ohmic, Hall, and ambipolar diffusivities are respectively proportional to
$\eta_{\mathrm{O}}\propto\mathrm{const.}$, $\eta_{\mathrm{H}}\propto B/\rho$,
and $\eta_{\mathrm{A}}\propto B^{2}/\rho^{2}$. Consequently, Ohmic resistivity
dominates in high-density regions near the inner disk and at the midplane,
while ambipolar diffusion becomes prominent in the outer disk. Hall drift is
most effective between these regions. Both numerical simulations and
analytical theory have underscored the significance of non-ideal MHD effects
in shaping disk structure and evolution, such as dead zone structures (e.g.,
Flock et al., 2017), magnetized wind kinematics (e.g., Wang et al., 2019),
dust dynamics and distribution (e.g., Lin & Hsu, 2022; Wu et al., 2024),
heating and disk temperature profiles (e.g., Béthune & Latter, 2020), and the
formation of annular substructures (e.g., Riols & Lesur, 2019; Cui & Bai,
2021, 2022).

Previously, the non-ideal MHD effects were incorporated in the linear studies
of several instabilities applicable to protoplanetary disks, for example, the
MRI (Balbus & Hawley, 1991) and the vertical shear instability (Goldreich &
Schubert, 1967; Nelson et al., 2013). These investigations unveiled
significant alterations in the unstable modes induced by non-ideal MHD
physics. In the context of MRI, both Ohmic resistivity and ambipolar diffusion
tend to stabilize the instability, although under specific conditions,
ambipolar diffusion shear instability can be excited (Blaes & Balbus, 1994;
Kunz & Balbus, 2004; Latter & Kunz, 2022). The impact of Hall drift is
twofold, based upon the sign of the Hall Elsässer number (Wardle, 1999; Pandey
& Wardle, 2012). Concerning the vertical shear instability, strong ionization
can substantially weaken the unstable modes, while subdued modes regain
strength under weak ionizations (Latter & Papaloizou, 2018; Cui & Bai, 2020;
Cui & Lin, 2021; Latter & Kunz, 2022).

To deepen our comprehension of RWI behavior in magnetized disks,
investigations must encompass both ideal and non-ideal MHD physics. Previous
studies have primarily focused on the ideal MHD regime. Yu & Li (2009)
utilized Lagrangian perturbation theory to probe the influence of toroidal
magnetic fields on the RWI. Their findings revealed a continuous reduction in
growth rates with increasing magnetization. Subsequently, Yu & Lai (2013)
explored the impact of poloidal fields using a vertically integrated disk
model. Interestingly, they observed a dichotomy in results: pure vertical
fields typically diminish RWI growth rates, whereas the presence of radial
fields tends to enhance them.

In this work, we explore the RWI unstable modes in the ideal and non-ideal MHD
limit, employing three-dimensional linear analysis with Eulerian
perturbations. The paper is organized as follows. In Section 2, we introduce
the governing dynamical equations and the perturbation equations that
delineate our theoretical framework. Section 3 elaborates on the numerical
methodologies employed to solve the set of ordinary differential equations
(ODEs) governing the magnetized RWI. In Section 4, we present numerical
solutions obtained and elucidate the observed RWI behaviors. We summarize the
main findings in Section 5 and compare the RWI linear growth to MRI in Section
6.

##  2 Theory

###  2.1 Dynamical equations

The stability of a three-dimensional, thin, magnetized disk with background
radial vortensity extrema is analyzed in cylindrical coordinates ($r,\phi,z$).
The gravitational potential is given by
$\Phi=-GM_{\star}/(r^{2}+z^{2})^{1/2}$, where $M_{\star}$ is the mass of
central star. Disk self-gravity is neglected. The governing equations for this
compressible, magnetized disk in Gaussian units are the continuity, momentum,
and entropy conservation equations,

$\frac{d\rho}{dt}+\rho\nabla\cdot\mbox{\boldmath{$v$}}=0,$ (1)
$\frac{d\mbox{\boldmath{$v$}}}{dt}+\frac{1}{\rho}\nabla\bigg{[}P+\frac{B^{2}}{8\pi}\bigg{]}+\nabla\Phi-\frac{1}{4\pi\rho}(\mbox{\boldmath{$B$}}\cdot\nabla)\mbox{\boldmath{$B$}}=0,$
(2) $\frac{dS}{dt}=0,$ (3)

where the material derivative is defined as $d/dt\equiv\partial/\partial
t+v\cdot\nabla$, and $S\equiv P/\rho^{\Gamma}$ is the entropy of the disk
matter. The induction equation is written as,

$\frac{\partial\mbox{\boldmath{$B$}}}{\partial
t}-\nabla\times(\mbox{\boldmath{$v$}}\times\mbox{\boldmath{$B$}}-c\mbox{\boldmath{$E$}}^{\prime})=0.$
(4)

The non-ideal MHD terms manifest in the electric field of the rest fluid
frame,

$\mbox{\boldmath{$E$}}^{\prime}=\frac{4\pi}{c^{2}}[\eta_{\mathrm{O}}\mbox{\boldmath{$J$}}+\eta_{\mathrm{H}}\mbox{\boldmath{$J$}}\times\mbox{\boldmath{$b$}}-\eta_{\mathrm{A}}(\mbox{\boldmath{$J$}}\times\mbox{\boldmath{$b$}})\times\mbox{\boldmath{$b$}}],$
(5)

where the unit vector of magnetic field is denoted by
$\mbox{\boldmath{$b$}}=\mbox{\boldmath{$B$}}/|B|$, and Ohmic, Hall, ambipolar
diffusivities are denoted by $\eta_{\mathrm{O}}$, $\eta_{\mathrm{H}}$, and
$\eta_{\mathrm{A}}$, respectively. The current density is
$\mbox{\boldmath{$J$}}=c\nabla\times\mbox{\boldmath{$B$}}/4\pi$. Using the
divergence free condition $\nabla\cdot\mbox{\boldmath{$B$}}=0$, the induction
equation can be cast into

$\frac{d\mbox{\boldmath{$B$}}}{dt}-(\mbox{\boldmath{$B$}}\cdot\nabla)\mbox{\boldmath{$v$}}+(\nabla\cdot\mbox{\boldmath{$v$}})\mbox{\boldmath{$B$}}+c\nabla\times\mbox{\boldmath{$E$}}^{\prime}=0.$
(6)

###  2.2 Equilibrium of the disk

The equilibrium disk model is stationary ($\partial/\partial t=0$),
axisymmetric ($\partial/\partial\phi=0$), and radially global. All background
quantities are independent of $z$. The steady-state physical quantities are
denoted by the subscript “0”. The equilibrium velocity field has only the
azimuthal component, $v_{\phi 0}=\Omega_{0}r$.

The necessary condition for Rossby wave instability is the presence of radial
vortensity extrema. This can be achieved by setting up a Gaussian bump
centered at $r=r_{0}$ in the density profile (e.g., Li et al., 2000),

$\frac{\rho_{\mathrm{0}}}{\rho_{00}}=1+(A-1)\exp[-\frac{1}{2}\bigg{(}\frac{r-r_{0}}{\Delta
r}\bigg{)}^{2}\bigg{]},$ (7)

where $\rho_{00}$ is the background density profile without the Gaussian bump.
Despite of the radially global nature of the disk model, $\rho_{00}$ is taken
to be a constant for simplicity. Setting it to a power-law distribution will
not qualitative alter the results as noted in Li et al. (2000).

To compute the background pressure $P_{0}$, we consider a barotropic flow, and
hene the pressure is only a function of density $P_{0}(\rho_{0})$, expressed
by

$\frac{P_{0}}{P_{0\ast}}=\bigg{[}\frac{\rho_{0}}{\rho_{0\ast}}\bigg{]}^{\Gamma},$
(8)

where subscript “$0\ast$” denotes background quantities evaluated at $r_{0}$,
and $\Gamma$ is the adiabatic index. The adiabatic sound speed is defined as
$c_{s0}\equiv(\Gamma P_{0}/\rho_{0})^{1/2}$. By specifying $c_{s0\ast}$, we
can obtain $P_{0\ast}$ and subsequently $P_{0}$. Throughout the paper, we set
$GM=\rho_{0}=r_{0}=1$, $\Delta r/r_{0}=0.05$, $\Gamma=5/3$, $A=1.5$, and the
disk aspect ratio $c_{s0\ast}/v_{\phi 0\ast}=0.06$ at $r_{0}$. The pressure
scale height $H\approx c_{s0\ast}/\Omega_{\mathrm{K0\ast}}$, where
$\Omega_{\mathrm{K0\ast}}$ is the Keplerian angular speed at $r_{0}$.

We construct a simplified equilibrium solution for this study, setting
$\mbox{\boldmath{$B$}}_{0}$ a constant vector. In equilibrium, the radial
momentum equation gives

$\frac{v_{\phi 0}^{2}(r,z)}{r}=\frac{1}{\rho_{0}}\frac{\partial
P_{0}(r)}{\partial r}+\frac{\partial\Phi(r,z)}{\partial r}.$ (9)

In the ideal MHD limit, the $\phi$-component of induction equation is

$\bigg{[}B_{r0}\frac{\partial}{\partial r}+B_{z0}\frac{\partial}{\partial
z}-\frac{B_{r0}}{r}\bigg{]}v_{\phi 0}=0.$ (10)

It is immediately seen that $B_{\phi 0}$ is not involved in these equations,
and hence is free to specify. By thin disk approximations,
$\Phi(r,z)\approx\Phi(r)=-GM_{\star}/r$, and $v_{\phi 0}(r,z)\approx v_{\phi
0}(r)$. Then for $B_{r0}=0$, we are completely free to specify $B_{z0}$, and
the constant equilibrium magnetic field is taken to be
$\mbox{\boldmath{$B$}}_{0}=(0,B_{\phi 0},B_{z0})$.

When considering non-ideal MHD physics to equilibrium solutions, only the
induction equation (4) is modified. The current density $J$ vanishes for a
pure $B_{z}$ field, and eq (10) is strictly satisfied. However, the existence
of curvature terms from cylindrical coordinates render non-zero $\nabla\times
c\mbox{\boldmath{$E$}}^{\prime}$ for a pure $B_{\phi}$ field. Hence, our
equilibrium magnetic field is only valid if we ignore the curvature terms in
curl operator.

The strengths of $B_{\phi 0}$ and $B_{z0}$ are parameterized by plasma
$\beta$, defined as the ratio of gas pressure to magnetic pressure,

$\beta=\frac{8\pi P_{0}}{B_{0}^{2}},$ (11)

and $\beta>1$ is commonly satisfied in protoplanetary disks. Finally, we
quantify the strength of non-ideal MHD effects by their respective Elsässer
numbers,

$\Lambda=\frac{v_{\textrm{A}}^{2}}{\eta_{\mathrm{O}}\Omega_{\text{K}}},$ (12)
$\mathrm{Ha}=\frac{v_{\mathrm{A}}^{2}}{\eta_{\mathrm{H}}\Omega_{\text{K}}},$
(13) ${\rm Am}=\frac{v_{\mathrm{A}}^{2}}{\eta_{\mathrm{A}}\Omega_{\text{K}}},$
(14)

where the Alfvén velocity is $v^{2}_{\mathrm{A}}=B_{0}^{2}/4\pi\rho$, and
$\Omega_{\text{K}}$ is the Keplerian angular speed.

###  2.3 Perturbations of the disk

Consider small perturbations to eqs (1)-(4), such that
$\mbox{\boldmath{$v$}}=\mbox{\boldmath{$v$}}_{0}+\delta\mbox{\boldmath{$v$}}(r,z,\phi,t)$,
$\mbox{\boldmath{$B$}}=\mbox{\boldmath{$B$}}_{0}+\delta\mbox{\boldmath{$B$}}(r,z,\phi,t)$,
…, . We linearize these equations by considering Eulerian perturbations
$\propto f(r)\exp(ik_{z}z+im\phi-i\omega t)$, where $k_{z}$ is the vertical
wavenumber, $m$ is the azimuthal mode number, and $\omega=\omega_{r}+i\gamma$
is the mode frequency, where $\gamma$ denotes the growth rate. We further
define the Doppler-shifted wave frequency $\Delta\omega=\omega-m\Omega$, the
azimuthal wavenumber $k_{\phi}=m/r$, and the radial epicyclic frequency
$\kappa=[r^{-3}d(r^{4}\Omega^{2})/dr]^{1/2}$.

We now drop subscript “$0$” for the background quantities throughout rest of
the paper. Our model encompasses eight perturbed quantities,
$\delta\mbox{\boldmath{$v$}},\delta\mbox{\boldmath{$B$}},\delta\rho,\delta\Psi$,
where

$\delta\Psi=\frac{\delta P}{\rho},$ (15)

and

$\frac{\partial\delta\Psi}{\partial r}=\frac{1}{\rho}\frac{\partial\delta
P}{\partial r}-\frac{1}{\rho}\frac{\partial\delta\rho}{\partial r}\delta\Psi.$
(16)

We follow Lovelace et al. (1999)’s and Li et al. (2000)’s original paper, and
define radial length scales of entropy, pressure, and density variations as

$\mathrm{L_{S}\equiv\frac{\Gamma}{d\ln S/dr}},$ (17)
$\mathrm{L_{P}\equiv\frac{\Gamma}{d\ln P/dr}},$ (18)
$\mathrm{L_{\rho}\equiv\frac{1}{d\ln\rho/dr}}.$ (19)

These length scales are related by

$\mathrm{\frac{1}{L_{P}}=\frac{1}{L_{S}}+\frac{1}{L_{\rho}}}.$ (20)

For a barotropic flow, the length scale of entropy approaches infinity,
$1/\mathrm{L_{S}}\rightarrow 0$.

To present the perturbation equations optimally, we separate them into two
cases: $k_{z}=0$ (§2.3) and $k_{z}\neq 0$ (Appendix A). We first show the set
of linearized equations in the ideal MHD limit. In the non-ideal MHD limit, we
split the equations into pure vertical and toroidal magnetic field regimes.

####  2.3.1 ideal MHD

We first derive the linearized equations in the ideal MHD limit. For
continuity equation (1) it is

$\frac{\partial\delta v_{r}}{\partial
r}+\bigg{[}\frac{1}{r}+\frac{1}{\mathrm{L_{\rho}}}\bigg{]}\delta
v_{r}+ik_{\phi}\delta v_{\phi}-i\Delta\omega\frac{\delta\Psi}{c_{s}^{2}}=0.$
(21)

The linearized momentum equations (2) are

$\displaystyle i\Delta\omega$ $\displaystyle\delta v_{r}+2\Omega\delta
v_{\phi}-\frac{\partial\delta\Psi}{\partial r}$
$\displaystyle+\frac{1}{4\pi\rho}\bigg{[}ik_{\phi}B_{\phi}\delta
B_{r}-\frac{2B_{\phi}}{r}\delta B_{\phi}-B_{z}\frac{\delta B_{z}}{\partial
r}-B_{\phi}\frac{\delta B_{\phi}}{\partial
r}+\frac{B_{\phi}^{2}}{r}\frac{\delta\rho}{\rho}\bigg{]}$ $\displaystyle=0,$
(22) $i\Delta\omega\delta v_{\phi}-\frac{\kappa^{2}}{2\Omega}\delta
v_{r}-ik_{\phi}\delta\Psi+\frac{1}{4\pi\rho}\bigg{[}\frac{B_{\phi}}{r}\delta
B_{r}-ik_{\phi}B_{z}\delta B_{z}\bigg{]}=0,$ (23) $i\Delta\omega\delta
v_{z}+\frac{ik_{\phi}B_{\phi}}{4\pi\rho}\delta B_{z}=0.$ (24)

The linearized induction equations (4) are

$i\Delta\omega\delta B_{r}+ik_{\phi}B_{\phi}\delta v_{r}=0,$ (25)
$i\Delta\omega\delta B_{\phi}+\bigg{[}\frac{\partial v_{\phi}}{\partial
r}-\frac{v_{\phi}}{r}\bigg{]}\delta B_{r}-B_{\phi}\frac{\partial\delta
v_{r}}{\partial r}=0,$ (26) $i\Delta\omega\delta B_{z}-\frac{B_{z}}{r}\delta
v_{r}-ik_{\phi}B_{z}\delta v_{\phi}+ik_{\phi}B_{\phi}\delta
v_{z}-B_{z}\frac{\partial\delta v_{r}}{\partial r}=0.$ (27)

Lastly, the entropy conservation (3) yields

$i\Delta\omega\frac{\delta
P}{P}-i\Delta\omega\Gamma\frac{\delta\rho}{\rho}-\frac{\Gamma}{\mathrm{L_{S}}}\delta
v_{r}=0.$ (28)

The barotropic assumption gives $1/\mathrm{L_{S}}\rightarrow 0$, and eq (28)
simplifies to

$\delta\Psi=c_{s}^{2}\frac{\delta\rho}{\rho}.$ (29)

It is the barotropic fluid assumption that transforms the set of perturbation
equations into standard eigenvalue problems (see §3). According to eq (15), it
is straightforward to see that the perturbed pressure is related to density by
$\delta P=c_{s}^{2}\delta\rho$.

####  2.3.2 non-ideal MHD limit: pure $B_{\phi}$

When non-ideal MHD effects apply, the continuity and entropy equations (eqs
(21) and (29)) remain unchanged, while the other equations are modified. We
now express the linearized equations for pure $B_{\phi}$. The perturbed
momentum equations are

$\displaystyle i\Delta\omega\delta v_{r}$ $\displaystyle+2\Omega\delta
v_{\phi}-\frac{\partial\delta\Psi}{\partial r}$
$\displaystyle+\frac{1}{4\pi\rho}\bigg{[}ik_{\phi}B_{\phi}\delta
B_{r}-\frac{2B_{\phi}}{r}\delta B_{\phi}-B_{\phi}\frac{\partial\delta
B_{\phi}}{\partial r}+\frac{B_{\phi}^{2}}{r}\frac{\delta\rho}{\rho}\bigg{]}$
$\displaystyle=0,$ (30) $i\Delta\omega\delta
v_{\phi}-\frac{\kappa^{2}}{2\Omega}\delta
v_{r}-ik_{\phi}\delta\Psi+\frac{1}{4\pi\rho}\bigg{[}\frac{B_{\phi}}{r}\delta
B_{r}\bigg{]}=0,$ (31) $i\Delta\omega\delta
v_{z}+\frac{ik_{\phi}B_{\phi}}{4\pi\rho}\delta B_{z}=0.$ (32)

The perturbed induction equations are

$\displaystyle i\Delta\omega\delta B_{r}$
$\displaystyle+ik_{\phi}B_{\phi}\delta v_{r}$
$\displaystyle+[\eta_{\mathrm{O}}+\eta_{\mathrm{A}}]\bigg{[}\frac{\partial^{2}\delta
B_{r}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{r}}{\partial
r}-k_{\phi}^{2}\delta B_{r}-\frac{\delta
B_{r}}{r^{2}}-\frac{2}{r}ik_{\phi}\delta B_{\phi}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}[k_{\phi}^{2}\delta B_{z}]$ $\displaystyle=0,$
(33) $\displaystyle i\Delta\omega\delta B_{\phi}$
$\displaystyle+\bigg{[}\frac{\partial v_{\phi}}{\partial
r}-\frac{v_{\phi}}{r}\bigg{]}\delta B_{r}-B_{\phi}\frac{\partial\delta
v_{r}}{\partial r}$
$\displaystyle+[\eta_{\mathrm{O}}+\eta_{\mathrm{A}}]\bigg{[}\frac{\partial^{2}\delta
B_{\phi}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{\phi}}{\partial
r}-k_{\phi}^{2}\delta B_{\phi}-\frac{\delta
B_{\phi}}{r^{2}}+\frac{2}{r}ik_{\phi}\delta B_{r}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}ik_{\phi}\bigg{(}\frac{\partial\delta
B_{z}}{\partial r}-\frac{\delta B_{z}}{r}\bigg{)}\bigg{]}$ $\displaystyle=0,$
(34) $\displaystyle i\Delta\omega\delta B_{z}$
$\displaystyle+ik_{\phi}B_{\phi}\delta v_{z}$
$\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{z}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{z}}{\partial
r}-k_{\phi}^{2}\delta B_{z}\bigg{]}$
$\displaystyle+\eta_{\mathrm{A}}\bigg{[}-k_{\phi}^{2}\delta B_{z}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}ik_{\phi}\bigg{(}-\frac{\partial\delta
B_{\phi}}{\partial r}-\frac{\delta B_{\phi}}{r}+ik_{\phi}\delta
B_{r}\bigg{)}\bigg{]}$ $\displaystyle=0.$ (35)

In the limit of pure $B_{\phi}$ and $k_{z}=0$, all three non-ideal MHD effects
are present. Ohmic resistivity and ambipolar diffusion exhibit identical terms
in both radial and azimuthal components of the induction equation. We set
$\delta B_{z}=\delta v_{z}=0$ because $\delta B_{z}$ and $\delta v_{z}$ solely
appear in eqs (32) and (35).

####  2.3.3 non-ideal MHD limit: pure $B_{z}$

For a pure $B_{z}$ field, the vertical component of the momentum equation
yields $\delta v_{z}=0$, while the remaining two linearized momentum equations
are

$\displaystyle i\Delta\omega\delta v_{r}+2\Omega\delta
v_{\phi}-\frac{\partial\delta\Psi}{\partial
r}-\frac{1}{4\pi\rho}\bigg{[}B_{z}\frac{\delta B_{z}}{\partial r}\bigg{]}=0,$
(36) $i\Delta\omega\delta v_{\phi}-\frac{\kappa^{2}}{2\Omega}\delta
v_{r}-ik_{\phi}\delta\Psi+\frac{1}{4\pi\rho}\bigg{[}-ik_{\phi}B_{z}\delta
B_{z}\bigg{]}=0.$ (37)

The non-ideal MHD effects manifest in the induction equation (4),

$\displaystyle i\Delta\omega\delta B_{r}$
$\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{r}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{r}}{\partial
r}-k_{\phi}^{2}\delta B_{r}-\frac{\delta
B_{r}}{r^{2}}-\frac{2}{r}ik_{\phi}\delta B_{\phi}\bigg{]}$ $\displaystyle=0,$
(38) $\displaystyle i\Delta\omega\delta B_{\phi}$
$\displaystyle+\bigg{[}\frac{\partial v_{\phi}}{\partial
r}-\frac{v_{\phi}}{r}\bigg{]}\delta B_{r}$
$\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{\phi}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{\phi}}{\partial
r}-k_{\phi}^{2}\delta B_{\phi}-\frac{\delta
B_{\phi}}{r^{2}}+\frac{2}{r}ik_{\phi}\delta B_{r}\bigg{]}$ $\displaystyle=0,$
(39) $\displaystyle i\Delta\omega\delta B_{z}$
$\displaystyle+B_{z}\bigg{[}\frac{1}{L_{\rho}}\delta
v_{r}-\frac{i\Delta\omega}{c_{s}^{2}}\delta\Psi\bigg{]}$
$\displaystyle+[\eta_{\mathrm{O}}+\eta_{\mathrm{A}}]\bigg{[}\frac{\partial^{2}\delta
B_{z}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{z}}{\partial
r}-k_{\phi}^{2}\delta B_{z}\bigg{]}$ $\displaystyle=0.$ (40)

In the limit of pure $B_{z}$ and $k_{z}=0$, only Ohmic resistivity and
ambipolar diffusion contribute to the linearized induction equations, but not
Hall drift. Ambipolar diffusion manifests only in the vertical component of
the induction equation, exhibiting identical terms to Ohmic resistivity.
Furthermore, $\delta B_{r}$ and $\delta B_{\phi}$ are involved only in eqs
(38) and (39), and hence can be set as $\delta B_{r}=\delta B_{\phi}=0$.

##  3 Numerical Methods

We solve the linearized equations presented in §2.3 numerically. Two numerical
methods are employed to solve ODE eigenvalue value problems: the
pseudospectral method (§3.1) and the finite difference method (§3.2).

We find that the spectral method, via the Python package dedalus, outperforms
the finite difference method. The finite difference method produced a number
of spurious modes with oscillation frequencies $w_{r}$ that closely resembled
those of RWI modes, particularly when incorporating Hall drift. Consequently,
we use the spectral method as the primary numerical approach for solving the
ODE systems and compare its results to those obtained from the finite
difference method for a selected range of parameters. The details of these two
methods are provided below.

###  3.1 Pseudospectral method

Pseudospectral methods, also known as orthogonal collocation methods,
approximate the solution of differential equations at selected collocation
points, by a weighted sum of orthogonal basis functions, which are often
chosen to be orthogonal polynomials up to a certain degree. Chebyshev
polynomials of the first kind $T_{n}$, where $n=0,1,2,...,N-1$, are chosen in
our problem, for which the eigenfunctions are expanded in. The radial domain,
that spans $r\in[0.4,1.6]$, is discretized into $N$ Chebyshev collocation
points. To minimize the interpolation error, these nodes are non-uniform and
are selected as the roots of $N$th degree Chebyshev polynomial $T_{N}$, which
cluster near the ends of the domain (Figure 1; Burns et al., 2020).

The differential equations described in §2 construct standard linear
eigenvalue problems. Written compactly in a generalized matrix form, it is

$A\vec{x}=\mathcal{L}\vec{x}+\omega\mathcal{M}\vec{x}=0,$ (41)

where $\omega$ is the eigenvalue, $\vec{x}=[\vec{\delta v}_{r},\vec{\delta
B}_{r},\vec{\delta\rho},...]^{\mathrm{T}}$ is a vector of eigenfunctions with
$M$ perturbed quantities, and $A$, $\mathcal{L}$, $\mathcal{M}$ are $MN\times
MN$ sized matrices, with $\mathcal{L}$ composed of linear operators. We employ
dedalus111https://dedalus-project.org/, a general purpose spectral code for
differential equations (Burns et al., 2020), to solve the linear eigenvalue
problem described above. We utilize the dense solver method by solve_dense in
dedalus, for which matrix $A$ is converted to dense arrays, and the
scipy.linalg.eig routine in Python is utilized to directly solve the
eigenvalue problem.

Generally, a numerical resolution of $N=256$ is adopted. Spurious solutions
can arise due to the truncation of differential equations to finite-
dimensional algebraic equations (Zebib, 1987; McFadden et al., 1990). Most
spurious modes are eliminated by increasing the numerical resolution by 1.5
times. By comparing solutions between low and high resolutions, we retain only
the valid modes within a tolerance of $10^{-6}$.

###  3.2 Finite difference method

In finite difference method, we approximate the ODEs by finite difference
equations, a technique previously adopted in RWI literature by Ono et al.
(2016); Liu & Bai (2023). We discretize the ODEs on a grid with $N=1001$
points and uniform spacing $h$, spanning the domain $r\in[0.4,1.6]$. The
system can be compactly expressed using matrices. The $MN\times MN$ sized
matrix $A$, corresponding to the set of ODEs, is written as

$A\vec{x}=\begin{vmatrix}A_{00}&...&A_{0,M-1}\\\ ...&...&...\\\
A_{M-1,0}&...&A_{M-1,M-1}\end{vmatrix}\vec{x}=0.$ (42)

Matrix $A$ consists of $M\times M$ submatrices $A_{ij}$, for $i,j=0,...,M-1$.
The value of $M$ is set by the number of perturbed quantities contained in the
system of ODEs. Each submatrix $A_{ij}$ has a size of $N\times N$. The vector
of perturbed quantities is expressed as $\vec{x}=[\vec{\delta
v}_{r},\vec{\delta B}_{r},\vec{\delta\rho},...]^{\mathrm{T}}$, where each
perturbed quantity, for example, $\delta\vec{v}_{r}=\delta\vec{v}_{r}(r_{k})$,
is evaluated at grid points $r_{k}$, for $k=0,...,N-1$.

The entries of $A_{ij}$ are filled in by constructing the linear operators of
the first derivative $d/dr$ and the second derivative $d^{2}/dr^{2}$. We adopt
the central differencing scheme to obtain the linear operators of the first
derivative (Liu & Bai, 2023),

$A_{ij,n,n-1}=-\frac{1}{2h},$ (43) $A_{ij,n,n+1}=\frac{1}{2h},$ (44)

and the second derivative,

$A_{ij,n,n-1}=\frac{1}{h^{2}},$ (45) $A_{ij,n,n}=-\frac{2}{h^{2}},$ (46)
$A_{ij,n,n+1}=\frac{1}{h^{2}},$ (47)

for $n=1,...,N-2$.

The differential equations described in §2 form standard linear eigenvalue
problems. The eigenvalue ($\omega$) and the corresponding eigenfunction
(perturbed quantities $\vec{x}$) are readily determined by utilizing
numpy.linalg.eig routine in Python.

###  3.3 Boundary conditions

To accomodate the density waves, WKBJ boundary conditions are imposed, and the
linear operator $d/dr-ik_{r}=0$ is applied to the boundary.

In the finite difference method, the boundary conditions involve in the zeroth
and $(N-1)$ th row of $A_{ij}$. Using the forward differencing scheme for the
outer boundary conditions and the backward differencing scheme for the inner
boundary conditions, the entries of $A_{ij}$ at boundaries can be written as

$A_{ij,00}=-\frac{1}{h},$ (48) $A_{ij,01}=\frac{1}{h}-ik_{r},$ (49)
$A_{ij,N-1,N-2}=-\frac{1}{h}-ik_{r},$ (50) $A_{ij,N-1,N-2}=\frac{1}{h}.$ (51)

To solve for $k_{r}$ at the boundaries, away from the RWI region where the
background quantities vary slowly over radius, we use WKBJ analysis by
assuming $\delta v,\delta B,\delta\rho...\propto\exp(ik_{r}r)$. Substituting
the derivatives by $k_{r}$, we obtain $B(k_{r},\omega)\vec{y}=0$ at the
boundaries, where $B$ is a matrix containing $k_{r}$ and $\omega$, and
$\vec{y}=[\delta v_{r},\delta B_{r},\delta\rho,...]^{\mathrm{T}}$. Both $B$
and $\vec{y}$ are evaluated at $r_{0}$ and $r_{N-1}$. We take $\omega$ to be
the value obtained from last $\beta$. This approach works if we gradually
decrease $\beta$ values in the domain of interest. Now the only unknown in
matrix $B$ is $k_{r}$, which is determined using the Newton-Raphson method by
setting the determinant of $B$ to zero, $\mathrm{det}\>B(k_{r})=0$. The
implementation of this method utilizes the source code of
scipy.optimize.newton as a reference. The iteration converges within five
steps for a tolerance of $10^{-10}$.

We choose azimuthal mode number $m=2,3,4,5$ for the calculations presented in
this work. For smaller or larger mode numbers, the growth rates can be
relatively small and are sensitive to the grid resolution. We note that when
$m=1$, the inner boundary condition differs from the WKB one described
earlier. In the hydrodynamic limit, the inner boundary lies within the
evanescent region, and the inner Lindblad resonance is absent. Instead, one
can assume a power-law behavior in $r$ and ensure regularity at $r=0$ (Ono et
al., 2016).

##  4 Results

The numerical solutions are presented for the ideal MHD limit (§4.1), Ohmic
resistivity (§4.2), ambipolar diffusion (§4.3), and the Hall drift (§4.4),
respectively. The analysis will be focused on the RWI growth rates, as the
eigenfunctions are similar to those in the pure hydrodynamic limit; two Rossby
waves are located on each side of the corotation radius, with density waves
excited at the inner and outer Lindblad resonances, propagating away from the
corotation radius (e.g., Lovelace & Romanova, 2014). Note that in the Figures
shown below, we select and present only the fastest growing mode for a
specific parameter, though a handful of slower growing modes also exist.

###  4.1 Ideal MHD

To analyze the numerical solutions in the ideal MHD limit, we begin with the
pure $B_{\phi}$ field. The top panel of Figure 1 depicts normalized RWI growth
rates $\gamma/\Omega_{\mathrm{K0}}$ as a function of plasma $\beta$, for
different azimuthal mode numbers $m=2,3,4,5$ at an infinite vertical
wavelength, $k_{z}=0$. At $\beta=100$, the growth rates nearly recover the
hydrodynamic RWI results, with the maximum growth rate observed at $m=4$ for
the chosen set of parameters. As $\beta$ exceeds $100$, the curves begin to
plateau. Conversely, as plasma $\beta$ decreases, the growth rates diminish
for all azimuthal mode numbers $m$ investigated. This suggests that increased
magnetization suppresses RWI growth. It can be interpreted as that the
magnetic fields are perfectly coupled with gas in the ideal MHD limit,
potentially hindering the free movement of perturbed gas. The top panel of
Figure 2 shows growth rates for different vertical wavenumbers
$k_{z}=0,1/10H,1/5H$ at $m=4$. It is evident that increasing vertical
wavenumbers tend to reduce RWI growth rates across the range of $\beta$ values
explored.

Figure 3 displays the RWI growth rates in the ideal MHD limit for the pure
$B_{z}$ field. Unlike in the pure $B_{\phi}$ field, where growth rates
monotonically decrease towards smaller $\beta$, different azimuthal mode
numbers $m$ exhibit distinct behaviors across the range of $\beta$ (top
panel). The parameter space of $\beta$ is now extended down to $0.01$, in
order to provide ample coverage to observe the trend in growth rates. For
relatively large $m=5$, the RWI growth rates decrease monotonically as $\beta$
decreases, similar to pure $B_{\phi}$ model. Conversely, smaller $m=2,3,4$
show an increase in growth rates as $\beta$ decreases, followed by a decline
as $\beta$ drops further. The peak of the growth rates occurs at higher
$\beta$ values for larger $m$. This is in agreement with Yu & Lai (2013), who
observed a peak in growth rate for $m=4$ around $\beta\approx 0.1$, as shown
in their Figure 2. As $\beta$ approaches $0.01$, the growth rates
asymptotically approach a stable value for each $m$ investigated. When there
is a finite vertical wavelength ($k_{z}\neq 0$), the growth rates diminish
with higher $k_{z}$ depicted in the top panel of Figure 4.

###  4.2 Ohmic resistivity

Next, we examine Ohmic resistive disks. The middle panel of Figure 1 shows RWI
growth rates as a function of the Ohmic Elsässer number $\Lambda$ at
$\beta=2$. As $\Lambda\rightarrow\infty$, the growth rates tend to approach
those of ideal MHD limit. For example, at $\beta\approx 2.212$ the growth
rates calculated for $m=2,3,4,5$ in the ideal MHD limit are
$\omega/\Omega_{\mathrm{K0}}\approx 1.974+0.0554i,\ 2.963+0.0708i,\
3.956+0.0726i,\ 4.949+0.0647i$, respectively (top panel; Figure 1).
Correspondingly, the growth rates calculated for resistive disks at
$\Lambda=100$ are $\omega/\Omega_{\mathrm{K0}}\approx 1.973+0.0528i,\
2.963+0.0666i,\ 3.956+0.0666i,\ 4.949+0.0526i$. As Ohmic Elsässer number
decreases, the growth rates for all $m$ models gradually increase. At
$\Lambda=10^{-2}$, they converge towards hydrodynamic results. For instance,
at $\beta=100$, the growth rates for $m=2,3,4,5$ in the ideal MHD limit are
$\omega/\Omega_{\mathrm{K0}}\approx 1.976+0.0892i,\ 2.966+0.118i,\
3.957+0.131i,\ 4.950+0.127i$, respectively (top panel; Figure 1). The
corresponding growth rates computed for the resistive disks at $\Lambda=0.01$
are $\omega/\Omega_{\mathrm{K0}}\approx 1.971+0.0881i,\ 2.960+0.103i,\
3.948+0.130i,\ 4.938+0.130i$, respectively.

The pure $B_{z}$ disk follows the similar trend as the pure $B_{\phi}$ disk.
As $\Lambda$ approaches infinity, the growth rates resemble those in the ideal
MHD limit. Specifically, the growth rates at $\beta\approx 2.212$ in the ideal
MHD limit are $\omega/\Omega_{\mathrm{K0}}\approx
1.973+0.0979i,2.963+0.126i,3.956+0.136i,4.949+0.128i$, respectively (top
panel; Figure 3). For comparison, at $\beta=2$ and $\Lambda=100$, the RWI
growth for $m=2,3,4,5$ are, respectively, $\omega/\Omega_{\mathrm{K0}}\approx
1.973+0.0984i,2.963+0.126i,3.956+0.136i,4.948+0.128i$. As the Ohmic Elsässer
number decreases, growth rates of different $m$ models all converge towards
hydrodynamic results. This can be readily seen by that at $\beta=100$, the
growth rates for $m=2,3,4,5$ in the ideal MHD limit are
$\omega/\Omega_{\mathrm{K0}}\approx 1.976+0.0903i,\ 2.966+0.120i,\
3.958+0.133i,\ 4.950+0.129i$, respectively (top panel; Figure 3). The
corresponding growth rates in the resistive disks at $\Lambda=0.01$ and
$\beta=2$ are $\omega/\Omega_{\mathrm{K0}}\approx 1.977+0.0878i,\
2.960+0.118i,\ 3.963+0.122i,\ 4.928+0.139i$, respectively. Among the explored
azimuthal mode numbers, only the $m=2$ curve rises with $\Lambda$, while the
others drop. This is due to the peak pattern observed for $m=3,4,5$ in the
ideal MHD limit (top panel; Figure 3).

The middle panels of Figure 2 and Figure 4 show variations in growth rates
with Ohmic Elsässer number $\Lambda$ at $m=4$ for different vertical
wavenumbers $k_{z}$. It is evident that the growth rates decrease steadily
with the vertical wavenumber $k_{z}$.

###  4.3 ambipolar diffusion

Ambipolar diffusion shares the exact perturbation equations to Ohmic
resistivity if $k_{z}=0$ (§2.3). Consequently, results obtained for
resistivity can directly apply to it. The curves obtained for Am overlap
exactly with those obtained for $\Lambda$ in the middle panels of Figure 1 and
Figure 3.

Nevertheless, for non-zero $k_{z}$ the induction equations differ between
resistivity and ambipolar diffusion. In the pure $B_{\phi}$ model, the
$\phi-$component of the induction equation shares the same terms for
resistivity and ambipolar diffusion, while the other two components differ.
Despite these differences, the growth rates in the ambipolar diffusion limit
closely approximate those in the resistivity limit when $k_{z}=1/10H,1/5H$
(middle panel; Figure 2). In the pure $B_{z}$ model, we suspect the modes
obtained are spurious and provide a detailed explanation in Appendix B.

###  4.4 The Hall drift

Unlike resistivity or ambipolar diffusion, Hall physics is sensitive to
rotation/radial shear (Latter & Kunz, 2022). This is evident by reversing the
sign of $B$ in the induction equation. In the limit of an axial field and
wavenumber, the Hall Elsässer number is written by
$\mathrm{Ha}=v_{Az}^{2}/v_{H}^{2}$, where $v_{H}^{2}\equiv\Omega B_{z}c/(2\pi
en_{e})$ is the square of the Hall velocity. It follows that $v_{H}^{2}$,
though squared, and Ha take negative (positive) values if $\Omega$ and $B_{z}$
are oriented oppositely (parallel). With more general field geometries and
wavenumbers, the sign of Ha depends on
$(\mbox{\boldmath{$\Omega$}}\cdot\mbox{\boldmath{$k$}})(\mbox{\boldmath{$B$}}\cdot\mbox{\boldmath{$k$}})$
(Balbus & Terquem, 2001). For a fixed direction of rotation, the toroidal
field reverses its sign along with vertical field. We explore positive and
negative Ha in pure $B_{\phi}$ and $B_{z}$ disk models.

In the case of pure $B_{\phi}$ and $k_{z}=0$, the Hall drift has a similar
influence as resistivity and ambipolar diffusion (bottom panel; Figure 1).
Moreover, the results obatined for Ha and -Ha overlap. As
Ha/-Ha$\rightarrow\infty$, the growth rate resembles those of ideal MHD, while
as Ha/-Ha$\rightarrow 0$, the growth rates revive and approach hydrodynamic
results. For example, at $\beta\approx 2$ for $m=2,3,4,5$ the growth rates in
the ideal MHD limit are provided in §4.2. In the Hall drift limit at
$\mathrm{Ha}=100$, the corresponding growth rates are
$\omega/\Omega_{\mathrm{K0}}\approx 1.973+0.0526i,\ 2.963+0.0675i,\
3.956+0.0683i,\ 4.949+0.0560i$, respectively. At $\beta=100$ and $m=2,3,4,5$,
the ideal MHD growth rates are also given in §4.2. In the Hall drift limit at
$\mathrm{Ha}=0.01$ and $\beta=2$, the corresponding growth rates are
$\omega/\Omega_{\mathrm{K0}}\approx 1.971+0.0896i,\ 2.960+0.121i,\
3.950+0.135i,\ 4.949+0.133i$, respectively. The bottom panel of Figure 2 shows
growth rates at $k_{z}\neq 0$. The variation of vertical wavelengths does not
significantly influence the RWI growth. The sign of Ha does not yield distinct
values of $\gamma$ when $k_{z}\neq 0$ either.

In the case of pure $B_{z}$, the Hall terms only appear in the induction
equations when the wavelength is finite ($k_{z}\neq 0$). As shown in bottom
panel of Figure 4, when vertical wave number is very small, for example
$k_{z}=1/50H$, the growth rates do not vary much across Ha/-Ha, as Hall terms
vanish at $k_{z}=0$. Increasing $k_{z}$ generically diminish the RWI growth.
Furthermore, positive and negative Ha exhibit distinct behaviors; positive Ha
results in a trough at $\mathrm{Ha}\approx 0.1$, whereas negative Ha results
in a peak at around the same location, slightly shifted towards larger |Ha|.

Figure 1: Pure $B_{\phi}$ and $k_{z}=0$. Top panel: normalized RWI growth
rates $\gamma/\Omega_{\mathrm{K0}}$ vs plasma $\beta$ in the ideal MHD limit
for different $m$. Middle panel: growth rates vs Ohmic ($\Lambda$) or
ambipolar (Am) Elsässer numbers at $\beta=2$. Bottom panel: growth rates vs
Hall Elsässer number Ha or -Ha at $\beta=2$.  Figure 2:  Pure $B_{\phi}$ and
$k_{z}\neq 0$. Top panel: normalized RWI growth rates
$\gamma/\Omega_{\mathrm{K0}}$ vs $\beta$ in the ideal MHD limit for different
$k_{z}=0,1/10H,1/5H$ at fixed $m=4$. Middle panel: growth rates vs Ohmic
($\Lambda$; solid) or ambipolar (Am; dotted) Elsässer numbers at $\beta=2$.
Bottom panel: growth rates vs Hall Elsässer number Ha or -Ha at $\beta=2$.

##  5 Summary

We studied the magnetized and lightly ionized RWI by Eulerian perturbations.
The framework is three-dimensional and radially global. Ideal MHD and non-
ideal MHD effects, including Ohmic resistivity, ambipolar diffusion, and the
Hall drift are considered. Taking advantage of the barotropic fluid
assumption, the perturbation equations form standard eigenvalue problems. The
spectral method via dedalus excels in resolving the eigenmodes. Our results
are summarized as follows.

For a pure $B_{\phi}$ field:

  * •

When gas and magnetic fields are perfectly coupled, RWI growth rates increase
with $\beta$. Strong magnetization tends to impede RWI. For all three non-
ideal MHD effects, as Elsässer numbers approach infinity, the results resemble
the ideal MHD limit; as Elsässer numbers approach zero, the results resemble
the hydrodynamic limit. In the limits of ideal MHD, Ohmic resistivity, and
ambipolar diffusion, non-zero vertical wavenumbers generically diminish RWI
growth compared to $k_{z}=0$ limit. In Hall-dominated disks, non-zero
wavenumbers do not significantly impact the results. The sign of Ha does not
yield different results.

For a pure $B_{z}$ field:

  * •

In the ideal MHD limit, RWI growth rates can either increase or decrease with
$\beta$, depending on the azimuthal mode number $m$. In the limit of
resistivity, similar to the $B_{\phi}$ model, as Elsässer numbers approach
infinity, the results resemble the ideal MHD limit. As Elsässer numbers
approach zero, the results resemble the hydrodynamic limit. In the limit of
ideal MHD and resistivity, similar to the $B_{\phi}$ model, vertical
wavenumbers generically diminish RWI growth. Hall drift only appears when
$k_{z}\neq 0$. The sign of Ha slightly complicates the growth rates.

Figure 3: Same as top and middle panels of Figure 1, but for pure $B_{z}$
disks. Figure 4: Same as Figure 2, but for pure $B_{z}$ disks. The middle
panel only shows growth rates vs Ohmic Elsässer numbers ($\Lambda$). The
bottom panel shows positive (solid) and negative (dotted) Ha at
$k_{z}=1/50H,1/10H,1/5H$, respectively.

##  6 Discussion

We compare the linear growth rates between RWI and MRI. We follow the
framework, derivations, and analyses described in Latter & Kunz (2022) for
MRI. In most cases of it, we consider channel modes with $k_{x}/k_{z}=0$, a
pure vertical field $B_{z}$ and $\beta_{z}$, Keplerian rotation, and no
vertical shear. We also consider azimuthal fields in the ambipolar diffusion
dominated regime.

We start with the simplest ideal MHD limit. Eq (29) in Latter & Kunz (2022)
shows the bi-quadratic dispersion relation, and the MRI growth rate $s$ is
found to be:

$\frac{2s^{2}}{\Omega^{2}}=-\bigg{(}\frac{2k^{2}_{z}v^{2}_{Az}}{\Omega^{2}}+1\bigg{)}+\bigg{(}\frac{16k^{2}_{z}v^{2}_{Az}}{\Omega^{2}}+1\bigg{)}^{1/2},$
(52)

where $v_{Az}=B_{z}/\sqrt{4\pi\rho}$ is the vertical Alfvén velocity. The
maximum growth rate $s=3/4\Omega$ occurs when
$k^{2}_{z}v^{2}_{Az}=15/16\Omega^{2}$ as expected. On the other hand, growth
rates vanish when $k^{2}_{z}v^{2}_{Az}=3\Omega^{2}$. In a realist
protoplanetary disk, the vertical wavelength should be shorter than the
pressure sclae height, or $k_{z}>1/H$. Therefore, the criteria for the
existence of MRI modes is $\beta_{z}>q_{R}^{-1}$, where
$q_{R}=-\partial\ln\Omega/\partial\ln R$ is the dimensionless orbital shear,
and has $q_{R}=3/2\sim O(1)$ for Keplerian rotation. The maximum growth can
occur when $\beta_{z}>32/15$. Both criteria can be readily satisfied in
protoplanetary disks, allowing linear MRI modes to easily surpass RWI modes.

For resistivity dominated disks, the instability criterion for MRI channel
modes is presented in eq (34) of Latter & Kunz (2022). Requiring the vertical
wavelengths to be longer than the thickness of the disk yields

$\beta_{z}>q_{R}^{-1}(1+\Lambda^{-2}).$ (53)

If $\Lambda\rightarrow\infty$, the system is almost in the ideal MHD regime,
and $\beta_{z}>q_{R}^{-1}$ is required for channel modes to emerge. Thus, MRI
surpasses RWI as shown above. If $\Lambda\ll 1$, the system is strongly
resistivity dominated, and $\beta_{z}\gtrsim\Lambda^{-2}$. In this regime, RWI
modes resemble hydrodynamic results, and can surpass MRI. For a zero
$B_{\phi}$ field, the instability criteria and physical behavior of ambipolar
diffusion is very similar to the Ohmic regime.

For a non-zero $B_{\phi}$ field, the MRI channel modes are diminished by small
Am, for which regime ambipolar diffusion shear instability (ADSI) can emerge
(Desch, 2004; Kunz & Balbus, 2004; Kunz, 2008). We start with the rather
simple pure MRI channel modes. These modes exist when

$\beta_{z}>q_{R}^{-1}(1+\mathrm{Am}^{-2}B^{2}/B_{z}^{2}).$ (54)

The presence of azimuthal fields contributes to stabilizing the MRI. In three-
dimensional global numerical simulations incorporating ambipolar diffusion, a
typical ratio of $B_{\phi}/B_{z}\sim O(10)$ has been observed (Cui & Bai,
2022). Using this ratio, channel modes occur under conditions where
$\beta_{z}\gtrsim 10^{4},10^{2},1$ for $\mathrm{Am}=0.1,1,10$, respectively.
Thereby, RWI may dominate over MRI for $\mathrm{Am}\lesssim 0.1$ in a
protoplanetary disk. On the other hand, ADSI becomes significant when channel
modes are diminished. ADSI modes have $|k_{x}/k_{z}|>0$. Although these modes
can always emerge for sufficiently large $|k_{x}/k_{z}|$, their growth rates
might be rather small (Latter & Kunz, 2022). To investigate the dominance
between MRI/ADSI and RWI, we compute the growth rates by the ambipolar
dispersion relation shown in eq (31) of Kunz & Balbus (2004). We utilize
$B_{\phi}/B_{z}\sim O(10)$ and set a critical growth rate threshold at
$0.01\Omega_{\text{K}}$ (see middle panels of Figure 1 and Figure 3). Under
these conditions, no modes satisfy the growth rate criterion for
$\mathrm{Am}<1$, leading us to conclude that RWI surpasses MRI/ADSI in this
regime.

Lastly, we examine the Hall dominated regime. Positive Ha exhibits fast growth
rates around $\sim O(0.75\Omega_{\text{K}})$, stemming from a blend of MRI and
Hall shear instability. There are three regimes associate with negative Ha.
For $\mathrm{Ha}<-0.25$, we are in the Hall modified and diffusive MRI
regimes, for which growth rates are rather fast $\sim
O(0.75\Omega_{\text{K}})$. For $\mathrm{Ha}>-0.25$, there is no instability
possible, rendering it ideal for RWI modes to grow.

## Acknowledgements

CC acknowledges funding from Natural Sciences and Engineering Research Council
of Canada and UK STFC grant ST/T00049X/1. AT acknowledges summer studentship
from the Centre for Mathematical Sciences, University of Cambridge. CY is
supported by the National SKA Program of China (grant 2022SKA0120101) and the
National Natural Science Foundation of China (grants 11873103 and 12373071).
MKL is supported by the National Science and Technology Council (grants
112-2112-M-001-064-, 113-2124-M-002-003-) and an Academia Sinica Career
Development Award (AS-CDA-110-M06).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Armitage (2011) Armitage P. J., 2011, ARA&A, 49, 195
  * Bai (2011) Bai X.-N., 2011, ApJ, 739, 51
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balbus & Terquem (2001) Balbus S. A., Terquem C., 2001, ApJ, 552, 235
  * Béthune & Latter (2020) Béthune W., Latter H., 2020, MNRAS, 494, 6103
  * Blaes & Balbus (1994) Blaes O. M., Balbus S. A., 1994, ApJ, 421, 163
  * Burns et al. (2020) Burns K. J., Vasil G. M., Oishi J. S., Lecoanet D., Brown B. P., 2020, Physical Review Research, 2, 023068
  * Cimerman & Rafikov (2023) Cimerman N. P., Rafikov R. R., 2023, MNRAS, 519, 208
  * Cui & Bai (2020) Cui C., Bai X.-N., 2020, ApJ, 891, 30
  * Cui & Bai (2021) Cui C., Bai X.-N., 2021, MNRAS, 507, 1106
  * Cui & Bai (2022) Cui C., Bai X.-N., 2022, MNRAS, 516, 4660
  * Cui & Lin (2021) Cui C., Lin M.-K., 2021, MNRAS, 505, 2983
  * Desch (2004) Desch S. J., 2004, ApJ, 608, 509
  * Flock et al. (2017) Flock M., Fromang S., Turner N. J., Benisty M., 2017, ApJ, 835, 230
  * Fuente et al. (2017) Fuente A., et al., 2017, ApJ, 846, L3
  * Fung & Ono (2021) Fung J., Ono T., 2021, ApJ, 922, 13
  * Galli & Shu (1993) Galli D., Shu F. H., 1993, ApJ, 417, 220
  * Girart et al. (2006) Girart J. M., Rao R., Marrone D. P., 2006, Science, 313, 812
  * Girart et al. (2009) Girart J. M., Beltrán M. T., Zhang Q., Rao R., Estalella R., 2009, Science, 324, 1408
  * Godon & Livio (1999) Godon P., Livio M., 1999, ApJ, 523, 350
  * Goldreich & Schubert (1967) Goldreich P., Schubert G., 1967, ApJ, 150, 571
  * Goldreich & Ward (1973) Goldreich P., Ward W. R., 1973, ApJ, 183, 1051
  * Güttler et al. (2010) Güttler C., Blum J., Zsom A., Ormel C. W., Dullemond C. P., 2010, A&A, 513, A56
  * Hammer et al. (2017) Hammer M., Kratter K. M., Lin M.-K., 2017, MNRAS, 466, 3533
  * Huang et al. (2018) Huang P., Isella A., Li H., Li S., Ji J., 2018, ApJ, 867, 3
  * Kunz (2008) Kunz M. W., 2008, MNRAS, 385, 1494
  * Kunz & Balbus (2004) Kunz M. W., Balbus S. A., 2004, MNRAS, 348, 355
  * Latter & Kunz (2022) Latter H. N., Kunz M. W., 2022, MNRAS, 511, 1182
  * Latter & Papaloizou (2018) Latter H. N., Papaloizou J., 2018, MNRAS, 474, 3110
  * Lesur (2021) Lesur G., 2021, Journal of Plasma Physics, 87, 205870101
  * Lesur & Papaloizou (2009) Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1
  * Li et al. (2000) Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023
  * Li et al. (2001) Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874
  * Li et al. (2020) Li Y.-P., Li H., Li S., Birnstiel T., Drążkowska J., Stammler S., 2020, ApJ, 892, L19
  * Lin (2012) Lin M.-K., 2012, MNRAS, 426, 3211
  * Lin & Hsu (2022) Lin M.-K., Hsu C.-Y., 2022, ApJ, 926, 14
  * Lin & Papaloizou (2011) Lin M.-K., Papaloizou J. C. B., 2011, MNRAS, 415, 1426
  * Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
  * Liu & Bai (2023) Liu H., Bai X.-N., 2023, MNRAS, 526, 80
  * Lovelace & Romanova (2014) Lovelace R. V. E., Romanova M. M., 2014, Fluid Dynamics Research, 46, 041401
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Lyra & Mac Low (2012) Lyra W., Mac Low M.-M., 2012, ApJ, 756, 62
  * McFadden et al. (1990) McFadden G., Murray B., Boisvert R., 1990, Journal of Computational Physics, 91, 228
  * Meheut et al. (2012) Meheut H., Keppens R., Casse F., Benz W., 2012, A&A, 542, A9
  * Miranda et al. (2016) Miranda R., Lai D., Méheut H., 2016, MNRAS, 457, 1944
  * Miranda et al. (2017) Miranda R., Li H., Li S., Jin S., 2017, ApJ, 835, 118
  * Muto et al. (2015) Muto T., et al., 2015, PASJ, 67, 122
  * Nelson et al. (2013) Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
  * Ono et al. (2016) Ono T., Muto T., Takeuchi T., Nomura H., 2016, ApJ, 823, 84
  * Pandey & Wardle (2012) Pandey B. P., Wardle M., 2012, MNRAS, 423, 222
  * Riols & Lesur (2019) Riols A., Lesur G., 2019, A&A, 625, A108
  * Rometsch et al. (2021) Rometsch T., Ziampras A., Kley W., Béthune W., 2021, A&A, 656, A130
  * Surville & Mayer (2019) Surville C., Mayer L., 2019, ApJ, 883, 176
  * Surville et al. (2016) Surville C., Mayer L., Lin D. N. C., 2016, ApJ, 831, 82
  * Tsang & Lai (2008) Tsang D., Lai D., 2008, MNRAS, 387, 446
  * Wang et al. (2019) Wang L., Bai X.-N., Goodman J., 2019, ApJ, 874, 90
  * Wardle (1999) Wardle M., 1999, MNRAS, 307, 849
  * Wardle (2007) Wardle M., 2007, Ap&SS, 311, 35
  * Weidenschilling (1977) Weidenschilling S. J., 1977, Ap&SS, 51, 153
  * Weiss et al. (2021) Weiss B. P., Bai X.-N., Fu R. R., 2021, Science Advances, 7
  * Wu et al. (2024) Wu Y., Lin M.-K., Cui C., Krapp L., Lee Y.-N., Youdin A. N., 2024, ApJ, 962, 173
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Yu & Lai (2013) Yu C., Lai D., 2013, MNRAS, 429, 2748
  * Yu & Li (2009) Yu C., Li H., 2009, ApJ, 702, 75
  * Zebib (1987) Zebib A., 1987, Journal of Computational Physics, 70, 521
  * Zhu & Stone (2014) Zhu Z., Stone J. M., 2014, ApJ, 795, 53
  * Zhu et al. (2014) Zhu Z., Stone J. M., Rafikov R. R., Bai X.-n., 2014, ApJ, 785, 122
  * Zsom et al. (2010) Zsom A., Ormel C. W., Güttler C., Blum J., Dullemond C. P., 2010, A&A, 513, A57
  * van der Marel et al. (2013) van der Marel N., et al., 2013, Science, 340, 1199
  * van der Marel et al. (2021) van der Marel N., et al., 2021, AJ, 161, 33

##  Appendix A Linearized equations for $k_{z}\neq 0$

###  A.1 ideal MHD

The linearized continuity equation is

$\frac{i\Delta\omega}{c_{s}^{2}}\delta\Psi-\frac{\partial\delta
v_{r}}{\partial
r}-\bigg{[}\frac{1}{r}+\frac{1}{\mathrm{L_{\rho}}}\bigg{]}\delta
v_{r}-ik_{\phi}\delta v_{\phi}-ik_{z}\delta v_{z}=0.$ (55)

The linearized momentum equations are

$\displaystyle i\Delta\omega\delta v_{r}$ $\displaystyle+2\Omega\delta
v_{\phi}-\frac{\partial\delta\Psi}{\partial
r}+\frac{1}{4\pi\rho}\bigg{[}i(k_{\phi}B_{\phi}+k_{z}B_{z})\delta
B_{r}-\frac{2B_{\phi}}{r}\delta B_{\phi}-B_{z}\frac{\delta B_{z}}{\partial
r}-B_{\phi}\frac{\delta B_{\phi}}{\partial
r}+\frac{B_{\phi}^{2}}{r}\frac{\delta\rho}{\rho}\bigg{]}=0,$ (56)
$\displaystyle i\Delta\omega\delta v_{\phi}$
$\displaystyle-\frac{\kappa^{2}}{2\Omega}\delta
v_{r}-ik_{\phi}\delta\Psi+\frac{1}{4\pi\rho}\bigg{[}\frac{B_{\phi}}{r}\delta
B_{r}-ik_{\phi}B_{z}\delta B_{z}+ik_{z}B_{z}\delta B_{\phi}\bigg{]}=0,$ (57)
$i\Delta\omega\delta
v_{z}-ik_{z}\delta\Psi+\frac{1}{4\pi\rho}\bigg{[}ik_{\phi}B_{\phi}\delta
B_{z}-ik_{z}B_{\phi}\delta B_{\phi}\bigg{]}=0,$ (58)

The linearized induction equations are

$i\Delta\omega\delta B_{r}+i(k_{\phi}B_{\phi}+k_{z}B_{z})\delta v_{r}=0,$ (59)
$i\Delta\omega\delta B_{\phi}+\bigg{[}\frac{\partial v_{\phi}}{\partial
r}-\frac{v_{\phi}}{r}\bigg{]}\delta B_{r}-B_{\phi}\frac{\partial\delta
v_{r}}{\partial r}+ik_{z}[B_{z}\delta v_{\phi}-B_{\phi}\delta v_{z}]=0,$ (60)
$i\Delta\omega\delta B_{z}-\frac{B_{z}}{r}\delta v_{r}-ik_{\phi}B_{z}\delta
v_{\phi}+ik_{\phi}B_{\phi}\delta v_{z}-B_{z}\frac{\partial\delta
v_{r}}{\partial r}=0.$ (61)

###  A.2 non-ideal MHD limit: pure $B_{\phi}$

The non-ideal MHD effects manifest in the induction equations,

$\displaystyle i\Delta\omega\delta B_{r}$
$\displaystyle+ik_{\phi}B_{\phi}\delta v_{r}$
$\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{r}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{r}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{r}-\frac{\delta
B_{r}}{r^{2}}-\frac{2}{r}ik_{\phi}\delta B_{\phi}\bigg{]}$
$\displaystyle+\eta_{\mathrm{A}}\bigg{[}-\frac{\delta
B_{\phi}}{r}-\frac{\partial\delta B_{\phi}}{\partial r}+ik_{\phi}\delta
B_{r}\bigg{]}ik_{\phi}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}-ik_{\phi}\delta B_{z}+ik_{z}\delta
B_{\phi}\bigg{]}ik_{\phi}$ $\displaystyle=0,$ (62) $\displaystyle
i\Delta\omega\delta B_{\phi}$ $\displaystyle+\bigg{[}\frac{\partial
v_{\phi}}{\partial r}-\frac{v_{\phi}}{r}\bigg{]}\delta
B_{r}-B_{\phi}\frac{\partial\delta v_{r}}{\partial r}-ik_{z}B_{\phi}\delta
v_{z}$
$\displaystyle+[\eta_{\mathrm{O}}+\eta_{\mathrm{A}}]\bigg{[}\frac{\partial^{2}\delta
B_{\phi}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{\phi}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{\phi}-\frac{\delta
B_{\phi}}{r^{2}}+\frac{2}{r}ik_{\phi}\delta B_{r}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}(-ik_{z}\delta
B_{r}+\frac{\partial\delta B_{z}}{\partial r})ik_{\phi}+(-ik_{\phi}\delta
B_{z}+ik_{z}\delta B_{\phi})\frac{1}{r}\bigg{]}$ $\displaystyle=0,$ (63)
$\displaystyle i\Delta\omega\delta B_{z}$
$\displaystyle+ik_{\phi}B_{\phi}\delta v_{z}$
$\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{z}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{z}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{z}\bigg{]}$
$\displaystyle+\eta_{\mathrm{A}}\bigg{[}-k_{\phi}^{2}\delta
B_{z}+k_{\phi}k_{z}\delta B_{\phi}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}-\frac{\partial\delta
B_{\phi}}{\partial r}-\frac{\delta B_{\phi}}{r}+ik_{\phi}\delta
B_{r}\bigg{]}ik_{\phi}$ $\displaystyle=0.$ (64)

###  A.3 non-ideal MHD limit: pure $B_{z}$

The linearized induction equations are written as

$\displaystyle i\Delta\omega$ $\displaystyle\delta B_{r}+ik_{z}B_{z}\delta
v_{r}$ $\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{r}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{r}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{r}-\frac{\delta
B_{r}}{r^{2}}-\frac{2}{r}ik_{\phi}\delta B_{\phi}\bigg{]}$
$\displaystyle+\eta_{\mathrm{A}}\bigg{[}-k_{z}^{2}\delta
B_{r}-ik_{z}\frac{\partial\delta B_{z}}{\partial r}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}-ik_{\phi}\delta B_{z}+ik_{z}\delta
B_{\phi}\bigg{]}ik_{z}$ $\displaystyle=0,$ (65) $\displaystyle
i\Delta\omega\delta B_{\phi}$ $\displaystyle+\bigg{[}\frac{\partial
v_{\phi}}{\partial r}-\frac{v_{\phi}}{r}\bigg{]}\delta B_{r}+ik_{z}B_{z}\delta
v_{\phi}$ $\displaystyle+\eta_{\mathrm{O}}\bigg{[}\frac{\partial^{2}\delta
B_{\phi}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{\phi}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{\phi}-\frac{\delta
B_{\phi}}{r^{2}}+\frac{2}{r}ik_{\phi}\delta B_{r}\bigg{]}$
$\displaystyle+\eta_{\mathrm{A}}\bigg{[}k_{z}k_{\phi}\delta
B_{z}-k_{z}^{2}\delta B_{\phi}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}-ik_{z}\delta
B_{r}+\frac{\partial\delta B_{z}}{\partial r}\bigg{]}ik_{z}$
$\displaystyle=0,$ (66) $\displaystyle i\Delta\omega\delta B_{z}$
$\displaystyle+B_{z}\bigg{[}\frac{1}{L_{\rho}}\delta
v_{r}-\frac{i\Delta\omega}{c_{s}^{2}}\delta\Psi\bigg{]}$
$\displaystyle+[\eta_{\mathrm{O}}+\eta_{\mathrm{A}}]\bigg{[}\frac{\partial^{2}\delta
B_{z}}{\partial r^{2}}+\frac{1}{r}\frac{\partial\delta B_{z}}{\partial
r}-(k_{\phi}^{2}+k_{z}^{2})\delta B_{z}\bigg{]}$
$\displaystyle+\eta_{\mathrm{H}}\bigg{[}-\frac{\partial\delta
B_{\phi}}{\partial r}-\frac{\delta B_{\phi}}{r}+ik_{\phi}\delta
B_{r}\bigg{]}ik_{z}$ $\displaystyle=0.$ (67)

##  Appendix B spurious modes

In a pure $B_{z}$ field, the ambipolar term in $z$-component of the induction
equation shares the same expression as the resistivity term. Unfortunately, we
found that the growth rates of RWI increase indefinitely with $k_{z}$ for
$\mathrm{Am}>0.2$. This trend persists even at a resolution of $N=512$ for
spectral method. We suspect these modes are spurious and cannot be eliminated
by high numerical resolutions. Future studies should focus on addressing this
issue.


# High temperature superconductivity: from complexity to simplicity

Sudip Chakravarty  Department of Physics and Astronomy, University of
California Los Angeles, Los Angeles, CA 90095, USA

###### Abstract

I discuss the recent quantum oscillation experiments in the underdoped high
temperature superconductors.

> HAMLET: Do you see yonder cloud that’s almost in shape of a camel ?
> POLONIUS: By th’ mass and ’tis, like a camel indeed.
> HAMLET: Methinks it is like a weasel.
> POLONIUS: It is backed like a weasel.
> HAMLET: Or like a whale.
> POLONIUS: Very like a whale.
> —William Shakespeare

More than 20 years ago, Bednorz and Müller discovered superconductivity in
copper oxides at remarkably high temperatures Bednorz:1986 . Since then,
physicists have struggled to understand the mechanisms at work. Recently, a
set of experiments on cuprates in high magnetic fields Doiron-Leyraud et al.
(2007); Bangura et al. (2007); LeBoeuf et al. (2007); Yelland et al. (2007);
Jaudet et al. (2007) has completely changed the landscape of research in high-
temperature superconductors (HTSs). In particular, the data suggest that the
current carriers are both electrons and holes, when in fact the materials are
hole doped —i.e., the current carriers should be positively charged. Moreover,
the data cannot be reconciled with an important theorem about how electrons
are organized in materials Luttinger (1960) unless one assumes that the
signals arise from a combination of both holes and electrons. Until now,
physicists have not been able to decide whether the cuprates, in Shakespeare’s
terms, are camels or whales; in fact, these experiments foreshadow a
remarkable degree of simplicity in these complex materials.

The cuprates start out as insulators and become superconductors when doped
with additional charge carriers. These so-called Mott insulators insulate by
virtue of strong repulsive Coulomb interaction and need not break any
symmetries in the lowest energy state, the ground state. A symmetry of a
system is a transformation, such as a translation or a rotation, that keeps it
unchanged. Such a symmetry is said to be broken, or spontaneously broken, if
the system does not obey the symmetry of the underlying fundamental physical
nature of the material; for example, a ferromagnet breaks the spin-rotational
symmetry with its magnetization pointing in a definite direction. The notion
of symmetry and broken symmetry finds many deep applications in physics.

Soon after the discovery of the cuprate superconductors, Anderson proposed
Anderson (1987) that their parent compounds begin as a featureless spin liquid
that does not break any symmetries, called the resonating valence bond (RVB)
state: “The preexisting magnetic singlet pairs of the insulating state become
charged superconducting pairs when the insu lator is doped sufficiently
strongly” Anderson (1987). Unfortunately, experiments show that the insulating
phase is a simple antiferromagnet in which the spins are arranged in
antiparallel manner, that is, with a broken symmetry. The materials remain
antiferromagnets for a range of doping, and then, after a sequence of not well
understood states as a function of doping, they become superconductors.

How this plays out experimentally can be understood by looking at the Fermi
surface, a fundamental concept in condensed matter physics. The Fermi surface
differentiates the occupied electronic states from the unoccupied states (in
coordinates of momentum rather than real space). Electrons fill the Fermi
surface (FS) up to some highest occupied energy called the Fermi energy (see
the figure). The excitations from the FS (e.g., when a current flows) are
called Landau quasiparticles (quasi, because they are combinations or
superpositions of real particles). The robustness of FS is due to its
topological invariance, one of the most basic invariances in mathematical
physics, which signifies stability with respect to “small deformations”
Volovik:2003 . Even when the quasiparticles are absent due to electron-
electron interactions, as in one-dimensional electronic systems, the FS is
still defined by the same topological invariant. A reconstruction of this
surface, such as a break up of a single surface into hole-like and electron-
like pockets, requires a global deformation, most likely a broken symmetry.
The new experimental work (2 6) yielded measurements of the oscillations that
arise from energy levels created by imposing a magnetic field on the material
(the Landau levels). As the magnetic field is increased, the highest fully
occupied levels sweep past the Fermi energy, and the system periodically
returns to itself, hence the oscillation in physical properties. The
oscillations of the Hall resistance Doiron-Leyraud et al. (2007); LeBoeuf et
al. (2007), capable of detecting the sign of the charge carriers, seem to show
the presence of electron and hole pockets in the Fermi surface, suggesting
that it undergoes some kind of reconstruction.

One might complain that these high field measurements are still considerably
below the upper critical field where superconductivity disappears (about 100 T
or more) and are affected by the complex motion of vortices generated by the
magnetic field. This may be true, but quantum oscillations in many
superconductors are observed at fields as small as half the critical field,
with the oscillation frequencies unchanged from the nonsuperconducting state
(with an increased damping, however). It is also known that the quasiparticles
of HTSs do not form Landau levels (9). Thus, it is very likely that the
quantum oscillation experiments are accessing the normal state beyond the
realm of superconductivity. But what kind of state? As the oscillations
definitively point to both electron and hole pockets, it cannot be a
conventional Fermi surface, rather one that has undergone a reconstruction due
to a broken symmetry at variance with the RVB picture Chakravarty and Kee
(2007).

We may be finally beginning to understand these superconductors after two
decades. The fly in the ointment is the lack of observation of electron and
hole pockets in other measurements in hole-doped superconductors (in angle
resolved photoemission spectroscopy, for instance) that are also capable of
measuring Fermi surfaces [see, however, the work on electron-doped materials
Armitage et al. (2001)]. Missing so far in experiments are also the higher
frequency oscillations that must arise from the hole pockets, not just the
electron pockets LeBoeuf et al. (2007). With further experimental work, we
should be able to tell just what kind of animal we are dealing with. [Note
added : A higher frequency is now seen and is believed to arise from an
incommensurate order—S. Sebastian, N. Harrison, and G. Lonzarich, private
communication.]

This work was supported by NSF under grant: DMR-0705092. Special thanks are
due to Hae-Young Kee for many discussions and collaborations. I would also
like to thank E. Abrahams, N. P. Armitage, R. B. Laughlin, Z. Tesanovic, and
J. Zaanen for important comments.

## References

  * (1) J. G. Bednorz and K. A. Müller, Z. Physik B 64, 189 (1986).
  * Doiron-Leyraud et al. (2007) N. Doiron-Leyraud, C. Proust, D. LeBoeuf, J. Levallois, J. B. Bonnemaison, R. X. Liang, D. A. Bonn, W. N. Hardy, and L. Taillefer, Nature 447, 565 (2007).
  * Bangura et al. (2007) A. F. Bangura, J. D. Fletcher, A. Carrington, J. Levallois, M. Nardone, B. Vignolle, P. J. Heard, N. Doiron-Leyraud, D. LeBoeuf, L. Taillefer, et al., arXiv:0707.4461 (2007).
  * LeBoeuf et al. (2007) D. LeBoeuf et al., Nature 450, 533 (2007).
  * Yelland et al. (2007) D. A. Yelland et al., arXiv:0707.0057 (2007).
  * Jaudet et al. (2007) C. Jaudet, D. Vignolles, A. Audouard, J. Levallois, D. LeBoeuf, N. Doiron-Leyraud, B. Vignolle, M. Nardone, A. Zitouni, R. Liang, et al., arXiv.org:0711.3559 (2007).
  * Luttinger (1960) J. M. Luttinger, Phys. Rev. 119, 1153 (1960).
  * Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
  * (9) G. E. Volovik, The universe in a Helium droplet (Cambridge University Press, Cambridge, 2003)
  * Franz and Tesanovic (2000) M. Franz and Z. Tesanovic, Phys. Rev. Lett. 84, 554 (2000).
  * Chakravarty and Kee (2007) S. Chakravarty and H.-Y. Kee, arXiv:0710.0608 (2007).
  * Armitage et al. (2001) N. P. Armitage, D. H. Lu, C. Kim, A. Damascelli, K. M. Shen, F. Ronning, D. L. Feng, P. Bogdanov, Z.-X. Shen, Y. Onose, et al., Phys. Rev. Lett. 87, 147003 (2001).

Figure 1: Fermi surface reconstruction: In a crystalline lattice of
periodicity $a$, the available quantum states are contained within the
Brillouin zone (BZ). In (a) we have shown a two-dimensional example,
$-\frac{\pi}{a}\leq k_{x}<\frac{\pi}{a}$, $-\frac{\pi}{a}\leq
k_{y}<\frac{\pi}{a}$, which is pertinent to high temperature superconductors
that consist of weakly coupled CuO-planes. For $N$-atoms, with unit cells
containing 1 atom in a square $a\times a$, there are exactly $2N$ available
states, as each state can be filled by electrons with two distinct spin
orientations. The diamond in (a) is called the reduced Brillouin zone (RBZ)
and contains exactly half the number of available states. In (a) the
unoccupied states are colored red; these, holes, in a otherwise filled BZ,
respond as charge carriers with positive sign in response to electric fields.
The constant energy contours are shown as the set of black curves. The filled
diamond corresponds to one electron per unit cell. The complementary sets
within the BZ, when reassembled together, will form an identical diamond,
which we may say is filled with holes. In (a) the red area corresponds to
$(1+x)$ holes per unit cell of the crystal lattice. The excess, $x$, is called
the doped holes. High temperature superconductors are schizophrenic. In some
regimes they behave as though they consist of $(1+x)$ charge carriers, holes,
while in the underdoped regime their properties are determined instead by $x$
doped holes. This is an important mystery. A class of theories posit that FS
reconstructs in the uderdoped regime. Consider shifting the FS in (a) by
vectors $(\pm\frac{\pi}{a},\pm\frac{\pi}{a})$, which will give rise to the
Figure (b), ignoring the shading for clarity. Interesting quantum mechanical
processes, about which there can be much debate, can result in reconnections
shown in (c), as in a kaleidoscope. However, if we continued to consider the
full BZ, we would double the number of states. All distinct states are
contained in the RBZ, but there are now two distinct sets of energy levels,
the upper band and the lower band. However, we continue to use the full BZ in
(c), as a better aid for visualization. Because the RBZ is the fundamental
unit in the wave vector space, the new unit cell of the crystal lattice is
doubled, given by a square $\sqrt{2}a\times\sqrt{2}a$, and the full
translational symmetry of the original lattice is broken. The FS now consists
of disconnected sheets of blue and red areas. The remarkable fact is that the
charge carriers in the blue region behave like electrons of fraction $n_{e}$
and in the red region like holes of fraction $n_{h}$. The doped holes is
easily seen to be $x=2n_{h}-n_{e}$, as there are two hole pockets and one
electron pocket in the RBZ. The broken symmetry invoked here is called
commensurate, as the translational invariance of the crystal of integer
multiples of the next nearest neighbor lattice vectors of the original lattice
is still preserved. The broken symmetry can also be incommensurate with the
original crystal lattice and can give rise to more complex FS reconstructions.