Ogilvie Magnetized Accretion Discs 1998

n a J

Printed 24 May 2021

(MN plain TEX macros v1.6)

Mon. Not. R. Astron. Soc. 000, 000-000 (1994)

G. I. Ogilvie1,2 1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA 2Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Silver Street, Cambridge CB3 9EW

Waves and instabilities in a diﬀerentially rotating disc containing a poloidal magnetic field

ABSTRACT The theory of waves and instabilities in a diﬀerentially rotating disc containing a poloidal magnetic field is developed within the framework of ideal magnetohydrody- namics. A continuous spectrum, for which the eigenfunctions are localized on indi- vidual magnetic surfaces, is identified but is found not to contain any instabilities associated with diﬀerential rotation. The normal modes of a weakly magnetized thin disc are studied by extending the asymptotic methods used previously to describe the equilibria. Waves propagate radially in the disc according to a dispersion relation which is determined by solving an eigenvalue problem at each radius. The dispersion relation for a hydrodynamic disc is re-examined and the modes are classified according to their behaviour in the limit of large wavenumber. The addition of a magnetic field introduces new, potentially unstable, modes and also breaks up the dispersion diagram by causing avoided crossings. The stability boundary to the magnetorotational insta- bility in the parameter space of polytropic equilibria is located by solving directly for marginally stable equilibria. For a given vertical magnetic field in the disc, bending of the field lines has a stabilizing eﬀect and it is shown that stable equilibria exist which are capable of launching a predominantly centrifugally driven wind.

There are good reasons for believing that magnetic fields are important to the physics of accretion discs. Magnetohydrodynamic (MHD) mechanisms provide the most convincing explanations for the anomalous transport of angular momentum that is required for accretion to proceed. One possibility is that angular momentum is removed from the disc by a rotating MHD wind (Blandford & Payne 1982). These ﬂows have the property of collimating into jets perpendicular to the disc (Heyvaerts & Norman 1989), which is attractive in view of the observed association of discs and jets. A rather more convincing mechanism for angular momentum transport is provided by the magnetorotational instability (Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991). A Keplerian shear ﬂow, while hydrodynamically stable, is destabilized by the presence of a magnetic field, provided that the magnetic field is suﬃciently weak and the disc is suﬃciently ionized. The instability is linear, dynamical and operates for toroidal as well as poloidal magnetic fields (Foglizzo & Tagger 1995; Ogilvie & Pringle 1996; Terquem & Papaloizou 1996). In local simulations, the non-linear development leads to turbulence, and the associated Reynolds and Maxwell stresses transport angular momentum radially outwards, as is required for accretion (Brandenburg et al. 1995, 1996; Hawley, Gammie & Balbus 1995, 1996; Stone et al. 1996).

In an earlier work (Ogilvie 1997; hereafter, Paper I) some idealized equilibrium models of accretion discs containing magnetic fields were presented, with the aim of studying the magnetorotational instability in more realistic geometry. The model system consists of a perfectly conducting, non-self-gravitating ﬂuid in diﬀerential rotation about a massive central object. The ﬂuid contains a purely poloidal magnetic field of dipolar symmetry, which bends as it passes through the disc, enforcing isorotation on magnetic surfaces. This rather general model is governed by a non-linear, elliptic partial diﬀerential equation in two dimensions, a version of the Grad-Shafranov equation. Two classes of special solutions of this general problem were described.

First, asymptotic solutions were obtained in the limit of a thin disc, using as the small parameter ǫ a characteristic value of H(r)/r, where H(r) is the height of the upper surface of the disc above the equatorial plane at radius r. It was shown that two families of solutions resembling accretion discs exist in this limit, with diﬀerent asymptotic scalings representing diﬀerent balances of forces in the radial and vertical directions. The weakly magnetized discs are the natural generalization of the standard, hydrodynamic thin discs (Pringle 1981). The sound speed and Alfv´en velocity are both O(ǫ) [given that the

Key words: accretion, accretion discs - hydrodynamics - instabilities - MHD - waves.

/ h p

o r t s a : v i X r a

INTRODUCTION

c(cid:13) 1994 RAS

2 G. I. Ogilvie

Secondly, solutions were obtained by assuming self-similarity in the spherical radial coordinate, as is often done in the analysis of winds and jets from accretion discs. This is a convenient method of studying thick equilibria and also of verifying the asymptotic results for thin discs.

Although the accretion ﬂow, the resistivity of the ﬂuid, and any toroidal magnetic field are neglected in this model, it is expected that these additional terms constitute only small perturbations to the internal equilibrium of the disc. It was shown in Paper I that a model of a wind-driven accretion disc can be built up by superimposing these additional features on the solution for a weakly magnetized thin disc without disturbing the equilibrium at leading order in ǫ.

Keplerian velocity is O(1)], while the fractional deviation from Keplerian rotation, caused by the radial Lorentz force, is O(ǫ). The strongly magnetized discs are not directly related to hydrodynamic thin discs but resemble models of solar prominences in which a sheet of matter is supported against gravity by a bending magnetic field (Kippenhahn & Schl¨uter 1957). The sound speed and Alfv´en velocity are both O(ǫ1/2), while the fractional deviation from Keplerian rotation is O(1). Previous approaches based on a vertical integration of the equations (e.g. Heyvaerts & Priest 1989) described only the strongly magnetized discs. It is, however, the weakly magnetized discs that are capable of launching a predominantly centrifugally driven wind if the magnetic field lines at the surface of the disc are inclined to the vertical by an angle greater than π/6.

The primary purpose of this paper is to determine the spectrum of waves and instabilities in a weakly magnetized thin disc, within the framework of ideal MHD, by extending the asymptotic methods used to describe the equilibria. In particular, the stability of thin discs to the magnetorotational instability is to be studied. There may, of course, be other types of instability in magnetized accretion discs. A non-axisymmetric interchange instability (Spruit, Stehle & Papaloizou 1995) or bending instability (Agapitou, Papaloizou & Terquem 1997) may be present if the magnetic field provides a significant amount of support against gravity in the radial direction. A poloidal magnetic field is also expected to alter the criterion for convective instability (Moss & Tayler 1969). Finally, a global, non-axisymmetric instability, related to the Papaloizou-Pringle instability (Papaloizou & Pringle 1984) may be anticipated (Curry & Pudritz 1996). However, not all of these instabilities will feature in the analysis that follows, either because they are not present in weakly magnetized discs, or because their growth rates are too small to be detected.

In an accretion disc the diﬀerential rotation is the most important feature of the dynamics and cannot be ignored or approximated by uniform rotation. The energy principle of Bernstein et al. (1958), which applies only to static equilibria, cannot be used. Nevertheless, Papaloizou & Szuszkiewicz (1992) showed that, for a diﬀerentially rotating, non-self-gravitating ﬂuid containing a poloidal magnetic field, a generalization of the energy principle exists in the form of a variational principle for the frequency eigenvalues of axisymmetric normal modes. Using the variational principle, they deduced some suﬃcient conditions for stability to axisymmetric perturbations. In this case, however, the problem does not reduce to a consideration of each magnetic surface separately. This suggests that the most unstable (or least stable) part of the spectrum of axisymmetric modes is discrete rather than continuous. This is consistent with the fact that the axisymmetric magnetorotational instability, although often described as a local instability, requires a relatively long wavelength in the direction perpendicular to the magnetic field, rather than being localized on a single magnetic surface. Nevertheless, it is possible to analyse the continuous spectrum directly without reference to an energy principle.

This paper is concerned principally with the weakly magnetized thin discs and will involve an extension of the asymptotic methods used in Paper I. Primary consideration is given to axisymmetric modes, but non-axisymmetric perturbations are also discussed brieﬂy. The continuous spectrum requires a separate analysis which can be made for a much more general equilibrium, and which is not restricted to axisymmetric modes. This is presented in Section 2. In Section 3 the equations and boundary conditions for axisymmetric modes in a weakly magnetized thin disc are derived. Some general properties of these equations are discussed in Section 4 and the numerical method of solution is described. In Section 5 some analytical results are obtained, principally for non-magnetized discs, which allow the modes to be enumerated and classified meaningfully in certain circumstances. For weakly magnetized discs, the most important result of this paper is the stability boundary in the parameter space of polytropic equilibria, which is obtained in Section 6. Non-axisymmetric modes are discussed brieﬂy in Section 7, and a concluding discussion is given in Section 8.

The energy principle of Bernstein et al. (1958) has played an important role in the analysis of the stability of magnetostatic equilibria relevant to astrophysics. Moss & Tayler (1969) examined the case of an axisymmetric, non-rotating star containing a poloidal magnetic field, and demonstrated that, if self-gravitation is neglected, the most unstable (or least stable) modes are those for which the azimuthal wavenumber m tends to infinity, and that the problem reduces to the consideration of a system of ordinary diﬀerential equations on each magnetic field line separately. Tayler (1973) considered the stability of a non-rotating star with a purely toroidal magnetic field, and found that the problem reduces to a system of algebraic equations at each separate point in the meridional plane. In each case, the most important displacements are those localized on a single magnetic field line or magnetic surface, and the equations obtained are closely related to those governing the continuous spectrum, which could be derived directly using the methods outlined in Section 2 below.

Throughout this paper, physical quantities are written in SI units with the permeability of free space, µ0, omitted for convenience. The coordinates used are cylindrical polar coordinates (r, φ, z) and the magnetic ﬂux coordinates (ψ, φ, χ) defined

c(cid:13) 1994 RAS, MNRAS 000, 000-000

(2.2)

(2.1)

δΠ = −(γp + 1

2 THE CONTINUOUS SPECTRUM

2 B2)∇·ξ − ξ·∇Π + B · (B·∇ξ)

Waves and instabilities in a diﬀerentially rotating disc

The linearized equation governing the Lagrangian displacement ξ(r, t) corresponding to a small departure from any state of ideal MHD may be written

D2ξ Dt2 = −∇δΠ − (∇·ξ)∇Π − ξ·∇∇Π − ρξ·∇∇Φ + B ·∇ [B·∇ξ − (∇·ξ)B] , where Π is the total pressure and

in Paper I. The magnetic ﬂux coordinates form a right-handed orthogonal coordinate system such that the magnetic field is B = ∇ψ × ∇φ = B(ψ, χ) eχ, where φ is the usual azimuthal angular coordinate. The Jacobian of the coordinate system is J = 1/|∇ψ||∇φ||∇χ|.

is its linearized Eulerian perturbation. Here ρ, p, B and γ are the density, pressure, magnetic field and adiabatic exponent, respectively, and D/Dt is the Lagrangian time derivative. This equation is equivalent to that given by Frieman & Rotenberg (1960), but is more general in that the gravitational potential Φ is included, and the equation holds for an arbitrary basic ﬂow. The self-gravitation of the ﬂuid is neglected here, as it was neglected throughout the construction of equilibria in Paper I. Although self-gravitation has, in general, a destabilizing inﬂuence on long-wavelength perturbations, its eﬀect on continuum modes, which are localized on a single magnetic surface, is nil. The equilibrium state under consideration is of the general class described in Section 2 of Paper I.

where ω ∈ C is the frequency eigenvalue, m ∈ ZZ the azimuthal wavenumber, and k ∈ IR a parameter which is allowed to tend to infinity. The function f is a unimodal ‘wavelet’1 of width w(k) centered on ψ = ψ0, and forms the envelope of the wave packet. In the limit k → ∞, the function f is to become infinitely localized at ψ = ψ0, but the width w(k) of the function should tend to zero more slowly than k−1, perhaps w ∝ k−1/2. Then the number of oscillations under the envelope tends to infinity, and a derivative of ξ with respect to ψ corresponds at leading order to multiplication by ik. The generalized function defined in the limit k → ∞ represents an infinitely localized wave packet. It is now possible to obtain the equations satisfied by the mode in this limit.

