A&A 467, 21-35 (2007)
DOI: 10.1051/0004-6361:20066979

Unstable magnetohydrodynamical continuous spectrum of accretion disks

A new route to magnetohydrodynamical turbulence in accretion disks

J. W. S. Blokland1 - R. Keppens1,2,3 - J. P. Goedbloed1,3


1 - FOM-Institute for Plasma Physics "Rijnhuizen'', Association Euratom-FOM, Trilateral Euregio Cluster, PO Box 1207, 3430 BE Nieuwegein, The Netherlands
2 - Centre for Plasma Astrophysics, K.U. Leuven, Belgium
3 - Astronomical Institute, Utrecht University, The Netherlands

Received 20 December 2006 / Accepted 25 February 2007

Abstract
Context. We present a detailed study of localised magnetohydrodynamical (MHD) instabilities occurring in two-dimensional magnetized accretion disks.
Aims. We model axisymmetric MHD disk tori, and solve the equations governing a two-dimensional magnetized accretion disk equilibrium and linear wave modes about this equilibrium. We show the existence of novel MHD instabilities in these two-dimensional equilibria which do not occur in an accretion disk in the cylindrical limit.
Methods. The disk equilibria are numerically computed by the FINESSE code. The stability of accretion disks is investigated analytically as well as numerically. We use the PHOENIX code to compute all the waves and instabilities accessible to the computed disk equilibrium.
Results. We concentrate on strongly magnetized disks and sub-Keplerian rotation in a large part of the disk. These disk equilibria show that the thermal pressure of the disk can only decrease outwards if there is a strong gravitational potential. Our theoretical stability analysis shows that convective continuum instabilities can only appear if the density contours coincide with the poloidal magnetic flux contours. Our numerical results confirm and complement this theoretical analysis. Furthermore, these results show that the influence of gravity can either be stabilizing or destabilizing on this new kind of MHD instability. In the likely case of a non-constant density, the height of the disk should exceed a threshold before this type of instability can play a role.
Conclusions. This localised MHD instability provides an ideal, linear route to MHD turbulence in strongly magnetized accretion disk tori.

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

1 Introduction

Accretion disks are very common objects in astrophysics. These objects can be found in stellar systems as well as in active galactic nuclei (AGN). In the first case, the disk accretes matter onto a protostar (young stellar object (YSO)), white dwarf, neutron star, or a black hole. The typical size of these disks is a few hundred astronomical units (AU). A massive black hole is the central object in the case of an accretion disk in AGN. Here, the typical size of the disk is a hundred parsecs (pc).

Much research on accretion disk dynamics focuses on occuring accretion processes in an MHD framework. This is done linearly as well as non-linearly. The linear studies aim to understand the drivers of the accretion mechanism. This can be done by looking at waves and instabilities about a disk equilibrium in the cylindrical limit (see e.g. Keppens et al. 2002; Blokland et al. 2005; van der Swaluw et al. 2005, and references therein). In this limit the disk equilibrium is essentially one-dimensional and one can use a self-similar model like the one of Spruit et al. (1987). When instabilities are found, one must compute the resulting evolution non-linearly to see if these instabilities give rise to turbulence. This turbulence may be the source of an enhanced effective viscosity mechanism which causes an increased transport of angular momentum outwards, as needed for accretion (Shakura & Sunyaev 1973). Balbus & Hawley (1991) discussed that the MHD turbulence resulting from the magneto-rotational instability (Velikhov 1959 and Chandrasekhar 1960) could provide this angular momentum transport.

In the last years, global two or even three-dimensional magnetohydrodynamical (MHD) simulations of accretion disks have been performed (see for example Matsumoto & Shibata 1997 and Armitage 1998). In these simulations, the non-linear dynamics is usually interpreted as a direct consequence of the magneto-rotational instability (see for example Hawley et al. 2001) leading to angular momentum transport outwards. Another possibility for the initial and evolving dynamics is that both the convective and the magneto-rotational instability play an important role in the angular momentum transport (see for example Igumenshchev et al. 2003). In the latter case, these disks are known as convection-dominated accretion flows (CDAF). However,especially those global simulations starting from initial axisymmetric accretion tori (see for example Hawley 2000) are at best loosely connected to the linear stability studies. A detailed catalogue of the MHD spectrum of such two-dimensional disks is completely lacking. Moreover, the usual identification of the cause of the turbulent dynamics with the linear magneto-rotational or convective instability is mostly based on the extrapolation of the linear stability analysis from a one-dimensional to a two or three-dimensional accretion disk. Such extrapolating ignores the fact that if one performs a detailed MHD stability study of a two-dimensional accretion disk, one may find many new types of instability that all could lead to effective angular momentum transport. A recent example of such new, poloidal flow-driven type of instability is presented by Goedbloed et al. (2004a), where the Trans-slow Alfvén Continuum (TSAC) mode is shown to occur inside a disk with toroidal and super-slow poloidal flow in the presence of a strong gravitational potential.

This paper has two aims. The first one is to present the equations and the numerical solutions of two-dimensional MHD disk equilibria with toroidal flows, but without poloidal flows. This case is fundamentally different from the previously discussed one of combined poloidal and toroidal flows (Goedbloed et al. 2004a) because the equilibria are described by quite different flux functions. By considering an axisymmetric equilibrium, we are able to model a thick accretion disk without any further approximations. All the equilibria presented below are computed using the code FINESSE (Beliën et al. 2002). The second aim is to present a detailed stability analysis of these two-dimensional equilibria. The analysis is done theoretically as well as numerically. For the stability computations we have used the recently published spectral code PHOENIX (Blokland et al. 2006). To our knowledge, this kind of analysis for MHD disk equilibria with purely toroidal flow has never been presented in the astrophysical literature until now.

We will present a sample of disk equilibria where the disk plasma typically varies from strongly magnetized up to close to equipartition. Also, the rotational speed of the plasma varies from sub-Keplerian up to Keplerian. The theoretical part of the spectral analysis of these equilibria shows that the convective continuum instabilities may appear inside the disk. These instabilities are part of the continuous branches that exist in the MHD spectrum of the linear eigenmodes of the disk tori and are localised on magnetic flux surfaces. We derive a stability criterion for this instability. This criterion looks similar to the Schwarzschild criterion. However, our criterion governs convective instabilities along the poloidal magnetic field lines while the Schwarzschild criterion applies in the direction perpendicular to the poloidal magnetic field lines. This theoretical analysis has been verified by our numerical stability calculations.

The paper is organised as follows. In Sect. 2, we present the equations which govern a two-dimensional accretion disk equilibrium. In Sect. 3, we recall the essential elements from spectral theory of MHD waves and instabilities, such as the Frieman and Rotenberg formalism (Frieman & Rotenberg 1960), straight field line coordinates and representation. In Sect. 4, these are used to derive the equations for the continuous MHD spectrum and the stability criterion for the convective instability. The numerical codes FINESSE and PHOENIX are discussed in Sect. 5. In Sect. 6, we present our numerical results on the disk equilibria and their stability analysis. Finally, in Sect. 7, we summarise and present our conclusions.

   
2 Accretion disk equilibrium

We consider an axisymmetric accretion disk. Because of this symmetry, cylindrical coordinates $(R,Z,\varphi)$ are used and the equilibrium quantities only depend on the poloidal coordinates R and Z. The disk equilibrium is modeled by the ideal MHD equations,

    
                      $\displaystyle \rho\displaystyle \frac{\partial {\vec{v}}}{\partial t}$ = $\displaystyle -\rho{\vec{v}}\cdot{\vec{\nabla}}{\vec{v}} - {\vec{\nabla}}p + {\vec{j}}\times{\vec{B}} - \rho{\vec{\nabla}}\Phi,$ (1)
