Parametric Survey of Nonaxisymmetric Accretion Disk Instabilities: Magnetorotational Instability to Super-Alfvénic Rotational Instability

Accretion disks are highly unstable to magnetic instabilities driven by shear flow, where classically, the axisymmetric, weak-field magnetorotational instability (MRI) has received much attention through local WKB approximations. In contrast, discrete nonaxisymmetric counterparts require a more involved analysis through a full global approach to deal with the influence of the nearby magnetohydrodynamic (MHD) continua. Recently, rigorous MHD spectroscopy identified a new type of ultralocalized, nonaxisymmetric instability in global disks with super-Alfvénic flow. These super-Alfvénic rotational instabilities (SARIs) fill vast unstable regions in the complex eigenfrequency plane with (near eigen)modes that corotate at the local Doppler velocity and are radially localized between Alfvénic resonances. Unlike discrete modes, they are utterly insensitive to the radial disk boundaries. In this work, we independently confirm the existence of these unprecedented modes using our novel spectral MHD code Legolas, reproducing and extending our earlier study with detailed eigenspectra and eigenfunctions. We calculate the growth rates of SARIs and MRI in a variety of disk equilibria, highlighting the impact of field strength and orientation, and find correspondence with analytical predictions for thin, weakly magnetized disks. We show that nonaxisymmetric modes can significantly extend instability regimes at high mode numbers, with maximal growth rates comparable to those of the MRI. Furthermore, we explicitly show a region filled with quasi-modes whose eigenfunctions are extremely localized in all directions. These modes must be ubiquitous in accretion disks, and play a role in local shearing box simulations. Finally, we revisit recent dispersion relations in the appendix, highlighting their relation to our global framework.


Introduction
Accretion disks are generally assumed to be in a turbulent state, as observed accretion rates are higher than expected from angular momentum transport through molecular viscosity alone (Pringle 1981;Balbus & Hawley 1998).The origin of this turbulence and the resulting enhanced viscosity has long been unclear: hydrodynamic turbulence is largely excluded in typical accretion disk flow profiles that satisfy the Rayleigh criterion for linear stability.Nonaxisymmetric hydrodynamic instability pathways, both linear and nonlinear, for differentially rotating tori were identified early on (Papaloizou & Pringle 1984;Blaes & Hawley 1988), but these more global modes are generally not expected to lead to fully developed turbulence in disks.
Only some 30 years ago, a convincing linear route to turbulence in disks was discovered through the weak-field, axisymmetric magnetorotational instability (MRI; Balbus & Hawley 1991), revisiting earlier analysis by Velikhov (1959) and Chandrasekhar (1960).Even a seed vertical magnetic field was shown to dramatically alter stability by acting as a tether between fluid elements, transporting angular momentum outward.Nonlinear simulations confirmed how, at unstable MRI wavelengths, the quickly developing field distortions lead to angular momentum interchanges (Hawley & Balbus 1991).Hawley & Balbus (1992) extended these nonlinear 2D simulations in a shearing sheet to address their long-term evolution and found magnetic energy growth and sustained angular momentum transport.The linearly unstable, quickly growing MHD eigenmodes induce the generation of turbulent eddies at multiple scales, whose interactions may explain the enhanced effective viscosity through the stress tensor.In the meantime, MRI has become an accepted pathway to disk turbulence through linear MHD instability in many fully 3D nonlinear MHD simulations, from fully global Newtonian (Mishra et al. 2020) and general relativistic settings (Begelman et al. 2022;Ripperda et al. 2022) to local 3D shearing boxes in ideal or resistive MHD (Lesaffre et al. 2009;Bai & Stone 2013;Hirai et al. 2018) or even kinetic particle-in-cell settings (Bacchini et al. 2022(Bacchini et al. , 2024)).
The linear dispersion relation for the standard weak-field MRI is obtained from a WKB analysis in a Cartesian shearing sheet or box, where the background curvature is ignored (Balbus & Hawley 1991).This local approach is in perfect correspondence with results from global cylindrical studies of the MRI where the full background variations are taken into account (Blokland et al. 2005;Latter et al. 2015).However, the linear MHD spectrum of eigenoscillations is intricately organized and truly sensitive to background variations, as exposed by MHD spectroscopy, which has its theoretical foundations in laboratory plasma contexts (Goedbloed et al. 2019).Early examples of MHD spectroscopy in cylindrical models relate to rotating screw pinches (Hameiri 1981), but it has also been applied to accretion disk theory (Keppens et al. 2002) to obtain the general equations governing all linear eigenmodes of any gravitating, rotating, magnetized disk with radial variation.
In contrast to the axisymmetric MRI, nonaxisymmetric modes fundamentally require such a global spectroscopic analysis because of singularities in the governing equations produced by the overlap of the continuous MHD spectra for radially varying, rotating equilibrium states.These nonaxisymmetric instabilities are of particular interest because they are related to dynamo activity, which is impossible under axisymmetric conditions (Moffatt 1978), and because their coupling to the toroidal field is essential for self-sustained instability (Hawley et al. 1996;Kitchatinov & Rüdiger 2010;Begelman et al. 2022).Furthermore, although global, they can become very localized with growth rates as high as those of the axisymmetric MRI, and they may still be unstable under conditions where the MRI is typically stabilized (Begelman et al. 2022;Goedbloed & Keppens 2022), as is shown in this work.Additionally, secondary parasitic nonaxisymmetric instabilities developing on channel flow solutions in shearing boxes have growth rates that can even exceed that of the MRI (Goodman & Xu 1994).Although the original MRI description was immediately followed by a local WKB description of the nonaxisymmetric MRI featuring transient growth (Balbus & Hawley 1992), a systematic MHD spectroscopic approach has shown the existence of unstable modes that cannot be found from local theory alone (Matsumoto & Tajima 1995;Ogilvie & Pringle 1996;Goedbloed & Keppens 2022).Additionally, growth rates for these "local" nonaxisymmetric MRIs are an order of magnitude underestimated as compared to those of their "global" versions (Curry & Pudritz 1996).
Research on the global linear treatment of nonaxisymmetric modes peaked in the decade after the MRI was discovered: Matsumoto & Tajima (1995) performed a global analysis of a local shearing sheet and found evanescent modes localized between two Alfvén resonances; Curry & Pudritz (1996) investigated a cylindrical equilibrium in the Boussinesq approximation with a constant vertical field and found two corotating Alfvénic modes, generalizing the hydrodynamic Papaloizou-Pringle instability (Papaloizou & Pringle 1984); Terquem & Papaloizou (1996) performed linear analysis and simulations of a cylindrical disk having both radial and vertical structure and a toroidal field, and in their simulations obtained unstable modes growing on increasingly smaller scales; Ogilvie & Pringle (1996) compared nonaxisymmetric modes in both a cylinder and a shearing sheet for a purely toroidal field, and found two branches of modes clustering toward the continua in the limit k → +∞ with arbitrarily localized eigenfunctions; Noguchi et al. (2000) further explored the evanescent eigenmodes found by Matsumoto & Tajima (1995), although the range of unstable fields in their shearing sheet is at odds with the earlier cylindrical models; and Keppens et al. (2002) gave a fully general theoretical and numerical description of the MHD eigenspectrum of a cylindrical model and found an abundance of unstable nonaxisymmetric modes.This general formalism has yet to be applied to perform a comprehensive linear study of nonaxisymmetric accretion disk instabilities, which is what we aim to accomplish in this work.In contrast to many of the above works, we will consider an arbitrary magnetic field orientation and strength to show how predictions from these limiting cases (purely vertical and azimuthal fields) hold up in general.
In recent years, nonaxisymmetric mode analysis has been sporadically revived, mostly in local shearing box models where the radial wavenumber changes in time and the instability is transient (Kitchatinov & Rüdiger 2010;Mamatsashvili et al. 2013) or in WKB and global treatments of cylindrical models with a dominant azimuthal field (Das et al. 2018;Begelman et al. 2022) as it becomes clear that these modes must play an important role in the dynamics of accretion disks.However, in nonlinear simulations, the development of turbulence is still often attributed solely to the axisymmetric MRI.Concurrently, there is claimed laboratory evidence for nonaxisymmetric MRI in an azimuthal field (Seilmayer et al. 2014), for the standard vertical field MRI (Wang et al. 2022a), and for unidentified nonaxisymmetric modes (Wang et al. 2022b), to which linear global analysis has also been applied (Ebrahimi & Pharr 2022).All this proves that these are still exciting times for linear studies of nonaxisymmetric instabilities, contrary to popular belief.
Indeed, Goedbloed & Keppens (2022) recently found a completely new type of linear ideal MHD instability in thin, weakly magnetized, cylindrical accretion disks that could provide a broadly accessible pathway to turbulence.These modes densely populate a 2D region in the complex eigenfrequency plane and have eigenfunctions that resemble radially localized wave packages.The instability is termed super-Alfvénic rotational instability (SARI), which reflects the fact that a disk has a rotational Doppler frequency mv θ /r dominating over the Alfvén frequency if it is assumed to be thin and weakly magnetized.These modes are promising because of the apparent ease with which they can be excited (they are not normal modes but require only a tiny addition of energy to be brought into resonance), their insensitivity to boundary conditions (the eigenfunctions are radially bound by Alfvén resonances), and the property that they exist all over the disk (each perturbation peaks at the local Doppler corotation radius).SARIs exist because of a delicate interplay between the MHD continua and the (near) singularities they introduce in the analysis of the governing ordinary differential equations (ODEs); hence they can only be obtained through a global analysis.Goedbloed & Keppens (2022) provided full analytic as well as numerical evidence for SARIs, using the powerful Spectral Web approach (Goedbloed 2018a(Goedbloed , 2018b) ) to locate complex eigenfrequencies of ideal, stationary MHD states.In addition, they rigorously generalized a criterion for instability and growth rate prediction.In this work, we aim to confirm the existence of SARIs with an independent numerical approach and to show the applicability of this instability criterion.
To distinguish between the transient growth of local nonaxisymmetric MRI as described in the literature (Balbus & Hawley 1992;Kitchatinov & Rüdiger 2010;Mamatsashvili et al. 2013) and the global eigenmode nature of the nonaxisymmetric SARIs, we will in the remainder of this work refer to any axisymmetric mode as the MRI and to nonaxisymmetric modes as SARIs in the regime of thin, weakly magnetized disks.We address some properties of nonaxisymmetric modes that are known from local theory or limiting cases in global analysis but are underexposed in the general global approach.In particular, we confirm analytical predictions for instability over a range of equilibria and show that the Alfvén frequency plays a central role, which opens up more regimes to instability.Furthermore, we show that global high-m modes are nevertheless very localized, with possible implications for shearing box simulations.
In this work, we first revisit, augment, and further validate the results of Goedbloed & Keppens (2022) using the Legolas code (Claes et al. 2020).The equilibrium setup and numerical approach are described in Section 2. Section 3 summarizes the properties of the novel quasi-continuum modes, as well as those of the standard MRI and discrete SARI modes in a global approach.It will become clear that resonances with the Alfvén continua play a crucial role in the behavior of nonaxisymmetric modes.Second, we perform a comprehensive study of the most unstable MRI and SARI growth rates in a range of equilibria.Sections 5.1-5.3compare growth rates for the MRI and SARIs and show that nonaxisymmetric modes can significantly extend regimes of instability to stronger fields.These results are shown to match classical and new analytical predictions in Section 5.4, and a comparison with local dispersion relation solutions highlights the necessity of a global approach.Finally, we again focus on the full MHD eigenspectra in Section 5.5 to explicitly show a quasicontinuum at a high azimuthal wavenumber, featuring an abundance of modes with significant growth rates that are global in nature but radially and vertically extremely localized.As additional material, Appendix A summarizes the mathematical tools required for a global and local analysis.Appendix B compares the general local dispersion relation of Blokland et al. (2005) and other recent dispersion relations that have been used for nonaxisymmetric modes (Das et al. 2018;Ebrahimi & Pharr 2022).

