A&A 444, 337-346 (2005)
DOI: 10.1051/0004-6361:20042591

Magneto-rotational overstability in accretion disks

J. W. S. Blokland^1,2 - E. van der Swaluw¹ - R. Keppens^1,2 - J. P. Goedbloed^1,2

1 - FOM-Institute for Plasma Physics Rijnhuizen, PO Box 1207, 3430 BE Nieuwegein, The Netherlands
2 - Association EURATOM-FOM, Trilateral Euregio Cluster

Received 22 December 2004 / Accepted 13 April 2005

Abstract
We present analytical and numerical studies of magnetorotational instabilities occuring in magnetized accretion disks. These calculations are performed for general radially stratified disks in the cylindrical limit. We elaborate on earlier analytical results and confirm and expand them with numerical computations of unstable eigenmodes of the full set of linearised compressible MHD equations. We compare these solutions with those found from approximate local dispersion equations from WKB analysis.

In particular, we investigate the influence of a nonvanishing toroidal magnetic field component on the growth rate and oscillation frequency of magnetorotational instabilities in Keplerian disks. These calculations are performed for a constant axial magnetic field strength. We find the persistence of these instabilities in accretion disks close to equipartition. Our calculations show that these eigenmodes become overstable (complex eigenvalue), due to the presence of a toroidal magnetic field component, while their growth rate reduces slightly.

Furthermore, we demonstrate the presence of magneto-rotational overstabilities in weakly magnetized sub-Keplerian rotating disks. We show that the growth rate scales with the rotation frequency of the disk. These eigenmodes also have a nonzero oscillation frequency, due to the presence of the dominant toroidal magnetic field component. The overstable character of the MRI increases as the rotation frequency of the disk decreases.

Key words: accretion, accretion disks - instabilities - magnetohydrodynamics (MHD) - plasmas

1 Introduction

Accretion disks are a rather common phenomenon in astrophysics. They vary in size from a few up to a hundred astronomical units (AU) in young stellar objects (YSO) to approximately a hundred parsecs (pc) in the centers of active galactic nuclei (AGN). The central accreting object is a protostar in YSOs, a white dwarf, neutron star or a black hole in a binary, and a massive black hole in the center of an AGN. Part of the gravitational energy released during the accretion process is radiated away over a wide range of frequencies. Multi-wavelength studies obtained from modern observational facilities have significantly increased our knowledge about the nature of accretion disks and their associated central objects.

On the other hand, insights into the nature of the accretion process itself have largely been obtained by theoretical and computational studies. In all the above mentioned astrophysical objects, one needs outward transport of angular momentum of the accreting material in order to sustain an accretion disk around the central object. This angular momentum transport can be sustained by a turbulent viscosity mechanism operating inside the disk material itself, where this mechanism excerts a torque on the accretion disk (Shakura & Sunyaev 1973). This turbulent viscosity, in turn, can originate from the development of fluid or magnetofluid instabilities occuring in the accretion disk. In the early 1990s it was realised by Balbus & Hawley (1991) that magneto-rotational instability (MRI) could provide the physical basis for this angular momentum transport in accretion disks. This instability was already known in the literature (Velikhov 1959 and Chandrasekhar 1960), but had not been applied in the context of accretion disks.

In more recent years, global magnetohydrodynamical (MHD) simulations of accretion disks have been performed, where much of the dynamics is interpreted as a direct consequence of the presence of the magneto-rotational instability (see for example Hawley et al. 2001). Other papers rather claim that both convective and magneto-rotational instabilities can play a dominant role in the transport of angular momentum (see for example Igumenshchev et al. 2003). These types of disks are referred to as convection-dominated accretion flows (CDAF). Less attention has been paid in recent years to the spectral analysis of instabilities occurring in accretion disks (see however Christodoulou et al. 2003).

The aim of this paper is to present a detailed linear analysis of magnetorotational instabilities present in a variety of global accretion disk configurations. We will consider a sample of accretion disk configurations that vary from sub-Keplerian (thick) to Keplerian (thin) rotating disks. In the case of Keplerian rotating disks, we will consider both weakly magnetized disks and disks that are close to equipartition. In particular, we investigate the influence of a nonvanishing toroidal magnetic field component on the growth rate and oscillation frequency of magnetorotational instabilities in Keplerian disks. To our knowledge, this has not extensively been done for compressible plasmas.

In this paper, we will limit ourselves to axisymmetric instabilities and the configurations are all taken in the cylindrical limit. The calculations continue the work presented by Keppens et al. (2002), advocating the need for a more elaborate magnetohydrodynamic spectroscopic analysis of all waves and instabilities in magnetized disks. We will use a semi-analytical approach, as well as numerical solutions to the full set of linearised, compressible MHD equations obtained with the code LEDAFLOW (Nijboer et al. 1997).

We will first show that magneto-rotational instabilities are present in both sub-Keplerian and Keplerian rotating disks and compare the growth rates obtained for these models. We find the presence of MRI even for cases where the disk is close to equipartition (i.e. $\beta\sim 1$), however the strength of the toroidal magnetic field should then be much larger than the axial magnetic field strength.

This paper is organised as follows: in Sect. 2, we recall the essential elements from spectral theory of MHD waves and instabilites. In Sect. 3, we present the model and the limitations of the accretion disk configuration we use. In Sect. 4, we explain the numerical strategies. In Sect. 5, we discuss the magnetorotational overstabilities which occur in our models and finally, in Sect. 6, we summarise and present our conclusions.

2 Spectral theory

We make use of the ideal MHD approximation to model an accretion disk in the presence of a magnetic field. This assumes that the disk matter consists of a plasma that is sufficiently ionized and that one treats the dynamics on a length scale such that the one-fluid approximation is valid. The ideal MHD equations are

    

