Quasi-periodic oscillations and the global modes of relativistic, MHD accretion discs

Dewberry, Janosz W; Latter, Henrik N; Ogilvie, Gordon I

doi:10.1093/mnras/sty385

Abstract

The high-frequency quasi-periodic oscillations that punctuate the light curves of X-ray binary systems present a window on to the intrinsic properties of stellar-mass black holes and hence a testbed for general relativity. One explanation for these features is that relativistic distortion of the accretion disc's differential rotation creates a trapping region in which inertial waves (r-modes) might grow to observable amplitudes. Local analyses, however, predict that large-scale magnetic fields push this trapping region to the inner disc edge, where conditions may be unfavourable for r-mode growth. We revisit this problem from a pseudo-Newtonian but fully global perspective, deriving linearized equations describing a relativistic, magnetized accretion flow, and calculating normal modes with and without vertical density stratification. In an unstratified model we confirm that vertical magnetic fields drive r-modes towards the inner edge, though the effect depends on the choice of vertical wavenumber. In a global model we better quantify this susceptibility, and its dependence on the disc's vertical structure and thickness. Our calculations suggest that in thin discs, r-modes may remain independent of the inner disc edge for vertical magnetic fields with plasma betas as low as β ≈ 100–300. We posit that the appearance of r-modes in observations may be more determined by a competition between excitation and damping mechanisms near the ISCO than by the modification of the trapping region by magnetic fields.

accretion, accretion discs, black hole physics, magnetic fields, MHD, waves, X-rays: binaries

1 INTRODUCTION

Understanding the variability observed in the emission from X-ray binaries remains an important task in astrophysics. So-called quasi-periodic oscillations (QPOs) frequently appear as 0.1–450 Hz features in the power density spectrum of stellar-mass black hole candidate systems. Historically, high-frequency quasi-periodic oscillations (HFQPOs) of ∼30–450 Hz, observed in anomalous states of high accretion and luminosity, have attracted interest for the comparability of their frequencies to the characteristic orbital frequencies close to a black hole. Because of their relative insensitivity to luminosity variations, HFQPOs are thought to issue from the black hole's imprint on the inner regions of its encircling accretion disc (Remillard & McClintock 2006; Motta 2016).

Okazaki, Kato & Fukue (1987) showed that in a purely hydrodynamic and isothermal model strong gravitational effects on the epicyclic frequencies κ and Ω_z create a narrow, annular, ‘self-trapping’ region where large-scale, global inertial waves (sometimes dubbed gravito-inertial modes or g-modes, but here called r-modes) might be constrained to oscillate as standing waves. Protected from dissipation at the inner disc edge by their confinement within this trapping region, hydrodynamic r-modes are subject to excitation by large-scale warping and eccentric deformations in the accretion flow, potentially growing to amplitudes sufficiently large to cause variations in luminosity detectable as HFQPOs (Kato 2004, 2008; Ferreira & Ogilvie 2008, 2009).

However, with temperatures of ≳1 keV providing sufficient ionization to support magnetic fields, discs around black holes are unlikely to be purely hydrodynamic. Magnetohydrodynamic (MHD) turbulence sustained by the magnetorotational instability (MRI) offers a widely accepted explanation for an effective viscosity driving accretion (Balbus & Hawley 1998). The survival of trapped inertial modes in the presence of such MHD turbulence is uncertain; in fact, Reynolds & Miller (2008) observed prominent r-mode signatures in hydrodynamic simulations, but, along with Arras, Blaes & Turner (2006), found them absent from MRI-turbulent simulations. It should be noted, though, that Arras et al. (2006) and Reynolds & Miller (2008) omitted any r-mode excitation mechanism other than the MHD turbulence itself, the latter authors concluding only that the MRI does not actively excite the oscillations. On the other hand, taking an analytical approach, Fu & Lai (2009) questioned the very existence of the self-trapping region when the accretion disc is threaded by large-scale poloidal magnetic fields of appreciable strength. However, Fu & Lai (2009) presented only a local analysis, omitting the radial structure of the inherently global oscillations, as well as vertical density stratification.

In this paper, we present a fully global, semi-analytical treatment of the problem, calculating r-mode oscillations within isothermal relativistic accretion discs threaded by azimuthal and vertical magnetic fields. We focus on two models: ‘cylindrical’ discs, which omit vertical density structure, and fully global discs, which include it. The wave modes are numerically calculated using pseudo-spectral and hybrid pseudo-spectral-Galerkin methods, respectively.

In support of the conclusions of Fu & Lai (2009), we find that increasingly strong vertical fields force the localisation of r-modes inwards, towards the innermost stable circular orbit (ISCO). In cylindrical models, azimuthal fields of equipartiton strength have little effect. Meanwhile, r-mode sensitivity to the vertical magnetic field depends on the parameterization of vertical structure through isothermal sound speed c_s and vertical wavenumber k_z, and is hence somewhat uncertain. In fully global models, the point at which r-mode confinement depends on reflection at the inner boundary depends on disc temperature and thickness. For isothermal sound speeds of c_s ∼ 0.001c, where c is the speed of light, we find critical values of the plasma beta (ratio of gas to magnetic pressure) as low as β ∼ 100–300, depending on the scale height's rate of increase with radius. These estimates are consistent with those of Fu & Lai (2009), but we stress that such field strengths are relatively high for large-scale, ordered magnetic fields; simulations of the MRI produce fields of this strength or greater, but they are small-scale and disordered. Finally, we demonstrate that with a choice of vertical wavenumber motivated by the basis functions used in our density stratified analysis, r-modes calculated in a cylindrical model reproduce the frequencies and localizations of fully global trapped inertial modes to within |$1\,\,\rm{per\,\,cent}$|⁠.

Our results indicate that if the accretion disc is sufficiently thin, the trapping of r-modes is minimally affected by large-scale magnetic fields (unless those fields are very strong). And indeed, QPOs only appear in emission states usually described by thin disc models. The frequencies, however, of the r-modes will be enhanced by the presence of large-scale fields, which may complicate their use when divining the properties of the central black hole. On the other hand, this frequency enhancement could be used to estimate the strength of the magnetic field itself, especially if the black hole mass and spin are constrained by other measurements. Finally, we stress that even if r-modes are pushed up against the ISCO this does not necessarily mean that they are destroyed. The oscillations will certainly be damped by radial inflow but could nonetheless achieve appreciable amplitudes if sufficiently excited by a strong disc eccentricity and/or warp. In short, very strong large-scale magnetic fields may not exterminate r-modes, just make them more difficult to excite.

The structure of the paper is as follows. In Section 2, we review the nature of global oscillation modes and the potential for their formation in relativistic accretion discs. Readers familiar with the subject are invited to skip to Section 3, where we present the linearized equations describing our magnetized, relativistic accretion disc model. In Section 4, we examine the effects of azimuthal and vertical magnetic fields on radially global modes calculated in the cylindrical approximation, and in Section 5 we present our vertically stratified results. Finally, in Section 6 we summarize and discuss our findings.

2 BACKGROUND

Analogous to helioseismology, discoseismology involves the study of accretion disc oscillations that are global in that they maintain their frequency and a coherent structure across a considerable radial extent. Examples might include warps and eccentricities, which can be thought of as non-axisymmetric waves with very low frequency (Ferreira & Ogilvie 2009), or large-scale spiral density waves. Subject to sufficient excitation, such oscillations might reach amplitudes large enough to cause observable variations in luminosity.

The establishment and persistence of global oscillations in a realistic accretion flow may require specific conditions, however, especially at the disc boundaries. For example, the growth of inertial-acoustic oscillations in relativistic discs via the co-rotation instability requires wave reflection at the inner radius and the transmission of wave energy at key resonances (Papaloizou & Pringle 1984; Narayan, Goldreich & Goodman 1987; Lai et al. 2012), while the r-modes considered in this paper may be damped by radial inflow at the inner disc edge (Ferreira 2010).

Accretion discs around stellar-mass black holes offer protection from such damping, the effects of strong gravity possibly providing a narrow region separated from the inner edge within which global r-mode oscillations could reside and grow. While the horizontal epicyclic frequency, κ, has a monotonic radial profile in Newtonian centrifugally supported discs in near Keplerian rotation, strong gravity introduces radii at which κ² attains a maximum and falls below zero. The latter radius defines the ISCO within which matter plunges towards the black hole, while the existence of the former has spurred the development of several theories of wave propagation aimed at explaining HFQPOs (Kato 2001). In this section, we review the key elements of these theories, treating hydrodynamic and magnetohydrodynamic models in turn.

2.1 Hydrodynamic waves in relativistic discs

Okazaki et al. (1987) used a pseudo-Newtonian treatment with a Paczynski–Wiita potential (see Section 3) in considering a hydrodynamic model of a thin, vertically isothermal, relativistic accretion disc. The authors argued that although a rigorous treatment of an accretion disc around a black hole would require a general relativistic model, the effects of strong gravity on wave propagation arise primarily through the relativistic modification of the epicyclic frequency. They showed that the linearized equations derived for perturbations of the form δ(r, z)exp [imϕ − iωt], for azimuthal wavenumber m and frequency ω, are approximately separable in r and z under the assumption of a scale height H(r) = c_s/Ω_z that varies slowly with radius. This separation allows for a description of the perturbations’ vertical structure by modified Hermite polynomials of vertical order n. Projecting on to a given vertical order leaves a decoupled, second-order system of ordinary differential equations in r which, with appropriate boundary conditions, can be solved numerically for the direct calculation of global normal modes.

Although we are interested in global oscillations, local analyses provide insight into the effect that the epicyclic frequency's non-monotonic radial profile has on wave propagation. Neglecting the radial variation of background quantities and assuming a radial dependence for perturbations of δ(r) ∝ exp [ik_rr], where k_r is a local radial wavenumber, Okazaki et al. (1987) derived the dispersion relation

\begin{eqnarray} k_r^2 =\frac{\left ({\tilde{\omega }}^2-\kappa ^2\right ) \left ({\tilde{\omega }}^2-n\Omega _z^2\right )}{\tilde{\omega }^2 c_s^2}, \end{eqnarray}

(1)

for |$\tilde{\omega }=\omega -m\Omega$|⁠. Since oscillatory behaviour requires real k_r, the radial profile for |$-k_r^2$| can be thought of as an effective ‘potential well’, allowing wave propagation wherever |$-k_r^2<0$| (Li, Goodman & Narayan 2003). The profiles for the horizontal and vertical epicyclic frequencies κ and Ω_z then define different regions of localization in the disc for three different families of waves. When n = 0, 2D ‘inertial-acoustic’ waves can propagate at radii where |$\tilde{\omega }^2>\kappa ^2$|⁠. Non-zero n gives wave propagation where either |$\tilde{\omega }^2>\max [\kappa ^2,n\Omega _z^2]=n\Omega _z^2$| or |$\tilde{\omega }^2<\min [\kappa ^2,n\Omega _z^2]=\kappa ^2$|⁠. The former, higher frequency waves are known as acoustic or p-modes, while the latter, lower frequency waves are the inertial r-modes of our interest. Fig. 1 shows the regions of propagation for such global modes in the axisymmetric case, as defined by the conditions on ω² (top) and |$-k_r^2$| (bottom).

$Axisymmetric (m = 0) wave propagation regions for fully relativistic characteristic frequencies κG and ΩGz (defined in Section 3) with spin parameter $a=0$, sound speed $c_s=0.01c$ and vertical quantum number $n=1$. R-mode trapping regions, as defined by bounding frequencies κ2 and $\Omega _z^2$ (top), and representative effective potential wells $-k_r^2$ (bottom) are plotted both with no magnetic field (solid lines) and with a vertical magnetic field with midplane β ∼ 500 (bottom, dash–dotted line).$

Figure 1.

Axisymmetric (m = 0) wave propagation regions for fully relativistic characteristic frequencies κ_G and Ω_Gz (defined in Section 3) with spin parameter |$a=0$|⁠, sound speed |$c_s=0.01c$| and vertical quantum number |$n=1$|⁠. R-mode trapping regions, as defined by bounding frequencies κ² and |$\Omega _z^2$| (top), and representative effective potential wells |$-k_r^2$| (bottom) are plotted both with no magnetic field (solid lines) and with a vertical magnetic field with midplane β ∼ 500 (bottom, dash–dotted line).

Open in new tab Download slide

The Lindblad, co-rotation and vertical resonances occur where |$\tilde{\omega }^2=\kappa ^2$|⁠, |$\tilde{\omega }=0$|⁠, and |$\tilde{\omega }^2=n\Omega _z^2$|⁠, respectively. They define the wave propagation regions in a hydrodynamic, isothermal disc, and depend closely on the radial profiles of the orbital and epicyclic frequencies. As illustrated in Fig. 1(top), the maximum in κ² (solid line) introduced by the prescription of relativistic frequencies provides a narrow region of confinement for modes with non-zero n. The two Lindblad resonances at radii where |$\tilde{\omega }^2=\kappa ^2$| define turning points, within which the inertial modes might be confined irrespective of conditions at the ISCO. The effective potential well defined by |$-k_r^2$| (solid line in Fig. 1, bottom) has a minimum at the radius where κ achieves its maximum, denoted as R_κ. This minimum is only local, however, as the vertical resonance occurring at the radius where |$\tilde{\omega }^2=n\Omega _z^2$| allows r-modes to ‘leak’ to the outer regions of the disc and propagate as p-modes. The extent of this leakage is limited by the effective potential barrier between the two propagation regions (Ferreira & Ogilvie 2008).

In isolation, trapped inertial waves are almost purely oscillatory, a small exponential decay coming only from wave leakage. Any explanation of HFQPOs involving r-modes then requires an excitation mechanism.¹ Kato (2004, 2008) estimated the growth rates of trapped inertial waves coupled to large-scale warps or eccentricities in the accretion flow, describing a global analogue to the local parametric instability introduced by Goodman (1993) and Ryu & Goodman (1994). This coupling involves the transfer of negative wave energy from the fundamental, axisymmetric r-mode, through the global disc deformations, to non-axisymmetric inertial modes propagating within their co-rotation radius. Ferreira & Ogilvie (2008) generalized the estimations of Kato (2004, 2008), calculating explicitly the growth rates of the fundamental trapped inertial mode when excited by warps and eccentricities, and proffering r-modes as a promising explanation for HFQPOs.

2.2 Magnetohydrodynamic waves in relativistic discs

The global calculations conducted by Okazaki et al. (1987) and Ferreira & Ogilvie (2008) bore out the predictions of the local analysis, bolstering the argument that trapped waves cause HFQPOs. However, as recognized by Reynolds & Miller (2008) and Fu & Lai (2009), the picture is much more complicated for magnetized accretion flows.