Model Equilibrium
Our basic model representing a cylindrical accretion disk is an ideal MHD equilibrium (Keppens et al. 2002;Blokland et al. 2005;Goedbloed & Keppens 2022) that radially varies according to self-similar profiles (Spruit et al. 1987) and that is vertically invariant.The disk is modeled as an annulus with radius r ä [r 1 , r 2 ], with a line source approximating the central gravitating object, The assumption of vertical invariance in the disk is valid if small vertical wavelengths (or large wavenumbers k) are considered, much smaller than the scale height H = c S /(Ω 1 r 1 ) = r 1 in a thin, Keplerian disk (Shakura & Sunyaev 1973), where g r = c p S is the sound speed and Ω 1 is the disk rotation frequency.From here on, we eliminate all proper units from the variables through normalizing by the equilibrium values at the inner boundary r 1 , all denoted by a subscript 1.The normalization is fixed by choosing new units for density ρ 1 , distance r 1 , and velocity = * v GM r K1 1 (the inner Keplerian rotation velocity).The nondimensionalized radial equilibrium profiles for the density, pressure, velocity, and magnetic field are then given by such that the stationary (∂/∂t = 0) MHD force balance equation is satisfied, Here, Λ = 0 if the flow is Keplerian, with modifications by pressure gradients and magnetic tension.The dimensionless radial domain is r ä [1, 1 + δ], where the parameter δ = (r 2 − r 1 )/r 1 describes the size of the domain.It is instructive to introduce three additional parameters that determine the constants at the inner disk edge in Equation (1) (Goedbloed & Keppens 2022): magnetic field strength or Alfvén speed , and inverse pitch angle μ 1 = B θ1 /B z1 .It should be noted that the general inverse pitch angle μ := B θ /(rB z ) is actually radially dependent, but that the plasma beta is uniform throughout the disk as p and B 2 have the same r variation.The latter also implies that b 2 is constant, so the varying sound speed is a fixed multiple of the varying Alfvén speed.We take γ = 5/3 for a monatomic gas.
In terms of these free parameters, which fully determine an equilibrium, the constants at the inner boundary become where the rotation frequency is obtained from the radial force balance equation.Note that, although the rotation power law is Keplerian, the actual rotation frequency Ω 1 is in general not Keplerian because of pressure and tension forces.For convenience, we denote j = q B B arctan z 1 1 for the angle between the equilibrium magnetic field and the z-axis.Variation of the field orientation from poloidal (vertical) to toroidal (azimuthal) is hence achieved through varying j from 0°to 90°.

Eigenvalue Problem
The force-balanced equilibrium given above can then be analyzed for its eigenoscillations, i.e., the complete set of complex eigenfrequencies and eigenfunctions of linear perturbations.Perturbed quantities f are assumed to be of the form corresponding to a Fourier decomposition in the spatial and temporal domains with complex Fourier coefficients ˆ( ) f r , real azimuthal and vertical wavenumbers m and k, respectively, and complex eigenfrequencies ω = σ + iν.Due to the θ-periodicity of the cylindrical equilibrium, m is quantized to integer values, whereas the vertical wavenumber k can take arbitrary values in this case.The set of eigenoscillations of this equilibrium is fully determined with the wavenumbers m and k.
We calculate all MHD eigenspectra in this work using the Legolas3 code developed by Claes et al. (2020) andDe Jonghe et al. (2022) and recently upgraded in terms of solver and memory performance (Claes & Keppens 2023).It uses a finite-element approach to solve the following generalized matrix eigenvalue problem obtained from the linearized MHD equations, where the Fourier coefficients become eigenfunctions: Here, X is an eight-component eigenvector before its Galerkin decomposition (dropping the hat on the Fourier coefficients), containing eigenfunctions (ρ, v r , v θ , v z , T, a r , a θ , a z ), with ρ denoting the density perturbation, T the plasma temperature perturbation, and a the magnetic vector potential, which fixes the perturbed magnetic field B = ∇ × a and enforces the monopole constraint ∇ • B exactly.Legolas constructs the matrices A and B from the equilibrium variations, the chosen finite-element base functions in the eigenmode representations, and the required physical effects like the viscosity, resistivity, Hall MHD, etc. (De Jonghe et al. 2022).As this work is limited to ideal MHD, the matrices A and B are modified from the static, linear MHD case with only flow and gravity as additional effects.The matrix system is then solved using various available solvers depending on the requirements.In this work, we use the QR-invert algorithm to produce complete eigenspectra, and the Arnoldi shift-invert method to search for eigenvalues in the neighborhood of some complex number (Claes & Keppens 2023).The mathematical formulation of the eigenvalue problem is closed by fixing the boundary conditions at radius 1 and 1 + δ.By default, Legolas assumes perfectly conducting, solid walls as boundaries, requiring v r = 0 and a θ = a z = 0 at those locations.As explained in Goedbloed & Keppens (2022), these solid-wall boundaries will affect the detailed variation of MRIs, while SARIs are fully agnostic of these boundary conditions, an important distinction between both mode categories.
To substantiate our Legolas results, we combine them in several figures with results obtained with the complementary ROC code, which employs a shooting method (Goedbloed 2018a(Goedbloed , 2018b)).Appendix A gives a brief description of the formulation of the spectral problem for stationary, cylindrical equilibria in terms of a second-order ODE, which forms the foundation for the ROC code and for the analytical results of Goedbloed & Keppens (2022).Key to our further discussion of the SARI are the MHD continua, which are singular frequency ranges that organize the spectrum.Because of the rotational equilibrium flow, all frequencies are Doppler-shifted according to the Doppler range: Note that, because v z = 0 in our setup, the Doppler shift vanishes for axisymmetric (m = 0) perturbations.The four singular ranges are dubbed forward and backward Alfvén and slow continua, as the frequencies lie in a band on the real axis.The continua are given by which are the Doppler-shifted Alfvén and slow continua.From here on, W  A,S denotes the closed range of continuum frequencies expressed by in Goedbloed & Keppens (2022).In accretion disks with subthermal magnetic fields, these Alfvén and slow frequency ranges almost coincide because of incompressibility (if γ → +∞, they exactly coincide).Our numerical results do always account for compressibility effects, as we adopt a ratio of specific heats γ = 5/3 for the disk plasma.Whereas ±ω A are always separated for m = 0 modes like the MRI, the continua can overlap for nonaxisymmetric modes.Exactly this phenomenon is responsible for many of the interesting properties of SARIs.For large values of m, Ω 0 dominates over ω A and the forward and backward continua overlap to a large degree.

