The Super-Alfvénic Rotational Instability in Accretion Disks about Black Holes

Accretion processes in disks are of central importance in astrophysics, as disks arise in the vicinity of black holes at both galactic and stellar scales (in active galactic nuclei (AGNs) or X-ray binary systems), but they also play a crucial role in the context of protoplanet formation. The accretion disk radiative properties relate to the accretion rate (Shakura & Sunyaev 1973), which is in turn connected to angular momentum transport. The literature distinguishes many accretion disk types; see, e.g., Armitage (2011) and Turner et al. (2014) for reviews on protoplanetary disk physics, or Abramowicz & Fragile (2013) and Yuan & Narayan (2014) for reviews of accretion onto black holes. Disks may be varying in geometry between thick and thin disks, varying in radiative aspects being hot versus cold (optically thin vs. thick), or varying in whether advection, gasdynamical, or radiation effects are dominant. Many disk structures can be described in terms of essentially one-dimensional (height-integrated) radially self-similar accretion solutions, as pioneered by Shakura & Sunyaev (1973). These self-similar profiles involve an angular frequency Ω = v_θ (r)/r ∼ r^−p , which can deviate from perfect Keplerian flow (where p = 1.5) owing to, e.g., gas pressure or magnetic field effects, and ultimately relate the radially inward accretion velocity v_r ∼ −α r^−a to the celebrated α-parameter collecting "turbulent viscous" processes. In this v_r (r) profile, the accretion inflow power index a is taken according to prevailing physics, and, e.g., Spruit et al. (1987) argue for radial profiles where v_r ∼ − α r^−0.5 while the sound speed varies as r^−0.5, a property shared by the viscous advection-dominated accretion flows (ADAFs) discussed in Narayan & Yi (1994). These purely hydrodynamic prescriptions can be generalized in various ways, e.g., by incorporating convection effects, by yielding convection-dominated accretion flows (CDAFs; Narayan et al. 2000), or by allowing for a magnetic field at a fixed plasma beta (Yuan et al. 2012). The latter really requires the use of magnetohydrodynamics (MHD), which has become a central working tool for accretion physics as a whole. Indeed, in various reviews on accretion disk theory (see, e.g., Balbus & Hawley 1998), the role of MHD turbulence in the disk is acknowledged as essential for accretion to proceed, as the fully developed turbulent state would realize the corresponding α-values to agree with observed disk properties. The ultimate trigger for disk turbulence is in the ideal MHD regime invariably ascribed to a linear MHD instability, the magnetorotational instability (MRI). Naturally, extensions beyond ideal MHD have rightfully gained attention, especially in the protoplanetary disk context (Turner et al. 2014), where resistive, ambipolar, and Hall MHD effects came into focus, introducing MRI-active versus "dead" zones (see, e.g., Lesur et al. 2014, Lesur 2021). Cartoon views on protoplanetary disks identify where ideal MHD, MRI-like turbulence is still applicable, or where only weak turbulence is expected, and how the α-parameter affects the growth of coupled gas-particle streaming instabilities (Umurhan et al. 2020). For both protoplanetary and black hole disks, the way the observationally inferred α-parameter (or the related central accretion rate; see Rafikov 2017) maps to turbulence-induced angular momentum transport is in reality also complicated by (magnetized) winds and/or jets in the accreting object + disk environment, which are known to influence the overall angular momentum budget (see, e.g., Casse & Keppens 2002, 2004).

The MRI was introduced as a "powerful local shear instability in weakly magnetized disks" in the seminal work by Balbus & Hawley (1991). Their work presented dispersion relations from linear ideal MHD for axisymmetric, linear perturbations about a weakly magnetized disk of finite vertical extent, and all it requires for instability is a radially decreasing angular velocity law where ${\rm{\Omega }}^{\prime} \lt 0$ and a weak poloidal B_z magnetic field component. Their linear MHD treatment of a disk, with added (but weak) toroidal B_θ (r, z) and poloidal B_z (r) components, assumed a Boussinesq approximation and analyzed $\exp [{\rm{i}}({k}_{r}r+{k}_{z}z-\omega t)]$ perturbations with radial and vertical mode numbers k_r and k_z . In an accompanying paper (Hawley & Balbus 1991), 2.5D axisymmetric nonlinear MHD simulations were performed of a magnetized Keplerian, cylindrical Couette flow, where the presence of a weak magnetic field was confirmed as essential to get growing linear mode behavior, and where its role in triggering turbulence and angular momentum redistribution was demonstrated for the first time. Further 2.5D nonlinear MHD simulations, employing the "shearing sheet" approximation as appropriate for disks, provided additional support (Hawley & Balbus 1992). The MRI, whose essential ingredients were analyzed earlier by Velikhov (1959) and Chandrasekhar (1960), has ever since been paramount to understanding and explaining accretion disk behavior. State-of-the-art research invokes MRI in countless studies. For example, in protoplanetary disks, Xu & Bai (2022) investigate how MRI turbulence in the outer disk zones interplays with gas-coupled particle prescriptions (dust), where gas−dust velocity differences affect dust clumping by streaming instabilities. The "gas" is then described by standard MHD equations, extended with ambipolar diffusion to account for the influence of partial ionization. In modern studies of accretion flows about black holes, general relativistic MHD (GRMHD) simulations serve to interpret observations by Event Horizon Telescope Collaboration et al. (2019), or observed multiwavelength variability in light curves (Chatterjee et al. 2021), whereas they employ "MRI-quality factors" that dictate the needed numerical resolutions to "resolve MRI turbulence." General relativistic, resistive MHD simulations nowadays address how plasmoids may form within current sheets in black hole accretion disks (Ripperda et al. 2020; Begelman et al. 2022), clearly differing for "Standard And Normal Evolution" (SANE) versus "Magnetically Arrested Disks" (MAD) states, whereby SANE disks are particularly prone to MRI turbulence throughout. Also in Newtonian settings, modern 3D resistive MHD runs address plasmoid reconnection aspects generated by "primary MRI" (Rosenberg & Ebrahimi 2021). Numerous authors performed MRI studies in shearing boxes, with, e.g., Simon & Hawley (2009) addressing visco-resistive modifications, or recently Held & Latter (2022) studying the effect of parameterized cooling prescriptions, where MRI interacts with thermal instability cycles.

The original axisymmetric MRI treatment (Balbus & Hawley 1991) presented a purely local stability analysis, akin to a WKB treatment, and various later studies highlighted important novel aspects. For example, Knobloch (1992) revisited the axisymmetric MRI using linear incompressible MHD in a more rigorous way, pointing out the role of boundary conditions (BCs) when performing a local stability analysis of shear flows, as well as the fact that any azimuthal field will turn the purely growing mode into an overstable oscillation. This was further elaborated on by Blokland et al. (2005), where axisymmetric modes in a stratified disk, as analyzed using the WKB approach, were confirmed by a complete numerical treatment of the linear compressible MHD equations. In terms of the traditional plasma beta parameter β ≡ 2p/B², this study identified how m = 0 MRIs could persist up to equipartition β ≈ 1 conditions, when dominant toroidal magnetic field components are present.

Although the initial 2.5D nonlinear MHD simulations (Hawley & Balbus 1991, 1992) gave clear hints for efficient angular momentum transport due to MRI, the assumption of axisymmetry precluded definitive statements on its importance for sustaining 3D MHD turbulence, as needed for, e.g., disk dynamo processes. Purely hydrodynamic 3D simulations of accretion tori did show nonaxisymmetric (low azimuthal mode number m) spirals forming (Hawley 1991), but no hints of actual turbulence. That linearly unstable m = 1 modes indeed exist for hydrodynamic tori with constant specific angular momentum was predicted analytically by Papaloizou & Pringle (1984). To make the case for the MRI as relevant for triggering 3D MHD turbulence, nonaxisymmetric modes for a thin disk were addressed in Balbus & Hawley (1992), where all radial structure was ignored except for the shear flow resulting from Ω(r). The analysis used time-dependent radial mode numbers where ${{dk}}_{r}/{dt}=-m{\rm{\Omega }}^{\prime}$ and demonstrated transient amplification of magnetic field perturbations, up to many orders of magnitude. Further strengthened by full 3D nonlinear MHD simulations, originally in the local "shearing box" approach (Hawley et al. 1995), but later (Hawley 2000) in global 3D MHD simulations of accretion tori, all these findings essentially appeared to close the case on MRIs, as relevant (1) to explain turbulence in disks, (2) for providing efficient outer angular momentum fluxes due to magnetic stresses to get accretion going, and (3) for their role in disk magnetic field amplification.

However, as this paper intends to show, many crucial aspects of especially nonaxisymmetric linear MHD instabilities in realistic accretion disks have been completely ignored thus far. This is mainly because—with few exceptions (e.g., Ogilvie 1998)—the astrophysical literature on MHD disk eigenmodes does not exploit the full power of MHD spectroscopy, where a key role is played by the continuous parts in the MHD spectrum. This "art" of quantifying the complete spectral eigenstructure of a given MHD equilibrium state has matured in the laboratory fusion context and has been transferred to astrophysical settings as documented in our textbook (Goedbloed et al. 2019). The variation of all MHD equilibrium quantities introduces both continuous and discrete sequences of eigenmodes, which have been shown to organize and link all known stable and unstable (overstable/damped) modes in intriguing ways. By using a fully general treatment of MHD eigenmodes in the cylindrical disk approach (see Section 1.2), we therefore aim to show (1) that the axisymmetric MRIs represent only a finite amount of actually unstable (i.e., overstable) modes, within an otherwise infinite clustering sequence of mostly stable modes accumulating toward the continuous spectra; (2) that a novel type of nonaxisymmetric modes exists, which we denote as super-Alfvénic rotational instabilities (SARIs), as a result of overlapping continuum frequency ranges, and these are truly infinite sequences consisting of unstable modes only; (3) the existence of "virtual walls" (due to skin currents at near singularities) in the localization of the SARI eigenfunctions, which makes them insensitive to adopted BCs, unlike the m = 0 MRIs; and (4) the completely uncharted role of extremely localized quasi-continuum, nonaxisymmetric instabilities, which seemingly occupy complete regions in the complex eigenfrequency plane. Together, these findings open up an entirely new window on linear MHD modes, far beyond the celebrated MRI mode, as likely triggers for MHD turbulence, accretion, and dynamo activity in disks.

Of course, even in purely hydrodynamical settings, the important role of nonaxisymmetric perturbations on (thin, viscous) disks is well established, with possibilities for transient growth as in Rebusco et al. (2009), or m ≠ 0 Rossby wave instabilities induced by local overdensities, as in Meheut et al. (2012), as well as various other pure hydro effects discussed in the protoplanetary context by Turner et al. (2014). We will focus here on magnetized disks, extending the ideal MHD theory way beyond the original m = 0 MRI, but point out that our findings have potential implications for purely hydro settings as well, as discussed in Section 6.2. Various previous studies investigated nonaxisymmetric modes in magnetized sheared flows. Ogilvie & Pringle (1996) presented a treatment in a cylindrical setting where gravity was ignored, the density was homogeneous, and a purely toroidal, potential B_θ ∼ r⁻¹ field was incorporated. That study did admit the role of forward and backward Alfvén continua and found discrete nonaxisymmetric modes whose limit points could be obtained for large axial wavenumber k_z → ∞. Nonaxisymmetric MRI modes, called azimuthal MRIs or AMRIs, were also found in a Taylor–Couette setting in visco-resistive MHD with a purely toroidal imposed field (Rüdiger et al. 2007a, 2007b; Hollerbach et al. 2010). The latter works extended earlier findings (Hollerbach & Rüdiger 2005) on m = 0 modes in helical magnetic fields in the same Taylor–Couette context, and this "helical MRI" is still a topic of concern for modern liquid sodium experiments (Mishra et al. 2021). How a uniform, axial magnetic field can open the way to nonaxisymmetric modes in dissipative disks is studied in Kitchatinov & Rüdiger (2010). Truly 2D disks with purely toroidal B_θ (r, z) magnetic fields were analyzed as well in Terquem & Papaloizou (1996), pointing out that such disks are always unstable. Curry & Pudritz (1996) studied nonaxisymmetric modes in incompressible cylinder settings and thereby generalized the hydrodynamic Papaloizou–Pringle (Papaloizou & Pringle 1984) mode to magnetic regimes, while finding additional pure Alfvén modes coupled to the local rotation frequency of a cylindrical shell. Our treatment will generalize and augment these findings considerably, since we will allow for arbitrary poloidal B_z (r) and toroidal B_θ (r), incorporating all subtleties due to the MHD continuous spectra, while we intentionally address the ideal MHD regime only, which we believe to be central to any further modification by nonideal effects. This will involve completely new spectral structures as organized in the Spectral Web, which allows us to reveal the surprising spread of modes throughout the complex eigenfrequency plane in unprecedented ways.

It must be emphasized that the countless nonlinear MHD simulations of turbulent disks, nowadays extended from Newtonian to fully general relativistic settings, covering aspects from magnetized protoplanetary disks, to disks in X-ray binaries, all the way up to galactic disks and AGNs, do indeed show clearly that MHD turbulence is virtually inevitable. In that respect, our general linear MHD treatment serves to seriously question the present attitude throughout the astrophysical literature, ascribing all turbulence in disks to MRI activity. In 3D shearing box simulations, it is clear that precisely the high-m modes are the ones of interest. In global disk settings, the 3D time-dependent (ideal) MHD state obtained at any individual moment may well be analyzed spectrally, and it has been shown (Keppens & Demaerel 2016) that the same mathematical operators (a generalized force operator and the Doppler–Coriolis operator) built up from these instantaneous (ideal) MHD profiles govern the state's stability at any time in its nonlinear evolution. It would be a major oversimplification to not discriminate mode properties depending on the state at hand, which may transit in its evolution from weak to strong field and has disks in thin to thick settings, with locally different partitions between rotational, magnetic, thermal, and radiative energies. In that respect, there is a clear dichotomy in the astrophysical literature, where the nonlinear, height-averaged equilibrium accretion flow of magnetized disks is being categorized in various states known as SANE, MAD, ADAFs, CDAFs, thin, thick, slim, radiatively efficient or inefficient, etc., while all turbulent activity is simply denoted as a direct consequence of the linear MRI.

In the established approach to diagnose all waves and instabilities about a stationary (i.e., time-independent) ideal MHD state (see Section 2.2), one must always specify the fully force-balanced equilibrium state up front. In the context of accretion disks, this equilibrium can be taken to be a magnetized, axisymmetric accretion torus, where the combination of pressure gradients, external gravity due to the central object, centrifugal and inertial effects, and Lorentz forces balance. This in essence renders the background equilibrium 2D in nature, with profiles ρ(r, z), p(r, z) and magnetic and flow vector quantities that in turn determine the linear eigenmode distribution. In the laboratory fusion context, magnetic flux surfaces forming nested tori are realized within tokamak configurations, and the same magnetic topology can be adopted for initializing magnetized accretion tori equilibrium states. Rigorous analysis of the MHD spectrum of such 2D equilibria has been undertaken and has already identified two new (i.e., distinct from MRI) likely routes to turbulence, related to the fact that purely flux surface localized modes, which define the continuous parts of the MHD spectrum, can actually be driven unstable. This invariably involves intricate coupling schemes between multiple linear modes of different poloidal mode number: the 2D poloidal (r,z) variation of the equilibrium itself now causes Fourier modes of different poloidal angular variation to couple. In Goedbloed et al. (2004a), as well as in Chapter 18 of Goedbloed et al. (2019), a new class of local instabilities called the trans-slow Alfvén continuum (TSAC) modes was identified, whenever the poloidal Alfvén Mach number M (found from ${M}^{2}\equiv \rho {v}_{p}^{2}/{B}_{p}^{2}$ using the poloidal vector components) of the equilibrium would exceed a critical value. In a later study, Blokland et al. (2007) identified unstable continuous spectra related to convective modes, aptly called convective continuum instabilities (CCIs), that even persist in strongly magnetized disks (a property known to suppress the standard MRI). The same is true for the TSAC modes, since in the strong-field limit β ≪ 1 of magnetically dominated thick accretion tori, this continuum instability switches on whenever the squared poloidal Alfvén Mach number exceeds a value of about $\tfrac{1}{2}\gamma \beta$, and this for all poloidal–toroidal mode number pairs. At the same time, Haverkort & de Blank (2012) clearly showed that in such 2D equilibrium settings as valid for accretion tori, a radially decreasing toroidal rotation frequency can be stabilizing for nonaxisymmetric MHD modes, by a Coriolis pressure effect. This finding is relevant both for accretion disks and for rotating, toroidally confined laboratory plasmas (which are known to feature transport barriers).

Vertical stratification of the disk also features in the analyses of Marcus et al. (2013) and Goodman & Xu (1994). In both papers the main thrust is on the nonlinear phase of the perturbations, grown either from a linearly stable state (like the hydrodynamic Keplerian disk analyzed in the first paper) or from a linear instability (like the MRI in magnetized accretion disks analyzed in the second paper). The paper by Marcus et al. (2013) considers a purely hydrodynamic vertical flow field, sheared in the horizontal direction. While linearly stable, it generates sequences of vortices that advect in the cross-stream directions. It is conjectured that this mechanism of finite-amplitude instability may destabilize hydrodynamic Keplerian rotations in protoplanetary disks with a negligible magnetic field. In the paper by Goodman & Xu (1994), on the other hand, a magnetized plasma disk with small magnetic fields subject to MRIs is considered. The problem posed is whether these instabilities survive at finite amplitudes. It is indicated that the nonlinear development of secondary parasitic instabilities will stop the growth of the primary MRIs. However, an "exact" finite-amplitude MRI is shown to occur for magnetic field perturbations of even equipartition strength! As admitted by the authors, this analysis is quite restricted in the number of degrees of freedom that can be handled. It gives additional significance to the MRI picture of the initial phase of turbulence in disks, though, which might equally well apply to the nonaxisymmetric instabilities discussed in the present paper.

For the time being, we will stick to the consideration of magnetic fields in the overall picture of accretion and the ensuing formation of jets. Also, we restrict the analysis to linear perturbations, awaiting a nonlinear approach of the nonaxisymmetric instabilities, which may show either saturation at low amplitude levels of the instabilities or further growth to large amplitudes, with the possible generation of parasitic instabilities, sequences of vortices, etc. As the history of the MRI illustrates, such analyses can only come after the linear version of the dynamics has been sufficiently documented. That is the implicit goal of the present paper.

Hence, we will rather adopt a purely 1D "cylindrical disk" model and analyze the full consequences of having a radially stratified Ω(r), density ρ(r), pressure p(r), and further arbitrary B_z (r) and B_θ (r). This cylindrical disk approach has the distinct advantage that the poloidal couplings mentioned as crucial to the TSAC and CCI unstable continuum mode types do not occur at all: a linear analysis can focus attention on a single mode $\exp [{\rm{i}}(m\theta +{kz}-\omega t)]$ at a time, with prechosen mode numbers m and k. The radial variation of the eigenmodes must always be computed as part of the eigenmode quantification, unless some WKB-type assumption can be justified. The basic governing equations for such cylindrical disk eigenmodes in ideal MHD have been known for some time and were presented in Keppens et al. (2002), building on all previous efforts that excluded the influence of a (cylindrical) gravitational potential. Due to the 1D nature of the equilibrium, the continuous—and always stable—parts of the MHD spectrum are known for each given set of mode numbers (m,k). In the mentioned paper, the axisymmetric MRI modes of a weakly magnetized disk were computed, showing their relation to the continua, and clear indications of a richer nonaxisymmetric eigenmode structure were given. The present paper revives this uncharted field of magnetoseismology of accretion disks, by using the powerful Spectral Web technique to locate all the eigenmodes. In this approach, we can choose the equilibrium properties according to the disk model believed to prevail, e.g., take profiles in accord with the ADAF, CDAF, or whichever accretion flow solution one realizes (e.g., using instantaneously height-averaged profiles from a fully 3D MHD simulation), and then show how different unstable modes appear owing to the intricate interplay of the various restoring forces that are present.

Consider the cylindrical slice model of an accretion disk rotating about a compact or protostellar object of mass M_* (Figure 1, adapted from Goedbloed et al. 2019). The vertical variation of the plasma equilibrium is neglected so that we obtain an annular cylindrically symmetric slice of height Δz and located within the radial range r₁ ≤ r ≤ r₂. The disk is confined by the gravitational field of the central object,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{gr}}=-\displaystyle \frac{{{GM}}_{* }}{\sqrt{{r}^{2}+{z}^{2}}}\approx -\displaystyle \frac{{{GM}}_{* }}{r}\quad \Rightarrow \quad g={{\rm{\Phi }}}_{\mathrm{gr}}^{\prime} ={{GM}}_{* }{/r}^{2}\,\left(^{\prime} \equiv d/{dr}\,\right),\end{eqnarray} \tag{ 1 }$$

where the approximation on the right-hand side (RHS) is justified by the assumption of short wavelengths in the vertical direction, kΔz ≫ 1. The Newtonian potential exploited could be replaced by the Paczyński–Wiita potential (Paczyński & Wiita 1980), to approximate general relativistic effects, but this complication is avoided here assuming that r₁ is far outside the Schwarzschild radius. Also, accretion is supposed to take place on a timescale much longer than that of the instabilities and ensuing turbulence that cause it, so that the radial velocity v_r is negligible with respect to the rotation velocity v_θ . The equilibrium is then described by rather arbitrary radial distributions of the following quantities. The disk has a density ρ(r), it exerts a pressure p(r), it is rotating with angular frequency Ω ≡ v_θ (r)/r, and it is immersed in a helical magnetic field with a toroidal component B_θ (r) and a vertical component B_z (r). The only restriction on these five profiles is that they should satisfy the radial equilibrium constraint

$$\begin{eqnarray}&&\rho r({{\rm{\Omega }}}^{2}-{{GM}}_{* }/{r}^{3})=(p+\tfrac{1}{2}{B}^{2})^{\prime} +{B}_{\theta }^{2}/r,\end{eqnarray} \tag{ 2 }$$

where the RHS pressure and magnetic field contributions produce deviations from the Keplerian rotation on the left-hand side. A vertical flow field v_z (r) would not modify this equilibrium, and it has also been incorporated in the general analysis of Section 13.3 of Goedbloed et al. (2019), on which this paper is based, but it is put to zero in the present paper. The kinetic pressure is assumed to be large with respect to the magnetic pressure, β ≡ 2p/B² ≫ 1, so that the perturbations of this equilibrium are assumed to be incompressible in the explicit analyses of Sections 3–5. The numerical results from the program ROC (Goedbloed 2018a, 2018b) presented throughout the paper are not restricted by this condition; the program exploits the fully general spectral differential equation solver.

Zoom In Zoom Out Reset image size

We will utilize the scale independence of the ideal MHD equations (see Section 4.1.2 of Goedbloed et al. 2019) by exploiting units of length, mass, and time, effectively provided by the distance from the origin to the reference position r = r₁, the density ρ₁, and the Keplerian velocity ${v}_{{\rm{K}}1}\equiv \sqrt{{{GM}}_{* }/{r}_{1}}$ at that position. This yields dimensionless quantities $\bar{r}\equiv r/{r}_{1}$, $\bar{\rho }\equiv \rho /{\rho }_{1}$, $\bar{{\rm{\Omega }}}\equiv {r}_{1}{\rm{\Omega }}/{v}_{{\rm{K}}1}$, ${\bar{B}}_{z,\theta }\equiv {B}_{z,\theta }/(\sqrt{{\mu }_{0}{\rho }_{1}}{v}_{{\rm{K}}1})$, etc. (In the plots, the coordinate $\bar{r}$ will be indicated by x, as shown in Figure 1.) From now on, we get rid of these trivial scaling parameters by dropping the bars and exploiting the resulting dimensionless quantities. This rescaling also gets rid of the awkward constant μ₀ of the magnetic variables, as anticipated in Equation (2). In contrast, four essential parameters will describe all our numerical examples, viz.,