$$\rho \frac{\partial {\vec{v}}}{\partial t} = -\rho{\vec{v}}\cdot\nabla {\vec{v}} - \nabla p + {\vec{j}} \times {\vec{B}} - \rho\nabla \Phi,$$ (1)

$$\frac{\partial p}{\partial t} = -{\vec{v}}\cdot\nabla p - \gamma p\nabla \cdot {\vec{v}},$$ (2)

$$\frac{\partial {\vec{B}}}{\partial t} = \nabla \times ({\vec{v}} \times {\vec{B}}),$$ (3)

$$\frac{\partial \rho}{\partial t} = -\nabla \cdot (\rho{\vec{v}}),$$ (4)

where the variables, $\rho$, p, ${\vec{v}}$, ${\vec{B}}$, $\Phi$, ${\vec{j}} = \nabla \times {\vec{B}}$ and $\gamma$are the density, pressure, velocity, magnetic field, gravitional potential, electric current and ratio of the specific heats, respectively. Furthermore, from Maxwell's theory, the equation  $\nabla \cdot {\vec{B}} = 0$ must also be satisfied. We have assumed a dimensionalization where the permeability of vacuum $\mu_{0}=1$. Equations (1)-(4) are the momentum, entropy, induction and mass conservation equation, respectively.

2.1 Frieman-Rotenberg formalism

To investigate the stability properties of accretion disks, we linearize the ideal MHD Eqs. (1)-(4). The linearisation is done by assuming a time-independent equilibrium. This assumption is justified at time scales which are much shorter than the dynamical time scale of the disk equilibrium. When linearizing, we write all variables of the ideal MHD equations in the following way:

$$f = f_{0}({\vec{r}}) + f_{1}({\vec{r}},t), \qquad \textrm{and} \ \vert f_{1}\vert \ll f_{0}.$$ (5)

Here, f₀ is the equilibrium quantity while f₁ is the time-dependent fluctuation about the equilibrium quantity. The resulting equations can be rewritten in terms of the Lagrangian displacement field,  $\vec{\xi}({\vec{r}},t)$, which is related to the perturbed velocity,  ${\vec{v}}_{1}$, in the following way

$${\vec{v}}_{1} = \left( \frac{\partial }{\partial t} + {\vec... ...t\nabla \right)\vec{\xi} - \vec{\xi}\cdot\nabla {\vec{v}}_{0},$$ (6)

where  ${\vec{v}}_{0}$ represents the equilibrium flow velocity. By introducing the Lagrangian displacement field, there is no longer confusion between equilibrium and perturbed quantities when we suppress the subscript 0 on all equilibrium quantities. The govering equation for the displacement field, called the Frieman-Rotenberg equation (Frieman & Rotenberg 1960), is

$$\rho \frac{\partial^{2}\vec{\xi}}{\partial t^{2}} + 2\rho{\... ...rtial \vec{\xi}}{\partial t} - {\bf\mathsf{F}}(\vec{\xi}) = 0,$$ (7)

where  ${\bf\mathsf{F}}(\vec{\xi})$ represents the force operator, defined by  

$${\bf\mathsf{F}}(\vec{\xi}) \equiv -\nabla \Pi + {\vec{B}}\cdot\na... ...{\vec{v}}\cdot\nabla {\vec{v}} - \rho{\vec{v}}{\vec{v}}\cdot\nabla \vec{\xi} ],$$ (8)

with the Eulerian perturbation of the total pressure,

$$\Pi \equiv -\gamma p\nabla \cdot \vec{\xi} - \vec{\xi}\cdot\nabla p + {\vec{B}}\cdot{\vec{Q}},$$ (9)

and the Eulerian perturbation of the magnetic field,

$${\vec{Q}} \equiv \nabla \times (\vec{\xi} \times {\vec{B}}).$$ (10)

By assuming an exponential time-dependence  ${\rm e}^{-{\rm i}\omega_{{\rm }}t}$ for the displacement field, we can distinguish two cases. In the case of no equilibrium flow, the force operator  ${\bf\mathsf{F}}(\vec{\xi})$ is self-adjoint, meaning that its eigenvalues $\omega_{{\rm }}$ are purely real or imaginary. This results in stable, damped and unstable modes. When there is flow, the force operator is no longer self-adjoint, meaning that its eigenvalues  $\omega_{{\rm }}$ are complex in general. This introduces the possiblity of damped stable waves as well as overstable modes.

2.2 System of first order differential equations

The Frieman-Rotenberg Eq. (7) will be applied in the case of a one-dimensional cylindrical plasma equilibrium. For this kind of equilibrium, the MHD Eqs. (1)-(4) reduce to the radial force balance equation,

$$\left( p + \frac{1}{2}B^{2} \right)' + \frac{B_{\theta}^{2}}{r} = \frac{\rho v_{\theta}^{2}}{r} - \rho g,$$ (11)

where the prime indicates the derivative with respect to r. The symbols B, $B_{\theta}$, and $v_{\theta}$ are the total magnetic field, toroidal magnetic field, and toroidal velocity, respectively. Furthermore, g represents the gravitational acceleration at the distance r,

$$g = \frac{GM_{*}}{r^{2}},$$ (12)

with G the gravitational constant and M_* the mass of the central object.

For the three-dimensional perturbations, we choose Fourier mode solutions of the form