$\displaystyle \displaystyle \frac{\partial p}{\partial t}$ = $\displaystyle -{\vec{v}}\cdot{\vec{\nabla}}p - \gamma p{\vec{\nabla}}\cdot{\vec{v}},$ (2)
$\displaystyle \displaystyle \frac{\partial {\vec{B}}}{\partial t}$ = $\displaystyle {\vec{\nabla}}\times\left({\vec{v}}\times{\vec{B}}\right),$ (3)
$\displaystyle \displaystyle \frac{\partial \rho}{\partial t}$ = $\displaystyle -{\vec{\nabla}}\cdot\left(\rho{\vec{v}}\right),$ (4)

where $\rho$, p, ${\vec{v}}$, ${\vec{B}}$, $\Phi$ and $\gamma$ are the density, pressure, velocity, magnetic field, gravitational potential, and the ratio of the specific heats, respectively. The current density  ${\vec{j}} = {\vec{\nabla}}\times{\vec{B}}$ and the equation ${\vec{\nabla}}\cdot{\vec{B}}=0$ has to be satisfied. Furthermore, the disk equilibrium is assumed to be time-independent. The relation between the density and the thermal pressure is expressed by the ideal gas law:  $p = \rho T$. From the equations  ${\vec{\nabla}}\cdot{\vec{B}}=0$ and  ${\vec{\nabla}}\cdot{\vec{j}}=0$ it follows that the magnetic field and the current density can be written as
              $\displaystyle {\vec{B}}$ = $\displaystyle \frac{1}{R}{\vec{e}}_{\varphi}\times{\vec{\nabla}}\psi + B_{\varphi}{\vec{e}}_{\varphi},$ (5)
$\displaystyle {\vec{j}}$ = $\displaystyle -\frac{1}{R}{\vec{e}}_{\varphi}\times{\vec{\nabla}}I + j_{\varphi}{\vec{e}}_{\varphi},$ (6)

respectively. Here, $2\pi\psi$ is the poloidal flux and the poloidal stream function  $I = RB_{\varphi}$. We restrict ourselves to disk equilibria with purely toroidal flow,

 \begin{displaymath}{\vec{v}} = v_{\varphi}{\vec{e}}_{\varphi}= R\Omega{\vec{e}}_{\varphi}.
\end{displaymath} (7)

In this case, Eqs. (2) and (4) are trivially satisfied. The angular velocity $\Omega$ is related to the electric field,

\begin{displaymath}{\vec{E}} = -{\vec{v}}\times{\vec{B}} = \Omega{\vec{\nabla}}\psi,
\end{displaymath} (8)

and making use of the induction Eq. (3) it can be easily shown that $\Omega = \Omega(\psi)$.

Three different projections can be applied on the momentum Eq. (1). The projections are in the toroidal direction and in the poloidal plane parallel and perpendicular to the poloidal magnetic field lines. From the toroidal projection one can show that the poloidal stream function is a flux function, i.e.  $I = I(\psi)$. The projection parallel to the poloidal magnetic field results in two equations,

 
$\displaystyle \left. \displaystyle \frac{\partial p}{\partial R} \right\vert _{\psi={\rm const}}$ = $\displaystyle \rho\left( R\Omega^{2} - \displaystyle \frac{\partial \Phi}{\partial R} \right),$  
$\displaystyle \left. \displaystyle \frac{\partial p}{\partial Z} \right\vert _{\psi={\rm const}}$ = $\displaystyle -\rho \displaystyle \frac{\partial \Phi}{\partial Z},$ (9)

which have to be satisfied simultaneously. From these two equations we conclude that the pressure  $p = p(\psi; R, Z)$. The last projection, perpendicular to the poloidal magnetic field, leads to the extended Grad-Shafranov equation,

 \begin{displaymath}R^{2}{\vec{\nabla}}\cdot\left( \frac{1}{R^{2}}{\vec{\nabla}}\...
...} \psi} - R^{2}\displaystyle \frac{\partial p}{\partial \psi}.
\end{displaymath} (10)

The two equations parallel to the poloidal magnetic field (9) can be solved analytically under the extra assumption that either the temperature T, or the density $\rho$, or the entropy  $S=p\rho^{-\gamma}$ is a flux function. The assumption that the temperature is a flux function can be justified due to the high thermal conductivity along the magnetic field lines. This is true at least up to the transport time scale, which is long compared to the Alfvén time. The resulting pressure can then be written as

 \begin{displaymath}p(\psi;R,Z) = p_{0}(\psi) \exp \left[ (R^{2} - R_{0}^{2}) \Lambda_{T}(\psi) - \frac{\Phi(R,Z)}{T(\psi)} \right],
\end{displaymath} (11)

where $\Lambda_{T} \equiv \Omega^{2} / (2T)$ is a flux function, p0 is the pressure for a static pure Grad-Shafranov equilibrium without gravity, and R0 is the geometric axis of the accretion disk. The extended Grad-Shafranov Eq. (10) reduces to

 \begin{displaymath}R^{2}{\vec{\nabla}}\cdot\left( \frac{1}{R^{2}}{\vec{\nabla}}\...
...^{2} - R_{0}^{2} \right) \Lambda_{T} - \frac{\Phi}{T} \right].
\end{displaymath} (12)

Another possibility, on MHD time scales, is to assume that the density is a flux function. In this case, the pressure reads

 \begin{displaymath}p(\psi;R,Z) = p_{0}(\psi) \left[ 1 + (R^{2}-R_{0}^{2})\Lambda_{\rho}(\psi) - \frac{\Phi(R,Z)}{T_{\rho}(\psi)} \right],
\end{displaymath} (13)

where  the quasi-temperature  $T_{\rho} \equiv p_{0} / \rho$ and $\Lambda_{\rho} \equiv \Omega^{2} / (2T_{\rho})$. Under this assumption, the extended Grad-Shafranov Eq. (10) can be written as

 \begin{displaymath}R^{2}{\vec{\nabla}}\cdot\left( \frac{1}{R^{2}}{\vec{\nabla}}\...
...^{2}-R_{0}^{2})\Lambda_{\rho} - \frac{\Phi}{T_{\rho}} \right].
\end{displaymath} (14)

The final option is to assume that the entropy S is a flux function. The advantage of this assumption is that it allows for a natural extension to include both toroidal and poloidal flows in the equilibrium, where the entropy has to be a flux function (Zehrfeld et al. 1972; Hameiri 1983). In this case, the pressure reads

 \begin{displaymath}p(\psi;R,Z) = p_{0}(\psi) \left\{ 1 + \frac{\gamma - 1}{\gamm...
...Phi(R,Z)}{T_{S}(\psi)} \right] \right\}^{\gamma/(\gamma - 1)},
\end{displaymath} (15)

and the extended Grad-Shafranov Eq. (10) reduces to

\begin{displaymath}R^{2}{\vec{\nabla}}\cdot\left( \frac{1}{R^{2}} {\vec{\nabla}}...
...ambda_{S} - \frac{\Phi}{T_{S}} \right) \right]^{-1} \right\}
\end{displaymath} (16)
   \begin{displaymath}\phantom{= -I\displaystyle \frac{{\rm d} I}{{\rm d} \psi} - R...
...} - \frac{\Phi}{T_{S}} \right] \right\}^{\gamma / (\gamma-1)},
\end{displaymath}