$$\begin{eqnarray}&&\epsilon \equiv {B}_{1}\equiv {\left({B}_{z1}^{2}+{B}_{\theta 1}^{2}\right)}^{1/2},\quad \beta \equiv 2{p}_{1}/{B}_{1}^{2},\quad {\mu }_{1}\equiv {B}_{\theta 1}/{B}_{z1},\quad \delta \equiv ({r}_{2}-{r}_{1})/{r}_{1},\end{eqnarray} \tag{ 3 }$$

where we have replaced the previous definitions of and δ used in Goedbloed et al. (2019) and Goedbloed (2018b) by the present, more relevant ones. In particular, the parameter now directly measures the magnitude of the magnetic field, which will be argued below to be the main dynamical quantity (next to rotation of course). The subscript on the constant μ₁ is essential since the inverse pitch of the field lines is not constant: μ ≡ B_θ /(rB_z ) = μ₁ r⁻¹. Adopting the standard self-similar radial dependence of accretion disks introduced by Spruit et al. (1987), we then obtain the following explicit equilibrium:

$$\begin{eqnarray}&&\rho ={r}^{-3/2},\quad {\rm{\Omega }}={{\rm{\Omega }}}_{1}{r}^{-3/2},\quad {B}_{z}={B}_{z1}{r}^{-5/4},\quad {B}_{\theta }={\mu }_{1}{B}_{z1}{r}^{-5/4},\quad p=\tfrac{1}{2}\beta {\epsilon }^{2}{r}^{-5/2},\end{eqnarray} \tag{ 4 }$$

where the angular rotation frequency Ω₁ follows from the equilibrium Equation (2):

$$\begin{eqnarray}&&{B}_{z1}=\displaystyle \frac{\epsilon }{{\left(1+{\mu }_{1}^{2}\right)}^{1/2}}\quad \Rightarrow \quad {{\rm{\Omega }}}_{1}={\left\{1-\tfrac{5}{4}{\epsilon }^{2}\left[\beta +\displaystyle \frac{1+\tfrac{1}{5}{\mu }_{1}^{2}}{1+{\mu }_{1}^{2}}\right]\right\}}^{1/2}\approx 1-\tfrac{5}{8}{\epsilon }^{2}\beta .\end{eqnarray} \tag{ 5 }$$

Approximating Ω₁ ≈ 1 (i.e., Keplerian rotation at r₁ and, hence, everywhere) would result from the assumptions ≪ 1 and ² β ≪ 1 (but β ∼ ⁻¹ ≫ 1) that are made for the examples in this paper. These assumptions also imply that the equilibrium rotation velocities considered will be extremely super-Alfvénic: v_θ ∼ 1 ≫ v_A ∼ . It is to be noted that the parameter β, even though large, disappears from the spectral analysis since it leads to incompressibility of the modes. Hence, somewhat paradoxically, the characteristic ordering of importance for these accretion disk problems is (1) dominant rotation (quantified by the Keplerian velocity v_K1), (2) small magnetic field magnitude (quantified by the parameter ), and, finally, (3) larger pressure effects (quantified by the parameter β ²) that are of no importance, however, for the dynamics of the modes described in this paper. The radial distribution of the equilibrium variables for a typical cylindrical accretion disk is shown in Figure 2.

Zoom In Zoom Out Reset image size

Many more self-similar equilibria could be considered. In particular, the functions ρ(r), ${B}_{\theta }^{2}(r)$, ${B}_{z}^{2}(r)$, and p(r) could be multiplied by an arbitrary power of r. For definiteness and comparison with the existing literature (Shakura & Sunyaev 1973; Spruit et al. 1987), we have set that power equal to 1 in expressions (4). Again, the numerical tools ROC (Goedbloed 2018a, 2018b), exploited in this paper, or Legolas (Claes et al. 2020), exploited in a forthcoming paper, are not restricted to these equilibria.

The most effective starting point for the study of waves and instabilities of stationary equilibria, such as the present one, is the seminal paper of Frieman & Rotenberg (1960). For modes exponentiating as $\exp (-{\rm{i}}\omega t)$, their spectral differential equation may be written as

$$\begin{eqnarray}&&{\boldsymbol{G}}({\boldsymbol{\xi }})-2\omega U{\boldsymbol{\xi }}+\rho {\omega }^{2}{\boldsymbol{\xi }}=0,\end{eqnarray} \tag{ 6 }$$

where ${\boldsymbol{G}}({\boldsymbol{\xi }})\equiv {\boldsymbol{F}}({\boldsymbol{\xi }})+{\rm{\nabla }}\cdot \left({\boldsymbol{\xi }}\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{v}}\right)-\rho {\left({\boldsymbol{v}}\cdot {\rm{\nabla }}\right)}^{2}{\boldsymbol{\xi }}$ is the force operator expression, generalizing the more well-known expression F ( ξ ) for static equilibria (Bernstein et al. 1958), and U ≡ − iρ v · ∇ is the gradient operator projected onto the velocity field. With suitable BCs on ξ, Equation (6) becomes a quadratic eigenvalue problem in terms of the two operators G and U, which are both self-adjoint. Hence, the two associated quadratic forms

$$\begin{eqnarray}&&W\equiv -\tfrac{1}{2}\int {{\boldsymbol{\xi }}}^{* }\cdot {\boldsymbol{G}}({\boldsymbol{\xi }})\,{dV}\quad \mathrm{and}\quad \ V\equiv \tfrac{1}{2}\int {{\boldsymbol{\xi }}}^{* }\cdot U{\boldsymbol{\xi }}\,{dV}\end{eqnarray} \tag{ 7 }$$

are real for solutions satisfying the BCs. The meaning of the ambiguous symbol V, referring to both the integral and the plasma volume, will be evident from the context. Introducing the normalization $I\equiv \tfrac{1}{2}\int \rho | {\boldsymbol{\xi }}{| }^{2}\,{dV}$, Equations (6) and (7) yield a quadratic equation for the complex variable ω ≡ σ + iν:

$$\begin{eqnarray}&&{\omega }^{2}-2\overline{V}\omega -\overline{W}=0,\qquad \overline{V}\equiv V/I,\qquad \overline{W}\equiv W/I,\end{eqnarray} \tag{ 8 }$$

where $\overline{V}$ and $\overline{W}$ are the solution averages of the Doppler–Coriolis shift and of the potential energy of the perturbations. The solutions of this quadratic equation for instabilities, with real frequency σ and growth rate ν, may be written as

$$\begin{eqnarray}&&\sigma =\overline{V},\qquad \nu =\pm {\left(-\overline{W}-{\overline{V}}^{2}\right)}^{1/2}.\end{eqnarray} \tag{ 9 }$$

It would look like: problem solved! However, these expressions cannot be evaluated a priori. In the standard approach, they can only be computed a posteriori when the final solution of the spectral Equation (6), consisting of the eigenvalue–eigenfunction pair {ω, ξ ( r ; ω)}, is known. This obstacle is overcome by the new method of the Spectral Web described in Section 2.4.

The complex eigenvalues of the waves and instabilities of the stationary cylindrical equilibria described in Section 2.1 are determined by the solutions of an ordinary differential equation (ODE) for the radial component χ ≡ r ξ_r of the plasma displacement, which may be obtained from the Frieman–Rotenberg spectral Equation (6) by straightforward reduction. We skip over all special cases of this ODE that have been exploited over the past 60 yr and, for reference, just indicate the general forms for different assumptions of the underlying equilibrium. For static isothermal (γ = 1) cylindrical equilibria, the relevant second-order ODE was derived by Hain & Lüst (1958) and generalized to adiabatic equilibria by Goedbloed (1971). The stationary counterpart was derived by Hameiri (1981) and numerically investigated by Bondeson et al. (1987). These studies were concerned with stability of magnetically confined laboratory plasmas for thermonuclear purposes. For astrophysical applications, the extension to the completely general cylindrical ODE with a gravitational field, applied to the various instabilities of an accretion disk rotating about a massive central object, e.g., the MRI (Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991, 1998), was obtained by Keppens et al. (2002).

We exploit the latter form of the cylindrical spectral equation, as given in Section 3.3.2 of Goedbloed et al. (2019) for normal modes of the form $\chi (r)\exp [{\rm{i}}(m\theta +{kz}-\omega t)]$:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left(\,\displaystyle \frac{N}{D}\,\displaystyle \frac{d\chi }{{dr}}\,\right)+\left[\,A+\displaystyle \frac{B}{D}+\left(\displaystyle \frac{C}{D}\right)^{\prime} \,\right]\,\chi =0,\end{eqnarray} \tag{ 10 }$$

where the modes are supposed to be localized within "rigid" boundaries at r = r₁ and r = r₂ represented by the BCs

$$\begin{eqnarray}&&\chi ({r}_{1})=0\ \ (\mathrm{left}),\qquad \chi ({r}_{2})=0\ \ (\mathrm{right}).\end{eqnarray} \tag{ 11 }$$

Of course, the actual boundaries of an accretion disk will not comply with the assumption of an annular cylindrical slice, so that we have to convince ourselves that, in the end, the eigenfunctions computed will be well within the range r₁ ≤ r ≤ r₂. For numerical integration it is expedient (Appert et al. 1974) to transform ODE (10) into a pair of complex first-order ODEs in terms of χ and the total pressure perturbation, ${\rm{\Pi }}\equiv {p}_{1}+{B}_{0}{B}_{1}\equiv -(N/D)\chi ^{\prime} -(C/D)\chi$:

$$\begin{eqnarray}&&N\displaystyle \frac{d}{{dr}}\left(\begin{array}{c}\chi \\ {\rm{\Pi }}\end{array}\right)+\left(\begin{array}{cc}C & D\\ E & -C\end{array}\right)\left(\begin{array}{c}\chi \\ {\rm{\Pi }}\end{array}\right)=0,\end{eqnarray} \tag{ 12 }$$

where the coefficient E is defined in terms of the other coefficients,

$$\begin{eqnarray}&&E\equiv -N\left(A+B/D\right)-{C}^{2}/D.\end{eqnarray} \tag{ 13 }$$

Of course, all physics now resides in the expressions for the coefficients N, D, A, B, and C, which are functions of the radial position r and of the Doppler-shifted complex frequency $\widetilde{\omega }$, which is also a function of r:

$$\begin{eqnarray}&&\widetilde{\omega }\equiv \omega -{{\rm{\Omega }}}_{0},\qquad {{\rm{\Omega }}}_{0}\equiv m{\rm{\Omega }}+{{kv}}_{z}.\end{eqnarray} \tag{ 14 }$$

For generality, we have kept the vertical velocity v_z , though it will be neglected in the rest of this paper. In contrast to the usual analysis of axisymmetric MRIs (m = 0), where $\widetilde{\omega }=\omega$ is constant, the radial profile Ω₀(r) plays a central role in the description of the nonaxisymmetric modes (m ≠ 0), so that $\widetilde{\omega }$ is not constant. This significantly complicates the analysis but, of course, also greatly enlarges the possible dynamics of accretion disks about compact objects. In particular, ODEs (10) and (12) now become a fourth-order system for the determination of the four components χ₁, χ₂, Π₁, and Π₂ of the essentially complex variables χ ≡ χ₁ + iχ₂ and Π ≡ Π₁ + iΠ₂.

Recall that the collection of real frequencies ω = {Ω₀(r)∣r₁ ≤ r ≤ r₂} forms the flow continuous spectrum in hydrodynamics (Case 1960). In magnetized plasmas, these frequencies still play an important role, as we will see, but their role as a continuum is superseded by the forward and backward Doppler-shifted Alfvén and slow magnetosonic continuum frequencies

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{A}}}^{\pm }\equiv {{\rm{\Omega }}}_{0}(r)\pm {\omega }_{{\rm{A}}}(r),\qquad {{\rm{\Omega }}}_{{\rm{S}}}^{\pm }\equiv {{\rm{\Omega }}}_{0}(r)\pm {\omega }_{{\rm{S}}}(r),\end{eqnarray} \tag{ 15 }$$

where ${\omega }_{{\rm{A}}}\equiv F/\sqrt{\rho }$, F ≡ mB_θ /r + kB_z , and ${\omega }_{{\rm{S}}}\equiv {[\gamma p/(\gamma p+{B}^{2})]}^{1/2}{\omega }_{{\rm{A}}}$ are the static expressions. The latter square root factor ≈ 1 − (4γ β)⁻¹ ≈ 1 for β ≫ 1, so that the modes are essentially incompressible and the difference between the slow and Alfvén frequencies disappears. Note that these frequencies are real! The definitions of the singularity coefficients N and D can then be written as

$$\begin{eqnarray}\begin{array}{rcl}N & \equiv & \widetilde{A}\,\widetilde{S}/r,\qquad D\,\equiv {\rho }^{2}{\widetilde{\omega }}^{4}-({m}^{2}/{r}^{2}+{k}^{2})\widetilde{S},\\ \widetilde{A} & \equiv & \rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})=\rho (\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{+})(\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{-}),\\ \widetilde{S}\, & \equiv & \rho (\gamma p+{B}^{2})({\widetilde{\omega }}^{2}-{\omega }_{{\rm{S}}}^{2})=\rho (\gamma p+{B}^{2})(\omega -{{\rm{\Omega }}}_{{\rm{S}}}^{+})(\omega -{{\rm{\Omega }}}_{{\rm{S}}}^{-}),\end{array}\end{eqnarray} \tag{ 16 }$$

whereas the definitions of the remaining coefficients, as well as the fourth-order ODEs exploited in the numerical analysis, may be found in Appendix A.1.

The expressions for the tangential components ξ_θ and ξ_z in terms of $\chi ^{\prime}$ and χ, explicitly given by Equations (13.93) and (13.94) of Goedbloed et al. (2019), enter the crucial relation of the real frequency part σ of the modes and the solution-averaged Doppler–Coriolis shift,

$$\begin{eqnarray}&&\sigma =\overline{V}=\int \rho \left[r{{\rm{\Omega }}}_{0}| {\boldsymbol{\xi }}{| }^{2}+{\rm{i}}{\rm{\Omega }}({\chi }^{* }{\xi }_{\theta }-\chi {\xi }_{\theta }^{* })\right]{dr}/\int \rho | {\boldsymbol{\xi }}{| }^{2}{rdr},\end{eqnarray} \tag{ 17 }$$

which is real! The first part of the integral is the solution-averaged Doppler shift, to be distinguished from the local Doppler shift Ω₀(r), and the second part is the Coriolis shift of the real part of the frequencies of the modes. The physical significance of equality (17) is that the solution-averaged Doppler–Coriolis-shifted real part of the frequency has to vanish, $\sigma -\overline{V}=0$, for every instability.

In Goedbloed (2009a, 2009b) and Goedbloed (2018a, 2018b), a new method of solving the quadratic eigenvalue problem (6) + BCs was developed based on the fact that the operator U is actually self-adjoint irrespective of whether ξ satisfies the BCs or not. Hence, the expression $\overline{V}$ is real throughout the complex ω-plane for any solution of Equation (6). In particular, one may compute the value of $\overline{V}$ for arbitrary ω by "shooting" from one of the boundaries, say, at r = r₁ in the present case, with solutions of Equation (6) that only satisfy the BCs there. For a given value ν = ν₀, the resulting nonlinear algebraic equation $\sigma -\overline{V}[{\boldsymbol{\xi }}({\boldsymbol{r}};\sigma +{\rm{i}}{\nu }_{0})]=0$ may then be solved by straightforward iteration on the zeros, which gives one or more solutions σ = σ₀(ν₀). The collection of all these points {σ₀(ν₀) + iν₀} of the ω-plane, where the real part of the Doppler–Coriolis-shifted frequency vanishes, has been called the solution path since all eigenvalues have to lie on that path. The actual eigenvalues may then be obtained by a second iteration, along the solution path, to also satisfy the remaining BCs, i.e., at r = r₂ in the present case.

The original solution path method (Goedbloed 2009a, 2009b) has been modified substantially (Goedbloed 2018a, 2018b) by combining the mentioned two steps into a single scheme involving the complementary energy W_com. This is a complex surface integral representing the energy to be injected into or extracted from the plasma such that energy conservation applies for the solution ξ of the spectral equation for any arbitrary value of ω. The zeros of the imaginary part of the complementary energy provide the solution path (where W_com is real), whereas the zeros of the real part provide a second path, called the conjugate path (where W_com is imaginary):

$$\begin{eqnarray}&&\begin{array}{r}{W}_{{\rm{c}}{\rm{o}}{\rm{m}}}[{\boldsymbol{\xi }}({\boldsymbol{r}};\omega )]=0\,\Rightarrow \,\left\{\begin{array}{l}{\rm{I}}{\rm{m}}({W}_{{\rm{c}}{\rm{o}}{\rm{m}}})=0\,(\text{solution path})\\ {\rm{R}}{\rm{e}}({W}_{{\rm{c}}{\rm{o}}{\rm{m}}})=0\,(\text{conjugate path})\end{array}\right..\end{array}\end{eqnarray} \tag{ 18 }$$

It is shown in Goedbloed (2018a) that Equation (18)(a) for the solution path is equivalent to $\sigma -\overline{V}=0$, whereas Equation (18)(b) for the conjugate path is equivalent to satisfying the remaining BC, i.e., at r = r₂ if one "shoots" from the left. The ideal MHD spectrum of unstable modes of rotating equilibria is then found by constructing the Spectral Web, consisting of a superposition of all the curves of the solution path and the conjugate path in the complex ω-plane. The eigenvalues then emerge by computing the intersections of those curves, i.e., where both real and imaginary parts of the complementary energy W_com vanish.

For the present investigation of the instabilities of the cylindrical slice model of an accretion disk, it is most expedient to exploit a mixed integration scheme of ODEs (12) by "shooting" both from the left and from the right and joining the solution χ in some mixing point r_mix in the middle of the interval (r₁, r₂), notably where χ has a large amplitude. One can then adapt the amplitude of the left solution χ^ℓ , satisfying χ^ℓ (r₁) = 0, to that of the right solution χ^r, satisfying χ^r(r₂) = 0, by renormalizing such that χ^ℓ (r_mix) = χ^r(r_mix). For an arbitrary value of ω, the derivative of χ will not be continuous at r_mix, i.e., the total pressure perturbation Π will exhibit a jump there, $[[{\rm{\Pi }}]]\equiv {\rm{\Pi }}({r}_{\mathrm{mix}}^{+})-{\rm{\Pi }}({r}_{\mathrm{mix}}^{-})$. This jump determines the complementary energy:

$$\begin{eqnarray}&&{W}_{\mathrm{com}}^{\mathrm{mix}}=-\tfrac{1}{2}\int {\xi }^{* }\,[[{\rm{\Pi }}]]\,{{dS}}_{\mathrm{mix}}=-\pi {\rm{\Delta }}z{\left({\chi }^{* }[[{\rm{\Pi }}]]\right)}_{{r}_{\mathrm{mix}}}.\end{eqnarray} \tag{ 19 }$$

The spectral web for the mixed solutions then follows by constructing the two paths from the two expressions in Equation (18) for "all" values of ω. Note that the numerical procedure involved is trivially parallel. Moreover, it may be restricted to the particular strip of the ω-plane where physical solutions are expected.

For the purpose of explicit analysis of the accretion disk instabilities, here we will simplify the basic ODE (10) by exploiting the smallness of the parameters of the equilibrium (4). The only place where the pressure, i.e., the parameter β, enters the coefficients of the spectral differential Equation (10) is in the definition of the slow continuum frequency (15), i.e., in the coefficient $\gamma p/{(\gamma p+{B}^{2})}^{1/2}$. We have already seen in Section 2.3 that this coefficient is approximately unity, so that the modes are essentially incompressible. The slow and Alfvén singularities then coalesce, ${{\rm{\Omega }}}_{{\rm{S}}}^{+}\to {{\rm{\Omega }}}_{{\rm{A}}}^{+}$, and similarly for the backward singularities, so that the singularity quotient from expressions (16) becomes

$$\begin{eqnarray}&&\displaystyle \frac{N}{D}\approx -\displaystyle \frac{r\rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})}{{m}^{2}+{k}^{2}{r}^{2}}\equiv -\displaystyle \frac{r\rho (\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{+})(\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{-})}{{m}^{2}+{k}^{2}{r}^{2}}.\end{eqnarray} \tag{ 20 }$$

We will indicate the coalesced continua by $\{{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)\}$ and $\{{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)\}$, i.e., we will not distinguish between slow and Alfvén anymore. This way, the more essential distinction between forward and backward continua is highlighted. We will show that the radial dependence of these continua determines the kind of instabilities that occur, with fundamental differences between the usual axisymmetric MRIs (where the tilde on $\widetilde{\omega }$ may be dropped since Ω₀ = 0) and the present nonaxisymmetric modes (where Ω₀ ≠ 0 is crucial).

In the incompressible approximation, the differential Equation (10) with expressions (A1)–(A6) for the coefficients simplifies to

$$\begin{eqnarray}&&r\displaystyle \frac{d}{{dr}}\left[\displaystyle \frac{\rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})}{{m}^{2}+{k}^{2}{r}^{2}}r\displaystyle \frac{d\chi }{{dr}}\right]-\left[\,\rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})+{\rm{\Delta }}-\displaystyle \frac{4{k}^{2}{\left({B}_{\theta }F+\rho r{\rm{\Omega }}\widetilde{\omega }\right)}^{2}}{({m}^{2}+{k}^{2}{r}^{2})\rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})}-2{mr}\left(\displaystyle \frac{{B}_{\theta }F+\rho r{\rm{\Omega }}\widetilde{\omega }}{r({m}^{2}+{k}^{2}{r}^{2})}\right)^{\prime} \,\right]\,\chi =0,\end{eqnarray} \tag{ 21 }$$

where the function Δ(r) is defined in Equation (A5). Introducing the squared epicyclic frequency ${\kappa }_{{\rm{e}}}^{2}\equiv r({{\rm{\Omega }}}^{2})^{\prime} +4{{\rm{\Omega }}}^{2}$ (which may be negative, in general, but it is positive for equilibrium profiles (4), where ${\kappa }_{{\rm{e}}}^{2}={{\rm{\Omega }}}_{1}^{2}{r}^{-3}$) and exploiting the smallness of the parameters and k⁻¹, the orders of magnitude of the different terms of Equation (21) are given by

$$\begin{eqnarray}&&{m}^{2}\sim 1\ll {k}^{2}{r}^{2},\quad {\rm{\Delta }}\approx -\rho r({{\rm{\Omega }}}^{2})^{\prime} \equiv -\rho ({\kappa }_{{\rm{e}}}^{2}-4{{\rm{\Omega }}}^{2})\sim 1,\quad {B}_{\theta }F\sim {\epsilon }^{2}\ll \rho r{\rm{\Omega }}\widetilde{\omega }\sim 1,\end{eqnarray} \tag{ 22 }$$

whereas the last term of Equation (21) (with the derivative) may be neglected. We then get

$$\begin{eqnarray}&&\displaystyle \frac{r}{\rho }\displaystyle \frac{d}{{dr}}\left[\displaystyle \frac{\rho }{r}({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})\displaystyle \frac{d\chi }{{dr}}\right]-{k}^{2}\left[{\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}-{\kappa }_{{\rm{e}}}^{2}-\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}}\right]\,\chi =0.\end{eqnarray} \tag{ 23 }$$

Since approximations (22) are extremely well satisfied for the accretion disk equilibria, Equation (23) can be considered as the basic ODE for all instabilities of these configurations in the cylindrical representation. For that reason, we will call it the accretion disk ODE in the following. Note that the special equilibrium distributions (4) have not been exploited yet in this equation.