It is convenient to adopt a normalization such that |ξ| = O(1) at ψ = ψ0. If one is to obtain a solution of equation (2.1) with a frequency eigenvalue that has a finite limit (and a finite value of m), then it is clear that the ordering must be such that ∇·ξ = O(1) [rather than O(k)] and δΠ = O(k−1) [rather than O(k) or O(1)]. Therefore the fast magnetoacoustic wave is filtered out in this limit, as in the magneto-Boussinesq approximation (Spiegel & Weiss 1982). It also follows that the component ξψ is O(k−1), so the displacement is confined within the magnetic surface (at leading order). One may write

The continuum modes may be described using a technique which was introduced by Papaloizou & Pringle (1982) and developed by Lin, Papaloizou & Kley (1993) and Terquem & Papaloizou (1996). The modes take the form of wave packets which are infinitely localized in the ψ-direction so that they are eﬀectively confined to a single magnetic surface. Localization in the χ-direction cannot be achieved because of the infinite restoring force that would result from bending the magnetic field lines. Using the magnetic ﬂux coordinates (ψ, φ, χ) of Paper I, one considers a sequence of functions of the form

In order to write equation (2.1) in magnetic ﬂux coordinates, it is convenient to introduce the angle of inclination i(ψ, χ)

(cid:19) deφ = − cos i eψ dφ − sin i eχ dφ,

Then the derivatives of the unit vectors are given by

of the magnetic field to the vertical, such that

c(cid:13) 1994 RAS, MNRAS 000, 000-000

where ∆ and ̟ are both O(1).

1 For example, f (x) = exp(−x2).

eχ dψ + cos i eφ dφ −

eψ dψ + sin i eφ dφ +

cos i = er · eψ = rB

exp(ikψ + imφ − iωt)

sin i = er · eχ =

ψ − ψ0 w(k)

ξ(r, t) = Re

deψ = −

∇·ξ = ∆

δΠ = k

−1̟,

eψ dχ,

eχ dχ,

deχ =

∂i ∂ψ

∂r ∂ψ

∂r ∂χ

∂i ∂χ

ξ(χ)f

(2.3)

(2.6)

(2.7)

(2.8)

(2.4)

(2.5)

and

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:20)

(cid:21)

(cid:17)

(cid:16)

−

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

ξχ

− ρ

and

B J

∆ −

with

(2.9)

ξφ r

− rB

(2.11)

(2.10)

− B∆

− sin i

∂ ∂χ

∂i ∂χ

∂ ∂χ

∂i ∂ψ

∂Φ ∂χ

∂Φ ∂ψ

∂Φ ∂χ

= −rB

∂Π ∂χ

∂Π ∂ψ

∂Π ∂χ

ξχ JB

∂ξφ ∂χ

∂ξχ ∂χ

(cid:19) ∂i ∂χ

∂(JB) ∂ψ

0 = −(γp + B2)∆ −

4 G. I. Ogilvie

− ρˆω2ξχ + 2iρˆωΩ sin i ξφ − ρΩ2 sin2 i ξχ = −

−ρˆω2ξφ − 2iρˆωΩ sin i ξχ − ρΩ2ξφ = − cos i rB

in terms of ξφ, ξχ and ∆, but is otherwise unimportant. The remaining equations at leading order are

The quantity ̟ appears at leading order only in the ψ-component of equation (2.1), which therefore serves to define ̟

where ˆω = ω − mΩ is the intrinsic frequency. These are the exact equations defining the continuous spectrum. It should be understood that all quantities are to be evaluated on the magnetic surface ψ = ψ0; the equations are essentially ordinary diﬀerential equations with ψ treated as a parameter. Once ∆ has been eliminated, the remaining equations may be written

If, instead, the field line is closed, then periodic boundary conditions should be applied, with δB being continuous on S. It can be shown that δBχ vanishes automatically on S (at leading order), and no boundary condition on ξχ is obtained. The explanation is that equation (2.14) is singular at the surface, where vs vanishes, and a regularity condition applies to ξχ there. In the absence of rotation, equations (2.13) and (2.14) are two uncoupled, second-order diﬀerential equations in Sturm- Liouville form, describing the Alfv´en continuum and the cusp continuum, respectively, for the polarization for which the displacement is confined to the magnetic surface. These equations have been given by Poedts, Hermans & Goossens (1985), who also showed that the addition of a toroidal magnetic field results in a coupling between the equations. Similarly, in this case, the equations are coupled by rotation, if the magnetic field is not purely vertical. A second eﬀect of rotation is to modify the eﬀective gravitational acceleration gχ.

while equation (2.11) becomes vacuous. If the field line crosses the surface S of the disc at points χ = χ1 and χ = χ2, say, and extends to infinity, then the relevant solution of equation (2.17) is δBφ = 0, precisely as if the exterior were a vacuum. Since δB must be continuous on S, the boundary condition on ξφ is

Boundary conditions must be supplied for these equations. The poloidal magnetic field line may either form a closed loop or extend to infinity, and in general will have segments both inside and outside the disc. The exterior region is to be treated as a force-free medium of zero density and pressure, in which equation (2.10) reduces to

where vs = (γp/ρ)1/2 and vA = (B2/ρ)1/2 are the sound speed and the Alfv´en velocity, respectively,

is a quantity analogous to the square of the Brunt-V¨ais¨al¨a frequency for displacements within the magnetic surface.

is the eﬀective gravitational acceleration parallel to the magnetic field, and

c(cid:13) 1994 RAS, MNRAS 000, 000-000

ˆω2ξχ − 2iˆωΩ sin i ξφ = −

ˆω2ξφ + 2iˆωΩ sin i ξχ = −

at χ = χ1 and χ = χ2.

v2 A s + v2 v2

v2 s s + v2 v2

rΩ2 sin i

∂ ln ρ ∂χ

ρJB

JB2

(rδBφ) =

ρrJ

gχ = −

χ = gχ

∂Φ ∂χ

B2 J

∂ ∂χ

(2.13)

(2.14)

(2.15)

(2.16)

(2.17)

(2.18)

(2.12)

gχ v2 s

r2 J

ξχ B

ξφ r

A (cid:19)

Bgχ

= 0,

(cid:26)(cid:18)

(cid:17)(cid:21)

χ −

and

N 2

= 0

(cid:21)(cid:27)

(cid:20)(cid:18)

ξχ,

−

(cid:16)

(cid:17)

(cid:20)

(cid:21)

(cid:20)

(cid:19)

(cid:18)

r2 J

∂ ∂χ

∂un ∂χ

nρr2Jun = 0,

Waves and instabilities in a diﬀerentially rotating disc

subject to the boundary condition ∂un/∂χ = 0 at χ = χ1 and χ = χ2 (or periodic boundary conditions for a closed field line), and the normalization condition