where the quasi-temperature  $T_{S} \equiv S\rho_{0}^{\gamma-1}$ and $\Lambda_{S} \equiv \Omega^{2} / (2T_{S})$. The extended Grad-Shafranov equation has been used previously by van der Holst et al. (2000b) for the first two cases without gravity.

The density can be easily derived for all three assumptions by inserting the corresponding equation for the pressure into the momentum equations parallel to the poloidal magnetic field lines (9). The resulting density is

\begin{displaymath}\rho(\psi; R,Z) = \rho_{0}(\psi) \times
\left\{\begin{array}{...
...t\}^{1/(\gamma-1)} & {\rm for}~~S=S(\psi).
\end{array}\right.
\end{displaymath} (17)

The flux function $\rho_{0}$ corresponds to the density of a static equilibrium without gravity. All these cases with the inclusion of gravity have been discussed in the appendix to Blokland et al. (2006).

   
3 Spectral formulation

3.1 Frieman-Rotenberg formalism

For the investigation of stability properties of stationary MHD equilibria, the formalism by Frieman & Rotenberg (1960) has been used. They derived from the linearised MHD equations, one second order differential equation for the Lagrangian displacement vector field  ${\vec{\xi}}$,

 \begin{displaymath}{\vec{F}}({\vec{\xi}}) - 2\rho{\vec{v}}\cdot{\vec{\nabla}}\di...
...al t} -\rho\frac{\partial^{2}{\vec{\xi}}}{\partial t^{2}} = 0,
\end{displaymath} (18)

where the force operator ${\vec{F}}({\vec{\xi}})$ is

 \begin{displaymath}
{\vec{F}}({\vec{\xi}}) = {\vec{F}}_{{\rm static}}({\vec{\xi...
...\rho{\vec{v}}{\vec{v}}\cdot{\vec{\nabla}}{\vec{\xi}} \right].
\end{displaymath} (19)

Here,

 \begin{displaymath}
{\vec{F}}_{{\rm static}}({\vec{\xi}}) =
-{\vec{\nabla}}\Pi...
...\vec{\nabla}}{\vec{Q}} + {\vec{Q}}\cdot{\vec{\nabla}}{\vec{B}}
\end{displaymath} (20)

is the force operator for static equilibria without gravity, which has been derived by Bernstein et al. (1958). Here, the Eulerian perturbation of the total pressure,

 \begin{displaymath}
\Pi = -\gamma p{\vec{\nabla}}\cdot{\vec{\xi}} - {\vec{\xi}}\cdot{\vec{\nabla}}p + {\vec{B}}\cdot{\vec{Q}},
\end{displaymath} (21)

and the Eulerian perturbation of the magnetic field,

 \begin{displaymath}
{\vec{Q}} = {\vec{\nabla}} \times \left( {\vec{\xi}} \times {\vec{B}} \right).
\end{displaymath} (22)

The time-dependence for the displacement field is assumed to be an exponential one with normal mode frequencies $\omega$, ${\vec{\xi}} = \hat{{\vec{\xi}}} \exp ( -{\rm i}\omega t)$. Using this assumption, the Frieman-Rotenberg equation can be written as

 \begin{displaymath}{\vec{F}}({\vec{\xi}}) + 2{\rm i}\rho\omega{\vec{v}}\cdot{\vec{\nabla}}{\vec{\xi}} + \rho\omega^{2}{\vec{\xi}} = 0.
\end{displaymath} (23)

In the derivations below, the toroidal symmetry of the accretion disk equilibrium is exploited. This is done by writing the toroidal dependence of the displacement field as  ${\vec{\xi}} \sim \exp ({\rm i}n\varphi)$, where n is the toroidal mode number.

3.2 Straight field line coordinates

The analytical and numerical analysis of MHD waves and instabilities are done in "straight field line'' coordinates. To convert from cylindrical  $(R,Z,\varphi)$ to straight field line coordinates $(x^{1}\equiv\psi, x^{2}\equiv\vartheta, x^{3}\equiv\varphi)$ one needs the metric tensor and the Jacobian associated with the non-orthogonal coordinates in which the equilibrium field lines appear straight. This is standard practice in MHD spectral studies for laboratory tokamak plasmas. The metric elements gij and the Jacobian J are

gij = $\displaystyle {\vec{\nabla}}x^{i} \cdot {\vec{\nabla}}x^{j},\hspace*{1cm} g_{ij...
...la}}\psi\times{\vec{\nabla}}\vartheta \cdot {\vec{\nabla}}\varphi \right)^{-1},$ (24)

respectively. Here, the poloidal angle $\vartheta$ is constructed such that the magnetic field lines are straight in the $(\vartheta,\varphi)$-plane. The slope of these lines is a flux function:

 \begin{displaymath}\left. \displaystyle \frac{{\rm d} \varphi}{{\rm d} \vartheta...
...B}}\cdot{\vec{\nabla}}\vartheta} = \frac{IJ}{R^{2}} = q(\psi),
\end{displaymath} (25)

where q is the safety factor. Furthermore, in straight field coordinates the poloidal and toroidal curvature of the magnetic surfaces are