Finally, that the pressure, i.e., compressibility, hardly affects the accretion disk instabilities does not imply that it may also be neglected in the equilibrium Equation (2). Accordingly, all results presented will involve the exact solutions of that equation. Also, whereas the accretion disk ODE (23) will be exploited for all analytical investigations, for the numerical results the solutions of the exact ODEs (10) or (12) will be exploited.

To set the stage for the analysis of the nonaxisymmetric instabilities to be discussed in Section 4, it is instructive to recall the major characteristics of the axisymmetric MRIs following from Equation (23). For the MRIs, since m = 0 so that Ω₀ = 0, the function $\widetilde{\omega }(r)$ reduces to the eigenvalue ω and ω² becomes real, as in static MHD, and waves and instabilities may be distinguished by the simple criterion ω² > 0, or ω² < 0. In that case, ODE (23) reduces to a second-order system to determine the real functions χ(r) and Π(r). (These enormous simplifications may explain the popularity of MRIs to model the anomalous dissipation required for accretion.) A rather general dispersion equation may then be derived by means of a local WKB approximation of the form $\chi (r)=p(r)\exp \left[\pm {\rm{i}}\int q(r)\,{dr}\right]$, where radial variations of the equilibrium variables over distances of the order q⁻¹ are neglected. For m = 0, where ${\omega }_{{\rm{A}}}={{kB}}_{z}/\sqrt{\rho }$, this yields the local dispersion equation,

$$\begin{eqnarray}&&({k}^{2}+{q}^{2}){\left({\omega }^{2}-{\omega }_{{\rm{A}}}^{2}\right)}^{2}-{k}^{2}{\kappa }_{{\rm{e}}}^{2}({\omega }^{2}-{\omega }_{{\rm{A}}}^{2})-4{k}^{2}{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}\,=\,0,\end{eqnarray} \tag{ 24 }$$

having the two solutions

$$\begin{eqnarray}&&{\omega }^{2}={\omega }_{{\rm{A}}}^{2}+\tfrac{1}{2}{\kappa }_{{\rm{e}}}^{2}{E}^{-1}\left[1\pm {\left(1+16E{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}/{\kappa }_{{\rm{e}}}^{4}\right)}^{1/2}\right],\qquad E\equiv 1+{q}^{2}/{k}^{2},\end{eqnarray} \tag{ 25 }$$

where the plus sign corresponds to the frequencies of the stable modes of the epicyclic subspectrum and the minus sign corresponds to the stable and unstable modes of the MRI subspectrum. Expression (25) provides reasonable estimates for the range of q (or the number of zeros of the eigenfunctions, n) from the largest growth rate (smallest n) to the transition to stability. The "exact" unstable eigenvalues obtained from the numerical solutions presented below are in general complex, though.

The numerically computed MRI part of the Spectral Web is illustrated in Figure 3 for a representative configuration, whereas the corresponding eigenfunctions are shown in Figure 4 (these figures are for incompressible MRIs; the compressible counterparts are presented in Goedbloed et al. 2019). This Spectral Web, consisting of the solution path (in red) and the conjugate path (in blue), clearly demonstrates the merits of this method: the 29 unstable eigenvalues are clearly aligned along the solution path that deviates from the imaginary axis (which would be the solution path of the analytical solutions (25)), whereas each of them is located on a kind of pancake of the conjugate path. The epicyclic modes would be situated to the left of the {−ω_A} continuum and to the right of the {+ω_A} continuum in the top panel of Figure 3. They are of no further interest here. On the other hand, only the most global modes (n ≤ 29) of the MRI subspectrum are unstable, whereas the two sequences clustering toward the continua are stable (for simplicity, we will permit the misnomer "MRIs" also for these modes).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The deviations shown in Figure 3 of the numerical solutions of the complex ODEs (12), or rather Equation (A7), from the real analytic solutions (25) have also been found from an extended WKB method with complex frequencies, with proper matching at the turning points (Blokland et al. 2005). This actually provides excellent approximations of the growth rates and eigenfunctions of the MRIs, but it fails to properly describe the approach to the continuum singularities. In general, for frequencies close to the continua, the two singularity coefficients ${\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}$ dominate the behavior of the solutions of Equation (23). For the m = 0 modes, only the static parts { ± ω_A(r)} of the continua enter. For the accretion disk equilibria of Section 2.1, the function ω_A(r) is decreasing so that the stable modes cluster toward the edges ± ω_A2 of the continua at r₂. Expansion for the modes localized about r₂, i.e., ${\omega }_{{\rm{A}}}(r)\approx {\omega }_{{\rm{A}}2}+{\omega }_{{\rm{A}}}^{\prime} (r-{r}_{2})$, transforms Equation (23) for $\widetilde{\omega }=\omega$ into the dominant differential equation

$$\begin{eqnarray}&&\displaystyle \frac{d}{d\tilde{x}}\left[(\tilde{x}+\lambda )\displaystyle \frac{d\chi }{d\tilde{x}}\right]+\displaystyle \frac{{p}^{2}}{\tilde{x}+\lambda }\chi =0,\qquad \tilde{x}\equiv \displaystyle \frac{{r}_{2}-r}{\delta },\quad \lambda \equiv \displaystyle \frac{\omega -{\omega }_{{\rm{A}}2}}{\delta {\omega }_{{\rm{A}}}^{\prime} },\quad p\equiv {\left|\displaystyle \frac{k{\rm{\Omega }}}{{\omega }_{{\rm{A}}}^{\prime} }\right|}_{{r}_{2}},\end{eqnarray} \tag{ 26 }$$

having the solutions

$$\begin{eqnarray}&&\chi ={C}_{1}\sin [p\mathrm{ln}(x+\lambda )]+{C}_{2}\cos [p\mathrm{ln}(x+\lambda )].\end{eqnarray} \tag{ 27 }$$

Satisfaction of the BCs at r₁ and r₂ yields the desired cluster spectra:

$$\begin{eqnarray}&&\lambda \approx \exp (-n\pi /p)\quad \Rightarrow \quad \omega \approx \pm [{\omega }_{{\rm{A}}2}+\delta {\omega }_{{\rm{A}}}^{\prime} \exp (-n\pi /p)].\end{eqnarray} \tag{ 28 }$$

Of course, these solutions only provide the asymptotic behavior for n ≫ 1 and small δ, i.e., for the stable modes. For the global unstable modes (lower n), in particular the MRIs, the approximate WKB expressions (25) are appropriate. Since the pertinent unstable part of the MRI spectrum is located far away from the continua, these instabilities are not affected by the details of the clustering at the continua. This is no longer the case for the nonaxisymmetric modes! For those, we will have to delve much deeper into the subtleties of the continua. This will turn out to have major physical consequences.

Turning now to the nonaxisymmetric instabilities, it would appear logical to present those in the same manner as the axisymmetric MRIs of the previous section, starting from the accretion disk ODE (23), and then to introduce some further approximations to derive an approximate dispersion equation. That approach is bound to fail! A local WKB analysis cannot be applied to these modes because the continuous spectra $\{{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)\}$ and $\{{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)\}$ now are never far away from the actual unstable eigenvalues, but they enter in a crucial way. Hence, these singularities need to be incorporated in any analytical theory of the nonaxisymmetric instabilities of accretion disks. This is done in Section 4.2, on the approximate asymptotic analysis of these instabilities.

4.1.1. Central Role of the Doppler Frequency Range

In the present subsection, we derive the Spectral Webs of the nonaxisymmetric instabilities by numerical means, as described in Section 2.4. Here the search can be restricted to a strip in the complex ω-plane centered about the real Doppler frequencies. As shown in Goedbloed et al. (2004b), the real Doppler frequency range

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{0}\,=\,m{\rm{\Omega }}(r)\equiv {{mv}}_{\theta }(r)/r\end{eqnarray} \tag{ 29 }$$

is not a continuous spectrum in MHD, but it does play a very central role in the description of the nonaxisymmetric instabilities. That range is necessarily located in between the forward and the backward Alfvén continua, since

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{A}}}^{\pm }\equiv {{\rm{\Omega }}}_{0}(r)\pm {\omega }_{{\rm{A}}}(r),\quad {\omega }_{{\rm{A}}}\equiv {{mB}}_{\theta }(r)/r+{{kB}}_{z}(r).\end{eqnarray} \tag{ 30 }$$

We will show that, for the nonaxisymmetric instabilities, the real frequency parts σ of the eigenvalues necessarily lie inside the Doppler range {Ω₀(r)}, so that σ − Ω₀(r) is much smaller than any of the frequencies Ω₀(r), whereas their growth rates ν are also of this smaller order of magnitude. Furthermore, the static Alfvén frequencies ω_A are generally much smaller than the Doppler frequencies Ω₀ so that the three frequency ranges $\{{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)\}$, {Ω₀(r)}, $\{{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)\}$ overlap, at least partly. The closeness of the complex eigenvalues to all three of these real frequency ranges implies that the analysis of the nonaxisymmetric instabilities becomes a very intricate problem with three near singularities. This particular terminology indicates that the two continuum singularities are not actual (which would require ν = 0) but just near because ν is small, whereas singularity of the Doppler frequency would require in addition that ω_A = 0, whereas ω_A is just small compared to Ω₀. Nevertheless, these three near singularities determine the behavior of the solutions completely, as we will see.

(In this discussion, we do not separately mention the slow magnetosonic continua $\{{{\rm{\Omega }}}_{{\rm{S}}}^{\pm }(r)\}$ anymore since they virtually coincide with the Alfvén continua and exactly coalesce in the limit β → ∞. That approximation is extremely well satisfied for the present modes, as well as for the MRIs, since β ≫1 for the equilibria. Of course, in the numerical analysis that approximation is not necessary, but it will be made anyway to simplify the presentation. The changes of the eigenvalues and the Spectral Web by compressibility are small enough to be neglected, unless explicitly mentioned.)

To prove that the real frequencies of the nonaxisymmetric instabilities lie in the Doppler frequency range, we recall Equation (17), which demands that the solution path of these modes is given by the condition that the solution-averaged Doppler–Coriolis-shifted real part of the frequency vanishes, $\sigma -\overline{V}=0$. This condition coincides with the solution path condition (18)(a), which is proved in Goedbloed (2018a). Clearly, without the contribution of the Coriolis term (the second term of expression (17) involving ξ_θ ), this demands that the eigenfunction should have at least two parts with a different sign of the expression σ − Ω₀(r), i.e., σ should be inside the range {Ω₀(r)}. The Coriolis term would spoil this argument, but, for the nonaxisymmetric instabilities, it is small enough to be negligible compared to the Doppler term. This follows by applying the ordering for large mode numbers,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{0}\sim m\gg {\omega }_{{\rm{A}}}\sim k\epsilon .\end{eqnarray} \tag{ 31 }$$

In this ordering, the order of magnitude of the Doppler term is ${ \mathcal O }(m)$, whereas an estimate from the explicit expression for ξ_θ shows that the order of magnitude of the Coriolis term is only ${ \mathcal O }({k}^{-1})$, i.e., a factor (mk)⁻¹ ≪ 1, QED.

For the MRIs, in contrast to the present nonaxisymmetric modes, the above argument does not apply because the Doppler frequencies vanish identically (m = 0 by definition). With that, the rotation profile Ω(r) does not influence the eigenvalues directly but only indirectly through the near-Keplerian equilibrium. (This is actually one of the few cases where compressibility makes a difference. The solution path of the Spectral Web shown here in Figure 3, for incompressible MRIs, deviates from the vertical imaginary axis by a factor of 4 compared to the much smaller deviation for the compressible MRIs that was shown in Figure 13.17 of Goedbloed et al. 2019.)

Consequently, the nonaxisymmetric modes are clearly distinguished from the MRIs by the dominance of the Doppler frequency range, with close "sidebands" of the backward and forward Alfvén frequency ranges, as expressed by Equations (30) and (31). This triple of near singularities is the determining feature of the nonaxisymmetric instabilities. The dominance of the Doppler frequency over the static Alfvén frequency, ∣Ω₀∣ ≫ ∣ω_A∣, makes it appropriate to call these modes SARIs.

Since the Doppler frequency is related to the rotation frequency through Ω₀ = mΩ, these instabilities come in two flavors, viz., counterrotating SARIs (for m < 0) with negative frequencies (σ < 0), and corotating SARIs (for m > 0) with positive frequencies (σ > 0). Furthermore, since ${{\rm{\Omega }}}_{0}^{{\prime} }=m{{\rm{\Omega }}}^{{\prime} }$, the radial frequency profiles of the forward and backward Alfvén singularities are different for the two flavors, viz., increasing for the first kind so that ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})\lt {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$, and decreasing for the second kind so that ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})\gt {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$. This is illustrated in the top panel of Figure 5 and in the top panel of Figure 8 below.

Zoom In Zoom Out Reset image size

Note that whereas the unstable wave packages move at super-Alfvénic speeds, either counter- or corotating with the disk, the disk itself also rotates at super-Alfvénic speeds. The two speeds are well to be distinguished. In a nonlinear context, the super-Alfvénic rotation of the disk will admit shocks, which demand symmetry breaking of the dynamics. This could be initiated by the nonaxisymmetric SARIs, but not by the axisymmetric MRIs. The complex nonlinear interplay of SARIs with rotation of the disk is a subject that deserves separate investigation, clearly far beyond the scope of the present paper.

4.1.2. Counterrotating SARIs

4.1.3. "Virtual Walls"

Three eigenfunctions of counterrotating SARIs, with radial intervals corresponding to the labels (a), (b), and (c) in the top panel of Figure 5, are shown in Figure 6. Anticipating the asymptotic numbering scheme, the orders of the functions in a particular sequence are given in the caption of this figure. The eigenfunctions of Figures 6(a) and (c) are taken from the "inner" sequence, whereas the eigenfunction of Figure 6(b) is taken from the "outer" sequence. The nomenclature "inner" versus "outer" now also becomes clear. The radial eigenfunctions of the "inner" modes, clustering (for ∣ν∣ ↓ 0) to the extremum ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})$ of the forward Alfvén continuum, are intersected by the backward continuum $\{{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)\}$ at some radius ${r}_{2}^{* }\leqslant {r}_{2}$. In the range ${r}_{2}^{* }\leqslant r\leqslant {r}_{2}$, the amplitudes of the "inner" modes are effectively "cut off" by the singularity, so that the radius $r={r}_{2}^{* }$ functions as a kind of "virtual wall." Similarly, the eigenfunctions of the "outer" modes are intersected by the forward continuum $\{{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)\}$ at some radius ${r}_{1}^{* }\geqslant {r}_{1}$, and they cluster to the extremum ${{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{2})$ of the backward Alfvén continuum. In the range ${r}_{1}\leqslant r\leqslant {r}_{1}^{* }$, the amplitudes of the "outer" modes are again "cut off" by the singularity, so that the radius $r={r}_{1}^{* }$ also functions as a kind of "virtual wall." The "cutoff" is not absolute, as illustrated in the blown-up parts of the eigenfunctions shown in Figures 7(a), (b1), and (c1), but their amplitudes are extremely small in the "cutoff" region and, of course, vanish exactly at the prescribed boundary as demanded by the BC. Hence, it is clear that, asymptotically (for small ν), the eigenfunctions and eigenvalues would hardly change if the outer boundary at r = r₂ were moved inward to the point $r={r}_{2}^{* }$ for the "inner" modes, or if the inner boundary were moved outward to the point $r={r}_{1}^{* }$ for the "outer" modes. We will exploit this fact in the approximate analysis of Section 4.2.

Zoom In Zoom Out Reset image size

Actually, the "cutoff" is only absolute in the limit ∣ν∣ ↓ 0, when the singularities produce a split in independent subintervals, very much like the ones introduced in the seminal paper by Newcomb (1960) on the stability of the diffuse linear pinch. The physics of this phenomenon is that the singularity facilitates the induction of skin currents at the position of the singularity, which effectively counteract any motion across that radius. For finite ν, the skin currents are spread out in the "cutoff" region, but they are still quite effective, as shown by the eigenfunction details of Figures 7(a), (b1), and (c1). The same holds for the influence of finite conductivity of the plasma when the induced currents become diffuse. Nevertheless, these "virtual wall" effects quite well explain, e.g., the localization of internal kink instabilities in tokamaks, replacing the actual wall by the virtual one at the singularity. For accretion disks, this effect is more than welcome since the introduction of rigid walls in the numerical calculation of MHD instabilities is anyway a rather questionable procedure that always needs some form of justification. On the other hand, the concept of a "virtual wall" can be justified easily since it is produced by a genuine physical effect, viz., the large conductivity of magnetized plasmas. We will see in Section 5 how two such "virtual walls" (at $r={r}_{1}^{* }$ and $r={r}_{2}^{* }$) effectively produce the localization of "quasi-continuum SARIs" at rather arbitrary locations in the disk.

In Section 5.3 on the Alfvén wave dynamics of the SARIs, we will discuss the relationship between the induction of the perpendicular current density perturbation ${\tilde{j}}_{\perp }$ and the appearance of "virtual walls" more extensively by means of the explicit expression for ${\tilde{j}}_{\perp }$ derived in Appendix A.3. Anticipating that analysis, the distribution of this current component for the two eigenfunctions of Figures 6(b) and (c) is presented in Figures 7(b2) and (c2). In fact, a sharp skin current-like distribution is evident in Figure 7(c2) for the 10th "inner" mode with the very small value of ν, very much like that following from the stability analysis by Newcomb (1960), which is still at the basis of much present tokamak stability theory. However, farther away from marginal stability, the broad oscillatory current distribution of Figure 7(c2) for the fourth "outer" mode shows no sign anymore of a skin current. This points to an important difference between tokamaks and accretion disks that is further discussed in Section 5.3.

Zoom In Zoom Out Reset image size

In summary, the SARIs are unstable waves rotating with super-Alfvénic velocities at the Doppler frequency, evaluated at approximately the middle of the reduced radial interval. The contributions of the forward and backward Alfvén waves, though localized at the end points of the reduced interval, are averaged such that the wave packages as a whole move in phase with that particular value of the Doppler frequency.

4.1.4. Corotating SARIs

The motion in phase with the Doppler frequency is either opposite to the direction of the rotation of the disk, for the counterrotating SARIs just discussed, or in the same direction, for the corotating SARIs. A representative Spectral Web of the latter modes is shown in Figure 8, again with the radial distributions of the three near singularities in the top panel. The continua now decrease as a function of σ so that the positions of the "inner" and "outer" modes in the Spectral Web are reversed with respect to that of Figure 5. With that in mind, there is no need to show the eigenfunctions of the corotating SARIs since they are completely analogous to the counterrotating ones shown in Figure 6.

Zoom In Zoom Out Reset image size

The mode numbers m and k for the counterrotating SARIs of Figure 5 and for the corotating SARIs of Figure 8 are chosen such that the values of the static Alfvén frequency ω_A are approximately the same for the two cases. Consequently, the maximum growth rates are of the same order of magnitude, in agreement with the overestimated limit on the maximum growth rate given by expression (A12) of Appendix A.2. (Incidentally, that expression is also valid for the MRIs of Figure 3, where the maximum growth rate 0.6208 is much closer to the estimated value 0.6988 because the neglected first integral of Equation (A11) is quite small in that case, due to the small average amplitude of ∣χ∣ shown in Figure 4(a).)

Except for the estimate of the maximum growth rate, given in Appendix A.2, not much explicit analysis can be presented on the SARIs. The four distinct cases of the asymptotic approach of the eigenvalues along straight lines to the two cluster points, viz., to ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{2})$ for the counterrotating SARIs of Figure 5 and to ${{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})$ for the corotating SARIs of Figure 8, are completely beyond the kind of local analysis usually exploited for the MRIs. This requires the systematic incorporation of the behavior at the cluster points, as given in the next section. We will take the example of the "inner" corotating SARIs to illustrate that analysis.

Shortcut #1: In summary, we now have showed numerically (using the Spectral Web technique) that nonaxisymmetric SARIs are quite distinct from axisymmetric MRIs: a local WKB analysis does not apply since overlapping Doppler-shifted forward and backward continua interplay. SARIs come in "outer" and "inner," co- and counterrotating flavors. They self-create virtual walls of skin currents, becoming more evident for the higher modes of the infinite sequences of unstable modes. This makes them insensitive to one of our artificial boundaries, and they rotate super-Alfvénically at the Doppler frequency. The next section, which gives a complete analytic treatment of the SARIs, could be skipped on first reading, to continue with our search for quasi-modes that are truly local wave packages described in Section 5.

For an analytic description of the SARIs, we again try to exploit a local approximation where radial variations of χ(r) occur over distances smaller than the scale length of equilibrium variations and consistency implies the assumption of a thin radial shell, δ ≪ 1. We also exploit the smallness of the deviations of the eigenvalues from the rotation frequencies, so that $| \widetilde{\omega }{| }^{2}\sim {\omega }_{{\rm{A}}}^{2}\sim {\epsilon }^{2}\ll {\kappa }_{{\rm{e}}}^{2}\sim 1$. Our basic differential Equation (23) then simplifies to

$$\begin{eqnarray}&&\frac{d}{{dr}}\left[({\tilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})\frac{d\chi }{{dr}}\right]+{k}^{2}\left[{\kappa }_{{\rm{e}}}^{2}-({\tilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})+\frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\tilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}}\right]\,\chi =0,\end{eqnarray} \tag{ 32 }$$

where the coefficients are assumed constant, except for the first and the last one, which produce the singularities, whereas the term −$({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})$ will later be assumed constant with an average magnitude derived below. From the quadratic form corresponding to this ODE, estimates of the eigenvalues may be obtained (see Appendix A.2).