Reynolds & Miller (2008) searched for trapped inertial waves in global simulations of relativistic accretion discs utilizing a Paczynski–Wiita potential. Peaks in the power density spectrum at the frequencies and radii appropriate for r-modes appeared in initially perturbed, hydrodynamic accretion discs, but were absent from the spectrum obtained from MHD-turbulent simulations. The authors concluded that turbulence caused by the MRI does not actively excite trapped inertial modes, but made no claims about active damping. Indeed, O'Neill, Reynolds & Miller (2009) saw trapped r-mode propagation in viscous simulations with α > 0.05, but found that the wave amplitudes were below the noise levels observed in MHD simulations by Reynolds & Miller (2008). Arras et al. (2006) observed a similar absence of inertial waves in MHD turbulent shearing box simulations. However, like Reynolds & Miller (2008), these authors did not include any form of excitation for r-modes other than the MRI turbulence itself. Henisey et al. (2009) considered such an excitation mechanism in simulations of MRI turbulent, relativistic discs with non-zero spin, including a tilt to examine numerically the excitation mechanism explored by Ferreira & Ogilvie (2008), and observed the emergence of modes ‘at least partially inertial in character.’

Separately, Fu & Lai (2009) put aside the question of MHD turbulence and argued that large-scale magnetic fields destroy the r-mode trapping region itself. They derived a local dispersion relation analogous to equation (1) for a relativistic accretion disc threaded by a purely vertical, constant magnetic field B_z, and found the self-trapping region modified by an Alfvénic restoring force. Rather than being constrained to propagate only where ω² < κ², their local analyses predict that axisymmetric, MHD r-modes are evanescent except in the region where

\begin{eqnarray} k_z^2V_{\rm{A}z}^2 <\omega ^2< \frac{1}{2} \left[ \kappa ^2+2k_z^2V_{\rm{A}z}^2 +\sqrt{\kappa ^4+16k_z^2V_{\rm{A}z}^2\Omega ^2} \right], \end{eqnarray}

(2)

where k_z is the vertical wavenumber corresponding to an exponential dependence exp [ik_zz], and |$V_{\rm{A}z}=B_z/\sqrt{\mu _0\rho }$| is the Alfvén speed.² The inclusion of a large-scale vertical magnetic field modifies the effective potential well defined by |$-k_r^2$| (see Section 4.2.1), an illustrative radial profile of which is plotted with a dash–dotted line in Fig. 1(bottom) for β ≈ 500 and the prescription of a vertical wavenumber k_z = 1/H(r) = Ω_z(r)/c_s. In the presence of a magnetic field, the potential well moves closer to the ISCO, as indicated. For much stronger fields (smaller β), the trapping region includes the inner edge of the disc, implying that r-mode oscillations would need to be confined by the outer turning point and the inner edge itself, where conditions may be uncertain and most likely unfavourable (e.g. Gammie 1999; Afshordi & Paczynski 2003; Ferreira 2010).

The self-trapping region's loss of distinction from the inner boundary does not necessarily rule out the existence of r-modes. Modes pushed to the inner edge of the disc would need to rely on excitation sufficient to overcome the damping effects of radial inflow (Ferreira 2010), and might require some reflection mechanism for their confinement. Such a reflection is required for explanations of HFQPOs that involve 2D, non-axisymmetric inertial-acoustic modes, which local analyses suggest are less susceptible to the effects of large-scale magnetic fields and subject to amplification at co-rotation (Lai & Tsang 2009; Fu & Lai 2011; Lai et al. 2012; Yu & Lai 2015).

In summary, Reynolds & Miller (2008) and Fu & Lai (2009) threw the efficacy of trapped inertial waves as an explanation for HFQPOs into some doubt. However, the simulations only indicated that the MRI does not actively excite r-modes, and in applying local approximations Fu & Lai (2009) neglected the global nature of discoseismic oscillations. Local analyses do give insight into the structure of wave propagation in accretion discs, but the WKBJ ansatz limits us to scales much smaller than those of the background flow, an inappropriate assumption when considering oscillations coherent across large radial and vertical extents. Additionally, although the scaling k_z ∝ 1/H is natural in the nearly separable, hydrodynamic problem (Okazaki et al. 1987), the inseparability of the MHD equations in r and z makes the exact proportionality unclear. Further, as noted by Kato (2017), the inverse dependence of the Alfvén speed on density enhances the importance of vertical density stratification. These considerations motivate our focus on fully global normal mode calculations, and our examination of the role that density stratification plays in localizing MHD trapped inertial waves.

3 DISC MODEL AND LINEARIZED EQUATIONS

We consider a thin, inviscid, non-self-gravitating, differentially rotating flow. The ideal MHD equations can be written as

\begin{eqnarray} \frac{\mathrm{\partial} {{\boldsymbol u}}}{\mathrm{\partial} t} +{{\boldsymbol u}\cdot \nabla {\boldsymbol u}} =-\frac{\nabla P}{\rho } -\nabla \Phi +\frac{1}{\mu _0\rho }{\left(\nabla \times {\boldsymbol B}\right)\times {\boldsymbol B}}, \end{eqnarray}

(3)

\begin{eqnarray} \frac{\mathrm{\partial} \rho }{\mathrm{\partial} t} =-\nabla \cdot {(\rho {\boldsymbol u})}, \end{eqnarray}

(4)

\begin{eqnarray} \frac{\mathrm{\partial} {{\boldsymbol B}}}{\mathrm{\partial} t} =\nabla \times \left({{\boldsymbol u\times B}}\right), \end{eqnarray}

(5)

\begin{eqnarray} \nabla \cdot {{\boldsymbol B}}=0, \end{eqnarray}

(6)

where |${{\boldsymbol u}}$|⁠, ρ, |${\boldsymbol B}$|⁠, and P are the fluid velocity, density, magnetic field, and gas pressure, respectively.

Equations (3)–(6) provide a Newtonian description of a plasma in the limit of ideal MHD, and are therefore not strictly applicable to relativistic accretion discs. However, as mentioned in Section 2, relativistic effects outside of the ISCO are frequently approximated, in both analysis and simulation (e.g. Okazaki et al. 1987; Reynolds & Miller 2008), with the prescription of a ‘pseudo-Newtonian’ Paczynski–Wiita potential, given in cylindrical coordinates (r, ϕ, z) as

\begin{eqnarray} \Phi =\frac{-GM}{\sqrt{r^2+z^2}-r_S}, \end{eqnarray}

(7)

where G is the gravitational constant, M the mass of the central black hole, and r_S = 2GM/c² is the Schwarzschild radius of the event horizon. In the absence of modification by pressure gradients and magnetic stresses, the introduction of a singularity at r = r_S gives midplane (z = 0) orbital and horizontal epicyclic frequencies

\begin{eqnarray} \Omega _P^2 =\frac{GM}{r(r-r_S)^2}, \end{eqnarray}

(8)

\begin{eqnarray} \kappa _P^2 =\frac{GM(r-3r_S)}{r(r-r_S)^3}, \end{eqnarray}

(9)

the second of which passes through zero at r = 3r_S, and achieves a maximum at r ≈ 3.73r_S. The point at which κ_P = 0 defines the ISCO, denoted by r_ISCO. The maximum replicates qualitatively that which appears in the exact frequencies for a particle in orbit around a Kerr black hole, given by

\begin{eqnarray} \Omega _G =\frac{1}{(r^{3/2}+a)}, \end{eqnarray}

(10)

\begin{eqnarray} \kappa _G =\Omega _G\sqrt{1-\frac{6}{r}+\frac{8a}{r^{3/2}}-\frac{3a^2}{r^2}}, \end{eqnarray}

(11)

\begin{eqnarray} \Omega _{Gz} =\Omega _G\sqrt{1-\frac{4a}{r^{3/2}}+\frac{3a^2}{r^2}}, \end{eqnarray}

(12)

where a ∈ ( − 1, 1) is the dimensionless spin parameter, and radii and frequencies are expressed in units of gravitational radius r_g = GM/c² and the gravitational frequency ω_g = c³/(GM). These expressions are commonly inserted in otherwise Newtonian analyses to approximate the effects of black hole spin on wave propagation in the accretion flow (Ferreira & Ogilvie 2008).

We consider a background equilibrium flow of the form |${{\boldsymbol u}}=r\Omega (r){\hat{\boldsymbol {\phi }}}$|⁠, in isorotation with a magnetic field |${{\boldsymbol B}}=B_{\phi }(r){\hat{\boldsymbol {\phi }}}+B_z(r){\hat{\boldsymbol {z}}}$|⁠. We neglect radial magnetic fields for their complication of the equilibrium. For a gravitational potential of the form given in equation (7), the radial component of equation (3) then yields

\begin{eqnarray} \Omega ^2 =\frac{1}{r}\bigg [ \frac{\mathrm{\partial} \Phi }{\mathrm{\partial} r} +\frac{1}{\rho }\frac{\mathrm{\partial} P}{\mathrm{\partial} r} +\frac{1}{\mu _0\rho }\bigg ( B_{\phi }\frac{\rm{d}B_{\phi }}{\rm{d}r} +B_z\frac{\rm{d}B_z}{\rm{d}r} +\frac{B_{\phi }^2}{r} \bigg )\bigg ]. \end{eqnarray}

(13)

For simplicity we assume the gas is globally isothermal with |$P=c_s^2\rho$|⁠. In a thin disc without magnetic fields, a barotropic equation of state with constant c_s guarantees an equilibrium angular velocity profile Ω that depends only on radius. As can be seen from equation (13), this remains true when our setup includes either a density profile ρ = ρ(r) independent of z (as considered in Section 4), or a purely uniform, vertical magnetic field (as considered in Section 5). In reality, the temperature, and hence the sound speed, ought to decrease with radius. We discuss the ramifications of a globally isothermal equation of state in Section 6; in short, we predict this assumption leads to an overestimate of the disc scale height as a function of radius, and consequently an overestimate of r-modes’ susceptibility to background magnetic fields.

Equation (13) implies that steep gradients in the pressure and magnetic fields, as well as the so-called hoop stress due to azimuthal magnetic fields, may cause deviations from the orbits of particles rotating with ‘pseudo-Keplerian’ angular velocities defined by |$\Omega _K^2=(1/r)\mathrm{\partial} _r\Phi$|⁠. We include these deviations in calculations utilizing a Paczynski–Wiita potential, but ignore them as |$\mathcal {O}(1/r^2)$| effects when utilizing the general relativistic formula for the characteristic frequencies. The vertical independence of |${{\boldsymbol u}}$| and |${{\boldsymbol B}}$| yields, in the classic thin disc approximation, the hydrostatic equilibrium |$c_s^2\mathrm{\partial} _z \ln \rho =-\mathrm{\partial} _z \Phi \approx -\Omega _z^2 z$|⁠, which gives a background density distribution ρ(r, z) = ρ₀(r)exp [ − z²/(2H²)].

The independence of the equilibrium from ϕ and t allows us to consider perturbations of the form δ(r, ϕ, z, t) = δ(r, z)exp [imϕ − iωt]. Because non-axisymmetric r-modes encounter a co-rotation resonance and may be strongly damped at the radius where |$\tilde{\omega }=\omega -m\Omega =0$| (Li et al. 2003), historical focus has been given to axisymmetric trapped inertial waves with m = 0. The co-rotation resonance might be avoided by non-axisymmetric modes with frequencies large enough that the radius where ω = mΩ(r) lies outside of the self-trapping region. However, such frequencies would be of the order of kHz, too high to explain HFQPO observations in stellar-mass black hole binaries. For this reason, we restrict our attentions to axisymmetric modes. In particular, we focus on the fundamental axisymmetric r-mode with the simplest radial and vertical structure, as it is the least likely to be disrupted, and the most likely to produce a net luminosity variation.

We disturb the equilibrium with axisymmetric fluid velocity, magnetic field and enthalpy (h ≡ δP/ρ) perturbations of the form |$[{{\boldsymbol v}},{{\boldsymbol b}},h](r,z)\exp [-\rm{i}\omega t]$|⁠, where ω is in general complex. Linearizing equations (3)–(6) then yields

\begin{eqnarray} -\rm{i}\omega v_r &=& 2\Omega v_{\phi } -\left[ \frac{\mathrm{\partial} }{\mathrm{\partial} r} -\frac{1}{\mu _0P} \left(B_z\frac{\rm{d}B_z}{\rm{d}r}+\frac{B_{\phi }}{r}\frac{\rm{d}(rB_{\phi })}{\rm{d}r}\right) \right]h\nonumber \\ &+&\frac{1}{\mu _0\rho } \Bigg [ B_z\frac{\mathrm{\partial} b_r}{\mathrm{\partial} z}-\left(B_{\phi } \frac{\mathrm{\partial} }{\mathrm{\partial} r} +\frac{\rm{d}B_{\phi }}{\rm{d}r}+\frac{2B_{\phi }}{r} \right)b_{\phi }\nonumber\\ &-&\left( B_z\frac{\mathrm{\partial} }{\mathrm{\partial} r}+\frac{\rm{d}B_z}{\rm{d}r} \right)b_z \Bigg ], \end{eqnarray}

(14)

\begin{eqnarray} -\rm{i}\omega v_{\phi } &=&-\frac{\kappa ^2}{2\Omega }v_r +\frac{1}{\mu _0\rho } \left( \frac{1}{r}\frac{\rm{d}(rB_{\phi })}{\rm{d} r}b_r +B_z\frac{\mathrm{\partial} b_{\phi }}{\mathrm{\partial} z} \right), \end{eqnarray}

(15)

\begin{eqnarray} -\rm{i}\omega v_z &=&-\frac{\mathrm{\partial} h}{\mathrm{\partial} z} +\frac{1}{\mu _0\rho } \left( \frac{\rm{d}B_z}{\rm{d} r}b_r -B_{\phi }\frac{\mathrm{\partial} b_{\phi }}{\mathrm{\partial} z} \right), \end{eqnarray}

(16)

\begin{eqnarray} -\rm{i}\omega h &=&-c_s^2 \left\lbrace \left[ \frac{\mathrm{\partial} }{\mathrm{\partial} r}+\frac{1}{r} \left( 1+\frac{\rm{d} \ln \rho _0}{\rm{d} \ln r} \right) \right] v_r\right.\nonumber \\ &&\left. +\left( \frac{\mathrm{\partial} }{\mathrm{\partial} z} +\frac{\mathrm{\partial} \ln \rho }{\mathrm{\partial} z} \right)v_z \right\rbrace , \end{eqnarray}

(17)

\begin{eqnarray} -\rm{i}\omega b_r &=&B_z\frac{\mathrm{\partial} v_r}{\mathrm{\partial} z}, \end{eqnarray}

(18)

\begin{eqnarray} \!-\!\rm{i}\omega b_{\phi }&=&\!-\!\left( B_{\phi }\!\frac{\mathrm{\partial} }{\mathrm{\partial} r} \!+\!\frac{\rm{d}B_{\phi }}{\rm{d}r} \right)v_r +B_z\frac{\mathrm{\partial} v_{\phi }}{\mathrm{\partial} z} \!-\!B_{\phi }\!\frac{\mathrm{\partial} v_z}{\mathrm{\partial} z} \!+\!\frac{\rm{d}\Omega }{\rm{d}\ln r}b_r, \end{eqnarray}

(19)