\begin{displaymath}\kappa_{{\rm p}}= -{\vec{n}} \cdot \left( {\vec{t}}\cdot{\vec...
...ial_{\vartheta}\frac{g_{12}}{g_{22}} \right) JB_{\vartheta},
\end{displaymath} (26)


\begin{displaymath}\kappa_{{\rm t}}= -{\vec{n}} \cdot \left( {\vec{e}}_{\varphi}...
...al_{\psi}- \frac{g_{12}}{g_{22}}\partial_{\vartheta}\right) R,
\end{displaymath} (27)

respectively. Here, ${\vec{n}} = {\vec{\nabla}}\psi / \vert{\vec{\nabla}}\psi\vert$ and ${\vec{t}} = {\vec{B}}_{\vartheta} / B_{\vartheta}$, where  $B_{\vartheta}$ is the poloidal magnetic field. It is important to realise that the straight field coordinates can only be constructed when the solution $\psi(R,Z)$ has been computed from the extended Grad-Shafranov Eq. (10).

3.3 Field line projection and representation

On each flux surface the distinction between two wave directions can be made: one parallel and the other one perpendicular to the magnetic field. These directions correspond to the short wavelength limit of slow and Alfvén waves in cylindrical geometry. Hence, it is useful to exploit a projection based on the magnetic surfaces and field lines using the triad of unit vectors,

 \begin{displaymath}{\vec{n}} \equiv \frac{{\vec{\nabla}}\psi}{\vert{\vec{\nabla}...
...{n}}, \hspace*{10mm} {\vec{b}} \equiv \frac{{\vec{B}}}{B}\cdot
\end{displaymath} (28)

Using these unit vectors, the components of the displacement field  ${\vec{\xi}}$ can be written as

 \begin{displaymath}X \equiv RB_{\vartheta}{\vec{\xi}}\cdot{\vec{n}}, \hspace*{10...
...*{10mm} Z \equiv {\rm i}\frac{1}{B} {\vec{\xi}}\cdot{\vec{b}},
\end{displaymath} (29)

and the projected gradient operators become
   
                      $\displaystyle \mathcal{D}$ $\textstyle \equiv$ $\displaystyle \frac{1}{RB_{\vartheta}}{\vec{n}} \cdot{\vec{\nabla}} = \partial_{\psi}- \frac{g_{12}}{g_{22}}\partial_{\vartheta},$ (30)
$\displaystyle \mathcal{G}$ $\textstyle \equiv$ $\displaystyle -{\rm i}RB_{\vartheta}B {\vec{\pi}}\cdot{\vec{\nabla}} = \frac{-{\rm i}I}{J}\partial_{\vartheta}- nB_{\vartheta}^{2},$ (31)
$\displaystyle \mathcal{F}$ $\textstyle \equiv$ $\displaystyle -{\rm i}{\vec{B}} \cdot{\vec{\nabla}} = \frac{-{\rm i}}{J}\partial_{\vartheta}+ \frac{nq}{J}\cdot$ (32)

The second equality in the expressions for $\mathcal{F}$ and $\mathcal{G}$ only holds when these operators act on a component of the displacement field  ${\vec{\xi}}$.

Using the straight field line coordinates and the projections (29) the Frieman-Rotenberg Eq. (18) can be written as

 \begin{displaymath}(\textbf{\sf A} + \textbf{\sf R}){\vec{x}} + 2\rho\widetilde{...
...left(
\begin{array}{c}
X \\
Y \\
Z
\end{array}\right),
\end{displaymath} (33)

where   $\widetilde{\omega}(\psi) \equiv \omega - n\Omega(\psi)$ is the Doppler shifted eigenfrequency. Here,  $\textbf{\sf A}$ and  $\textbf{\sf B}$ are  $3 \times 3$ matrix operators which are also present in the case of a static equilibrium without gravity. The matrix operators  $\textbf{\sf C}$ and  $\textbf{\sf R}$ contain the elements due to the toroidal rotation of the plasma and of an external gravitational potential. The matrix elements of  $\textbf{\sf A}$ and  $\textbf{\sf B}$ are
 
                                     $\displaystyle {\rm A}_{11}$ $\textstyle \equiv$ $\displaystyle -\mathcal{D}\frac{\gamma p + B^{2}}{J}\mathcal{D}^{\dagger}J + F\...
...ft(\frac{1}{J} \frac{B_{\varphi}\kappa_{{\rm t}}}{B_{\vartheta}}\right)\right],$  
$\displaystyle {\rm A}_{12}$ $\textstyle \equiv$ $\displaystyle -\mathcal{D}\gamma p\mathcal{G}\frac{1}{B^{2}} - \mathcal{D}\math...
...i}\kappa_{{\rm t}}}{B_{\vartheta}}\frac{{\rm i}}{J}\partial_{\vartheta}\right),$  
$\displaystyle {\rm A}_{13}$ $\textstyle \equiv$ $\displaystyle -\mathcal{D}\gamma p\mathcal{F},$  
$\displaystyle {\rm A}_{21}$ $\textstyle \equiv$ $\displaystyle \frac{1}{B^{2}}\mathcal{G}\frac{\gamma p}{J}\mathcal{D}^{\dagger}...
...}\partial_{\vartheta}\frac{B_{\varphi}\kappa_{{\rm t}}}{B_{\vartheta}} \right),$  
$\displaystyle {\rm A}_{22}$ $\textstyle \equiv$ $\displaystyle \frac{1}{B^{2}}\mathcal{G}\gamma p\mathcal{G}\frac{1}{B^{2}} + \m...
...B^{2}}\mathcal{G}
+ \mathcal{F}\frac{R^{2}B_{\vartheta}^{2}}{B^{2}}\mathcal{F},$  
$\displaystyle {\rm A}_{23}$ $\textstyle \equiv$ $\displaystyle \frac{1}{B^{2}}\mathcal{G}\gamma p\mathcal{F},$  
$\displaystyle {\rm A}_{31}$ $\textstyle \equiv$ $\displaystyle \mathcal{F}\frac{\gamma p}{J}\mathcal{D}^{\dagger}J,$  
$\displaystyle {\rm A}_{32}$ $\textstyle \equiv$ $\displaystyle \mathcal{F}\gamma p\mathcal{G}\frac{1}{B^{2}},$  
$\displaystyle {\rm A}_{33}$ $\textstyle \equiv$ $\displaystyle \mathcal{F}\gamma p\mathcal{F},$ (34)

and

 \begin{displaymath}{\rm B}_{11} \equiv \frac{1}{R^{2}B_{\vartheta}^{2}},
{\rm ...
...c{R^{2}B_{\vartheta}^{2}}{B^{2}},
{\rm B}_{33} \equiv B^{2}.
\end{displaymath} (35)

Here, the operators

 \begin{displaymath}\mathcal{D}^{\dagger}\equiv \mathcal{D}- \left[ \partial_{\va...
...partial_{\vartheta}\left(\frac{g_{12}}{g_{22}}\right) \right],
\end{displaymath} (36)

are related to the normal gradient operator  $\mathcal{D}$, which has been introduced by Goedbloed (1997) for the spectral analysis of static tokamak equilibria. The square brackets in the expressions (34) for the matrix elements of  $\textbf{\sf A}$ and the gradient operators defined by (36) indicate that the differential operator only acts on the term inside the bracket. This notation is also used in the expressions below. The matrices  $\textbf{\sf A}$ and  $\textbf{\sf B}$ were derived by Goedbloed (19751997) for static equilibria.

The matrices  $\textbf{\sf R}$ and  $\textbf{\sf C}$ enter if there is toroidal flow and/or external gravity. The expression for these matrices are

 \begin{displaymath}\textbf{\sf R} =
\left(\begin{array}{lcc}
-\rho\left\{ \lef...
...frac{1}{J}\partial_{\vartheta}\mu \right]
\end{array}\right),
\end{displaymath} (37)


 \begin{displaymath}\textbf{\sf C} =
\left(\begin{array}{lcc}
0&-\displaystyle\...
...c{R}{J}\partial_{\vartheta}R \right] &
0
\end{array}\right),
\end{displaymath} (38)

where

 \begin{displaymath}
\mu \equiv \left[ \frac{R}{J}\partial_{\vartheta}R \right] ...
...p = \frac{1}{\rho}{\vec{B}}_{\vartheta}\cdot{\vec{\nabla}}p,
\end{displaymath} (39)


\begin{displaymath}\lambda \equiv \frac{R\kappa_{{\rm t}}}{B_{\vartheta}}\Omega^{2} - \mathcal{D}\Phi.
\end{displaymath} (40)

Notice that the term $\mu$ represents the pressure variation on a flux surface. The matrix  $\textbf{\sf C}$ represents the Coriolis effect due to the rotation while the matrix  $\textbf{\sf R}$ contains rotational as well as gravitational effects. The case without an external gravitational potential has been discussed by van der Holst et al. (2000b).

Two kinds of cylindrical limits can be obtained from the spectral Eq. (33). The first one is the infinitely slender torus, meaning that the radial position of accretion disk is taken to be at infinity. In that case the matrices $\textbf{\sf R}$ and $\textbf{\sf C}$ disappear. This is due to the fact that all equilibrium quantities will be become independent of the angle $\vartheta$, the gravitational potential at infinity is zero, and the toroidal curvature  $\kappa_{{\rm t}}$ becomes zero. The flow enters only as a Doppler shift,  $-n\Omega(\psi)$, in the Doppler shifted eigenfrequency  $\widetilde{\omega}(\psi)$.

The other limit is the cylindrical limit of a thin ( $\vert Z\vert/R \ll 1$) slice of plasma at the equatorial plane of the accretion disk. In this region all equilibrium quantities only depend on the radius R, approximately. The matrices  $\textbf{\sf R}$ and  $\textbf{\sf C}$ do not disappear as in the previous limit. Also the toroidal curvature  $\kappa_{{\rm t}}$ is non-zero but instead the poloidal curvature  $\kappa_{{\rm p}}= 0$. In this case the spectral Eq. (33) reduces to the matrix equation for a one-dimensional accretion disk presented by Keppens et al. (2002).

   
4 Continuous MHD spectrum

In the previous section, the spectral Eq. (33) governing all MHD waves and instabilities in axisymmetric accretion disks has been derived. This equation is the starting point to all following MHD spectral computations. In particular, we can derive the equations for the continuous MHD spectrum by considering modes localised on a particular flux surface  $\psi = \psi_{0}$. Van der Holst et al. (2000a) showed that the toroidal flow can drive the continuous spectrum unstable or overstable for non-gravitating plasma tori. Here, we show that the combination of toroidal flow and gravity can also drive the MHD continua overstable.

4.1 General formalism

To derive the equations for the continuous spectrum, the normal derivative  $(\partial / \partial \psi)$ of the eigenfunctions is considered to be large compared to the eigenfunctions themselves and the Doppler shifted eigenfrequency  $\widetilde{\omega}(\psi)$ is assumed finite. In this case, the first row of the spectral Eq. (33) can be solved approximately,

 \begin{displaymath}\frac{1}{J}\mathcal{D}^{\dagger}JX \approx \frac{-\gamma p}{\...
...rho\mu}{\gamma p + B^{2}} \left( \frac{I}{B^{2}}Y + Z \right).
\end{displaymath} (41)

Notice, that this solution implies that $\partial X / \partial \psi$, Y, and Z are of the same order, but more importantly that X is small compared to Y and Z. This means that the continuum modes are mainly tangential to a particular flux surface.

Inserting the expression (41) into the second and third row of the spectral Eq. (33) results in an eigenvalue problem which is independent of the normal derivative. Hence, the reduced problem becomes non-singular. Exploiting this property, we can write the projected displacement field components Y and Z as follows:

$\displaystyle Y \approx \delta (\psi - \psi_{0}) \eta (\vartheta),$      
$\displaystyle Z \approx \delta (\psi - \psi_{0}) \zeta(\vartheta),$     (42)

where $\delta (\psi - \psi_{0})$ is the Dirac delta function. Here, Y and Z have been split in an improper ($\psi$-dependence) and proper part ($\vartheta$-dependence). This kind of splitting has been introduced by Goedbloed (1975) for static equilibria without gravity. The reduced, non-singular eigenvalue problem becomes

 \begin{displaymath}%
\left( \textbf{\sf a} + \textbf{\sf r} + 2\rho\widetilde{\o...
...t( \begin{array}{l}
\eta \\
\zeta
\end{array}\right) = 0,
\end{displaymath} (43)

where

 \begin{displaymath}\textbf{\sf a} =
\left(
\begin{array}{cc}
\mathcal{F}\dis...
...amma pB^{2}}{\gamma p + B^{2}}\mathcal{F}
\end{array}\right),
\end{displaymath} (44)


 \begin{displaymath}\textbf{\sf r} =
\left(\begin{array}{lc}
-4\displaystyle \fra...
...} \right] + \rho B^{2} N_{{\rm m,pol}}^{2}
\end{array}\right),
\end{displaymath} (45)


 \begin{displaymath}\textbf{\sf c} =