Accretion Disk Instabilities
By means of the tools above, we can now perform an MHD spectroscopic study of accretion disk instabilities.We begin by revisiting the spectra typically associated with MRI and SARI as given by Goedbloed & Keppens (2022) to calibrate our results.We find exact correspondence in terms of eigenvalues and eigenfunctions for both MRI and discrete SARIs.
The MHD spectrum of a typical weakly magnetized, near-Keplerian disk for axisymmetric m = 0 modes contains the MRI, among other modes.The structure of such a spectrum is well known (Keppens et al. 2002;Blokland et al. 2005;Das et al. 2018;Goedbloed & Keppens 2022), and it is shown in Figure 1(a) for a weak field with equal poloidal and toroidal components and k = 70.Vertically, we find a finite sequence of overstable discrete MRI modes (and their damped counterparts), and horizontally, we find oscillating Alfvénic/slow modes clustering to either edge of the Alfvén continua, indicated by red and black segments on the real axis.In fact, the oscillating modes clustering to the edges ±ω A2 belong to the sequence of MRI modes when listed according to the increasing number of zeros of their eigenfunctions but they are stable (in a local MRI analog, the radial wavenumber q becomes larger, stabilizing the MRI).Two sequences of fast modes cluster toward infinity, but they are not shown in the figure.The up-down symmetry of the spectrum is a standard feature of ideal MHD where all operators are self-adjoint (Goedbloed et al. 2019), but the damped MRI variant is of course of limited physical significance for accretion disk turbulence (but is significant in the initial value problem).The spectrum is, however, not left-right symmetric as a result of background flow.In particular, the axisymmetric MRI is not exponentially unstable (note the slight deviation from the imaginary axis), contrary to local WKB predictions (Balbus & Hawley 1991), although there exists a local dispersion relation that does reproduce this behavior (Blokland et al. 2005).A typical eigenfunction, here the real part of v r for the most unstable eigenvalue, is shown in the inset (the imaginary part looks similar).The perturbation is localized at the inner boundary and has no zeros inside the domain.Moving down through the sequence of MRIs, the number of zeros increases and the perturbations become global over the entire domain.They are hence sensitive to the boundary conditions.
Although a complete spectroscopic picture of nonaxisymmetric accretion disk instabilities has been missing until recently, it has been known for some time that a wavenumber m ≠ 0 produces an additional branch of discrete unstable modes if the Alfvén continua overlap.These modes resonate with both Alfvén continua and corotate with the disk at the local Doppler frequency (Matsumoto & Tajima 1995;Curry & Pudritz 1996;Ogilvie & Pringle 1996).Since a weakly magnetized, thin disk rotates super-Alfvénically (v A1 = 1 ≈ Ω 1 ), Goedbloed & Keppens (2022) termed them super-Alfvénic rotational instabilities in their general analysis, where they showed that the two branches consist of infinitely many discrete modes clustering toward the internal edges of the Doppler-shifted continua W + A and W - A .These modes are all unstable, in contrast with the finite number of unstable MRIs in  A (red and black) overlap in this super-Alfvénic regime.For these counterrotating SARIs (since m < 0), both Ω 0 and W  A are decreasing, and the inner SARI eigenfrequencies cluster to W -( ) r A 1 whereas the outer SARIs cluster to W + ( ) r A 2 (and vice versa for corotating modes with m > 0).The most unstable SARIs are indicated in the figure, and we note that the most unstable inner SARI has the highest growth rate -a general feature of disks with outwardly decreasing equilibrium profiles (for an explanation, refer to Section 5.4).The fact that both branches have oscillation frequencies σ inside the Doppler range (gray) implies that these modes corotate with (or against) the local flow at the Doppler corotation radius r * , where Ω 0 (r * ) = σ.Since σ is also in one of the continua W  A , each mode has one Alfvénic resonance, which has direct bearing on the eigenfunctions shown in the two insets of panel (b), with vertical lines denoting the resonances: the inner v r eigenfunction (orange) is localized between the inner boundary and the resonance with the backward Alfvén continuum, the outer eigenfunction (green) between the resonance with the forward Alfvén continuum and the outer boundary.Both peak at the location of resonance with the local Doppler frequency.Clearly, the names "inner" and "outer" also reflect the localization of the eigenfunctions (Ogilvie & Pringle 1996).Skin currents at the resonant location shield the mode from one boundary, acting as a virtual wall (Goedbloed & Keppens 2022).The number of zeros of the eigenfunctions (the "radial wavenumber") again increases down the sequence, but it is less straightforward to number them accordingly (Goedbloed & Keppens 2022).