\begin{eqnarray} -\rm{i}\omega b_z&=&-\left( B_z\frac{\mathrm{\partial} }{\mathrm{\partial} r} +\frac{\rm{d}B_z}{\rm{d}r} +\frac{B_z}{r} \right)v_r, \end{eqnarray}

(20)

where |$\nabla \cdot {{\boldsymbol b}}=0$| has been substituted into the induction equation, so that the perturbations satisfy the solenoidal condition by construction. In the following sections, we solve equations (14)–(20) under various approximations and in full.

4 CYLINDRICAL CALCULATIONS

To set the scene and isolate clearly the radially global aspects of MHD r-modes, we first present calculations made in the cylindrical approximation, in which the vertical lengthscale of the perturbations is assumed to be much smaller than that of the equilibrium flow. Strictly speaking, this assumption is inappropriate; the r-modes of most interest are those with the simplest vertical structure. (These are the least affected by large-scale magnetic fields.) However, the cylindrical model is particularly attractive from a numerical standpoint, providing a less expensive framework for global, non-linear simulations. Furthermore, we show in Sections 5.3.3 and 5.4.2 that with a specific choice of vertical wavenumber, cylindrical r-modes can be very closely identified with trapped inertial waves calculated from a fully global and self-consistent model.

Neglecting the vertical dependence of Φ and ρ in equations (14)–(20), we Fourier transform and prescribe a plane-wave z-dependence exp [ik_zz], where the vertical wavenumber k_z is assumed for simplicity and self-consistency to be constant (revisited in Section 5.4.2). We do consider radial variation in |${{\boldsymbol B}}$| and ρ in making calculations with a Paczynski–Wiita potential, writing

\begin{eqnarray} {{\boldsymbol B}} &=&B_{\phi 0}\left(\frac{r}{r_0}\right)^q{\hat{\boldsymbol {\phi }}}+B_{z0}\left(\frac{r}{r_0}\right)^p{\hat{\boldsymbol {z}}}, \end{eqnarray}

(21)

\begin{eqnarray} \rho &=& \rho_{00} \left(\frac{r}{r_0}\right)^{\sigma }, \end{eqnarray}

(22)

where r₀ is the radius of the ISCO in the absence of modification by gas pressure and magnetic forces, B_ϕ0 and B_z0 are constant magnetic field strengths, |$\rho_{00}$| is a constant density, and q, p, and σ are constant power-law indices. For convenience, we follow Fu & Lai (2009) in defining the midplane Alfvénic Mach number |$\mathcal {M}_{\rm{A}}=|{{\boldsymbol V_{\rm{A}}}}|/c_s=\sqrt{2}\beta ^{-1/2}$|⁠, where |${{\boldsymbol V_{\rm A}}}={{\boldsymbol B}}/\sqrt{\mu _0\rho }$| is the midplane Alfvén velocity and β is the midplane plasma beta. We then write |$\mathcal {M}_{\rm{A}z}(r)=V_{\rm{A}z}/c_s$| and |$\mathcal {M}_{\rm{A}\phi }(r)=V_{\rm{A}\phi }/c_s$|⁠, noting that both are functions of radius. With these prescriptions, the equilibrium flow subject to a Paczynski–Wiita potential follows the midplane angular velocity

\begin{eqnarray} \Omega ^2=\frac{GM}{r(r-r_S)^2} +\frac{c_s^2}{r^2}\big [\sigma +\left (q+1\right )\mathcal {M}_{\rm{A}\phi }^2+p\mathcal {M}_{\rm{A}z}^2\big ]. \end{eqnarray}

(23)

Note that we only use this expression when p, q, σ ≠ 0. When omitting all background radial gradients in |${\boldsymbol B}$| and ρ we employ the correct general relativistic frequencies, setting Ω ≈ Ω_G and κ ≈ κ_G, from equations (10)–(12).

Trading |${{\boldsymbol b}}$| for the Alfvén velocity perturbation |${{\boldsymbol v}}_{\rm A}={{\boldsymbol b}}/\sqrt{\mu _0\rho (r)}$|⁠, equations (14)–(20) can be written as

\begin{eqnarray} -\rm{i}\omega v_r &=&2\Omega v_{\phi } -\left( \frac{\rm{d} }{\rm{d}r} -\frac{(q+1)}{r}\mathcal {M}_{\rm{A}\phi }^2-\frac{p}{r}\mathcal {M}_{\rm{A}z}^2 \right)h \nonumber \\ &&+\, c_s \Bigg [ \rm{i}k_z\mathcal {M}_{\rm{A}z}v_{\rm{A}r} -\mathcal {M}_{\rm{A}\phi } \left( \frac{\rm{d} }{\rm{d} r}+\frac{\sigma +2q+4}{2r} \right)v_{\rm{A}\phi } \nonumber \\ &&-\, \mathcal {M}_{\rm{A}z} \left( \frac{\rm{d} }{\rm{d} r}+\frac{\sigma +2p}{2r} \right)v_{\rm{A}z} \Bigg ], \end{eqnarray}

(24)

\begin{eqnarray} -\rm{i}\omega v_{\phi } &=&-\frac{\kappa ^2}{2\Omega }v_r +c_s \left( \frac{(q+1)}{r}\mathcal {M}_{\rm{A}\phi }v_{\rm{A}r} +\rm{i}k_z\mathcal {M}_{\rm{A}z}v_{\rm{A}\phi } \right), \nonumber\\ \end{eqnarray}

(25)

\begin{eqnarray} -\rm{i}\omega v_z &=&-\rm{i}k_zh +c_s \big( \frac{p}{r}\mathcal {M}_{\rm{A}z}v_{\rm{A}r} -\rm{i}k_z\mathcal {M}_{\rm{A}\phi }v_{\rm{A}\phi } \big), \end{eqnarray}

(26)

\begin{eqnarray} -\rm{i}\omega h &=&-c_s^2 \left[ \left( \frac{\rm{d}}{\rm{d} r}+\frac{1+\sigma }{r} \right) v_r +\rm{i}k_zv_z \right], \end{eqnarray}

(27)

\begin{eqnarray} -\rm{i}\omega v_{\rm{A}r} &=&\rm{i}k_zc_s \mathcal {M}_{\rm{A}z}v_r, \end{eqnarray}

(28)

\begin{eqnarray} -\rm{i}\omega v_{\rm{A}\phi }\!&=&\! -c_s \left[ \mathcal {M}_{\rm{A}\phi } \left( \frac{\rm{d}}{\rm{d}r}+\frac{q}{r} \right)v_r -\rm{i}k_z\mathcal {M}_{\rm{A}z}v_{\phi } +\rm{i}k_z\mathcal {M}_{\rm{A}\phi }v_z \right] \nonumber \\ &&+\, \frac{\rm{d}\Omega }{\rm{d}\ln r}v_{\rm{A}r}, \end{eqnarray}

(29)

\begin{eqnarray} -\rm{i}\omega v_{\rm{A}z} &=&-c_s\mathcal {M}_{\rm{A}z} \left( \frac{\rm{d}}{\rm{d} r}+\frac{p+1}{r} \right)v_r. \end{eqnarray}

(30)

Non-dimensionalizing radii, frequencies, and velocities in units of r_g, ω_g, and r_gω_g = c (respectively), equations (24)–(30) can be written as a generalized eigenvalue problem of the form |${{\boldsymbol{\sf A}}} {U}=\omega {{{\boldsymbol{\sf B}} U}}$|⁠, where |${{\boldsymbol{\sf B}}}$| specifies the boundary conditions placed on the eigenvector |${{\boldsymbol U}}$|⁠, and solved using a pseudo-spectral method. With derivatives approximated by Chebyshev spectral derivative matrices, this requires only the use of standard numerical packages for finding the eigenvalues and eigenvectors of matrices (Boyd 2001). All normal mode calculations made in the cylindrical approximation were performed using a Gauss–Lobatto grid on r ∈ [r_ISCO, r_ISCO + 14r_g] with a discretization of N = 200.

4.1 Cylindrical, hydrodynamic r-modes

We point out that the hydrodynamic dispersion relation in the cylindrical approximation is similar to equation (1), but with |$n\Omega _z^2$| replaced by |$k_z^2c_s^2.$| In our simplified model with constant k_z, the vertical resonance has vanished: |$k_z^2c_s^2$| does not decrease with radius like |$n\Omega _z^2$|⁠, so p-modes with |$\omega ^2>k_z^2c_s^2>\kappa ^2$| can propagate freely, while r-modes with ω² < κ² are evanescent everywhere except in the narrow region near the radius of maximal κ, R_κ. As a result, the r-modes calculated with the cylindrical model experience no decay from wave leakage, and ω is entirely real. This justifies the use of the rigid wall boundary condition v_r = 0 at the outer as well as the inner boundaries in hydrodynamic calculations using the cylindrical approximation.

Apart from the absence of wave leakage, our hydrodynamic trapped inertial modes resemble those of Ferreira & Ogilvie (2008). The solid curves in Fig. 2 show example radial profiles of Re[v_r] for normal modes. Here, the calculations employ the correct general relativistic formulas for the frequencies Ω_G and κ_G. Though hydrodynamic calculations are not the focus of this paper, we briefly summarize the effects of different parameters as a point of comparison with previous work. Increasing k_z causes the fundamental trapped inertial mode's confinement to narrow, and its frequency to grow. Increasing c_s, which causes the potential well |$-k_r^2$| to become wider, expands the radial extent of the r-modes. Using a Paczynski–Wiita potential to investigate radial density variation, we find that in the hydrodynamic model, physically reasonable power-law indices 0 > σ ≥ −1.5 give a pressure gradient that acts only minimally through Ω and κ.

$Radial profiles of Re[vr] for hydrodynamic and magnetohydrodynamic r-modes calculated with non-zero Bz, with inner and outer boundary conditions given as vr = 0 (top), dvϕ/dr = 0 (middle), and dbr/dr = 0 (bottom). General relativistic characteristic frequencies were used (⟹p, q, σ = 0), with a = 0 and $k_z=\Omega _z(r_0)/c_s\approx 22.7r_g^{-1}$ for cs = 0.003c. The filled black and unfilled red dots give the radii of maximal vr (denoted rmax) and inner radii where $-k_r^2=0$ for the calculated frequencies ω (denoted rTR).$

Figure 2.

Radial profiles of Re[v_r] for hydrodynamic and magnetohydrodynamic r-modes calculated with non-zero B_z, with inner and outer boundary conditions given as v_r = 0 (top), dv_ϕ/dr = 0 (middle), and db_r/dr = 0 (bottom). General relativistic characteristic frequencies were used (⟹p, q, σ = 0), with a = 0 and |$k_z=\Omega _z(r_0)/c_s\approx 22.7r_g^{-1}$| for c_s = 0.003c. The filled black and unfilled red dots give the radii of maximal v_r (denoted r_max) and inner radii where |$-k_r^2=0$| for the calculated frequencies ω (denoted r_TR).

Open in new tab Download slide

4.2 Cylindrical, magnetohydrodynamic r-modes

4.2.1 Magnetohydrodynamic boundary conditions

Although nonzero magnetic fields introduce more derivatives to equations (24)–(30), with appropriate choices of variable the system can be reduced to a single ODE of second order (e.g. Appert, Gruber & Vaclavik 1974; Ogilvie & Pringle 1996). We then require exactly two boundary conditions, and considering a local dispersion relation for a radial wavenumber k_r (derived from equations 24–30 with the assumptions |$\delta(r)\propto\exp[{\rm i}k_rr],\,\, k_r,k_z\gg1/r$|⁠) provides insight into which may be appropriate (for the case with |$B_{\phi} = 0$|⁠, see equation 28 in Fu & Lai 2009).³

For frequencies within the range defined by equation (2) and constant |$k_z\gtrsim\Omega_z(r_0)/c_s$|⁠, we find that effective potential wells defined by |$-k_r^2$| remain positive at |$r_{\rm{out}} = r_{\rm{ISCO}}+14r_g$| for the magnetic field strengths considered. This is because our assumption of a constant k_z still excludes a vertical resonance. The choice of outer boundary condition then makes little difference. Sufficiently strong vertical magnetic fields can drive the profiles for |$-k_r^2$| to negative values at the ISCO, however, suggesting that MHD r-mode confinement may become dependent on reflection at the inner boundary. To assess the degree of this dependence, we also consider the boundary conditions |$db_r/dr = 0$| and |$dv_{\phi}/dr = 0$| at |$r_{\rm{in}} = r_{\rm{ISCO}}$|⁠, utilizing the former unless otherwise stated.

4.2.2 Purely vertical magnetic field, |$B_{\phi }(r)=0,B_z(r)\not=0$|

We use the same pseudo-spectral method as earlier to calculate solutions with a non-zero vertical magnetic field. As in the hydrodynamical case, each calculation at a given k_z produces a discrete set of modes, each with radial structures of increasing complexity. The modes may be distinguished by quantum numbers l, where l = 0 corresponds to the fundamental, which exhibits the simplest structure.

Fig. 2 shows plots of Re[v_r] for fundamental MHD r-modes calculated with a purely constant vertical magnetic field of strengths up to |$\mathcal {M}_{\rm{A}z}=0.15$| (β ≲ 100), using general relativistic formulas for the characteristic frequencies, all three boundary conditions, c_s = 0.003c and the naive assumption |$k_z=\Omega _z(r_0)/c_s\approx 22.7r_g^{-1}$|⁠. The frequencies Re[ω] can be used to construct profiles for |$-k_r^2$|⁠, which give a local estimate, here called r_TR, of the r-mode turning point, or inner edge of the trapping region (plotted as red, unfilled dots in Fig. 2). The black dots give the modes’ radius of maximal radial velocity, r_max. While r_max does not measure the evanescence of the r-modes, unlike r_TR it provides a global measure of their localization, and qualitatively illustrates the effect of the vertical magnetic field.

Fig. 2 indicates that the inclusion of a purely constant vertical magnetic field forces r-modes towards the plunging region within the ISCO, as predicted by the local analysis. This migration is accompanied by changes in mode frequency, illustrated in Fig. 3. These alterations to the mode behaviour occur because the addition of a vertical magnetic field amplifies the restoring force to horizontal disturbances, causing the r-modes to become hybrid epicyclic-Alfvénic, and their frequencies to increase with B_z as well as k_z. If the Alfvénic restoring force is too large, the modes cease to be confined without reflection at the inner boundary (compare the middle and bottom panels of Fig. 2 with the top, although a non-zero value of v_r at the inner boundary does not necessarily mean the mode is not evanescent there). We find that the frequency increases shown in Fig. 3 follow the dependence on Alfvén frequency ω_Az = k_zV_Az given by equation (2), but diverge by as much as |$10\,\,\rm{per\,\,cent}$| for larger values of k_z and V_Az.

$Re[ω] for the fundamental l = 0 r-mode calculated in the cylindrical approximation, as in Fig. 2(bottom), for varying kz and Alfvénic Mach number $\mathcal {M}_{\rm{A}z}$. Modes found to have rmax ≤ rISCO are excluded.$

Figure 3.