\left(\begin{array}{cc}
0 &
i \left[ \dis...
...ac{R}{J}\partial_{\vartheta}R \right] & 0
\end{array}\right),
\end{displaymath} (46)


 \begin{displaymath}\textbf{\sf b} =
\left(\begin{array}{ll}
\displaystyle \fra...
...\vartheta}^{2}}{B^{2}} & 0 \\
0 & B^{2}
\end{array}\right).
\end{displaymath} (47)

In these expressions  $\kappa_{{\rm g}}$ is the geodesic curvature,

\begin{displaymath}\kappa_{{\rm g}}= {\vec{\kappa}}\cdot{\vec{\pi}}= \frac{I}{JRB_{\vartheta}B^{2}}\partial_{\vartheta}B.
\end{displaymath} (48)

Here, ${\vec{\kappa}} = {\vec{b}}\cdot{\vec{\nabla}}{\vec{b}}$ is the field line curvature and $N_{{\rm m,pol}}^{2}$ is the magnetically modified Brunt-Väisälä frequency projected on a flux surface,

\begin{displaymath}N_{{\rm m,pol}}^{2} = \frac{\mu}{B^{2}}\left[ \frac{1}{J\rho}...
...frac{\rho}{\gamma p + B^{2}} {\vec{\nabla}}p \right) \right].
\end{displaymath} (49)

This magnetically modified Brunt-Väisälä frequency is similar as the one presented by van der Holst et al. (2000b) for tokamaks with purely toroidal flow and as the one for static gravitational plasmas published by Poedts et al. (1985) and Beliën et al. (1997). Furthermore, notice that the poloidal variation of the pressure presented by the matrix  $\textbf{\sf r}$ shows up in all its matrix elements. This destroys the diagonal dominance of the matrix  $\textbf{\sf a}$, which corresponds to the separation of Alfvén and slow continuum modes. Similar as the matrix  $\textbf{\sf C}$ in spectral Eq. (33), the matrix  $\textbf{\sf c}$ represents the Coriolis effect.

The non-singular eigenvalue problem (43) is solved for poloidally periodic boundary conditionsfor the eigenfunctions $\eta$, $\zeta$. Solving this problem for a given toroidal mode number n on each flux surface separately results in a set of discrete eigenvalues. All these discrete sets together map out the continuous spectra. In the MHD formulation with the primitive variables $\rho$, ${\vec{v}}$, p, ${\vec{B}}$, an additional entropy continuum $n\Omega$ is found. This Eulerian entropy continuum does not couple to any of the other continua or stable and unstable modes (Goedbloed et al. 2004b).

4.2 Spectral properties and stability criterion

In this subsection a stability criterium for the continua will be derived. For this derivation we need to construct a Hilbert space with appropriate inner product. The parallel gradient operator  $\mathcal{F}$ is a Hermitian operator under the inner product

 \begin{displaymath}\oint J v^{*} \left( \mathcal{F}w \right) {\rm d}^{}\vartheta = \oint J \left( \mathcal{F}v \right)^{*} w {\rm d}^{}\vartheta
\end{displaymath} (50)

for poloidally periodic functions v and w. Inspired by this property a Hilbert space can be defined with the inner product

 \begin{displaymath}\langle {\vec{v}},{\vec{w}} \rangle \equiv \oint {\vec{v}}^{*}\cdot{\vec{w}} J {\rm d}^{}\vartheta .
\end{displaymath} (51)

In this Hilbert space, the matrix operators  $\textbf{\sf a}$, $\textbf{\sf b}$, $\textbf{\sf c}$, and  $\textbf{\sf r}$ are Hermitian operators. By setting ${\vec{v}}$ and ${\vec{w}}$ equal to the eigenfunctions  $(\eta,\zeta)^{{\rm T}}$, one can derive a quadratic polynomial for the eigenfrequency  $\widetilde{\omega}$ from the spectral Eq. (43) for the continua.

 \begin{displaymath}a\widetilde{\omega}^{2} - 2b\widetilde{\omega}- c = 0.
\end{displaymath} (52)

In this equation, the coefficients are
   
                                a $\textstyle \equiv$ $\displaystyle \oint \rho\left( \frac{R^{2}B_{\vartheta}^{2}}{B^{2}}\vert\eta\vert^{2} + B^{2}\vert\zeta\vert^{2} \right) J{\rm d}^{}\vartheta,$ (53)