Radially Extended Disks
Even though the radial domain in Figure 1(b) is chosen to be small, the SARI eigenfunctions are not global because they are limited by an Alfvén resonance.In some sense, each mode actually has a second Alfvén resonance outside of the domain at a fictitious radius if the continua are extended.Intuitively, these additional resonances are encapsulated when the domain is enlarged, and the result might be hybrid modes that have properties of both inner and outer SARIs: they are localized between the two resonances, and hence shielded from the boundaries.It is clear that for a real accretion disk, these boundaries are unphysical.Physically relevant modes should hence preferentially be insensitive to these artificial constraints, like the modes in radially extended disks proposed here.
In reality, the matter is more complicated, but indeed, overlapping forward-and backward-shifted Alfvén continua can produce modes with two resonances, forming an unprecedented 2D region in the complex eigenfrequency plane of overstable quasi-modes, as recently reported by Goedbloed & Keppens (2022).In regimes with a large overlap between W  A (if δ is not small, or k and m are large), a wealth of modes appears above the Doppler range.These modes have a nonzero complementary energy W com , a complex-valued quantity that should vanish for true normal modes (Goedbloed 2018a(Goedbloed , 2018b; see Appendix A).This complementary energy quantifies the addition of energy required to excite these modes; they are almost eigenoscillations like normal modes, but not exactly.However, for these modes, W com is very low, sometimes even below machine precision, rendering the resulting near eigenmodes numerically indistinguishable from actual normal modes -hence Goedbloed & Keppens (2022) termed them quasimodes.Allowing a low but finite W com causes dense regions to appear in the complex plane containing such modes, from here on named continua of quasi-modes or quasi-continua.In practice, nonzero W com implies that the total pressure perturbation shows a minute jump at the location where a shooting method links up two numerical solutions starting from the left and right boundaries, but this is numerically virtually unnoticeable at the small values of W com = 10 −12 that we have here.Physically, W com is the additional energy required to bring the mode into resonance with the time dependency w - . Such tiny excitation energies as for the quasi- modes in this work could easily be provided by ambient energy fluctuations in turbulent disks, so that there is no observable difference from a normal mode.Moreover, in actual nonlinear numerical MHD simulations, one always has roundoff and truncation errors, related to the spatiotemporal discretizations used.Deviations from machine precision can be sufficient to introduce nonzero complementary energy at any time during the run.
An unstable quasi-continuum has oscillation frequencies inside the Doppler range and near the Alfvén continua.An analytic treatment of the SARI quasi-modes must hence deal with three near-singular resonance locations (Goedbloed & Keppens 2022), which are only accessible through global analysis.In Goedbloed & Keppens (2022), analytic and shooting-based numerical results were presented using the Spectral Web technique (Goedbloed 2018a(Goedbloed , 2018b)), emphasizing the distinctive properties of these quasi-continuum SARIs.We revisit these results in the next section using a different and independent approach with Legolas.

Quasi-continuum Modes
The above discussion on quasi-modes and their low but finite complementary energy represents an unusual numerical challenge to a finite-element-based solver like Legolas.Because of the Galerkin decomposition of eigenmodes, the boundary conditions and eigenfunction continuity (even differentiability) will be satisfied per construction, and will not show a minute jump associated with nonzero W com .If a quasi-mode is found, it is a priori truly indistinguishable from any other mode in the spectrum.However, a change in resolution will always slightly shift the quasi-mode eigenvalues found previously, whereas convergence toward a solution generally requires that an eigenvalue/eigenvector pair change minimally when increasing the resolution.This regular convergence behavior is indeed true for discrete modes like MRI or the inner and outer, discrete SARI variants.
Despite these difficulties, it is possible to validate the quasicontinuum results from Goedbloed & Keppens (2022) with Legolas with the combined effort of multiple solvers implemented in the code.Figure 2 overlays some of our results on Figure 15 of Goedbloed & Keppens (2022) for the same equilibrium.Panel (a) shows the m = 1 co-propagating SARI spectrum found using the QR-invert algorithm for increasing resolution n grid ä {150, 200, 300, 500}.The spectra are overplotted on the W com 10 −12 (black) and 10 −8 regions calculated with the ROC code and show a good match: the QRinvert algorithm locates the left and right edges of the quasicontinuum, as well as additionally finding the discrete inner and outer SARI modes clustering to the continua and disappearing in a sea of quasi-modes.For increasing grid resolution, the upper edge of the quasi-continuum as calculated by Legolas remains more or less stationary, but the lower edge shifts downward, thereby increasingly resolving the lateral edges.Since the independent Spectral Web approach can identify the entire black shaded region as having a minute complementary energy, we must conclude that the lower edge of the quasi-continuum obtained by the QR algorithm is hence an artificial feature of the limited resolution.The correspondence between the quasi-continuum edges and the black region for W com = 10 −12 (small compared to the normalized eigenfunctions) indicates that Legolas operates on a similar accuracy.It should be noted that low-resolution versions of Figure 1(b) also show a quasi-continuum-like structure between the two SARI branches, but it becomes smaller as the resolution is increased and its upper edge approaches the real axis.On the other hand, the upper edges of true quasicontinua quickly converge to their final values when the resolution is increased.We use this observation to define the first criterion for the detection of quasi-continua with Legolas: a QR solver locates the edges of the quasicontinuum, and with increasing resolution, the upper edge converges whereas the lateral edges are increasingly resolved.
To show that the regions delimited by the QR-invert solver actually contain infinitely many unstable modes, we use the Arnoldi shift-invert solver (ARPACK; Lehoucq et al. 1998), which finds a predefined number of eigenvalues closest to some complex target frequency ω t .We define a grid of target shift frequencies that sample the 2D region of interest, where we use the quasi-continuum edges found by the QR algorithm.Thanks to the recent upgrades to the Legolas code (Claes & Keppens 2023), this can now easily be done at high radial resolutions: we choose a resolution of n = 1000 gridpoints for the equilibria and look for 100 eigenvalues in the neighborhood of each shift.depending on the location, resolution, and number of eigenvalues.However, it is clear that the annuli are only artificial: the eigenvalue clouds found by nearby shifts overlap, allowing for the 2D quasi-continua to be filled completely (as in Figure 3).This observation can be used as a second criterion for defining quasi-continua with Legolas: when applying shift-invert within the bounds determined by the QR algorithm, arbitrarily many modes can be located in the direct neighborhood of the shift.Panels A is large, as it is for high-m modes, the separation is very small and the modes are hence radially very localized (but still global in nature).As for the two quasi-modes in panel (d), the orange and green eigenfunctions only resonate with the forward (red) and backward (black) Alfvén continua, respectively.In this sense, modes like these are the quasi-continuum equivalent of discrete SARI modes, where the second Alfvén resonance can now be either inside the domain or not.An inspection of the other eigenfunctions (not shown) reveals that the Lagrangian perturbations along and perpendicular to the magnetic field are the largest, highlighting the slow/Alfvénic nature of the instability.
The huge potential of quasi-mode SARIs to govern instability in accretion disks becomes clear when plotting composite spectra for several wavenumbers, where the Doppler range moves for every m.In Figure 3, we show spectra for modes with wavenumbers m ä {−3, −2, K, 3} and wavenumbers k/|m| = 50 for the SARIs and k = 50 for the MRI, recreating Figure 20 of Goedbloed & Keppens (2022) using Legolas.We use the combined efforts of the QR algorithm to find the quasi-continuum edges and the Arnoldi shift-invert algorithm to fill the 2D regions.The color code adopted here will be used in the remainder of this paper-the m = 0 MRI in black and the m < 0 counterrotating and m > 0 corotating SARIs in red and blue, respectively-to give the Doppler shift its usual connotation.We find an excellent correspondence with the W com 10 −12 regions delimited by ROC.The Alfvén continua (red and black) and Doppler range (gray) are indicated for each value of m, and it is clear that each quasi-continuum is limited to the Doppler range.The furthest edges of the m = ±3 continua fall out of this frame at around σ = −3.5 and 3.8, respectively.With Legolas, we additionally obtain the discrete SARI modes clustering toward the continua, as well as stable Alfvénic modes (visible on the real axis for m = ±1).Note again the up-down symmetry of the spectrum, a feature of ideal MHD (Goedbloed et al. 2019).This figure underlines that a typical equilibrium has a huge range of unstable perturbations over many wavenumbers, with perturbations that can be radially localized anywhere in the disk.If a linear MHD view suffices to explain the cause of turbulence, this abundance of entire 2D regions of near-eigenmode SARIs that grow exponentially and experience differential rotation in the disk is a much more convincing argument than the existence of a few isolated unstable MRIs.

SARI and MRI in Accretion Disk Equilibria
The previous section made it clear that allowing for modes with m ≠ 0 has important implications for accretion disk stability.We independently confirmed the existence of continuous regions of quasi-modes using Legolas, and pointed out how the axisymmetric MRI differs from its nonaxisymmetric counterparts.In this section, we set out to systematically explore the stability of inner/outer SARIs and MRI in various disk equilibria, thereby comparing the fastestgrowing modes over a range of wavenumbers m.In Section 5.1, we first look at the influence of the vertical wavenumber on MRI/SARI growth.We then explore the influence of the magnetic field orientation (Section 5.2), and field strength and plasma beta (Section 5.3).These results are then compared to predictions from analytical theory in Section 5.4.In Section 5.5, we use this knowledge to explicitly obtain a quasi-continuum of extremely localized modes.
Before exploring the instabilities further numerically, it is useful to state the MRI/SARI instability criterion and analytical maximum growth rate predictions derived by Goedbloed & Keppens (2022).This upper bound for growth, obtained from the global problem (their Equation (84)), improves on an upper bound derived under the WKB approximation (Balbus & Hawley 1992): which yields the well-known instability criterion (Balbus & Hawley 1991, 1998;Ogilvie & Pringle 1996;Blokland et al. 2005), Here, the squared epicyclic frequency is W 2 for our Keplerian rotation profile.Note that the Alfvén contains both wavenumbers k and m, and it can now go through zero and become negative for m < 0, in contrast to the situation for the MRI.If ω A vanishes and the lower bound is violated, the modes should stabilize.The same should happen if the upper bound is violated because the field gets too strong or the wavenumbers too large-physically, magnetic tension then stabilizes modes with wavelengths that are sufficiently small.Equation (9) can also be differentiated with respect to w A 2 to obtain the maximum growth rate over all wavenumbers at a fixed radius: which exactly matches the well-known local MRI prediction (Balbus & Hawley 1998) but for a completely general field orientation.We will compare the above predictions to our numerical results in Section 5.4.The dimensionless parameters ò, β, μ 1 , and δ defined in Section 2.1 completely determine an ideal MHD equilibrium with cylindrical flow and gravity, representing an accretion disk.Goedbloed & Keppens (2022) take v A1 = B 1 = ò = 1 and p 1 ∼ ò 2 β = 1, and k 2 r 2 ?m 2 so that the equilibrium state represents a near-Keplerian (Ω 1 ≈ 1), weakly magnetized, thin accretion disk (vertical scale height 2 ; Pringle 1981).Super-Alfvénic rotation then follows, since v A1 = 1 ∼ mΩ 1 .These assumptions allow for an analytical treatment of the global eigenvalue problem, but a numerical code like Legolas is not limited to those regimes.Nevertheless, we will find that in most cases, the analytical growth rate predictions are still sufficient to understand the stability of the modes (Section 5.4).
The values of the parameters are fixed to one of the two following sets throughout this section: 1. Case (a): ò = 0.0141, μ 1 = 1, β = 100, δ = 1, and k = 50, corresponding to a weak field and pitch angle j = 45°( see Section 2.1); or 2. Case (b): ò = 0.045, μ 1 = 10, β = 10, δ = 1, and k = 50, corresponding to a slightly stronger, predominantly toroidal field (j ≈ 85°), where one or more parameters are varied.Both sets of parameters have a sizable B θ component-indeed, there is a growing set of evidence that significant toroidal fields may be generated in accretion disks (Das et al. 2018;Begelman et al. 2022).The main differences between the two cases are then an equipartition between the poloidal and toroidal fields versus a dominant poloidal field, and a weak field versus a slightly less weak field.In both cases, Ω 1 ≈ 0.987.
Legolas works in dimensionless units, which are fixed in all runs by the following default choices: B u = 10 G for the magnetic field, L u = 10 9 cm as the length unit, and T u = 10 6 K for the temperature.The spectra of Figures 1, 2, 3, and 11 are given in these dimensionless values.The growth rates in Figures 4-8 are additionally normalized by the inner disk rotation frequency Ω 1 .All spectra in those figures are calculated at a resolution of 200 base gridpoints.

Vertical Localization
It is well known from local linear theory that the MRI has a cutoff value for the vertical wavenumber k (Balbus & Hawley 1991, 1998): the strength of the poloidal field component determines the smallest vertical wavelength not stabilized by magnetic tension.We first investigate if this generally applies to SARIs in a global eigenfunction description, and at which vertical wavelengths the growth is maximal by quantifying the most unstable modes from numerical spectra obtained with Legolas.
Figure 4 shows how the growth rates vary with k for Case (b) (j = 85°), which is more illustrative here than Case (a).The growth rates of the m = 0 MRI are shown in black, and those of the m > 0 and m < 0 SARIs in blue and red, respectively, with the inner SARIs in solid lines and the outer branch dotted.One immediate observation is that all instabilities have a clear cutoff wavelength and reach maximum growth for intermediate wavelengths.Both the cutoff wavenumber k cutoff and the wavenumber of maximal growth k max are ordered by m and significantly increased for counterrotating m < 0 SARIs, which allow for much smaller unstable vertical wavelengths, as was also found by Begelman et al. (2022).Maximum growth rates are comparable for the MRI and inner SARI, but the outer SARIs consistently have lower growth rates.The latter growth rate curves largely follow those of the MRI and inner SARIs, but are scaled down by a factor of 2 horizontally and a factor of around 3 vertically.The reason for this scaling will be explained in Section 5.4.Strikingly, for m < 0 SARIs, there is an intermediate value of k for which the growth rates vanish before the maximum is reached.This is most prominent for m = −10, where the inner branch is stabilized at k = 100 and the outer branch at k = 50.At these wavenumbers, the Alfvén frequency vanishes locally, i.e., ω A (r * ) = 0, at the corotation location r * , which violates the SARI instability criterion given by Equation (10).This highlights the fact that the magnetic field is essential for instability, and instability does not prevail for those wavenumbers where the wave is directed perpendicularly to the field (k • B = 0).It is important to bear in mind that the range of allowed k is limited in a realistic, 2D disk with a vertical extent and hence a quantization on k or at least an upper limit for vertical wavelengths.Hence, the range of k where the modes are physical is determined by H ∼ c S , which gives in this case k  30.Case (a) (not shown) has much smaller values of k max and k cutoff because of the stronger vertical field component.The variation with m is similar to Case (b) but less pronounced: because the field has a larger vertical component (B z1 = B θ1 ), the influence of the azimuthal wavenumber m is reduced and one needs to go to more negative m to significantly increase k cutoff .

Field Orientation
Most analytical treatments of accretion disk instabilities consider the limiting cases of a dominant poloidal background field (Balbus & Hawley 1991, 1998;Mamatsashvili et al. 2013) or a dominant toroidal background field (Ogilvie & Pringle 1996;Terquem & Papaloizou 1996;Das et al. 2018).However, a global spectroscopic approach can treat any completely general equilibrium with radial variation (Keppens et al. 2002;Blokland et al. 2005), where in our case the orientation and strength of the magnetic field can be smoothly varied by the parameters μ 1 and ò.In what follows, we consider inverse pitch angles j = q B B : arctan z 1 1 between 0°and 90°, where the former corresponds to a purely poloidal (vertical) field and the latter to a purely toroidal (azimuthal) field, and calculate SARI and MRI growth rates at a fixed vertical wavenumber k = 50.
Figure 5(a) shows the growth rates for Case (a) (weak field and β = 100).The MRI and low-m SARI are the fastestgrowing modes whenever there is a significant B z component present, with inner and outer SARI growth rates decreasing with |m|.All growth rates drop significantly as the field becomes more toroidal, but the SARI growth rates remain finite whereas the MRI is stabilized completely as j → 90°: with B z → 0 and m = 0 there is no way for the instability to couple to the toroidal component, and ω A = 0 everywhere (Balbus & Hawley 1992).For a purely toroidal field, the growth rates increase with |m|.This is understood from Figure 4: an equivalent figure for fixed k and varying m would show an equivalent growth rate curve for m with an m max and m cutoff .Hence, for values of |m| below m max , the growth rates increase with |m|.As in Figure 4, for every m < 0 there is a value of j at which ω A (r * ) = 0 near the inner or outer boundary and as a result the respective inner or outer mode that resonates with this frequency is locally stabilized.as m becomes increasingly negative: the relative strength of the poloidal field allowed for instability increases.Comparing with Case (a), we conclude that for ò = 0.014, the vertical component is able to support instability at k = 50, but ò = 0.045 is too strong unless m is very negative.The cutoff's being less sharp for m = ±10 SARIs is a resolution effect.For a purely toroidal field, the MRI is again suppressed and SARIs reach comparably high growth rates, again increasing with |m| in this range of azimuthal wavenumbers.Again, m < 0 SARIs can be stabilized locally.The maximum growth rates of all instabilities are comparable to those of Case (a), but the stronger azimuthal field increases the SARI growth rates in the limit j → 90°.
Ogilvie & Pringle (1996) performed a global analysis of a cylindrical accretion disk model with a purely toroidal background field.They found that instability prevails for arbitrary values of the wavenumbers k and m and that the fastest growth happens on the smallest scales, implying that there is no cutoff vertical wavelength k.Terquem & Papaloizou (1996) also came to this conclusion when finding modes growing on increasingly localized scales, eventually limited only by dissipative scales.Figure 6 the growth rate variation with k for six field orientations and fixed m = 1.If the poloidal component dominates (j 45°), the k cutoff varies only slightly below k = 200.As the relative strength of the toroidal component increases, k cutoff moves further away to values above k = 1000, and eventually to arbitrarily large values as j → 90°.A similar observation can be made for k max .This confirms that indeed, for all practical purposes the vertical scale of the fastest-growing modes becomes arbitrarily small for a toroidal field.However, this extreme localization immediately breaks down when the field gets a small poloidal component, as was also noted by Balbus & Hawley (1992).Additionally, the growth rates remain finite in the limit of a purely toroidal field even for low k because of the nonaxisymmetric coupling to the toroidal component.For an axisymmetric mode like the MRI, this is not the case.

Magnetization
Two of the essential assumptions for SARI modes (and for the standard MRI in Balbus & Hawley 1991) are a weak field, i.e., v A1 = 1, and a thin disk with scale height , so that near-Keplerian disks rotate super-Alfvénically and the background is close to incompressible.We now investigate fields from weak to dynamically relevant, modeling disks with near-Keplerian to strongly sub-Keplerian rotation speeds.
Figure 7 shows the MRI and SARI growth rates as a function of Alfvén speed v A1 = ò normalized by the inner rotation speed v θ1 = r 1 Ω 1 , with the most unstable inner and outer SARIs again denoted by solid and dotted lines, respectively.Panels  Figure 6.Growth rates of the m = 1 inner SARI peak and cutoff at increasingly large k when the field becomes nearly toroidal but this extreme localization breaks down for even a small poloidal component.However, growth rates remain finite for small k in the limit j → 90°.
corresponds roughly with that of the MRI.This pattern will be explained from the SARI instability criterion in Section 5.4.Note that in both cases the outer SARIs again behave very similarly to the inner SARIs, but with lower growth rates and a stricter cutoff.In Case (b), the outer m = −10 SARI is stabilized, as k = 50 is the vertical wavenumber where ω A (r 2 ) = 0 (this is also clear from Figure 4).Note that the normalization factor Ω 1 is itself dependent on ò (Section 2.1).This parameter can become much smaller than its Keplerian value Ω 1 = 1 for the strongest fields in Case (b), going down to Ω 1 < 0.3.For such sub-Keplerian disks, the inner disk is magnetically supported instead of rotationally.
The field strength can also be considered relative to the gas pressure through varying the plasma beta.Although the standard MRI was first obtained for disks threaded by a weak, subthermal field (Balbus & Hawley 1991), it has been shown to prevail in suprathermal disks as well (Das et al. 2018), motivated by simulations where a strong suprathermal toroidal field is generated (β ∼ 0.1-1; Bai & Stone 2013).Figure 8(a) shows how the MRI and SARI growth rates depend on the plasma beta for our Case (a) with equal poloidal and toroidal field components.The MRI is the most unstable mode, closely followed by the inner SARI, with growth rates decreasing with |m|.The ±m SARI growth rates almost coincide for β < 1 but are separated for subthermal fields, which becomes more apparent as |m| increases.For these equilibria, suprathermal magnetic fields (β > 1) produce slightly higher growth rates than subthermal fields (β < 1).The outer SARIs show a similar trend, but the ±m growth rates almost coincide for subthermal fields.They diverge, however, in the limit β → 0, where the corotating modes retain a finite growth rate but the counterrotating modes have a vanishing growth rate.Note that the m = −10 outer mode is stabilized completely for β → 0, while the m = 10 outer mode is briefly stabilized and then reappears.However, a closer inspection of the spectra for these equilibria (not shown) reveals that these modes are not SARIs: the slow continua are separated from the Alfvén continua and have their proper overlap, and the outer eigenfunctions are not confined to the region between the resonances with the Alfvén continuum.They could be related to the hybrid Alfvén/slow modes described in Das et al. (2018).Figure 8(b) shows quite a different behavior for Case (b) with a slightly stronger dominant toroidal field.All modes show a large variation in growth rates and the MRI is not the most unstable mode over most of the β range.All growth rates become low for β → 0. Note the sudden drop in growth of the m = 10 modes at large values of β, while other growth rates drastically increase.This is a result of the normalization, where Ω 1 depends critically on β.When the disk becomes strongly sub-Keplerian, Ω 1 = 1.It is expected that for values of β where Ω 1 → 0, the growth rates  will remain finite and suddenly vanish like those for the m = 10 modes.In contrast to most of our results shown above, the growth rates of the outer m = −3 SARI are slightly higher than those of the inner SARI over almost the entire parameter range.Indeed, the wavenumber k = 50 is close to that where the m = −3 inner mode is locally stabilized (Figure 4), and hence it has a reduced growth rate.The m = −10 outer mode is again fully stabilized at this wavenumber.

Analytical Predictions
We now revert back to the MRI/SARI instability criterion and maximum growth rate to explain the results of the previous sections and test the analytical predictions.

Instability Criterion
Much of the SARI behavior can be explained by the instability criterion (Equation ( 10)).Compared to the case of the m = 0 MRI, the toroidal field plays a significant role in the instability of nonaxisymmetric modes.An important consequence is that the perturbations can couple to the toroidal field: for any given value of k, there exists a negative m that makes ω A low enough to satisfy the criterion.A second consequence is that for certain wavenumber combinations, ω A (r * ) = 0 locally, hence violating the lower end of the instability criterion.Both of these possibilities are completely absent for the MRI.
A direct consequence of the first observation is a shift in the cutoff vertical wavenumber with m as found in Figure 4. From Equation (10), it is clear that the cutoff vertical wavenumber shifts to larger values compared to the axisymmetric MRI for m < 0 and to smaller values for m > 0. For a sizable poloidal field, this effect is small: for Case (a), which is not shown, all cutoffs are in the interval [150,175] for |m| 10.For a dominant poloidal field, the effect is more noticeable: according to the analytical predictions, k max and k cutoff shift to large values when m becomes arbitrarily negative.Furthermore, the wavenumber k where ω A vanishes locally shifts to larger values as m becomes increasingly negative.For increasingly positive m, the range of allowed k decreases due to stabilization by the increasing magnetic tension.
The validity of the instability criterion is made strikingly clear when we plot the growth rates as a function of ω A (r 1 ), as in Figure 9.It contains the same results of Figure 4, where k is varied, but without the outer SARIs.Stability is indeed approximately limited to the region where w < W ( ) r 3 A 1 1 and it is clear how the m < 0 SARIs have access to the whole unstable range, and are stabilized when the Alfvén frequency moves through zero.Note how the m > 0 SARIs cannot access the whole range of unstable Alfvén frequencies.The maximum growth rates are indeed approximately reached for w = ( ) r A 1 W 15 16 1 , as predicted by Equation (11), but this value shifts depending on m.The theoretical maximum growth of 0.75Ω 1 is, however, not attained.The slight differences for w Amax and the maximum growth rate depending on m are partly explained by the fact that these inner SARI growth rates are evaluated not at the exact Doppler radius where their eigenfunctions peak, but at the inner radius instead.Finally, note that the equilibrium for Figure 9 has v A1 ∼ ò = 1 ∼ r 1 Ω 1 and hence is in the super-Alfvénic regime.This is not in contradiction with Figure 9: indeed, although instability extends to wavenumbers k where ω A1 > Ω 1 , the rotation frequency in the local Doppler frame ω A1 < mΩ 1 is still super-Alfvénic.
The outer SARI growth rates (not shown) obey a similar dependency on ω A (r 2 ) when scaled by Ω(r 2 ).Note also that the condition w 2 for outer SARI instability is more restrictive than that for the inner SARIs, as the right-hand side of the inequality decreases faster than the left-hand side.This explains why outer SARI modes are stabilized more easily compared to inner SARIs.We also noted in the previous sections that outer SARIs have growth rates three times lower than those of inner SARIs.When the outer SARIs are, however, scaled by the local rotation frequency ≈Ω(r 2 ) instead of Ω 1 , this discrepancy disappears.Indeed, Relative to the local frame of reference, the outer modes are hence equally fast-growing over one rotation period, with maximum growth rates up to 0.7Ω in our numerical results.The fact that the outer SARIs are more easily stabilized than the inner modes gives two possibilities for the evolution of the entire spectrum when an equilibrium parameter is varied: (1) The quasi-continuum first collapses to a point as the overlap between the Alfvén continua decreases, where it might detach from the real axis.The inner and outer branches then merge into one branch of discrete modes, or the outer branch is stabilized before that happens.(2) The outer branch is stabilized before the quasicontinuum disappears, and the quasi-continuum still exists as long as the overlap of the continua is large enough.In either case, the remaining branch of inner modes is reminiscent of the unstable sequence of MRI modes, but it contains an infinite number of modes clustering to an internal edge of one of the continua W  A (Goedbloed & Keppens 2022).It, too, eventually disappears to the real axis when it is stabilized.
To compare the predictions to the actual values of k cutoff in Figures 4-8, we evaluate Equation (10) at the local Doppler radius r * , where σ = Ω 0 (r * ), and the MRI growth rates at the inner boundary r 1 .These predictions are reasonably accurate when compared to the actual growth rates of Figure 4, e.g., the m = −3 SARI is predicted to have k cutoff = 411 while in reality k cutoff = 400, and the value of k where the inner mode is locally stabilized is predicted to be k stab = 33 while in reality it is k stab = 30.When keeping k fixed in Figures 5 and 7, it is the orientation or strength of the magnetic field that determines whether Equation (10) is satisfied.The Alfvén frequency in terms of the equilibrium parameters is given by w Hence, the instability criterion predicts that, once the field gets stronger, large values of k are only allowed when the toroidal component is significant; this is exactly the difference between Figures 5(a) and (b), which both have k = 50.The instability criterion again explains the ordering of the cutoff vertical field obtained for Case (b): if m < 0, B θ can become smaller compared to B z , and vice versa for m > 0. For Case (a), there is no cutoff as k is below the critical value k cutoff for a purely poloidal field.For a fixed orientation and wavenumber, instability is only possible in a range of field strengths where the upper limit of Equation (10) is satisfied.The ordering of m < 0 and m > 0 mode critical field strengths in Figures 7(a) and (b) is again explained from this criterion, except for the high v A1 in Case (b).Interestingly, Equation (10) leaves the option for MRI/SARI in a range of wavenumbers to be "switched off" at intermediate pitch angles j, since ω A as a function of j might first increase above the critical threshold of Equation (10) and then again decrease.This implies that some high-m corotating SARIs are stabilized for j ä [5°, 25°], for example, which means that they are unstable for an almost purely vertical field or when a significant azimuthal field is present.Finally, a visualization of ω A as μ 1 or ò is varied gives almost the same results as those in Figure 9.

Growth Rate Estimates
We now investigate the upper bound to the SARI growth rates by evaluating Equation (9) at the local Doppler corotation radius, where σ = Ω 0 (r * ), which is also where the eigenfunctions peak.Expression (9) for the maximal growth rate varies with r and is usually radially decreasing for our accretion disk equilibria.This implies that inner SARIs (localized around r 1 ) indeed have higher growth rates than outer SARIs (localized around r 2 ), as is found in most of our spectra and parameter scans.There are, however, specific instances where n max is nondecreasing and the outer SARI can have higher growth rates.Recall for example the cases with vanishing ω A (r 1 ) above, or the m = −3 SARI in Figure 8.In the second case, n max increases from its value at r 1 up to a global maximum, after which it again decreases to its value in r 2 , which is larger than that in r 1 .However, there is no mode corresponding to the maximum of n max with a growth rate higher than that of the outer mode, so the outer SARI is the most unstable mode.This also points to a second observation: the predicted upper bound for growth rates is never strict.It usually overestimates the growth rates, for both discrete and quasi-continuum modes, but one can nevertheless derive the shape of the quasi-continuum upper edge from it, as in Figure 10.A possible reason for this discrepancy is that the derivative of the eigenfunctions, c¢, is not taken into account in the derivation of Equation (9), as noted in Goedbloed & Keppens (2022).We will now use Equations ( 9) and (10) to explain several observations made earlier.
To fully confirm the finding of Figure 6 that a purely toroidal field allows for arbitrarily large unstable wavenumbers k, we extend the numerical results with predictions using Equation (9).Indeed, as j → 90°, the values of k max and of k cutoff increase dramatically: the growth rates for j = 89°.99 are predicted to peak at k = 2 × 10 5 , and those for j = 88°.999 at k = 1.5 × 10 6 ; the most unstable modes can become arbitrarily localized in the vertical direction.Here, the distinction with the m = 0 MRI becomes clear: we perform the same Legolas runs as in Figure 6 but now for the MRI, extended with analytical predictions for very large k.Unsurprisingly, the m = 0 MRI also shows increased localization as j → 90°b ecause the poloidal component gets weaker and hence larger k are allowed.However, whereas the SARI retains finite growth rates for small k as j → 90°, the MRI is stabilized in the limit: k max moves to arbitrarily large values, but for values of k 1000, the growth rates become arbitrarily low.In contrast to the SARI modes, the m = 0 MRI has no way to couple to the toroidal field, and once the poloidal component becomes negligible, it is stabilized.
Expression (9) for n max also explains why ±m SARIs have coalescing growth rates for j → 0°or j → 90°in Figure 5 for both Case (a) and Case (b).Indeed, for a purely toroidal field, n max only depends on m 2 and for a purely poloidal field, m completely disappears from the expression (but evaluating n max at the proper resonance location r * for each mode does produce the observed dependence on m).Physically, a purely vertical or toroidal field induces no preferential (coor counterrotating) direction for waves to grow in.As mentioned in Section 5.3, Figure 7(b) shows a peculiar ordering of the critical field strengths v A1cutoff , which decrease as m increases from m = −3 to 10, but the m = −10 SARI has a smaller range of unstable field strengths.From Equation (9), we find the following behavior of the instability as m varies: for m > 0, the range of unstable field strengths decreases and disappears in the limit m → ∞, since then ω A becomes high.For m < 0, the critical field strength first increases to a maximum, then decreases and eventually vanishes too in the limit m → −∞.When the field has a dominant toroidal component, this switchback occurs for lower m than when the field is dominantly poloidal.
We conclude that, in general, Equation (9) correctly predicts the general trend of growth rates in weakly magnetized, thin disks, the parameter regime for which it is derived.Outside of this regime, it breaks down: when the assumption k 2 r 2 ?m 2 is violated in Figure 4 all growth rates nearly vanish as k → 0 but the predicted growth rates remain finite; for β < 1 in Figure 8, growth rates are not at all correctly predicted, and the m = ±10 SARIs behave in a distinct way.It is clear that in the compressible regime, one has to resort to numerical means as done in this work.Equation (9) does not replicate the decreasing maximum growth rates as |m| increases, which is observed above in most figures.

Global versus Local Approach
The global eigenvalue problem as in Keppens et al. (2002) has often been in the shadow of the local WKB approximation because of the involved methods required for its analysis.However, the interaction of nearby singularities of the governing ODE is exactly what is responsible for the special properties of nonaxisymmetric modes.We briefly show in this section that local dispersion relations are unable to correctly predict both discrete and quasi-continuum SARI modes.In Appendix B, we show how such a dispersion relation can be derived from the global formulation using a WKB approach (Keppens et al. 2002;Blokland et al. 2005), and how some recent dispersion relations obtained directly from the linearized MHD equations also follow from the global formulation (Das et al. 2018;Ebrahimi & Pharr 2022).
Solving a local WKB dispersion relation requires a radius r where all equilibrium parameters should be evaluated, and a radial wavenumber q in addition to k and m.Blokland et al. (2005) showed that a sixth-order local dispersion relation can exactly replicate MRI eigenfrequencies if r and q are derived self-consistently from the peak radius and radial wavelength of an exact eigenfunction of the global formulation.We now evaluate this dispersion relation over the whole radial domain to obtain the growth rates of the most unstable modes, and consider several values of q.In Figure 10, a range of predicted eigenfrequencies from the sixth-order dispersion relation (Equation (B4)) is plotted for three fixed values of q (orange).For the most unstable inner SARI, this leads to a good approximation (q ≈ 50 is the value obtained from the eigenfunction).The same cannot be said for the quasicontinuum SARI modes: since the eigenfunctions for modes with growth rate ν are very similar in shape but shifted horizontally according to the local Doppler resonance, a constant value of q would have to correspond with a nearly horizontal band of quasi-continuum modes.However, as Figure 10 shows, the maximal growth rates are quickly decreasing toward the outer edge of the disk and the dispersion relation predicts a stabilization of all modes for all three radial wavenumbers.Note that the larger wavenumbers q, which should correspond to radially localized modes, perform worse in predicting ranges of unstable modes, and even start missing out the entire Doppler range.Indeed, also in the local dispersion relation of Terquem & Papaloizou (1996), the MRI is stabilized for q ?k.The WKB analysis assumes radially oscillating solutions over the whole domain, whereas SARI eigenfunctions are radially localized.MRI eigenfunctions can fill the entire radial domain and hence are better suited for such an approximation.In particular, the dispersion relation gives no correct predictions for the outer SARI growth rates.The local WKB approximation is hence unable to replicate the complex interactions of the two singularities, which are crucial in the analysis of SARI modes.In contrast, the green line in Figure 10, representing n max (Equation ( 9)) obtained through global means, provides an upper bound for the discrete SARIs as well as for the quasi-continuum modes, and correctly predicts the radial decrease of the maximum growth rates.

Extremely Localized Modes in Global Spectra
In Sections 3 and 4, we show typical MRI and SARI spectra (Figure 1), as well as unstable 2D continua filled with quasimodes (Figure 3).We then proceed in Section 5 to calculate the growth rates of the most unstable MRI and SARIs over a range of parameters, showing how the inner SARIs usually dominate over the outer SARIs and how modes can be stabilized by locally vanishing ω A .However, not just the fastest-growing modes are of interest in this study: the presence of quasicontinua provides a source of unstable modes with localized eigenfunctions all over the disk, as shown in Section 4. We set forth to calculate a quasi-continuum of large-m SARIs that are extremely localized (but global in nature).
In principle, a SARI quasi-continuum can occur whenever W + A and W - A overlap, but the height that it might reach can be severely limited, requiring a high resolution (n grid 1000) to adequately resolve its edges.This depends on the ratio |k|/|m|, which appears in the exponent of exact quasi-continuum solutions (Goedbloed & Keppens 2022).As this ratio increases, the associated complementary energy greatly decreases and as a result, the quasi-continuum "grows."These two conditions are summarized in Equation (91) of Goedbloed & Keppens (2022), which we repeat here: where the upper bound makes sure that the Alfvén continua overlap for the self-similar radial profiles in this setup.
We use Equation (12) and the results of Section 5 to obtain a SARI quasi-continuum for a high azimuthal mode number m = −10.The two conditions need to be balanced: the ratio |k/m| must be high, so a k ?|m| is needed-for Figure 3, k/|m| = 50-but at the same time the instability has a cutoff with k when the magnetic tension of the vertical field component stabilizes the mode, as seen in Section 5.1.Additionally, the overlap between W  A strongly depends on ò −1 , which needs to be large.A high-m quasi-continuum hence prevails in weak fields, where large k are allowed for instability.Alternatively, a stronger toroidal component moves k cutoff to larger values as shown in Section 5.2 and thus can also allow large k.An additional difficulty is that, for strongly overlapping W  A , SARI resonances are close by and the modes are radially very localized between those resonances, requiring a high resolution for resolving the eigenfunctions.All in all, we use the equilibrium parameters μ 1 = 10, β = 10, ò = 0.01, δ = 1, and k = 400 so that k/|m| = 40.
Figure 11 shows the results of four QR runs with increasing resolution that gradually find the true edges of a quasicontinuum for m = −10, overlaid on a quasi-continuum calculated with ROC at two values of W com .Both the left and upper edges of the quasi-continuum are only fully resolved at a Legolas resolution of n grid = 1000.The latter has a value of ν ≈ 0.2, which is comparable to the results of Figure 3.The scattered modes found within the quasi-continuum bounds are part of dense clouds of eigenvalues when applying Arnoldi shift-invert to locations inside the 2D region (not shown).Furthermore, the eigenfunctions, shown in Figure 11(b), again resemble wave packages localized between the Alfvén resonances, which are now minimally separated (by only ∼1% of the radial domain).Hence, the criteria from Section 4 to classify this Legolas spectrum as a quasi-continuum have been met.Using a modified version of ROC that enables shooting between the very close-by Alfvén resonances rather than over the whole domain, we see that the eigenfunctions oscillate heavily for ν → 0 and that they fill the whole region between the resonances.All together, we have modes that are highly localized in all directions: azimuthally with m = −10, longitudinally with k = 400, and strongly oscillating radially in a small portion of the accretion disk, far away from the boundaries.Hence, these modes are truly localized, even though they are global in nature.It is hence not unthinkable that they have shearing box analogs, where only a tiny radial extent still produces overlapping continua that generate SARIlike modes (Matsumoto & Tajima 1995;Noguchi et al. 2000).

Discussion and Conclusions
The first part of this work considered the SARI quasicontinua discovered by Goedbloed & Keppens (2022).We independently confirmed their existence, noting that they are truly indistinguishable from true normal modes by any numerical code-and by physical plasma, for that matter.Although the mathematical treatment of these modes is quite novel, hints of their existence were already present in earlier works: the nonaxisymmetric spectra of Keppens et al. (2002) show large collections of modes that are reminiscent of the spectra in Figure 2, where the QR algorithm finds many scattered modes inside a quasi-continuum region.An earlier example is the modes excited in 2D simulations by Terquem & Papaloizou (1996), which have some resemblance with the quasi-mode eigenfunctions and move at the local Doppler speed.More analysis is needed to determine whether they are truly like the wave packages seen in Section 4. In any case, it is clear that a global treatment of the problem is essential to capture these effects: the nearness of singularities in the overlapping continua lies at the origin of the problem, and this cannot be obtained through a local analysis.With the definition of quasi-modes having a tiny but finite complementary energy on a firm physical and mathematical basis, we presented how quasi-modes can be identified using modern MHD spectroscopy with a tool like Legolas.The ROC and Legolas codes can be operated in tandem, so that the predictory value of the latter can be tested against the explicit W com calculations of the former, as was demonstrated in Figure 11.It will be interesting to establish whether a shearing box equivalent of SARIs exists in our formalism, resembling the radially localized nonaxisymmetric modes found by Matsumoto & Tajima (1995) and Noguchi et al. (2000).Global high-m modes with eigenfunctions that are adequately localized could perhaps even survive in the small radial variation of a local box.
In the second part, we performed a comprehensive study of instability growth rates in a wide range of equilibria, encompassing accretion disks with weak to strong, vertical to azimuthal, and subthermal to suprathermal magnetic fields.The stabilization of modes at various values of m was investigated for the most unstable inner and outer SARIs and the MRI.It was shown that stability is indeed governed by a condition (Equation ( 10)) on the Alfvén frequency that follows from analytical theory, which was used to explain growth rates for weakly magnetized, thin disks.A recent upper bound for the growth rates (Equation ( 9)) explained why the outer SARI growth rates behave like a scaled-down version of their inner equivalents.We also concluded that the m < 0 counterrotating SARIs play an important role in opening up more regimes for instability, significantly extending the range of unstable equilibria at high k, as shown in Figures 4 and 7(b).Our findings for nonaxisymmetric modes and their implications are summarized below.In particular, the maximal values of the numerically calculated growth rates are ∼0.6Ω 1 , in correspondence with the global analysis of Terquem & Papaloizou (1996) and Begelman et al. (2022), while Balbus & Hawley (1992) find growth rates of only several percent of Ω 1 in a local shearing sheet approximation.This again highlights the importance of using a global method when dealing with nonaxisymmetric modes.
Localization.In a thin, weakly magnetized disk, the fastestgrowing mode is found at intermediate to small wavelengths, with both k max and k cutoff increasing with increasingly negative m as instability depends on ω A , a conclusion also drawn by Begelman et al. (2022).For equilibria with a sizable toroidal field, the highest growth rates are found at large k (Balbus & Hawley 1992;Ogilvie & Pringle 1996).This implies that the essence of the instability is still captured in this approximation of an infinitely long cylinder, which requires the assumption of small vertical wavelengths compared to the scale height, and the modes found represent modes that are localized around the midplane of a 3D accretion disk (Balbus & Hawley 1998).In contrast, in models of accretion tori and 3D shearing boxes, unstable modes near the vertical edges of the disk are found (Balbus & Hawley 1998), and the strongest growth for the MRI occurs at wavelengths about the size of the vertical scale height H (Mamatsashvili et al. 2013).
Additionally, we showed in Section 5.5 that large-m modes have significant maximum growth rates, usually lower than those of the MRI but sometimes comparable.However, given the vast range of accessible wavelengths for nonaxisymmetric instabilities, even those with growth close to optimal can dynamically outperform the axisymmetric MRI (Mamatsashvili et al. 2013).In any case, the obtained growth rates are an order of magnitude higher than those applicable to the shearing sheet model by Balbus & Hawley (1992), as was also remarked by Curry & Pudritz (1996).The possibility of a locally vanishing ω A locally stabilizes SARI modes.We also explicitly showed that quasi-continua do exist in configurations with large m and k (Figure 11).Quasi-continuum SARI modes with such large wavenumbers are hence also very relevant to shearing box simulations, where the simulation domain is radially, vertically, and azimuthally localized (Lesaffre et al. 2009;Bai & Stone 2013;Hirai et al. 2018).The quasi-modes resulting from these spectra can corotate anywhere in the radial domain, making them ideal candidates for instability in both local Cartesian and global simulations.
In practice, however, large-k and large-m modes can only be realized as long as the scales involved are above the resistive and viscous length scales.The effect of these nonideal terms will generally be to quench the instability, particularly on smaller scales (Balbus & Hawley 1991, 1998;Kitchatinov & Rüdiger 2010).We plan to further investigate the effect of these diffusive terms on the MRI/SARI, especially on localized quasi-continuum SARI modes.Since the inclusion of nonideal MHD effects breaks the up-down symmetry of MHD spectra (Goedbloed et al. 2019), it will be interesting to see how the quasi-continuum regions are modified in the structure of the spectra.We note that Legolas can analyze a number of nonideal effects in combination, and thus can explore viscoresistive or Hall and ambipolar effects on the ideal MHD spectrum.
Magnetic field.We can conclude from Sections 5.2 and 5.3 that a strong poloidal field inhibits large-k instabilities, while a strong toroidal field inhibits large-m > 0 instabilities: the magnetic tension produced by these small-wavelength modes would be too high (Balbus & Hawley 1998).As a result, if the field is predominantly poloidal and strong, the MRI will be stabilized at small k while negative-m modes can still have significant growth rates for large k if ω A satisfies the criterion of Equation (10) (Begelman et al. 2022;Begelman & Armitage 2023), as was also confirmed by Figure 5.In particular, counterrotating SARIs have one value of m where the unstable range of field strengths is maximally extended, as in Figure 7(b).This implies that strong, poloidal fields could be unstable to extremely high-m modes beyond the capabilities of our linear (and nonlinear) codes: their eigenfunctions would be constrained to a tiny fraction of the domain.We also confirmed that for a purely toroidal field, a "complete breakdown of the most rapidly growing local modes" (Balbus & Hawley 1998) does indeed occur: the  ¥ k max for a purely toroidal field predicted by Ogilvie & Pringle (1996) quickly moves to much smaller values for even a seed poloidal field, as was noted in Figure 6.Finally, in our exploration of the plasma beta, we found that a dominant toroidal field has low growth rates for β < 1.This was also found by Das et al. (2018)  ).Not only the fastest-growing modes are of importance for stability: the linear evolution is governed by an entire MHD spectrum of eigenmodes.We found that the outer SARI branch is usually stabilized first, implying that the quasi-continuum between the discrete SARI branches can detach from the real axis and disappear before that, or continue to exist until the inner SARIs are also stabilized.It should be noted that we did not cover SARI quasi-continua for a dominant poloidal field, but these do exist.We confirmed this for m = 1 at weak field strength.The appendix covers the basic theory needed for the global formulation of the eigenvalue problem, from which a local dispersion relation can be derived.We show the connection between the general dispersion relation from Blokland et al. (2005) and recent work by Das et al. (2018) and Ebrahimi & Pharr (2022).By comparing solutions to this dispersion relation with a SARI spectrum in Figure 10, we highlighted that the local approximation is insufficient to describe outer and quasi-continuum SARI eigenmodes.
with nonorthogonal eigenfunctions (Goedbloed & Keppens 2022).In the Spectral Web method, the notion of a complementary energy is used to define exact modes and quasi-modes.This quantity is calculated from the total energy of an open system (Goedbloed et al. 2019): where W is the total energy, W p is the potential energy, and W com is the complementary energy.W p is real and calculated from the symmetric part of the integrand, whereas W com can be complex and represents a surface integral over the boundaries.For normal modes, the complementary energy is exactly zero, and so has vanishing real and imaginary parts.These two conditions can be used to construct two paths in the spectral plane on which eigenvalue solutions must lie.Physically, the complementary energy represents the additional energy required to bring a mode into resonance with the time dependence w -( ) i t exp . True normal modes are natural resonances and hence have vanishing W com .In the Spectral Web method, nonzero W com implies a discontinuity in the total pressure perturbation where the left and right solutions obtained from shooting from the boundaries should match.For quasi-continuum SARIs, this jump is so small that the modes can almost be regarded as normal modes.
The vector Equation (A4) in ξ can be reduced to a scalar ODE by replacing gradients in the azimuthal and vertical directions by their Fourier-transformed representations, reducing the equation to one for the Fourier coefficients x ˆ.The second and third components of x ˆcan then be written in terms of the first component ξ r .These expressions can then be used to write the equation of motion in terms of χ = rξ r alone as a second-order ODE: The singular solutions of Equation (A6) are produced by the coefficients N and D. Actually, the latter function only leads to apparent singularities, which becomes clear after rewriting Equation (A6) as a system of first-order ODEs (Blokland et al. 2005).The genuine singularities, resulting from the zeros of N, can be obtained explicitly: the functions c ˜and q can be determined from the real and imaginary parts of the equation and by using the WKB assumption (Blokland et al. 2005): Finally, a sixth-degree dispersion relation in w ˜is obtained from the first equality (Keppens et al. 2002), which is of the form where the coefficients can be found in Blokland et al. (2005).This equation describes the relation between w ˜, the wavenumbers (q, m, and k), and the equilibrium quantities at every location r ä [r 1 , r 2 ].This dispersion relation contains various local dispersion relations obtained for specific approximate cases, including the original standard weak-field MRI dispersion relation by Balbus & Hawley (1991) as shown in Blokland et al. (2005), and fourth-and sixth-order relations by Kim & Ostriker (2000) and by Balbus & Hawley (1998), respectively, as demonstrated in Keppens et al. (2002).The sixth-order dispersion relation Equation (B4) describes any general cylindrical disk equilibrium for any (m, k) combination.Its six solutions consist of four oscillating modes (the Alfvén and fast modes) and a conjugate pair of unstable modes (MRIor SARIlike modes).It was originally applied the overstability of the m = 0 mode, which is well described by this WKB approach.