Re[ω] for the fundamental l = 0 r-mode calculated in the cylindrical approximation, as in Fig. 2(bottom), for varying k_z and Alfvénic Mach number |$\mathcal {M}_{\rm{A}z}$|⁠. Modes found to have r_max ≤ r_ISCO are excluded.

Open in new tab Download slide

The degree to which the mode's localization is affected by |$\mathcal {M}_{\rm{A}z}$| depends on both the vertical wavenumber k_z and the sound speed c_s, the combination of which parametrizes vertical structure in the cylindrical model. Fig. 4 shows a heatmap of r_max for the modes whose frequencies are shown in Fig. 3. The two figures show a similar dependence of the mode frequency and localization on |$\mathcal {M}_{{\rm A}z}$| and k_z. The sound speed chosen has a more significant impact on r_max than on ω, however. Naively choosing k_z = Ω_z(r₀)/c_s and varying c_s and |$\mathcal {M}_{\rm{A}z}$|⁠, we find that lower values of c_s make the r-modes marginally less susceptible to the vertical magnetic field, as measured by r_max. More importantly, larger values of c_s increase the width of the trapping regions for the r-modes, causing them to encounter the inner boundary at lower magnetic field strengths. Including a non-zero spin parameter also prolongs the r-modes’ march towards the ISCO.

$Radius of maximal vr, rmax , for the fundamental l = 0 r-mode calculated in the cylindrical approximation, as in Fig. 2(bottom), for varying kz and Alfvénic Mach number $\mathcal {M}_{\rm{A}z}$. Modes found to have rmax ≤ rISCO are excluded.$

Figure 4.

Radius of maximal v_r, r_max, for the fundamental l = 0 r-mode calculated in the cylindrical approximation, as in Fig. 2(bottom), for varying k_z and Alfvénic Mach number |$\mathcal {M}_{\rm{A}z}$|⁠. Modes found to have r_max ≤ r_ISCO are excluded.

Open in new tab Download slide

Performing calculations with a Paczynski–Wiita potential yields qualitatively the same results, and allows for the convenient evaluation of the effects of radial density and magnetic field variation, which modify both the rotation profile (23) and equations (24)–(30). Radial density variation with power laws 0 ≥ σ ≥ −1.5 impact MHD r-modes more than in the hydrodynamic case, acting through the Alfvén speed to delay the point at which r_max ∼ r_ISCO, and to marginally increase their frequencies further. Radially decreasing B_z acts also through the background Alfvén velocity V_Az, but instead forces the fundamental r-mode closer to the ISCO at a given value of B_z0. This can be understood through equation (2). As the r-modes migrate inwards, the restoring force from the magnetic field becomes larger and larger.

4.2.3 Purely toroidal magnetic field, |$B_{\phi }(r)\not=0,B_z(r)=0$|

Trapped inertial mode calculations made with |$B_{\phi }\not=0,B_z=0$| and rigid wall boundary conditions confirm the predictions of Fu & Lai (2009) that large-scale azimuthal magnetic fields modify r-modes only negligibly. For the ranges of sound speed and vertical wavenumber considered in Section 4.2.2, an equipartition field strength with β = 1 (⁠|$\mathcal {M}_{{\rm A}\phi }=\sqrt{2}$|⁠) is required for even a |$0.5\,\,\rm{per\,\,cent}$| change in frequency. Changes in r-mode localization are imperceptible.

5 VERTICALLY STRATIFIED, FULLY GLOBAL CALCULATIONS

Our calculations of trapped inertial waves made in the ideal MHD, cylindrical approximation confirm the qualitative behaviour predicted by the local analyses of Fu & Lai (2009). Although our simple cylindrical model provides an inexpensive framework for non-linear simulations, it nevertheless has too many free parameters, and is strictly inappropriate for the case of most interest with k_z ∼ 1/H. To find accurate measurements of the critical field strengths at which the mode confinement becomes dependent on the inner boundary, then, it is essential to consider the vertical structure of the disc. In this section, we present a fully general, vertically stratified model for axisymmetric oscillations in an accretion disc threaded by a uniform vertical magnetic field in Sections 5.1 and 5.2, and calculations of critical magnetic field strengths both without and with the coupling of vertical modes provided by radial scale height variation in Sections 5.3 and 5.4 (respectively).

In this section, we also seek to validate the use of an unstratified, cylindrical model to study the behaviour of trapped inertial waves and other global oscillations. The ansatz of an exponential dependence |$\delta (r,z)\propto \tilde{\delta }(r)\exp [\rm{i}k_zz]$| relies on the assumption that the background flow's scale of vertical variation is much larger than that of the perturbations. Cylindrical models may then be understood to target phenomena taking place at the midplane of the disc. For this approach to be valid, we must find a correspondence between r-modes calculated with both the cylindrical model, and a fully general treatment. In particular, we must identify an appropriate vertical wavenumber. The choice of k_z greatly influences cylindrical r-modes’ response to increasing vertical magnetic field strength (see Section 4), but, as discussed in Section 2, is made unclear by the non-separability introduced by such a field. In Sections 5.3.3 and 5.4.2, we present calculations demonstrating a clear correspondence both without and with (respectively) the inclusion of radial scale height variation.

The importance of vertical density stratification has been recognized by Kato (2017), who noted that the Alfvén speed associated with a purely constant vertical magnetic field would have a much steeper vertical than horizontal gradient in an isothermal disc. Kato (2017) also considered a vertically stratified, isothermal accretion disc threaded by a purely constant vertical magnetic field, but used radially local analyses, and asymptotic expansions with |$\mathcal {M}_{\rm{A}z}=V_{\rm{A}z}/c_s$| as a small parameter. Further, the author used the argument of a rarified corona to impose a rigid lid. Here, we instead present fully global calculations made without explicit restrictions on the magnetic field strength, and relax this imposition of vertical disc truncation.

5.1 Stratified equations: a change of variables

Since our cylindrical analyses suggest that both azimuthal magnetic fields and radial variation in ρ and |${\boldsymbol B}$| have very little effect on the trapped inertial modes, we consider a purely constant |${{\boldsymbol B}}=B_z{\hat{\boldsymbol {z}}}$| and set p = q = σ = 0. Thus, the density profile is radially constant: ρ = ρ₀₀g(η), where η(r, z) ≡ z/H(r) is a new vertical coordinate, and g = exp [ − η²/2] is the Gaussian density profile associated with an isothermal disc.

In addition to this change of coordinate, we modify equations (14)–(20) with a change of dependent variable. In the case of a purely constant |${{\boldsymbol B}}=B_z{\hat{\boldsymbol {z}}},$| the perturbation to the Lorentz force |${{\boldsymbol J}}\times {\boldsymbol B}$| is proportional to |$(\nabla \times {{\boldsymbol b}})\times {{\boldsymbol B}}$|⁠, prompting the definition

\begin{eqnarray} {{\boldsymbol L}} \equiv B_z\left[ \frac{1}{H}\frac{\mathrm{\partial} b_r}{\mathrm{\partial} \eta } -\left( \frac{\mathrm{\partial} }{\mathrm{\partial} r} +\frac{\mathrm{\partial} \eta }{\mathrm{\partial} r}\frac{\mathrm{\partial} }{\mathrm{\partial} \eta } \right) b_z \right]{\hat{\boldsymbol {r}}}+\frac{B_z}{H}\frac{\mathrm{\partial} b_{\phi }}{\mathrm{\partial} \eta }{\hat{\boldsymbol {\phi }}}. \end{eqnarray}

(31)

Eliminating b_ϕ and b_z in favour of |${{\boldsymbol L}}$|’s r and ϕ components, L_r and L_ϕ, and trading the enthalpy perturbation |$h=c_s^2\delta \rho /\rho$| for the non-dimensional Γ = δρ/ρ for convenience, equations (14)–(20) can be re-written as

\begin{eqnarray} -\rm{i}\omega v_r &=&2\Omega v_{\phi } -c_s^2\left( \frac{\mathrm{\partial} }{\mathrm{\partial} r} -\frac{\rm{d}\ln H}{\rm{d} r}\eta \frac{\mathrm{\partial} }{\mathrm{\partial} \eta } \right)\Gamma +\frac{L_r}{\mu _0\rho }, \end{eqnarray}

(32)

\begin{eqnarray} -\rm{i}\omega v_{\phi } &=&-\frac{\kappa ^2}{2\Omega }v_r +\frac{L_{\phi }}{\mu _0\rho }, \end{eqnarray}

(33)

\begin{eqnarray} -\rm{i}\omega v_z &=&-\frac{c_s^2}{H}\frac{\mathrm{\partial} \Gamma }{\mathrm{\partial} \eta }, \end{eqnarray}

(34)

\begin{eqnarray} -\rm{i}\omega \Gamma &=&-\mathcal {L}_rv_r +\frac{\rm{d} \ln H }{\rm{d} r}\frac{\eta }{g}\frac{\mathrm{\partial} (gv_r)}{\mathrm{\partial} \eta } -\frac{1}{gH}\frac{\mathrm{\partial} (gv_z)}{\mathrm{\partial} \eta }, \end{eqnarray}

(35)

\begin{eqnarray} -\rm{i}\omega b_r &=&\frac{B_z}{H}\frac{\mathrm{\partial} v_r}{\mathrm{\partial} \eta }, \end{eqnarray}

(36)

\begin{eqnarray} \rm{i}\omega L_r = B_z^2\left\lbrace \mathcal {L}_{rr} \!-\!\mathcal {L}_{HH}\eta \!\frac{\mathrm{\partial} }{\mathrm{\partial} \eta } \!+\!\left[ \frac{1}{H^2} \!+\!\left(\frac{\rm{d}\ln H}{\rm{d}r}\right)^2\eta ^2 \right]\!\frac{\mathrm{\partial} ^2}{\mathrm{\partial} \eta ^2} \right\rbrace v_r, \nonumber\\ \end{eqnarray}

(37)

\begin{eqnarray} -\rm{i}\omega L_{\phi } =\frac{B_z^2}{H^2}\frac{\mathrm{\partial} ^2 v_{\phi }}{\mathrm{\partial} \eta ^2} +\frac{B_z}{H}\frac{\rm{d}\Omega }{\rm{d}\ln r}\frac{\mathrm{\partial} b_r}{\mathrm{\partial} \eta }, \end{eqnarray}

(38)

where we have defined the differential operators

\begin{eqnarray} \mathcal {L}_r \equiv \frac{\mathrm{\partial} }{\mathrm{\partial} r}+\frac{1}{r}, \end{eqnarray}

(39)

\begin{eqnarray} \mathcal {L}_{rr} \equiv \frac{\mathrm{\partial} ^2 }{\mathrm{\partial} r^2} +\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r} -\frac{1}{r^2}, \end{eqnarray}

(40)

\begin{eqnarray} \mathcal {L}_{HH} \equiv \frac{\rm{d}\ln H}{\rm{d}r} \left(2\frac{\mathrm{\partial} }{\mathrm{\partial} r} +\frac{1}{r} -2\frac{\rm{d}\ln H}{\rm{d}r} \right) +\frac{1}{H}\frac{\rm{d}^2H}{\rm{d}r^2}. \end{eqnarray}

(41)

5.2 Boundary conditions and series expansions

In this section, we describe the series expansions and numerical method utilized in solving equations (32)–(38). Readers interested only in the results of our calculations may skip to Sections 5.3 and 5.4. For numerical efficiency we follow Okazaki et al. (1987) and Ogilvie (2008) in making the ansatz that for a given perturbation δ(r, η), we can write

\begin{eqnarray} \delta (r, \eta ) =\sum _n^{\infty }\delta _n(r)\mathcal {B}_n(\eta ), \end{eqnarray}

(42)

where |$\lbrace \mathcal {B}_n\rbrace$| is a set of orthogonal basis functions satisfying the boundary conditions imposed on δ(r, η) as η → ±∞. Finding the exact solution would require solving for an infinite number of radially variable coefficients δ_n(r), but we find that with the appropriate choice of |$\mathcal {B}_n$|⁠, truncating the series expansions at a finite n = M can produce excellent approximations to the eigenvalues ω.

In the hydrodynamic case, Okazaki et al. (1987) showed that the fluid variables |${{\boldsymbol v}}$| and h are well described by modified Hermite polynomials He_n of order n. However, in the MHD problem, the r and ϕ components of the induction equation imply that velocities increasing polynomially with η would result in an unbounded bending of the field lines high above the disc. Although bending may in reality occur for large-scale poloidal field loops with footpoints in the outer disc, our focus on the inner radii of the black hole accretion disc prompts us to impose the condition that the magnetic field lines become purely vertical at large η. This implies that |$b_r,b_{\phi }, \mathrm{\partial} _zv_r, \mathrm{\partial} _zv_{\phi }\rightarrow 0$| as |η| → ∞. A more appropriate choice, then, is to expand v_r, v_ϕ and Γ as

\begin{eqnarray} v_r(r, \eta ) =\sum _{n=1}^{\infty }u_n(r)F_n(\eta ), \end{eqnarray}

(43)

\begin{eqnarray} v_{\phi }(r, \eta ) =\sum _{n=1}^{\infty }v_n(r)F_n(\eta ), \end{eqnarray}

(44)

\begin{eqnarray} \Gamma (r, \eta ) =\sum _{n=1}^{\infty }\Gamma _n(r)F_n(\eta ), \end{eqnarray}

(45)

where u_n, v_n, and Γ_n are functions to be determined. The basis functions F_n(η) are solutions to the equation

\begin{eqnarray} \frac{\rm{d}^2F}{\rm{d}\eta ^2}+K^2gF=0, \end{eqnarray}

(46)

where K is a constant, dimensionless eigenvalue. Subject to the boundary condition that |$\mathrm{\partial} _{\eta }F\rightarrow 0$| as η → ∞ this equation (derived by Gammie & Balbus 1994; Ogilvie 1998; Latter, Fromang & Gressel 2010) is in Sturm–Liouville form, with guaranteed discrete sets {F_n} and {K_n} (written as {K_nH} in Latter et al. 2010). The F_n describe the horizontal velocity components of MRI channel modes in the local, anelastic approximation, and are related to another set of basis functions, {G_n(η)}, by

\begin{eqnarray} \frac{\rm{d}F_n}{\rm{d}\eta }=K_nG_n. \end{eqnarray}

(47)

Both sets of functions (cf. Fig. A1 and fig. 1 in Latter et al. 2010) are orthogonal, and can be normalized such that

\begin{eqnarray} \int _{-\infty }^{\infty }gF_nF_m\, \mathrm{d}\eta &=\delta _{nm}, \end{eqnarray}

(48)

\begin{eqnarray} \int _{-\infty }^{\infty }G_nG_m\, \mathrm{d}\eta &=\delta _{nm}. \end{eqnarray}

(49)

The functions {G_n} provide an appropriate basis set for the radial magnetic field perturbation, which we expand as

\begin{eqnarray} b_r(r, \eta )=\sum _{n=1}^{\infty }b_r^n(r)G_n(\eta ). \end{eqnarray}

(50)