2b $\textstyle \equiv$ $\displaystyle 2i\oint \rho\Omega\left[ \frac{R}{J}\partial_{\vartheta}R \right] \left( \eta^{*}\zeta - \eta\zeta^{*} \right) J{\rm d}^{}\vartheta,$ (54)
c $\textstyle \equiv$ $\displaystyle \oint \left\{ \frac{R^{2}B_{\vartheta}^{2}}{B^{2}} \left\vert \ma...
...vert \frac{I}{B^{2}}\eta + \zeta \right\vert^{2} \right\} J{\rm d}^{}\vartheta,$ (55)

where $N_{{\rm BV,pol}}^{2}$ is the Brunt-Väisälä frequency projected on a flux surface,
                    $\displaystyle N_{{\rm BV,pol}}^{2}$ = $\displaystyle \frac{\mu}{B^{2}}\left[ \frac{1}{J\rho}\partial_{\vartheta}\rho - \frac{\rho\mu}{\gamma p} \right]$  
  = $\displaystyle \left[ \frac{{{\vec{B}}}_{\vartheta}\cdot{\vec{\nabla}}p}{\rho B}...
...eft( {\vec{\nabla}}\rho - \frac{\rho}{\gamma p} {\vec{\nabla}}p \right) \right]$  
  = $\displaystyle -\left[ \frac{{{\vec{B}}}_{\vartheta}\cdot{\vec{\nabla}}p}{\rho B...
...]
\left[ \frac{{{\vec{B}}}_{\vartheta}\cdot{\vec{\nabla}}S}{\gamma BS} \right].$ (56)

Note that the coefficient a is always non zero but, more importantly, that the coefficients a, b, c are real due to the Hermitian property of the inner product. Solutions of this polynomial are  $\widetilde{\omega}= (b \pm \sqrt{b^{2} + ac}) / a$. This means that if $b^{2} + ac \ge 0$, the solutions are waves with real frequencies. But if b2 + ac < 0, the solutions contain a damped stable wave and an overstable mode. For the complex solution  $\widetilde{\omega}$ one derives its absolute value  $\vert\widetilde{\omega}\vert=\sqrt{-c / a}$ by taking the linear combination of the polynomial with its conjugate. Furthermore, ${\rm Re}(\omega)=n\Omega + b/a$. It can easily be shown that the continuum is stable if $c \ge 0$. The coefficient c is always equal or greater than zero if

 \begin{displaymath}N_{{\rm BV,pol}}^{2} \ge 0
\end{displaymath} (57)

is satisfied everywhere in the plasma. From this stability criterion, we conclude that equilibria with $S = S(\psi)$ and equilibria with $T = T(\psi)$, $\gamma \ge 1$ are always stable. Equilibria where the density is a flux function violate this criterion, which may result in the appearance of damped stable and overstable modes. Notice that the stability criterion (57) is analogous to the Schwarzschild criterion for convective instability. The only important difference is that in the Schwarzschild criterion one deals with normal derivatives while this criterion deals with tangential ones. This has also been noticed by Hellsten & Spies (1979), Hameiri (1983), and Poedts et al. (1985).

  
5 Numerical codes

For the stability analysis of axisymmetric accretion tori two numerical codes, the equilibrium code FINESSE (Beliën et al. 2002) and the spectral code PHOENIX (Blokland et al. 2006), have been used. In the next two subsections their algorithmic details will be briefly discussed.

5.1 The equilibrium code FINESSE


  \begin{figure}
\par\includegraphics[width=10cm,clip]{6979-01.ps} \end{figure} Figure 1: An example of the geometry of an accretion disk with a circular cross-section. Here, R0 and a are the major radius of the geometric axis and the radius of the last closed flux surface, respectively.
Open with DEXTER

The FINESSE code developed by Beliën et al. (2002) can compute a stationary axisymmetric, gravitating MHD equilibrium described by the extended Grad-Shafranov Eqs. (12), (14), or (2) for a given poloidal cross-section and inverse aspect ratio  $\epsilon \equiv a / R_{0}$. Here, a and R0 are the minor radius of the last closed flux surface and the major radius of the geometric axis, respectively. Figure 1 shows an example of the geometry of an accretion disk with a circular cross-section as seen in a poloidal plane with the central gravitational object at the origin. The extended Grad-Shafranov equation has been discretized using a finite element method in combination with the standard Galerkin method. The elements used are isoparametric bicubic Hermite ones. These bicubic elements ensure that the computed solution has the desired high accuracy (fourth order) needed for the stability analysis. FINESSE solves the elliptic PDE problem using the Picard iteration scheme. The imposed boundary conditions assume that the fixed boundary represents the last closed flux surface. The same FINESSE code can actually handle the more general case of stationary MHD equilibrium with non-vanishing poloidal flows as well, but here we restrict our discussion to equilibria with purely toroidal flows.

5.2 The spectral code PHOENIX

The spectral code PHOENIX developed by van der Holst (Blokland et al. 2006) is used for the stability analysis of a two dimensional accretion disk. PHOENIX can compute the complete MHD spectrum for a given two dimensional equilibrium, toroidal mode number n, and a range of poloidal mode numbers m. A range of poloidal mode numbers is needed due to the mode coupling caused by allowing for a non-circular cross-section, gravitational stratification, and the toroidal geometry of the accretion disk. PHOENIX does not directly solve the spectral Eq. (33), but instead employs the linearised version of the full set of MHD Eqs. (1)-(4). These linearised equations are then discretized using a finite element method in the normal $\psi$-direction, a spectral method in the poloidal direction and a standard Galerkin method is used to obtain a linear generalised eigenvalue problem for the eigenfrequencies $\omega$. The elements used in the $\psi$-direction for the perturbed quantities are a combination of quadratic and cubic Hermite elements to prevent the creation of spurious eigenvalues (Rappaz 1977). This results in a non-Hermitian eigenvalue problem, which is solved using the Jacobi-Davidson method (Sleijpen & van der Vorst 1996). The boundary conditions used for the perturbations are those of a perfect conducting wall. Strictly speaking, these conditions can not be applied to waves and instabilities that appear inside an accretion disk. However, if the wave phenomena are sufficiently localised or not significantly affected by the position of the wall, these boundary conditions have hardly any influence. This is particularly true for the continuous branches of the MHD spectrum on which we concentrate in what follows. These modes are purely localised on a flux surface, and are not affected by the wall. Like FINESSE, PHOENIX can handle poloidal flows as well, but we have exploit it here for equilibria with toroidal flow only.

The continuous spectrum is computed by a less expensive method introduced by Poedts & Schwarz (1993). In this method, on each individual flux surface the cubic Hermite and the quadratic elements are replaced by $\log (\epsilon)$and  $1/\epsilon$, respectively. Due to this replacement, the singular behaviour of the continuous spectrum has been approximated by the perturbed quantities. This singular behaviour has been described by Goedbloed (1975) and Pao (1975) for static, axisymmetric and toroidal plasma. Recently, the analogous treatment has been used by Goedbloed et al. (2004a) for plasmas with toroidal and poloidal flows and gravity. The resulting eigenvalue problem is then a small algebraic problem of a size set by the range in poloidal mode numbers m on each flux surface, and is then solved using a direct QR method and a scan over all flux surfaces. This yields detailed information on all MHD continua.

  
6 Accretion disks with density as a flux function

We present our numerical results on disk equilibria and stability using the codes described in the previous section. For the equilibrium computations the ratio of the specific heats $\gamma = 5/3$ and the gravitational potential is an external Newtonian one, i.e.  $\Phi = -GM_{*} / \sqrt{R^{2} + Z^{2}}$, where G and M* are the gravitational constant and the mass of the central object, respectively. Furthermore, we assume that the density is a flux function for all equilibrium computations. From our theoretical analysis above, this is the relevant case which may give rise to damped stable and overstable continuum mode pairs due to the presence of toroidal flow and gravity.

6.1 Thin accretion disk with constant density


  \begin{figure}
\par\includegraphics[width=6cm,clip]{6979-02.ps} \end{figure} Figure 2: The two-dimensional pressure (gray-scale) and plasma beta  $\beta = 2p/B^{2}$ (contours) profile for an accretion disk with constant density. The cross-section is circular and the inverse aspect ratio  $\epsilon =0.1$. The central object is to the left of the figure.
Open with DEXTER

As the first case, we take an accretion disk with constant density as also presented by Blokland et al. (2006). The equilibrium is described by the following flux functions:

 
$\displaystyle I^{2}(\psi) = A (1 - 0.0385\psi + 0.02\psi^{2} + 0.00045\psi^{3}), \quad \rho(\psi) = 1,$      
$\displaystyle p_{0}(\psi) = AB(1 - 0.9\psi), \quad \Omega(\psi) = C (1 - 0.9\psi^{2}),$     (58)

where the values of the coefficients A, B, and C are given by 112, 0.01, and 0.1. Furthermore, the considered poloidal cross-section is circular and the inverse aspect ratio  $\epsilon =0.1$. In these calculations, GM* = 1 after scaling with respect to a reference density, the vacuum magnetic field, and the minor radius of the accretion disk. Figure 2 shows the computed two-dimensional pressure and plasma beta  $\beta (= 2p/B^{2})$ profile. The pressure decreases monotonically outwards due to the fact that the term proportional to the gravitational potential $\Phi$ dominates in the pressure Eq. (13). In contrast, the plasma beta increases monotonically outwards from 0.63 at the inner part to 0.76 at the outer part of the disk. The ratio  $v_{\varphi} / v_{{\rm Kepler}} = [ 0.085 , 1.049]$, where  $v_{{\rm Kepler}}$ is the Keplerian velocity.


  \begin{figure}
\par\includegraphics[width=6cm,clip]{6979-03a.ps}\includegraphics...
...=6cm,clip]{6979-03b.ps}\includegraphics[width=6cm,clip]{6979-03c.ps}\end{figure} Figure 3: a) Real, and b) imaginary parts of the sub-spectrum of the MHD continua as functions of the radial flux coordinate  $s \equiv \sqrt {\psi }$ and c) represented in the complex plane for toroidal mode number n=-1 and poloidal mode numbers m=[-3,5]. The corresponding disk equilibrium shown in Fig. 2 has a constant density distribution.
Open with DEXTER