B.1. Recent Dispersion Relations
With the above approach, a local dispersion relation is selfconsistently derived from the global framework of the generalized Hain-Lüst equation.However, most dispersion relations in the literature are derived by directly manipulating the linearized MHD equations and applying WKB arguments throughout.The two approaches are mutually consistent, since these relations are all contained in the general description by Blokland et al. (2005).In this section, we explicitly show this for two recent examples.Das et al. (2018) derived a general fourth-order dispersion relation from the linearized MHD equations, incorporating a suprathermal poloidal field, curvature, compressibility, and nonaxisymmetric modes (their Equation (40)).It is indeed contained in the description by Blokland et al. (2005) and can be obtained by applying the following approximations to Equation (B4): 1. assuming the fast frequencies are far away, i.e.,  w | ˜| qc hc ,

Figure
Figure 1(a).Figure 1(b) features a typical SARI spectrum for m = −10, considering a smaller domain with δ = 0.1 and the same equilibrium as that in panel (a).Note the large Doppler shift that moves the oscillation frequencies to negative values, making the continua W A (red and black) overlap in this super-Alfvénic regime.For these counterrotating SARIs (since m < 0), both Ω 0 and W  Figure 1(a).Figure 1(b) features a typical SARI spectrum for m = −10, considering a smaller domain with δ = 0.1 and the same equilibrium as that in panel (a).Note the large Doppler shift that moves the oscillation frequencies to negative values, making the continua W A (red and black) overlap in this super-Alfvénic regime.For these counterrotating SARIs (since m < 0), both Ω 0 and W 
Figure 2(b) shows the result of several such runs for m = 1.An annulus of eigenvalues is found around each shift in the quasi-continuum, with the inner and outer radii
(c) and (d) of Figure 2 show three quasi-mode eigenfunctions, sampled at various locations in the quasicontinuum as indicated in panel (b).All three eigenfunctions ( ) R v r resemble wave packages having a different number of oscillations and shifted to various locations in the disk.Note how the number of oscillations increases as w ( ) I decreases, but the eigenfunctions remain radially localized-in contrast to the MRI, which becomes more global as the number of oscillations increases.Like the discrete SARI eigenfunctions, each quasimode eigenfunction peaks at the local Doppler corotation radius (dashed, gray) and has one or two resonances with the Alfvén continua (red and black, dotted).The eigenmode in panel (c) has resonances with both continua.At these locations, large skin currents appear that effectively shield the mode from the radial boundaries (Goedbloed & Keppens 2022).This eigenmode then represents a truly radially localized wave package moving at the local Doppler speed, insensitive to the boundaries.How local a quasi-continuum mode is, is determined by the separation of these two resonances: if the overlap of W 
Figure 5(b) of Case (b) (slightly stronger field and β = 10) shows quite a different situation: at k = 50, both MRI and SARI are suppressed when the vertical field component dominates, and they are only activated for a sizable toroidal component.The "cutoff" inverse pitch angle is seen to decrease