Equation (32) still yields the correct MRI dispersion Equation (24) of the axisymmetric modes (where $\widetilde{\omega }=\omega$) in the WKB approximation, but we now focus on the expansion about the continuum frequencies ${{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)$ when $\widetilde{\omega }\ne \omega$. The basic differential Equation (32) is now essentially complex, hence of fourth order, where the companion system (12), or rather Equation (A7), determines the real and imaginary parts χ₁, Π₁ and χ₂, Π₂ of the eigenfunctions. In this case, the continua usually overlap and the unstable discrete modes have eigenvalues that are close to those continua, as shown by the Spectral Webs of Figures 5 and 8 and the eigenfunctions of Figures 6 and 7. These pictures show that the radial variation of the singularity terms dictates the behavior of the SARIs, where it is essential that both singularities, $\widetilde{\omega }={\omega }_{{\rm{A}}}$ and $\widetilde{\omega }=-{\omega }_{{\rm{A}}}$ , contribute.

4.2.1. Reduction to the Legendre Equation

For definiteness, we restrict the analysis to the sequence of "inner" corotating SARIs, where the various geometrical quantities are illustrated in Figure 9. For a given real frequency component σ of the complex eigenvalue ω, indicated by the vertical line, the intersection of the continuum function ${{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)$ at ${r}_{1}^{* }$ is outside the physical interval, whereas the continuum function ${{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)$ is intersected at ${r}_{2}^{* }$. Hence, the relevant radial interval is reduced to $({r}_{1},{r}_{2}^{* })$. We consider the sequence of modes tending to the cluster point, $\sigma \to {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$, so that the solutions are effectively determined by the two near singularities ${{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$. The left BC is unchanged, χ(r₁) = 0, whereas the right one should be replaced by $\chi ({r}_{2}^{* })=0$ to get the proper boundary value problem for the reduced interval. We introduce a unit radial variable for that interval:

$$\begin{eqnarray}&&\tilde{x}\equiv (r-{r}_{1})/{\delta }^{* },\qquad {\delta }^{* }\equiv {r}_{2}^{* }-{r}_{1},\end{eqnarray} \tag{ 33 }$$

where δ^* is the dimensionless size of the interval (recall that distances have been made dimensionless by dividing through r₁), so that $\tilde{x}=0$ corresponds to r = r₁ and $\tilde{x}=1$ corresponds to $r={r}_{2}^{* }$. From Figure 9 it is clear that the eigenvalue problem would involve a variable interval size δ^* since it would depend on the horizontal distance $\sigma -{{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ to the cluster point. To avoid this complication, ${r}_{2}^{* }$ and δ^* are fixed at the asymptotic value corresponding to $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$ (analogously for the "outer" modes, with interval $({r}_{1}^{* },{r}_{2})$ and clustering toward ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})$).

Zoom In Zoom Out Reset image size

From definition (16) of the Alfvén factor $\tilde{A}$, the singular expression may be expanded about the end points $\tilde{x}=0$ (r = r₁) and $\tilde{x}=1$ ($r={r}_{2}^{* }$) of the interval, i.e.,

$$\begin{eqnarray}&&\begin{array}{c}-\mathrm{for}\,r\,\approx \,{r}_{1}:\qquad {\tilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}\approx \left[\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)\right]\left[{{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})-{{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})\right]\\ \,\approx \,2{\delta }^{* }{\omega }_{{\rm{A}}}({r}_{1}){{\rm{\Omega }}}_{{\rm{A}}}^{-\mbox{'}}({r}_{1})(\tilde{x}-\lambda ),\quad \mathrm{where}\quad \lambda \equiv \frac{\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})}{{\delta }^{* }{{\rm{\Omega }}}_{{\rm{A}}}^{-\mbox{'}}({r}_{1})},\end{array}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\begin{array}{l}-\mathrm{for}\,r\,\approx \,{r}_{2}^{* }:\qquad {\tilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}\approx \left[\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)\right]\left[{{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })-{{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{2}^{* })\right]\\ \,\approx \,-\,2{\delta }^{* }{\omega }_{{\rm{A}}}({r}_{2}^{* }){{\rm{\Omega }}}_{{\rm{A}}}^{+\mbox{'}}({r}_{2}^{* })(\tilde{x}-1-\hat{\lambda }),\quad \mathrm{where}\quad \hat{\lambda }\equiv \frac{\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })}{{\delta }^{* }{{\rm{\Omega }}}_{{\rm{A}}}^{+\mbox{'}}({r}_{2}^{* })}.\end{array}\end{eqnarray} \tag{ 35 }$$

The derivatives ${{\rm{\Omega }}}_{{\rm{A}}}^{\pm ^{\prime} }$ of the continuum frequencies may be transformed into derivatives of the Doppler frequency and, hence, of the rotation frequency, ${{\rm{\Omega }}}_{{\rm{A}}}^{\pm ^{\prime} }\approx {{\rm{\Omega }}}_{0}^{{\prime} }=m{{\rm{\Omega }}}^{{\prime} }$, since ∣ω_A∣ ≪ ∣Ω₀∣. Also, we will assume ∣λ∣ and $| \hat{\lambda }| \ll 1$. Of course, the eigenvalue parameters λ and $\hat{\lambda }$ are related since there is only one eigenvalue ω:

$$\begin{eqnarray}&&\lambda \approx \displaystyle \frac{\sigma -{{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})+{\rm{i}}\nu }{{\delta }^{* }{{\rm{\Omega }}}_{0}^{{\prime} }({r}_{1})},\qquad \hat{\lambda }\approx \displaystyle \frac{{\rm{i}}\nu }{{\delta }^{* }{{\rm{\Omega }}}_{0}^{{\prime} }({r}_{2}^{* })}\quad \Rightarrow \quad {\hat{\lambda }}_{2}={\lambda }_{2}/h,\quad h\equiv \displaystyle \frac{{{\rm{\Omega }}}^{{\prime} }({r}_{2}^{* })}{{{\rm{\Omega }}}^{{\prime} }({r}_{1})}\approx 1-\tfrac{5}{2}{\delta }^{* }.\end{eqnarray} \tag{ 36 }$$

Defining the average of equilibrium quantities at the two end points, $\langle f\rangle \equiv \tfrac{1}{2}\left[f({r}_{1})+f({r}_{2}^{* })\right]$, the two expressions (34) and (35) may be combined into a single approximate expression,

$$\begin{eqnarray}&&{\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}\approx 2{\delta }^{* }\langle | {\omega }_{{\rm{A}}}{{\rm{\Omega }}}_{0}^{{\prime} }| \rangle (\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda }),\end{eqnarray} \tag{ 37 }$$

which is considered to be valid over the whole interval. The absolute signs are introduced in this combined expression (not in the definitions of λ and $\hat{\lambda }$!) to get positive multiplication factors in the basic parameters of the ODE presented below. This implies that that ODE is valid for the present corotating SARIs, where ${{\rm{\Omega }}}_{0}^{{\prime} }\lt 0$, as well as for the counterrotating SARIs, where ${{\rm{\Omega }}}_{0}^{{\prime} }\gt 0$. However, the defining relations (34) and (35) between the complex variables λ and ω would differ in sign for the two cases. Note that 〈f〉〈g〉 ≈ 〈fg〉 since the equilibrium variations are assumed small. With this approximation, the differential Equation (32) transforms into an explicit representation involving all basic parameters:

$$\begin{eqnarray}&&\displaystyle \frac{d}{d\tilde{x}}\left[(\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda })\displaystyle \frac{d\chi }{d\tilde{x}}\right]+\left[\displaystyle \frac{{k}^{2}{\delta }^{* }}{2| m| }\left\langle \displaystyle \frac{{\kappa }_{{\rm{e}}}^{2}}{| {\omega }_{{\rm{A}}}{{\rm{\Omega }}}^{{\prime} }| }\right\rangle -{k}^{2}{\delta }^{* 2}(\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda })+\displaystyle \frac{({k}^{2}/{m}^{2})\langle {\rm{\Omega }}/{{\rm{\Omega }}}^{{\prime} }{\rangle }^{2}}{(\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda })}\right]\chi =0,\end{eqnarray} \tag{ 38 }$$

where the mode number m now appears since ${{\rm{\Omega }}}_{0}^{{\prime} }=m{{\rm{\Omega }}}^{{\prime} }$. Because of approximation (37), the second term of Equation (38) is negligible at the end points and has its maximum contribution in the middle of the interval. Hence, its line average,

$$\begin{eqnarray}&&\overline{(\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda })}\approx {\int }_{0}^{1}\tilde{x}(\tilde{x}-1)\,d\tilde{x}=-\tfrac{1}{6},\end{eqnarray} \tag{ 39 }$$

gives a contribution that can be included in the first term.

Introducing the complex variable $z\equiv (2\tilde{x}-1-\lambda -\hat{\lambda })/(1-\lambda +\hat{\lambda })$, where $\tilde{x}=0$ corresponds to z ≈ − 1 − 2λ and $\tilde{x}=1$ corresponds to $z\approx 1-2\hat{\lambda }$, the singular expression becomes ${\left(\tilde{x}-\lambda )(\tilde{x}-1-\hat{\lambda })\equiv \tfrac{1}{4}(1-\lambda +\hat{\lambda }\right)}^{2}({z}^{2}-1)\approx \tfrac{1}{4}({z}^{2}-1)$. Equation (38) then transforms into a differential equation in terms of the complex variable z:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dz}}\left[({z}^{2}-1)\displaystyle \frac{d\chi }{{dz}}\right]+\left[{u}^{2}+\tfrac{1}{4}-{w}^{2}({z}^{2}-\tfrac{1}{3})+\displaystyle \frac{{v}^{2}}{{z}^{2}-1}\right]\,\chi =0,\end{eqnarray} \tag{ 40 }$$

where the physical parameters u, v, and w are defined by

$$\begin{eqnarray}&&u\equiv {\left[\displaystyle \frac{{k}^{2}{\delta }^{* }}{2| m| }\left\langle \displaystyle \frac{{\kappa }_{{\rm{e}}}^{2}}{| {\omega }_{{\rm{A}}}{{\rm{\Omega }}}^{{\prime} }| }\right\rangle +\tfrac{1}{6}{k}^{2}{\delta }^{* 2}-\tfrac{1}{4}\right]}^{1/2},\qquad v\equiv \left|\displaystyle \frac{2k}{m}\left\langle \displaystyle \frac{{\rm{\Omega }}}{{{\rm{\Omega }}}^{{\prime} }}\right\rangle \right|,\qquad w\equiv \tfrac{1}{2}k{\delta }^{* }.\end{eqnarray} \tag{ 41 }$$

Since we have assumed δ ∼ δ^* ≪ 1, the term ${w}^{2}\left({z}^{2}-\tfrac{1}{3}\right)$, representing the deviation from the average given by Equation (39), is formally negligible with respect to the term ${u}^{2}+\tfrac{1}{4}$. (We will return to this at the end of this section.) Hence, Equation (40) transforms into Legendre's equation solely involving the coefficients u and v, which are related to the canonical Legendre exponents ν_1,2 and μ_1,2 through

$$\begin{eqnarray}&&\nu (\nu +1)=-{u}^{2}-\tfrac{1}{4}\ \ \Rightarrow \ \ {\nu }_{\mathrm{1,2}}\equiv -\tfrac{1}{2}\pm {\rm{i}}u,\quad \mathrm{and}\quad {\mu }^{2}=-{v}^{2}\ \ \Rightarrow \ \ {\mu }_{\mathrm{1,2}}\equiv \pm {\rm{i}}v.\end{eqnarray} \tag{ 42 }$$

(Note the tiny but crucial difference in shape between the Greek symbol ν and the Latin symbol v.) The solutions are the associated Legendre functions ${P}_{-(1/2)\pm {\rm{i}}u}^{\pm {\rm{i}}v}(z)$ and ${Q}_{-(1/2)\pm {\rm{i}}u}^{\pm {\rm{i}}v}(z)$, where the different combinations may be chosen to suit the requirements of a particular boundary value problem. For our case, the pair ${P}_{-(1/2)+{\rm{i}}u}^{{\rm{i}}v}(z)$ and ${P}_{-(1/2)+{\rm{i}}u}^{-{\rm{i}}v}(z)$ turns out to be the most expedient choice, so that a general solution of Equation (40) can be written as

$$\begin{eqnarray}\begin{array}{rcl}\chi & = & {F}_{1}(z)+{{CF}}_{2}(z),\\ {F}_{1} & \equiv & {\rm{\Gamma }}(1-{\rm{i}}v){P}_{-\tfrac{1}{2}+{\rm{i}}u}^{{\rm{i}}v}(z),\quad {F}_{2}\equiv {\rm{\Gamma }}(1+{\rm{i}}v){P}_{-\tfrac{1}{2}+{\rm{i}}u}^{-{\rm{i}}v}(z),\end{array}\end{eqnarray} \tag{ 43 }$$

where the Γ-factors are inserted for (later) algebraic convenience. The behavior of this solution, subjected to BCs (11) with r₂ replaced by ${r}_{2}^{* }$, is to be investigated for values of ω approaching the two singularities z = − 1 and z = 1, as illustrated in Figure 10.

Zoom In Zoom Out Reset image size

The convenience of the Legendre solution (43) is that it provides the full connection between the boundaries of the physical interval at $\tilde{x}=0$ and $\tilde{x}=1$, so that we just need to impose the BCs there to determine the free parameters of the problem. Those are the complex eigenvalue parameter λ and the complex eigenfunction constant C,

$$\begin{eqnarray}&&\lambda =| \lambda | {e}^{{\rm{i}}\alpha },\qquad C=| C| {e}^{{\rm{i}}\varphi },\end{eqnarray} \tag{ 44 }$$

involving four so far unknown constants: ∣λ∣, α, ∣C∣, and φ.

4.2.2. Right Boundary Condition

It is expedient to start with the BC at $\tilde{x}=1$ ($r={r}_{2}^{* }$) since a ready-to-use representation of the associated Legendre functions is given in Olver et al. (2010, Equation 14.3.6) in terms of an algebraic function multiplied by a hypergeometric function that converges for ∣z − 1∣ < 2 (i.e., in the right circle of Figure 10):

$$\begin{eqnarray}&&{[{F}_{1}(z)]}_{{\rm{R}}}\equiv {\left(\displaystyle \frac{z+1}{z-1}\right)}^{\tfrac{1}{2}{\rm{i}}v}F\left(\tfrac{1}{2}+{\rm{i}}u,\tfrac{1}{2}-{\rm{i}}u;1-{\rm{i}}v;\tfrac{1}{2}-\tfrac{1}{2}z\right),\end{eqnarray} \tag{ 45 }$$

whereas ${[{F}_{2}(z)]}_{{\rm{R}}}$ is obtained from ${[{F}_{1}(z)]}_{{\rm{R}}}$ by replacing v by − v. The dominant behavior approaching the singularity, z → 1 for $\tilde{x}=1$, follows from these expressions by evaluating the algebraic factor for $z\approx 1-2\hat{\lambda }$, with $| \hat{\lambda | }\ll 1$, and eliminating the hypergeometric function accordingly using the property $F(a,b;c;\hat{\lambda })\to F(a,b;c;0)\equiv 1$. For the algebraic factors, it is to be noted that the parameters λ₁, λ₂, and ${\hat{\lambda }}_{2}$ are negative, as follows from Equation (36) for the "inner" corotating SARIs, so that the approach of the singularities is from above the branch cut, as illustrated in Figure 10. (Recall that they are positive for the "inner" counterrotating SARIs, so that the singularities would be approached from below the branch cut in that case.) This yields

$$\begin{eqnarray}&&{[{F}_{1}(1-2\hat{\lambda })]}_{{\rm{R}}}\approx {\left({\rm{i}}| {\hat{\lambda }}_{2}| \right)}^{-\tfrac{1}{2}{\rm{i}}v}\approx {{\rm{e}}}^{\tfrac{1}{4}\pi v}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}v\mathrm{ln}| {\hat{\lambda }}_{2}| },\quad {[{F}_{2}(1-2\hat{\lambda })]}_{{\rm{R}}}\approx {\left({\rm{i}}| {\hat{\lambda }}_{2}| \right)}^{\tfrac{1}{2}{\rm{i}}v}\approx {{\rm{e}}}^{-\tfrac{1}{4}\pi v}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}v\mathrm{ln}| {\hat{\lambda }}_{2}| },\end{eqnarray} \tag{ 46 }$$

so that the BC at $\tilde{x}=1$ determines two of the four constants:

$$\begin{eqnarray}&&\begin{array}{l}{[{F}_{1}(1-2\hat{\lambda })]}_{{\rm{R}}}+C{[{F}_{2}(1-2\hat{\lambda })]}_{{\rm{R}}}\approx {{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}\varphi }\left[{{\rm{e}}}^{\tfrac{1}{4}\pi v}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| {\hat{\lambda }}_{2}| +\varphi )}+| C| {{\rm{e}}}^{-\tfrac{1}{4}\pi v}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| {\hat{\lambda }}_{2}| +\varphi )}\right]=0\\ \quad \Rightarrow \quad | C| ={{\rm{e}}}^{\tfrac{1}{2}\pi v},\quad \mathrm{and}\quad \cos (v\mathrm{ln}| {\hat{\lambda }}_{2}| +\varphi )=-1.\end{array}\end{eqnarray} \tag{ 47 }$$

The latter condition yields $v\mathrm{ln}| {\hat{\lambda }}_{2}| +\varphi =(2n^{\prime} +1)\pi$, where $n^{\prime}$ is any (positive or negative) integer. However, since $| {\hat{\lambda }}_{2}|$ is the scaled growth rate of the instabilities, we should demand that it monotonically decreases to zero for positive integers n = 1, 2, ... referring to increasingly oscillatory behavior of the solutions. Writing $n^{\prime} ={n}_{0}-n$, where n₀ is to be determined yet, and noting that − π < α < 0, this gives

$$\begin{eqnarray}&&| {\hat{\lambda }}_{2}| \equiv -h| \lambda | \sin \alpha =h| \lambda | \sin (\pi +\alpha )={{\rm{e}}}^{-[(2n-2{n}_{0}-1)\pi +\varphi ]/v}\quad (n=1,2,\ldots ),\end{eqnarray} \tag{ 48 }$$

where the negative angle α in the complex λ-plane has been replaced by the positive angle π + α (later called $\alpha ^{\prime}$) of the locus of eigenvalues in the ω-plane. There is some arbitrariness in where to start counting. For n = 1, the exponent should be negative to ensure monotonicity for the fastest-growing mode, so that the integer n₀ should at least satisfy

$$\begin{eqnarray}&&{n}_{0}\leqslant \tfrac{1}{2}(1+\varphi /\pi ),\end{eqnarray} \tag{ 49 }$$

whereas a slightly sharper criterion is obtained by applying inequality (A12) of Appendix A.2. Equation (48) implies that the logarithmic sequence of eigenvalue moduli decreases inversely proportional to the singularity exponent v:

$$\begin{eqnarray}&&\mathrm{ln}(| {\lambda }_{n+1}| /| {\lambda }_{n}| )=-2\pi /v,\end{eqnarray} \tag{ 50 }$$

i.e., assuming that the constants α and φ (to be determined yet by the other BC) do not depend on n, which turns out to be the case.

4.2.3. Left Boundary Condition

As stated above, solution (43) provides the full connection between the two boundaries. However, to apply the BC at $\tilde{x}=0$ (r = r₁), representation (45) of the basic solution cannot be used since it is only valid in the right circle of convergence shown in Figure 10, i.e., for ∣z − 1∣ < 2. For the left BC we need an expression that converges in the left circle, i.e., for ∣z + 1∣ < 2. Such a representation may be obtained from Olver et al. (2010, Equation 15.8.4) or Abramowitz & Stegun (1964, Equation 15.3.6), by converting the two hypergeometric functions presented there in terms of z and 1 − z to functions in terms of $\tfrac{1}{2}-\tfrac{1}{2}z$ and $\tfrac{1}{2}+\tfrac{1}{2}z$. This linear transformation from the representation for the right circle of convergence to the one for the left circle is valid for $| \arg (1-z)| \lt \pi$ and $| \arg (1+z)| \lt \pi$, i.e., above as well as below but not on the branch cut. We then obtain the following expression for the basic solution for ∣z + 1∣ < 2:

$$\begin{eqnarray}\begin{array}{rcl}{[{F}_{1}(z)]}_{{\rm{L}}} & \equiv & \displaystyle \frac{{\rm{i}}\pi }{\sinh (\pi v)}{\left(\displaystyle \frac{z+1}{z-1}\right)}^{\tfrac{1}{2}{\rm{i}}v}\left[\displaystyle \frac{{\rm{\Gamma }}(1-{\rm{i}}v)F(\tfrac{1}{2}+{\rm{i}}u,\tfrac{1}{2}-{\rm{i}}u;1+{\rm{i}}v;\tfrac{1}{2}+\tfrac{1}{2}z)}{{\rm{\Gamma }}(1+{\rm{i}}v){\rm{\Gamma }}(\tfrac{1}{2}+{\rm{i}}u-{\rm{i}}v){\rm{\Gamma }}(\tfrac{1}{2}-{\rm{i}}u-{\rm{i}}v)}\right.\\ & & -\left.{\left(\tfrac{1}{2}+\tfrac{1}{2}z\right)}^{-{\rm{i}}v}\,\displaystyle \frac{F(\tfrac{1}{2}+{\rm{i}}u-{\rm{i}}v,\tfrac{1}{2}-{\rm{i}}u-{\rm{i}}v;1-{\rm{i}}v;\tfrac{1}{2}+\tfrac{1}{2}z)}{{\rm{\Gamma }}(\tfrac{1}{2}+{\rm{i}}u){\rm{\Gamma }}(\tfrac{1}{2}-{\rm{i}}u)}\right],\end{array}\end{eqnarray} \tag{ 51 }$$

whereas ${[{F}_{2}(z)]}_{{\rm{L}}}$ is obtained from ${[{F}_{1}(z)]}_{{\rm{L}}}$ by replacing v by − v.

The dominant behavior approaching the singularity, z → − 1 for $\tilde{x}=0$, follows from these expressions by evaluating the algebraic factors for z ≈ − 1 − 2λ, with ∣λ∣ ≪ 1, and eliminating the hypergeometric functions as above. It is to be noted that the −λ occurring inside the brackets of the second algebraic factor is to be evaluated above the branch cut, so that ${\left(\tfrac{1}{2}+\tfrac{1}{2}z\right)}^{-{\rm{i}}v}\approx {(-\lambda )}^{-{\rm{i}}v}={{\rm{e}}}^{\pi v}\cdot {{\rm{e}}}^{\alpha v}{{\rm{e}}}^{-{\rm{i}}v\mathrm{ln}| \lambda | }$. The first solution at $\tilde{x}=0$ then becomes

$$\begin{eqnarray}&&{[{F}_{1}(-1-2\lambda )]}_{{\rm{L}}}\approx {\rm{i}}{f}^{-1}\left[\sqrt{1+{f}^{2}}{{\rm{e}}}^{-\tfrac{1}{2}\alpha v}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}v\mathrm{ln}| \lambda | }{{\rm{e}}}^{-{\rm{i}}\tau }-{{\rm{e}}}^{(\tfrac{1}{2}\alpha +\pi )v}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}v\mathrm{ln}| \lambda | }\right],\end{eqnarray} \tag{ 52 }$$

where f is a function of u and v,

$$\begin{eqnarray}&&f\equiv \sinh (\pi v)/\cosh (\pi u),\end{eqnarray} \tag{ 53 }$$

and the angle τ involves the arguments γ_v and δ_u±v of the Γ-functions,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }}(1+{\rm{i}}v) & = & \sqrt{\pi v/\sinh (\pi v)}\,{{\rm{e}}}^{{\rm{i}}{\gamma }_{v}},\\ {\rm{\Gamma }}(\tfrac{1}{2}+{\rm{i}}u) & = & \sqrt{\pi /\cosh (\pi u)}\,{{\rm{e}}}^{{\rm{i}}{\delta }_{u}},\quad {\delta }_{u}={\gamma }_{2u}-{\gamma }_{u}-2u\mathrm{ln}2,\\ \tau & \equiv & 2{\gamma }_{v}-{\delta }_{u+v}+{\delta }_{u-v}=2{\gamma }_{v}-{\gamma }_{2(u+v)}+{\gamma }_{u+v}+{\gamma }_{2(u-v)}-{\gamma }_{u-v}+4v\mathrm{ln}2.\end{array}\end{eqnarray} \tag{ 54 }$$

The second basic solution ${[{F}_{2}(-1-2\lambda )]}_{{\rm{L}}}$ follows again by replacing v by − v.

With the given approximations, the BC at $\tilde{x}=0$,

$$\begin{eqnarray}&&{[{F}_{1}(-1-2\lambda )]}_{{\rm{L}}}+{{\rm{e}}}^{\tfrac{1}{2}\pi v}{{\rm{e}}}^{{\rm{i}}\varphi }{[{F}_{2}(-1-2\lambda )]}_{{\rm{L}}}=0,\end{eqnarray} \tag{ 55 }$$

may be symmetrized by judicious manipulation of the coefficients to yield

$$\begin{eqnarray}&&\begin{array}{l}\sqrt{1+{f}^{2}}\left[S{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| \lambda | +\varphi )}{{\rm{e}}}^{{\rm{i}}(\varphi +\tau )}-{S}^{-1}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| \lambda | +\varphi )}{{\rm{e}}}^{-{\rm{i}}(\varphi +\tau )}\right]\\ \quad +\ T{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| \lambda | +\varphi )}-{T}^{-1}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| \lambda | +\varphi )}=0,\qquad S\equiv {{\rm{e}}}^{\tfrac{1}{2}(\alpha +\tfrac{1}{2}\pi )v},\quad T\equiv {{\rm{e}}}^{\tfrac{1}{2}(\alpha +\tfrac{3}{2}\pi )v}.\end{array}\end{eqnarray} \tag{ 56 }$$