The computation of the continuous spectrum of this equilibrium is done for toroidal mode number n=-1 and poloidal mode numbers m=[-3,5]. Figure 3 shows the sub-spectrum of the MHD continua. The plots clearly indicate the existence of stable and overstable modes in the MHD continua. The plots (a) and (b) show the real and imaginary part of the eigenfrequencies $\omega$ that belong to the continuous spectrum as a function of the radial flux coordinate  $s \equiv \sqrt {\psi }$. The frequencies of the continua in the complex plane are shown in the plot (c). It should be noted that it is difficult to label the different branches of the MHD continua by a "dominant'' poloidal mode number m. This could be done by making use of analytical expressions for an infinitely slender torus ( $\epsilon \to 0$), but then the gravity of the central object would be neglected. In our case, these expressions can hardly be used to identify the different branches of the MHD continua, because gravity plays a dominant role in the equilibrium as well as in the spectral analysis. Notice that in plot (a) two branches merge at $s \approx 0.59$, while at the same location the imaginary part of the eigenfrequency starts to become non-zero. Intricate mode couplings thus give rise to the emergence of pairs of overstable and damped stable modes.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{6979-04.ps} \end{figure} Figure 4: The growth rate of the most unstable continuum mode as a function of the value of gravitational potential on the magnetic axis  $R_{{\rm M}}$. The gravitational potential is scaled with respect to the Alfvén speed on the magnetic axis.
Open with DEXTER


  \begin{figure}
\par\begin{tabular}{ll}
\includegraphics[width=5.5cm,clip]{6979-...
...
\includegraphics[width=5.5cm,clip]{6979-05c.ps} &
\end{tabular} \end{figure} Figure 5: The Brunt-Väisälä frequency projected on a flux surface for the points a), b), and c) shown in Fig. 4 shown as a contour plot in a poloidal cross-section.
Open with DEXTER

Next, we investigate the influence of the external gravitational potential on the overstable modes presented in the previous figure. This is done by varying the mass of the central object GM*. The result is shown in Fig. 4, where the most overstable mode of the MHD continua is plotted as function of the gravitational potential on the magnetic axis  $R_{{\rm M}}$. Notice that each point in this plot corresponds to a different equilibrium calculated with FINESSE. PHOENIX is used to analyse the stability of the computed equilibrium. The gravitational potential is stabilizing from 0 to -0.01 and destabilizing from -0.01 to lower values. The zero value corresponds to a toroidally rotating, non-gravitating "tokamak'' plasma. Its instability is in agreement with the results by van der Holst et al. (2000a). Figure 5 shows the projected Brunt-Väisälä frequency for the points (a), (b), and (c) indicated in Fig. 4. Plot (a) shows that the minimum value of this frequency is reached inside the plasma, while plot (c) shows that this value is reached at the edge of the accretion disk. Point (b) is the transition from the minimum value inside to the edge of the disk as seen in plot (b).

6.2 Thin accretion disks with non-constant density

To model a more realistic accretion disk equilibrium, we allow a density gradient across the poloidal magnetic surfaces. All other flux functions are kept as in the previous case, i.e.

 
                                     $\displaystyle I^{2}(\psi)$ = $\displaystyle A (1 - 0.0385\psi + 0.02\psi^{2} + 0.00045\psi^{3}), \quad \rho(\psi) = 1 - 0.4\psi,$  
$\displaystyle p_{0}(\psi)$ = $\displaystyle AB(1 - 0.9\psi), \quad \Omega(\psi) = C (1 - 0.9\psi^{2}),$ (59)

where the coefficients A, B, C are given by 91.5, 0.01, and 0.1. Also in this case, the poloidal cross-section is circular, the inverse aspect ratio  $\epsilon =0.1$, and the scaled mass of the central object is GM*=1. The two-dimensional pressure and plasma beta profile are shown in Fig. 6. This profile shows that the pressure reaches its maximum within the disk, instead of at the inner boundary of the disk, as in the case of a constant density. This difference can be explained as follows. Here, also the term  $p_{0} \Phi(R,Z)/T_{\rho} = \rho \Phi(R,Z)$ dominates in the pressure Eq. (13), but, due to the small inverse aspect ratio, the variation in the gravitational potential is small compared to the one in the density. This implies that the density, which is a flux function, mainly determines the shape of the dominant term in the Eq. (13). The temperature (not shown) decreases monotonically outwards. The range of plasma beta  $\beta=[0.082 , 0.158]$ and its maximum value is reached within the accretion disk. These low values for the plasma beta imply that the disk is strongly magnetized. Furthermore, the range of the ratio  $v_{\varphi}/v_{{\rm Kepler}}=[0.085 , 1.044]$.
  \begin{figure}