Figure 4 .
Figure 4. Normalized growth rates for MRI (black, dashed-dotted) and SARI (inner: solid; outer: dotted) over a range of vertical wavenumbers k in Case (b) (dominant azimuthal field).Both k cutoff and k max increase as m decreases.The outer growth rates are scaled-down versions of the inner curves.Cusps in the growth rate curves, e.g., at k = 100 for the inner m = −10 SARI, signify local stabilization if ω A vanishes.
(a) and (b) again show Case (a) with j = 45°and Case (b) with a dominant toroidal field, respectively, both with k = 50 and β > 10.An immediate result is that for every k and m there is a bounded range of field strengths at which instability operates.For Case (a), all growth rates are maximal for v A1 ≈ 0.02v θ1 , well in the extremely super-Alfvénic regime, but for Case (b) the range of v A1 in which instability prevails for the different m is much larger, with v A1 /v θ1 close to unity for m < 0. The m = −3 SARI almost doubles the range of instability as compared to the MRI at this value of k.In both panels, the values of v A1 max and v A1cutoff shift with the value of m in a distinctive way.In panel (a), both v A1 max and v A1cutoff decrease as the wavenumbers go from m < 0 to m > 0. In panel (b), the shifts are more pronounced and do not obey the same ordering: remarkably, the m = −10 SARI is cut off at significantly weaker fields than other m < 0 SARIs, and its v A1 max