Substituting relationship (48) from the BC at $\tilde{x}=1$, the exponential factors involving ∣λ∣ may be transformed to eliminate the dependence on n,

$$\begin{eqnarray}&&{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}(v\mathrm{ln}| \lambda | +\varphi )}={{\rm{e}}}^{-{\rm{i}}(n-{n}_{0}-\displaystyle \frac{1}{2})\pi }{{\rm{e}}}^{-{\rm{i}}\eta }={\rm{i}}{\left(-1\right)}^{n-{n}_{0}}{{\rm{e}}}^{-{\rm{i}}\eta },\qquad \eta \equiv \tfrac{1}{2}v\mathrm{ln}[h\sin (\pi +\alpha )],\end{eqnarray} \tag{ 57 }$$

giving a complex equation for the determination of the angles α and φ:

$$\begin{eqnarray}&&\sqrt{1+{f}^{2}}\left[S{{\rm{e}}}^{{\rm{i}}(\varphi +\tau +\eta )}+{S}^{-1}{{\rm{e}}}^{-{\rm{i}}(\varphi +\tau +\eta )}\right]+T{{\rm{e}}}^{{\rm{i}}\eta }+{T}^{-1}{{\rm{e}}}^{-{\rm{i}}\eta }=0.\end{eqnarray} \tag{ 58 }$$

Separating real and imaginary parts of this equation yields two expressions,

$$\begin{eqnarray}\begin{array}{rcl}\sin (\varphi +\tau ) & = & -\displaystyle \frac{\sinh (\tfrac{1}{2}\pi v)\sin (2\eta )}{\sqrt{1+{f}^{2}}\sinh [(\alpha +\tfrac{1}{2}\pi )v]}\,\equiv s(\alpha ),\\ \cos (\varphi +\tau ) & = & -\displaystyle \frac{\sinh [(\pi +\alpha )v]-\sinh (\tfrac{1}{2}\pi v)\cos (2\eta )}{\sqrt{1+{f}^{2}}\sinh [(\alpha +\tfrac{1}{2}\pi )v]}\,\equiv c(\alpha ),\end{array}\end{eqnarray} \tag{ 59 }$$

which, together, determine α and φ in terms of the parameters u and v.

The final equation determining α is obtained by adding the squares of these expressions:

$$\begin{eqnarray}\begin{array}{rcl}P(\alpha ) & \equiv & {\sinh }^{2}(\tfrac{1}{2}\pi v)+{\sinh }^{2}[(\pi +\alpha )v]-2\sinh (\tfrac{1}{2}\pi v)\sinh [(\pi +\alpha )v]\cos [v\mathrm{ln}(h\sin (\pi +\alpha ))]\\ & & -(1+{f}^{2}){\sinh }^{2}[(\alpha +\tfrac{1}{2}\pi )v]=0.\end{array}\end{eqnarray} \tag{ 60 }$$

The angle α should have a value in the range $-\pi \lt \alpha \lt -\tfrac{1}{2}\pi$ to correspond to the range $0\lt \alpha ^{\prime} \equiv \pi +\alpha \lt \tfrac{1}{2}\pi$ of the locus of eigenvalues of the "inner" corotating modes in the ω-plane (see Figure 8). A zero of P(α) in that range is guaranteed since $P\left(-\pi \right)\lt 0$ and $P(-\tfrac{1}{2}\pi )\gt 0$ (for h ≠ 1). Once this zero is found, the final equation determining the parameter φ follows directly from relations (59):

$$\begin{eqnarray}&&\varphi =-\tau +\arctan \left[s(\alpha )/c(\alpha )\right].\end{eqnarray} \tag{ 61 }$$

This completes the asymptotic analysis of the SARIs. The eigenvalue parameters ∣λ∣ and α are determined by Equations (48), or (57), and (60), whereas the eigenfunction parameters ∣C∣ and φ are determined by Equations (47) and (61).

4.2.4. Comparison with the Spectral Web Results

For the "inner" corotating SARIs, corresponding to the right branch of the Spectral Web depicted in Figure 8, the numerical values of the input and the output parameters of the above analytic expressions are as follows:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{input}:\quad \epsilon & = & 0.0141,\quad \beta =100,\quad {\mu }_{1}\,=\,1,\quad \delta =0.1,\quad m\,=\,10,\quad k\,=\,50\\ & \Rightarrow & \quad u\,=\,3.8764,\quad v\,=\,6.9584,\quad h\,=\,0.81078,\quad {\delta }^{* }=0.087526;\\ \mathrm{output}:\quad | C| & = & 5.5841\,\times \,{10}^{4},\quad \tau =5.7643\pi \ \ (\to -0.2357\pi ),\quad \varphi =-1.7747\pi ,\\ | {\lambda }_{1}| & = & 0.35152,\quad \alpha =-0.56194\pi .\end{array}\end{eqnarray} \tag{ 62 }$$

To compare these numbers with the numerical Spectral Web results, the right branch of the Spectral Web of Figure 8 is replotted in double logarithmic coordinates $\mathrm{log}(| {\rm{\Delta }}\sigma | ),\mathrm{log}(\nu )$ in Figure 11, as appropriate for the exponential approach (48) of the cluster point. In the logarithmic representation, ${\rm{\Delta }}\sigma \equiv \sigma -{{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})$ is the real frequency shift with respect to the cluster point. A sequence of nine analytic eigenvalues is shown in green on a dashed line through those points. Note that the identical slope of 45° for the locus of both the analytical and the numerical eigenvalues is accidental, just due to having chosen the same scale for the two logarithmic axes; it has nothing to do with the angle $\alpha ^{\prime}$. However, it is evident from the figure that the approximate eigenvalues (green dots) differ significantly from the "exact" (i.e., numerical) eigenvalues (black dots). In particular, the asymptotic value of the slope of the locus of the numerical "inner" eigenvalues shown in Figure 8 in regular coordinates is $\alpha ^{\prime} \equiv \pi +\alpha =0.36063\pi$, as opposed to the much larger analytical value $\alpha ^{\prime} =0.43806\pi$ following from the numbers (62). This is not surprising in view of the rather crude assumptions that had to be made in approximations (37) and (39) to reduce the problem to the solution of the Legendre equation.

Zoom In Zoom Out Reset image size

The influence of approximation (39), which was necessary to neglect the contribution of the term multiplied by the parameter w ≡ k δ^* in the "extended" Legendre Equation (40), can be estimated by solving that equation numerically for the pertinent value of w. That value is 2.1881 for the present example; the corresponding eigenvalues (green diamonds) shown in Figure 11 are shifted significantly toward the exact eigenvalues (black dots). This suggests that the present approximate asymptotic analysis, based on the assumption of small variation of the equilibrium variables, is justified because the contribution of w vanishes faster for δ^* → 0 than the other contributions according to the orders of the three parameters of the "extended" Legendre equation: v² ∼ 1, ${u}^{2}+\tfrac{1}{4}\sim {\delta }^{* 2}$, w² ∼ δ^*2.

Incidentally, note that the values of u and v, when multiplied by π, produce large values of the hyperbolic functions and the arguments of the Γ functions occurring in the auxiliary angle τ. That angle then has to be reduced to the principal value by adding or subtracting the appropriate number of multiples of 2π. Extreme care had to be taken to maintain the needed accuracy by proper canceling of the oscillatory integrands of the Γ function integrals.

Finally, two qualitative observations on the presented analytic properties of the SARIs are in order. First, the asymmetry of the equilibrium between the two end points of the interval (i.e., the two singularities) was partly accounted for by exploiting the ratio $h\equiv {{\rm{\Omega }}}_{\mathrm{top}}^{{\prime} }/{{\rm{\Omega }}}_{\mathrm{bot}}^{{\prime} }\ne 1$ of the rotation shears, defined in Equation (36). However, the mentioned approximation (37) implies that the index v was forced to be equal for the two singularities, whereas for this particular equilibrium the two indices actually differ, viz., by a factor v_top/v_bot = 1.087. The index v is the only parameter occurring in relation (50), which describes the logarithmic decay of the eigenvalues toward the cluster point. From Figure 11 it is clear that this decay is quite well described by the approximate analytic expression. There is a very surprising feature to this. Recall that Equation (50) was obtained by just applying the BC at $r={r}_{2}^{* }$, so that, in a sense, the logarithmic sequence would be described by the local value of v at that point. But the same expression would also have been obtained by applying the BC at r = r₁, resulting in a logarithmic sequence with another local value of v. From the "exact" (numerical) results shown in Figure 11 this is not the case! Apparently, the two singularities communicate in a nonlocal manner to establish just one sequence with an averaged value of v, as in the approximate analysis. This is another manifestation of the intricate interaction of the three near singularities in the dynamics of the SARIs mentioned at the end of Section 4.1.

Second, the alignment of the numbered sequence of approximate eigenvalues with the "exact" Spectral Web eigenvalues shown in Figure 11 has important consequences for the classification of the different solutions. Whereas numbering of the eigenvalues by counting the number of zeros of the associated eigenfunctions is a standard tool in the spectral analyses of static equilibria (which still applies for the MRIs, as illustrated in Figure 4), this is impossible for the complex eigenfunctions of the SARIs. This is already evident from the fact that the number of zeros (n_z1, n_z2) of the real and the imaginary components χ₁ and χ₂ of χ are usually different, so that there appears to be no unique way of labeling the modes with one number of zeros. For example, for the modes shown in Figure 6 those numbers are (2, 2) for the first "inner" mode, (10, 8) for the fourth "outer" mode, and (23, 20) for the 10th "inner" mode shown. Here "number of zeros" is counted for the full interval, excluding the zero of the BC at r₁ but including the one at r₂, so that this number corresponds with the intuitive notion of "number of lobes" of the eigenfunction (whereas the "number of nodes" is just one less). However, the Legendre analysis provides the unique numbering of eigenvalues and associated eigenfunctions for the SARIs, which was already used in Figure 6. Note that, quite in contrast to the standard spectra with real functions and one singularity, there is no zero-node (one zero) solution! Both sequences start off with solutions with at least two zeros. Moreover, all successive eigenfunctions have about two additional zeros. Apparently, each of the two singularities simultaneously contributes one zero at a time, so that the number of zeros of the eigenfunctions of both sequences is incremented with at least two as n is incremented with one: quite a peculiar aspect of the SARIs!

Of course, counting based on the Legendre analysis will only apply as long as the approximations made in this section are justified. In the following section, we will move far away from these assumptions so that this remnant of spectral counting also will break down.

The SARIs discussed in Section 4 appear to need an edge of one of the Alfvén continua so that they localize about the inner or outer boundary of the accretion disk (see the eigenfunctions of Figure 6). Actually, localization at the inner boundary also happens, for different reasons, for the axisymmetric MRIs discussed in Section 3 (see the eigenfunctions of Figure 4). The analysis of Section 4.2 was necessarily restricted to modes localized in narrow layers of the accretion disk (δ ≪ 1) so that the urgent question to be addressed is, are the SARIs for larger values of δ still localized at the boundaries of the region investigated? Since that region is a kind of necessary artifact of the numerical method used, it would be nice if we also found instabilities that are localized in the middle of the interval, far away from the boundaries of the disk. Hence, we now investigate modes for values of δ that are not small. This analysis has a number of surprises in store for us. It will surpass anything that has been investigated before in the study of the linear MRIs: the best is yet to come!

In Figure 12, the Spectral Web for the corotating SARIs is shown for the same equilibrium and mode numbers investigated in Section 4.2, except that the layer is widened to δ = 1. In fact, concentration of the Spectral Web contours is found for a small range σ ≈ Ω₀(r₀) about the Doppler frequency associated with the fixed matching point r = r_mix = r₀ of the mixed integration scheme described in Section 2.4. If that point is moved to the left where σ ≈ Ω₀(r₂), the "outer modes" clustering toward the forward Alfvén singularity ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})$ appear. If that point is moved to the right where σ ≈ Ω₀(r₁), the "inner modes" clustering toward the backward Alfvén singularity ${{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ appear. This checks with the above remark on the localization of the modes at the edges of the layer. Unfortunately, in the whole range in between those points, the Spectral Webs found just consist of densely packed solution paths and conjugate paths, but they never cross. This appears to suggest that there are simply no instabilities to be found in the middle of the layer!

Zoom In Zoom Out Reset image size

However, for each value of ω in the region of the "condensed" Spectral Web, solutions of the accretion disk eigenvalue problem (10) and (11) are found (as, e.g., illustrated in Figure 13) that do not correspond to discrete eigenvalues, but they do have quite small values of the complementary energy, typically W_com ∼ 10⁻⁵ for this case. They also have the required property of localization in the middle of the interval: they are "confined" by the two singularities, so that they may be characterized as being both "inner" and "outer." Also note that the Spectral Web of Figure 12 is for fixed r₀ = 1.5. Varying that value of r₀, the whole range of the Doppler frequency, from σ = 3.49 to σ = 9.87, is covered with complex values of ω that all are close to genuine eigenvalues. We will now choose parameters such that this behavior is maximized so that modes appear for values of ω that cannot be distinguished from genuine eigenvalues since the measure of distinction, viz., the value of the complementary energy W_com, will be extremely small (less than machine accuracy in most numerical results). For that reason, they may be termed quasi-discrete "eigenvalues." It is important to notice that, in the dynamical environment of an accretion disk, excitation energies of this order of magnitude are readily available so that the distinction between actual eigenvalues and quasi-discrete "eigenvalues" becomes irrelevant.

Zoom In Zoom Out Reset image size

In the next example, the equilibrium and mode numbers are chosen such that the Spectral Web degenerates into a completely fragmented configuration with myriads of small-scale structures of the solution and conjugate paths that never appear to cross though (Figure 14). The values of and ${p}_{1}\equiv \tfrac{1}{2}\beta {\epsilon }^{2}$ are chosen about a factor of 2 smaller than in the equilibrium exploited for Figure 12. The ratio μ₁ of the toroidal over the longitudinal magnetic field is chosen a factor of 10 larger so that a completely different magnetic structure of the accretion disk equilibrium is obtained. Most importantly, the ratio k/m of the longitudinal and toroidal mode numbers is chosen a factor of 10 larger by keeping the value of k the same but decreasing the value of m to 1. The rationale behind this choice will become clear in the analysis of exponential factors in Section 5.2.

Zoom In Zoom Out Reset image size

The fragmentation of the Spectral Web shown in Figure 14 appears to represent a disastrous collapse of the concepts we have presented in this paper. However, quite the opposite is true. The very reason of the fragmentation is not the disappearance of crossings of the two paths, but the fact that whether or not they cross becomes irrelevant because the absolute value ∣W_com∣ of the complementary energy becomes extremely small there. Since the complementary energy distinguishes between quasi-modes (small ∣W_com∣) and actual eigenvalues (vanishing ∣W_com∣), that distinction becomes irrelevant when the value of ∣W_com∣ becomes of the order of the machine accuracy of the computer exploited to derive these results. This is the case for the quasi-modes that are found anywhere in the fragmented region.

At this point, a significant simplification of the method of the Spectral Web is in order. Since the separate solution and conjugate paths and their crossings make no sense anymore for this type of SARI, it is expedient to replace their plots with the contours of ∣W_com∣ instead. Those contours delineate regions of the complex ω-plane where quasi-modes are to be found that cannot be distinguished from genuine eigenmodes when reasonable measures of accuracy, or, better, reasonable measures for a closed system, are exploited. Remember that W_com is the amount of energy that needs to be provided or extracted to bring the system in resonance. That energy vanishes for eigenmodes (when the system is closed), and it is ridiculously small for quasi-modes (when the system is nearly closed). As already mentioned, such small amounts of the complementary energy are readily available from fluctuations in any dynamical system. In Figure 15, the resulting diagram is shown of contours of ∣W_com∣ = 10⁻⁸ and ∣W_com∣ = 10⁻¹² for the same case that produced the Spectral Web of Figure 14. One extension was exploited, though: the matching radius r₀ joining left and right solutions of the ODEs (which was kept fixed for Figure 14) was varied to correspond to the position of the Doppler "resonance," σ − Ω₀(r₀) = 0. This way, a significantly larger area of the ω-plane with quasi-modes is found. That extension could have been exploited for Figure 14 as well, but the corresponding figure is omitted here because the resulting Spectral Web is even more convoluted than the one shown. Moreover, the quasi-continuum representation of the complementary energy contours is just as instructive as the Spectral Web contours.

Zoom In Zoom Out Reset image size

In Figure 16 the quasi-continuum "eigenfunctions" of the SARIs are shown for values of ω corresponding to the labels (a)–(f) of Figure 15. The left column shows that, for fixed growth rate ν, the localization of the modes shifts from the right to the left position in the disk following the value of σ corresponding to the Doppler frequency in the top panel. The right column shows that, for fixed frequency σ, the "eigenfunctions" become ever more oscillatory as the continua on the real axis are approached. Note that, with the exception of the one shown in Figure 16(d), all "eigenfunctions" of Figure 16 have values of ∣W_com∣ ∼ 10⁻¹⁴: the omitted functions Π₁ and Π₂ would show a discontinuity of this order of magnitude at the matching radius r₀ that is completely invisible! For the same reason, the position and the amplitude of that discontinuity, which depend on the choice of r_mix, do not influence the visible shape of those functions (as long as r_mix is not close to ${r}_{1}^{* }$ or ${r}_{2}^{* }$). Consequently and quite conveniently, all these quasi-continuum "eigenfunctions" are obtained by a single shot solution of the accretion disk eigenvalue problems (10) and (11) for arbitrary values of ω in the quasi-continuum range. This is in sharp contrast to all genuine eigenfunctions shown in the previous sections that are obtained by time-consuming iteration for isolated values of ω to get ∣W_com∣ to become small enough to qualify for an eigenvalue. How the single shot "miracle" comes about is described in the analysis of the next section.

Zoom In Zoom Out Reset image size

In summary: We have found a genuine 2D continuum (without any holes!) of quasi-discrete modes in the lower unstable part of the complex ω-plane. This continuum will be called "quasi-continuum" for short. The modes have small jumps in the total pressure perturbation such that the complementary energy is extremely small. This continuum is well to be distinguished from the classical continuous spectrum, where the modes have nonsquare integrable singularities, like in the 2D continuous spectrum found by Lifschitz (1997) for hyperbolic flows. As will be clarified in the next section, the only condition, in addition to the general instability criteria for the SARIs formulated in the previous sections, is that the ratio of the mode numbers is large, ∣k∣/∣m∣ ≫ 1 (recall that all variables have been made dimensionless, so that k is actually $\bar{k}\equiv {{kr}}_{1}$). The superposition of all quasi-continua for all values of k and m satisfying this condition will occupy a vast area of the complex ω-plane!

Shortcut #2: Again, the next section, which provides a full asymptotic analysis of the 2D quasi-continua, could be skipped on first reading. The physical implications of the quasi-continua are picked up in Section 5.3, where we explain that these local instabilities, completely unaware of any artificial BCs, set up intricate electric current distributions in between the forward and backward resonant locations, fully aware of the local magnetic field orientation. This makes them resemble truly local, growing wave packages that surf super-Alfvénically along the disk.

For the analysis of the quasi-continuum SARIs a considerable revision of the analysis of Section 4.2 is in order. The Legendre equation is no longer applicable since the variation of the coefficients of the accretion disk ODE (23) over the interval cannot be approximated anymore as was done in that section. However, the dominant contribution of the Alfvén singularities remains, and with it, some striking similarities of the solutions. To demonstrate that, we will write the solutions of the general spectral equation in the same form as in the first part of Equation (43):

$$\begin{eqnarray}&&\chi ={F}_{1}(z)+{{CF}}_{2}(z),\end{eqnarray} \tag{ 63 }$$

where the complex variable z of Section 4.2 is exploited again. Now, F₁ and F₂ are not Legendre functions, but the solutions of the accretion disk ODE (23), or rather the general spectral differential Equation (10), subject to BCs (11). As we have seen, the latter effectively may be replaced by the BCs

$$\begin{eqnarray}&&\chi ({r}_{1}^{* })=0\ \ (\mathrm{left}),\qquad \chi ({r}_{2}^{* })=0\ \ (\mathrm{right})\end{eqnarray} \tag{ 64 }$$

at the "virtual walls" at $r={r}_{1}^{* }$ and $r={r}_{2}^{\star }$. The positions of these "walls" are determined by the backward and forward Alfvén continuum singularities $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1}^{* })$ and $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$ for corotating modes, similar to that shown in Figure 9 (but with two intersections with the continua, whereas the overlap region ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})\leqslant \sigma \leqslant {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ is now significantly larger than depicted in that figure). For the contra-rotating modes, the positions of ${{\rm{\Omega }}}_{{\rm{A}}}^{-}$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{+}$ are interchanged.

5.2.1. Large and Small Solutions

We again expand the singular factor ${\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}\equiv [\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)][\omega -{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)]$, like in Equations (34) and (35), with ${r}_{1}\to {r}_{1}^{* }$, but of course without the averaging procedure exploited in Section 4.2 that gave the Legendre equation. For this case, we introduce the scaled radial variable $\tilde{x}\equiv (r-{r}_{1}^{* })/{\delta }^{* }$, where ${\delta }^{* }\equiv {r}_{2}^{* }-{r}_{1}^{* }$, generalizing Equation (33), and again relate that variable to the complex variable $z\equiv (2\tilde{x}-1-\lambda -\hat{\lambda })/(1-\lambda +\hat{\lambda })$. Since now $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1}^{* })={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$, the parameters λ and $\hat{\lambda }$ are both imaginary (the parameter α of Equation (44) no longer appears):

$$\begin{eqnarray}&&\lambda =-{\rm{i}}| \lambda | \equiv \displaystyle \frac{{\rm{i}}\nu }{{\delta }^{* }{{\rm{\Omega }}}_{{\rm{A}}}^{-^{\prime} }({r}_{1}^{* })},\qquad \hat{\lambda }=-{\rm{i}}| \hat{\lambda }| \equiv \displaystyle \frac{{\rm{i}}\nu }{{\delta }^{* }{{\rm{\Omega }}}_{{\rm{A}}}^{+^{\prime} }({r}_{2}^{* })}=-{\rm{i}}| \lambda | /h,\quad h\equiv \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{A}}}^{+^{\prime} }({r}_{2}^{* })}{{{\rm{\Omega }}}_{{\rm{A}}}^{-^{\prime} }({r}_{1}^{* })}\approx {\left(\displaystyle \frac{{r}_{2}^{* }}{{r}_{1}^{* }}\right)}^{-5/2}.\end{eqnarray} \tag{ 65 }$$

Hence, the images of the points $\tilde{x}=0$ and $\tilde{x}=1$ will be situated vertically above the points z = − 1 and z = 1 in the z-plane depicted in Figure 10. As before, the approximations are for Alfvén frequency small compared to the Doppler frequency, ∣ω_A∣ ≪ ∣Ω₀∣, but that assumption is not really necessary for the analysis.

We exploit Frobenius expansions about the singular points z = −1 and z = +1. Close to the backward singularity on the left, where $z\approx -1+2(\tilde{x}-\lambda$) and $\tilde{x}\ll 1$, the dominant part of the accretion disk ODE then yields a left solution χ_L built of two independent solutions ${[{F}_{1}]}_{{\rm{L}}}$ and ${[{F}_{2}]}_{{\rm{L}}}$ with imaginary indices:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{dz}}\left[(z+1)\displaystyle \frac{d\chi }{{dz}}\right]+\left[\,\ldots +\displaystyle \frac{{v}_{1}^{2}}{4(z+1)}\right]\chi =0,\qquad {v}_{1}\equiv \left|\displaystyle \frac{2k{\rm{\Omega }}({r}_{1}^{* })}{{{\rm{\Omega }}}_{{\rm{A}}}^{-^{\prime} }({r}_{1}^{* })}\right|\approx \left|\displaystyle \frac{4{{kr}}_{1}^{* }}{3m}\right|\\ \quad \Rightarrow {\chi }_{{\rm{L}}}:\quad {[{F}_{\mathrm{1,2}}]}_{{\rm{L}}}\equiv {\left[-\tfrac{1}{2}(z+1)\right]}^{\pm \tfrac{1}{2}{\rm{i}}{v}_{1}}{f}_{\mathrm{1,2}}(z).\end{array}\end{eqnarray} \tag{ 66 }$$