It is clear that the azimuthal wavenumber m appears in equations (2.13), (2.14) and (2.18) only in the combination ˆω = ω − mΩ. Moreover, ω, Ω and ˆω are all constant on the magnetic surface. It follows that the eigenfunctions are the same for all values of m, while the frequency eigenvalues of non-axisymmetric modes are related to those of axisymmetric modes simply by a Doppler shift. In the same way that Papaloizou & Szuszkiewicz (1992) derived a variational principle for axisymmetric modes, so it is possible to derive a variational principle for modes of arbitrary m, provided that attention is n(ψ) : n ∈ ZZ+} and {un(ψ, χ) : n ∈ ZZ+} to be the restricted to the continuous spectrum. One may proceed by defining {ω2 ordered, real eigenvalues and normalized, real eigenfunctions of the ‘Alfv´en operator’ that appears in equation (2.13). This is a Sturm-Liouville problem,

It should be noted, however, that there is a significant diﬀerence in the continuous spectrum between the case of no rotation (or uniform rotation) and that of diﬀerential rotation. In the former case, one may consider functions of the form (2.3) and take the simultaneous limits k → ∞ and m → ∞. A continuous range of polarizations is obtained by choosing diﬀerent limiting values of the ratio k/m, and the corresponding eigenfunctions need not have displacements confined within the magnetic surface. This leads naturally to the conclusions of Moss & Tayler (1969). The equations given by Poedts et al. (1985) in the case of a purely poloidal magnetic field describe only that part of the continuous spectrum with the polarization for which k/m → ∞, since m is assumed finite. When the ﬂuid is in diﬀerential rotation, however, the limit m → ∞ cannot be taken, because an infinite diﬀerential Doppler shift would result. The presence of a toroidal magnetic field would also prevent the limit m → ∞ from being taken.

and exists provided that the appropriate coeﬃcient an vanishes should ˆω2 be equal to any of the eigenvalues of the Sturm- Liouville problem. The spectrum of the Alfv´en operator has the following characteristics. The eigenvalues {ω2 n(ψ)} for each magnetic surface form an ordered, denumerably infinite sequence of distinct, non-negative real numbers with limit +∞. The lowest eigenvalue is always ω2 0(ψ) = 0, the corresponding eigenfunction u0(ψ, χ) being independent of χ. This zero-frequency mode is a trivial displacement in which the entire magnetic surface is rotated through an infinitesimal angle about the z-axis.2 A standard method (Courant & Hilbert 1953) leads to the asymptotic expression

Z for each value of ψ. Let {an} and {bn} be the components of (sin i ξχ/r) and (ξφ/r) with respect to this set of eigenfunctions, such that

2 In the absence of a magnetic field, the infinitesimal rotation of any ring of ﬂuid about the axis is a trivial displacement, but these trivial displacements are easily eliminated from the analysis by discarding the root ˆω = 0.

The eigenfunction expansions (2.21) may be substituted into equation (2.14), the equation multiplied by ρJξ∗

Then the solution of the inhomogeneous equation (2.13) is given by

for the eigenvalues in the limit n → ∞, where

is the Alfv´en time of the magnetic field line.

integrated with respect to χ to yield

c(cid:13) 1994 RAS, MNRAS 000, 000-000

ρr2Jumun dχ = δmn,

L[ξχ; ˆω2] = R[ξχ],

n2π2 τ 2 A(ψ)

ρ1/2J dχ =

hχ dχ vA

bnun(ψ, χ).

anun(ψ, χ)

ˆω2 − ω2

τA(ψ) =

2iˆωΩ(ψ)

n(ψ) =

bn = −

O(1)

χ and

(2.19)

(2.20)

(2.21)

(2.22)

(2.25)

(2.23)

(2.24)

n=0 X

ξχ r

ξφ r

n(ψ)

sin i

and

an,

ω2

χ1

χ2

χ1

χ2

∞

χ1

χ2

χ1

χ2

χ1

(cid:16)

(cid:12) (cid:12) (cid:12) (cid:12)

(cid:20)(cid:18)

(cid:21)(cid:27)

N 2

and

χ −

Bgχ

A (cid:19)

dχ +

ξχ B

n=1 X

(2.28)

(2.27)

(2.26)

|an|2

(ˆω2)i

∂ ∂χ

χ1 (cid:18)

B2 J

(cid:17)(cid:12) (cid:12) (cid:12) (cid:12)

χ1 (cid:26)(cid:18)

JB2

ρ|ξχ|2 J dχ

ρ|ξχ|2 J dχ +

ρ|ξχ|2 J dχ −

L[ξχ; ˆω2] = ˆω2

v2 A v2 s + v2

v2 s v2 s + v2

4ˆω2Ω2 ˆω2 − ω2 n

4ω2 nΩ2 |ˆω2 − ω2

nΩ2 4ω2 (ˆω2 − ω2

ρ|ξχ|2 J dχ + (ˆω2)i

R[ξχ] = 4Ω2|a0|2 +

6 G. I. Ogilvie

n|2 |an|2 = 0,

n)2 |an|2 > 0,

where the two functionals

∂ ∂ ˆω2 L[ξχ; ˆω2] =

have been introduced. Note that the infinite sum on the left-hand side starts at n = 1, the term n = 0, which is independent of ˆω2, having been transferred to the right-hand side. If ˆω2 = (ˆω2)r + i(ˆω2)i, then the imaginary part of equation (2.25) is

from which it follows that the eigenvalues ˆω2 are real. The eigenfunctions ξχ may also be taken to be real. Now consider L[ξχ; ˆω2] as a function of ˆω2 for a given trial function ξχ. Its derivative with respect to ˆω2 is

and so the function is monotonic wherever the derivative is defined. The function has a pole, however, at each of the values ˆω2 = ω2 n (other than n = 0) for which the coeﬃcient an does not vanish. The function therefore increases once through all real values in each of the intervals (−∞, p1), (p1, p2), (p2, p3), … , where p1, p2, p3, … are the abscissae of the poles. Equation (2.25) may therefore be understood as defining a multiple-valued functional ˆω2[ξχ]. Moreover, the stationary values of this functional may be demonstrated by standard techniques to be identical to the true eigenvalues ˆω2. Corresponding to each eigenfunction ξχ is one eigenvalue in each of the intervals (−∞, p1), (p1, p2), (p2, p3), … .3 Now L[ξχ; 0] = 0, which is the reason for transferring the term n = 0 to R[ξχ]. A necessary condition for the existence of a negative eigenvalue ˆω2 is that the functional R[ξχ] admit negative values for trial functions satisfying the boundary conditions. The variational principle ensures that this condition is also suﬃcient. Therefore a necessary and suﬃcient condition for instability in the continuous spectrum is that the functional R[ξχ] admit negative values for trial functions satisfying the boundary conditions. In particular, stability is assured if v2 s s + v2 v2

The principal aim of this paper is to determine the spectrum of waves and instabilities in a magnetized thin disc of the type described in Paper I. Although asymptotic analysis in the parameter ǫ underlies the mathematics, calculations are made to leading order only, and therefore few explicit references to ǫ are required. Of the two families of thin discs, the weakly magnetized discs have the more interesting spectrum, for they are potentially unstable to the magnetorotational instability.

Significantly, the confinement of the displacements within the magnetic surface implies that no instabilities associated with diﬀerential rotation are found. In particular, the magnetorotational instability does not appear in the continuous spec- trum. This is in contrast with the case of a purely toroidal magnetic field (Terquem & Papaloizou 1996). Neither is a magnetoconvective instability associated with the ψ-component of gravity found.

The only type of instability that can be present in the continuous spectrum is a magnetoconvective (Parker) instability modified by rotation. The instability is due to the release of gravitational potential energy by displacements parallel to the magnetic field. The first term in R[ξχ], which may be expressed in the form

represents the stabilizing eﬀect of rotation, but only on modes for which a0 6= 0. In particular, if the disc is symmetrical about the equatorial plane, then this stabilizing eﬀect is zero for all modes with odd symmetry.

3 Although these modes share the same displacement ξχ, they have diﬀerent displacements ξφ, according to the relation (2.22).

(cid:18) throughout χ1 < χ < χ2, for each magnetic surface.

3 AXISYMMETRIC WAVES AND INSTABILITIES IN THIN DISCS

c(cid:13) 1994 RAS, MNRAS 000, 000-000

4Ω2 sin2 i ρ|ξχ|2 J dχ,

v2 A s + v2 v2

4Ω2|a0|2 =

JB2

∂ ∂χ

|an|2

(2.29)

(2.30)

(2.31)

! Z

|a0|2

n=0 X

A (cid:19)

Bgχ

χ −

N 2

> 0

(cid:20)(cid:18)

(cid:30)

(cid:21)

χ2

χ1

∞

1/2

Ω0 =

GM r3

equations for a weakly magnetized disc are

3.1 Review of weakly magnetized discs

Waves and instabilities in a diﬀerentially rotating disc

It is assumed that the pressure and density satisfy a polytropic relation on each magnetic surface. Then the leading-order

Indeed, a major objective is to determine the stability boundary in the parameter space of weakly magnetized, polytropic discs. Strongly magnetized discs are stable to the magnetorotational instability and so are not considered in detail.

The analysis of thin discs in Paper I is based on a small parameter ǫ which may be defined as either the maximum value, or a characteristic value, of H(r)/r, where z = H(r) is the location of the upper surface of the disc at radius r. The internal structure of the disc is resolved by introducing a stretched vertical coordinate ζ = ǫ−1z whose value at the surface of the disc is ζs = ǫ−1H. There are two aspects to obtaining a solution at leading order. First, there are ordinary diﬀerential equations in ζ that determine the vertical equilibrium of the disc at each radius. Secondly, there is an integral relation that determines the global magnetic structure; however, this is not directly relevant to the analysis that follows.

The analysis presented here may be seen as an extension of the work of Lubow & Pringle (1993; hereafter, LP) and Korycansky & Pringle (1995; hereafter, KP) on the equivalent problem for hydrodynamic thin discs. Previously, Ruden, Papaloizou & Lin (1988) had examined convectively unstable modes in a thin disc using a similar method. This analysis also extends the work of Gammie & Balbus (1994), who examined the magnetorotational instability of a stratified disc containing a uniform magnetic field. The present analysis is concerned with a polytropic disc containing a poloidal magnetic field that bends, in general, as it passes through the disc.

The asymptotic method used to construct the equilibria can be applied to study their normal modes. For a weakly magnetized disc, the angular velocity, the buoyancy frequency, and the characteristic frequencies of acoustic and Alfv´en waves with wavelengths comparable to H are all O(1), which implies that this is the appropriate scaling for the frequency eigenvalue of a normal mode. Only if a mode is found with eigenvalue 0 or ∞ in this scaling need a diﬀerent scaling be considered. In

where Br0s(r) denotes the value of Br0 on the upper surface, and is supposed to be known from the global magnetic structure. The solution in −ζs < ζ < 0 is inferred from symmetry, ρ0, p0 and Ω1 being even functions of ζ, while Br0 is odd.

Here s is a free parameter which should be positive if self-gravitation is to be unimportant.

These equations are to be solved on 0 < ζ < ζs(r), with boundary conditions

where the angular velocity, density, pressure and magnetic field are given by

B(r, z) = ǫs+1 [Br0(r, ζ) + O(ǫ)] er + ǫs+1 [Bz0(r) + O(ǫ)] ez.

(cid:3) p0(r, ζ) + ǫp1(r, ζ) + O(ǫ2)

Ω(r, z) = Ω0(r) + ǫΩ1(r, ζ) + O(ǫ2),

ρ0(r, ζ) + ǫρ1(r, ζ) + O(ǫ2)

c(cid:13) 1994 RAS, MNRAS 000, 000-000

2ρ0rΩ0Ω1Br0 Bz0

(cid:16) (cid:17) = −ρ0Ω2

(cid:2) p(r, z) = ǫ2s+2

Br0(r, ζs(r)) = Br0s(r),

3.2 Normal modes

2ρ0rΩ0Ω1 Bz0

3Ω0Br0 2rBz0

ρ0(r, ζs(r)) = 0

p0 = K0ρΓ 0 ,

ρ(r, z) = ǫ2s

Br0(r, 0) = 0,

∂Br0 ∂ζ

∂Ω1 ∂ζ

∂p0 ∂ζ

(3.10)

0ζ +

(3.1)

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

(3.8)

(3.9)

= −

and

(cid:2)

(cid:3)

(cid:20)

−1

(cid:26)

(cid:21)(cid:27)

(3.12)

(3.11)

k(r) dr

−iωt + iǫ

ξ(r, t) ∼ Re

ξ0(r, ζ) exp

8 G. I. Ogilvie

ξ(r, t) ∼ Re [ξ0(r, ζ) exp(−iωt)] .

When the expression (3.12) is substituted into equation (2.1), the following diﬀerential system is obtained at leading

deciding the form of the eigenfunctions of axisymmetric normal modes at leading order, it is natural to assume that they take the form

where k(r) is a radial wavenumber. This is, of course, a WKB function, but it should be emphasized that no additional approximation or limit is involved here; this form arises simply from the fact that the eigenfunction varies on a radial length scale that is much shorter than that of the equilibrium disc.

If this is substituted into equation (2.1), there results at leading order a sixth-order diﬀerential system4 on −ζs < ζ < ζs at each radius separately, which, when supplemented by appropriate boundary conditions, would constitute an eigenvalue problem for ω. The eigenvalues would be discrete, but would depend on r in general, which means that a global solution, connecting a range of radii, could not be obtained other than in exceptional circumstances. One is led to conclude that the eigenfunction must depend more strongly on r. The correct form of the eigenfunction is

Strictly speaking, equation (3.12) represents a wave propagating in the r-direction (when ω is real), whereas a normal mode is formed from the superposition of two such waves propagating in opposite directions. However, the local dispersion relation derived below does not depend on the sign of k, and so it is suﬃcient to consider travelling waves of this form. When ω is imaginary, equation (3.12) represents an exponentially growing or decaying normal mode rather than a travelling wave. The scaling of ξ with ǫ is arbitrary, since this is a linear problem, but is chosen to be O(1) for convenience. For simplicity of notation, the subscript ‘0’ on all quantities is omitted hereafter.

This system is of sixth order and requires six boundary conditions (apart from the arbitrary normalization condition that applies to all linear eigenvalue problems). There are infinitely many discrete eigenvalues ω for each value of k, and these form a local dispersion relation ω(k; r) for the disc which has infinitely many branches. The result of Papaloizou & Szuszkiewicz (1992) that ω2 must be real for any global normal mode is reﬂected in the fact, demonstrated below, that the eigenvalues ω of the system (3.13)-(3.15) are either purely real or purely imaginary when k is itself real. As in a standard WKB problem, a global solution here is trapped in a ‘wave region’ between two points, each of which may be either a turning point (at which k goes to zero) or a boundary point (either the inner or outer radius of the disc), and within which k is real and non-zero. The variation of ξ0 with r is determined by higher-order terms which could be calculated in principle.

The boundary conditions for equations (3.13)-(3.15) are complicated because there is a free surface with an inclined magnetic field. In the ‘corona’ above the disc, where ζ > ζs and ρ = 0, the magnetic field is uniform on the spatial scale H and equations (3.13)-(3.15) reduce to the two equations

The equations for an incompressible ﬂuid (γ → ∞) are obtained by imposing the constraint ∆ = 0 on ξ, and making δΠ

a distinct variable in its own right, rather than using the indeterminate equation (3.17).

ω2 + Ω2 ∂ ln ρ ∂ ln ζ (cid:18) where

δΠ = −(γp + B2)∆ + ρΩ2ζξz + BrDξr + BzDξz

4 Identical to equations (3.13)-(3.15) below, but with k = 0.

ρ(ω2ξφ + 2iωΩξr) = −D2ξφ,

− ρΩ2ζ∆ − D(Dξz − Bz∆),

∂2ξψ ∂ζ 2 = k2B2

c(cid:13) 1994 RAS, MNRAS 000, 000-000

(ω2 + 3Ω2)ξr − 2iωΩξφ

= ik δΠ − D(Dξr − Br∆),

D = ikBr + Bz

∆ = ikξr +

∂ δΠ ∂ζ

∂ξz ∂ζ

order:

(3.13)

(3.14)

(3.15)

(3.16)

(3.17)

(3.18)

(3.19)

∂ ∂ζ

ξz =

z ξψ

and

B2 z

(cid:19)

(cid:2)

(cid:3)

(cid:19)

(cid:18)

(cid:19)

and

∂ ∂ζ

− Brs

(3.23)

(3.22)

(3.21)

(3.20)

δBφ =

ξφ = 0,

ξφ = 0.

∂ξr ∂ζ

ikBrs + Bz

ξψ ∝ exp (−|k|ζ) ,

= −|k|(Bzξr − Brsξz)

ξφ ∝ exp [−i(Brs/Bz)kζ] ,

since this corresponds to a vanishing Eulerian perturbation

when k is real. The appropriate solution of equation (3.20) is

Waves and instabilities in a diﬀerentially rotating disc

(cid:18) where ξψ is the component of ξ in the meridional plane perpendicular to the magnetic field. The appropriate solution of equation (3.19) is

In each case the solution is such that the components of the Eulerian perturbation of the magnetic field go to zero at infinity. The corresponding boundary conditions for the interior solution at ζ = ζs are easily obtained by applying the continuity of ξ and its first derivatives. Thus ∂ξz ∂ζ

applies, where r1 and r2 are the limits of the wave region, n is an integer and δ is a phase constant. This implies that n = O(ǫ−1) for all modes with the scalings considered, and therefore the fact that it is an integer is irrelevant at this level of approximation. If at any radius the local dispersion relation has a root ω corresponding to some real value of k, then a global normal mode must exist with that frequency eigenvalue at leading order. The spectrum is dense in the asymptotic limit under consideration. This also means that a necessary and suﬃcient condition for the existence of an unstable global normal mode is that the local dispersion relation should have an imaginary root ω for some real value of k. This is true, at least, for axisymmetric modes with dynamical [i.e. O(1)] growth rates.

The three remaining boundary conditions are the equivalents of the first three at the lower surface ζ = −ζs. Alternatively, these may be replaced by symmetry conditions at ζ = 0, since the solutions are either even or odd. It is convenient to define an ‘even’ mode as one which preserves the symmetry of the equilibrium about the equatorial plane. For such a mode, the Eulerian perturbation of any quantity must have the same symmetry as that quantity has in the equilibrium state. This implies that ξr and ξφ are even functions of ζ, while ξz is odd. An ‘odd’ mode is one that breaks the symmetry of the equilibrium, implying that ξr and ξφ are odd, while ξz is even. The boundary conditions at ζ = 0 are therefore

For completeness, the equations for axisymmetric waves in a strongly magnetized thin disc should be mentioned. In such a disc, the angular velocity and the buoyancy frequency are both O(1), but the characteristic frequencies of acoustic and Alfv´en waves with wavelengths comparable to H are both O(ǫ−1/2). This implies that an appropriate scaling for the frequency eigenvalue of a normal mode is O(ǫ−1/2). A set of equations is then obtained, equivalent to equations (3.13)-(3.15) but with Ω set to zero. The modes are either Alfv´en waves with ξr = ξz = 0 or magnetoacoustic waves with ξφ = 0; neither rotation nor buoyancy has any eﬀect, and there is no instability. However, there is a mode with exactly zero frequency in this scaling,

The third boundary condition at ζ = ζs is a regularity condition, which is necessary because the sound speed goes to zero at the surface. The equation for ξχ, the component of ξ parallel to the magnetic field, has a singular point there and only the regular solution is acceptable.

The relation between the local dispersion relation and the spectrum of global normal modes is quite straightforward. As

for an odd mode. Note that this convention is opposite to that used by LP and KP.

in any WKB problem, a quantization condition of the form

c(cid:13) 1994 RAS, MNRAS 000, 000-000

for an even mode, and

k(r) dr = nπ + δ

= −ikBrsξφ.

ξr = ξφ =

= ξz = 0

∂ξφ ∂ζ

∂ξz ∂ζ

∂ξr ∂ζ

(3.24)

(3.25)

(3.26)

(3.27)

(3.28)

−1 ǫ

= 0

4.1 Overview

10 G. I. Ogilvie

4.2 A variational principle

4 PROPERTIES OF THE LOCAL DISPERSION RELATION

having an eigenfunction corresponding to a uniform displacement. This allows a class of modes to exist that have frequency eigenvalues ω = O(1) and vary on a long radial length scale, as in equation (3.11). These are the bending modes, which can be described using a two-dimensional treatment, and which can be unstable if the magnetic field provides a significant amount of support against gravity (Agapitou et al. 1997).

Before deriving the variational principle it may be noted that the dispersion relation ω(k) has reﬂectional symmetry in both ω- and k-axes. This can be seen from the fact that equations (3.13)-(3.15) have the symmetries (ξφ, ω) 7→ (−ξφ, −ω) and either (ξ, ω, k) 7→ (ξ∗, −ω, −k), if ω is real, or (ξ, ω, k) 7→ (ξ∗, ω, −k), if ω is imaginary. When presenting the numerical results it is therefore suﬃcient to restrict attention to positive values of ω (or ω/i) and of k.

The sixth-order eigenvalue problem defined by equations (3.13)-(3.15) and the boundary conditions is considerably more complicated than the hydrodynamic equivalent. An attempt is made in Section 5 below to classify the hydrodynamic modes, but little progress can be made when a magnetic field is included. It is useful, however, to derive a variational principle for the frequency eigenvalues. This not only shows that ω2 must be real when k is real, but also helps to clarify under what circumstances unstable modes are expected. The method here is an adaptation of the analysis of Papaloizou & Szuszkiewicz (1992) and closely follows the method of Section 2.

Z for each r. This is entirely analogous to the situation in Section 2. Physically, the eigenfunctions represent Alfv´en waves in a fictitious, non-rotating disc with the same density and vertical magnetic field, but without a radial magnetic field. The eigenvalues {ω2

n(r) : n ∈ ZZ+} and the real eigenfunctions {un(r, ζ) : n ∈ ZZ+} of the Sturm-Liouville equation ∂2un ∂ζ 2 + ω2

To invert equation (3.14), which is subject to the boundary condition (3.25), one may define the ordered, real eigenvalues {ω2

is related to the Sturm-Liouville equation by a unitary transformation. It has the same eigenvalues and its eigenfunctions are simply

(cid:20) Equation (3.14) may be written in the form

and is conveniently solved using the eigenfunction expansions

n} are the corresponding squared frequency eigenvalues.

subject to the boundary conditions

D2ξφ + ω2ρξφ = −2iωΩρξr,

and the normalization condition

c(cid:13) 1994 RAS, MNRAS 000, 000-000

˜un(r, ζ) = un(r, ζ) exp

2Ω1(r, ζ) 3Ω(r)

an(r)un(r, ζ) exp

bn(r)un(r, ζ) exp

ρumun dζ = δmn

= un(r, ζ) exp

D2 ˜un + ω2

The equation

nρun = 0,

nρ˜un = 0,

ξφ(r, ζ) =

ξr(r, ζ) =

D˜un = 0

∂un ∂ζ

ζ = ±ζs

Br Bz

Z (cid:16)

n=0 X

(4.1)

(4.2)

(4.3)

(4.4)

(4.8)

(4.9)

(4.5)

(4.6)

(4.7)

and

B2 z

−ik

= 0

ikr

−ζs

dζ

−

(cid:17)

(cid:21)

(cid:20)

(cid:21)

(cid:20)

(cid:21)

(cid:20)

∞

ζs

(cid:3)

(cid:2)

∞

(cid:1)

(cid:0)

(cid:19)

−ζs

dζ,

and

n(r)

dζ −

n=1 X

ζs −ζs

(4.12)

(4.11)

(4.10)

|an|2

∂ρ ∂ζ

an(r),

−ζs (cid:20)

2iωΩ(r)

ω2 − ω2

bn(r) = −

(cid:18) |ξz|2

|ξr|2 + |ξz|2

The solution is

γp + B2

4ω2Ω2 ω2 − ω2 n

|k||Bzξr − Brξz|2

−3ρΩ2|ξr|2 − Ω2ζ

L[ξr, ξz; ω2, k] = ω2

involving the two functionals

R[ξr, ξz; k] = 4Ω2|a0|2 +

L[ξr, ξz; ω2, k] = R[ξr, ξz; k]

|BrDξr + BzDξz + ρΩ2ζξz|2

|δΠ|2 γp + B2 + |Dξr|2 + |Dξz|2 −

Waves and instabilities in a diﬀerentially rotating disc

where it has been assumed that k is real. The imaginary part of equation (4.11) is

and exists at any particular value of r provided that the appropriate coeﬃcient an(r) vanishes should ω2 be equal to any of the eigenvalues {ω2 n(r)}. One can now substitute for ξφ in equation (3.13) and form the integral relation [cf. equation (2.25)]

(cid:19) If both the specific entropy and the (radial) magnetic field strength are non-decreasing functions of ζ for 0 < ζ < ζs, then this coeﬃcient is non-negative and it follows that only the term −3ρΩ2|ξr|2 in the integrand of equation (4.15) can lead to instability. This implies that, in a disc with a convectively stable stratification, any instability is due to the diﬀerential rotation.

where ω2 = (ω2)r + i(ω2)i, and it follows that ω2 is real. Moreover, it can be shown that equation (4.11) has a variational property which implies that the necessary and suﬃcient condition for the existence of an unstable mode at any radius is that there exists a trial displacement ξ, satisfying the boundary conditions, which makes the functional R[ξr, ξz; k] negative for some real value of k.

etc., in which case the trial displacement must be defined on −∞ < ζ < ∞ and go to zero as |ζ| → ∞, but need not satisfy any boundary conditions at ζ = ±ζs.

R[ξr, ξz; k] = 4Ω2|a0|2 ∞ |δΠ|2 γp + B2 + |Dξr|2 + |Dξz|2 −

The functional R[ξr, ξz; k] can be expressed a number of diﬀerent ways which are useful for some purposes. One alternative

BrDξr + BzDξz − (cid:12) (cid:12) ∂ρ (cid:12) Ω2ζ (cid:12) ∂ζ

A second useful variation is to extend the integral in R[ξr, ξz; k] over the interval −∞ < ζ < ∞ and write

(cid:20) The coeﬃcient of |ξz|2 in the integrand is

B2 |BzDξr − BrDξz|2

(cid:19) −3ρΩ2|ξr|2 −

|BrDξr + BzDξz + ρΩ2ζξz|2

c(cid:13) 1994 RAS, MNRAS 000, 000-000

R[ξr, ξz; k] =4Ω2|a0|2 +

|δΠ|2 (cid:2) γp + B2 +

n|2 |an|2 = 0,

−3ρΩ2|ξr|2 − Ω2ζ

|k||Bzξr − Brξz|2

nΩ2 4ω2 |ω2 − ω2

(ρΩ2ζ)2 γp

γp γp + B2

γp + B2

(cid:18) |ξz|2

|ξr|2 + |ξz|2

Ω2ζ (cid:20)

dζ + (ω2)i

∂ ln ρ ∂ζ

∂ ln p ∂ζ

is to write

∂Br ∂ζ

= ρΩ2ζ

ρΩ2ζξz

−ζs (cid:26)

B2 γp

−∞ (cid:20)

Br γp

(ω2)i

∂ρ ∂ζ

(4.13)

(4.14)

(4.15)

(4.17)

(4.16)

ζs −ζs

n=1 X

dζ.

dζ,

(cid:12) (cid:12) (cid:12) (cid:12)

−ζs

(cid:27)

(cid:18)

(cid:19)

(cid:18)

(cid:19)

−

(cid:21)

(cid:0)

(cid:1)

∞

(cid:3)

ζs

(cid:2)

(cid:19)

−ζs

dζ,

(ω2

ρξ∗

ζs −ζs

|ξz|2

(4.19)

(4.18)

∂ρ ∂ζ

−ζs (cid:18) Z

1 − ω2 2)

1 · ξ2 dζ = 0,

|k||Bzξr − Brξz|2

12 G. I. Ogilvie

4.3 Numerical method

R[ξr, ξz; k] = 4Ω2|a0|2 +

|Dξr|2 + |Dξz|2 − 3ρΩ2|ξr|2 − Ω2ζ

Finally, the version for an incompressible ﬂuid is

(cid:3) where ξ is subject to the constraint ∆ = 0.

Despite the unconventional location of the frequency eigenvalue in the variational principle, the orthogonality principle is straightforward. If ξ1 and ξ2 are two eigenfunctions at the same radius, and for the same real value of k, which correspond to frequency eigenvalues ω1 and ω2, then

which implies that the eigenfunctions are orthogonal provided that the squared eigenvalues are distinct. Note that this orthogonality principle involves all three components of ξ, not just the meridional part. It is conjectured that the eigenfunctions at each radius, and for each real value of k, form a complete set in the space of continuous, vector-valued functions of ζ on −ζs < ζ < ζs.

The sixth-order eigenvalue problem (3.13)-(3.15) is solved numerically by the following method. First, the equilibrium itself must be computed, as in Paper I. When Brs and Bz are specified, equations (3.1)-(3.5) constitute a third-order, non-linear eigenvalue problem for Ω1s (the value of Ω1 at ζ = ζs). This is solved using the shooting method: a value of Ω1s is guessed and the equations are integrated from ζ = ζs to ζ = 0. The amount by which the boundary condition Br0(r, 0) = 0 fails to be satisfied defines a mismatch function whose derivative with respect to Ω1s is obtained by simultaneously integrating the equations diﬀerentiated with respect to this parameter. Newton-Raphson convergence is then applied to obtain the eigenvalue. A similar technique is applied to equations (3.13)-(3.15). However, the coeﬃcients in these equations are known only numerically, and so to avoid the use of interpolation one must integrate equations (3.1)-(3.5) simultaneously (with the previously determined eigenvalue). When k is specified, equations (3.13)-(3.15) contain only one parameter, ω. However, since only three boundary conditions apply at ζ = ζs, one must also guess two boundary values (the sixth being given by a normalization condition). The full problem requires shooting in a three-dimensional complex space. However, Newton-Raphson iteration can still be applied successfully.

where the numerators Nr and Nz are regular and are linear functions of ξr, ξφ, ξz and their first derivatives. The condition that Nr and Nz both go to zero as fast as γp as ζ → ζs defines a regularity condition. However, to start the numerical integration at ζ = ζs, these quotients must be evaluated. This is done by expanding all quantities in power series (not Taylor series) about the singular point, a procedure that involves much tedious algebra.

Once a mode has been computed for some value of k, it can be followed quasi-continuously as k is varied, so tracing out a single branch of the dispersion relation. This is then repeated for all branches within some limited region of the dispersion diagram. Fortunately there are certain limits in which the modes can be enumerated using analytical methods.

Similarly, frequencies and growth rates are expressed in units of the local Keplerian angular velocity, and the radial wavenumber in units of H −1.

All numerical calculations are carried out in units such that GM = K0 = r = 1, and ǫ is defined such that ζs = 1 at the

The most diﬃcult part of the numerical method is the treatment of the singular point at ζ = ζs. One can re-write

Taking into account the factors of ǫ, this means that the true unit of magnetic field strength is

value of r under consideration. In particular, this means that B0 is expressed in units of

equations (3.13) and (3.15) in the form

c(cid:13) 1994 RAS, MNRAS 000, 000-000

−3Γ/2(Γ−1)H Γ/(Γ−1).

−3Γ/2(Γ−1)ζ Γ/(Γ−1)

(GM )Γ/2(Γ−1)K

∂2ξr ∂ζ 2 =

∂2ξz ∂ζ 2 =

−1/2(Γ−1)

−1/2(Γ−1)r

Nz γp

Nr γp

(4.20)

(4.22)

(4.23)

(4.21)

(cid:17)

h(cid:16)

(5.3)

(5.2)

(5.1)

ρ0 =

0(ζ 2

1/(Γ−1)

2 Ω2

Γ/(Γ−1)

h0 K0

h0 = 1

p0 = K0

s − ζ 2).

Γ − 1 Γ

(cid:17) Γ − 1 Γ

5.1 Hydrodynamic discs revisited

5 CLASSIFICATION OF MODES

Waves and instabilities in a diﬀerentially rotating disc

i h0 K0 h(cid:16) where the enthalpy is

The spectrum of axisymmetric waves in a polytropic, hydrodynamic thin disc has been described by KP. When the magnetic field is absent, the solution for the equilibrium is

Equations (3.13)-(3.15) then become equivalent to those solved by KP in terms of Eulerian perturbations. In their analysis, which was for a convectively stable disc, they identified three distinct classes of modes: p modes (acoustic modes, due to compressibility and driven mainly by pressure forces) and g modes (gravity modes, due to buoyancy and driven mainly by gravitational forces), both with ω2 ≥ Ω2, and also r modes5 (inertial modes, due to rotation and driven mainly by inertial forces), with ω2 ≤ Ω2. Examples of the eigenfunctions for these modes can be found in KP. They found no modes analogous to the f mode (the fundamental, surface gravity mode) of stellar oscillations, but reported certain peculiarities concerning the ordering of the eigenvalues. The aim of this subsection is to clarify the situation and demonstrate that there are in fact two f modes.

The f modes are the only modes with ω2 ≥ Ω2 to survive both limits (ii) and (iii). There are two of them, one of each parity. The classification of the modes is particularly clear in the limit kζs → ∞. All modes with ω2 > Ω2 are then trapped near the surfaces of the disc. As explained by KP, the loss of contact between the two surfaces implies that the frequency eigenvalues of adjacent even and odd modes must coalesce so that a single surface wave of neither symmetry can exist. The WKB method used by KP gives the asymptotic form of the dispersion relation for large vertical mode number, but is not appropriate for studying modes of small vertical mode number such as the f modes. Instead, one can use an asymptotic method based on the large parameter kζs. It is found numerically that all modes with ω2 > Ω2 are trapped in a layer of extent O near each surface of the disc, and have ω2 = O(kζs) in this limit. This is consistent with the interpretation of (i) f modes as surface gravity modes (‘ω2 ∼ gk’), since g = O(1) in the layer; (ii) p modes as acoustic modes (‘ω2 ∼ v2 (iii) g modes as buoyancy modes (‘ω2 ∼ N 2’), since N 2 = O(kζs) in the layer. The scalings imply that neither rotation nor the non-uniformity of gravity is significant to these modes in the limit kζs → ∞, and the modes therefore obey the same equations at leading order as in a static, polytropic atmosphere with uniform gravity. This problem has been solved by Lamb (1932) in terms of conﬂuent hypergeometric functions. It was noted by Christensen- Dalsgaard (1980) that Lamb’s analysis is valid in an asymptotic sense for modes of large degree ℓ in an arbitrary stellar model. The same is true of modes of large kζs in an accretion disc. To prove this, one should introduce a stretched coordinate

An example of the local dispersion relation for a hydrodynamic disc is shown in Fig. 1. The disc is convectively stable, with Γ = 4/3 and γ = 5/3, and has a spectrum that is qualitatively similar to that calculated by KP using slightly diﬀerent parameters. The diﬀerent classes of modes may be distinguished (in a convectively stable disc) by the following characteristics: (i) f, p and g modes have ω2 ≥ Ω2, while r modes have ω2 ≤ Ω2; (ii) in the limit γ → Γ, in which the disc becomes adiabatically stratified, the g modes all collapse to ω2 = Ω2, while all other

(iii) in the limit γ → ∞, in which the ﬂuid becomes incompressible, the p modes have ω2 → ∞, while all other modes have

which is O(1) in the layer in which the modes are trapped. Equations (3.13)-(3.15) are then solved asymptotically with

c(cid:13) 1994 RAS, MNRAS 000, 000-000

modes have non-trivial limits;

ω2 = λΩ2kζs + O(1),

ξz(r, ζ) = ξz0(r, x1) + O

ξr(r, ζ) = ξr0(r, x1) + O

5 Called g modes by LP.

s k2’), since v2

non-trivial limits.

x1 = k(ζs − ζ)

in the layer;

(kζs)−1

s = O

(kζs)

(5.6)

(5.7)

(5.4)

(5.5)

and

−1

(cid:0)

(cid:1)

(cid:0)

(cid:1)

(cid:0)

(cid:1)

(cid:0)

14 G. I. Ogilvie

Figure 1. Left: local dispersion diagram for a polytropic, hydrodynamic thin disc. Right: an expanded view of the lower-frequency branches. The disc is convectively stable, with Γ = 4/3 and γ = 5/3. The frequency eigenvalues of various branches of modes, in units of the local angular velocity, are plotted against the radial wavenumber, in units of H −1. The two f modes, and the first five p, g and r modes of each parity are shown. In order of increasing frequency these are ro 3, ge 2, 3, pe 2, ge go

It is exponentially small as x1 → +∞, thereby justifying the layer analysis, if and only if n is a positive integer, in which case the conﬂuent hypergeometric function is proportional to a generalized Laguerre polynomial of degree n − 1. Equation (5.11) gives the corresponding eigenvalues λ± n as the roots of a quadratic equation. Provided that γ > Γ > 1, these roots are real and ordered such that

as n → ∞. The modes corresponding to λ+ additional mode which corresponds to the trivial solution ∆0 = 0 of equation (5.8). This has

4, ge 2, ro 5, although even and odd f, p and g modes coalesce for suﬃciently large k.

(kζs)−1/2 while ξφ = O Lamb (1932), can be combined to obtain

and λ = λ0 = 1. It is a surface gravity mode and corresponds to the f mode of stellar oscillations.

. The quantity λ is the dimensionless eigenvalue of the leading-order equations, which, following

n are called pn and gn respectively in stellar oscillations. There is an

is the leading-order part of ∇·ξ, while c and n are given by

The solution regular at x1 = 0 is (e.g. Erd´elyi et al. 1953a)

Γλ2 − γ [2(Γ − 1)n + 1] λ + (γ − Γ) = 0.

c(cid:13) 1994 RAS, MNRAS 000, 000-000

[2(n − 1) + c − x1] ∆0 = 0,

2γ(Γ − 1) Γ

γ − Γ 2γ(Γ − 1)

1F1 (1 − n; c; 2x1) .

1 < 1 < λ+

2Γ − 1 Γ − 1

(cid:0) ∂∆0 ∂x1

n and λ−

∂2∆0 ∂x2

ξr0 ∝ ξz0 ∝ e

1, f e, f o, pe

∆0 = iξr0 −

1 < λ+

∂ξz0 ∂x1

2 < … ,

− 2 < λ

∆0 ∝ e

… < λ

λ+ n ∼

− n ∼

where

(5.10)

(5.15)

(5.14)

(5.11)

(5.12)

(5.13)

4, po

1, po

5, po

2, po

3, po

3, go

4, go

5, go

1, go

4, pe

1, pe

5, ge

2, pe

5, ge

1, ro

4, ro

3, ro

3, re

1, re

5, re

2, re

4, re

(5.8)

(5.9)

with

and

c =

−x1

(cid:1)

−

(cid:1)

(cid:0)

while

(5.17)

(5.16)

n, ge

n, po

− n →

(kζs)−1

− n → 0.

n and go

and scalings

(kζs)−1/2

ξr(r, ζ) ∼ (kζs)

−1/2ξr0(r, x2),

ξφ(r, ζ) ∼ ξφ0(r, x2),

λ+ n → ∞ while

Conversely, when γ → Γ,

λ+ n → 2(Γ − 1)n + 1

2(Γ − 1)n + 1

(cid:0) x2 = (kζs)1/2(ζ/ζs)

centred on the equatorial plane, and have ω2 = O

Waves and instabilities in a diﬀerentially rotating disc

a layer of extent O in a similar way, with a stretched coordinate

mode has ∇·ξ = 0 at leading order. Nevertheless, the g modes do survive the incompressible limit. When γ → ∞,

Somewhat surprisingly, the p modes and g modes are equally compressive, sharing the same function ∆0, and only the f

The r modes behave quite diﬀerently in the limit kζs → ∞. In a convectively stable disc, the r modes become trapped in . The asymptotic analysis proceeds

The notation ‘f’, ‘p’, ‘g’, etc., for stellar oscillations can be readily applied to thin accretion discs provided that an additional distinction is made between even and odd modes. This distinction is lost in the limit kζs → ∞, as described above. For general values of k, one should refer to the modes as f e, f o, pe n, where n ∈ IN is a vertical mode number and the superfix denotes the parity of the mode.6

and HN is an Hermite polynomial (e.g. Erd´elyi et al. 1953b). If N is odd, the mode is even and may be called re n, where N = 2n − 1. If N is even, the mode is odd and may be called ro n, where N = 2n − 2. The r modes have ∇·ξ = 0 at leading order, and their limiting behaviour depends as much on buoyancy as on inertial forces. In the marginally stable case γ = Γ, the r modes have ω2 = O

6 An alternative classification scheme is possible in which pe n are renamed p2n and p2n−1, respectively, and similarly for other classes of modes. This is closer to the notation of KP, and has the advantage that the superfix is not required. However, the frequency eigenvalues do not then fall in sequence for all classes of modes. Also, to refer to f1 and f2 rather than f o and f e is perhaps less helpful since it obscures the uniqueness of the f mode.

These asymptotic results are all verified numerically, as shown in Fig. 2. It should be emphasized that this classification scheme is based entirely on the behaviour of the modes for large values of kζs and may not accurately reﬂect the properties of the modes for smaller kζs.

where N is a non-negative integer; the solution is then

and the eigenfunctions are not localized.

A bounded solution exists if and only if

The leading-order equations are

c(cid:13) 1994 RAS, MNRAS 000, 000-000

and ∂2ξz0 ∂x2

2(γ − Γ) γ(Γ − 1)

ξz(r, ζ) ∼ ξz0(r, x2)

ξz0 ∝ exp(− 1

λ = (2N + 1)

ω2 ∼ (kζs)

2)HN (αx2),

ξφ0 = −2iλ

−1/2ξr0

n and po

∂ξz0 ∂x2

−1Ω2λ.

(kζs)−2

2 α2x2

ξz0 = 0.

iξr0 +

where

(5.25)

(5.26)

(5.19)

(5.20)

(5.21)

(5.22)

(5.23)

(5.24)

(5.27)

(5.28)

(5.18)

= 0,

α =

and

λ −

1/4

(cid:21)

(cid:20)

(cid:21)

(cid:20)

(cid:0)

(cid:1)

16 G. I. Ogilvie

The presence of a poloidal magnetic field in the disc has profound implications for the spectrum of waves and instabilities. As well as modifying the modes that occur in a hydrodynamic disc, the magnetic field gives rise to additional modes which are potentially unstable. Furthermore, the branches are no longer separated in the dispersion diagram, but undergo avoided cross- ings. This makes the classification of modes diﬃcult, although there is one limit which can be investigated semi-analytically. This is the case k = 0 in a disc with a purely vertical magnetic field, as was considered by Gammie & Balbus (1994).

When the magnetic field is purely vertical, the Lorentz force vanishes and equations (5.1)-(5.3) for the equilibrium are valid. If also k = 0, equations (3.13)-(3.15) simplify considerably because the horizontal and vertical components of ξ become decoupled. There exist purely horizontal modes, with

Figure 2. Asymptotic behaviour of the local dispersion relation for large radial wavenumber k, for the same disc as in Fig. 1. The f modes and p modes (top left) and the g modes (top right) have ω2 = O(k) as k → ∞, while the r modes (bottom) have ω2 = O(k−1). In each case the dotted lines indicate the asymptotic limits derived in the text. Note that the vertical axis is diﬀerent in each plot.

c(cid:13) 1994 RAS, MNRAS 000, 000-000

5.2 Magnetized discs

ξφ(r, ζ) = aφ(r)uN (r, ζ),

ξr(r, ζ) = ar(r)uN (r, ζ)

(5.29)

(5.30)

and

(cid:1)

(cid:0)

(cid:21)

(cid:20)

(cid:21)

1/2

1 ±

viz.

ar aφ

2 Ω2

(5.33)

(5.32)

(5.31)

N (ω2

N + 1

N /Ω2

ω2 = ω2

ω2 = Ω2

1 + 16ω2

ω4 − (2ω2

N − 3Ω2) = 0,

n if N = 2n, ro

N + Ω2)ω2 + ω2

n if N = 2n, or mo

N −2iωΩ ω2 − ω2

ω2 + 3Ω2 − ω2 2iωΩ

N < 3Ω2. Since the eigenvalues {ω2

16 Ω2. In the limit Bz → 0, the eigenvalues ω2

Waves and instabilities in a diﬀerentially rotating disc

where uN is one of the eigenfunctions of equation (4.1), and the components ar and aφ satisfy the equation

N (cid:21) (cid:20) (cid:20) The frequency eigenvalues are the roots of the equation

N } are ordered, the N = 1 mode 4 Ω, which is achieved if N → 0 (although not uniformly with respect to N ), and equation (5.33)

so that there is an unstable mode if and only if 0 < ω2 is unstable if any mode is unstable. The maximum possible growth rate is the Oort parameter A = 3 ω2 N = 15 reduces to

Accordingly, the lesser solution of equation (5.33) may be called an m mode, since it is due entirely to the magnetic field. Specifically, it is me n if N = 2n − 1. If N = 0 the mode is trivial. The greater solution of equation (5.33) may be called re n if N = 2n − 1, or f e if N = 0. However, the r and g modes really lose their identity when a magnetic field is introduced.

An example of the dispersion diagram for a disc containing a purely vertical magnetic field is shown in Fig. 3. In this case the strength of the magnetic field is such that one unstable magnetorotational mode of each parity exists. These modes are unstable only for suﬃciently small k, and their growth rates are greatest for k = 0. The stable modes form a complicated dispersion diagram because of the avoided crossings that occur between diﬀerent branches. Broadly speaking, the pattern consists of branches which, like the f and p modes in a hydrodynamic disc, sweep upwards in the diagram, and other branches which are almost ﬂat. Where these would cross over each other, a closer examination reveals that each branch is in fact continuous; the two branches approach and then diverge without crossing. The character of the eigenfunctions is exchanged smoothly in the neighbourhood of these avoided crossings. This means that a classification of the modes based on continuity with k = 0 would not be meaningful. The reason for plotting even and odd modes separately is that avoided crossings occur only between modes of equal parity. The fact that parity is a discrete property which cannot be exchanged smoothly implies that modes of opposite parity cannot interact in this way, and so branches of even and odd modes cross freely over one another. Similar avoided crossings in the spectrum of a magnetized isothermal atmosphere (without the symmetry of reﬂection in a

N is a Gegenbauer polynomial (e.g. Erd´elyi et al. 1953b). Equation (5.37) then leads to the dispersion relation

1 modes in a Keplerian disc without a magnetic field involve both horizontal and vertical displacements even when

Appropriate solutions, regular at ζ = ±ζs, are obtained when N is a non-negative integer, in which case

There are also purely vertical modes, which satisfy the equation

These modes may be identified as pe

7 In fact the f o and ro k = 0.

∂2ξz ∂ζ 2 − (1 + 2q)ζ

n if N = 2n, or f o if N = 0.7

c(cid:13) 1994 RAS, MNRAS 000, 000-000

n if N = 2n − 1, po

N (N + 2q)ξz = 0,

ω2 − Ω2 Ω2

Γ + 1 2(Γ − 1)

ω2/Ω2 = 1 + γ

N (N + 2q) =

Γ − 1 2Γ

Γ + 1 Γ − 1

2Γ Γ − 1

where C q

N (ζ/ζs),

ξz ∝ C q

s − ζ 2)

∂ξz ∂ζ

where

(5.35)

(5.36)

(5.37)

(5.38)

(5.39)

(5.34)

(cid:17) (cid:18)

N +

and

(ζ 2

q =

(cid:17)i

(cid:19)

(cid:16)

(cid:17)

(cid:16)

18 G. I. Ogilvie

Figure 3. Part of the local dispersion relation for a thin disc containing a purely vertical magnetic field. The parameters are Γ = 4/3, γ = 5/3 and Bz = 0.01. Modes with even and odd symmetry are shown in the left and right panels, respectively. The frequency eigenvalues of various branches of modes, in units of the local angular velocity, are plotted against the radial wavenumber, in units of H −1. A dashed line indicates an unstable mode for which ω is imaginary. In this case one unstable magnetorotational mode of each parity exists. A maximum growth rate of 0.7151Ω is achieved at k = 0 by the mo 1 mode. The frequency eigenvalues of the continuous spectrum, which are all real, are marked on the scale at the right of each panel; these are the eventual limits of the branches as k → ∞.

When the magnetic field is weaker, more unstable modes exist but their branches fit neatly without mode interactions, as shown in Fig. 4. It is well known that the addition of a weak magnetic field to a hydrodynamic disc constitutes a highly singular perturbation if ideal MHD is assumed. This is reﬂected in the fact that the convergence of the frequency eigenvalues of mn modes to zero as B → 0 is not uniform with respect to n. Modes of large n clutter the dispersion diagram even when the magnetic field is exceedingly weak, and can be unstable with large growth rates.

Figure 4. Part of the local dispersion relation for a thin disc containing a purely vertical magnetic field. The parameters are Γ = 4/3, γ = 5/3 and Bz = 0.002. Modes with even and odd symmetry are shown in the left and right panels, respectively. In this case seven unstable magnetorotational modes of each parity exist. A maximum growth rate of 0.7495Ω is achieved at k = 0 by the me 4 mode. The frequency eigenvalues of the continuous spectrum, which are all real, are marked on the scale at the right of each panel.

Finally, an example of the dispersion diagram for a disc containing a bending poloidal magnetic field is shown in Fig. 5.

horizontal plane) have been analysed by Hasan & Christensen-Dalsgaard (1992).

c(cid:13) 1994 RAS, MNRAS 000, 000-000

Waves and instabilities in a diﬀerentially rotating disc

Figure 5. Part of the local dispersion relation for a thin disc containing a bending poloidal magnetic field. The parameters are Γ = 4/3, γ = 5/3, Bz = 0.015 and Brs = 0.005, with Ω1s ≈ 0.3376. Modes with even and odd symmetry are shown in the left and right panels, respectively. In this case only the mo 1 mode is unstable, achieving a maximum growth rate of 0.6019Ω at k = 0. The frequency eigenvalues of the continuous spectrum, which are all real, are marked on the scale at the right of each panel.

Before describing the marginal curves, it is appropriate to give a more detailed account of the solutions of the equilibrium equations for polytropic discs. When Brs and Bz are specified, equations (3.1)-(3.5) constitute a non-linear eigenvalue problem for Ω1s. The solutions of these equations lie on a two-dimensional manifold in the three-dimensional parameter space of (Brs, Bz, Ω1s). Some of these solutions are physically acceptable and may be called ‘regular’ equilibria: the magnetic field lines bend only once when passing through the disc, and the density and pressure decrease monotonically from the equatorial plane to the surface. In the remaining solutions, which may be called ‘irregular’ equilibria, the field lines bend more than once; some of these solutions also have density inversions. In Paper I only the class of regular equilibria (described as the ‘principal branch’) for the case Γ = 5/3 was mentioned. It is now to be shown that all irregular equilibria are unstable, so that the marginal curves need be drawn only on the principal branch.

The stability of a weakly magnetized thin disc to axisymmetric perturbations can be decided, in principle, by computing the local dispersion relation at each radius separately; the disc is unstable if and only if an imaginary frequency eigenvalue exists, for any real value of k, at any radius. However, a more eﬃcient algorithm is required in practice. In the case of polytropic discs, the parameter space divides into stable and unstable regions separated by a curve of marginal stability. The location of the marginal curve depends on a further parameter which is the adiabatic exponent γ. Equilibria on the marginal curve possess a mode with ω = 0, but this is not suﬃcient to determine the curve, since all unstable equilibria also possess modes with ω = 0 (cf. Figs 3-5). In principle, a marginal curve can be drawn for each m mode and for every value of k.

The behaviour of branches in the limit k → ∞ is quite diﬀerent from the hydrodynamic case. If any one branch is followed continuously, it undergoes a finite number of avoided crossings and then ω approaches a finite limit. It can be shown that these limits are the frequency eigenvalues of the continuous spectrum discussed in Section 2. In this limit the Lagrangian displacement is confined entirely within the magnetic surface and the magnetorotational instability, which is associated with the gradient of angular velocity perpendicular to the magnetic field, is lost.

Proof. When k = 0, the form of the functional R[ξr, ξz; k] for purely horizontal trial displacements of odd symmetry is, from equation (4.13),

The avoided crossings are wider than for a purely vertical magnetic field, which implies that the coupling between diﬀerent branches of modes is stronger when the field lines bend.

Theorem 1. All irregular equilibria possess an unstable mode of odd symmetry at k = 0.

STABILITY CRITERIA FOR MAGNETIZED DISCS

c(cid:13) 1994 RAS, MNRAS 000, 000-000

− 3ρΩ2|ξr|2

R[ξr, 0; 0] =

∂ξr ∂ζ

−ζs

(6.1)

B2 z

dζ,

(cid:12) (cid:12) (cid:12) (cid:12)

ζs

= 0

B2 z

(6.4)

(6.3)

(6.2)

∂y ∂ζ

f (ζ) =

ζ = ±ζ1.

20 G. I. Ogilvie

y1(ζ), sgn(ζ)y1(ζ1),

∂2y ∂ζ 2 + λρy = 0

∂2Ω1 ∂ζ 2 + 3ρΩ2Ω1 = 0,

on −ζ1 < ζ < ζ1, subject to the boundary conditions

and the trial function must satisfy ∂ξr/∂ζ = 0 at ζ = ±ζs. Given that the equilibrium is irregular, Br must have zeros, other than at ζ = 0, at points ζ = ±ζn, where 0 < ζ1 < ζ2 < · · · < ζN < ζs and N ≥ 1. Then consider the Sturm-Liouville equation

Since the equation is symmetrical about ζ = 0, the eigenfunctions y0, y1, y2, … alternate between even and odd symmetry, and the eigenvalues λ0, λ1, λ2, … form an ordered, increasing sequence of distinct, non-negative real numbers. The lowest eigenvalue is λ0 = 0, corresponding to y0 = constant. Now it follows from equations (3.1)-(3.4) for the equilibrium that Ω1 satisfies the equation

and also the boundary conditions (6.3). Moreover, since ζ1 is the first zero of Br in 0 < ζ < ζs, Br has exactly two extrema, and Ω1 has exactly two zeros, in the domain of the Sturm-Liouville equation. It follows that Ω1 is the eigenfunction y2, and has eigenvalue λ2 = 3Ω2. The properties of the Sturm-Liouville equation then ensure that the first odd eigenfunction y1 has an eigenvalue satisfying the (strict) inequalities 0 < λ1 < 3Ω2. Now construct the trial function

In Fig. 6, the diﬀerent branches of equilibria, now for the case Γ = 4/3, are shown in a number of suitably chosen sections through the parameter space. It is found that, as Bz is reduced, an increasingly large number of branches of irregular equilibria appear as the solution manifold folds and divides into many leaves. There is a critical value of Bz, approximately 0.01283, below which the principal branch becomes separated from the ‘vertical’ solution (which has Br = 0 and Ω1s = 0) by a stretch of irregular equilibria. (The vertical solution, however, is technically always regular.) As Bz is reduced further, a bifurcation disconnects the principal branch from the vertical solution. An infinite number of further bifurcations occur as Bz → 0. It is observed in the numerical calculations that, for regular equilibria, the last mode to be stabilized is the mo

where the trial displacement ξ is defined on −∞ < ζ < ∞ and is subject to the constraint ∆ = 0. Given that there exists an unstable eigenfunction ˆξ(ζ) at some wavenumber ˆk 6= 0, it follows that

which is continuous throughout −ζs < ζ < ζs and satisfies the boundary conditions for the variational principle. After an integration by parts one obtains

When k = 0, the constraint ∆ = 0 can be satisfied by using a purely horizontal trial displacement ξr = f , where the function f (ζ) is defined by

While this conjecture has not yet been proved in the most general circumstances envisaged, the following two theorems support it strongly.

8 When γ < Γ there may be a (magneto-) convective instability and the result cannot be expected to hold. In fact, a condition diﬀering slightly from γ ≥ Γ may be required in order to prove the conjecture.

Theorem 2. The conjecture holds in the case of an incompressible ﬂuid, provided there is no density inversion.

If a regular equilibrium possesses an unstable mode at any value of k, then it possesses an unstable mode at

Proof. An appropriate version of the functional R[ξr, ξz; k] for an incompressible ﬂuid is

and it is stabilized last at k = 0. This supports the following conjecture.

Conjecture. k = 0, assuming that γ ≥ Γ.8

and, since this is clearly negative, the theorem is proved.

|Dξr|2 + |Dξz|2 − 3ρΩ2|ξr|2 − Ω2ζ

c(cid:13) 1994 RAS, MNRAS 000, 000-000

R[ξr, ξz; k] = 4Ω2|a0|2 +

(3ρΩ2 − λ1)|f |2 dζ − 2

ˆξr(ζ) = f (ζ) exp

R[ ˆξr, ˆξz; ˆk] < 0.

|ζ| < ζ1, |ζ| ≥ ζ1,

−(2Ω1/3Ω)iˆkr

3ρΩ2|f |2 dζ,

R[f, 0; 0] = −2

−∞ (cid:18) Z

1 mode,

∂ρ ∂ζ

|ξz|2

(6.5)

(6.6)

(6.9)

(6.7)

(6.8)

dζ,

(cid:26)

(cid:19)

∞

(cid:2)

(cid:3)

ζ1

ζs

Waves and instabilities in a diﬀerentially rotating disc

Figure 6. Solution curves for weakly magnetized, polytropic thin discs with Γ = 4/3. The ‘eigenvalue’ Ω1s is plotted against Brs for various values of Bz (as indicated in each panel). All solutions in the range −0.025 < Brs < 0.025 are shown. Regular and irregular equilibria are indicated by solid and dotted lines, respectively. As Bz is decreased, the solution manifold folds and divides into many leaves, and the branch of regular equilibria becomes separated from the solution at the origin which has a purely vertical magnetic field.

c(cid:13) 1994 RAS, MNRAS 000, 000-000

ˆξr = iˆkBr ˆξr + Bz

−(2Ω1/3Ω)iˆkr

∂ ˆξr ∂ζ

so that

(6.10)

D(0)f

∂f ∂ζ

= Bz

D(ˆk)

exp

(cid:18)

(cid:19)

(cid:1)

(cid:0)

(cid:3)

(cid:2)

(cid:3)

∞

−

(cid:19)

dζ

(cid:12) (cid:12) (cid:12) (cid:12)

and

| ˆξz|2

(6.13)

(6.12)

(6.11)

|D(ˆk)

∂ρ ∂ζ

|D(ˆk)

−∞ (

Proof.

∂ξr ∂ζ

∂ξz ∂ζ

−∞ (cid:18) Z

ρΩ2ζ γp

ˆξz|2 − Ω2ζ

| ˆξr|2 = |f |2

−3ρΩ2|ξr|2 −

γp γp + B2

ˆξr|2 = |D(0)f |2.

is non-negative, and so

In this case one may take

22 G. I. Ogilvie

R[ ˆξr, ˆξz; ˆk] − R[f, 0; 0] =

R[f, 0; 0] ≤ R[ ˆξr, ˆξz; ˆk] < 0,

|δΠ|2 γp + B2 + B2

R[ξr, ξz; k] =4Ω2|a0|2 ∞

which proves that an unstable mode must exist with k = 0.

where the suﬃx on the operator D indicates the value of k used. Then

Theorem 3. The conjecture holds for a compressible ﬂuid when the magnetic field is purely vertical.

Also, the coeﬃcients {an} in the eigenfunction expansion (4.8) are the same when ξr = f and k = 0 as they are when ξr = ˆξr and k = ˆk. Thus the diﬀerence

If the conjecture is accepted on the basis of these lemmata and the numerical evidence, it implies that the overall stability boundary in the parameter space is the marginal curve for a mode (in fact, the mo 1 mode) at k = 0, drawn on the principal branch. This curve has been computed directly by solving the equations for an equilibrium (with a given value of Bz) and a mode (with ω = 0 and k = 0) as an eigenvalue problem, yielding a value for Brs. The stability boundaries for equilibria with Γ = 4/3 and Γ = 5/3 are shown in Fig. 7, using γ = 5/3 for the adiabatic exponent. In comparing the two panels it is important to note that the unit of magnetic field strength depends on Γ, according to equation (4.22). It is more convenient, therefore, to compare a dimensionless quantity such as the plasma beta β = 2p/B2. The plasma beta of the marginal equilibria, evaluated on the equatorial plane, is plotted against Bz in Fig. 8. It is seen that all marginal equilibria have values of β reasonably close to unity. However, β is not generally a reliable guide to stability, especially for equilibria in which the angle of inclination of the magnetic field exceeds π/6.

The emergence and disappearance of equilibrium solutions as the parameters are varied is not necessarily governed by conventional bifurcation theory. Typically, an equilibrium ceases to exist when the parameters are changed in such a way that would make the enthalpy become negative at some point. This would correspond to a branch point of the polytropic relation (3.5), and non-analytic behaviour is to be expected. In fact, the equilibrium may continue to exist in a complex-valued sense, even though this has no physical meaning. This means that there is little hope of using the bifurcations of the solution manifold to infer the stability of the solutions. Indeed, two further obstacles to such a method present themselves: first, only the stability to modes of even parity could be considered, since the equilibria are all symmetric; and secondly, the eﬀect of the adiabatic exponent γ on stability could not be taken into account.

When k = 0, the purely horizontal displacement ξr = ˆξr should be used, and δΠ then vanishes. The coeﬃcient of |ξz|2 in the integrand of equation (6.14),

The magnetorotational modes are not aﬀected greatly by buoyancy or compressibility. In Fig. 9 the stability boundaries for equilibria with Γ = 4/3 and Γ = 5/3 are plotted for diﬀerent values of the adiabatic exponent. The unstable region reduces

where δΠ = −(γp + B2)∆ + ρΩ2ζξz + B2 ∂ξz ∂ζ

Given that there exists an unstable eigenfunction ˆξ(ζ) at some wavenumber ˆk 6= 0, it follows that

(cid:16) is non-negative, since it is assumed that γ ≥ Γ. Thus

which proves that an unstable mode must exist with k = 0.

R[ ˆξr, 0; 0] ≤ R[ ˆξr, ˆξz; ˆk] < 0,

c(cid:13) 1994 RAS, MNRAS 000, 000-000

(cid:19) (ρΩ2ζ)2 γp

(cid:12) (cid:12) (cid:12) Ω2ζ (cid:12) (cid:20)

R[ ˆξr, ˆξz; ˆk] < 0.

(ρΩ2ζ)2 γp

Ω2ζ (cid:20)

(cid:12) (cid:12) (cid:12) |ξz|2 (cid:12)

γ − Γ Γ

(cid:18) ∂ρ ∂ζ

∂ρ ∂ζ

(6.14)

(6.15)

(6.16)

(6.17)

dζ,

(cid:12) (cid:12) (cid:12) (cid:12)

(cid:27)

(cid:18)

(cid:19)

ξz

−

(cid:17)

(cid:21)

Waves and instabilities in a diﬀerentially rotating disc

Figure 7. Stability boundaries in the parameter spaces for weakly magnetized, polytropic thin discs with Γ = 4/3 (left) and Γ = 5/3 (right). In each case, the dotted line indicates an angle of inclination i = π/6. The dashed line is the critical curve for the existence of the principal solution branch. The solid line is the marginal curve, at k = 0, for the mo 1 mode, when γ = 5/3. This is the overall stability boundary of the equilibria to axisymmetric perturbations. In the regions R, which extend indefinitely beyond the upper right of each figure, there exist equilibria which are capable of driving a wind and are also stable to the magnetorotational instability. Note the diﬀerent scales used in each plot.

Figure 8. The plasma beta of the marginal equilibria, evaluated on the equatorial plane, plotted against the vertical magnetic field strength in the interval in which a marginal equilibrium exists. The two panels correspond to the two panels of Fig. 7, i.e. equilibria with Γ = 4/3 (left) and with Γ = 5/3 (right), where γ = 5/3 in each case. Marginal equilibria typically have values of β that are reasonably close to unity.

slightly in size as γ is increased from Γ to ∞, demonstrating the mildly stabilizing eﬀect of a sub-adiabatic stratification. The eﬀect is most pronounced for equilibria in which the magnetic field bends significantly, and is zero for those with a purely vertical magnetic field.

The analysis of Section 3 may be extended to include non-axisymmetric waves and instabilities of azimuthal wavenumber m. A distinction must be made between modes with m = O(ǫ−1) and those with m = O(1). In the former case, the azimuthal

7 NON-AXISYMMETRIC MODES

c(cid:13) 1994 RAS, MNRAS 000, 000-000

24 G. I. Ogilvie

Figure 9. The eﬀect of sub-adiabatic stratification on the magnetorotational instability. Left: stability boundaries for weakly magnetized, polytropic thin discs with Γ = 4/3 are plotted for four diﬀerent values of the adiabatic exponent γ. The four solid lines are the stability boundaries for γ = 4/3, 5/3, 3 and ∞ (from right to left). Right: similar results for equilibria with Γ = 5/3. The three solid lines are the stability boundaries for γ = 5/3, 3 and ∞ (from right to left).

There remains the possibility that the dispersion relation yields imaginary values of ω for real values of k, as a result of the magnetorotational instability or a convective instability. In that case one cannot simply replace ω with ˆω, for the following reason. Since ω must be constant for a normal mode, it follows that the real part of ˆω can vanish at (at most) one radius, the corotation radius of the mode. Therefore ˆω is imaginary at the corotation radius, but anywhere else it cannot be a solution of the dispersion relation for any real value of k. The conclusion is that, for a non-axisymmetric unstable mode, k is real only at the corotation radius.

wavenumber is comparable to the radial and vertical wavenumbers, and the equations of Section 3 are not valid in any sense. These modes are expected to resemble the global non-axisymmetric instabilities of thick tori studied by Papaloizou & Pringle (1984). The strong diﬀerential Doppler shift means that these modes are probably localized in a region of small radial extent about their corotation radius, and require a boundary of the disc to be present in this small region if they are to have dynamical [O(1)] growth rates.

which is an analytic function apart from a branch cut along the imaginary axis. It then follows that each branch of the dispersion relation is an analytic function, except on Re(k) = 0. The derivative ∂ω/∂k is known for real values of k from the real dispersion relation, and is imaginary for an unstable branch. Let k = kc be the real value of the radial wavenumber at the corotation radius r = rc of such a mode. Then, at a neighbouring radius r, the wavenumber is (to first order)

and the frequency eigenvalue ω must be replaced with the intrinsic frequency ˆω = ω − mΩ. In this way LP were able to discuss the propagation of non-axisymmetric waves in an isothermal accretion disc. Provided that the dispersion relation yields real values of ˆω for real values of k, there is no diﬃculty because ω, ˆω and k are all real within the wave region of a mode.

assuming that the derivative (∂ω/∂k)c does not vanish. The quantity in square brackets is purely imaginary; depending on the sign of the various terms, Im(k) is either an increasing function or a decreasing function of (r − rc). If it is an increasing

For modes with m = O(1), however, the equations of Section 3 require very few changes. The eigenfunctions have the

When k is not real, the variational principle of Section 4.2 is not valid. Otherwise, the only change that need be made in

Section 3 is to replace |k| with the function

+k, Re(k) > 0, −k, Re(k) < 0,

c(cid:13) 1994 RAS, MNRAS 000, 000-000

−iωt + imφ + iǫ

ξ0(r, ζ) exp

3mΩc 2rc

ξ(r, t) ∼ Re

k ≈ kc +

(r − rc),

k(r) dr

p(k) =

∂ω ∂k

c (cid:21) (cid:17)

(cid:30)(cid:16)

(7.1)

(7.3)

(7.2)

form

(cid:21)(cid:27)

(cid:26)

−1

(cid:20)

8 DISCUSSION

Waves and instabilities in a diﬀerentially rotating disc

function, then the mode will take the form of a localized Gaussian wave packet around r = rc, according to equation (7.1). If it is a decreasing function, the function (7.1) would increase greatly in magnitude away from r = rc and cannot satisfy the radial boundary conditions.

In conclusion, non-axisymmetric unstable modes are localized around the corotation radius, and obey the real dispersion relation only at that point. The modes have either positive or negative values of m depending on the sign of the imaginary group velocity. At any given point on an unstable branch of the real dispersion relation, non-axisymmetric modes of this kind can be found.

This paper has addressed a number of issues relating to the spectrum of waves and instabilities in an accretion disc containing a poloidal magnetic field. The general analysis of the continuous spectrum in Section 2 demonstrates that the only type of possible instability that is truly localized on a single magnetic surface is an essentially axisymmetric magnetoconvective instability associated with the component of gravity parallel to the magnetic field. Unlike the case of uniform (or zero) rotation, the interchange instability, if present, is not manifest in the continuous spectrum. Neither is the magnetorotational instability or any other instability associated with diﬀerential rotation. These must instead be sought in the normal modes of the system. The extension of the asymptotic methods of Paper I leads to a WKB description of axisymmetric waves and instabilities in a weakly magnetized thin disc in terms of a local dispersion relation which generalizes earlier work by LP and KP. The waves propagate radially as in a slowly varying waveguide. The modes in a hydrodynamic disc can be classified as f, p, g and r modes according to their behaviour in the limit of large radial wavenumber. When a magnetic field is introduced, the dispersion diagram is complicated by a large number of avoided crossings which make a systematic classification diﬃcult. If the magnetic field is relatively weak, unstable magnetorotational modes also occur for suﬃciently small values of the radial wavenumber. The overall stability boundary in the parameter space of polytropic equilibria is obtained by computing the marginal curve for the first magnetorotational mode of odd symmetry. For a given value of the vertical magnetic field, the addition of a radial component (which makes the field lines bend) has a stabilizing inﬂuence. A sub-adiabatic stratification also has a mild stabilizing inﬂuence if the magnetic field is not purely vertical.

It is also important to note that the analysis is valid only for axisymmetric modes and non-axisymmetric modes of small azimuthal wavenumber [m = O(1)], and that it can only detect instabilities with dynamical [O(1)] growth rates. It is expected that all models are subject to global non-axisymmetric instabilities which depend strongly on the radial boundaries of the disc, resembling the Papaloizou-Pringle instability, but these probably have very small growth rates in a thin disc. The radial interchange instability (Spruit et al. 1995), if it is not entirely stabilized by the diﬀerential rotation, would also have a sub-dynamical growth rate in a weakly magnetized thin disc. The equilibria described here as ‘stable’ may therefore not exist as truly laminar ﬂows but could be subject to weak turbulence or at least ﬂuctuations. However, that would be quite diﬀerent from equilibria that are locally unstable to the magnetorotational instability, which are bound to degenerate into strong MHD turbulence.

There are many ways in which this analysis could be improved and extended. In order to study the radial propagation of waves, an explicit global equilibrium model of a magnetized disc should be constructed by choosing the functional forms of ζs(r), ψ0(r) and K0(r), say, and solving for the equilibrium at each radius. The variation of the radial wavenumber with radius could then be determined for any mode by following the dispersion relation. The dependence of the amplitude of the mode on radius could be found either by solving the linearized equations at the next order in ǫ, or, more simply, by appealing to a wave-action conservation relation (cf. LP). Other areas to be explored are the possible stabilization of a super-adiabatically stratified disc by a suﬃciently strong magnetic field, the propagation of the m = 1 ‘tilt’ mode in a magnetized disc, and the eﬀects of self-gravitation on the equilibria and spectra of both hydrodynamic and magnetized discs.

This analysis shows for the first time that it is possible to construct stable equilibria which are capable of driving a wind. Indeed, increasing the magnetic field strength not only tends to stabilize the equilibria but also makes it easier to construct equilibria in which the magnetic field lines at the surface of the disc are inclined to the vertical at angles significantly greater than π/6. It should be emphasized, however, that the stability analysis applies only to the disc and not to any wind solution that may be superimposed on it.

I would like to thank Jim Pringle, Douglas Gough and Ulf Torkelsson for helpful discussions. A research studentship from the Particle Physics and Astronomy Research Council is acknowledged.

Agapitou V., Papaloizou J. C. B., Terquem C., 1997, MNRAS, 292, 631

c(cid:13) 1994 RAS, MNRAS 000, 000-000

ACKNOWLEDGMENTS

REFERENCES

26 G. I. Ogilvie

Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214 Bernstein I. B., Frieman E. A., Kruskal M. D., Kulsrud R. M., 1958, Proc. R. Soc. Lond. A, 244, 17 Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883 Brandenburg A., Nordlund ˚A., Stein R. F., Torkelsson U., 1995, ApJ, 446, 741 Brandenburg A., Nordlund ˚A., Stein R. F., Torkelsson U., 1996, ApJ, 458, L45 Chandrasekhar S., 1960, Proc. Natl Acad. Sci., 46, 253 Christensen-Dalsgaard J., 1980, MNRAS, 190, 765 Courant R., Hilbert D., 1953, Methods of Mathematical Physics, vol. 1. Interscience, New York Curry C., Pudritz R. E., 1996, MNRAS, 281, 119 Erd´elyi A., Magnus W., Oberhettinger F., Tricomi F. G., 1953a, Higher Transcendental Functions, vol. 1. McGraw-Hill, New York Erd´elyi A., Magnus W., Oberhettinger F., Tricomi F. G., 1953b, Higher Transcendental Functions, vol. 2. McGraw-Hill, New York Foglizzo T., Tagger M., 1995, A&A, 301, 293 Frieman E., Rotenberg M., 1960, Rev. Mod. Phys., 32, 898 Gammie C. F., Balbus S. A., 1994, MNRAS, 270, 138 Hasan S. S., Christensen-Dalsgaard J., 1992, ApJ, 396, 311 Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742 Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 464, 690 Heyvaerts J. F., Norman C., 1989, ApJ, 347, 1055 Heyvaerts J. F., Priest E. R., 1989, A&A, 216, 230 Kippenhahn R., Schl¨uter A., 1957, Z. Astrophys., 43, 36 Korycansky D. G., Pringle J. E., 1995, MNRAS, 272, 618 (KP) Lamb H., 1932, Hydrodynamics, 6th edn. Cambridge Univ. Press, Cambridge Lin D. N. C., Papaloizou J. C. B., Kley W., 1993, ApJ, 416, 689 Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360 (LP) Moss D. L., Tayler R. J., 1969, MNRAS, 145, 217 Ogilvie G. I., 1997, MNRAS, 288, 63 (Paper I) Ogilvie G. I., Pringle J. E., 1996, MNRAS, 279, 152 Papaloizou J. C. B., Pringle J. E., 1982, MNRAS, 200, 49 Papaloizou J. C. B., Pringle J. E., 1984, MNRAS, 208, 721 Papaloizou J. C. B., Szuszkiewicz E., 1992, Geophys. Astrophys. Fluid Dyn., 66, 223 Poedts S., Hermans D., Goossens M., 1985, A&A, 151, 16 Pringle J. E., 1981, ARA&A, 19, 137 Ruden S. P., Papaloizou J. C. B., Lin D. N. C., 1988, ApJ, 329, 739 Spiegel E. A., Weiss N. O., 1982, Geophys. Astrophys. Fluid Dyn., 22, 219 Spruit H. C., Stehle R., Papaloizou J. C. B., 1995, MNRAS, 275, 1223 Stone J. M., Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 463, 656 Tassoul J.-L., 1978, Theory of Rotating Stars. Princeton Univ. Press, Princeton Tayler R. J., 1973, MNRAS, 161, 365 Terquem C., Papaloizou J. C. B., 1996, MNRAS, 279, 767 Velikhov E. P., 1959, Sov. Phys. JETP, 9, 995

This paper has been produced using the Royal Astronomical Society/Blackwell Science TEX macros.

c(cid:13) 1994 RAS, MNRAS 000, 000-000