For the Lorentz force perturbation, we appeal to the models of Gammie & Balbus (1994) and Sano & Miyama (1999), who assumed that far from the disc the magnetic field goes to a force free configuration, such that |${{\boldsymbol L}}\rightarrow 0$|⁠. Finally, for the vertical velocity perturbation, we note that the assumption that Γ goes to a constant implies that v_z → 0 as η → ±∞ as well. In theory, then, v_z, L_r, and L_ϕ can be expanded in any set of orthogonal basis functions that go to zero at infinity, but we find the fastest convergence for

\begin{eqnarray} v_z(r, \eta )&=&\sum _{n=1}^{\infty }w_n(r)g^{1/4}(\eta )F_{n-1}(\eta ), \end{eqnarray}

(51)

\begin{eqnarray} L_{r}(r, \eta )&=&\sum _{n=1}^{\infty }L_{r}^n(r)g(\eta )F_n(\eta ), \end{eqnarray}

(52)

\begin{eqnarray} L_{\phi }(r, \eta )&=&\sum _{n=1}^{\infty }L_{\phi }^n(r)g(\eta )F_n(\eta ), \end{eqnarray}

(53)

where w_n, |$L_r^n$|⁠, and |$L_\phi ^n$| are functions to be determined.

It should be noted that G₀ = 0, K₀ = 0, and F₀ is a constant. Since F₋₁ is undefined, the expansion for v_z must start from n = 1. For the vertically structured inertial waves, the n = 0 basis functions do not describe any of the other perturbations, and their exclusion has no impact on the solutions.

The basis functions for v_z and |${{\boldsymbol L}}$| are also orthogonal with the weights g^1/2 and g⁻¹ (respectively). Substituting our expansions into equations (32)–(38), we gain an infinite number of coupled systems of seven equations, but this orthogonality allows us to greatly reduce the dimensionality of our problem. First scaling r and H by r_g, frequencies by ω_g, |${{\boldsymbol v}}$| and c_s by r_gω_g = c, b_r by B_z, and |${{\boldsymbol L}}$| by |$B_z^2/r_g$|⁠, we multiply equations (32)–(38) by suitable combinations of basis functions of arbitrary order m and integrate over |$\int _{-\infty }^{\infty }\, \mathrm{d}\eta$| to gain (after swapping indices)

\begin{eqnarray} \!-\!\rm{i}\omega u_n &=&2\Omega v_n \!-\!c_s^2\frac{\rm{d}\Gamma _n }{\rm{d} r} \!+\!c_s^2\!\frac{\rm{d}\ln H}{\rm{d} r}\sum _{m=1}^{\infty }\nu _{nm}\Gamma _m \!+\!\mathcal {M}_{\rm{A}z}^2c_s^2L_r^n, \end{eqnarray}

(54)

\begin{eqnarray} -\rm{i}\omega v_n &=&-\frac{\kappa ^2}{2\Omega }u_n +\mathcal {M}_{\rm{A}z}^2c_s^2L_{\phi }^n, \end{eqnarray}

(55)

\begin{eqnarray} -\rm{i}\omega w_n &=&-\frac{c_s^2}{H}\sum _{m=1}^{\infty }\theta _{nm}\Gamma _m, \end{eqnarray}

(56)

\begin{eqnarray} -\rm{i}\omega \Gamma _n &=&-\mathcal {L}_ru_n -\frac{\rm{d}\ln H}{\rm{d} r}\sum _{m=1}^{\infty }\lambda _{nm}u_m +\frac{1}{H}\sum _{m=1}^{\infty }\alpha _{nm}w_m, \end{eqnarray}

(57)

\begin{eqnarray} -\rm{i}\omega b_r^n &=&\frac{K_n}{H}u_n, \end{eqnarray}

(58)

\begin{eqnarray} &&{-\rm{i}\omega \sum _{m=1}^{\infty }\frac{\gamma _{nm}}{K_m^2}L_r^m =\mathcal {L}_{rr}u_n }\nonumber \\ &&\quad-\sum _{m=1}^{\infty }\left[ \nu _{nm}\mathcal {L}_{HH} +\frac{1}{H^2}\gamma _{nm} +\left(\frac{\rm{d}\ln H}{\rm{d}r}\right)^2\mu _{nm} \right]u_m, \end{eqnarray}

(59)

\begin{eqnarray} -\rm{i}\omega L_{\phi }^n =-\frac{K_n^2}{H^2}v_n -\frac{K_n}{H}\frac{\rm{d}\Omega }{\rm{d}\ln r}b_r^n, \end{eqnarray}

(60)

where the coefficients θ_nm, α_nm, γ_nm, λ_nm, ν_nm, and μ_nm are defined by finite integrals involving various combinations of η, g, F_n, F_m, G_n, and G_m (see Appendix A1). These coupling coefficients connect the equation sets of different orders, and thus illustrate the non-separability of the problem, though as we shall see the coupling is relatively weak. A given set of equations of nth order couples to other sets of equations with orders n + 2, n + 4, … of the same parity, but the integrals quickly go to zero away from n ≈ m. This suggests that we may lose little in accuracy by truncating the series in equations (54)–(60) at a relatively small m = M. A danger with this numerical method is that inappropriately chosen basis functions can lead to unconverged solutions. In general, however, we find that with these basis expansions, calculating a converged trapped inertial mode of a given vertical order n only requires extending the series expansion to M ≥ n + 4. With truncation at larger M, energetic contributions to the solution from coefficients of higher order are negligible, as are changes in frequency (further discussion in Appendix A2).

Discretizing our radial domain on a Gauss–Lobatto grid with N gridpoints, equations (54)–(60) can once again be written as a generalized eigenvalue problem |${{\boldsymbol{\sf A}} {\boldsymbol U}}=\omega {{\boldsymbol{\sf B}} {\boldsymbol U}}$|⁠, where |${\boldsymbol{\sf B}}$| now carries the coupling coefficients γ_nm in addition to encoding the boundary conditions, of which we require two for each order n included in the truncated series. The solution takes the form |${{\boldsymbol U}} =({{\boldsymbol U}}_1,{{\boldsymbol U}}_2,\ldots ,{{\boldsymbol U}}_M)^T$| with nth order radial coefficients |${{\boldsymbol U}}_n =({{\boldsymbol u}}_n,{{\boldsymbol v}}_n,{{\boldsymbol w}}_n,{{\boldsymbol h}}_n,{{\boldsymbol b}}_r^n,{{\boldsymbol L}}_r^n,{{\boldsymbol L}}_{\phi }^n)^T$|⁠. We find that 10 grid points per r_g are sufficient to resolve the least radially complicated fundamental r-mode of physical interest.

Once the truncated series of coupled eigenvalue problems has been solved, the fully global, (r, η)-dependent fields can be reconstructed for each perturbation variable δ(r, η) with the summations |$\sum _n^M\delta _n(r)\mathcal {B}_n(\eta )$|⁠. In the following sections, we present fully global r-mode solutions calculated without (Section 5.3) and with (Section 5.4) the coupling provided by radial scale height variation. We use our calculations to find approximations to the critical magnetic field strengths at which r-modes require reflection at the inner boundary for confinement, as well as a correspondence with, and k_z prescription for, the cylindrical model presented in Section 4.

5.3 Global r-modes in a flat disc

In this section, we make the simplifying approximation that the scale height of the disc is purely constant: H = c_s/Ω_z(r_ISCO). This is motivated by the relative proximity of the trapping region to the inner edge of the disc and H’s weak radial variation within this region. With this approximation |$\mathrm{\partial} _r\ln H=0$| and |$\mathcal {L}_{HH}=0$|⁠, so the couplings involving the coefficients λ_nm, ν_nm, and μ_nm drop out. The absence of a vertical resonance again makes the outer boundary condition irrelevant, since the r-modes should be evanescent at r_out as long as it is placed beyond the outer turning point.

5.3.1 The global r-mode spectrum

Solving equations (54)–(60) with series truncation at a vertical order n = M produces M sets of global trapped inertial modes. Each set consists of modes with similar vertical structure but differing radial structure. As with our cylindrical calculations, we order the members of each set according to the number of radial nodes in v_r, e.g. l = 0, 1, 2, …. The members of each set are dominated by a single component |${{\boldsymbol U}}_k$| of order n = k, testifying to the weak separability of the problem. As a consequence, all the members of that set express a vertical profile closely resembling |${{\mathcal B}}_k$|⁠, and for that reason we assign a quantum number (or vertical order) ‘k’ to the set of modes dominated by the kth order radial coefficients. There are hence two numbers denoting modes: l and k, representing their radial and vertical quantization, respectively.

For example, Fig. 5 shows a representative fundamental r-mode. The (r, η) heat maps show the full solution |${{\boldsymbol U}}(r, \eta )$|⁠, and alongside them we plot the radial profiles of the various coefficients |${{\boldsymbol U}}_n(r)$| that constitute the full solution. As is clear from the latter, the mode is completely dominated by the n = 1 component. We hence assign the vertical quantum number k = 1 to this mode, and its radial structure tells us that l = 0. In addition, Fig. 6 shows the (exaggerated) effect of the mode on the background magnetic field lines, close to the trapping region.

$Left: Real parts of the radially variable coefficients $[u_n,v_n,w_n, \Gamma _n,b_r^n,L_r^n,L_{\phi }^n](r)$ for the radially and vertically fundamental (l = 0, k = 1) trapped inertial mode calculated with cs = 0.003c, a = 0, r/rg ∈ [6, 20], N = 200, M = 7, and midplane $\mathcal {M}_{\rm{A}z}=0.04$ (β ≈ 1250). Right: Heatmaps of the fully reconstructed solutions Re[vr, vϕ, vz, Γ, br, Lr, Lϕ](r, η). The induction equation has been used to find the fields for bϕ(r, η) (bottom left) and bz(r, η) (bottom right). Note that the magnetic field and Lorentz force perturbations are proportionately larger in magnitude because they are scaled by much smaller quantities than the velocity perturbations, which are scaled by the speed of light.$

Figure 5.

Left: Real parts of the radially variable coefficients |$[u_n,v_n,w_n, \Gamma _n,b_r^n,L_r^n,L_{\phi }^n](r)$| for the radially and vertically fundamental (l = 0, k = 1) trapped inertial mode calculated with c_s = 0.003c, a = 0, r/r_g ∈ [6, 20], N = 200, M = 7, and midplane |$\mathcal {M}_{\rm{A}z}=0.04$| (β ≈ 1250). Right: Heatmaps of the fully reconstructed solutions Re[v_r, v_ϕ, v_z, Γ, b_r, L_r, L_ϕ](r, η). The induction equation has been used to find the fields for b_ϕ(r, η) (bottom left) and b_z(r, η) (bottom right). Note that the magnetic field and Lorentz force perturbations are proportionately larger in magnitude because they are scaled by much smaller quantities than the velocity perturbations, which are scaled by the speed of light.

Open in new tab Download slide

$Magnetic field geometry in the r − z plane, close to the trapping region for the l = 0, k = 1 r-mode shown in Fig. 5. Magnetic field lines are plotted as evenly spaced contours of the poloidal magnetic flux function. Given the axisymmetry of the field, this can be calculated as $\Psi (r,z)=\int r(\mathcal {A}b_z+B_z)\, \mathrm{d}r$, where the amplifying factor $\mathcal {A}=1/\max [b_z]$ normalizes the linear perturbation (which is already scaled by Bz) so as to make its effect discernible.$

Figure 6.

Magnetic field geometry in the r − z plane, close to the trapping region for the l = 0, k = 1 r-mode shown in Fig. 5. Magnetic field lines are plotted as evenly spaced contours of the poloidal magnetic flux function. Given the axisymmetry of the field, this can be calculated as |$\Psi (r,z)=\int r(\mathcal {A}b_z+B_z)\, \mathrm{d}r$|⁠, where the amplifying factor |$\mathcal {A}=1/\max [b_z]$| normalizes the linear perturbation (which is already scaled by B_z) so as to make its effect discernible.

Open in new tab Download slide

To illustrate the structure of the oscillations within the disc, the component of the linearized momentum perturbation p_r = ρ₀₀gv_r is shown for similarly calculated modes with l = 0, 1, 2 radial and k = 1, 2, 3 vertical structures in Fig. 7. In the weakly magnetized case with small |$\mathcal {M}_{\rm{A}z},$| the r-modes possessing the same radial node structure but dominated by different n have nearly the same frequencies. With increasing magnetic field strength, however, these groups of solutions separate in both the frequency and physical domains. As foreshadowed by the dependence on k_z observed in our cylindrical calculations, we find that the more vertically complicated modes with larger quantum number k are more strongly affected by the vertical magnetic field. This can be seen in the radial momentum fields shown in Fig. 7. The k = 2 modes (middle row) have been forced further inwards by the background magnetic field (β = 1250) than the k = 1 modes (top row), and the k = 3 modes (bottom row) even further. Additionally, with more complicated vertical structure, increases in Re[ω] are larger for modes of a given radial structure.

Real parts of the radial momentum perturbation pr = ρ00gvr for trapped inertial waves with radial quantum numbers l = 0, 1, 2 (from left to right) and vertical orders k = 1, 2, 3 (from top to bottom), calculated as in Fig. 5. The black dashed line indicates Rκ, while the white dashed lines mark z = ±H and z = ±2H.

Figure 7.

Real parts of the radial momentum perturbation p_r = ρ₀₀gv_r for trapped inertial waves with radial quantum numbers l = 0, 1, 2 (from left to right) and vertical orders k = 1, 2, 3 (from top to bottom), calculated as in Fig. 5. The black dashed line indicates R_κ, while the white dashed lines mark z = ±H and z = ±2H.

Open in new tab Download slide

As mentioned, the radially and vertically fundamental mode (Fig. 5, top left in Fig. 7) is the least likely to be disrupted and most likely to cause an observable luminosity variation, so we focus on its response to the vertical magnetic field.

5.3.2 Critical magnetic field strengths

We seek a convenient metric for evaluating the critical magnetic field strength at which trapped inertial mode confinement becomes dependent on conditions at the inner boundary. At low magnetic field strengths and sound speeds the frequencies measured (Fig. 8) are extremely insensitive to the boundary conditions imposed at r_in and r_out, agreeing to the precision allowed by our numerical method. However, for sufficiently strong magnetic fields, the frequencies do become dependent on the inner boundary, with different boundary conditions producing different frequencies (Fig. 9). We identify the point at which this happens as one (conservative) measure of the critical magnetic field strength for r-mode independence from the ISCO. As a point of comparison, we have also considered the first derivatives of the radial velocity for modes calculated with rigid boundary conditions. The magnetic field strength at which these derivatives become non-zero at the inner boundary provides a very similar, slightly more stringent measure.

$Frequencies for the fundamental l = 0, k = 1 r-mode, calculated using a purely constant scale height (and the boundary condition $\mathrm{\partial} _rb_r=0$), with varying cs and midplane $\mathcal {M}_{\rm{A}z}$, for a = 0 (circles), a = 0.5 (triangles), and a = 0.8 (stars). Sound speed is indicated with the same colour-map used in Fig. 9, black corresponding to cs = 0.001c and yellow to cs = 0.01c (r ∈ [rISCO, rISCO + 14rg], N = 300, M = 7).$