Close to the forward singularity on the right, where $z\approx 1-2(1-\tilde{x}+\hat{\lambda })$ and $1-\tilde{x}\ll 1$, the accretion disk ODE then yields a right solution χ_R built of two other independent solutions ${[{F}_{1}]}_{{\rm{R}}}$ and ${[{F}_{2}]}_{{\rm{R}}}$ with imaginary indices:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{dz}}\left[(z-1)\displaystyle \frac{d\chi }{{dz}}\right]+\left[\,\ldots +\displaystyle \frac{{v}_{2}^{2}}{4(z-1)}\right]\chi =0,\qquad {v}_{2}\equiv \Space{0ex}{3.0ex}{0ex}\displaystyle \frac{2k{\rm{\Omega }}({r}_{2}^{* })}{{{\rm{\Omega }}}_{{\rm{A}}}^{+^{\prime} }({r}_{2}^{* })}\left|\Space{0ex}{3.8ex}{0ex}\approx \Space{0ex}{3.8ex}{0ex}\right|\displaystyle \frac{4{{kr}}_{2}^{* }}{3m}\Space{0ex}{3.0ex}{0ex}| \\ \qquad \Rightarrow \quad {\chi }_{{\rm{R}}}:\quad {[{F}_{\mathrm{1,2}}]}_{{\rm{R}}}\equiv {\left[\tfrac{1}{2}(z-1)\right]}^{\mp \tfrac{1}{2}{\rm{i}}{v}_{2}}{g}_{\mathrm{1,2}}(z).\end{array}\end{eqnarray} \tag{ 67 }$$

Here f_1,2(z) and g_1,2(z) are regular analytic functions with power series expansions that converge in overlapping regions of the z-plane, similar to the hypergeometric functions exploited in Section 4.2 as illustrated in Figure 10. We do not need any information about these functions other than that they are regular and normalized by imposing the conditions f_1,2( − 1) = 1 and g_1,2(1) = 1. In principle, they could be obtained from the numerical solution of the ODE that is exploited to construct the Spectral Web.

To avoid confusion, it is important to spend a few words on the notation used to distinguish the different quantities. In the parameters ${r}_{1,2}^{* }$ and v_1,2, the subscripts 1 and 2 refer to the left and right end points of the interval $0\leqslant \tilde{x}\leqslant 1$, respectively. For the functions F_1,2, f_1,2, and g_1,2, they refer to the two independent solutions of the ODE, where the upper sign of the powers in Equations (66) and (67) corresponds to the subscript 1 and the lower sign corresponds to the subscript 2. The use of the notation F_1,2 for the left and for the right solutions will be justified shortly.

At this point, it should be observed that the distinguishing feature of the present analysis, compared to that of Section 4.2 exploiting the Legendre functions, is the magnitude of the parameters v₁ and v₂, which is now much larger: v₁ ≫ 1 and v₂ ≫ 1. In terms of the physical variable $\tilde{x}$, this implies that the above solutions may be termed "large" (indicated by the upper sign) and "small" (indicated by the lower sign), exploiting the terminology introduced by Newcomb (1960) for similar singular functions. For example, from Equation (66),

$$\begin{eqnarray}&&{[{F}_{\mathrm{1,2}}]}_{{\rm{L}}}\approx {\left[-{\rm{i}}(| \lambda | -{\rm{i}}\tilde{x})\right]}^{\pm \tfrac{1}{2}{\rm{i}}{v}_{1}}{f}_{\mathrm{1,2}}(z)={{\rm{e}}}^{\pm \tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{\pm \tfrac{1}{2}{v}_{1}\arctan (\tilde{x}/| \lambda | )}{{\rm{e}}}^{\pm \tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}\sqrt{| \lambda {| }^{2}+{\tilde{x}}^{2}}}{f}_{\mathrm{1,2}}(z),\end{eqnarray} \tag{ 68 }$$

where, away from $\tilde{x}=0$, the second exponential factor gives another boost (up or down) of the first factor:

$$\begin{eqnarray}&&{{\rm{e}}}^{\pm \tfrac{1}{2}{v}_{1}\arctan (\tilde{x}/| \lambda | )}=\left\{\,1\ (\mathrm{for}\,\tilde{x}=0)\ \ \to \ \ {{\rm{e}}}^{\pm \tfrac{1}{4}\pi {v}_{1}}\ (\mathrm{for}\,| \lambda | \ll \tilde{x}\ll 1)\,\right\}.\end{eqnarray} \tag{ 69 }$$

Analogously, from Equation (67) for the right solutions at $\tilde{x}=1$,

$$\begin{eqnarray}&&{[{F}_{\mathrm{1,2}}]}_{{\rm{R}}}\approx {\left[{\rm{i}}(| \hat{\lambda }| +{\rm{i}}(1-\tilde{x}))\right]}^{\pm \tfrac{1}{2}{\rm{i}}{v}_{2}}{g}_{\mathrm{1,2}}(z)={{\rm{e}}}^{\pm \tfrac{1}{4}\pi {v}_{2}}{{\rm{e}}}^{\pm \tfrac{1}{2}{v}_{2}\arctan [(1-\tilde{x})/| \hat{\lambda }| ]}{{\rm{e}}}^{\mp \tfrac{1}{2}{\rm{i}}{v}_{2}\mathrm{ln}\sqrt{| \hat{\lambda }{| }^{2}+{\left(1-\tilde{x}\right)}^{2}}}{g}_{\mathrm{1,2}}(z),\end{eqnarray} \tag{ 70 }$$

where, away from $\tilde{x}=1$, the second exponential factor again gives a boost of the first factor,

$$\begin{eqnarray}&&{{\rm{e}}}^{\pm \tfrac{1}{2}{v}_{2}\arctan [(1-\tilde{x})/| \hat{\lambda }| ]}=\left\{\,1\ (\mathrm{for}\,\tilde{x}=1)\ \ \to \ \ {{\rm{e}}}^{\pm \tfrac{1}{4}\pi {v}_{2}}\ (\mathrm{for}\,| \hat{\lambda }| \ll 1-\tilde{x}\ll 1)\,\right\}.\end{eqnarray} \tag{ 71 }$$

Note that, with the parameters used to construct the Spectral Web and quasi-continuum of Figures 14 and 15, the indices themselves are large, v₁ = 84.42 and v₂ = 100.8, but the exponential factors are huge: ${{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{1}}=3.9\times {10}^{57}$ and ${{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{2}}=6.2\times {10}^{68}$, i.e., perfect for expansions! This implies that the solution ${[{F}_{1}]}_{{\rm{L}}}$ is large at $\tilde{x}=0$ and huge away from it, whereas ${[{F}_{2}]}_{{\rm{L}}}$ is small at $\tilde{x}=0$ and tiny away from it. Similarly, the solution ${[{F}_{1}]}_{{\rm{R}}}$ is large at $\tilde{x}=1$ and huge away from it, whereas ${[{F}_{2}]}_{{\rm{R}}}$ is small at $\tilde{x}=1$ and tiny away from it. The crucial conclusion is that the solution that is large at the left end point is also large at the right end point, and similarly for the small solution. Hence, apart from a constant multiplicator, the functions ${[{F}_{1}]}_{{\rm{L}}}$ and ${[{F}_{1}]}_{{\rm{R}}}$ refer to the same independent (large) solution of the ODE, and the functions ${[{F}_{2}]}_{{\rm{L}}}$ and ${[{F}_{2}]}_{{\rm{R}}}$ refer to the other independent (small) solution. This is the justification for choosing the same symbol for the solutions left and right, inspired by the nomenclature of Section 4.2, where the different representations (45) and (51) of the same function F₁(z) already illustrate the crucial property just discussed.

5.2.2. No Discrete Modes

With F₁(z) and F₂(z) two independent solutions of the accretion disk ODE (23), we can now construct the solutions that satisfy the left and the right BC, respectively. For the left solution,

$$\begin{eqnarray}&&{\chi }_{{\rm{L}}}(\tilde{x})=A\left\{{[{F}_{1}(z)]}_{{\rm{L}}}+{C}_{1}{[{F}_{2}(z)]}_{{\rm{L}}}\right\},\quad {C}_{1}=| {C}_{1}| {{\rm{e}}}^{{\rm{i}}{\varphi }_{1}},\end{eqnarray} \tag{ 72 }$$

satisfaction of the BC at $\tilde{x}=0$ yields

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{{\rm{L}}}(0) & = & A\left\{{{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}| \lambda | }+| {C}_{1}| {{\rm{e}}}^{{\rm{i}}{\varphi }_{1}}{{\rm{e}}}^{-\tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}| \lambda | }\right\}=0,\\ & \Rightarrow & | {C}_{1}| ={{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{1}},\qquad {{\rm{e}}}^{{\rm{i}}{\varphi }_{1}}=-{{\rm{e}}}^{{\rm{i}}{v}_{1}\mathrm{ln}| \lambda | },\end{array}\end{eqnarray} \tag{ 73 }$$

i.e., the large contribution is balanced by the small contribution at $\tilde{x}=0$ by means of the huge factor ∣C₁∣. However, away from that point, for $\tilde{x}\gg | \lambda |$, the "large" contribution completely dominates over the "small" one because of the second exponential factor:

$$\begin{eqnarray}&&{\chi }_{{\rm{L}}}(\tilde{x})=A{{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}| \lambda | }\left\{{{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}(\tilde{x}/| \lambda | )}{f}_{1}(z)-{{\rm{e}}}^{-\tfrac{1}{4}\pi {v}_{1}}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}{v}_{1}\mathrm{ln}(\tilde{x}/| \lambda | )}{f}_{2}(z)\right\}.\end{eqnarray} \tag{ 74 }$$

Hence, in the middle of the interval, in particular about the matching point $\tilde{x}={\tilde{x}}_{\mathrm{mix}}$, the left solution is dominated by the independent solution that is "large" at $\tilde{x}=0$. Completely analogously for the right solution,

$$\begin{eqnarray}&&{\chi }_{{\rm{R}}}(\tilde{x})=B\left\{{[{F}_{1}(z)]}_{{\rm{R}}}+{C}_{2}{[{F}_{2}(z)]}_{{\rm{R}}}\right\},\quad {C}_{2}=| {C}_{2}| {{\rm{e}}}^{{\rm{i}}{\varphi }_{2}},\end{eqnarray} \tag{ 75 }$$

satisfaction of the BC at $\tilde{x}=1$ yields

$$\begin{eqnarray}&&| {C}_{2}| ={{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{2}},\qquad {{\rm{e}}}^{{\rm{i}}{\varphi }_{2}}=-{{\rm{e}}}^{-{\rm{i}}{v}_{2}\mathrm{ln}| \hat{\lambda | }},\end{eqnarray} \tag{ 76 }$$

and the solution away from that point, for $1-\tilde{x}\gg | \hat{\lambda }|$, becomes

$$\begin{eqnarray}&&{\chi }_{{\rm{R}}}(\tilde{x})=B{{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{2}}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}{v}_{2}\mathrm{ln}| \hat{\lambda | }}\left\{{{\rm{e}}}^{\tfrac{1}{4}\pi {v}_{2}}{{\rm{e}}}^{\tfrac{1}{2}{\rm{i}}{v}_{2}\mathrm{ln}[(1-\tilde{x})/| \hat{\lambda }| ]}{g}_{1}(z)-{{\rm{e}}}^{-\tfrac{1}{4}\pi {v}_{2}}{{\rm{e}}}^{-\tfrac{1}{2}{\rm{i}}{v}_{2}\mathrm{ln}[(1-\tilde{x})/| \hat{\lambda }| ]}{g}_{2}(z)\right\},\end{eqnarray} \tag{ 77 }$$

so that in the middle of the interval, the right solution is dominated by the independent solution that is "large" at $\tilde{x}=1$.

It remains to match the two solutions (74) and (77) at the matching point $\tilde{x}={\tilde{x}}_{\mathrm{mix}}$, which we will choose to be the Doppler point r = r₀ so that $\sigma ={{\rm{\Omega }}}_{0}({r}_{0})={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1}^{* })={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$. For an eigenvalue, this implies that we have to impose the conditions $[[\chi ]]\equiv {\chi }_{{\rm{R}}}({\tilde{x}}_{\mathrm{mix}})-{\chi }_{{\rm{L}}}({\tilde{x}}_{\mathrm{mix}})=0$ and $[[\chi ^{\prime} ]]\equiv {\chi }_{{\rm{R}}}^{{\prime} }({\tilde{x}}_{\mathrm{mix}})-{\chi }_{{\rm{L}}}^{{\prime} }({\tilde{x}}_{\mathrm{mix}})=0$. Separating the large (exponential) contributions from the small (order unity) contributions, these jumps are determined by

$$\begin{eqnarray}&&[[\chi ]]=B\left({b}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{2}}+{b}_{2}\right)-A\left({a}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{1}}+{a}_{2}\right),\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&[[\chi ^{\prime} ]]=B\left({p}_{1}{b}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{2}}+{p}_{2}{b}_{2}\right)-A\left({p}_{1}{a}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{1}}+{p}_{2}{a}_{2}\right),\quad {p}_{1}\equiv \frac{{F}_{1}^{\mbox{'}}}{{F}_{1}},\quad {p}_{2}\equiv \frac{{F}_{2}^{\mbox{'}}}{{F}_{2}},\end{eqnarray} \tag{ 79 }$$

where we do not need any information about the constants a_1,2, b_1,2, and p_1,2 other than that they are complex and of order unity. Formally, they are supposed to be known through the solutions of the ODE. It is important that, in contrast to the pairs {a₁, a₂} and {b₁, b₂}, which are different owing to the different left and right representations, the logarithmic derivatives p₁ and p₂ are the same left and right since they refer to the same function (apart from the constant amplitude factor). It is trivial now to impose continuity of χ at $\tilde{x}={\tilde{x}}_{\mathrm{mix}}$. From expression (78), this just determines the right amplitude B in terms of the left amplitude A. Inserting this into the jump expression (79) for $\chi ^{\prime}$ and exploiting the smallness of ${{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{1}}$ and ${{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{2}}$, the leading-order term cancels but a smaller part remains:

$$\begin{eqnarray}&&[[\chi ^{\prime} ]]\approx {{Aa}}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{1}}({p}_{1}-{p}_{2})\left[({a}_{2}/{a}_{1}){{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{1}}-({b}_{2}/{b}_{1}){{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{2}}\right].\end{eqnarray} \tag{ 80 }$$

This jump is small compared to the order of magnitude of the function χ itself, $\chi \sim {{Aa}}_{1}{{\rm{e}}}^{\tfrac{1}{2}\pi {v}_{1}}$, because of the small term in square brackets. In general, this jump does not vanish, as would be required for a discrete eigenvalue!

Imposing continuity of the two jumps would be similar to the procedure applied in Section 4.2, except that it was hidden there because the Legendre functions automatically transport the solution continuously from the right to the left. A more fundamental difference is that the eigenvalue was complex there, λ = ∣λ∣e^iα with two free parameters, whereas here λ = − i∣λ∣ with only one free parameter. Whereas continuity of [[χ]] can always be satisfied, in general (except maybe for coincidental values of the parameters), continuity of expression (79) for $[[\chi ^{\prime} ]]$ cannot be satisfied because there are not enough free parameters in this problem. Stated differently, the original eigenvalue problem (63)–(64) will not have a solution, as is also evident from the fact that it contains only one complex constant C, whereas the present analysis contains two essentially different complex constants C₁ and C₂. We wind up with a very negative conclusion: In general, there are no eigenvalues in the region we are looking for, viz., away from the edges of the continua.

5.2.3. Instead, a Continuum of Quasi-modes Emerges

As has been illustrated already by Figures 14–16, the Spectral Web method provides a powerful way out of this conundrum. So far, eigenvalues were determined by the condition that the complementary energy vanishes, W_com = 0. This is associated with the intersections of the solution paths with the conjugate paths, which conspicuously are absent now. For the mixed integration scheme, exploiting the left and right solutions of the accretion disk ODE, the complementary energy is given by expression (19), involving the jump [[Π]] of the total pressure perturbation, defined above Equation (12). This expression assumes that the left and right solutions have been matched by imposing [[χ]] = 0, so that

$$\begin{eqnarray}\begin{array}{rcl}[[{\rm{\Pi }}]] & \equiv & {Q}_{0}[[\chi ^{\prime} ]],\qquad {Q}_{0}\equiv -{(N/D)}_{{r}_{0}}\approx {\left[\displaystyle \frac{\rho ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})}{r({m}^{2}/{r}^{2}+{k}^{2})}\right]}_{{r}_{0}}=-\displaystyle \frac{{\rho }_{0}({\omega }_{{\rm{A}},\ 0}^{2}+{\nu }^{2})}{{r}_{0}({m}^{2}/{r}_{0}^{2}+{k}^{2})}\\ & \Rightarrow & {W}_{\mathrm{com}}=-\pi {\rm{\Delta }}z| {\chi }_{0}{| }^{2}\displaystyle \frac{[[{\rm{\Pi }}]]}{{\chi }_{0}}=\pi {\rm{\Delta }}z| {\chi }_{0}{| }^{2}{Q}_{0}\displaystyle \frac{[[\chi ^{\prime} ]]}{{\chi }_{0}}.\end{array}\end{eqnarray} \tag{ 81 }$$

Inserting $[[\chi ^{\prime} ]]$ from Equation (80) and χ₀ from the second term of Equation (78) then yields

$$\begin{eqnarray}\begin{array}{rcl}{W}_{\mathrm{com}} & \approx & \pi {\rm{\Delta }}z| {\chi }_{0}{| }^{2}{Q}_{0}({p}_{1}-{p}_{2})\left[({a}_{2}/{a}_{1}){{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{1}}-({b}_{2}/{b}_{1}){{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{2}}\right]\\ & \approx & \pi {\rm{\Delta }}z| {\chi }_{0}{| }^{2}{Q}_{0}({p}_{1}-{p}_{2})({a}_{2}/{a}_{1}){{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{1}}\quad \Rightarrow \quad | {W}_{\mathrm{com}}| \sim {{\rm{e}}}^{-\tfrac{1}{2}\pi {v}_{1}},\end{array}\end{eqnarray} \tag{ 82 }$$

where the largest (the one with v₁) of the two small exponentials survives. The very last expression exhibits the key result: asymptotically, the complementary energy is exponentially small! For example, ∣W_com∣ = 6.7 × 10⁻⁵⁰ when ν ≪ 1 in Figures 14 and 15.

In conclusion, the asymptotic analysis confirms the numerical results found in Section 5.1: The Spectral Web "condenses" into a 2D continuum of quasi-discrete modes with eigenfunctions that cannot be distinguished from genuine eigenfunctions because ∣W_com∣ is extremely small, i.e., the discontinuity of $\chi ^{\prime}$ at the matching radius is completely determined by the solution that is exponentially small everywhere.

It should be noted that the asymptotic analysis is for ν → 0, whereas the numerical results shown in Section 5.1 are for finite values of ν, away from the continua on the real axis. For these values of ν, the quasi-modes have their maximum amplitudes in the middle of the interval, as illustrated by all "eigenfunctions" shown in Figures 16. Hence, their oscillatory properties are rather described by expanding the accretion disk ODE about the Doppler point r = r₀ in the middle than about the Alfvén frequencies at ${r}_{1}^{* }$ and ${r}_{2}^{* }$. This yields

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{d\tilde{x}}\left[({\omega }_{{\rm{A}}}^{2}+{\nu }^{2})\displaystyle \frac{d\chi }{d\tilde{x}}\right]-{k}^{2}{\delta }^{* 2}\left[{\kappa }_{{\rm{e}}}^{2}+{\omega }_{{\rm{A}}}^{2}+{\nu }^{2}-\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{A}}}^{2}+{\nu }^{2}}\right]\chi =0\\ \quad \Rightarrow \chi \sim {{\rm{e}}}^{{\rm{i}}q\tilde{x}},\quad q\,=\,k{\delta }^{* }{\left[\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\left({\omega }_{{\rm{A}}}^{2}+{\nu }^{2}\right)}^{2}}-\displaystyle \frac{{\kappa }_{{\rm{e}}}^{2}}{{\omega }_{{\rm{A}}}^{2}+{\nu }^{2}}-1\right]}^{1/2}.\end{array}\end{eqnarray} \tag{ 83 }$$

In the middle of the interval, the wavelength 2π/(q δ^*) corresponds to that of the quasi-mode eigenfunctions depicted in Figure 16. From this expression, oscillatory behavior is limited to

$$\begin{eqnarray}&&\nu \leqslant {\nu }_{\max }={\left[\tfrac{1}{2}\sqrt{{\kappa }_{{\rm{e}}}^{4}+16{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}-\tfrac{1}{2}{\kappa }_{{\rm{e}}}^{2}-{\omega }_{{\rm{A}}}^{2}\right]}^{1/2}\approx | {\omega }_{{\rm{A}}}| \sqrt{4{{\rm{\Omega }}}^{2}/{\kappa }_{{\rm{e}}}^{2}-1}\approx | {\omega }_{{\rm{A}}}| \sqrt{3},\end{eqnarray} \tag{ 84 }$$

which should be an upper limit to the quasi-mode continuum. Not surprisingly, this agrees with expression (A12) of Appendix A.2, except that the parameters are now evaluated at the Doppler point where $\sigma ={{\rm{\Omega }}}_{0}({\tilde{x}}_{0})$. Formally, this yields the same approximate criterion for instability as for the MRIs:

$$\begin{eqnarray}&&0\lt {\omega }_{{\rm{A}}}^{2}\lt 4{{\rm{\Omega }}}^{2}-{\kappa }_{{\rm{e}}}^{2}\approx 3{{\rm{\Omega }}}^{2},\end{eqnarray} \tag{ 85 }$$

where the left inequality guarantees nonvanishing growth rate and the approximation on the right is for the accretion disk equilibrium described in Section 2.1. The important difference with the MRIs is that, for the SARIs, the Alfvén frequency ${\omega }_{{\rm{A}}}\equiv ({{mB}}_{\theta }/r+{{kB}}_{z})/\sqrt{\rho }$ contains both the mode number m and the toroidal field component B_θ .

Expression (84) explains why the upper limit of the quasi-continuum shown in Figure 15 goes up with increasing frequency σ. At ν = 0, the quasi-continuum is limited by the singularity $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})$ on the left (black dash) and by the singularity $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ on the right (red dash). For larger values of ν, it widens in the directions of the two vertical dashed lines at σ = Ω₀(r₂) and σ = Ω₀(r₁), but it is limited from above by the curve that is increasing because ${\nu }_{\max }$ from Equation (84) is smaller on the left than on the right: ω_A(r₂) < ω_A(r₁) for the equilibria investigated.

To our knowledge, the present 2D continuum of quasi-modes has not been described before. It should be considered as a subclass of modes outside the spectrum itself. Quite fundamentally, the spectrum consists of the collection of proper (discrete) and improper (continuum) modes, and any complex frequency outside is part of the so-called resolvent set (see, e.g., Section 6.3.1 of Goedbloed et al. 2019). In a physical context, the resolvent set may be considered as the collection of complex frequencies for which the initial value problem can be solved or, from another perspective, where the driven problem has a finite response. In terms of the differential equations involved, a frequency belongs either to the spectrum, corresponding with solutions to the homogeneous equation, or to the resolvent set, corresponding with solutions to the inhomogeneous equation. This dichotomy is called the Fredholm alternative. It is exactly expressed by the concept of the complementary energy W_com of the Spectral Web, which should be considered as the energy that is needed to bring the open system (described by the inhomogeneous version of the differential equation) into resonance. That energy vanishes for discrete modes, which are natural resonators. What we have found here is a subclass of modes that, strictly speaking, are not discrete modes, but which require just a tiny amount of energy to bring them into resonance. Consequently, for all practical purposes, they cannot be distinguished from genuine discrete modes.