Figure 5 .
Figure 5. Normalized growth rates for MRI (black, dashed-dotted) and SARI (inner: solid; outer: dotted) for purely vertical (j = 0°) to purely azimuthal (j = 90°) background fields for Case (a) and Case (b) at k = 50.Note how a slight increase in field strength in panel (b) greatly reduces the allowed poloidal component.MRI stabilizes in a purely toroidal field; SARIs become more unstable.

Figure 7 .
Figure 7. Normalized growth rates of MRI (black, dashed-dotted) and SARI (inner: solid; outer: dotted) as a function of normalized inner Alfvén velocity (v A1 = B 1 ).A helical field (a) provides only a small range for instability at k = 50, in contrast to a dominant toroidal field (b), where especially the m < 0 SARIs are maximally unstable at much larger field strengths.

Figure 8 .
Figure 8. Normalized growth rates of MRI (black, dashed-dotted) and SARI (inner: solid; outer: dotted) varying with plasma beta at k = 50.In suprathermal fields, growth remains finite for a helical field (a) but converges to 0 for a dominant toroidal field (b).In panel (b), the MRI is not the most unstable mode at this value of k.

Figure 9 .
Figure 9. Alternative view of Figure 4, where the inner SARI growth rates now depend on ω A varying through k.Instability is perfectly governed by the criterion w < W  0 3 A 2 1 2 .Most growth rates peak at around w » W 15 16 A