Figure 8.

Frequencies for the fundamental l = 0, k = 1 r-mode, calculated using a purely constant scale height (and the boundary condition |$\mathrm{\partial} _rb_r=0$|⁠), with varying c_s and midplane |$\mathcal {M}_{\rm{A}z}$|⁠, for a = 0 (circles), a = 0.5 (triangles), and a = 0.8 (stars). Sound speed is indicated with the same colour-map used in Fig. 9, black corresponding to c_s = 0.001c and yellow to c_s = 0.01c (r ∈ [r_ISCO, r_ISCO + 14r_g], N = 300, M = 7).

Open in new tab Download slide

$Heatmap of percent differences between frequencies calculated as in Fig. 8, using the boundary conditions vr = 0 (ωV) and $\mathrm{\partial} _rb_r=0,$ (ωB), with varying sound speeds and magnetic field strengths (a = 0).$

Figure 9.

Heatmap of percent differences between frequencies calculated as in Fig. 8, using the boundary conditions v_r = 0 (ω_V) and |$\mathrm{\partial} _rb_r=0,$| (ω_B), with varying sound speeds and magnetic field strengths (a = 0).

Open in new tab Download slide

Fig. 10 shows estimates of the critical magnetic field strengths attained with the former metric. We plot the field strengths at which the percentage difference between frequencies calculated with different boundary conditions reaches |$0.01\,\,\rm{per\,\,cent}$| (i.e. when 100 × (ω_V − ω_B)/ω_B reaches 10⁻², where ω_V and ω_B are the frequencies calculated with the boundary conditions v_r = 0 and |$\mathrm{\partial} _rb_r=0,$| respectively). Fig. 10 illustrates the inverse relationship between disc temperature and r-mode resilience to vertical magnetic fields. This is to be expected, since trapped inertial modes are less well-confined in hotter, thicker discs even without magnetic fields (Ferreira & Ogilvie 2008). For sound speeds ranging from c_s = 0.001c to c_s = 0.01c, we find critical midplane Alfvénic Mach numbers of |$\mathcal {M}_{\rm{A}z}\sim 0.14{\rm -}0.06$| (β ∼ 100–550). We stress, however, that these values arise from the rather stringent condition that the frequencies ω_V and ω_B differ by one-hundredth of a percent, a somewhat arbitrary value. If we allow frequency discrepancies as large as |$0.1\,\,\rm{per\,\,cent}$| the critical Mach numbers rises to values above 0.1, and β values lower than 100.

$Critical Alfvénic Mach number $\mathcal {M}_{\rm{Ac}}$ found from the constant scale height mode calculations presented in Fig. 8 by considering the percent variations in frequency between modes calculated with the inner boundary conditions vr = 0 and $\mathrm{\partial} _rb_r=0$ (see Fig. 9).$

Figure 10.

Critical Alfvénic Mach number |$\mathcal {M}_{\rm{Ac}}$| found from the constant scale height mode calculations presented in Fig. 8 by considering the percent variations in frequency between modes calculated with the inner boundary conditions v_r = 0 and |$\mathrm{\partial} _rb_r=0$| (see Fig. 9).

Open in new tab Download slide

The only other free parameter, the spin of the black hole, has a weaker effect on the critical field strengths. Larger values of the parameter a allow r-modes to remain independent of the inner boundary in the presence of marginally stronger magnetic fields, and result in slightly steeper increases in frequency with |$\mathcal {M}_{\rm{A}z}$|⁠.

5.3.3 Cylindrical correspondence

In calculating trapped inertial modes in both cylindrical and stratified disc models, it is natural to ask how solutions found in the two regimes relate to one another. Latter, Fromang & Faure (2015) showed that local MRI channel modes calculated in the shearing box correspond to sections in the evanescent regions of much larger, radially global MRI modes. We consider an analogous relationship between r-modes in MHD discs with and without vertical density stratification.

With respect to an appropriate choice of vertical wavenumber, Fu & Lai (2009) adopted the relationship |$k_z\sim \sqrt{\epsilon }/H,$| with ε being of order unity. We find that in general this relationship holds both with and without the inclusion of radial variation in H. However, we go further, finding a precise value for the proportionality |$\sqrt{\epsilon }$|⁠. The rapidity of the convergence we find with our expansion in the basis functions F_n and G_n suggests the following identification, for k_z in units of |$r_g^{-1}$|⁠:

\begin{eqnarray} k_z =\frac{K_n}{H}. \end{eqnarray}

(61)

In particular we choose k_z = K₁/H ≈ 1.16/H to follow our interest in the vertically fundamental, k = 1 r-mode (cf. Fig. A1 and fig. 1 in Latter et al. 2010).

With this prescription, the frequencies for cylindrical modes closely follow those of fully global modes. Fig. 11(top) shows the percent difference between frequencies calculated for stratified r-modes and cylindrical r-modes with k_z = K₁/H, at sound speeds c_s = 0.001–0.01c, and a = 0 (similar results are obtained for a = 0.5, 0.8). In all cases, the frequencies differ by less than |$0.5\,\,\rm{per\,\,cent}$|⁠, but the deviation increases from |$0.01{\rm -}0.03\,\,\rm{per\,\,cent}$| to |$0.1{\rm -}0.3\,\,\rm{per\,\,cent}$| for larger sound speeds. For comparison, frequencies calculated with the prescription k_z = 1/H give percent differences as large as |$1{\rm -}10\,\,\rm{per\,\,cent}$|⁠, increasing with larger magnetic field strengths. Fig. 11(bottom) shows the percent differences between measurements of the radius of maximal v_r, r_max, between cylindrical and stratified calculations, which are also less than |$1\,\,\rm{per\,\,cent}$|⁠.

Figure 11.

Percent differences between frequencies (top) and radii of maximal v_r (bottom) calculated for fully global r-modes with a purely constant scale height (Figs 8 and 10) and cylindrical r-modes calculated with the same spin (a = 0), sound speeds and magnetic field strengths, and k_z = K₁/H (here Δ indicates the difference between quantities calculated using the two models). Marginally closer relationships were found for larger values of a.

Open in new tab Download slide

The close agreement between frequencies and localization found through fully global and unstratified calculations made with k_z = K₁/H suggests that the cylindrical model presented in Section 4 can accurately describe the behaviours of fully global, MHD r-modes. However, the correct choice of k_z is essential to exploit this property.

5.4 The effects of radial scale height variation

In this section, we include the radial variation of the isothermal scale height H = c_s/Ω_z, and thus present fully global, self-consistent calculations of MHD r-modes. We now must contend with non-zero profiles for |$\mathrm{\partial} _r\ln H$| and |$\mathcal {L}_{HH}$|⁠, which in turn provide further couplings in equations (54)–(60) through the coefficients ν_nm, λ_nm, and μ_nm. The integrals defining these coefficients (see Appendix A1) are larger than those defining θ_nm, α_nm, and γ_nm, but they couple with the equations only in concert with factors that are |$\lesssim\!\mathcal{O}(1/r)$|⁠, and we find that their inclusion has little effect on the r-modes. The largest effect comes from radial variation in the coefficients involving H⁻¹, which exacerbates the impact of the background magnetic field. However, the extent of this exacerbation is likely exaggerated by our simplifying assumption of a globally constant sound speed.

Introducing radial variation in disc thickness also complicates our choice of boundary conditions, this time at r_out. In hydrodynamic disc models, a vertical resonance occurs at the radius where |$\omega ^2=n\Omega _z^2$| (Ferreira & Ogilvie 2008), beyond which lies the propagation region for p-modes with frequency ω. Fu & Lai (2009) identified the analogous radius for an unstratified disc threaded by a purely constant vertical field as that where |$\omega ^2=k_z^2c_s^2$|⁠. Given the correspondence found in Section 5.3.3, we place our outer boundary beyond the radius at which |$\omega ^2=K_1^2\Omega _z^2$| and implement a wave propagation boundary to allow for ‘tunnelling’ of r-modes to the outer regions of the disc. At r_out, we require ∀n that |$\mathrm{\partial} _ru_n=\rm{i}k_ru_n,$| with k_r determined using k_z = K₁/H(r_out) and a frequency calculated for an r-mode using the modified cylindrical model described in Section 5.4.2. The frequencies calculated for these modes are then complex. However, we find that, as observed by Ferreira & Ogilvie (2008) in their hydrodynamic calculations (cf. their table 1), the decay rates are negligible for the range of sound speeds considered here.

5.4.1 Eigenmodes and critical field strengths when H = H(r)

The discrete spectrum of global r-modes calculated with the fully self-consistent inclusion of radial scale height variation is similar to that found with the approximation of a purely constant H. Truncating our series expansions at vertical order n = M results in M sets of r-modes, each with their own discrete spectrum of radial quantum numbers l. Qualitatively, the response of the r-modes to increasing background magnetic field strength is the same as that observed with a purely constant H. Once again solutions dominated by coefficients of larger vertical order n = k are pushed further towards the ISCO by a given increase in B_z (Fig. 12), and increases in frequency are comparable (Fig. 13).

$Real parts of the radial momentum perturbation pr = ρ00gvr (interpolated from r, η to r, z coordinates) for trapped inertial waves with radial quantum numbers l = 0, 1, 2 (from left to right) and vertical orders k = 1, 2, 3 (from top to bottom), calculated with cs = 0.001c, a = 0, r/rg ∈ [6, 20], N = 300, M = 7, midplane $\mathcal {M}_{\rm{A}z}=0.04$ (β ≈ 1250) and a radially variable scale height. The black dashed line indicates Rκ, while the white dashed lines mark z = ±H and z = ±2H.$

Figure 12.

Real parts of the radial momentum perturbation p_r = ρ₀₀gv_r (interpolated from r, η to r, z coordinates) for trapped inertial waves with radial quantum numbers l = 0, 1, 2 (from left to right) and vertical orders k = 1, 2, 3 (from top to bottom), calculated with c_s = 0.001c, a = 0, r/r_g ∈ [6, 20], N = 300, M = 7, midplane |$\mathcal {M}_{\rm{A}z}=0.04$| (β ≈ 1250) and a radially variable scale height. The black dashed line indicates R_κ, while the white dashed lines mark z = ±H and z = ±2H.

Open in new tab Download slide

Figure 13.

Frequencies for the fully general fundamental r-modes, calculated as in Fig. 8 but with the inclusion of radial scale height variation, as a function of Alfvénic Mach number, sound speed (colour-scale goes from black for c_s = 0.0005c to yellow for c_s = 0.004c), and black hole spin parameter a.

Open in new tab Download slide

However, the critical field strengths at which frequencies diverge for modes calculated with the boundary conditions v_r = 0 and |$\mathrm{\partial} _rb_r=0$| at r_in are lower by nearly a factor of 2 (Fig. 14). This arises through the dependence on disc thickness discussed in Section 5.3. For a given sound speed, the r-modes observe a thicker disc at the trapping region with scale height variation than without. Alternatively, one might look to equation (2). For the prescription k_z ∝ H⁻¹, the contribution to the restoring force provided by the Alfvén frequency k_zV_Az increases as the mode is forced inwards by the vertical magnetic field.

$Critical Alfvénic Mach number $\mathcal {M}_{\rm{Ac}}$ found from the variable scale height mode calculations presented in Fig. 13 by considering the percent variations in frequency between modes calculated with the inner boundary conditions vr = 0 and $\mathrm{\partial} _rb_r=0$.$

Figure 14.

Critical Alfvénic Mach number |$\mathcal {M}_{\rm{Ac}}$| found from the variable scale height mode calculations presented in Fig. 13 by considering the percent variations in frequency between modes calculated with the inner boundary conditions v_r = 0 and |$\mathrm{\partial} _rb_r=0$|⁠.

Open in new tab Download slide

The estimates of critical magnetic field strengths obtained in this section are likely exaggerated by the assumption of global isothermality, though, which results in more rapid flaring of the disc than might be found with more realistic temperature profiles. As mentioned in Section 3, however, radial temperature gradients significantly complicate the equilibrium flow, and we do not consider them here.

5.4.2 Cylindrical correspondence when H = H (r)

Fully global trapped inertial modes calculated with the inclusion of radial scale-height variation do not show a direct correspondence with r-modes calculated in the cylindrical model as presented in Section 4, due to the assumption of a purely constant vertical wavenumber k_z associated with the latter. Suppose now that |$k_z = k_z(r)$|⁠. For perturbations |$\delta$| with the dependence |$\delta(r,z)\propto \tilde{\delta}(r)\exp[\textrm{i}k_z(r)z]$|⁠, we then have

\begin{eqnarray} \frac{\mathrm{\partial} \delta }{\mathrm{\partial} r} ={\rm e}^{\rm{i}k_zz}\left( \frac{\rm{d}\tilde{\delta }}{\rm{d}r} +iz\frac{\rm{d}k_z}{\rm{d}r}\tilde{\delta } \right). \end{eqnarray}

(62)

The dependence on z of the second term on the right-hand side eliminates the advantage of separability for which one might assume a plane wave dependence in the first place. However, if in our cylindrical calculations we retain radial variation in using the prescription k_z = K₁/H(r) but assume that the scale height H varies slowly enough for dk_z/dr to be negligible, we find that the coupling of vertical modes is weak enough that the correspondence is still very close, although not as precise as that presented in Section 5.3.3.

As in Fig. 11, Fig. 15 shows the percent difference between frequencies calculated with increasing magnetic field strengths and varying sound speeds c_s/c = 0.0005–0.004 for stratified modes and cylindrical modes calculated with the k_z = K₁/H(r). Again the percent changes in frequency are less than |$1\,\,\rm{per\,\,cent}$|⁠, and percent changes in radii of maximal radial velocity less than |$2\,\,\rm{per\,\,cent}$|⁠, but the former increase more rapidly with sound speed than with the approximation of a purely constant scale height. This is perhaps due to the increasing relevance of the radial derivative ignored in equation (62). In terms of the scaled coordinate η, this prescription for the vertical wavenumber of a perturbation δ(r, η) might alternatively be written as |$\delta (r, \eta )=\tilde{\delta }(r)\exp [\rm{i}K_1\eta ],$| yielding

\begin{eqnarray} \frac{\mathrm{\partial} \delta }{\mathrm{\partial} r} ={\rm e}^{\rm{i}K_1\eta } \left( \frac{\rm{d}\tilde{\delta }}{\rm{d}r} -iK_1\eta \frac{\rm{d}\ln H}{\rm{\rm{d}}r}\tilde{\delta } \right). \end{eqnarray}

(63)

In implementing this prescription for a radially variable k_z(r) in our cylindrical, unstratified model, then, we are omitting the coupling of vertical modes that is included self-consistently in our fully global calculations, and this omission becomes more relevant in hotter, thicker discs.

Figure 15.

Percent differences between frequencies (top) and radii of maximal v_r (bottom) calculated for fully global r-modes with a radially variable scale height (Figs 13 and 14) and cylindrical r-modes calculated with the same spin (a = 0), sound speeds and magnetic field strengths, and k_z(r) = K₁/H(r) (a similar wave propagation boundary condition to that described in Section 5.4 was used for the cylindrical modes).