An important physical question remains to be addressed: what causes the extreme localization of the quasi-continuum SARIs about the Doppler point brought about by what we have called "virtual walls" at the positions of the Alfvén/slow singularities? We have argued in Section 4.1, discussing the similar behavior of the eigenfunctions of the regular SARIs shown in Figure 7, that induced currents at the positions of the singularities cause this, in analogy to the stabilization of internal kink modes in tokamaks. To analyze this behavior, it is necessary to compute the distribution of the perpendicular current density perturbation ${\tilde{j}}_{\perp }(r)$ since this yields radial confinement by the Lorentz force. This function may be constructed from the solution pair {χ, Π} of the accretion disk ODE. This requires rather involved algebra, presented in Appendix A.3. Exploiting the resulting expression (A16) for ${\tilde{j}}_{\perp }$, the current density distribution of Figure 7(c2) for a regular SARI with an eigenvalue ω close to the continuum singularities exhibited, in fact, skin currents at the positions of the singularities. However, in the broad current density distribution of the SARI with a larger growth rate shown in Figure 7(b2) skin currents were absent. Nevertheless, that mode is still "confined" by a "virtual wall" at $r={r}_{1}^{* }$. Evidently, this requires explanation beyond that given in Section 4.1.

Turning to the quasi-continuum SARIs, the extreme behavior of these modes, caused by the presence of the near singularities $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1}^{* })$ and $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2}^{* })$, is exhibited, once more, in Figure 17. These singularities are not even actual, or close, since the imaginary growth rate component of the "eigenfrequency" ω is finite, ν ≠ 0. Their importance is shown in the blow-up with a factor of 10¹⁰ in Figure 17(b), showing the "explosions" of the amplitude of the solution at those positions. The additional blow-up, with a factor of 10⁶², in Figure 17(c) shows that the BCs at r = r₁ and r = r₂ are actually satisfied, but also that the amplitude of the solution is completely negligible in the two regions $({r}_{1},{r}_{1}^{* })$ and $({r}_{2}^{* },{r}_{2})$ outside the singularities. How could small amplitudes ∼ 10⁻¹⁰ of the solution at the near singularities, compared to the large ∼1 amplitude at the Doppler position, cause the "confinement" demonstrated by all eigenfunctions of the SARIs shown in this paper?

To answer this question, we study the behavior of the perpendicular current distributions presented in Figures 18 and 19. These figures correspond to the same equilibrium and mode numbers as for the mode shown in Figure 17, with σ = 0.6 and ν = 0.01 for Figure 18, but ν = 0.001 for Figure 19. They demonstrate, more clearly than Figures 7(b2) and (c2) for the regular SARIs, the effective localization of both the "eigenfunction" χ and the perpendicular current density ${\tilde{j}}_{\perp }$ for the quasi-continuum SARIs. The confinement is now well within the singularities at the "virtual walls" at ${r}_{1}^{* }$ and ${r}_{2}^{* }$!

The first panel of Figure 18 shows the perpendicular current density perturbation ${\tilde{j}}_{\perp }$, according to expression (A16). It appears to follow the spatial oscillations of the "eigenfunction" of Figure 17(a) rather closely, implying the excitation of a large number of current filaments of opposite directions. The blow-up with only a factor of 10³ of the same current distribution, shown in Figure 18(b), shows that the dependence on the second derivative of χ (represented by the first term of expression (A16) for ${\tilde{j}}_{\perp }$) causes a significant increase of the magnitude of the current density close to the singularities. However, skin currents at the singularities are absent. Clearly, an additional effect is needed to bridge the gap between the positions of the singularities (at ${r}_{1}^{* }$ and ${r}_{2}^{* }$) and the large-amplitude perturbation within. The opposite directions of the perpendicular current in the filaments imply that each pair of currents will correspond to radial Lorentz forces ${{\boldsymbol{e}}}_{r}\cdot (\tilde{{\boldsymbol{j}}}\times {\boldsymbol{B}})\equiv {\tilde{j}}_{\perp }B$ of opposite directions. Maybe there is a systematic imbalance between the positive and the negative current channels that would produce a net inward Lorentz force? To decide on this, a careful analysis of the average current densities in the current channels was undertaken.

The period-averaged components ${[{\tilde{j}}_{\perp 1}]}_{\mathrm{pav}}$ and ${[{\tilde{j}}_{\perp 2}]}_{\mathrm{pav}}$ shown in Figures 18(c) and 19(c) are constructed as follows. Consider the N spatial periods of oscillation of, e.g., the component ${\tilde{j}}_{\perp 1}$, numbered by n. Each period begins at x = x_b,n where ${\tilde{j}}_{\perp 1}$ becomes positive, it becomes negative in the middle of the interval at x = x_m,n , and it ends at x = x_e,n , where it becomes positive again. Evaluating the positive and negative contributions separately, in search of a systematic difference indicative of localized sheet currents, this yields

$$\begin{eqnarray}&&[{\tilde{j}}_{\perp 1}{]}_{\mathrm{pav}}^{+,n}={\left({x}_{{\rm{m}},n}-{x}_{{\rm{b}},n}\right)}^{-1}{\int }_{{x}_{{\rm{b}},n}}^{{x}_{{\rm{m}},n}}{\left({\tilde{j}}_{\perp 1}\right)}^{+}\,{dx},\qquad [{\tilde{j}}_{\perp 1}{]}_{\mathrm{pav}}^{-,n}={\left({x}_{{\rm{e}},n}-{x}_{{\rm{m}},n}\right)}^{-1}{\int }_{{x}_{{\rm{m}},n}}^{{x}_{{\rm{e}},n}}{\left({\tilde{j}}_{\perp 1}\right)}^{-}\,{dx},\end{eqnarray} \tag{ 86 }$$

where the intervals are centered at, respectively, ${x}_{n}^{+}\equiv \tfrac{1}{2}({x}_{{\rm{b}},n}+{x}_{{\rm{m}},n})$ and ${x}_{n}^{-}\equiv \tfrac{1}{2}({x}_{{\rm{m}},n}+{x}_{{\rm{e}},n})$. Analogous expressions may be derived for $[{\tilde{j}}_{\perp 2}{]}_{\mathrm{pav}}^{\pm ,n}$. Expressions (86) are depicted in Figures 18(c) and 19(c), in black for the positive period-averaged components and in blue for the negative ones. Note that the black curves are defined only at the points $x={x}_{n}^{+}$ and the blue curves only at the points $x={x}_{n}^{-}$, whereas the points in between are obtained by simple interpolation. The curves exactly follow the envelopes of the oscillating current density components depicted in Figures 18(a) and 19(b). More importantly, reflecting the curve of the negative period-averaged components with respect to the x-axis, it precisely coincides with the curve of the positive period-averaged components: there is no systematic difference indicative of localized sheet currents! Moreover, the total positive and negative current perturbations flowing in the plasma exactly cancel (i.e., to machine accuracy): there is no net perpendicular current ${\tilde{j}}_{\perp }$ flowing in the plasma either.

So far, we have only considered the spatial dependence of the different modes. The viewpoint of spatial-period-averaging the perpendicular currents in the different channels yields perfect balance of the radial Lorentz forces. These forces are alternating between inward and outward, depending on the direction of the current in a particular channel, but they appear to yield a rather static picture. However, if we include the temporal behavior of the modes, dictated by their wave dependence e^{iσ t} , the picture becomes quite different. Since the directions of all currents in the different channels reverse sign in the second half-period of the temporal oscillation, automatically running waves are generated that propagate away from the singularities toward the Doppler point. For the corotating SARIs, these waves may be considered to be backward Alfvén waves traveling from $r={r}_{1}^{* }$ to the right and forward Alfvén waves traveling from $r={r}_{2}^{* }$ to the left, respectively. Pounding on the central part with the large current filaments, these Alfvén waves produce a wave pressure that effectively confines the perturbation within the "virtual walls." This resembles the similar mechanism of ponderomotive stabilization of external kink modes in tokamaks (D'Ippoloto 1988), where impinging ion cyclotron waves effectively produce the stabilization corresponding to a close-fitting conducting wall on the plasma, whereas the wall is actually far away. Similarly, the dynamics of the Alfvén waves associated with the SARIs produce "virtual walls" confining the perturbation within the actual boundaries.

It might be objected that singular Alfvén waves only propagate along the magnetic field, not across. For example, Alfvén wave heating occurs by means of excitation by fast magnetosonic waves propagating toward the singularity. Here, the situation is reversed. Since the plasma is approximated to be incompressible, fast waves degenerate into "waves" that propagate the Lorentz forces instantaneously. In a sense, the fast contribution becomes part of the Alfvén waves, as evidenced by the sizeable radial component of all eigenfunctions shown in this paper.

Of course, separating the mode into a central part with the large current filaments and inward-propagating Alfvén waves on the periphery is just a figure of speech since there is no separate plasma column and waves impinging on it, but the mode just keeps itself together by its wavy character represented by e^{iσ t} while exponentially growing with the factor e^{ν t} . Nevertheless, the analogy with stabilization of internal kink modes at marginal stability (ν = 0) and ponderomotive stabilization of external kink instabilities (ν ≠ 0) by waves helps to clarify the different perpendicular current distributions shown in Figures 7(b2) and (c2) for the regular SARIs and in Figures 18(a) and (c) and Figures 19(b) and (c) for the quasi-continuum SARIs. In fact, as shown in Figure 7(c2) for SARIs with a very small growth rate ν, the current distribution is very close to skin currents at the singularities, producing the separation into independent subintervals (Newcomb 1960) responsible for stability of internal kink modes. The wavy character of these currents is already visible in the small-scale oscillations of the peaked current distribution for that case, but it is more evident in the current distribution for the quasi-continuum modes illustrated in Figures 19(b) and (c). For the latter case, though rather close to marginal stability (ν = 0.001), the current distribution is still largest close to the singularities, but it is already spreading out by the oscillations of the current channels. Farther away from the singularities, for ν = 0.01 illustrated in Figure 18, the current distribution has become peaked at the Doppler point, i.e., farthest away from the "virtual walls."

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the same equilibrium parameters as exploited in the Spectral Web of Figure 14 and the corresponding quasi-continuum of Figure 15, a composition is shown in Figure 20 of three quasi-continua of counterrotating SARIs (m < 0) and three quasi-continua of corotating SARIs (m > 0), together with the complete discrete spectrum of MRIs (m = 0) for a fixed value of k. For ∣k∣ ≫ ∣m∣, there are probably also infinite sequences of discrete SARIs clustering at the tips of the continua, but they are swamped by the "sea" of quasi-continua: there is no way to tell them apart. This picture provides a glimpse of how the lower unstable part of the complex ω-plane would become filled with quasi-continua above the real continua $\{{{\rm{\Omega }}}_{{\rm{A}}}^{\pm }\}$ when all mode numbers would be incorporated. Note that the vertical scale of Figure 20 is of the same order of magnitude for the SARIs and the MRIs. This is in agreement with expressions (84) and (A12) for the maximum growth rate, which is similar for the two kinds of instability. Hence, the MRIs may be considered to be degenerate SARIs where the assumption of axisymmetry of the perturbations has eliminated the most important feature for the onset of turbulence, viz., the possible excitation of modes from the two-dimensional "sea" of nonaxisymmetric quasi-continuum SARIs illustrated in Figure 20. The latter modes are localized in all three directions: vertically (large k), toroidally (m ≠ 0 and possibly large), and radially about the Doppler radius. Moreover, because of the occurrence of "virtual walls," they are radially detached from the actual boundaries of the disk so that any narrow annular band of the disk is unstable with respect to SARIs. Except for the vertical localization, none of these features are applicable to the MRIs.

Zoom In Zoom Out Reset image size

We found that the analytical solution of the accretion disk ODE (23) for nonaxisymmetric SARIs is much more intricate than the usual elementary analysis of the axisymmetric MRIs. The reason is that the Doppler-shifted Alfvén/slow continua,

$$\begin{eqnarray}&&\{{{\rm{\Omega }}}_{{\rm{A}}}^{\pm }(r)\}\equiv \{{{\rm{\Omega }}}_{0}(r)\pm {\omega }_{{\rm{A}}}(r)\},\end{eqnarray} \tag{ 87 }$$

completely control the analysis of the SARIs. Recall that the slow continua $\{{{\rm{\Omega }}}_{{\rm{S}}}^{\pm }(r)\}$ are tacitly included since they coalesce with the Alfvén continua in the approximation of incompressibility, which is quite adequate for these modes. Whereas the equilibrium itself is super-Alfvénic for reasonable values of the mode numbers,

$$\begin{eqnarray}&&{\rm{\Omega }}(r)\gt {\omega }_{{\rm{A}}}(r)\qquad (\mathrm{for}\,\mathrm{all}\,r),\end{eqnarray} \tag{ 88 }$$

for MRIs as well as SARIs, the distinguishing feature of SARIs is that the Doppler frames are rotating at super-Alfvénic velocities:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{0}(r)\equiv m{\rm{\Omega }}(r)\gt {\omega }_{{\rm{A}}}(r)\qquad (\mathrm{for}\,\mathrm{all}\,r).\end{eqnarray} \tag{ 89 }$$

Hence, on average, the modes are also propagating with the Doppler velocity, either corotating (m > 0) or counterrotating (m < 0) with respect to the background rotation. This implies breaking of the axisymmetry of the configuration, of importance for the occurrence of shocks. This might be initiated by the SARIs, but that is clearly beyond the scope of the present paper. It also implies that the forward and backward continua overlap, which is instrumental for the occurrence of two infinite sequences of discrete SARIs emanating from the two tips of those continua.

The analysis of the discrete SARIs in Section 4, even for minor inhomogeneity (δ ≪ 1) to stay as closely as possible to the local analysis of the MRIs in Section 3.2, immediately runs into the dominance of the two sets of continuum singularities $\{\sigma -{{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)=0\}$ and $\{\sigma -{{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)=0\}$ on the real axis of the ω-plane. For finite values of the growth rate ν, they are situated relatively far away in the complex ω-plane, yet they turn out to have a major effect on the localization of the modes. At those values of σ, the Alfvén waves emitted from the tips of the continua affect the modes in such a way as to behave like they are confined by a virtual wall inside the actual boundaries. This could be documented by a detailed analysis of the accretion disk ODE, which was transformed into the Legendre equation by analytic continuation of the real radial variable x (≡ r) into the complex z-plane. This produced completely explicit solutions in terms of the Legendre functions (43), propagating the solution from one boundary to the other. For the discrete SARIs, e.g., for a mode from the "inner" sequence emanating from $\sigma ={{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})$ (see Figure 6(a) or Figure 6(c)), one of those boundaries was the actual one, say, r = r₁, whereas the other one corresponded to a virtual wall at $r={r}_{2}^{* }$ (<r₂). The crucial parameter describing the behavior at the singularities is the imaginary exponent iv, where v is defined in Equation (41), which produces solutions that rapidly oscillate at the singularities and grow explosively away from them (all referring to the spatial domain). This behavior is not restricted to the Legendre equation, but it is characteristic for all solutions of the accretion disk Equation (23), and also for the general ODEs (12), at the Alfvén singularities.

Quite unexpectedly, the analysis of the discrete SARIs for large inhomogeneity (δ ∼ 1) in Section 5 runs into the problem of fragmentation of the Spectral Web (Figure 14). Our method of computing eigenvalues by constructing the intersections of the solution path, from the condition $\mathrm{Im}({W}_{\mathrm{com}})=0$, and the conjugate path, from the condition $\mathrm{Re}({W}_{\mathrm{com}})=0$, appears to break down for certain values of the parameters. The two paths fragment into countless little islands, but there are no intersections anymore! The generalization in Section 5.2 of the Legendre analysis of Section 4.2 provides the answer to why this happens. Again, the two independent solutions of the accretion disk ODE can be distinguished as "small" and "large" according to their behavior at the virtual wall end points ${r}_{1}^{* }$ (>r₁) and ${r}_{2}^{* }$ (< r₂), but they have the peculiar property that a solution that is "small" at the left end point is still "small" at the right end point. This implies that the BCs at both end points can be satisfied for arbitrary values of ω in a wide neighborhood of the continuous spectra with just an exponentially small closing error. That error is the absolute magnitude ∣W_com∣ of the complementary energy, which is the energy that would have to vanish if the mode was to qualify as a discrete mode. Hence, these modes form a genuine two-dimensional (extending in frequency σ as well as growth rate ν) continuum of quasi-discrete modes that just require a tiny amount of energy to be brought into resonance. These extended regions in the complex eigenfrequency plane belong to the resolvent set of the operator, and their physical role in the actual solution of initial value calculations (as done by any numerical time integration of the MHD equations) is thus unquestionable. To our knowledge, these quasi-modes have not been described before in the literature, neither in the physics on instabilities nor in the mathematics on spectral theory.

The parameters mentioned for which the quasi-continuum SARIs occur are the generalizations of the imaginary exponent iv of the singular expansions of the Legendre equation, viz., the two exponents iv₁ and iv₂ defined in Equations (66) and (67). For super-Alfvénic Doppler frequencies, these exponents become quite simple and transparent:

$$\begin{eqnarray}&&{v}_{i}\approx \left|\displaystyle \frac{4{{kr}}_{i}^{* }}{3m}\right|\qquad (i=1,2).\end{eqnarray} \tag{ 90 }$$

Consequently, for values of ∣k∣ ≫ ∣m∣, these numbers become large, and the exponential factors in which they occur even become exponentially large, which implies exponentially small complementary energies according to Equation (82). All this occurs according to asymptotic analysis close to the continua (ν → 0). Farther away in the complex ω-plane, numerical analysis provides contour plots of the boundaries of the quasi-continua for ∣W_com∣ < c, where c is typically a small number dictated by what one considers to be small enough to correspond to "eigenfunctions" of quasi-modes that cannot be distinguished from genuine eigenfunctions of discrete modes. This way, for the quasi-continuum SARIs, Spectral Web plotting, as in Figure 14, is superseded by contour plotting of the absolute values of the complementary energy, as in Figure 15 for a single value and in Figure 20 for six values of m.

The quasi-continuum SARIs occur when the mode numbers satisfy the condition ∣k∣ ≫ ∣m∣ and the forward and backward Alfvén continuum frequency ranges ${{\rm{\Omega }}}_{{\rm{A}}}^{+}(r)$ and ${{\rm{\Omega }}}_{{\rm{A}}}^{-}(r)$ overlap. The latter requirement implies distinct conditions on the magnitudes of the mode numbers m and k and the equilibrium parameters , μ₁, and β, viz., ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{1})\lt {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{2})$ for the counterrotating SARIs (m < 0) and ${{\rm{\Omega }}}_{{\rm{A}}}^{+}({r}_{2})\lt {{\rm{\Omega }}}_{{\rm{A}}}^{-}({r}_{1})$ for the corotating SARIs (m > 0). Inserting the explicit equilibria of Section 2.1, these can be combined into the following composite condition for the occurrence of quasi-continuum SARIs:

$$\begin{eqnarray}&&1\ll \displaystyle \frac{| k| }{| m| }\lt \displaystyle \frac{{\epsilon }^{-1}\sqrt{1+{\mu }_{1}^{2}}(1-\tfrac{5}{8}{\epsilon }^{2}\beta )[1-{\left({r}_{1}/{r}_{2}\right)}^{3/2}]-\mathrm{sgn}(m){\mu }_{1}\,[1+{\left({r}_{1}/{r}_{2}\right)}^{3/2}]}{1+{\left({r}_{1}/{r}_{2}\right)}^{1/2}}.\end{eqnarray} \tag{ 91 }$$

When this condition is satisfied, 2D continua of quasi-discrete SARIs emanate from the overlapping continuous spectra $\{{{\rm{\Omega }}}_{{\rm{A}}}^{\pm }\}$ on the real axis of the ω-plane. How far the quasi-continua protrude into the complex ω-plane depends on the value of ∣W_com∣ that one considers small enough to qualify for quasi-modes, and eventually (when ∣W_com∣ is no longer small) on inequality (84). The RHS of inequality (91) is dominated by the factor ⁻¹, which should be large for quasi-modes to be possible. This implies rather small values of the vertical magnetic field B_z , unless there is also a sizable toroidal field component B_θ (i.e., μ₁ ≠ 0 and large). The influence of the parameter β, measuring the magnitude of the kinetic pressure p over the magnetic pressure $\tfrac{1}{2}{B}^{2}$, is rather modest. Even at equipartition (β = 1), quite large regions of quasi-continuum SARIs occur, which has been verified by running the compressible version of ROC. Clearly, the explicit RHS inequality of condition (91) only applies for the particular class of equilibria defined in Section 2.1. Other equilibria will require a modification of this condition, but it is important to notice that all major effects of nonaxisymmetric modes on an equilibrium with both vertical and toroidal components have been incorporated in the present analysis. Hence, the occurrence of quasi-continua will not be invalidated by the ramifications of all the possible flows in accretion disks (Simon & Hawley 2009; Hollerbach et al. 2010), if only they have overlapping continua. This might even apply to pure hydrodynamic disks (Umurhan et al. 2016; Lyra & Umurhan 2019) if the velocity profiles are nonmonotonic enough to yield a flow continuum {Ω₀} folding over onto itself.

For the old question on which linear modes are responsible for the turbulent interaction producing the dissipation that is needed for accretion, it is evident that these 2D continua of nonaxisymmetric quasi-modes, which just require a tiny amount of exciting energy, are extremely relevant. For ever larger ratio of ∣k∣/∣m∣ that exciting energy becomes ever smaller, so that the "sea" of quasi-continuum SARIs occupies ever larger portions of the complex ω-plane (i.e., as long as the above mode number restrictions (91) are respected). Clearly, the dynamical processes in the disk automatically will tune into the mode numbers and associated frequencies for which the response is maximum, thus initiating nonaxisymmetric turbulence. We conclude with the conjecture that, to all probability, the onset of 3D turbulence in accretion disks is not governed by the excitation of discrete axisymmetric MRIs but by the excitation of modes from the two-dimensional continua of quasi-discrete nonaxisymmetric SARIs. The final proof of this conjecture would require nonlinear or quasi-linear analysis beyond this paper.

We hope that this work helps to initiate a new research program of MHD spectroscopy for astrophysical plasma configurations. For example, it is well known that merely including sheared flow in hydro or MHD configurations introduces Kelvin–Helmholtz pathways to instability, while (external or any effective) gravity field introduces Rayleigh–Taylor, Parker, quasi-Parker, and quasi-interchange modes (Goedbloed et al. 2019). Since the basic reference for accretion disk studies is the disk-height-averaged hydrodynamic viscous description from Shakura & Sunyaev (1973), a thorough MHD spectroscopic study of that configuration is called for. The full complexity of including magnetic fields in accretion disk configurations on the modes they support is thus far more complex than previously realized, and besides the well-known axisymmetric MRIs, we have now identified both discrete SARI and quasi-continuum SARI modes as essential new ingredients. We look forward to future spectroscopic studies that should make a causal link to further nonlinear MHD behavior, by solving the full MHD equations with initial conditions targeted to excite a well-known single mode. For the particular context of accretion disks, these follow-up simulations must also quantify the role of the various modes in the angular momentum transport they might realize.

Perhaps the most important new aspect that came into focus through this study is that not all relevant information is to be found in pure eigenmodes of the linearized MHD system. The concept of complementary energy, which can be defined for any complex frequency ω, i.e., the energy required to bring the stationary equilibrium state into resonance with this particular frequency, is herewith established as more central to linear theory than previously realized. Exact eigenmodes, like the countable, discrete, unstable m = 0 MRI, happen to be individual frequencies where this complementary energy vanishes: a perfect resonance. In the cylindrical disk limit studied here, the exact eigenmodes also show infinite sequences of discrete unstable modes, like the m ≠ 0 co- or counterrotating SARIs, in addition to the stable continuous Alfvén (and slow) singular eigenmodes. However, any frequency with a nonzero but otherwise negligible complementary energy should be easily excited and is thus of equal relevance. This was shown to lead to entire 2D regions in the unstable eigenfrequency plane, where, strictly speaking, no modes exist but, for all practical purposes, a tiny perturbation suffices to get unlimited growth of localized wave packages that rotate with essentially the local Doppler frequency. These quasi-continua of SARIs are thus the central new concept elucidated in this work, which may bring a completely new view on turbulence, from the linear eigenmode perspective.

Another viewpoint providing information beyond that obtained from the evolution of the individual eigenmodes comes from nonmodal analysis. In typical hydrodynamical problems, like the solution of the Orr–Sommerfeld equation, it is usually stressed that the relevant linear operator is non-self-adjoint and, hence, that the eigenfunctions are nonorthogonal; see, e.g., Schmid (2007). This permits the construction of linear combinations of eigenmodes that, for a limited time period, may grow faster than the fastest-growing eigenmode proper. This is, rightly, considered to be extremely relevant for the excitation of turbulent motion. The viewpoint of nonmodal analysis is given additional impetus by the extension with the concept of pseudospectrum, where the linear operator is extended with an arbitrary operator of small norm ; see Trefethen & Embree (2005). This yields contours in the complex ω-plane that shrink to the spectrum proper in the limit → 0. This is quite relevant for the excitation of turbulence: the pseudospectrum may even extend into the unstable (upper) part of the ω-plane when the spectrum proper is restricted to the stable (lower) part. All this, evidently, evokes two questions: (1) Can the present theory for the SARIs be extended with a nonmodal analysis? (2) Are the contours of the absolute value ∣W_com∣ of the complementary energy, delineating the quasi-continuum SARIs, just another way of representing a pseudospectrum?

Before addressing these questions, it is useful to highlight the essential differences between the mentioned theory of the hydrodynamical problems and our analysis of the SARIs. The first one, in general, refers to nonconservative systems described by a non-self-adjoint linear operator having nonorthogonal eigenmodes. In our analysis, the system is conservative, and the occurring operators G and U are self-adjoint; nevertheless, the eigenmodes are also nonorthogonal! Introducing eigenmodes ξ _α and ξ _β corresponding to eigenvalues ω_α and ω_β of the eigenvalue problem posed by the spectral differential Equation (6) + BCs, this property has been demonstrated in Goedbloed et al. (2019), Equations (12.97) and (12.98), which we here reproduce for the convenience of the reader:

$$\begin{eqnarray}&&\langle {{\boldsymbol{\xi }}}_{\beta },{\rho }^{-1}{\boldsymbol{G}}({{\boldsymbol{\xi }}}_{\alpha })\rangle ={\omega }_{\beta }^{* }{\omega }_{\alpha }\langle {{\boldsymbol{\xi }}}_{\beta },{{\boldsymbol{\xi }}}_{\alpha }\rangle ,\qquad \langle {{\boldsymbol{\xi }}}_{\beta },U{{\boldsymbol{\xi }}}_{\alpha }\rangle =\tfrac{1}{2}({\omega }_{\beta }^{* }+{\omega }_{\alpha })\langle {{\boldsymbol{\xi }}}_{\beta },{{\boldsymbol{\xi }}}_{\alpha }\rangle .\end{eqnarray} \tag{ 92 }$$

Clearly for nonvanishing Doppler–Coriolis operator U, the RHS inner products do not vanish, so that the eigenmodes are nonorthogonal. However, the crucial feature of our approach is that the eigenvalue problem is nonlinear: the solutions are eigenfunctions neither of G nor of U, but of the specific nonlinear (with respect to the eigenvalue parameter ω) combination occurring in Equation (6). For the same reason, the question (usually dictated from quantum mechanical contexts) whether the operators G and U commute or not is irrelevant: the two operators have different roles to play in the problem. In particular, the most significant distinction with all other approaches of the spectral problem is probably the focus it provides on the central importance of the Doppler–Coriolis operator U describing the shear flow that occurs in almost all astrophysical plasmas (e.g., on Rossby waves in the Sun: Zaqarashvili et al. 2007, 2010; Dikpati et al. 2020; on different extensions of the MRI: Rüdiger et al. 2007a, 2007b; Kitchatinov & Rüdiger 2010; Hollerbach et al. 2010; on shear flow instabilities: Heifetz et al. 2015; finally, on turbulence in protoplanetary disks associated with the critical layers where the Doppler-shifted frequency of a Rossby wave coincides with the Brunt–Vaisälä frequency: Umurhan et al. 2016). Having established the important difference with respect to self-adjointness, it remains to notice though that the nonmodal analysis does not really depend on that, but only on the fact that the eigenmodes are nonorthogonal. Hence, the answer to the first question is affirmative: in principle, it should be possible to extend the present theory of the SARIs also with a nonmodal analysis. As a pertinent example, Squire & Bhattacharjee (2014) and also Singh Bhatia & Mukhopadhyay (2016) actually present a nonmodal analysis of the MRIs, the spectral analysis of which we have shown to be one of the solutions of Equation (6). However, these nonmodal analyses crucially depend on WKB approximations and a local dispersion equation within a shearing box model, which we have shown in Section 4.1 to be inadequate to describe the SARIs. To properly describe the discrete SARIs, at least the rapid oscillations approaching the continuum singularities should be incorporated, whereas for the quasi-continuum SARIs also the extreme connection between the two singularities, permitting continua of modes with large amplitudes at the Doppler frequencies, should be represented. Yet the physics of the nonmodal growth of the MRIs undoubtedly carries over to the SARIs, and nature will know how to exploit that in the excitation of turbulence, irrespective of our inability to solve such a complicated problem.

With respect to the second question, on the similarity of the contours of ∣W_com∣ and the pseudospectrum contours, it is clear that both methods are extensions to escape the narrow confines of the spectrum of discrete eigenvalues. They do that in entirely different directions, though. The ∣W_com∣ contours enlarge the class of permitted solutions by relaxing the condition of continuity of the normal component of the displacement vector ξ, whereas the pseudospectrum contours enlarge the differential operator itself. Both ∣W_com∣ and represent a kind of measure of the distance to the actual eigenvalues, justified on physical grounds by the observation that the effects of external excitation and/or dissipation of the system (no matter how small) eventually are to be taken into account. This is actually what the nonmodal analysis also is about, so that it is not surprising that the review on nonmodal stability theory by Schmid (2007) contains an extensive discussion on external forcing, like time-dependent flows, stochastic forcing, flows in complex geometries, etc. For the quasi-continuum modes, here reported for the first time, it is not so clear that the pseudospectral method could also find those since then W_com is not literally a distance to a discrete mode (there are none in a wide neighborhood) but just an insignificant difference with such modes. On the other hand, it should also be pointed out that the present analysis, based on the Frieman–Rotenberg spectral Equation (6), is strictly limited to ideal MHD, so that all the different dissipation mechanisms that occur in astrophysical plasmas cannot be incorporated. Fortunately, a multiplicity of such extensions has already been implemented in our finite-element code Legolas (Claes et al. 2020). This code has verified most results on the SARIs reported here, including the fragmentation of the Spectral Web in the quasi-continuum range. We hope to report in the near future on the powerful broadening of scope by means of the ROC and Legolas codes, operated in tandem.

In conclusion: Any small perturbation of the equilibrium will, in the first instance, excite a huge collection of the quasi-continuum SARIs at each point of the accretion disk since the condition for instability is satisfied everywhere, whereas the modes are truly local in all three directions. Moreover, the breaking of the axisymmetry by these modes also implies possible dynamo action and excitation of MHD shocks. Thus, the crucial presence of a magnetic field, as correctly highlighted in the MRI analysis of Balbus & Hawley (1991), opens up an enormous potential of additional new dynamical pathways to turbulence, necessary to explain why accretion occurs at all and why the accreted magnetized plasma is accelerated to powerful jets. Which of these pathways is chosen by nature is open to future research.

The singularity coefficients N and D of the cylindrical spectral Equations (10) and (12) were given in Equation (16). For the convenience of the reader we here reproduce the definitions of the remaining coefficients from Goedbloed (2018b),

$$\begin{eqnarray}&&A\,\equiv \,\displaystyle \frac{1}{r}\left(\widetilde{A}+{\rm{\Delta }}\right),\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&B\,\equiv \,-\displaystyle \frac{4}{r}\left\{(\rho {\widetilde{\omega }}^{2}-{k}^{2}\gamma p){P}^{2}-\rho {\widetilde{\omega }}^{2}R+\left(m\rho {\widetilde{\omega }}^{2}P-{r}^{2}{h}^{2}Q\right){\rm{\Lambda }}+\tfrac{1}{4}{r}^{2}{h}^{2}\widetilde{A}{{\rm{\Lambda }}}^{2}\right\},\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&C\,\equiv \,\displaystyle \frac{2}{{r}^{2}}\left\{m\widetilde{S}P-{r}^{2}\rho {\widetilde{\omega }}^{2}\left(Q-\tfrac{1}{2}\widetilde{A}{\rm{\Lambda }}\right)\right\},\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&E\,\equiv \,-\displaystyle \frac{1}{{r}^{2}}\left\{\widetilde{A}\widetilde{S}\left(\widetilde{A}+{\rm{\Delta }}\right)-4\widetilde{S}{P}^{2}+4{r}^{2}{\left(Q-\tfrac{1}{2}\widetilde{A}{\rm{\Lambda }}\right)}^{2}\right\},\end{eqnarray} \tag{ A4 }$$

where h² ≡ m²/r² + k² and abbreviations have been introduced for two equilibrium functions, Δ(r) and Λ(r), and three perturbation functions, $P(r;\widetilde{\omega })$, $Q(r;\widetilde{\omega })$, and $R(r;\widetilde{\omega })$, as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}\,\equiv \,r({B}_{\theta }^{2}/{r}^{2}-\rho {{\rm{\Omega }}}^{2})^{\prime} +\rho ^{\prime} {{GM}}_{* }/{r}^{2},\qquad {\rm{\Lambda }}\equiv \rho ({{\rm{\Omega }}}^{2}-{{GM}}_{* }/{r}^{3}),\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&P\,\equiv \,({B}_{\theta }/r)F+\rho {\rm{\Omega }}\widetilde{\omega },\qquad Q\equiv ({B}_{\theta }/r)\left[{FP}+({B}_{\theta }/r)\widetilde{A}\right],\qquad R\equiv {\left({B}_{\theta }/r\right)}^{2}\left[2{mP}-{h}^{2}({B}_{\theta }^{2}+\rho {r}^{2}{{\rm{\Omega }}}^{2})\right].\end{eqnarray} \tag{ A6 }$$

All coefficients are complex through $\widetilde{\omega }\equiv \sigma -{{\rm{\Omega }}}_{0}+{\rm{i}}\nu$, which is also a function of r through the Doppler shift Ω₀(r).

For complex ω, the system of first-order ODEs (12) splits into equations for the real and imaginary components:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\chi }_{1}^{\prime} \\ {\chi }_{2}^{\prime} \\ {{\rm{\Pi }}}_{1}^{\prime} \\ {{\rm{\Pi }}}_{2}^{\prime} \end{array}\right)+\left(\begin{array}{cccc}{\hat{C}}_{1} & -{\hat{C}}_{2} & {\hat{D}}_{1} & -{\hat{D}}_{2}\\ {\hat{C}}_{2} & {\hat{C}}_{1} & {\hat{D}}_{2} & {\hat{D}}_{1}\\ {\hat{E}}_{1} & -{\hat{E}}_{2} & -{\hat{C}}_{1} & {\hat{C}}_{2}\\ {\hat{E}}_{2} & {\hat{E}}_{1} & -{\hat{C}}_{2} & -{\hat{C}}_{1}\end{array}\right)\left(\begin{array}{c}{\chi }_{1}\\ {\chi }_{2}\\ {{\rm{\Pi }}}_{1}\\ {{\rm{\Pi }}}_{2}\end{array}\right)=0,\end{eqnarray} \tag{ A7 }$$