\par\includegraphics[width=6.4cm,clip]{6979-06.ps} \end{figure} Figure 6: The pressure (gray-scale) and plasma beta  $\beta = 2p/B^{2}$ (contours) distribution in a poloidal cross-section for an accretion disk with density  $\rho \propto (1-0.4\psi )$ and an inverse aspect ratio  $\epsilon =0.1$.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=6cm,clip]{6979-07a.ps}\includegraphics[width=6cm,clip]{6979-07b.ps}\end{figure} Figure 7: For the accretion disk equilibrium shown in Fig. 6, a) the real parts of the sub-spectrum of MHD continua as function of the radial flux coordinate  $s \equiv \sqrt {\psi }$ and b) sub-spectrum of the MHD continua in the complex plane are shown for toroidal mode number n=-1 and poloidal mode numbers m=[-3,5]. For these mode numbers, no unstable or overstable modes are found.
Open with DEXTER

For this case, the continuous spectrum has been computed for toroidal mode number n=1 and poloidal mode numbers m=[-3,5]. The results are shown in Fig. 7, which contains two plots. In the left one, the real part of the eigenvalues is plotted against the "radial'' coordinate  $s \equiv \sqrt {\psi }$, while the right one shows the eigenvalues in the complex plane. For this case, there are neither purely exponential unstable or overstable modes. Due to the presence of toroidal flow, gravity and a non-constant density, it is possible to create new or wider gaps between the different MHD continuum branches, with avoided crossings as a result of poloidal mode couplings. This has been shown by van der Holst et al. (2000b) for tokamak equilibria with toroidal flow. In these gaps, there may exist global modes similar to the Toroidal Flow-induced Alfvén Eigenmode (TFAE) for tokamaks found by van der Holst et al. (2000a) or to the Stratification-induced Alfvén Eigenmode (SAE) for coronal flux tubes found by Beliën et al. (1996). We will investigate this possibility for novel global modes in accretion disk context in future work.

6.3 Thick accretion disk


  \begin{figure}
\par\includegraphics[width=6cm,clip]{6979-08.ps} \end{figure} Figure 8: The pressure (gray-scale) and plasma beta  $\beta = 2p/B^{2}$ (contours) for a thick accretion disk with density  $\rho \propto (1-0.4\psi )$ and an inverse aspect ratio  $\epsilon =0.5$.
Open with DEXTER

Next, we increase the inverse aspect ratio $\epsilon$ from 0.1 to 0.5, while keeping all other quantities the same. Again, the poloidal cross-section is chosen circular. The resulting pressure profile is presented in Fig. 8. Here, The pressure decreases monotonically outward, similar to the equilibrium with constant density. Again, the term  $\rho \Phi(R,Z)$ dominates in the pressure Eq. (13). The inverse aspect ratio is now not small so that both the gravitational potential and the density determines the shape of the pressure profile. Also in this case, the temperature decreases in the outward direction. The plasma beta $\beta$ of this accretion disk ranges from 0.237 up to 0.970 and reaches its maximum within the disk. The deviation from a Keplerian disk is expressed by the ratio  $v_{\phi}/v_{{\rm Kepler}} = [0.007, 0.282]$. Thus the disk completely rotates sub-Keplerian.


  \begin{figure}
\par\begin{tabular}{l}
\includegraphics[height=0.2\textheight,cl...
...ludegraphics[height=0.2\textheight,clip]{6979-09c.ps} \end{tabular} \end{figure} Figure 9: For the thick disk equilibrium shown in Fig. 8, a) the real and b) imaginary parts of the sub-spectrum of the MHD continua as a function of the radial flux coordinate  $s \equiv \sqrt {\psi }$ and c) sub-spectrum of the MHD continua in the complex plane are shown for toroidal mode number n=-1 and poloidal mode numbers m=[-3,5].
Open with DEXTER

For this equilibrium, the continuous MHD spectrum is computed, again for toroidal mode number n=-1 and poloidal mode numbers m=[-3,5]. A part of the spectra is plotted in Fig. 9. Plot (a) shows the real part of the continuous MHD spectrum, while plot (c) shows the imaginary part as a function of the radial flux coordinate  $s \equiv \sqrt {\psi }$. Plot (c) contains the eigenfrequencies plotted in the complex plane. This shows the existence of overstable modes in a thick accretion disk due to the presence of toroidal flow and gravity. If one carefully examines the plot (a) one sees that one of the continuum branches splits in two at  $s \approx 0.87$, where the imaginary part becomes zero and at  $s \approx 0.91$ that these two continuum branches merge again, whereas the imaginary part becomes non-zero.

6.4 Maximum growth rate of the MHD continua of thin and thick accretion disks

The last two cases only differ by the value of the inverse aspect ratio $\epsilon$. We have systematically investigated the influence of this ratio on the existence of overstable modes by varying its value, while keeping all other equilibria quantities constant. In this parametric study, we essentially move the disk torus radially inwards while keeping the radius of the last closed flux surface a constant. We then gradually evolve from the equilibrium shown in Fig. 6 ( $\epsilon =0.1$) to the one shown in Fig. 8 ( $\epsilon =0.5$). The growth rate of the most overstable continuum mode is plotted as a function of the inverse aspect ratio $\epsilon$ in Fig. 10. Each point of the figure represents an equilibrium computed by FINESSE, from which the MHD continua is computed by PHOENIX. In this figure, one can clearly see that no overstable modes are present in these disks for  $\epsilon=[0.10 , 0.24]$. From  $\epsilon=0.24$ up to 0.32, there is a strong increase of the growth rate of the most overstable mode, while from $\epsilon=0.32$ up to at least 0.50 the growth rate remains approximately the same.

  
7 Conclusions

We have presented the equations describing an axisymmetric MHD accretion disk equilibrium and used FINESSE to find numerical solutions. We have considered thin as well as thick disks with a plasma beta of the order 0.1 where the largest part of the disk rotates sub-Keplerian. The numerical solutions show that the pressure profile can only decrease in the outward direction if the gravitational potential dominates in the pressure Eq. (13).


  \begin{figure}
\par\includegraphics[width=9cm,clip]{6979-10.ps} \end{figure} Figure 10: Growth rate of the most overstable continuum mode as a function of the inverse aspect ratio  $\epsilon = a / R_{0}$.
Open with DEXTER

The stability of the disk equilibria have been analysed theoretically as well as numerically. From theory, a stability criterion has been derived which closely resembles the Schwarzschild criterion. It is shown that in a disk where the density is a flux function, the MHD continua can turn overstable due to the presence of toroidal rotation and gravity. This instability is not possible for disks where the temperature or the entropy is a flux function. The numerical results support our theoretical conclusion. The effect of the gravity on the overstable modes can either be stabilizing or destabilizing. In the case of non-constant density the inverse aspect ratio $\epsilon$ has been varied. This study shows that the continua turn overstable if the inverse aspect ratio is larger than a certain critical value, approximately 0.24. These localised instabilities will lead to small-scale, turbulent dynamics when the equilibria are used in initial conditions for non-linear simulations. This dynamics is unrelated the usual magneto-rotational instability. Indeed, these instabilities, as well as the poloidal flow driven unstable continuum modes investigated by Goedbloed et al. (2004a), provide a new route to MHD turbulence in accretion disks. Furthermore, due to the presence of gravity, it is possible to widen or create new gaps in the MHD continua. It is possible to show this theoretically by expanding the eigenfunctions in poloidal harmonics. In these gaps new types of global ideal MHD eigenmodes can exist. The investigation of these global eigenmodes is left for future work.

Acknowledgements
This work was carried out within the framework of the European Fusion Programme, supported by the European Communities under contract of the Association EURATOM/FOM. Views and opinions expressed herein do not necessarily reflect those of the European Commission. 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)''.

References

 

Copyright ESO 2007