Open in new tab Download slide

6 DISCUSSION AND CONCLUSIONS

We have reconsidered the effects of poloidal and toroidal magnetic fields on trapped inertial waves in magnetised relativistic accretion discs by undertaking cylindrical and fully global linear eigenvalue calculations. Our work expands on and more thoroughly investigates a prediction made by local analyses (Fu & Lai 2009): that a purely constant, vertical magnetic field forces r-modes inward, possibly removing their independence from the inner edge of the disc.

We have verified that poloidal magnetic fields do force r-modes to migrate inwards, and established more realistic estimates for the critical magnetic field strength at which the modes are pushed up against the inner disc edge. However, we find that at lower temperatures MHD inertial modes can remain independent from the inner boundary, even in the presence of relatively strong magnetic fields; the cooler and thinner the disc, the less susceptible the modes are to the fields. We also find that r-mode frequencies can nearly double because of the magnetic field (see Figs 8 and 13). This frequency enhancement should be kept in mind when using QPOs to estimate properties of the black hole-disc system.

Normal modes calculated with a cylindrical model confirm that an azimuthal field of even equipartition strength has little effect on r-mode frequencies or localizations. Radial gradients in magnetic fields or density, for reasonable power-law indices, also have negligible effects. Finally, we have shown that with an appropriate choice of vertical wavenumber, trapped inertial modes calculated in a simplified cylindrical framework (ideal for global non-linear simulations) closely reproduce the frequencies and localizations of fully global r-modes (to within 1 per cent).

For an inner accretion disc temperature of ∼1 keV, isothermal sound speed estimates of

\begin{eqnarray} c_s=\sqrt{\frac{k_BT}{\mu m_p}} \approx 0.001c \end{eqnarray}

(64)

lead to the survival of the r-mode trapping region in the presence of vertical magnetic fields with midplane Alfvénic Mach numbers less than |$\mathcal {M}_{\rm{A}z}\approx 0.08{\rm -}0.14$| (see Fig 10 and 14). It should be emphasized that these critical field strengths should be regarded as lower bounds, given the stringency of the criterion used to find them. Moreover, a loss of isolation from the inner boundary need not preclude the existence of trapped inertial waves, given the possibility of wave reflection by steep density gradients or other features. All one can say is that r-modes that are forced up against the ISCO require greater excitation to counteract enhanced energy losses through the inner boundary.

We note that beyond assuming an isothermal equation of state, we have applied other simplifications to our model, in order to more feasibly investigate r-mode propagation in a fully global, MHD context. One such simplification is the omission of a background radial inflow, which intensifies near the ISCO. Ferreira (2010) considered the effects of radial inflow on trapped inertial modes, and found that decay rates were amplified by the presence of a transonic accretion flow in a viscous disc model. However, the author found that the severity of the damping depends on the location of the sonic point, a location r_sonic < r_ISCO still allowing r-mode excitation by warps and eccentricities in a hydrodynamic disc. The movement of the trapping region towards the inner boundary by a poloidal magnetic field would make damping by radial inflow more likely, but not a certainty for sonic points far enough within the ISCO.

Additionally, although the range of sound speeds considered is motivated by observations, radiation pressure is expected to play a large role in the black hole accretion discs associated with X-ray binaries in the emission states in which HFQPOs are primarily observed. The thickening of the disc associated with this radiation pressure might nullify the increase in critical magnetic field strength we have observed for discs with lower sound speeds, although the inclusion of radiative transfer may modify the dynamics in other ways.

We have also excluded the impact of radial magnetic fields, which some have suggested might counteract the disruptive effects of a purely vertical one (Ortega-Rodriguez et al. 2015), as might disc truncation above and below by a hot corona (Kato 2017). Further, we have considered only the effects of large-scale, ordered magnetic fields, for which midplane plasma betas of β ≲ 500 are actually rather strong. Although the MRI may saturate near equipartition on small scales, the net flux associated with large-scale poloidal fields with the strengths considered here strongly affects MHD turbulence and outflows in accretion disc simulations (e.g. Lesur, Ferreira & Ogilvie 2013; Bai & Stone 2014; Salvesen et al. 2016).

Our final take away point is the following. The appearance or not of trapped inertial waves in observations may be determined more by a competition between excitation by disc deformations, and damping by radial inflow and wave leakage, rather than the effect of a large-scale magnetic field on the self-trapping region alone. Models and non-linear simulations including more complicated flow dynamics, thermodynamics, and magnetic field configurations will certainly shed more light on the issue and form the basis for future work.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous reviewer for a positive and useful report. JWD thanks the Cambridge Commonwealth European and International Trust and the Vassar College De Golier Trust for funding this work.

Footnotes

1

Ortega-Rodriguez & Wagoner (2000) concluded incorrectly (because of a sign error in their analysis) that many hydrodynamic oscillation modes are destabilized by viscosity, but in fact they are damped.

2

The modes subject to confinement within the regions defined by equation (2) are epicyclic-Alfvénic in nature, but we refer to both them and purely hydrodynamic trapped inertial waves as r-modes throughout this paper.

3

We note that the profiles for |$-k_r^2$| found from this dispersion relation can deviate based on whether or not the classic equality |$d\Omega/d\ln r = \kappa^2/(2\Omega)-2\Omega$| is used in conjunction with the general relativistic versions of the characteristic frequencies. Although this equality was assumed in deriving equations (14)–(20), it is not accurate for characteristic frequencies |$\kappa_G,$ $\Omega_G$| derived from the Kerr metric (equations 10–12). In solving equations (14)–(20) with the general relativistic frequencies, we equate the direct formulae for |$\kappa_G^2$| and |$d\Omega_G/d\ln r$| with the epicyclic and shearing terms appearing in equations (15) and (19) (resp.).

REFERENCES

Afshordi

N.

,

Paczynski

B.

,

2003

,

ApJ

,

592

,

354

10.1086/375559

Crossref

Search ADS

Appert

K.

,

Gruber

R.

,

Vaclavik

J.

,

1974

,

Phys. Fluids

,

17

,

1471

Crossref

Search ADS

Arras

P.

,

Blaes

O.

,

Turner

N.

,

2006

,

ApJ

,

645

,

65

10.1086/505966

Crossref

Search ADS

Bai

X.

,

Stone

J.

,

2014

,

ApJ

,

796

,

14

10.1088/0004-637X/796/1/31

Crossref

Search ADS

Balbus

S.

,

Hawley

J.

,

1998

,

Rev. Mod. Phys.

,

70

,

1

10.1103/RevModPhys.70.1

Crossref

Search ADS

Boyd

J.

,

2001

,

Chebyshev and Fourier Spectral Methods

.

Dover Publications

,

Mineola, NY

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Ferreira

B.

,

2010

,

PhD thesis

,

Univ. Cambridge

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Ferreira

B.

,

Ogilvie

G.

,

2008

,

MNRAS

,

386

,

2297

10.1111/j.1365-2966.2008.13207.x

Crossref

Search ADS

Ferreira

B.

,

Ogilvie

G.

,

2009

,

MNRAS

,

392

,

428

10.1111/j.1365-2966.2008.14070.x

Crossref

Search ADS

Fu

W.

,

Lai

D.

,

2009

,

ApJ

,

690

,

1386

10.1088/0004-637X/690/2/1386

Crossref

Search ADS

Fu

W.

,

Lai

D.

,

2011

,

MNRAS

,

410

,

399

.

10.1111/j.1365-2966.2010.17451.x

Crossref

Search ADS

Gammie

C.

,

1999

,

ApJ

,

522

,

57

10.1086/312207

Crossref

Search ADS

Gammie

C.

,

Balbus

S.

,

1994

,

MNRAS

,

270

,

138

10.1093/mnras/270.1.138

Crossref

Search ADS

Goodman

J.

,

1993

,

ApJ

,

406

,

596

10.1086/172472

Crossref

Search ADS

Henisey

K.

,

Blaes

O.

,

Fragile

P.

,

Ferreira

B.

,

2009

,

ApJ

,

706

,

705

10.1088/0004-637X/706/1/705

Crossref

Search ADS

Kato

S.

,

2001

,

PASJ

,

53

,

1

10.1093/pasj/53.1.1

Crossref

Search ADS

Kato

S.

,

2004

,

PASJ

,

56

,

905

10.1093/pasj/56.5.905

Crossref

Search ADS

Kato

S.

,

2008

,

PASJ

,

60

,

889

10.1093/pasj/60.4.889

Crossref

Search ADS

Kato

S.

,

2017

,

MNRAS

,

472

,

1119

10.1093/mnras/stx2075

Crossref

Search ADS

Lai

D.

,

Tsang

D.

,

2009

,

MNRAS

,

393

,

979

10.1111/j.1365-2966.2008.14218.x

Crossref

Search ADS

Lai

D.

,

Fu

W.

,

Tsang

D.

,

Horak

J.

,

Yu

C.

,

2013

, in

Zhang

C. M.

et al. , eds,

IAU Symp. Vol. 290, Feeding Compact Objects: Accretion on all Scales

.

Cambridge Univ. Press

,

Cambridge

, p.

57

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Latter

H.

,

Fromang

S.

,

Gressel

O.

,

2010

,

MNRAS

,

406

,

848

Latter

H.

,

Fromang

S.

,

Faure

J.

,

2015

,

MNRAS

,

453

,

3257

10.1093/mnras/stv1890

Crossref

Search ADS

Lesur

G.

,

Ferreira

J.

,

Ogilvie

G. I.

,

2013

,

A&A

,

550

,

13

10.1051/0004-6361/201220395

Crossref

Search ADS

Li

L.

,

Goodman

J.

,

Narayan

R.

,

2003

,

ApJ

,

593

,

980

10.1086/376695

Crossref

Search ADS

Motta

S.

,

2016

,

Astron. Nachr.

,

337

,

398

10.1002/asna.201612320

Crossref

Search ADS

Narayan

R.

,

Goldreich

P.

,

Goodman

J.

,

1987

,

MNRAS

,

228

,

1

.

10.1093/mnras/228.1.1

Crossref

Search ADS

O'Neill

S.

,

Reynolds

C.

,

Miller

C.

,

2009

,

ApJ

,

693

,

1100

10.1088/0004-637X/693/2/1100

Crossref

Search ADS

Ogilvie

G. I.

,

1998

,

MNRAS

,

297

,

291

10.1046/j.1365-8711.1998.01507.x

Crossref

Search ADS

Ogilvie

G. I.

,

2008

,

MNRAS

,

388

,

1372

Ogilvie

G. I.

,

Pringle

J.

,

1996

,

MNRAS

,

279

,

152

10.1093/mnras/279.1.152

Crossref

Search ADS

Okazaki

A. T.

,

Kato

S.

,

Fukue

J.

,

1987

,

PASJ

,

39

,

457

Ortega-Rodriguez

M.

,

Wagoner

R.

,

2000

,

ApJ

,

537

,

922

10.1086/309060

Crossref

Search ADS

Ortega-Rodriguez

M.

,

Solis-Sanchez

H.

,

Wagoner

R.

,

Arguedas-Leiva

A.

,

Levine

A.

,

2015

,

ApJ

,

809

,

15

10.1088/0004-637X/809/1/15

Crossref

Search ADS

Papaloizou

J. C. B.

,

Pringle

J. E.

,

1984

,

MNRAS

,

208

,

721

10.1093/mnras/208.4.721

Crossref

Search ADS

Remillard

R. A.

,

McClintock

J. E.

,

2006

,

ARA&A

,

44

,

49

10.1146/annurev.astro.44.051905.092532

Crossref

Search ADS

Reynolds

C.

,

Miller

C.

,

2009

,

ApJ

,

692

,

869

10.1088/0004-637X/692/1/869

Crossref

Search ADS

Ruy

D.

,

Goodman

J.

,

1994

,

ApJ

,

422

,

269

10.1086/173725

Crossref

Search ADS

Salvesen

G.

,

Simon

J. B.

,

Armitage

P. J.

,

Begelman

M. C.

,

2016

,

MNRAS

,

457

,

857

10.1093/mnras/stw029

Crossref

Search ADS

Sano

T.

,

Miyama

S.

,

1999

,

ApJ

,

515

,

776

10.1086/307063

Crossref

Search ADS

Yu

C.

,

Lai

D.

,

2015

,

MNRAS

,

450

,

2466

10.1093/mnras/stv803

Crossref

Search ADS

APPENDIX A: NUMERICAL METHOD FOR STRATIFIED CALCULATIONS

A1 Coupling coefficients

Substituting the series expansions discussed in Section 5.2 in equations (32)–(38), multiplying by suitable combinations of basis functions F_n and G_n (see Fig. A1) of arbitrary order and integrating over |$\int _{-\infty }^{\infty }\, \mathrm{d}\eta$|⁠, the orthogonality relations (48)–(49) allow for the elimination of many of the infinite sums. This projection on to a given vertical order n is marred only by sums over the other vertical orders m that involve the coupling integrals

\begin{eqnarray} \theta _{nm} =K_m\int _{-\infty }^{\infty }g^{3/4}F_{n-1}G_m\, \mathrm{d}\eta , \end{eqnarray}

(A1)

\begin{eqnarray} \lambda _{nm} =\int _{-\infty }^{\infty } gF_m\left(F_n+\eta K_nG_n\right)\, \mathrm{d}\eta , \end{eqnarray}

(A2)

\begin{eqnarray} \alpha _{nm} =K_n\int _{-\infty }^{\infty }g^{5/4}F_{m-1}G_n\, \mathrm{d}\eta , \end{eqnarray}

(A3)

\begin{eqnarray} \gamma _{nm} =K_m^2\int _{-\infty }^{\infty }g^2F_nF_m\, \mathrm{d}\eta , \end{eqnarray}

(A4)

\begin{eqnarray} \nu _{nm} =K_m\int _{-\infty }^{\infty }\eta gF_nG_m\, \mathrm{d}\eta , \end{eqnarray}

(A5)

\begin{eqnarray} \mu _{nm} =K_m^2\int _{-\infty }^{\infty }\eta ^2 g^2F_nF_m\, \mathrm{d}\eta . \end{eqnarray}

(A6)

These integrals are evaluated by first calculating the basis functions F_n, G_n and eigenvalues K_n, then integrating numerically over the range of η within which the normalized integrands are >10⁻⁴. The numerical values for the integrals computed with n, m ∈ [1, 30] are shown in the heatmaps in Fig. A2. For small n, the integrals go to zero very quickly with increasing m, such that the couplings are very weak for equations with very different vertical order. Additionally, the heatmaps in Fig. A2 illustrate the symmetry in equations (32)–(38) ensuring that mode couplings are restricted to either even or odd parity.

$The n = 1, 2, 3, 4, 5 basis functions Gn (top) and Fn (bottom), normalized such that $\int _{-\infty }^{\infty }G_nG_m\, \mathrm{d}\eta = \delta _{nm}$ and $\int _{-\infty }^{\infty }gF_nF_m\, \mathrm{d}\eta = \delta _{nm}$. The latter are found from the former, which are calculated using pseudo-spectral methods utilizing Whittaker cardinal functions (compare with fig. 1 in Latter et al. 2010).$