where $\hat{C}\equiv C/N$, $\hat{D}\equiv D/N$, and $\hat{E}\equiv E/N$. The real and imaginary components ${\hat{C}}_{\mathrm{1,2}}$, ${\hat{D}}_{\mathrm{1,2}}$, and ${\hat{E}}_{\mathrm{1,2}}$ follow by straightforward expansion of expressions (16), (A3), and (A4).

Estimates of the complex frequencies of the SARIs may be obtained by deriving a quadratic form corresponding to the accretion disk ODE (23), or the approximated form (32),

$$\begin{eqnarray}&&\int ({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})| {\chi }^{{\prime} }{| }^{2}{dr}-{k}^{2}\int \left[{\kappa }_{{\rm{e}}}^{2}-({\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2})+\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\widetilde{\omega }}^{2}-{\omega }_{{\rm{A}}}^{2}}\right]\,| \chi {| }^{2}{dr}=0,\end{eqnarray} \tag{ A8 }$$

where the boundary term resulting from integration by parts has been canceled by applying the BCs. Writing $\widetilde{\omega }\equiv \tilde{\sigma }+{\rm{i}}\nu$, where $\tilde{\sigma }\equiv \sigma -{{\rm{\Omega }}}_{0}$, and splitting this equation into real and imaginary parts yields

$$\begin{eqnarray}&&\int ({\tilde{\sigma }}^{2}-{\nu }^{2}-{\omega }_{{\rm{A}}}^{2})| {\chi }^{{\prime} }{| }^{2}{dr}-{k}^{2}\int \left[{\kappa }_{{\rm{e}}}^{2}-({\tilde{\sigma }}^{2}-{\nu }^{2}-{\omega }_{{\rm{A}}}^{2})+\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}({\tilde{\sigma }}^{2}-{\nu }^{2}-{\omega }_{{\rm{A}}}^{2})}{{\left({\tilde{\sigma }}^{2}-{\nu }^{2}-{\omega }_{{\rm{A}}}^{2}\right)}^{2}+4{\tilde{\sigma }}^{2}{\nu }^{2}}\right]\,| \chi {| }^{2}{dr}=0,\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&\nu \left\{\int \tilde{\sigma }| {\chi }^{{\prime} }{| }^{2}{dr}+{k}^{2}\int \tilde{\sigma }\left[1+\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\left({\tilde{\sigma }}^{2}-{\nu }^{2}-{\omega }_{{\rm{A}}}^{2}\right)}^{2}+4{\tilde{\sigma }}^{2}{\nu }^{2}}\right]\,| \chi {| }^{2}{dr}\right\}=0.\end{eqnarray} \tag{ A10 }$$

From the latter equation it follows directly that, for instabilities (ν ≠ 0), the real Doppler-shifted frequency $\tilde{\sigma }$ should change sign on the interval, as discussed in Section 4.1. From the first equation it follows directly that there is no solution for ν → ∞ , so that there should be a maximum growth rate ${\nu }_{\max }$. An estimate is obtained by noting that, from the previous argument, the approximation ${\tilde{\sigma }}^{2}\ll {\nu }^{2}+{\omega }_{{\rm{A}}}^{2}$ should hold for the integrands for the most global modes corresponding to ${\nu }_{\max }$. The first equation then simplifies to

$$\begin{eqnarray}&&\int ({\nu }^{2}+{\omega }_{{\rm{A}}}^{2})| {\chi }^{{\prime} }{| }^{2}{dr}+{k}^{2}\int \left[{\kappa }_{{\rm{e}}}^{2}+({\nu }^{2}+{\omega }_{{\rm{A}}}^{2})-\displaystyle \frac{4{{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}}{{\nu }^{2}+{\omega }_{{\rm{A}}}^{2}}\right]\,| \chi {| }^{2}{dr}=0,\end{eqnarray} \tag{ A11 }$$

where the integrand of the second integral should be negative, on average, to have solutions at all. This yields

$$\begin{eqnarray}&&{\nu }^{2}\leqslant {\nu }_{\max }^{2}=\tfrac{1}{2}\sqrt{\langle {\kappa }_{{\rm{e}}}^{4}\rangle +16\langle {{\rm{\Omega }}}^{2}{\omega }_{{\rm{A}}}^{2}\rangle }-\tfrac{1}{2}\langle {\kappa }_{{\rm{e}}}^{2}\rangle -\langle {\omega }_{{\rm{A}}}^{2}\rangle .\end{eqnarray} \tag{ A12 }$$

For the explicit example of the SARIs discussed in Section 4.1, this gives ${\nu }_{\max }=0.6307$ and 0.6277, in agreement with the actual growth rates of the uppermost modes of the Spectral Webs shown in Figures 5 and 8, which are much lower (0.3913 and 0.3246). This is understood since the first integral of Equation (A11) is not incorporated in expression (A12).

Distinguishing perturbations with a tilde from background equilibrium quantities, the current density perturbation may be written as

$$\begin{eqnarray}&&\tilde{{\boldsymbol{j}}}={\rm{\nabla }}\times \widetilde{{\boldsymbol{B}}},\qquad \widetilde{{\boldsymbol{B}}}={\rm{\nabla }}\times ({\boldsymbol{\xi }}\times {\boldsymbol{B}}).\end{eqnarray} \tag{ A13 }$$

Projecting $\tilde{{\boldsymbol{j}}}$ onto the radial, the perpendicular and parallel (to B ) directions yield

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{j}}_{r} & = & -{\rm{i}}(m/r)({B}_{z}\chi )^{\prime} +{\rm{i}}k\left(({B}_{\theta }/r)\chi \right)^{\prime} -{\rm{i}}(m/{r}^{2}+{k}^{2})B\eta ,\\ {\tilde{j}}_{\perp } & = & \displaystyle \frac{1}{{rB}}\left\{-{F}^{2}\chi +{{rB}}_{z}\left[(1/r)({B}_{z}\chi )^{\prime} +({mB}/r)\eta \right]^{\prime} +{B}_{\theta }\left[r\left(({B}_{\theta }/r)\chi \right)^{\prime} -{krB}\eta \right]^{\prime} \right\},\\ {\tilde{j}}_{\parallel } & = & \displaystyle \frac{1}{{rB}}\left\{{FG}\chi +{{rB}}_{\theta }\left[(1/r)({B}_{z}\chi )^{\prime} +({mB}/r)\eta \right]^{\prime} -{B}_{z}\left[r\left(({B}_{\theta }/r)\chi \right)^{\prime} -{krB}\eta \right]^{\prime} \right\},\end{array}\end{eqnarray} \tag{ A14 }$$

where η ≡ i(B_z ξ_θ − B_θ ξ_z )/B is the perpendicular component of the displacement vector ξ, and G ≡ mB_z /r − kB_θ and F ≡ mB_θ + kB_z are proportional to the perpendicular and parallel components of the "wavevector," respectively. One easily checks that the three current density components satisfy the constraint ${\rm{\nabla }}\cdot \tilde{{\boldsymbol{j}}}=(1/r)(r{\tilde{j}}_{r})^{\prime} +{\rm{i}}(G/B){\tilde{j}}_{\perp }+{\rm{i}}(F/B){\tilde{j}}_{\parallel }=0$.

For our purposes, the radial and parallel current density perturbations are not important since only ${\tilde{j}}_{\perp }$ produces a radially directed Lorentz force ${{\boldsymbol{e}}}_{r}\cdot (\tilde{{\boldsymbol{j}}}\times {\boldsymbol{B}})\equiv {\tilde{j}}_{\perp }B$. Note, from Equation (A14), that this expression involves the second derivative of χ, which is dominant for the rapidly (spatially) oscillating solutions involved, so that the contribution ${{\boldsymbol{e}}}_{r}\cdot ({\boldsymbol{j}}\times \tilde{{\boldsymbol{B}}})$ to the radial Lorentz force may be neglected. The expression for ${\tilde{j}}_{\perp }$ needs to be transformed to one in terms of the variables χ and Π that are available from the ODE solver used in the program ROC. The expression of the variable η in terms of χ and $\chi ^{\prime}$ may be found in Equation (12.93) of Goedbloed et al. (2019). It is convenient, though, to transform it again into one involving ${\rm{\Pi }}\equiv -(N/D)\chi ^{\prime} -(C/D)\chi$, exploiting Equation (20) for (N/D) and Equations (16) and (A3) for (C/D) in the incompressible approximation, giving

$$\begin{eqnarray}&&\eta \approx \displaystyle \frac{1}{{mB}}\left(-{B}_{z}\chi ^{\prime} +\displaystyle \frac{{krF}}{\rho {\widetilde{\omega }}^{2}-{F}^{2}}{\rm{\Pi }}\right).\end{eqnarray} \tag{ A15 }$$

This yields, after some algebra, the final expression for ${\tilde{j}}_{\perp }$:

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{j}}_{\perp } & = & \displaystyle \frac{1}{{rB}}\left\{a\chi ^{\prime\prime} +b\chi ^{\prime} +c\chi +\displaystyle \frac{d}{\widetilde{A}}{\rm{\Pi }}^{\prime} +\left(\displaystyle \frac{e}{\widetilde{A}}+\displaystyle \frac{f}{{\widetilde{A}}^{2}}+\displaystyle \frac{g\widetilde{\omega }}{{\widetilde{A}}^{2}}\right){\rm{\Pi }}\right\},\qquad \widetilde{A}\equiv \rho {\widetilde{\omega }}^{2}-{F}^{2},\\ a & \equiv & {{rB}}_{\theta }F/m,\qquad b\equiv {B}_{z}{B}_{z}^{\prime} +2{B}_{\theta }{B}_{\theta }^{\prime} -{B}_{\theta }^{2}/r+({{kB}}_{\theta }/m)({{rB}}_{z})^{\prime} ,\\ c & \equiv & {B}_{z}{B}_{z}^{\prime\prime} -{B}_{z}{B}_{z}^{\prime} /r+{B}_{\theta }{B}_{\theta }^{\prime\prime} -{B}_{\theta }{B}_{\theta }^{\prime} /r+{B}_{\theta }^{2}/{r}^{2}-{F}^{2},\qquad d\equiv {{kr}}^{2}{FG}/m,\\ e & \equiv & [{{kr}}^{2}G(F^{\prime} -\rho ^{\prime} F/\rho )-2{k}^{2}{{rB}}_{\theta }F]/m,\qquad f\equiv -{{kr}}^{2}{F}^{2}G(\rho ^{\prime} F/\rho -2F^{\prime} )/m,\qquad g\equiv 2{{kr}}^{2}{FG}\rho {\rm{\Omega }}^{\prime} .\end{array}\end{eqnarray} \tag{ A16 }$$

The real coefficients a, ..., g are explicitly known through the equilibrium solutions (4). The complex functions $\chi ^{\prime}$, χ, ${\rm{\Pi }}^{\prime}$, and Π are directly provided by the solution of ODEs (12), whereas χ'' is produced from $\chi ^{\prime}$ by fourth-order-accurate interpolation on a fine grid. Note that the Alfvén factor $\widetilde{A}$ and the Doppler-shifted frequency $\widetilde{\omega }$ are complex, so that the contributions of the real and imaginary components of ${\rm{\Pi }}^{\prime}$ and Π to expression (A16) for the perpendicular current density ${\tilde{j}}_{\perp }$ are mixed in a complicated way. Because of the rapid spatial oscillations of the "eigenfunctions," those contributions are rather small compared to the first term involving the second derivative χ''.