ar Xiv:astro-ph/9801306v1 29 Jan 1998 Mon. Not. R. Astron. Soc.000, 000–000 (1994)Printed 24 May 2021(MN plain T E X macros v1.6) Waves and instabilities in a differentially rotating disc containing a poloidal magnetic field G. I. Ogilvie 1,2 1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA 2 Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Silver Street, Cambridge CB3 9EW The theory of waves and instabilities in a differentially 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 differential rotation. The normal modes of a weaklymagnetized 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 thedispersion diagram by causing avoided crossings. The stability boundary to the magnetorotational insta- bility in the parameter space of polytropic equilibria is located by solvingdirectly for marginally stable equilibria. For a given vertical magnetic field in the disc, bending of the field lines has a stabilizing effect and it is shown that stable equilibria exist which are capable of launching a predominantly centrifugally driven wind. Key words:accretion, accretion discs – hydrodynamics – instabilities – MHD – waves. 1 INTRODUCTION There are good reasons for believing that magnetic fields areimportant to the physics of accretion discs. Magnetohydrodynamic (MHD) mechanisms provide the most convincing explanationsfor 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 flows 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 flow, while hydrodynamically stable, is destabilized by the presence of a magnetic field, provided that the magnetic field is sufficiently weak and the disc is sufficiently 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 fluid in differential rotation about a massive central object. The fluid contains a purely poloidal magnetic field ofdipolar 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 differential equation in two dimensions, a version of the Grad–Shafranovequation. 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 radiusr. It was shown that two families of solutions resembling accretion discs exist in this limit, with different asymptotic scalings representing different balances of forces in the radial and vertical directions. Theweakly magnetized discsare the natural generalization of the standard, hydrodynamic thin discs (Pringle 1981). The sound speed and Alfv ́en velocity are both O(ǫ) [given that the c 1994 RAS 2G. I. Ogilvie Keplerian velocity is O(1)], while the fractional deviation from Keplerian rotation, caused by the radial Lorentz force, is O(ǫ). Thestrongly magnetized discsare 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 capableof launching a predominantly centrifugally driven wind if the magnetic field lines at the surface of the disc are inclined tothe vertical by an angle greater thanπ/6. 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 flow, the resistivity of the fluid, andany toroidal magnetic field are neglected in this model, it is expected that these additional terms constitute only smallperturbations 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ǫ. 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 & Papaloizou1995) 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. The energy principle of Bernstein et al. (1958) has played animportant 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 wavenumbermtends to infinity, and that the problem reduces to the consideration of a system of ordinary differential equations on each magnetic field line separately. Tayler (1973) considered the stability of a non-rotating star with a purely toroidal magnetic field, andfound 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 methodsoutlined in Section 2 below. In an accretion disc the differential 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 differentially rotating, non-self-gravitating fluid containing a poloidal magnetic field, a generalizationof 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 sufficient 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 arelatively long wavelength in the direction perpendicularto 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 briefly. 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, whichallow the modes to be enumerated and classified meaningfullyin 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 brieflyin 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 flux coordinates (ψ,φ,χ) defined c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc3 in Paper I. The magnetic flux coordinates form a right-handedorthogonal 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/|∇ψ||∇φ||∇χ|. 2 THE CONTINUOUS SPECTRUM The linearized equation governing the Lagrangian displacementξ(r,t) corresponding to a small departure from any state of ideal MHD may be written ρ D 2 ξ Dt 2 =−∇δΠ−(∇·ξ)∇Π−ξ·∇∇Π−ρξ·∇∇Φ +B·∇[B·∇ξ−(∇·ξ)B],(2.1) where Π is the total pressure and δΠ =−(γp+ 1 2 B 2 )∇·ξ−ξ·∇Π +B·(B·∇ξ)(2.2) is its linearized Eulerian perturbation. Hereρ,p,Bandγare the density, pressure, magnetic field and adiabatic exponent, respectively, and D/Dtis 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 flow. The self-gravitation of the fluid is neglected here, as it wasneglected throughout the construction of equilibria in Paper I. Although self-gravitation has, in general, a destabilizing influence on long-wavelength perturbations, its effect 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. The continuum modes may be described using a technique whichwas 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 effectively 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 flux coordinates (ψ,φ,χ) of Paper I, one considers a sequence of functions of the form ξ(r,t) = Re [ ξ(χ)f ( ψ−ψ 0 w(k) ) exp(ikψ+ imφ−iωt) ] ,(2.3) whereω∈ C is the frequency eigenvalue,m∈ZZthe azimuthal wavenumber, andk∈IR a parameter which is allowed to tend to infinity. The functionfis a unimodal ‘wavelet’ 1 of widthw(k) centered onψ=ψ 0 , and forms the envelope of the wave packet. In the limitk→ ∞, the functionfis to become infinitely localized atψ=ψ 0 , but the widthw(k) of the function should tend to zero more slowly thank −1 , perhapsw∝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 limitk→∞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 (anda finite value ofm), 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 ∇·ξ= ∆andδΠ =k −1 ̟,(2.4) where ∆ and̟are both O(1). In order to write equation (2.1) in magnetic flux coordinates, it is convenient to introduce the angle of inclinationi(ψ,χ) of the magnetic field to the vertical, such that cosi=e r ·e ψ =r B ∂r ∂ψ andsini=e r ·e χ = 1 JB ∂r ∂χ .(2.5) Then the derivatives of the unit vectors are given by de ψ =− ( ∂i ∂ψ ) e χ dψ+ cosie φ dφ− ( ∂i ∂χ ) e χ dχ,(2.6) de φ =−cosie ψ dφ−sinie χ dφ,(2.7) de χ = ( ∂i ∂ψ ) e ψ dψ+ sinie φ dφ+ ( ∂i ∂χ ) e ψ dχ,(2.8) 1 For example,f(x) = exp(−x 2 ). c 1994 RAS, MNRAS000, 000–000 4G. I. Ogilvie with ∂i ∂ψ = 1 JB ∂ ∂χ ( 1 r B ) and ∂i ∂χ =−r B ∂(JB) ∂ψ .(2.9) The quantity̟appears at leading order only in theψ-component of equation (2.1), which therefore serves to define̟ in terms ofξ φ ,ξ χ and ∆, but is otherwise unimportant. The remaining equations at leading order are −ρˆω 2 ξ φ −2iρˆωΩ siniξ χ −ρΩ 2 ξ φ =−cosir B ( ∂Π ∂ψ +ρ ∂Φ ∂ψ ) ξ φ r −sini 1 JB ( ∂Π ∂χ +ρ ∂Φ ∂χ ) ξ φ r + 1 J ∂ ∂χ ( 1 J ∂ξ φ ∂χ ) ,(2.10) −ρˆω 2 ξ χ + 2iρˆωΩ siniξ φ −ρΩ 2 sin 2 iξ χ =− 1 JB ∂Π ∂χ ∆− ξ χ JB ∂ ∂χ ( 1 JB ∂Π ∂χ ) −ρ ξ χ JB ∂ ∂χ ( 1 JB ∂Φ ∂χ ) −r B ( ∂Π ∂ψ −ρ ∂Φ ∂ψ ) ξ χ JB ( ∂i ∂χ ) + 1 J ∂ ∂χ ( 1 J ∂ξ χ ∂χ −B∆ ) − ( 1 J ∂i ∂χ ) 2 ξ χ (2.11) and 0 =−(γp+B 2 )∆− ξ χ JB ∂Π ∂χ + B J ∂ξ χ ∂χ ,(2.12) 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 differential equations withψtreated as a parameter. Once ∆ has been eliminated, the remaining equations may be written ˆω 2 ξ φ + 2iˆωΩ siniξ χ =− 1 ρr J ∂ ∂χ [ r 2 J ∂ ∂χ ( ξ φ r ) ] (2.13) and ˆω 2 ξ χ −2iˆωΩ siniξ φ =− 1 ρJB ∂ ∂χ [( v 2 s v 2 s +v 2 A ) B 2 J ∂ ∂χ ( ξ χ B ) ] + {( v 2 s v 2 s +v 2 A ) N 2 χ − 1 JB 2 ∂ ∂χ [( v 2 A v 2 s +v 2 A ) Bg χ ]} ξ χ ,(2.14) wherev s = (γp/ρ) 1/2 andv A = (B 2 /ρ) 1/2 are the sound speed and the Alfv ́en velocity, respectively, g χ =− 1 JB ∂Φ ∂χ +rΩ 2 sini(2.15) is the effective gravitational acceleration parallel to themagnetic field, and N 2 χ =g χ [ 1 JB ∂lnρ ∂χ − g χ v 2 s ] (2.16) is a quantity analogous to the square of the Brunt–V ̈ais ̈al ̈a frequency for displacements within the magnetic surface. 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 whichequation (2.10) reduces to ∂ ∂χ (rδB φ ) = ∂ ∂χ [ r 2 J ∂ ∂χ ( ξ φ r ) ] = 0,(2.17) while equation (2.11) becomes vacuous. If the field line crosses the surface Sof 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δBmust be continuous on S, the boundary condition onξ φ is ∂ ∂χ ( ξ φ r ) = 0atχ=χ 1 andχ=χ 2 .(2.18) If, instead, the field line is closed, then periodic boundaryconditions should be applied, withδBbeing 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, wherev s vanishes, and a regularity condition applies toξ χ there. In the absence of rotation, equations (2.13) and (2.14) are two uncoupled, second-order differential equations in Sturm– Liouville form, describing the Alfv ́en continuum and the cusp continuum, respectively, for the polarization for whichthe 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 fieldresults in a coupling between the equations. Similarly, in this case, the equations are coupled by rotation, if the magneticfield is not purely vertical. A second effect of rotation is to modify the effective gravitational accelerationg χ . c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc5 It should be noted, however, that there is a significant difference in the continuous spectrum between the case of no rotation (or uniform rotation) and that of differential rotation. In the former case, one may consider functions of the form (2.3) and take the simultaneous limitsk→ ∞andm→ ∞. A continuous range of polarizations is obtained by choosing different limiting values of the ratiok/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 whichk/m→∞, sincemis assumed finite. When the fluid is in differential rotation, however, the limitm→∞cannot be taken, because an infinite differential Doppler shift would result. The presence of a toroidal magnetic field would also prevent the limitm→∞from being taken. It is clear that the azimuthal wavenumbermappears 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 theeigenfunctions are the same for all values ofm, 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 arbitrarym, provided that attention is restricted to the continuous spectrum. One may proceed by defining{ω 2 n (ψ) :n∈ZZ + }and{u n (ψ,χ) :n∈ZZ + }to be the ordered, real eigenvalues and normalized, real eigenfunctions of the ‘Alfv ́en operator’ that appears in equation (2.13). This is a Sturm–Liouville problem, ∂ ∂χ ( r 2 J ∂u n ∂χ ) +ω 2 n ρr 2 Ju n = 0,(2.19) subject to the boundary condition∂u n /∂χ= 0 atχ=χ 1 andχ=χ 2 (or periodic boundary conditions for a closed field line), and the normalization condition ∫ χ 2 χ 1 ρr 2 Ju m u n dχ=δ mn ,(2.20) for each value ofψ. Let{a n }and{b n }be the components of (siniξ χ /r) and (ξ φ /r) with respect to this set of eigenfunctions, such that sini ξ χ r = ∞ ∑ n=0 a n u n (ψ,χ) and ξ φ r = ∞ ∑ n=0 b n u n (ψ,χ).(2.21) Then the solution of the inhomogeneous equation (2.13) is given by b n =− 2iˆωΩ(ψ) ˆω 2 −ω 2 n (ψ) a n ,(2.22) and exists provided that the appropriate coefficienta n vanishes should ˆω 2 be equal to any of the eigenvalues of the Sturm– Liouville problem. The spectrum of the Alfv ́en operator hasthe 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 eigenfunctionu 0 (ψ,χ) being independent ofχ. This zero-frequency mode is a trivial displacement in which the entire magnetic surface is rotated through an infinitesimal angle about thez-axis. 2 A standard method (Courant & Hilbert 1953) leads to the asymptotic expression ω 2 n (ψ) = n 2 π 2 τ 2 A (ψ) +O(1)(2.23) for the eigenvalues in the limitn→∞, where τ A (ψ) = ∫ χ 2 χ 1 ρ 1/2 Jdχ= ∫ χ 2 χ 1 h χ dχ v A (2.24) is the Alfv ́en time of the magnetic field line. The eigenfunction expansions (2.21) may be substituted into equation (2.14), the equation multiplied byρJξ ∗ χ and integrated with respect toχto yield L[ξ χ ; ˆω 2 ] =R[ξ χ ],(2.25) 2 In the absence of a magnetic field, the infinitesimal rotationof any ring of fluid about the axis is a trivial displacement, but these trivial displacements are easily eliminated from the analysis by discarding the root ˆω= 0. c 1994 RAS, MNRAS000, 000–000 6G. I. Ogilvie where the two functionals L[ξ χ ; ˆω 2 ] = ˆω 2 ∫ χ 2 χ 1 ρ|ξ χ | 2 Jdχ− ∞ ∑ n=1 4ˆω 2 Ω 2 ˆω 2 −ω 2 n |a n | 2 (2.26) and R[ξ χ ] = 4Ω 2 |a 0 | 2 + ∫ χ 2 χ 1 ( v 2 s v 2 s +v 2 A ) B 2 J ∣ ∣ ∣ ∣ ∂ ∂χ ( ξ χ B ) ∣ ∣ ∣ ∣ 2 dχ+ ∫ χ 2 χ 1 {( v 2 s v 2 s +v 2 A ) N 2 χ − 1 JB 2 ∂ ∂χ [( v 2 A v 2 s +v 2 A ) Bg χ ]} ρ|ξ χ | 2 Jdχ (2.27) have been introduced. Note that the infinite sum on the left-hand side starts atn= 1, the termn= 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 (ˆω 2 ) i ∫ χ 2 χ 1 ρ|ξ χ | 2 Jdχ+ (ˆω 2 ) i ∞ ∑ n=1 4ω 2 n Ω 2 |ˆω 2 −ω 2 n | 2 |a n | 2 = 0,(2.28) 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 ∂ ∂ˆω 2 L[ξ χ ; ˆω 2 ] = ∫ χ 2 χ 1 ρ|ξ χ | 2 Jdχ+ ∞ ∑ n=1 4ω 2 n Ω 2 (ˆω 2 −ω 2 n ) 2 |a n | 2 >0,(2.29) 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 thann= 0) for which the coefficienta n does not vanish. The function therefore increases once through all real values in each of the intervals (−∞,p 1 ), (p 1 ,p 2 ), (p 2 ,p 3 ),..., wherep 1 ,p 2 ,p 3 ,...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 beidentical to the true eigenvalues ˆω 2 . Corresponding to each eigenfunctionξ χ is one eigenvalue in each of the intervals (−∞,p 1 ), (p 1 ,p 2 ), (p 2 ,p 3 ),.... 3 Now L[ξ χ ; 0] = 0, which is the reason for transferring the termn= 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 sufficient. Therefore a necessaryand sufficient 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 ( v 2 s v 2 s +v 2 A ) N 2 χ − 1 JB 2 ∂ ∂χ [( v 2 A v 2 s +v 2 A ) Bg χ ] >0(2.30) throughoutχ 1 < χ < χ 2 , for each magnetic surface. 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 ofgravitational potential energy by displacements parallelto the magnetic field. The first term in R[ξ χ ], which may be expressed in the form 4Ω 2 |a 0 | 2 = ( |a 0 | 2 / ∞ ∑ n=0 |a n | 2 ) ∫ χ 2 χ 1 4Ω 2 sin 2 iρ|ξ χ | 2 Jdχ,(2.31) represents the stabilizing effect of rotation, but only on modes for whicha 0 6= 0. In particular, if the disc is symmetrical about the equatorial plane, then this stabilizing effect is zero for all modes with odd symmetry. Significantly, the confinement of the displacements within the magnetic surface implies that no instabilities associated with differential 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. 3 AXISYMMETRIC WAVES AND INSTABILITIES IN THIN DISCS 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. 3 Although these modes share the same displacementξ χ , they have different displacementsξ φ , according to the relation (2.22). c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc7 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 presented here may be seen as an extension of thework 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 unstablemodes 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. 3.1 Review of weakly magnetized discs 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, wherez=H(r) is the location of the upper surface of the disc at radiusr. The internal structure of the disc is resolved by introducing a stretchedvertical coordinateζ=ǫ −1 zwhose value at the surface of the disc isζ s =ǫ −1 H. There are two aspects to obtaining a solution at leading order. First, there are ordinary differential 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. It is assumed that the pressure and density satisfy a polytropic relation on each magnetic surface. Then the leading-order equations for a weakly magnetized disc are Ω 0 = ( GM r 3 ) 1/2 ,(3.1) ∂p 0 ∂ζ =−ρ 0 Ω 2 0 ζ+ 2ρ 0 rΩ 0 Ω 1 B r0 B z0 ,(3.2) ∂B r0 ∂ζ =− 2ρ 0 rΩ 0 Ω 1 B z0 ,(3.3) ∂Ω 1 ∂ζ = 3Ω 0 B r0 2r B z0 (3.4) and p 0 =K 0 ρ Γ 0 ,(3.5) where the angular velocity, density, pressure and magneticfield are given by Ω(r,z) = Ω 0 (r) +ǫΩ 1 (r,ζ) +O(ǫ 2 ),(3.6) ρ(r,z) =ǫ 2s [ ρ 0 (r,ζ) +ǫρ 1 (r,ζ) +O(ǫ 2 ) ] ,(3.7) p(r,z) =ǫ 2s+2 [ p 0 (r,ζ) +ǫp 1 (r,ζ) +O(ǫ 2 ) ] (3.8) and B(r,z) =ǫ s+1 [B r0 (r,ζ) +O(ǫ)]e r +ǫ s+1 [B z0 (r) +O(ǫ)]e z .(3.9) Heresis 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 B r0 (r,0) = 0, ρ 0 (r,ζ s (r)) = 0and B r0 (r,ζ s (r)) =B r0s (r),(3.10) where B r0s (r) denotes the value of B r0 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 ,p 0 and Ω 1 being even functions ofζ, while B r0 is odd. 3.2 Normal modes 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 thecharacteristic frequencies of acoustic and Alfv ́en waves with wavelengths comparable to Hare 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 different scaling be considered. In c 1994 RAS, MNRAS000, 000–000 8G. I. Ogilvie deciding the form of the eigenfunctions of axisymmetric normal modes at leading order, it is natural to assume that they take the form ξ(r,t)∼Re [ξ 0 (r,ζ) exp(−iωt)].(3.11) If this is substituted into equation (2.1), there results atleading order a sixth-order differential system 4 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 onrin general, which means that a global solution, connecting a range of radii, could not be obtained other thanin exceptional circumstances. One is led to conclude that the eigenfunction must depend more strongly onr. The correct form of the eigenfunction is ξ(r,t)∼Re { ξ 0 (r,ζ) exp [ −iωt+ iǫ −1 ∫ k(r) dr ]} ,(3.12) wherek(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. Strictly speaking, equation (3.12) represents a wave propagating in ther-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 ofk, and so it is sufficient to consider travelling waves of this form. When ωis imaginary, equation (3.12) represents an exponentiallygrowing or decaying normal mode rather than a travelling wave. The scaling ofξwithǫis arbitrary, since this is a linear problem, but is chosen tobe O(1) for convenience. For simplicity of notation, the subscript ‘0’ on all quantities is omitted hereafter. When the expression (3.12) is substituted into equation (2.1), the following differential system is obtained at leading order: ρ [ (ω 2 + 3Ω 2 )ξ r −2iωΩξ φ ] = ik δΠ−D(Dξ r −B r ∆),(3.13) ρ(ω 2 ξ φ + 2iωΩξ r ) =−D 2 ξ φ ,(3.14) ρ ( ω 2 + Ω 2 ∂lnρ ∂lnζ ) ξ z = ∂ δΠ ∂ζ −ρΩ 2 ζ∆−D(Dξ z −B z ∆),(3.15) where ∆ = ikξ r + ∂ξ z ∂ζ ,(3.16) δΠ =−(γp+B 2 )∆ +ρΩ 2 ζξ z +B r Dξ r +B z Dξ z (3.17) and D= ik B r +B z ∂ ∂ζ .(3.18) 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 ofk, 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 reflected in the fact, demonstrated below, that the eigenvaluesω of the system (3.13)–(3.15) are either purely real or purelyimaginary whenkis 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 kgoes to zero) or a boundary point (either the inner or outer radius of the disc), and within whichkis real and non-zero. The variation ofξ 0 withris determined by higher-order terms which could be calculated in principle. The equations for an incompressible fluid (γ→∞) 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). 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 B 2 z ∂ 2 ξ ψ ∂ζ 2 =k 2 B 2 z ξ ψ (3.19) 4 Identical to equations (3.13)–(3.15) below, but withk= 0. c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc9 and ( ik B rs +B z ∂ ∂ζ ) 2 ξ φ = 0,(3.20) whereξ ψ is the component ofξin the meridional plane perpendicular to the magnetic field. The appropriate solution of equation (3.19) is ξ ψ ∝exp (−|k|ζ),(3.21) whenkis real. The appropriate solution of equation (3.20) is ξ φ ∝exp [−i(B rs /B z )kζ],(3.22) since this corresponds to a vanishing Eulerian perturbation δB φ = ( ik B rs +B z ∂ ∂ζ ) ξ φ = 0.(3.23) 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 B z ∂ξ r ∂ζ −B rs ∂ξ z ∂ζ =−|k|(B z ξ r −B rs ξ z )(3.24) and B z ∂ξ φ ∂ζ =−ik B rs ξ φ .(3.25) 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 andonly the regular solution is acceptable. 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 ∂ξ r ∂ζ = ∂ξ φ ∂ζ =ξ z = 0(3.26) for an even mode, and ξ r =ξ φ = ∂ξ z ∂ζ = 0(3.27) for an odd mode. Note that this convention is opposite to thatused by LP and KP. The relation between the local dispersion relation and the spectrum of global normal modes is quite straightforward. As in any WKB problem, a quantization condition of the form ǫ −1 ∫ r 2 r 1 k(r) dr=nπ+δ(3.28) applies, wherer 1 andr 2 are the limits of the wave region,nis 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 ofk, 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 sufficient 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 ofk. This is true, at least, for axisymmetric modes with dynamical [i.e. O(1)] growth rates. For completeness, the equations for axisymmetric waves in astrongly magnetized thin disc should be mentioned. In such a disc, the angular velocity and the buoyancy frequencyare both O(1), but the characteristic frequencies of acoustic and Alfv ́en waves with wavelengths comparable to Hare 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 effect, and there is no instability. However, there is a mode with exactly zero frequency in this scaling, c 1994 RAS, MNRAS000, 000–000 10G. I. Ogilvie 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 canbe unstable if the magnetic field provides a significant amount of support against gravity (Agapitou et al. 1997). 4 PROPERTIES OF THE LOCAL DISPERSION RELATION 4.1 Overview 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 principlefor the frequency eigenvalues. This not only shows thatω 2 must be real whenkis real, but also helps to clarify under what circumstances unstable modes are expected. The method hereis an adaptation of the analysis of Papaloizou & Szuszkiewicz (1992) and closely follows the method of Section 2. Before deriving the variational principle it may be noted that the dispersion relationω(k) has reflectional symmetry in bothω- andk-axes. This can be seen from the fact that equations (3.13)–(3.15) have the symmetries (ξ φ ,ω)7→(−ξ φ ,−ω) andeither(ξ,ω,k)7→(ξ ∗ ,−ω,−k), ifωis real,or(ξ,ω,k)7→(ξ ∗ ,ω,−k), ifωis imaginary. When presenting the numerical results it is therefore sufficient to restrict attention to positive values ofω(orω/i) and ofk. 4.2 A variational principle To invert equation (3.14), which is subject to the boundary condition (3.25), one may define the ordered, real eigenvalues {ω 2 n (r) :n∈ZZ + }and the real eigenfunctions{u n (r,ζ) :n∈ZZ + }of the Sturm–Liouville equation B 2 z ∂ 2 u n ∂ζ 2 +ω 2 n ρu n = 0,(4.1) subject to the boundary conditions ∂u n ∂ζ = 0atζ=±ζ s (4.2) and the normalization condition ∫ ζ s −ζ s ρu m u n dζ=δ mn (4.3) for eachr. 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 }are the corresponding squared frequency eigenvalues. The equation D 2 ̃u n +ω 2 n ρ ̃u n = 0,(4.4) subject to the boundary conditions D ̃u n = 0atζ=±ζ s (4.5) is related to the Sturm–Liouville equation by a unitary transformation. It has the same eigenvalues and its eigenfunctions are simply ̃u n (r,ζ) =u n (r,ζ) exp [ −ik ∫ ( B r B z ) dζ ] =u n (r,ζ) exp [ − 2Ω 1 (r,ζ) 3Ω(r) ikr ] .(4.6) Equation (3.14) may be written in the form D 2 ξ φ +ω 2 ρξ φ =−2iωΩρξ r ,(4.7) and is conveniently solved using the eigenfunction expansions ξ r (r,ζ) = ∞ ∑ n=0 a n (r)u n (r,ζ) exp [ − 2Ω 1 (r,ζ) 3Ω(r) ikr ] (4.8) and ξ φ (r,ζ) = ∞ ∑ n=0 b n (r)u n (r,ζ) exp [ − 2Ω 1 (r,ζ) 3Ω(r) ikr ] .(4.9) c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc11 The solution is b n (r) =− 2iωΩ(r) ω 2 −ω 2 n (r) a n (r),(4.10) and exists at any particular value ofrprovided that the appropriate coefficienta n (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)] L[ξ r ,ξ z ;ω 2 ,k] =R[ξ r ,ξ z ;k](4.11) involving the two functionals L[ξ r ,ξ z ;ω 2 ,k] =ω 2 ∫ ζ s −ζ s ρ ( |ξ r | 2 +|ξ z | 2 ) dζ− ∞ ∑ n=1 4ω 2 Ω 2 ω 2 −ω 2 n |a n | 2 (4.12) and R[ξ r ,ξ z ;k] = 4Ω 2 |a 0 | 2 + [ |k||B z ξ r −B r ξ z | 2 ] ζ s −ζ s + ∫ ζ s −ζ s [ |δΠ| 2 γp+B 2 +|Dξ r | 2 +|Dξ z | 2 − ( 1 γp+B 2 ) |B r Dξ r +B z Dξ z +ρΩ 2 ζξ z | 2 −3ρΩ 2 |ξ r | 2 −Ω 2 ζ ∂ρ ∂ζ |ξ z | 2 ] dζ,(4.13) where it has been assumed thatkis real. The imaginary part of equation (4.11) is (ω 2 ) i ∫ ζ s −ζ s ρ ( |ξ r | 2 +|ξ z | 2 ) dζ+ (ω 2 ) i ∞ ∑ n=1 4ω 2 n Ω 2 |ω 2 −ω 2 n | 2 |a n | 2 = 0,(4.14) 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 sufficient 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 ofk. The functional R[ξ r ,ξ z ;k] can be expressed a number of different ways which are useful for some purposes. One alternative is to write R[ξ r ,ξ z ;k] =4Ω 2 |a 0 | 2 + [ |k||B z ξ r −B r ξ z | 2 ] ζ s −ζ s + ∫ ζ s −ζ s { |δΠ| 2 γp+B 2 + 1 B 2 |B z Dξ r −B r Dξ z | 2 + ( γp γp+B 2 ) 1 B 2 ∣ ∣ ∣ ∣ B r Dξ r +B z Dξ z − ( B 2 γp ) ρΩ 2 ζξ z ∣ ∣ ∣ ∣ 2 −3ρΩ 2 |ξ r | 2 − [ Ω 2 ζ ∂ρ ∂ζ + (ρΩ 2 ζ) 2 γp ] |ξ z | 2 } dζ.(4.15) The coefficient of|ξ z | 2 in the integrand is − [ Ω 2 ζ ∂ρ ∂ζ + (ρΩ 2 ζ) 2 γp ] =ρΩ 2 ζ ( 1 γ ∂lnp ∂ζ − ∂lnρ ∂ζ + B r γp ∂B r ∂ζ ) .(4.16) If both the specific entropy and the (radial) magnetic field strength are non-decreasing functions ofζfor 0< ζ < ζ s , then this coefficient 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 differential rotation. A second useful variation is to extend the integral in R[ξ r ,ξ z ;k] over the interval−∞< ζ <∞and write R[ξ r ,ξ z ;k] = 4Ω 2 |a 0 | 2 + ∫ ∞ −∞ [ |δΠ| 2 γp+B 2 +|Dξ r | 2 +|Dξ z | 2 − ( 1 γp+B 2 ) |B r Dξ r +B z Dξ z +ρΩ 2 ζξ z | 2 −3ρΩ 2 |ξ r | 2 −Ω 2 ζ ∂ρ ∂ζ |ξ z | 2 ] dζ,(4.17) 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 . c 1994 RAS, MNRAS000, 000–000 12G. I. Ogilvie Finally, the version for an incompressible fluid is R[ξ r ,ξ z ;k] = 4Ω 2 |a 0 | 2 + [ |k||B z ξ r −B r ξ z | 2 ] ζ s −ζ s + ∫ ζ s −ζ s ( |Dξ r | 2 +|Dξ z | 2 −3ρΩ 2 |ξ r | 2 −Ω 2 ζ ∂ρ ∂ζ |ξ z | 2 ) dζ,(4.18) 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 ofk, which correspond to frequency eigenvaluesω 1 andω 2 , then (ω 2 1 −ω 2 2 ) ∫ ζ s −ζ s ρξ ∗ 1 ·ξ 2 dζ= 0,(4.19) 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 ofk, form a complete set in the space of continuous, vector-valued functions ofζon −ζ s < ζ < ζ s . 4.3 Numerical method 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 B rs and B z 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 B r0 (r,0) = 0 fails to be satisfied defines a mismatch function whose derivative with respect to Ω 1s is obtained by simultaneously integrating the equations differentiated 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 coefficients in these equations areknown only numerically, and so to avoid the use of interpolation one must integrate equations (3.1)–(3.5) simultaneously (with the previously determined eigenvalue). Whenkis 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. The most difficult part of the numerical method is the treatment of the singular point atζ=ζ s . One can re-write equations (3.13) and (3.15) in the form ∂ 2 ξ r ∂ζ 2 = N r γp ,(4.20) ∂ 2 ξ z ∂ζ 2 = N z γp ,(4.21) where the numerators N r and N z are regular and are linear functions ofξ r ,ξ φ ,ξ z and their first derivatives. The condition that N r and N z both go to zero as fast asγpasζ→ζ 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 involvesmuch tedious algebra. Once a mode has been computed for some value ofk, it can be followed quasi-continuously askis 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. All numerical calculations are carried out in units such that GM=K 0 =r= 1, andǫis defined such thatζ s = 1 at the value ofrunder consideration. In particular, this means that B 0 is expressed in units of (GM) Γ/2(Γ−1) K −1/2(Γ−1) 0 r −3Γ/2(Γ−1) ζ Γ/(Γ−1) s .(4.22) Taking into account the factors ofǫ, this means that the true unit of magnetic field strength is (GM) Γ/2(Γ−1) K −1/2(Γ−1) r −3Γ/2(Γ−1) H Γ/(Γ−1) .(4.23) Similarly, frequencies and growth rates are expressed in units of the local Keplerian angular velocity, and the radial wavenumber in units of H −1 . c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc13 5 CLASSIFICATION OF MODES 5.1 Hydrodynamic discs revisited 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 ρ 0 = [( Γ−1 Γ ) h 0 K 0 ] 1/(Γ−1) ,(5.1) p 0 =K 0 [( Γ−1 Γ ) h 0 K 0 ] Γ/(Γ−1) ,(5.2) where the enthalpy is h 0 = 1 2 Ω 2 0 (ζ 2 s −ζ 2 ).(5.3) 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 gmodes (gravity modes, due to buoyancy and driven mainly by gravitational forces), both withω 2 ≥Ω 2 , and also r modes 5 (inertial modes, due to rotation and driven mainly by inertial forces), withω 2 ≤Ω 2 . Examples of the eigenfunctions for these modes can be foundin 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 toclarify the situation and demonstrate that there are in facttwo f modes. 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 different parameters. The different 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 modes have non-trivial limits; (iii) in the limitγ→ ∞, in which the fluid becomes incompressible, the p modes haveω 2 → ∞, while all other modes have non-trivial limits. 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 limitkζ 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 parameterkζ s . It is found numerically that all modes withω 2 >Ω 2 are trapped in a layer of extent O ( (kζ s ) −1 ) 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’), sinceg=O(1) in the layer; (ii) p modes as acoustic modes (‘ω 2 ∼v 2 s k 2 ’), sincev 2 s =O ( (kζ s ) −1 ) in the layer; (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 limitkζ 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 confluent 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 largekζ s in an accretion disc. To prove this, one should introduce a stretched coordinate x 1 =k(ζ s −ζ)(5.4) which is O(1) in the layer in which the modes are trapped. Equations (3.13)–(3.15) are then solved asymptotically with ξ r (r,ζ) =ξ r0 (r,x 1 ) +O ( (kζ s ) −1 ) ,(5.5) ξ z (r,ζ) =ξ z0 (r,x 1 ) +O ( (kζ s ) −1 ) (5.6) and ω 2 =λΩ 2 kζ s +O(1),(5.7) 5 Called g modes by LP. c 1994 RAS, MNRAS000, 000–000 14G. 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 r o 1 , r e 1 , r o 2 , r e 2 , r o 3 , r e 3 , r o 4 , r e 4 , r o 5 , r e 5 , g e 5 , g o 5 , g e 4 , g o 4 , g e 3 , g o 3 , g e 2 , g o 2 , g e 1 , g o 1 , f e , f o , p e 1 , p o 1 , p e 2 , p o 2 , p e 3 , p o 3 , p e 4 , p o 4 , p e 5 , p o 5 , although even and odd f, p and g modes coalesce for sufficiently largek. whileξ φ =O ( (kζ s ) −1/2 ) . The quantityλis the dimensionless eigenvalue of the leading-order equations, which, following Lamb (1932), can be combined to obtain x 1 ∂ 2 ∆ 0 ∂x 2 1 +c ∂∆ 0 ∂x 1 + [2(n−1) +c−x 1 ] ∆ 0 = 0,(5.8) where ∆ 0 = iξ r0 − ∂ξ z0 ∂x 1 (5.9) is the leading-order part of∇·ξ, whilecandnare given by c= 2Γ−1 Γ−1 (5.10) and Γλ 2 −γ[2(Γ−1)n+ 1]λ+ (γ−Γ) = 0.(5.11) The solution regular atx 1 = 0 is (e.g. Erd ́elyi et al. 1953a) ∆ 0 ∝e −x 1 1 F 1 (1−n;c; 2x 1 ).(5.12) It is exponentially small asx 1 →+∞, thereby justifying the layer analysis, if and only ifnis a positive integer, in which case the confluent hypergeometric function is proportional to a generalized Laguerre polynomial of degreen−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 ... < λ − 2 < λ − 1 <1< λ + 1 < λ + 2 < ...,(5.13) with λ + n ∼ 2γ(Γ−1) Γ nandλ − n ∼ γ−Γ 2γ(Γ−1) 1 n (5.14) asn→ ∞. The modes corresponding toλ + n andλ − n are called p n and g n respectively in stellar oscillations. There is an additional mode which corresponds to the trivial solution ∆ 0 = 0 of equation (5.8). This has ξ r0 ∝ξ z0 ∝e −x 1 (5.15) andλ=λ 0 = 1. It is a surface gravity mode and corresponds to the f mode of stellar oscillations. c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc15 Somewhat surprisingly, the p modes and g modes are equally compressive, sharing the same function ∆ 0 , and only the f mode has∇·ξ= 0 at leading order. Nevertheless, the g modes do survive theincompressible limit. Whenγ→∞, λ + n →∞whileλ − n → 1 2(Γ−1)n+ 1 .(5.16) Conversely, whenγ→Γ, λ + n →2(Γ−1)n+ 1whileλ − n →0.(5.17) 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 limitkζ s →∞, as described above. For general values ofk, one should refer to the modes as f e , f o , p e n , p o n , g e n and g o n , wheren∈IN is a vertical mode number and the superfix denotes the parity of the mode. 6 The r modes behave quite differently in the limitkζ s →∞. In a convectively stable disc, the r modes become trapped in a layer of extent O ( (kζ s ) −1/2 ) centred on the equatorial plane, and haveω 2 =O ( (kζ s ) −1 ) . The asymptotic analysis proceeds in a similar way, with a stretched coordinate x 2 = (kζ s ) 1/2 (ζ/ζ s )(5.18) and scalings ξ r (r,ζ)∼(kζ s ) −1/2 ξ r0 (r,x 2 ),(5.19) ξ φ (r,ζ)∼ξ φ0 (r,x 2 ),(5.20) ξ z (r,ζ)∼ξ z0 (r,x 2 )(5.21) and ω 2 ∼(kζ s ) −1 Ω 2 λ.(5.22) The leading-order equations are iξ r0 + ∂ξ z0 ∂x 2 = 0,(5.23) ξ φ0 =−2iλ −1/2 ξ r0 (5.24) and ∂ 2 ξ z0 ∂x 2 2 + [ λ− 2(γ−Γ) γ(Γ−1) x 2 2 ] ξ z0 = 0.(5.25) A bounded solution exists if and only if λ= (2N+ 1) √ 2(γ−Γ) γ(Γ−1) ,(5.26) where Nis a non-negative integer; the solution is then ξ z0 ∝exp(− 1 2 α 2 x 2 2 )H N (αx 2 ),(5.27) where α= [ 2(γ−Γ) γ(Γ−1) ] 1/4 ,(5.28) and H N is an Hermite polynomial (e.g. Erd ́elyi et al. 1953b). If Nis odd, the mode is even and may be called r e n , where N= 2n−1. If Nis even, the mode is odd and may be called r o 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 ( (kζ s ) −2 ) and the eigenfunctions are not localized. 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 ofkζ s and may not accurately reflect the properties of the modes for smallerkζ s . 6 An alternative classification scheme is possible in which p e n and p o n are renamed p 2n and p 2n−1 , respectively, and similarly for other classes of modes. This is closer to the notation of KP, and hasthe advantage that the superfix is not required. However, thefrequency eigenvalues do not then fall in sequence for all classes of modes. Also, to refer to f 1 and f 2 rather than f o and f e is perhaps less helpful since it obscures the uniqueness of the f mode. c 1994 RAS, MNRAS000, 000–000 16G. I. Ogilvie Figure 2.Asymptotic behaviour of the local dispersion relation for large radial wavenumberk, 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) ask→ ∞, while the r modes (bottom) haveω 2 =O(k −1 ). In each case the dotted lines indicate the asymptotic limitsderived in the text. Note that the vertical axis is different in each plot. 5.2 Magnetized discs 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 difficult, although there is one limit which can be investigated semi-analytically. This is the casek= 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 forcevanishes and equations (5.1)–(5.3) for the equilibrium are valid. If alsok= 0, equations (3.13)–(3.15) simplify considerably because the horizontal and vertical components ofξbecome decoupled. There exist purely horizontal modes, with ξ r (r,ζ) =a r (r)u N (r,ζ)(5.29) and ξ φ (r,ζ) =a φ (r)u N (r,ζ),(5.30) c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc17 whereu N is one of the eigenfunctions of equation (4.1), and the componentsa r anda φ satisfy the equation [ ω 2 + 3Ω 2 −ω 2 N −2iωΩ 2iωΩω 2 −ω 2 N ][ a r a φ ] = [ 0 0 ] .(5.31) The frequency eigenvalues are the roots of the equation ω 4 −(2ω 2 N + Ω 2 )ω 2 +ω 2 N (ω 2 N −3Ω 2 ) = 0,(5.32) viz. ω 2 =ω 2 N + 1 2 Ω 2 [ 1± ( 1 + 16ω 2 N /Ω 2 ) 1/2 ] ,(5.33) so that there is an unstable mode if and only if 0< ω 2 N <3Ω 2 . Since the eigenvalues{ω 2 N }are ordered, the N= 1 mode is unstable if any mode is unstable. The maximum possible growth rate is the Oort parameter A= 3 4 Ω, which is achieved if ω 2 N = 15 16 Ω 2 . In the limit B z →0, the eigenvaluesω 2 N →0 (although not uniformly with respect to N), and equation (5.33) reduces to ω 2 = Ω 2 or 0.(5.34) 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 m e n if N= 2n, or m o n if N= 2n−1. If N= 0 the mode is trivial. The greater solution of equation (5.33) may be called r e n if N= 2n, r o n if N= 2n−1, or f e if N= 0. However, the r and g modes really lose their identity whena magnetic field is introduced. There are also purely vertical modes, which satisfy the equation (ζ 2 s −ζ 2 ) ∂ 2 ξ z ∂ζ 2 −(1 + 2q)ζ ∂ξ z ∂ζ +N(N+ 2q)ξ z = 0,(5.35) where q= Γ + 1 2(Γ−1) (5.36) and N(N+ 2q) = 1 γ ( 2Γ Γ−1 ) ( ω 2 −Ω 2 Ω 2 ) ,(5.37) Appropriate solutions, regular atζ=±ζ s , are obtained when Nis a non-negative integer, in which case ξ z ∝C q N (ζ/ζ s ),(5.38) where C q N is a Gegenbauer polynomial (e.g. Erd ́elyi et al. 1953b). Equation (5.37) then leads to the dispersion relation ω 2 /Ω 2 = 1 +γ ( Γ−1 2Γ ) N [ N+ ( Γ + 1 Γ−1 )] .(5.39) These modes may be identified as p e n if N= 2n−1, p o n if N= 2n, or f o if N= 0. 7 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 sufficiently smallk, and their growth rates are greatest fork= 0. The stable modes form a complicated dispersion diagram because of the avoided crossings that occur between different 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 flat. 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 withk= 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 reflection in a 7 In fact the f o and r o 1 modes in a Keplerian disc without a magnetic field involve both horizontal and vertical displacements even when k= 0. c 1994 RAS, MNRAS000, 000–000 18G. 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 B z = 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 atk= 0 by the m o 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 ask→∞. 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 B z = 0.002. Modes with even and odd symmetry are shown in the left andright panels, respectively. In this case seven unstable magnetorotational modes of each parity exist. A maximum growth rate of 0.7495Ω is achieved atk= 0 by the m e 4 mode. The frequency eigenvalues of the continuous spectrum, which are all real, are marked on the scale at the right of each panel. horizontal plane) have been analysed by Hasan & Christensen-Dalsgaard (1992). 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 reflected in the fact that the convergence of the frequency eigenvalues of m n modes to zero as B→0is not uniform with respect ton. Modes of largenclutter the dispersion diagram even when the magnetic field is exceedingly weak, and can be unstable with large growth rates. Finally, an example of the dispersion diagram for a disc containing a bending poloidal magnetic field is shown in Fig. 5. c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc19 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,B z = 0.015 and B rs = 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 m o 1 mode is unstable, achieving a maximum growth rate of 0.6019Ω atk= 0. The frequency eigenvalues of the continuous spectrum, which are all real, are marked onthe scale at the right of each panel. The avoided crossings are wider than for a purely vertical magnetic field, which implies that the coupling between different branches of modes is stronger when the field lines bend. The behaviour of branches in the limitk→∞is quite different 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. 6 STABILITY CRITERIA FOR MAGNETIZED DISCS The stability of a weakly magnetized thin disc to axisymmetric perturbations can be decided, in principle, by computingthe local dispersion relation at each radius separately; the disc is unstable if and only if an imaginary frequency eigenvalue exists, for any real value ofk, at any radius. However, a more efficient algorithm is required in practice. In the case of polytropic discs, the parameter space divides into stable and unstableregions separated by a curve of marginal stability. The location of the marginal curve depends on a further parameter which isthe adiabatic exponentγ. Equilibria on the marginal curve possess a mode withω= 0, but this is not sufficient to determine the curve, since allunstable equilibria also possess modes withω= 0 (cf. Figs 3–5). In principle, a marginal curve can be drawnfor each m mode and for every value ofk. Before describing the marginal curves, it is appropriate togive a more detailed account of the solutions of the equilibrium equations for polytropic discs. When B rs and B z 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 (B rs ,B z ,Ω 1s ). Some of these solutions are physically acceptable and maybe 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 asthe ‘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. Theorem 1. All irregular equilibria possess an unstable mode of odd symmetry atk= 0. Proof. Whenk= 0, the form of the functional R[ξ r ,ξ z ;k] for purely horizontal trial displacements of odd symmetryis, from equation (4.13), R[ξ r ,0; 0] = ∫ ζ s −ζ s ( B 2 z ∣ ∣ ∣ ∣ ∂ξ r ∂ζ ∣ ∣ ∣ ∣ 2 −3ρΩ 2 |ξ r | 2 ) dζ,(6.1) c 1994 RAS, MNRAS000, 000–000 20G. I. Ogilvie and the trial function must satisfy∂ξ r /∂ζ= 0 atζ=±ζ s . Given that the equilibrium is irregular,B r 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 B 2 z ∂ 2 y ∂ζ 2 +λρy= 0(6.2) on−ζ 1 < ζ < ζ 1 , subject to the boundary conditions ∂y ∂ζ = 0atζ=±ζ 1 .(6.3) Since the equation is symmetrical aboutζ= 0, the eigenfunctionsy 0 ,y 1 ,y 2 ,...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 toy 0 = constant. Now it follows from equations (3.1)–(3.4) for the equilibrium that Ω 1 satisfies the equation B 2 z ∂ 2 Ω 1 ∂ζ 2 + 3ρΩ 2 Ω 1 = 0,(6.4) and also the boundary conditions (6.3). Moreover, sinceζ 1 is the first zero of B r in 0< ζ < ζ s ,B r has exactly two extrema, and Ω 1 has exactly two zeros, in the domain of the Sturm–Liouville equation. It follows that Ω 1 is the eigenfunctiony 2 , and has eigenvalueλ 2 = 3Ω 2 . The properties of the Sturm–Liouville equation then ensure that the first odd eigenfunctiony 1 has an eigenvalue satisfying the (strict) inequalities 0< λ 1 <3Ω 2 . Now construct the trial function f(ζ) = { y 1 (ζ),|ζ|< ζ 1 , sgn(ζ)y 1 (ζ 1 ),|ζ|≥ζ 1 , (6.5) which is continuous throughout−ζ s < ζ < ζ s and satisfies the boundary conditions for the variational principle. After an integration by parts one obtains R[f,0; 0] =−2 ∫ ζ 1 0 (3ρΩ 2 −λ 1 )|f| 2 dζ−2 ∫ ζ s ζ 1 3ρΩ 2 |f| 2 dζ,(6.6) and, since this is clearly negative, the theorem is proved. In Fig. 6, the different 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 B z 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 B z , approximately 0.01283, below which the principal branch becomes separated from the‘vertical’ solution (which has B r = 0 and Ω 1s = 0) by a stretch of irregular equilibria. (The vertical solution, however,is technically always regular.) As B z is reduced further, a bifurcation disconnects the principal branch from the vertical solution. An infinite number of further bifurcations occur as B z →0. It is observed in the numerical calculations that, for regular equilibria, the last mode to be stabilized is the m o 1 mode, and it is stabilized last atk= 0. This supports the following conjecture. Conjecture. If a regular equilibrium possesses an unstable mode at any value ofk, then it possesses an unstable mode at k= 0, assuming thatγ≥Γ. 8 While this conjecture has not yet been proved in the most general circumstances envisaged, the following two theorems support it strongly. Theorem 2. The conjecture holds in the case of an incompressible fluid,provided there is no density inversion. Proof. An appropriate version of the functional R[ξ r ,ξ z ;k] for an incompressible fluid is R[ξ r ,ξ z ;k] = 4Ω 2 |a 0 | 2 + ∫ ∞ −∞ ( |Dξ r | 2 +|Dξ z | 2 −3ρΩ 2 |ξ r | 2 −Ω 2 ζ ∂ρ ∂ζ |ξ z | 2 ) dζ,(6.7) where the trial displacementξis defined on−∞< ζ <∞and is subject to the constraint ∆ = 0. Given that there existsan unstable eigenfunction ˆ ξ(ζ) at some wavenumber ˆ k6= 0, it follows that R[ ˆ ξ r , ˆ ξ z ; ˆ k]<0.(6.8) Whenk= 0, the constraint ∆ = 0 can be satisfied by using a purely horizontal trial displacementξ r =f, where the function f(ζ) is defined by ˆ ξ r (ζ) =f(ζ) exp [ −(2Ω 1 /3Ω)i ˆ kr ] ,(6.9) 8 Whenγ <Γ there may be a (magneto-) convective instability and the result cannot be expected to hold. In fact, a condition differing slightly fromγ≥Γ may be required in order to prove the conjecture. c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc21 Figure 6.Solution curves for weakly magnetized, polytropic thin discs with Γ = 4/3. The ‘eigenvalue’ Ω 1s is plotted against B rs for various values of B z (as indicated in each panel). All solutions in the range−0.025< B rs <0.025 are shown. Regular and irregular equilibria are indicated by solid and dotted lines, respectively. As B z 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. so that D ( ˆ k) ˆ ξ r = i ˆ k B r ˆ ξ r +B z ∂ ˆ ξ r ∂ζ =B z ( ∂f ∂ζ ) exp [ −(2Ω 1 /3Ω)i ˆ kr ] = ( D (0) f ) exp [ −(2Ω 1 /3Ω)i ˆ kr ] ,(6.10) c 1994 RAS, MNRAS000, 000–000 22G. I. Ogilvie where the suffix on the operator Dindicates the value ofkused. Then | ˆ ξ r | 2 =|f| 2 and|D ( ˆ k) ˆ ξ r | 2 =|D (0) f| 2 .(6.11) Also, the coefficients{a n }in the eigenfunction expansion (4.8) are the same whenξ r =fandk= 0 as they are whenξ r = ˆ ξ r andk= ˆ k. Thus the difference R[ ˆ ξ r , ˆ ξ z ; ˆ k]−R[f,0; 0] = ∫ ∞ −∞ ( |D ( ˆ k) ˆ ξ z | 2 −Ω 2 ζ ∂ρ ∂ζ | ˆ ξ z | 2 ) dζ(6.12) is non-negative, and so R[f,0; 0]≤R[ ˆ ξ r , ˆ ξ z ; ˆ k]<0,(6.13) which proves that an unstable mode must exist withk= 0. Theorem 3. The conjecture holds for a compressible fluid when the magnetic field is purely vertical. Proof. In this case one may take R[ξ r ,ξ z ;k] =4Ω 2 |a 0 | 2 + ∫ ∞ −∞ { |δΠ| 2 γp+B 2 +B 2 ∣ ∣ ∣ ∣ ∂ξ r ∂ζ ∣ ∣ ∣ ∣ 2 + ( γp γp+B 2 ) B 2 ∣ ∣ ∣ ∣ ∂ξ z ∂ζ − ( ρΩ 2 ζ γp ) ξ z ∣ ∣ ∣ ∣ 2 −3ρΩ 2 |ξ r | 2 − [ Ω 2 ζ ∂ρ ∂ζ + (ρΩ 2 ζ) 2 γp ] |ξ z | 2 } dζ,(6.14) where δΠ =−(γp+B 2 )∆ +ρΩ 2 ζξ z +B 2 ∂ξ z ∂ζ .(6.15) Given that there exists an unstable eigenfunction ˆ ξ(ζ) at some wavenumber ˆ k6= 0, it follows that R[ ˆ ξ r , ˆ ξ z ; ˆ k]<0.(6.16) Whenk= 0, the purely horizontal displacementξ r = ˆ ξ r should be used, andδΠ then vanishes. The coefficient of|ξ z | 2 in the integrand of equation (6.14), − [ Ω 2 ζ ∂ρ ∂ζ + (ρΩ 2 ζ) 2 γp ] = ( γ−Γ Γ ) (ρΩ 2 ζ) 2 γp ,(6.17) is non-negative, since it is assumed thatγ≥Γ. Thus R[ ˆ ξ r ,0; 0]≤R[ ˆ ξ r , ˆ ξ z ; ˆ k]<0, which proves that an unstable mode must exist withk= 0. 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 m o 1 mode) atk= 0, drawn on the principal branch. This curve has been computed directly by solving theequations for an equilibrium (with a given value of B z ) and a mode (withω= 0 andk= 0) as an eigenvalue problem, yielding a value for B rs . 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/B 2 . The plasma beta of the marginal equilibria, evaluated on the equatorial plane, is plotted against B z in Fig. 8. It is seen that all marginal equilibria have valuesofβ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 effect of the adiabatic exponentγon stability could not be taken into account. The magnetorotational modes are not affected greatly by buoyancy or compressibility. In Fig. 9 the stability boundaries for equilibria with Γ = 4/3 and Γ = 5/3 are plotted for different values of the adiabatic exponent. The unstable region reduces c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc23 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 inclinationi=π/6. The dashed line is the critical curve for the existence of the principal solution branch. The solid line is the marginal curve, atk= 0, for the m o 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 different scales used in each plot. Figure 8.The plasma beta of the marginal equilibria, evaluated on theequatorial 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 effect of a sub-adiabatic stratification. The effect is most pronounced for equilibria in which the magnetic field bends significantly, and is zero for those with a purely vertical magnetic field. 7 NON-AXISYMMETRIC MODES The analysis of Section 3 may be extended to include non-axisymmetric waves and instabilities of azimuthal wavenumberm. A distinction must be made between modes withm=O(ǫ −1 ) and those withm=O(1). In the former case, the azimuthal c 1994 RAS, MNRAS000, 000–000 24G. I. Ogilvie Figure 9.The effect of sub-adiabatic stratification on the magnetorotational instability. Left: stability boundaries for weakly magnetized, polytropic thin discs with Γ = 4/3 are plotted for four different 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 equilibriawith Γ = 5/3. The three solid lines are the stability boundaries forγ= 5/3, 3 and∞(from right to left). 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 differential 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 thedisc to be present in this small region if they are to have dynamical [O(1)] growth rates. For modes withm=O(1), however, the equations of Section 3 require very few changes. The eigenfunctions have the form ξ(r,t)∼Re { ξ 0 (r,ζ) exp [ −iωt+ imφ+ iǫ −1 ∫ k(r) dr ]} ,(7.1) 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 isothermalaccretion disc. Provided that the dispersion relation yields real values of ˆωfor real values ofk, there is no difficulty becauseω, ˆωandkare all real within the wave region of a mode. There remains the possibility that the dispersion relationyields imaginary values ofωfor real values ofk, 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 ofk. The conclusion is that, for a non-axisymmetric unstable mode,kis real only at the corotation radius. Whenkis not real, the variational principle of Section 4.2 is not valid. Otherwise, the only change that need be made in p(k) = { +k,Re(k)>0, −k,Re(k)<0, (7.2) 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∂ω/∂kis known for real values ofkfrom the real dispersion relation, and is imaginary for an unstable branch. Letk=k c be the real value of the radial wavenumber at the corotation radiusr=r c of such a mode. Then, at a neighbouring radiusr, the wavenumber is (to first order) k≈k c + [ 3mΩ c 2r c / ( ∂ω ∂k ) c ] (r−r c ),(7.3) 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−r c ). If it is an increasing c 1994 RAS, MNRAS000, 000–000 Waves and instabilities in a differentially rotating disc25 function, then the mode will take the form of a localized Gaussian wave packet aroundr=r c , according to equation (7.1). If it is a decreasing function, the function (7.1) would increase greatly in magnitude away fromr=r c 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 ofmdepending 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 thiskind can be found. 8 DISCUSSION 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 inthe continuous spectrum. Neither is the magnetorotationalinstability or any other instability associated with differential 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 largeradial wavenumber. When a magnetic field is introduced, the dispersion diagram is complicated by a large number of avoided crossings which make a systematic classification difficult. If the magnetic field is relatively weak, unstable magnetorotational modes also occur for sufficiently 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 influence. A sub-adiabatic stratification also has a mild stabilizing influence if the magnetic field is not purely vertical. 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 ofthe disc are inclined to the vertical at angles significantlygreater 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. 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 differential 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 flows but could be subject to weak turbulence or at least fluctuations. However, that would be quite different 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 magnetizeddisc should be constructed by choosing the functional formsof ζ s (r),ψ 0 (r) and K 0 (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 tobe explored are the possible stabilization of a super-adiabatically stratified disc by a sufficiently strong magnetic field, the propagation of them= 1 ‘tilt’ mode in a magnetized disc, and the effects of self-gravitation on the equilibria and spectra ofboth hydrodynamic and magnetized discs. ACKNOWLEDGMENTS I would like to thank Jim Pringle, Douglas Gough and Ulf Torkelsson for helpful discussions. A research studentship fromthe Particle Physics and Astronomy Research Council is acknowledged. REFERENCES Agapitou V., Papaloizou J. C. B., Terquem C., 1997, MNRAS, 292, 631 c 1994 RAS, MNRAS000, 000–000 26G. I. Ogilvie Balbus S. A., Hawley J. F., 1991, Ap J, 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, Ap J, 446, 741 Brandenburg A., Nordlund ̊ A., Stein R. F., Torkelsson U., 1996, Ap J, 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. Mc Graw-Hill,New York Erd ́elyi A., Magnus W., Oberhettinger F., Tricomi F. G., 1953b, Higher Transcendental Functions, vol. 2. Mc Graw-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, Ap J, 396, 311 Hawley J. F., Gammie C. F., Balbus S. A., 1995, Ap J, 440, 742 Hawley J. F., Gammie C. F., Balbus S. A., 1996, Ap J, 464, 690 Heyvaerts J. F., Norman C., 1989, Ap J, 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, Ap J, 416, 689 Lubow S. H., Pringle J. E., 1993, Ap J, 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, Ap J, 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, Ap J, 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 T E X macros. c 1994 RAS, MNRAS000, 000–000