Figure A1.

The n = 1, 2, 3, 4, 5 basis functions G_n (top) and F_n (bottom), normalized such that |$\int _{-\infty }^{\infty }G_nG_m\, \mathrm{d}\eta = \delta _{nm}$| and |$\int _{-\infty }^{\infty }gF_nF_m\, \mathrm{d}\eta = \delta _{nm}$|⁠. The latter are found from the former, which are calculated using pseudo-spectral methods utilizing Whittaker cardinal functions (compare with fig. 1 in Latter et al. 2010).

Open in new tab Download slide

Figure A2.

Heatmaps showing the values computed numerically for the coupling integrals θ_nm, α_nm, γ_nm, λ_nm, ν_nm, and μ_nm. The heatmaps are saturated because the amplitudes of the coupling integrals on the diagonal do not matter (i.e. the n = 1 equations clearly ought to couple very strongly with the m = 1 equations).

Open in new tab Download slide

A2 Numerical convergence

To check for convergence with increasing series truncation order M, we track both the normal mode frequencies and the relative energetic contributions from each |${{\boldsymbol U}}_n$| to the full solution |${{\boldsymbol U}}.$| For |f|² = ff*, we quantify the latter through norms defined for each variable, at each vertical order n, as

\begin{eqnarray} v_r:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}\rho _{00}g|u_nF_n|^2\, \mathrm{d}r\, \mathrm{d}\eta , \end{eqnarray}

(A7)

\begin{eqnarray} v_{\phi }:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}\rho _{00}g|v_nF_n|^2\, \mathrm{d}r\, \mathrm{d}\eta , \end{eqnarray}

(A8)

\begin{eqnarray} v_z:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}\rho _{00}g|w_ng^{1/4}F_{n-1}|^2\, \mathrm{d}r\, \mathrm{d}\eta , \end{eqnarray}

(A9)

\begin{eqnarray} \Gamma :\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}|c_s^2\Gamma _nF_n|^2\, \mathrm{d}r\mathrm{d}\eta , \end{eqnarray}

(A10)

\begin{eqnarray} b_r:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}|b_r^nG_n|^2\, \mathrm{d}r\mathrm{d}\eta , \end{eqnarray}

(A11)

\begin{eqnarray} L_r:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}|L_r^ngF_n|^2\, \mathrm{d}r\mathrm{d}\eta , \end{eqnarray}

(A12)

\begin{eqnarray} L_{\phi }:\qquad \int _{-\infty }^{\infty }\int _{r_{in}}^{r_{\rm{out}}}|L_{\phi }^ngF_n|^2\, \mathrm{d}r\mathrm{d}\eta . \end{eqnarray}

(A13)

With the expansions discussed in Section 5.2, the energetic contributions from the non-dominant n are at least two orders of magnitude smaller than that from the dominant n, and quickly fall off with increasing series truncation to values smaller than the error introduced by our Chebyshev collocation method. Beyond truncation at M = 5, changes in frequency for the fundamental k = 1 r-mode are less than |$10^{-4}\,\,\rm{per\,\,cent}$|⁠, and continue to fall with increasing M. Frequencies found with increasing M and |$\mathcal {M}_{\rm{A}z}$| are shown for in Tables A1 and A2 for constant and variable scale height calculations (respectively). Like Table A2, Table A3 shows frequencies calculated with variable H(r) but with the coupling provided by the radial derivatives of H excluded.

Table A1.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated with a purely constant scale height and the boundary condition |$\mathrm{\partial} _rb_r=0$|⁠, for varying magnetic field strength and series truncation order M (a = 0, c_s = 0.001, r/r_g ∈ [6, 20], N = 200).

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.04	0.06	0.08	0.1	0.12	0.14
∼β	5000	1250	556	312	200	139	102
M = 1	0.022 789	0.024 658	0.027 08	0.029 722	0.032 439	0.035 166	0.037 882
M = 3	0.022 791	0.024 66	0.027 082	0.029 726	0.032 443	0.035 17	0.037 886
M = 5	0.022 791	0.024 661	0.027 083	0.029 726	0.032 443	0.035 171	0.037 887
M = 7	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 9	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 11	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 172	0.037 887
M = 13	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888
M = 15	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.04	0.06	0.08	0.1	0.12	0.14
∼β	5000	1250	556	312	200	139	102
M = 1	0.022 789	0.024 658	0.027 08	0.029 722	0.032 439	0.035 166	0.037 882
M = 3	0.022 791	0.024 66	0.027 082	0.029 726	0.032 443	0.035 17	0.037 886
M = 5	0.022 791	0.024 661	0.027 083	0.029 726	0.032 443	0.035 171	0.037 887
M = 7	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 9	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 11	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 172	0.037 887
M = 13	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888
M = 15	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888

Table A1.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated with a purely constant scale height and the boundary condition |$\mathrm{\partial} _rb_r=0$|⁠, for varying magnetic field strength and series truncation order M (a = 0, c_s = 0.001, r/r_g ∈ [6, 20], N = 200).

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.04	0.06	0.08	0.1	0.12	0.14
∼β	5000	1250	556	312	200	139	102
M = 1	0.022 789	0.024 658	0.027 08	0.029 722	0.032 439	0.035 166	0.037 882
M = 3	0.022 791	0.024 66	0.027 082	0.029 726	0.032 443	0.035 17	0.037 886
M = 5	0.022 791	0.024 661	0.027 083	0.029 726	0.032 443	0.035 171	0.037 887
M = 7	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 9	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 11	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 172	0.037 887
M = 13	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888
M = 15	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.04	0.06	0.08	0.1	0.12	0.14
∼β	5000	1250	556	312	200	139	102
M = 1	0.022 789	0.024 658	0.027 08	0.029 722	0.032 439	0.035 166	0.037 882
M = 3	0.022 791	0.024 66	0.027 082	0.029 726	0.032 443	0.035 17	0.037 886
M = 5	0.022 791	0.024 661	0.027 083	0.029 726	0.032 443	0.035 171	0.037 887
M = 7	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 9	0.022 791	0.024 661	0.027 083	0.029 727	0.032 444	0.035 172	0.037 887
M = 11	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 172	0.037 887
M = 13	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888
M = 15	0.022 791	0.024 661	0.027 084	0.029 727	0.032 444	0.035 173	0.037 888

Table A2.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated with a variable scale height and the boundary condition |$\mathrm{\partial} _rb_r=0$|⁠, for varying magnetic field strength and series truncation order M (a = 0, c_s = 0.001, r/r_g ∈ [6, 20], N = 200).

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 174	0.025 184	0.026 431	0.027 975	0.029 662
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 187	0.026 433	0.027 977	0.029 667
M = 5	0.022 347	0.022 773	0.023 379	0.024 178	0.025 188	0.026 434	0.027 977	0.029 668
M = 7	0.022 347	0.022 773	0.023 379	0.024 179	0.025 188	0.026 434	0.027 977	0.029 669
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 67
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 978	0.029 67

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 174	0.025 184	0.026 431	0.027 975	0.029 662
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 187	0.026 433	0.027 977	0.029 667
M = 5	0.022 347	0.022 773	0.023 379	0.024 178	0.025 188	0.026 434	0.027 977	0.029 668
M = 7	0.022 347	0.022 773	0.023 379	0.024 179	0.025 188	0.026 434	0.027 977	0.029 669
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 67
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 978	0.029 67

Table A2.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated with a variable scale height and the boundary condition |$\mathrm{\partial} _rb_r=0$|⁠, for varying magnetic field strength and series truncation order M (a = 0, c_s = 0.001, r/r_g ∈ [6, 20], N = 200).

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 174	0.025 184	0.026 431	0.027 975	0.029 662
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 187	0.026 433	0.027 977	0.029 667
M = 5	0.022 347	0.022 773	0.023 379	0.024 178	0.025 188	0.026 434	0.027 977	0.029 668
M = 7	0.022 347	0.022 773	0.023 379	0.024 179	0.025 188	0.026 434	0.027 977	0.029 669
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 67
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 978	0.029 67

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 174	0.025 184	0.026 431	0.027 975	0.029 662
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 187	0.026 433	0.027 977	0.029 667
M = 5	0.022 347	0.022 773	0.023 379	0.024 178	0.025 188	0.026 434	0.027 977	0.029 668
M = 7	0.022 347	0.022 773	0.023 379	0.024 179	0.025 188	0.026 434	0.027 977	0.029 669
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 669
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 434	0.027 978	0.029 67
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 978	0.029 67

Table A3.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated as in Table A2, but with the coupling provided by the radial derivatives of H excluded.

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 175	0.025 184	0.026 431	0.027 977	0.029 665
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 188	0.026 434	0.027 978	0.029 671
M = 5	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 7	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 175	0.025 184	0.026 431	0.027 977	0.029 665
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 188	0.026 434	0.027 978	0.029 671
M = 5	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 7	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673

Table A3.

Open in new tab

Frequencies Re[ω]/ω_g for the fundamental l = 0, k = 1 trapped inertial mode, calculated as in Table A2, but with the coupling provided by the radial derivatives of H excluded.

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 175	0.025 184	0.026 431	0.027 977	0.029 665
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 188	0.026 434	0.027 978	0.029 671
M = 5	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 7	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673

\|$\mathcal {M}_{\rm{A}z}$\|	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09
∼β	5000	2222	1250	800	555	408	312	246
M = 1	0.022 342	0.022 768	0.023 375	0.024 175	0.025 184	0.026 431	0.027 977	0.029 665
M = 3	0.022 346	0.022 772	0.023 379	0.024 178	0.025 188	0.026 434	0.027 978	0.029 671
M = 5	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 7	0.022 347	0.022 773	0.023 38	0.024 179	0.025 188	0.026 435	0.027 979	0.029 672
M = 9	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 11	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 13	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673
M = 15	0.022 347	0.022 773	0.023 38	0.024 179	0.025 189	0.026 435	0.027 979	0.029 673

It should be noted that we have more trouble resolving the coefficients |${{\boldsymbol U}}_n$| and frequencies at lower magnetic field strengths (⁠|$\mathcal {M}_{\rm{A}z}\lesssim 0.04$|⁠), when larger sound speeds (c_s ≳ 0.007c for constant H and c_s ≳ 0.003c for variable H) are used. Beyond the largest values of c_s considered in Sections 5.3 and 5.4 (0.01c and 0.004c, respectively), these resolution issues extend throughout the full range of magnetic field strengths considered. This indicates that the series expansions given in Section 5.2 are less appropriate in thicker discs, and highlights the importance of disc thinness to r-mode trapping in magnetized contexts.

Download all slides

Month:	Total Views:
February 2018	22
March 2018	13
April 2018	39
May 2018	13
June 2018	4
July 2018	16
August 2018	9
September 2018	12
October 2018	21
November 2018	10
December 2018	14
January 2019	9
February 2019	5
March 2019	2
April 2019	13
May 2019	6
June 2019	3
August 2019	2
October 2019	3
November 2019	10
December 2019	2
January 2020	1
May 2020	4
June 2020	8
July 2020	14
August 2020	1
September 2020	2
December 2020	1
March 2021	1
May 2021	3
June 2021	11
July 2021	5
August 2021	5
September 2021	7
October 2021	4
November 2021	5
December 2021	7
January 2022	5
February 2022	5
March 2022	4
April 2022	4
May 2022	5
June 2022	8
July 2022	14
August 2022	14
September 2022	13
October 2022	9
November 2022	7
December 2022	5
January 2023	11
February 2023	8
March 2023	18
April 2023	8
May 2023	5
June 2023	6
July 2023	9
August 2023	15
September 2023	8
October 2023	8
November 2023	8
December 2023	14
January 2024	9
February 2024	6
March 2024	12
April 2024	20
May 2024	12
June 2024	10
July 2024	22
August 2024	8
September 2024	9
October 2024	5
November 2024	5
December 2024	10
January 2025	4
February 2025	9
March 2025	19
April 2025	7
May 2025	13

Article Contents

Quasi-periodic oscillations and the global modes of relativistic, MHD accretion discs

Abstract

1 INTRODUCTION

2 BACKGROUND

2.1 Hydrodynamic waves in relativistic discs

2.2 Magnetohydrodynamic waves in relativistic discs

3 DISC MODEL AND LINEARIZED EQUATIONS

4 CYLINDRICAL CALCULATIONS

4.1 Cylindrical, hydrodynamic r-modes

4.2 Cylindrical, magnetohydrodynamic r-modes

4.2.1 Magnetohydrodynamic boundary conditions

4.2.2 Purely vertical magnetic field, |$B_{\phi }(r)=0,B_z(r)\not=0$|

4.2.3 Purely toroidal magnetic field, |$B_{\phi }(r)\not=0,B_z(r)=0$|

5 VERTICALLY STRATIFIED, FULLY GLOBAL CALCULATIONS

5.1 Stratified equations: a change of variables

5.2 Boundary conditions and series expansions

5.3 Global r-modes in a flat disc

5.3.1 The global r-mode spectrum

5.3.2 Critical magnetic field strengths

5.3.3 Cylindrical correspondence

5.4 The effects of radial scale height variation

5.4.1 Eigenmodes and critical field strengths when H = H(r)

5.4.2 Cylindrical correspondence when H = H (r)

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: NUMERICAL METHOD FOR STRATIFIED CALCULATIONS

A1 Coupling coefficients

A2 Numerical convergence

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Quasi-periodic oscillations and the global modes of relativistic, MHD accretion discs Free

Abstract

1 INTRODUCTION

2 BACKGROUND

2.1 Hydrodynamic waves in relativistic discs

2.2 Magnetohydrodynamic waves in relativistic discs

3 DISC MODEL AND LINEARIZED EQUATIONS

4 CYLINDRICAL CALCULATIONS

4.1 Cylindrical, hydrodynamic r-modes

4.2 Cylindrical, magnetohydrodynamic r-modes

4.2.1 Magnetohydrodynamic boundary conditions

4.2.2 Purely vertical magnetic field, |$B_{\phi }(r)=0,B_z(r)\not=0$|

4.2.3 Purely toroidal magnetic field, |$B_{\phi }(r)\not=0,B_z(r)=0$|

5 VERTICALLY STRATIFIED, FULLY GLOBAL CALCULATIONS

5.1 Stratified equations: a change of variables

5.2 Boundary conditions and series expansions

5.3 Global r-modes in a flat disc

5.3.1 The global r-mode spectrum

5.3.2 Critical magnetic field strengths

5.3.3 Cylindrical correspondence

5.4 The effects of radial scale height variation

5.4.1 Eigenmodes and critical field strengths when H = H(r)

5.4.2 Cylindrical correspondence when H = H (r)

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: NUMERICAL METHOD FOR STRATIFIED CALCULATIONS

A1 Coupling coefficients

A2 Numerical convergence

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Quasi-periodic oscillations and the global modes of relativistic, MHD accretion discs