$$\vec{\xi}(r,\theta,z,t) = \left( \begin{array}{c} \xi_{r,... ...nd{array} \right) {\rm e}^{{\rm i}(m\theta + kz - \omega t)},$$ (13)

where m and k are the toroidal and axial wavenumber, respectively. This choice can be made because of the symmetry of the equilibrium. Also by exploiting a projection based on the magnetic field lines, the Frieman-Rotenberg Eq. (7) can be reduced to a system of first order differential equations, which reads

$$\frac{AS}{r} \left(\begin{array}{r} \chi \\ \Pi \end{ar... ...\left( \begin{array}{r} \chi \\ \Pi \end{array}\right) = 0,$$ (14)

where

$$\chi \equiv r\xi_{r},$$ (15)

$$\equiv \rho\widetilde{\omega}^{2} - F^{2},$$ A (16)

$$\equiv (\gamma p + B^{2})\rho\widetilde{\omega}^{2} - \gamma p F^{2},$$ S (17)

$$\equiv \rho^{2}\widetilde{\omega}^{4} - k_{{\rm0}}^{2} S,$$ D (18)

$$\equiv -\frac{\rho\widetilde{\omega}^{2}}{r} \Biggl\{ \frac{\rho}{r} \le... ...widetilde{\omega}^{2} + (B_{\theta}\widetilde{\omega}+ Fv_{\theta})^{2} \right]$$ C

$$\quad \qquad \qquad + \rho gA \Biggr\} + 2 \frac{m}{r}\frac{S}{r^{2}} (FB_{\theta} + \rho\widetilde{\omega}v_{\theta}),$$ (19)

and

$$\equiv -\frac{AS}{r^{2}} \left[ A + r\left( \frac{B_{\theta}^{2} - \rho v_{\theta}^{2}}{r^{2}} \right)' + \rho' g \right]$$                                       E

$$\quad + 4\frac{S}{r^{4}} (B_{\theta}F + \rho\widetilde{\omega}v_{\theta})^{2}$$

$$\quad - \frac{1}{r^{2}} \Biggl\{ \frac{\rho}{r} \left[ (B_{\theta... ...} + (B_{\theta}\widetilde{\omega}+ Fv_{\theta}) \right] + \rho gA \Biggr\}^{2}.$$ (20)

In these expressions we use the following definitions: ${\vec{k_{{\rm0}}}} = \left[ 0,m/r,k \right]$, $k_{{\rm0}}=\vert{\vec{k_{{\rm0}}}}\vert$, the Doppler shifted frequency $\widetilde{\omega}= \omega - {\vec{k_{{\rm0}}}}\cdot{\vec{v}}$, and $F = {\vec{k_{{\rm0}}}}\cdot{\vec{B}}$. Notice that we have suppressed the subscripts m and k. The system (14) has been derived by Hameiri (1981) and Bondeson et al. (1987) for a cylindrical plasma with flow but without gravity. Recently, this system was obtained by Keppens et al. (2002).

The differential Eqs. (14) become singular when A = 0 or S = 0, which results in the MHD continua. When A=0, the eigenfrequency $\omega$ is equal to the local Doppler shifted Alfvén continuum frequency,

$$\Omega_{{\rm A}}^{\pm} \equiv {\vec{k_{{\rm0}}}}\cdot{\vec{v}} \pm \omega_{{\rm A}},$$ (21)

where $\omega_{{\rm A}} \equiv F / \sqrt{\rho}$ is the local Alfvén continuum frequency for static equilibrium. If S=0, the eigenfrequency $\omega$ is equal to the local Doppler shifted slow continuum frequency:

$$\Omega_{{\rm S}}^{\pm} \equiv {\vec{k_{{\rm0}}}}\cdot{\vec{v}} \pm \omega_{{\rm S}},$$ (22)

where $\omega_{{\rm S}} \equiv \sqrt{\gamma p / (\gamma p + B^{2})} \ F / \sqrt{\rho}$ is the local static slow continuum frequency. These MHD continua form the basic organizing structure of the full MHD spectrum, and are schematically shown in Fig. 1, where the equilibrium is assumed to be weakly inhomogeneous. The Eulerian entropy continuum (Goedbloed et al. 2004)  $\Omega_{{\rm0}} \equiv {\vec{k_{{\rm0}}}}\cdot{\vec{v}}$ is also included in the figure.

	Figure 1: Continuous parts of the MHD spectrum for a weakly inhomogeneous equilibrium.
Open with DEXTER

2.3 Local dispersion equation

One can study many of the occuring instabilities by means of the local dispersion equation,

$$q^{2}\frac{A^{2}S^{2}}{r^{2}} + C^{2} + DE - \frac{ASD}{r}\left( \frac{C}{D} \right)' = 0,$$ (23)

where q is the local radial "wavenumber''. This local dispersion equation can be derived from the differential Eqs. (14) using the WKB method as discussed in appendix A. The equation will be solved neglecting the term proportional to ASD/r. This approximation is valid if one assumes  $kc_{{\rm s}} \gg \widetilde{\omega}$, where  $c_{{\rm s}} \equiv \sqrt{\gamma p / \rho}$ is the sound speed.

The remaining local dispersion equation can be rewritten as a sixth-order polynomial in  $\widetilde{\omega}$ which governs all discrete local modes. This dispersion equation reads:

$$\widetilde{\omega}^{6} + a_{4}\widetilde{\omega}^{4} + a_{3... ...2}\widetilde{\omega}^{2} + a_{1}\widetilde{\omega}+ a_{0} = 0,$$ (24)

where the coefficients a_i are defined as

     

$$\equiv - \Biggl[ (q^{2} + k_{{\rm0}}^{2})\frac{\gamma p + B^{2}}{\rho} + \frac{F^{2}}{\rho}$$                                        a₄

$$\quad \qquad - \frac{r}{\rho} \left( \frac{B_{\theta}^{2}}{r^{2}} \right)' + \kappa^{2} + \frac{\rho'}{\rho}V_{g} \Biggr],$$ (25)

$$\equiv - 4\Omega \left[ \frac{m}{r}V_{g} + 2\frac{kB_{z}B_{\theta}}{r\rho} \right],$$ a₃ (26)

$$\equiv (q^{2} + k_{{\rm0}}^{2})\frac{2\gamma p + B^{2}}{\rho} \frac{F^{2}}{\rho}$$ a₂

$$\quad + k_{{\rm0}}^{2}\frac{\gamma p + B^{2}}{\rho} \Biggl[ - \frac{r}{\rho} \left( \frac{B_{\theta}^{2}}{r^{2}} \right)' + r{\Omega^{2}}'$$

$$\quad + \frac{\rho'}{\rho}V_{g} - \frac{\rho V_{g}^{2}}{\gamma p + B^{2}} \Biggr]$$

$$\quad + 4k^{2}\Omega^{2}\frac{\gamma p + B^{2}}{\rho} - 4k\left( \frac{m}{r}B_{z} - kB_{\theta} \right) \frac{B_{\theta}}{r\rho}V_{g}$$

$$\quad + 4\frac{\Omega^{2}F}{\rho} \left( \frac{m}{r}B_{\theta} - kB_{z} \right) - 4\frac{k^{2}}{r^{2}} \frac{B_{\theta}^{2}B^{2}}{\rho^{2}},$$ (27)

$$\equiv 4 \frac{\Omega FB_{\theta}}{\rho} \left[ k_{{\rm0}}^{2}V_{g} + 2k^{2}\frac{\gamma p}{r\rho} \right],$$ a₁ (28)

$$\equiv - \frac{\gamma p F^{2}}{\rho^{2}} \Biggl\{ (q^{2} + k_{{\rm0}}^{2})\frac{F^{2}}{\rho}$$ a₀

$$\quad + k_{{\rm0}}^{2} \left[ - \frac{r}{\rho} \left( \frac{B_{\t... ...Omega^{2}}' + \frac{\rho'}{\rho}V_{g} - \frac{\rho V_{g}^{2}}{\gamma p} \right]$$

$$\quad - 4k^{2}\frac{B_{\theta}^{2}}{r^{2}\rho} \Biggr\},$$ (29)

with rotation frequency  $\Omega \equiv v_{\theta}/r$, deviation from a Keplerian disk  $V_{\rm g} \equiv v_{\theta}^{2}/r - g$ and epicyclic frequency  $\kappa^{2} \equiv 2v_{\theta}(rv_{\theta})'/r^{2}$. The local dispersion Eq. (23) reduces to the local dispersion equation derived by Dubrulle & Knobloch (1993) if one considers axisymmetric perturbations, a constant density, a zero radial "wavenumber'' q and takes the incompressible limit.

2.4 Magneto-rotational instability

By making some extra assumptions, we can reduce the local dispersion Eq. (23) further. Consider axisymmetric perturbations and the situation where $qc_{\rm s}$, $kc_{\rm s}$, $qv_{{\rm A}}$ and $kv_{{\rm A}} \gg \widetilde{\omega}$ and a purely axial magnetic field. Here, $v_{{\rm A}} \equiv B_{z}/\sqrt{\rho}$ is the Alfvén speed. In that case, the local dispersion Eq. (24) reduces to the following 4th order polynomial:

$$b_{4}\widetilde{\omega}^{4} + b_{2}\widetilde{\omega}^{2} + b_{0} = 0,$$ (30)

where the coefficients are

$$(q^{2} + k^{2}) \frac{\gamma p + B_{z}^{2}}{\rho} + \frac{F^{2}}{\rho},$$                                      b₄ = (31)

$$- \Biggl\{ (q^{2} + k^{2})\frac{2\gamma p + B_{z}^{2}}{\rho} \frac{F^{2}}{\rho}$$ b₂ =

$$+ k^{2}\frac{\gamma p+B_{z}^{2}}{\rho} \left[ r{\Omega^{2}}' + \frac{\rho'}{\rho}V_{g} - \frac{\rho V_{g}^{2}}{\gamma p + B_{z}^{2}} \right]$$

$$+ 4k^{2}\Omega^{2} \frac{\gamma p}{\rho} \Biggr\},$$ (32)

$$\frac{\gamma p F^{2}}{\rho^{2}} \Biggl\{ (q^{2} + k^{2})\frac{F^{2}}{\rho}$$ b₀ =

$$+ k^{2}\left[ r{\Omega^{2}}' + \frac{\rho'}{\rho}V_{g} - \frac{\rho V_{g}^{2}}{\gamma p} \right] \Biggr\}.$$ (33)

This 4th order polynomial can be solved analytically, which results in the stability criterion

$$\frac{F^{2}}{\rho} \ge \frac{-k^{2}}{q^{2}+k^{2}} \left[ ... ...c{\rho'}{\rho}V_{g} - \frac{\rho V_{g}^{2}}{\gamma p} \right].$$ (34)

This is the stability criterion for the axisymmetric MRI for disks with a purely axial magnetic field. This criterion is more generally applicable than the one derived by Balbus and Hawley (1991) to which it reduces for a weakly magnetized disk ( $p \gg B_{z}^{2}$).

3 The accretion disk model in the cylindrical limit

In order to quantify instabilities using a linear analysis, one has to consider the equilibrium state of an accretion disk. The model we consider uses power-law scalings for the different flow variables, following the self-similar models of Spruit et al. (1987). In order to have a model that is in an equilibrium state, these profiles have to satisfy Eq. (11). We use the following profiles for the density $\rho$, thermal pressure p, toroidal magnetic field $B_{\theta}$, axial magnetic field B_z and the toroidal velocity  $v_{\theta}$:

$$\rho = r^{-3/2},$$ (35)

$$p = \epsilon^2 \;\; r^{-5/2},$$ (36)

$$B_{\theta} = -\alpha_1 \sqrt{{2\epsilon^2\over\beta(\alpha_1^2+\alpha_2^2)}}\;\; r^{-5/4},$$ (37)

$$B_{z} = \alpha_2 \sqrt{{2\epsilon^2\over\beta(\alpha_1^2+\alpha_2^2)}}\;\; r^{-5/4},$$ (38)

$$v_{\theta} = V_0\;\; r^{-1/2},$$ (39)

where the $\alpha$-parameters express the ratio of toroidal or axial to the total magnetic field,

$$\frac{B_{\theta}}{B} = \frac{-\alpha_{1}}{\sqrt{\alpha_{1}^{2} + \alpha_{2}^{2}}},$$ (40)

$$\frac{B_{z}}{B} = \frac{ \alpha_{2}}{\sqrt{\alpha_{1}^{2} + \alpha_{2}^{2}}},$$ (41)

and

$$V_0^2 = GM_{*} - {\epsilon^2\over 2\beta\left(\alpha_1^2 + \a... ...igl[ 5(1+\beta)(\alpha_1^2 + \alpha_2^2) - 4\alpha_1^2 \bigr].$$ (42)

Finally, the parameter $\beta$ measures the radially constant ratio between the thermal pressure and the total magnetic pressure, i.e. $\beta\equiv 2p/B^2$. Notice that the density and pressure profiles are such that the entropy is constant throughout the disk, thereby excluding convective modes.

One of the key parameters of our model is the parameter $\epsilon$, which is taken as a constant free parameter. For values of $\epsilon \ll 1$, the physical interpretation is the same as in the model of Shakura & Sunyaev (1973), i.e. $\epsilon=c_{\rm s}/v_{\theta}\simeq H/r$, where His the scale height of the disk. These cases correspond to thin Keplerian rotating accretion disks. For values of $\epsilon\sim 1$, one is in the regime of sub-Keplerian rotating disks (see e.g. Narayan & Yi 1994) and for these cases we stick to the interpretation $\epsilon\simeq H/r$, noting that the identification $\epsilon=c_{\rm s}/v_{\theta}$ is not valid anymore. Typically $c_{\rm s}/v_{\theta}$ will be larger, up to a maximum factor of order $\sim$10.

Models of accretion disks which only depend on the radius r are sometimes referred to as accretion disks in the cylindrical limit (see for example Armitage 1998 and Hawley 2001). The model we use is an example of such a model, where the results from these models approximate the interior equatorial region of the accretion disk. This approximation is correct as long as the axial wavenumber k is much larger than the inverse of the scale height H. Therefore the identity $k \gg 2\pi /(\epsilon r)$ has to be satisfied in order to connect the results from our spectral analysis with MHD stability properties of the interior of an accretion disk.

4 Numerical codes

For the investigation of the stability properties two numerical codes, LEDAFLOW and a Local Dispersion Equation Solver (LODES), have been used. In the next two subsections we briefly explain their algorithmic details.

4.1 The spectral code LEDAFLOW

The code LEDAFLOW developed by Nijboer et al. (1997) can compute the complete MHD spectrum for a given one-dimensional stationary equilibrium, axial "wavenumber'' k, and azimuthal wavenumber m. Instead of numerically solving the system (14) the full set of MHD Eqs. (1)-(4) are linearized and then discretized using finite elements in the inhomogeneous direction, and solved using the Galerkin method (Schwarz 1988). For the elements, a combination of cubic Hermite and quadratic elements for the perturbed quantities is used to prevent the creation of spurious eigenvalues. The result is a non-Hermitian eigenvalue problem, which is solved using a QR method or an inverse vector iteration method (Press et al. 1988). The latter one can be used to calculate the eigenfunctions for a single eigenvalue. The used boundary conditions treat the edges of the disk as perfect conducting walls. Strictly speaking, these boundary conditions are inappropriate for considering waves modes in an accretion disk. However, interior local modes that do not depend strongly on the precise variations in a boundary region are barely affected by these conditions.

4.2 The local dispersion equation solver

The local dispersion equation solver LODES directly computes the roots of the local dispersion Eq. (24) at a given radial position. This is done by making use of Laguerre's method (Press et al. 1988).

The code calculates the local value of coefficients a_i (25)-(29) depending on the given equilibrium, the radial "wavenumber'' q, the azimuthal wavenumber m and the axial wavenumber k. In principle the code can compute all six roots, but there is a switch to calculate only the mode with the largest growth rate. Because the code uses Laguerre's method (Press et al. 1988), it needs an initial guess for this root. As an initial guess, we use the analytical solution of the most unstable mode of the 4th order polynomial (30).

   
4.3 Coupling LEDAFLOW with LODES

As described in the previous subsection, LODES needs a value of the radial "wavenumber'' q and a radial position $r_{\rm i}$ to calculate the roots. In this subsection, we show how we extract these input parameters from the results of LEDAFLOW.

	Figure 2: Complete MHD spectrum of a weakly magnetized thin disk for wavenumbers m=0 and k=200.
Open with DEXTER

	Figure 3: The eigenfunction $\chi$ of the mode with the largest growth rate. The growth rate is $0.728\Omega$.
Open with DEXTER

Figure 2 shows the complete MHD spectrum of a calculation performed by LEDAFLOW with the following set of parameters: m=0, k=200, $\beta = 1000$, $\epsilon = 0.1$ and a purely axial magnetic field ( $\alpha_{1} = 0$ and $\alpha_{2} = 1$). In this calculation, we use 101 grid points in the radial direction on the domain r=[1,2]. It is clear from the resulting MHD spectrum that the disk is not stable. We identify these unstable modes as MRI modes. Notice that, for a purely axial magnetic field, the eigenvalues  $\omega^{2}$ are real (i.e. no overstability in this case).

The eigenfunction $\chi$ of the mode with the largest growth rate is shown in Fig. 3. From this eigenfunction the location $r_{\rm i}$ and the radial wavenumber q of the instability can be extracted. The location $r_{\rm i}$ is the radius where $\chi$ has its extremum. The radial wavenumber q follows from  $q = 2\pi / \lambda_{r}$, where $\lambda_{r}$ is chosen as the radial wavelength from the inner boundary to the point when the value of the eigenfunction is 0.3% of its extremum value. The justification of this heuristic method to calculate the radial `wavenumber' q can be found in the appendix, where the WKB solutions are matched to the analytical solutions in the turning point regions.

The mode shown in Fig. 3 has also been investigated on the enlarged radial domains r=[0.966,2], r=[0.938,2] and r=[0.915,2]. The inner boundary of these domains is such that the mode has its extremum at the same radial location for the different radial domains. The resulting growth rates are all almost equal to each other, which confirms our earlier statement that the mode does not strongly depend on the boundary conditions used.

5 The magneto-rotational instability

We present results on magneto-rotational instabilities in an accretion disk in the cylindrical limit using the methods described in the previous section. The calculations using LEDAFLOW are performed on a radial domain r=[1,2]. We only consider axisymmetric perturbations (m = 0) and take a ratio of specific heats  $\gamma = 5/3$. We discuss seven calculations of which the different sets of parameters are shown in Tables 1 and 2.

5.1 Influence of the toroidal magnetic field

First, we investigate the influence of the toroidal magnetic field on both the growth rate and the oscillation frequency of the MRI. This is done by varying the toroidal magnetic field strength, while keeping the axial magnetic field strength constant. This means that the plasma beta varies between the different calculations. The parameters of the four calculations are shown in Table 1. In these calculations we consider a thin disk by taking $\epsilon = 0.1$.

The results are shown in Figs. 4-7. In each figure, the left panel shows the growth rate and the right panel shows the oscillation frequency of the most unstable axisymmetric mode, both as a function of the axial wavenumber. We have scaled the growth rate and the oscillation frequency with respect to the rotation frequency. Instead of the axial wavenumber k, we use the corresponding local Alfvén frequency  $\omega_{{\rm A}}$scaled with respect to the local rotation frequency. This scaling is similar to the one used by Narayan et al. (2002).

In our figures, the diamonds represent the calculations performed with LEDAFLOW. The solid and dashed line correspond with the results from respectively the analytical solution of the approximate 4th order polynomial (30) and LODES. Recall that the last two calculations make use of LEDAFLOW results, in the manner described in Sect. 4.3.

In the calculations A, B, and C (shown in Figs. 4-6, respectively) we change from a purely axial (A) to a predominantly toroidal (C) magnetic field configuration. One sees that the growth rate of the analytical solution of 4th order polynomial (30) and LODES matches with the LEDAFLOW calculations. The differences to LEDAFLOW are less than 4% and 1%, respectively. These calculations (A, B and C) also show that the oscillation frequency increases away from zero as the toroidal magnetic field strength increases. A similar result was obtained by Dubrulle & Knobloch (1993). These authors considered an incompressible plasma with a similar toroidal velocity profile as used in this paper, but a different toroidal magnetic field profile. Again, there is a perfect match between the LEDAFLOW and LODES results. However, the 4th order polynomial (30) cannot be used to compute the oscillation frequency since it only yields a value for the growth rate. The differences between the LEDAFLOW and LODES results are again less than 1%.

	Figure 4: The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode as a function of the scaled Alfvén frequency for disk model A, where $\beta = 1000$. In this purely axial magnetic field configuration, the axisymmetric modes are purely exponentially growing in time.
Open with DEXTER

	Figure 5: The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode as a function of the scaled Alfvén frequency for disk model B, where $\beta = 500$.
Open with DEXTER

Figure 7 shows results from a calculation for which the accretion disk is close to equipartition in a strong toroidal  $\vert B_{\theta}\vert=30\vert B_{z}\vert$ magnetic field configuration (disk model D). Again, there is a perfect match between the results from LODES and LEDAFLOW (the discrepancy is less than 1%), but the 4th order polynomial (30) no longer predicts a correct approximation of the growth rate. This is due to the dynamical importance of the toroidal magnetic field component. Comparing Figs. 7 with 6, it is clear that the oscillation frequency remains roughly the same and the growth rate decreases slightly. The analytical work performed by Dubrulle & Knobloch (1993) also indicates a decrease in the growth rate, as a result of the increasing tension of the toroidal magnetic field. Hence, the cases A-D clearly demonstrate that the MRI remains active when the toroidal magnetic field increases. This remains true at least up to equipartition (Fig. 7).

5.2 Keplerian to sub-Keplerian disk

In this subsection, we investigate the properties of the axisymmetric MRI in sub-Keplerian disks. Three calculations have been performed, all with $\beta = 1000$, $\alpha_{1}=10$ and $\alpha_{2} = 1$ (see Table 2). A sub-Keplerian disk can be obtained by increasing the value of the $\epsilon$ parameter beyond $\epsilon =0.3$. All calculations have been performed such that the domain of the Alfvén frequency over the rotation frequency remains the same: the corresponding domain of k-values can be found in Table 2.

Table 1: Parameters of thin, Keplerian accretion disks.

Table 2: Parameters of Keplerian to sub-Keplerian accretion disks.

The results are shown in Figs. 8-10. Again, the growth rate (left panel) and oscillation frequency (right panel) are scaled with respect to the rotation frequency. These scaled quantities have been plotted as a function of the ratio of the Alfvén frequency to the rotation frequency.

Again, there is a perfect match (less than 1% discrepancy) for the growth rate obtained with the three methods: LEDAFLOW, LODES and the 4th order polynomial (30). Notice that the scaled growth rates in Figs. 8-10 coincide although the deviation from a Keplerian disk is significant in cases F and G. This can be explained by investigating the solution of the 4th order polynomial in the case of a weakly magnetized disk. The solution, scaled with respect to the rotation frequency, reads:

 

$$\frac{\widetilde{\omega}^{2}}{\Omega_{{\rm }}^{2}} \frac{\omega_{{\rm A}}^{2}}{\Omega_{{\rm }}^{2}} + \frac{k^{2}}{2... ...left( \frac{\rho'}{\rho}V_{\rm g} - \frac{\rho V_{\rm g}^{2}}{\gamma p} \right)$$ =

$$\pm \sqrt{ \left[ \frac{\kappa^{2}}{\Omega_{{\rm }}^{2}} + \frac{... ...^{2}+k^{2}}{k^{2}} \frac{\omega_{{\rm A}}^{2}}{\Omega_{{\rm }}^{2}} } \Biggr\}.$$ (43)

For the chosen power-law of the toroidal velocity, the epicyclic frequency  $\kappa^{2}=\Omega_{{\rm }}^{2}$. In the case of a weakly magnetized disk and a constant entropy, the term  $\rho' V_{\rm g}/\rho - \rho V_{\rm g}^{2}/(\gamma p)$ can be neglected. Furthermore, notice that for a pure hydrodynamical case this term is the Brunt-Väisala frequency. For the calculations in this paper, this frequency is equal to zero because we consider disks with constant entropy.

This results in a equation that only depends on the ratio $\omega_{{\rm A}}/\Omega_{{\rm }}$. Remember that the domain of this ratio has been kept the same in calculations E, F and G.

The oscillation frequencies in Figs. 8-10 obtained from LEDAFLOW and LODES match perfectly (less than 1%). These figures show an overall increase in the value of the scaled oscillation frequency as one increases the value of $\epsilon$. Hence, the overstable nature of the MRI is even more significant in sub-Keplerian disks with a dominant toroidal magnetic field component.

	Figure 6: The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode as a function of the scaled Alfvén frequency for disk model C, where $\beta = 9.901$.
Open with DEXTER

	Figure 7: The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode for a disk close to equipartition with a dominant toroidal magnetic field component (disk model D).
Open with DEXTER

	Figure 8: The growth rate and the oscillation frequency of the most unstable mode for disk model E, where the disk rotation is Keplerian.
Open with DEXTER

	Figure 9: The growth rate and the oscillation frequency of the most unstable mode for disk model F, where the disk rotation is weaklysub-Keplerian.
Open with DEXTER

	Figure 10: The growth rate and the oscillation frequency of the most unstable mode for disk model G, where the disk rotation is strongsub-Keplerian.
Open with DEXTER

6 Conclusions

We have considered the magnetorotational overstability in magnetized accretion disks with a toroidal magnetic field component. We have used results from the spectral code LEDAFLOW and compared them with more approximate solutions of both the fourth and sixth order polynomial local dispersion equation. In our calculations, we have considered axisymmetric perturbations (m=0) within a disk in the cylindrical limit. The most important results are summarised below:

• Magnetorotational instabilities are present in both sub-Keplerian and Keplerian rotating accretion disks.
• The magnetorotational instability is also present in Keplerian rotating disks close to equipartition, when the toroidal magnetic field strength dominates the total magnetic field strength.
• The oscillation frequency of the instability is nonzero due to the presence of a toroidal magnetic field component for both sub-Keplerian and Keplerian rotating disks.

In an accompanying paper by van der Swaluw et al. (2005), we will consider accretion disk models that also allow for convective instabilities, and analyse their relation with respect to MRI modes. All our solvers can be used to study non-axisymmetric perturbations $(m\ne 0)$. This is left to future work.

Acknowledgements

J. W. S. Blokland, R. Keppens and J. P. Goedbloed carried out this work within the framework of the European Fusion Programme, is supported by the European Communities under the contract of the Association between EURATOM/FOM. Views and opinions expressed herein do not necessarily reflect those of the European Commission. E. van der Swaluw did this research in the FOM projectruimte on "Magnetoseismology of accretion disks'', a collaborative project between R. Keppens (FOM Institute Rijnhuizen, Nieuwegein) and N. Langer (Astronomical Institute Utrecht). This work is part of the research programme of the "Stichting voor Fundamenteel Onderzoek der Materie (FOM)'', which is financially supported by the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)''.

   
Appendix A: Radial "wavenumber''

To justify how we calculate the radial "wavenumber'' q as discussed in Sect. 4.3, we discuss the WKB approximation in detail. Instead of using the system of first order differential Eqs. (14) we use the equivalent generalized Hain-Lüst equation,

$$(f \chi')' + \hat{g}\chi = 0,$$ (A.1)

where

$$\equiv \frac{AS}{rD},$$                               f (A.2)

$$\hat{g} \equiv \frac{-r}{ASD} \left[ C^{2} + DE - \frac{ASD}{r} \left( \frac{C}{D} \right)' \right].$$ (A.3)

This equivalent equation is preferred over system (14) for a more straightforward WKB analysis. In standard WKB analysis for the generalized Hain-Lüst Eq. (A.1) we write the solution in the following form:

$$\chi (r) = p(r) \exp \left[ {\rm i} \int_{r_{0}}^{r} q(r) {\rm d}r \right],$$ (A.4)

where r₀ is the radial location of the inner boundary. Inserting this solution in the generalized Hain-Lüst Eq. (A.1) gives

$$-fq^{2}p + \hat{g}p + fp'' + f'p' + {\rm i} \left( 2fqp' \!+\! f'qp \!+\! fq'p \right) = 0.$$ (A.5)

Introducing a length scale L for the radial background variation, we can neglect the second order terms  fp'' + f'p' compared to the first order ones if $qL \gg 1$. This kind of ordering yields the following expressions for the radial "wavenumber'' q and the amplitude p:

$$q \approx \pm \sqrt{\frac{\hat{g}}{f}} \quad \textrm{and} \quad p \approx \frac{1}{ (f\hat{g})^{1/4} }\cdot$$ (A.6)

Using the expression for q, we can plot this quantity for a given radial domain. This is shown in Fig. A.1 for the eigenfunction $\chi$ plotted in Fig. 3.

Figure A.1: The radial "wavenumber'' q of the eigenfunction $\chi$ shown in Fig. 3. The growth rate is $0.728\Omega$.

Here, we applied the same radial domain as used in the LEDAFLOW calculation. The plot shows that at a certain radius r₁ the radial "wavenumber'' becomes zero. Close to this radius the WKB approximation fails, because $qL \gg 1$ is no longer satisfied. However, we can find an analytical solution close to this radius. To find the analytical solution we write the generalized Hain-Lüst Eq. (A.1) in the following form

$$\frac{{\rm d}^{2}\chi}{{\rm d}R^{2}} + G \chi = 0\cdot$$ (A.7)

where

$$R' \equiv \frac{1}{f}, \quad \textrm{and} \quad G \equiv f\hat{g}.$$ (A.8)

Both the radial "wavenumber'' q and the function G become zero at the radius r₁ as expected. We apply a Taylor expansion of G about $R_{1} ( \equiv$R(r₁)) and insert this in the generalized Hain-Lüst Eq. (A.7). This yields

$$\frac{{\rm d}^{2}\chi}{{\rm d}z^{2}} - z \chi = 0,$$ (A.9)

where

$$z \equiv \left( \left\vert \frac{{\rm d}G}{{\rm d}R} \right\vert _{R=R_{1}} \right)^{1/3} \left( R - R_{1} \right).$$ (A.10)

The analytical solution of this differential equation is the Airy function,

$$Ai(z) = \frac{1}{\pi} \int_{0}^{\infty} \cos \left( \frac{1}{3}s^{3} + sz \right) {\rm d}s,$$ (A.11)

which for large |z| has the asymptotic form

$$\sim \frac{1}{2\sqrt{\pi} z^{1/4} } \exp \left[ -\frac{2}{3} z^{3/2} \right] \textrm{for } z > 0,$$                                  Ai(z)

$$\sim \frac{1}{ \sqrt{\pi}(-z)^{1/4}} \sin \left( \frac{2}{3}(-z)^{3/2} + \frac{1}{4} \pi \right) \textrm{for } z < 0.$$ Ai(z) (A.12)

We now calculate the integral

$$\int_{r}^{r_{1}} q(r) {\rm d}r = \frac{2}{3} (-z)^{3/2},$$ (A.13)

which, for R<0, results in the solution

$$\chi (r) = \frac{A_{0}}{ (f\hat{g})^{1/4} } \sin \left( \int_{r}^{r_{1}} q(r){\rm d}r + \frac{1}{4}\pi \right),$$ (A.14)

where A₀ is constant. We apply to this function the same boundary condition as in LEDAFLOW at the inner boundary. This means that $\chi(r)$ is zero at the inner boundary, which yields

$$\int_{r_{0}}^{r_{1}} q(r){\rm d}r = -\frac{1}{4}\pi + n\pi,$$ (A.15)

where n is a natural number. If we numerically calculate this integral for the radial "wavenumber'' q plotted in Fig. A.1, we find that it is approximately 0.76$\pi$. For numerical integration we used Simpson's rule (Press et al. 1988). Thus determining the radial "wavenumber'' as discussed in Sect. 4.3 gives approximately the right value for q. Hence, even for this most unstable mode, the local WKB analysis is in excellent agreement with the numerical solution.

References

• Armitage, P. J. 1998, ApJ, 501, L189 [NASA ADS] [CrossRef] (In the text)
• Balbus, S. A., & Hawley, F. H. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] (In the text)
• Bondeson, A., Iacono, R., & Bhattacharjee, A. 1987, Phys. Fluids, 30, 2167 [NASA ADS] [CrossRef] (In the text)
• Chandrasekhar, S. 1960, Proc. Nat. Acad. Sci. USA, 46, 53 (Oxford University Press) (In the text)
• Christodoulou, D. M., Contopoulos, J., & Kazanas, D. 2003, ApJ, 586, 372 [NASA ADS] [CrossRef] (In the text)
• Dubrulle, B., & Knobloch, E. 1993, A&A, 274, 667 [NASA ADS] (In the text)
• Frieman, E., & Rotenberg, M. 1960, Rev. Mod. Phys., 32, 898 [NASA ADS] [CrossRef] (In the text)
• Goedbloed, J. P., Beliën, A. J. C., van der Holst, B., & Keppens, R. 2004, Phys. Plasmas, 11, 4332 [NASA ADS] [CrossRef] (In the text)
• Hameiri, E. 1981, J. Math. Phys., 22, 2080 [NASA ADS] [CrossRef] (In the text)
• Hawley, F. H. 2001, ApJ, 554, 534 [NASA ADS] [CrossRef] (In the text)
• Hawley, F. H., & Balbus, S. A. 1991, ApJ, 376, 223 [NASA ADS] [CrossRef] (In the text)
• Hawley, J. F., Balbus, S. A., & Stone, J. M. 2001, ApJ, 554, L49 [NASA ADS] [CrossRef] (In the text)
• Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042 [NASA ADS] [CrossRef] (In the text)
• Keppens, R., Casse, F., & Goedbloed, J. P. 2002, ApJ, 569, L121 [NASA ADS] [CrossRef] (In the text)
• Narayan, R., Quataert, E., Igumenshchev, I. V., & Abramowiczm, M. A. 2002, ApJ, 577, 295 [NASA ADS] [CrossRef] (In the text)
• Narayan, R., Yi, I. 1994, ApJ, 428, L13 [NASA ADS] [CrossRef] (In the text)
• Nijboer, R. J., Holst, B. van der, Poedts, S., & Goedbloed, J. P. 1997, Comp. Phys. Com., 106, 39 [NASA ADS] (In the text)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1988, Numerical Recipes (Cambridge University Press) (In the text)
• Schwarz, H. R. 1988, Finite Element Methods (Academic Press) (In the text)
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
• Spruit, H. C., Matsuda, T., Inoue, M., & Sawada, K. 1987, MNRAS, 229, 517 [NASA ADS] (In the text)
• van der Swaluw, E., Blokland, J. W. S., & Keppens, R. 2005, A&A, 444, 347 [EDP Sciences] [CrossRef] (In the text)
• Velikhov, E. P. 1959, Sov. Phys. JETP, 36(9), 995 (In the text)