Figure 10 .
Figure 10.The m = 1 quasi-continuum of Figure 2 filled using shift-invert.The prediction of n max of Equation (9) (green) is compared to the fastest-growing modes from a WKB dispersion relation (Equation (B4)) at three values for the radial wavenumber q.Clearly, a local model fails to predict the existence of quasi-modes, and that of the outer SARI.

0
-Lüst equation was first obtained in the context of accretion disks byKeppens et al. (2002) in a different form, and was rewritten byBlokland et al. (2005) in the form above.It contains function coefficients that are in general quite complicated expressions involving equilibrium quantities, the wavenumbers m and k, and the frequency ω.They can be found inGoedbloed et al. (2019) andGoedbloed & Keppens (2022), but we will repeat the most important definitions here.As noted before, the operator U produces a Doppler shift in the frequency.In cylindrical coordinates, the Doppler range is given by Note that the Doppler range has a radial dependence and vanishes for axisymmetric perturbations in our equilibrium without vertical flow.
The four singular ranges are dubbed forward and backward Alfvén and slow continua, as the frequencies lie in a continuum band of Alfvén frequencies w the Doppler range Ω 0 .The eigenfunctions corresponding to these continuum frequencies are nonsquare integrable functions with a singularity at the location of resonance with the continuum, i.e., where w w h 2 = m 2 /r 2 + k 2 .For brevity of notation, we use the same definitions of the equilibrium functions asGoedbloed et al. (2019): van der Swaluw et al. (2005)ths, but they identified a new hybrid Alfvén/slow instability at stronger fields.It is unclear if a similar instability appears for our m = 10 SARI in Case (a).According tovan der Swaluw et al. (2005), these instabilities must be more MRI-like than convection-like, since the parameter range explored here corresponds to thin disks (