Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany

Abstract
We perform global three-dimensional simulations of accretion disks integrating the compressible, non-viscous, but diffusive MHD equations. The disk is supposed to be isothermal. We make use of the ZEUS-3D code integrating the MHD equations and added magnetic diffusivity. We measure the efficiency of the angular-momentum transport. Various model simulations delivered transport parameters of $\alpha_{\rm SS}=0.01$ to 0.05 which are consistent with several local numerical investigations. Two of the models reach a highly turbulent state at which $\alpha _{\rm SS}$ is of the order of 0.1. After a certain stage of saturating of the turbulence, Reynolds stress is found to be negative (inward transport) in many of the models, whereas Maxwell stresses dominate and deliver a positive (outward) total transport. Several of the models yield strongly fluctuating Reynolds stresses, while Maxwell stresses are smooth and always transport outwards. Dynamo action is found in the accretion disk simulations. A positive dynamo- $\alpha$ is indicated in the northern hemisphere of the most prominent run, coming along with negative kinetic and current helicities (all having the opposite sign on the southern side). The dipolar structure of the magnetic field is maintained throughout the simulations, although indication for a decay of antisymmetry is found. The simulations covered relatively thick disks, and results of thin-disk dynamo models showing quadrupolar fields may not be compatible with the results presented here.

Key words: accretion disks - instabilities - magnetic fields - MHD - turbulence

Accretion processes in astrophysical disks lead to enormous luminosity and huge changes in disk structure during astronomically short times. Efficient transport mechanisms are necessary to achieve such short time-scales. Anisotropic turbulence appears to be a major physical condition to provide astrophysical disks with strong transport. As these disks generally exhibit increasing specific angular momentum towards larger radii and thus fulfill the stability criterion of Rayleigh, they do not give rise to an instability leading to turbulence by themselves. Searches for instabilities in disks with rotation profiles similar to a Keplerian one unveiled several ways to turbulence being more or less favorable with respect to their prerequisites for the disk configuration. Gravitational instability needs the disk to be either cool or massive. Nonlinear and nonaxisymmetric perturbations require a severe additional perturber near the disk; conditions for instability in a purely hydrodynamical disk were derived by Dubrulle (1993). Hydrodynamic instability essentially comes down to violating the Rayleigh criterion saying that a rotation profile with an increasing specific angular momentum with radius is hydrodynamically stable, that is for $\partial l^2/\partial r>0$ . Since the typical length of perturbations caused be external forces is supposed to be very large, the required amplitude of the perturbations has to be considerable, too, in order to violate the Rayleigh criterion locally. If the length scale of the perturbation is comparable to the radius, a strong alteration of the Keplerian velocity of several per cent is needed. Convection was shown to deliver either negligible transport (Stone & Balbus 1996) or inward angular momentum transport (Kley et al. 1993).

The requirements for the magnetic shear-flow instability (Balbus & Hawley 1991) do match astrophysical conditions in accretion disks in many configurations. All it needs is a radially decreasing angular velocity and a weak magnetic field threading the rotating object. It can even be shown that the temperature range applicable to the magnetic shear-flow concept is very broad; even very small ionization fractions are sufficient to magnetize a disk in many cases (Balbus & Hawley 1998).

First numerical approaches to the magnetic shear-flow instability tackled the local problem; the linear analysis as described by Balbus & Hawley (1991) were immediately followed by 2D simulations (Hawley & Balbus 1991) of a small box-shaped domain which was cut out of the disk. These computations confirmed the relation between magnetic-field strength and wavenumber derived from the linear analysis earlier, they showed that the maximum growth rate of any wavenumber is independent from the field strength and that the system is capable of transporting angular momentum. Because of being restricted to axisymmetric problems, they could not provide self-sustained turbulence which needs dynamo-action, and the slow decay of the turbulence is an indirect effect of the Cowling theorem. Improved computations dealt with a three-dimensional even though local configuration, and particular care was taken for the radial boundary conditions which are not simply periodic, but account for the shear due to differential rotation (Hawley et al. 1995). These comprehensive computations indeed resulted in magnetically sustained turbulence whose anisotropy causes efficient outward angular momentum transport. This work was followed up by stratified models (Stone et al. 1996) covering more than 50 orbital periods of the box cut-out. As the computational domain covered 2 density scale heights, this work was a first step towards global simulations, followed up by similar approaches such as Ziegler & Rüdiger (2000a,b).

Linear studies of global configurations of disks threaded by magnetic fields in various directions were carried out. Curry & Pudritz (1995) investigated the stability for vertical and azimuthal fields threading the disk. They found in detail that the actual initial field geometry does not strongly depend on field topology as was suggested by the numerical simulations of Hawley & Balbus (1991). Rüdiger et al. (1999) particularly addressed the angular momentum transport in their linear study. The first nonlinear global approach by numerical simulations with full azimuthal coverage was presented by Armitage (1998) who omitted the density stratification for the sake of a large radial extent of the disk. A global approach with stratification was then presented by Hawley (2000) following the evolution of a thick torus under the influence of an external magnetic field threading parts of the computational domain. The magnetic shear-flow instability was found to set in quickly causing enough turbulent viscosity to soon form a Keplerian velocity profile. All the above mentioned numerical studies integrate the ideal MHD equations omitting magnetic diffusivity. Local simulations including diffusivity were performed by Fleming et al. (2000) which proved the onset of instability even for low conductivity.

The full understanding of accretion disks implies self-excited dynamo action as well as angular momentum transport. Can positive angular momentum transport be coupled with a suitable kinetic helicity providing the expected dynamo action according to the $\alpha$ -effect principle? As the kinetic helicity in stratified, rotating disks is expected to be negative, a wrong sign (positive) would follow for the dynamo- $\alpha$ -effect according to the conventional $\alpha$ -theories. We present global three-dimensional diffusive simulations and study the angular momentum transport as well as dynamo action and the sign of $\alpha^{\rm dyn}$ as a consequence of the correlation with the flow.

2 The simulation setup

The computations presented here make use of the ZEUS-3D code developed for astrophysical problems of magnetohydrodynamics (see the key papers Stone & Norman 1992a,b, Stone et al. 1992 for numerical; Clarke et al. 1994 for technical details). We use cylindrical coordinates and computational domains covering various annuli with a a vertical extension of z=-1 to +1 and the full azimuthal range of $\phi=0$ to $2\pi$ . See Table 1 for a list of models presented in this Paper. The notation $\Delta z$ , $\Delta r$ , and $\Delta \phi$ will refer to the extensions of the computational domain in the vertical, radial and azimuthal directions, respectively. In this approach, we assume an isothermal disk to save computation time on the energy equation. The remaining system for integration is

$\displaystyle % \frac{\partial\rho}{\partial t}+\mathop{\rm div}\nolimits(\rho {\vec u})=0$			(1)
$\displaystyle \frac{\partial{\rho \vec u}}{\partial t}+\mathop{\rm div}\nolimit... ...d}\nolimits p-\rho \mathop{\rm grad}\nolimits\Phi+{\vec J}\times {\vec B}+\dots$			(2)
$\displaystyle \frac{\partial {\vec B}}{\partial t}=\mathop{\rm curl}\nolimits({\vec u} \times {\vec B}) +\eta\triangle \vec B,$			(3)

Table 1: Model parameter used for this study. $T_{\rm orb}$ is the revolution period of the disk at the inner edge. The z-range is always from -1.0 to +1.0; the $\phi$ -range always from 0 to $2\pi$ . If accretion boundary conditions occur, they refer to the outer boundary; the inner boundary is "outflow'' then.
Model Resolution (z, r, $\phi$ ) Radial boundary condition r-range $T_{\rm orb}$ $\eta$

Ia $31\times61\times351$ inner outflow 5.0-6.0 0.222 0.001

Ib $31\times61\times351$ all outflow 5.0-6.0 0.222 0.001

II $31\times61\times351$ accretion $u_{\rm in}=-0.001c_{\rm ac}$ 4.0-6.0 0.159 0.001

III $31\times61\times351$ accretion $u_{\rm in}=-0.01c_{\rm ac}$ 4.0-6.0 0.159 0.001

IV $31\times31\times351$ accretion $u_{\rm in}=-0.001c_{\rm ac}$ 5.0-6.0 0.222 0.001

V $31\times61\times351$ Gaussian accretion $u_{\rm in}=-0.001c_{\rm ac}$ 4.0-6.0 0.159 0.01

VI $31\times61\times351$ Gaussian accretion $u_{\rm in}=-0.01c_{\rm ac}$ 4.0-6.0 0.159 0.01

VII $31\times31\times151$ Gaussian accretion $u_{\rm in}=-0.001c_{\rm ac}$ 4.0-6.0 0.159 0.01

VIII $31\times61\times351$ accretion $u_{\rm in}=-0.001c_{\rm ac}$ 3.0-7.0 0.103 0.001

**Table 1:** Model parameter used for this study. $T_{\rm orb}$ is the revolution period of the disk at the inner edge. The z-range is always from -1.0 to +1.0; the $\phi$ -range always from 0 to $2\pi$ . If accretion boundary conditions occur, they refer to the outer boundary; the inner boundary is "outflow'' then.
Model	Resolution (z, r, $\phi$ )	Radial boundary condition	r-range	$T_{\rm orb}$	$\eta$
Ia	$31\times61\times351$	inner outflow	5.0-6.0	0.222	0.001
Ib	$31\times61\times351$	all outflow	5.0-6.0	0.222	0.001
II	$31\times61\times351$	accretion $u_{\rm in}=-0.001c_{\rm ac}$	4.0-6.0	0.159	0.001
III	$31\times61\times351$	accretion $u_{\rm in}=-0.01c_{\rm ac}$	4.0-6.0	0.159	0.001
IV	$31\times31\times351$	accretion $u_{\rm in}=-0.001c_{\rm ac}$	5.0-6.0	0.222	0.001
V	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.001c_{\rm ac}$	4.0-6.0	0.159	0.01
VI	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.01c_{\rm ac}$	4.0-6.0	0.159	0.01
VII	$31\times31\times151$	Gaussian accretion $u_{\rm in}=-0.001c_{\rm ac}$	4.0-6.0	0.159	0.01
VIII	$31\times61\times351$	accretion $u_{\rm in}=-0.001c_{\rm ac}$	3.0-7.0	0.103	0.001

The gravitational potential is spherically symmetric, whence the z-component of the gravitation is retained within the disk. We therefore obtain a density stratification unlike the computations by Armitage (1998) who omits the z-component of the gravitation and applies periodic boundary conditions for the upper and lower boundaries. His approach allows much larger radial extents of the computational domain as the numerical restrictions from large density contrasts do not emerge.

The sound speed is always $c_{\rm ac}=10$ which is roughly 7% of the average Keplerian velocity in the simulated ring. The initial density distribution is Gaussian with r-dependent density scale-height H, thus $\rho(r,z)= \rho_{\rm c}\exp(-z^2/H(r)^2)$ with a density in the equatorial plane of $\rho_{\rm c}=1000$ . The equilibrium hydrodynamic model settles such that the z-extent of the domain covers 3Hat r=4 and 1.5H at r=6.

The magnetic field threads the disk vertically with a mere z-component. We tried a homogeneous initial field $B_z(r)={\rm const}$ for the z-boundaries as well as a field

The initial velocity field is a merely Keplerian motion following $u_\phi=u_{\rm K}=\sqrt{GM/r}$ , where G is the gravitational constant and M is the central mass which is 10⁵ in our computational units. Times are henceforth measured in orbital periods which convert by $T_{\rm orb}$ as given in Table 1 from the arbitrary units of the code.

The number of cells in each coordinate direction was typically $31\times61\times351$ for the z-, r-, and $\phi$ -directions. The aspect ratio of the grid cells is not unity. Finest resolution is achieved in radial direction, lowest in azimuthal direction: $\delta z/\delta r=2$ , $r\delta\phi/\delta r=2.18$ for the most often used inner radius r=4 and $r\delta\phi/\delta r=3.28$ for r=6 (Models II, III, V, VI, and VII), where $\delta r$ refers to the mesh width in radial direction and so on. Tests with up to 621 cells in $\phi$ -direction have been carried out, but the increased computation time does not allow reasonable periods to be covered by the simulation.

The fastest growing mode of the magneto-rotational instability has a wavelength of $\lambda_{\rm inst}=2\pi\sqrt{16/15}\,u_{\rm A}/\Omega$ . With the angular velocity of $\Omega=28.2$ for r=5 which is the middle of most of the radial ranges used, we obtain wavelengths between 0.4 and 0.9 with the above given Alfvén velocities of the initial field strength. These wavelengths are upper limits as they refer to the sites of maximum field and minimum density which need not coincide necessarily. Yet, they are well resolved by the computation domain covering z-ranges of 2.0 and two different r-intervals of 1.0 and 2.0.

The integration of the magnetic fields was done with the evolution of the electromotive forces by the Constraint Transport scheme (CT, Evans & Hawley 1988) which preserves a divergence-free configuration. The advection scheme of velocity and density involves the second-order van-Leer interpolation.

Two terms adding a source of viscosity to the hydrodynamic equations extend the above Navier-Stokes equation as indicated by dots in Eq. (2). The von Neumann-Richtmyer artificial viscosity depends on the square of velocity gradients and acts only on decreasing velocity in the direction of propagation, i.e. compression. The strength of this term is denoted by $q_{\rm con}$ in the implementation of the ZEUS-3D code. The second term depends linearly on the velocity gradient and acts on both compression and expansion; the strength is represented by $q_{\rm lin}$ . Tests with a stable Keplerian flow found the choice of $q_{\rm con}=0.1$ and $q_{\rm lin}=1.0$ to be suitable. These values were used throughout all the runs described here. The Courant number determining the "safety'' of a certain maximum allowed time step derived from the velocities, magnetic fields, and the artificial viscosities, is set to 0.5.

Occasionally, pockets of extremely low density may emerge coupled with extraordinarily high speeds exceeding the Keplerian velocity. We suppress the existence of such pockets by a mass replenishment as soon as the density of a certain cell drops below 10^-4 the central density. The actual mass being thus added to the total disk mass is found to be negligible. The mass replenishment thus acts like a lower limit for the time step.

2.1 Boundary conditions

The vertical boundary condition on the $z={\rm const}$ -faces of our computational domain are closed for the flow, i.e. u_z=0. The boundary condition is stress-free in the sense of $\partial u_r/\partial z=0$ and $\partial u_\phi/\partial z=0$ . The magnetic field has to fulfill the "opposite'' boundary condition as it is only allowed to cross the boundary perpendicularly, which is actually implemented affecting the electromotive force ${\vec{\cal E}}=\vec u\times \vec B - \eta \mathop{\rm curl}\nolimits\vec B$ such that ${\cal E}_z=0$ , $\partial {\cal E}_{r}/\partial z=0$ , and $\partial {\cal E}_\phi/\partial z=0$ .

Computations were carried out with two choices of radial boundary conditions regarding the velocity. The first does not allow flow through the top, bottom and outer radial boundaries. The inner boundary was chosen to be "outflow'', that is, the physical arrays are constantly extrapolated into the ghost zones of the computational domain. The only exception is an inflow, in our case a u_r > 0 at the inner boundary. This velocity component is reset to zero then. All radial velocities u_r<0 are copied to the ghost zones beyond the computational domain. The same holds for the density and the magnetic-field components. We will henceforth refer to these models as Model Ia (only inner boundary has "outflow'') and Ib (inner and outer radial and lower and upper vertical boundaries have "outflow'').

Since the outflow condition is likely to empty the disk on the viscous time-scale, the second choice tries to account for the accreted matter and feeds the disk at the outer radial boundary. We sum up the mass loss rate at the inner boundary $\dot M=r_{\rm min}\Delta\phi\Delta z \sum_i\sum_k \rho_{ik} (u_{r})_{ik}$ for the total mass loss. The indices i and k run along the vertical and azimuthal directions, respectively.

A constant, small inflow velocity $u_{\rm in}$ is assumed for the outer boundary. Accordingly, the density at the outer boundary is determined to account for the mass flux at the inner boundary. Since we use the average influx, its z-dependence (nor even the $\phi$ -dependence) is not used for a z- (nor even $\phi$ -) dependent inflow from far radii. The boudary conditions for the outer radial face are thus $u_{r}=u_{\rm in}$ and $\rho=\dot M/(A_{\rm out}u_{\rm in})$ where $A_{\rm out}$ is the surface of the outer boundary. The magnetic-field conditions are B_r=0, $\partial B_z/\partial r$ , and $\partial B_\phi/\partial r=0$ . Note that this is not consistent with the actual inflow of matter, but the $u_{\rm in}$ is extremely small, less than $10^{-4}u_{\rm K}(r_{\rm o})$ for the slow-inflow models and less than $10^{-3}u_{\rm K}(r_{\rm o})$ for the fast-infall models (see third column of Table 1). Upon expanding $\mathop{\rm curl}\nolimits(\vec u\times \vec B)$ we find $B_\phi\partial u_{r}/\partial r$ to be the relevant term generating a contribution to B_z. The field in the last three zones in Fig. 2 is an effect of this.

We tested inflow conditions with homogeneous and Gaussian and density distribution. It was supposed that an average influx of mass will establish according to the current viscosity (either numerical or turbulent) in the disk. These models are referred to as Model II-VIII in the following.

The azimuthal velocity, $u_\phi$ , is extrapolated into the radial boundary zones by a power law r^-1.5 based on the last zone of the computational domain. The velocities are thus not Keplerian, they just follow the same radial power law. Since we cover the full azimuthal range in $\phi$ , the boundaries at the $\phi={\rm const}$ surfaces are periodic, naturally.

$\begin{figure} \par {\epsfig{file=h2593f1.ps, width=8.8cm} } \par \end{figure}$	Figure 1: Vertical cut through the disk of Model II after 6.3 orbital revolutions, excluding the high-velocity coronal component of the disk. The shading represents the azimuthal component, the arrows are the velocity vectors projected onto the (r,z)-plane.
Open with DEXTER

$\begin{figure} \par {\epsfig{file=h2593f2.ps, width=8.8cm} } \par \end{figure}$	Figure 2: Vertical cut through the disk of Model II after 6.3 orbital revolutions with an azimuthal-component shading and poloidal magnetic-field vectors.
Open with DEXTER

3 Patterns, transport, spectra

3.1 Glancing at the solutions

Typical vertical slices from a run of Model II are shown in Figs. 1 and 2. The central mass is located on the left side, outside the computational domain. The shadings represent the azimuthal components of the magnetic and velocity fields, while the arrows are the (r,z)-projections of the actual vectors. Shaded areas near the equatorial plane of the disk show the typical Keplerian velocity profile. However, at high altitudes above the equator, a strong super-rotation emerges coupled with relatively high infall velocities. (These high radial velocities have been cut in Fig. 1 to better show the velocities in the disk.)

The computational domain thus covers part of the corona of the disk where the density falls below 10^-3 times the equatorial density. The infall velocity in the corona is of the same order of magnitude as the Keplerian velocity, whereas the non-orbital velocities within the disk are much smaller. The actual flow in the disk as shown in Fig. 1 is limited by $0.35 c_{\rm ac}$ .

After saturation of the turbulence, the temporal behavior of $\alpha _{\rm SS}$ changes. The fairly smooth function turns into an oscillatory behavior. The typical feature of this moment is the zero-crossing of the Reynolds stress which is essentially negative all through the saturation parts of the simulations. It is expected that the infall of matter at the outer radial boundary tends to perform relaxation motions, in particular for Models II-IV where matter falls in homogeneously. The exclusion of the outermost six cell planes from the computation of the spatially average $\alpha _{\rm SS}$ does not alter the result (the actual relaxation motions are observed only in the outermost three cell planes). Relaxation motions will consist mainly of u_z components which do not affect radial angular-momentum transport directly.

A problem connected with the homogeneous feeding of the disk from the outer boundary is the appearance of a "density front''. This is the reason why the density at the outer boundary is highest which may look odd in the context of an accretion disk. This incoming density enhancement can be partly leveled by the presumed infall velocity at the outer boundary.

3.2 Angular momentum transport

According to the standard model of a turbulent disk as essentially worked out by Shakura & Sunyaev (1973), the $r\phi$ -element of the stress tensor of velocity and magnetic-field fluctuations,

$\begin{figure} \par {\epsfig{file=h2593f3.ps, width=8.8cm} } \par \end{figure}$	Figure 3: Angular momentum transport by Reynolds (dashed) and Maxwell (dotted) stresses as derived from a run of Model Ia, that is with outflow boundary condition at the inner radial face and a closed boundary for flow and magnetic fields at the outer radial face. The total transport (thick solid line) is positive and of the order of 0.0005.
Open with DEXTER

$\begin{figure} \par {\epsfig{file=h2593f5.ps, width=8.8cm} } \par \end{figure}$	Figure 5: Angular momentum transport by Reynolds (dashed) and Maxwell (dotted) stresses as derived from a run of Model II. The total transport (thick solid line) is positive and of the order of 0.01.
Open with DEXTER

$\begin{figure} \par {\epsfig{file=h2593f7.ps, width=8.8cm} } \par \end{figure}$	Figure 7: Angular momentum transport by Reynolds (dashed) and Maxwell (dotted) stresses as derived from a run of Model III. The homogeneous inflow velocity is ten times higher in this simulation that in Fig. 5.
Open with DEXTER

The time series of the outflow model Ia is shown in Fig. 3. Model Ia in particular shows a short-lived occurrence of strong angular-momentum transport; a similar but less pronounced behavior is found in Model Ib. This period is the only time in both Models Ia and Ib when Reynolds stresses dominate the Maxwell stresses. A similarity to the transient channel solution as described by Hawley et al. (1995) is likely, although that was found for uniform initial B_z.

Time series for long runs of Models II and Model III are shown in Figs. 5 and 7. Generally, the accretion model shows stronger angular momentum transport than the the simpler inner-outflow Model Ia. Models Ia and Ib run out of matter after about 8 orbital revolutions with enabled magnetic fields, and the $\alpha _{\rm SS}$ values are no longer meaningful.

A striking feature of the accretion Models II and III is the change of sign in the Reynolds stress after about 4-6 orbits. We find that, in a long run, the outward angular momentum transport is a sole consequence of the magnetic stresses. Additionally, the Reynolds part shows strong amplitude variations unlike the Maxwell stress.

Figures 4 and 6 show the respective temporal evolution of the fluctuative kinetic energies. For the z- and r-components, the zero mode is subtracted such as $\langle u_z'^2\rangle=\langle u_z^2\rangle-\langle u_z\rangle^2$ ; the azimuthal fluctuative energy is found by $\langle u_\phi'^2\rangle=\langle (u_\phi - u_{\rm K})^2\rangle$ . All three models show a clear energy "sorting'' of the z-, r-, and $\phi$ -components in this order.

A particularly long run was performed with a computational domain only half as large as Models II and III and is denoted by IV. The initial and boundary conditions are like in Model II. The long-term behavior is not too similar to that of Model II with very strong oscillations of Reynolds versus Maxwell stresses. A behavior comparable with Model II may be spotted during the first 8 orbital periods. After that time, the energy sorting of the radial and vertical components is given up. All models from Ia to IV appear to be not fully saturated with respect to the ever-growing kinetic energies.

We should note that a homogeneous infall might be suitable for a realistic disk in the sense that the infalling matter does not "know'' the structure of the disk. Nevertheless, numerical problems may occur when the incoming homogeneous "matter front'' meets the near-Gaussian density structure. Models V and VI apply an inflow with pre-defined Gaussian infall profile over z. The respective $\alpha _{\rm SS}$ -profiles are shown in Figs. 8 and 10, the first having an infall velocity of $u_{\rm in}=-0.01$ , the second $u_{\rm in}=-0.1$ . The initial behavior is similar to the respective Models II and III including the change of sign of the Reynolds stress after a couple of orbits. After as many as 10 orbital periods, the disk turns into a significantly different regime with enhanced outward transport in both models.

The respective kinetic energies are given in Figs. 9 and 11. In contrast to the above studied models, energies do not develop in an ever-growing way, but appear to fluctuate about an average. The only exception is the regime change near $t=12.5\,T_{\rm orb}$ (Model V) and $t=10.5\,T_{\rm orb}$ (Model VI) when the average (essentially of $\langle u_z'^2\rangle$ ) changes, but the stationary appearance remains.

$\begin{figure} \par {\epsfig{file=h2593f8.ps, width=8.8cm} } \par \end{figure}$	Figure 8: Angular momentum transport by Reynolds (dashed) and Maxwell (dotted) stresses as derived from a run of Model V. Density is not falling in homogeneously, but with a predefined Gaussian vertical distribution; the (homogeneous) inflow velocity is $-0.001c_{\rm ac}$ .
Open with DEXTER

$\begin{figure} \par {\epsfig{file=h2593f10.ps, width=8.8cm} } \par \end{figure}$	Figure 10: Angular momentum transport by Reynolds (dashed) and Maxwell (dotted) stresses as derived from a run of Model VI. The inflow velocity is ten times higher than in model V, Fig. 8.
Open with DEXTER

Table 2: Results of the models given in Table 1. The duration is the total period of the run with enabled magnetic fields. The last column gives the period from which the $\alpha _{\rm SS}$ -average was obtained. Times are in units of orbital revolutions of the inner edge.
Mod. Resolution Duration $\langle\alpha_{\rm SS}\rangle$ during time

Ia $31\times61\times351$ 8.3 0.0003 4.0-8.9

Ib $31\times61\times351$ 12.2 0.01-0.02 3.5-9.0

II $31\times61\times351$ 14.7 0.014 3.6-14.6

III $31\times61\times351$ 9.7 0.041 2.8-9.7

IV $31\times31\times351$ 41.5 0.014 8.2-43.0

V $31\times61\times351$ 11.9 0.038 3.4-11.9

'' '' 0.12 12.4-15.9

VI $31\times61\times351$ 16.1 0.037 3.1-10.1

'' '' 0.14 11.1-16.1

VII $31\times31\times151$ 36.8 0.080 15.9-27.9

VIII $31\times61\times351$ 22.4 0.07 4.3-22.4

**Table 2:** Results of the models given in Table 1. The duration is the total period of the run with enabled magnetic fields. The last column gives the period from which the $\alpha _{\rm SS}$ -average was obtained. Times are in units of orbital revolutions of the inner edge.
Mod.	Resolution	Duration	$\langle\alpha_{\rm SS}\rangle$	during time
Ia	$31\times61\times351$	8.3	0.0003	4.0-8.9
Ib	$31\times61\times351$	12.2	0.01-0.02	3.5-9.0
II	$31\times61\times351$	14.7	0.014	3.6-14.6
III	$31\times61\times351$	9.7	0.041	2.8-9.7
IV	$31\times31\times351$	41.5	0.014	8.2-43.0
V	$31\times61\times351$	11.9	0.038	3.4-11.9
''	''		0.12	12.4-15.9
VI	$31\times61\times351$	16.1	0.037	3.1-10.1
''	''		0.14	11.1-16.1
VII	$31\times31\times151$	36.8	0.080	15.9-27.9
VIII	$31\times61\times351$	22.4	0.07	4.3-22.4

The distribution of $\alpha _{\rm SS}$ versus time and the z-direction is shown in Fig. 12 for Models V. Outward angular momentum transport is marked with white areas, inward transport by dark areas. Strong transport is observed at large distances from the equator shortly after the onset of the instability. Regions of outward transport migrate towards the equator. This is particularly marked in the Model V plot for which we found a development into a strong-transport regime (cf. Fig. 8). Positive angular momentum transport is not restricted to high altitudes anymore.

$\begin{figure} \par {\epsfig{file=h2593f12.ps, width=8.8cm} } \par \end{figure}$	Figure 12: Distribution of the total stress $\alpha _{\rm SS}$ versus time and vertical coordinate for Model V. Light areas represent inward transport, dark areas outward transport.
Open with DEXTER

The temporal development of the accretion rate measured at the inner radial boundary is given in Fig. 13. As we were free to choose a central disk density, the computation of the accretion rate in physical units will not provide much information, whence the arbitrary units in Fig. 13. Choosing the rough dimensions of an inner protostellar system of 100 AU radius for scaling the mass flux through our models, one would obtain an accretion rate of about $10^{-5}~M_\odot$ /yr for Model II shown in the figure.

3.3 Limits of the simulations

The application of accreting boundary conditions allows relatively long durations of the simulations at constant total mass. Yet, the influence of the inconsistent outer radial boundary becomes important after approximately 15 orbital periods. The major part of the mass concentrates near the outer boundary, while the rest of the disk appears to run empty with strong inward radial velocities of $u_{r}=2c_{\rm ac}$ to $3c_{\rm ac}$ . At this point, our adjustment of the inflow density to sustain the total mass does not match the actual process going on, and a strong pile-up of density at the inlow boundary will be observed. The disruption is particularly sudden for Models V, VI, and VII which are specific in that they have a higher magnetic diffusivity than the other models.

$\begin{figure} \par {\epsfig{file=h2593f14.ps, width=8.8cm} } \par \end{figure}$	Figure 14: Azimuthal Fourier decomposition of the magnetic-field components after 8.15 orbital periods in Model V. The dotted line shows a Kolmogorov power law with exponent -5/3.
Open with DEXTER

3.4 Spectral decomposition

Power spectra resulting from Fourier transforms of the magnetic field of Model V over the the azimuthal and radial directions are shown in Figs. 14 and 15, respectively. The azimuthal decomposition of the velocity field is shown in Fig. 16. The individual azimuthal spectra for all points in the (z,r)-plane were averaged, as were the points of the ( $\phi,z$ )-plane for the radial spectra. In the azimuthal spectra, we observe slight maxima in power at k=5 and k=8, followed by a strong, roughly power-law decline between $k\sim 20$ and $k\sim 100$ . The magnetic spectrum is steeper than the kinetic one. The underlying power-laws would be $E\propto k^{-3}$ to $E\propto k^{-4}$ which is significantly steeper than a Kolmogorov spectrum, and is similar to what Armitage (1998) found for a global simulation in the far-k range. The Kolmogorov spectrum of isotropic turbulence would exhibit a wavenumber exponent of -5/3; here we face MHD turbulence which implies an effect of Alfvén waves impeding the transfer of energy towards smaller scales. The energy spectrum of isotropic MHD turbulence is $E\propto k^{-3/2}$ . Contrasting with the azimuthal decomposition, the radial spectrum matches the Kolmogorov spectrum over the entire range of modes as can be seen in Fig. 15.

3.5 Comparison with other studies

Among stratified local simulations, we find a number of very similar results as regards the magnitude of $\alpha _{\rm SS}$ , such as in Stone et al. (1996) who obtained 0.01, Hawley et al. (1996) where an average of 0.016 was found, with typical values of 0.01. The largest box used delivered 0.02, the smallest, 0.003. Upon studying the dependence of the turbulence on the rotation law $\Omega\propto r^{-q}$ , Hawley et al. (1999) found 10^-2 for q=1.5. A most recent local-box simulation by Ziegler & Rüdiger (2000b) covering 288 orbital periods resulted in a space-time average of 0.015, where the Reynolds part was $2.8\times 10^{-3}$ and the Maxwell part was $1.2\times 10^{-2}$ . The stress was normalized with the averaged equatorial-plane pressure and thus provides a lower limit for $\alpha _{\rm SS}$ . The early stratified-box simulations by Brandenburg et al. (1995) provided $\alpha _{\rm SS}$ an order of magnitude lower between 0.001 and 0.005. All the models agree upon the dominance of Maxwell stresses over Reynolds stresses.

It is most interesting to note that other global simulations such as performed by Armitage (1998) or Hawley (2000) result in an order of magnitude higher values of 0.2 to 0.3. Contrasting with this, the global spherical, but unstratified simulation by Drecker (2000) gives again an $\alpha _{\rm SS}$ which is an order of magnitude lower than shown here, i.e. 10^-3. Those simulations alone included a physical viscosity. The simulations of Hawley (2000) - and with somewhat different objectives, those of Machida et al. (2000) - had open boundaries in radial and vertical directions. This fact limits the run-time of the models which are prone to running out of matter.

4 Dynamo effect

The amplification of the total magnetic energy by a factor of roughly 10³ from the initial magnetic-field perturbation to the state of strong angular-momentum transport indicates considerable dynamo action in the disk. In the following we will investigate the characteristics of magnetic-field generation in our simulations.

The local-box simulations by Hawley et al. (1996) impose severe constraints on the growth of an average magnetic field because of the periodicity of the boundary conditions. A non-zero term of the average magnetic-field growth is due to the shear and requires the presence of an average radial component, in the local nomenclature $\langle B_x\rangle$ . We follow their derivation of the averages by replacing the volume curl integral by a vector-product surface integral of (3). The average

$\displaystyle % \frac{\partial\langle B_z\rangle}{\partial t} = -\left\langle{\frac{\partial}{\partial r}u_{r} B_z}\right\rangle$			(13)
$\displaystyle \frac{\partial\langle B_{r}\rangle}{\partial t} = \phantom{-}\lef... ...angle{\frac{\partial}{\partial z}\frac{\partial B_z}{\partial r} }\right\rangle$			(14)
$\displaystyle \frac{\partial\langle B_\phi\rangle}{\partial t}= \phantom{-}\lef... ...}{\partial r}\frac{1}{r} \frac{\partial rB_\phi}{\partial r}}\right\rangle\cdot$			(15)

4.1 Sketching mean-field dynamos

Kinematic dynamos have been extensively examined in numerous previous publications concerning various types of geometries. A widely used mean-field approach circumnavigates the difficulties with resolving the small scales in a model. We address the applicability of the mean-field theory to accretion disk dynamos by extracting characteristic quantities from the simulations.

The mean-field dynamo applies the likely amplification of a magnetic field through a fluctuative electromotive force parallel to the actual large-scale magnetic field. This effect is generally expressed by the term " $\alpha$ -effect'' and can be thoroughly studied in e.g. Krause & Rädler (1981). In a model splitting large and small scales, the induction equation extends beyond the large-scale electromotive force such as

Second-order correlation approximation leads to a relation of $\alpha^{\rm dyn}$ with the helicity of the flow,

4.2 Symmetry of magnetic fields

Special attention is pid to the symmetry of the magnetic fields of the solutions. In general, the preference of a certain symmetry is a tool to check the consistency with mean-field dynamo models. Additionally, dipolar (antisymmetric) fields are supposed to support the launch of winds forming jets from the disk better than quadrupolar (symmetric) fields, with the latter rather providing closed field lines within the disk. The mean-field disk dynamo of Rekowski et al. (2000) produced dipolar magnetic fields for negative $\alpha^{\rm dyn}$ in the northern hemisphere.

A number of mean-field dynamo models in one to three dimensions delivered varying symmetry of the solutions with a major influence of the boundary conditions - vacuum versus perfectly conducting. Three-dimensional investigations of stability by Meinel et al. (1990) and Elstner et al. (1992) found quadrupolar solutions for the condition of low-conductivity surroundings of the disk, whereas perfectly conducting boundaries delivered dipolar solutions. A local model with periodic boundaries will thus not be able to answer the question of the final symmetry of the solutions. An ideal global model would comprise the entire investigated object and is less depending on the surroundings (which will be vacuum in very good approximation).

We should emphasize that the initial perturbation of the magnetic field had dipolar symmetry. Since the total flux through the vertical surfaces, however, vanishes, field configurations with either symmetry may arise from the evolution of the model. The dipolar initial perturbation may not have been fully re-organized as the diffusive time-scale

$\begin{figure} \par {\epsfig{file=h2593f17.ps, width=8.8cm} } \par \end{figure}$	Figure 17: Total $\alpha _{\rm SS}$ and parity for Model II with $\eta =0.001$ . A parity of -1 denotes fully antisymmetric. (dipolar) solutions.
Open with DEXTER

4.3 The dynamo of the disk model

Average kinetic helicities were computed based on a layer slightly below one density scale-height in both hemispheres. The average is taken in $\phi$ - and r-direction. There is a considerable scatter in the time series of the helicity, but temporal averages are all negative on the northern side for all models. The only exception is the northern part of Model II with a near-zero value. Negative sign also holds for the two test models with mixed-parity initial perturbation. Average current helicities are also negative on the northern side (and positive on the southern one) throughout almost all models. The only exceptions are the mixed-parity models which show a current helicity with opposite sign compared with the kinetic helicity in both hemispheres.

Since the angular momentum transport is dominated by magnetic stresses, we also ask about the energy in the flow fluctuations and the magnetic field. The temporal evolutions of these two quantities in Figs. 19 and 20 show a clear difference between the Models II and V. Kinetic-energy dominance holds for Model II, whereas magnetic dominance is found for Model V.

$\begin{figure} \par {\epsfig{file=h2593f19.ps, width=8cm} } \par \end{figure}$	Figure 19: Temporal evolution of kinetic (solid line) and magnetic (dashed line) energies of Model II with $\eta =0.001$ . The energy is measured as a radial and azimuthal total in a horizontal layer at z=+0.39 which is - measured in density scale-heights H - near 1.2H at the inner radial boundary and near 0.6H at the outer boundary, due to the variable disk height along radius.
Open with DEXTER

$\begin{figure} \par {\epsfig{file=h2593f20.ps, width=8cm} } \par \end{figure}$	Figure 20: Temporal evolution of kinetic (solid line) and magnetic (dashed line) energies of Model V with $\eta =0.01$ . The energy is measured as a radial and azimuthal total in a horizontal layer at z=+0.39.
Open with DEXTER

In order to evaluate the applicability of the $\alpha^{\rm dyn}$ approach, we compute the average toroidal field $\langle B_\phi\rangle$ and compare it with the actual average electromotive force ${\cal E}_\phi=(\langle\vec{u}\rangle\times \langle\vec{B}\rangle)_\phi$ derived directly from the simulated vector fields. A correlation between the two quantities may justify the $\alpha$ -effect for dynamo generation of magnetic fields. As a meaningful $\alpha^{\rm dyn}$ must change its sign at the equator, we plot the correlation for both hemispheres of the disk separately; the result is given in Fig. 22. As the temporal fluctuations are strong, we plot the time-averages of $B_\phi$ and $(\vec{u}\times \vec{B})_\phi$ in Fig. 21 only with their resulting sign. The quantities are also averaged in azimuthal direction and plotted in the (r,z)-plane. White areas denote a negative sign; the dominance of negative averaged electromotive forces indicates a positive $\alpha^{\rm dyn}$ in the northern hemisphere in accordance with Fig. 22.

The classical understanding of the Coriolis force giving preference to left-handed helicity on the northern hemisphere (whence positive $\alpha^{\rm dyn}$ ) appears not applicable for the correlation plots of Models II and III. The fluctuations were strong, and no meaningful sign of $\alpha^{\rm dyn}$ can be derived. However, Models Ia, V, and VI show a mostly negative averaged EMF for both hemispheres, while $\langle B_\phi\rangle$ changes its sign. An indication for positive $\alpha^{\rm dyn}$ is thus found from these simulations.

$\begin{figure} \par {\epsfig{file=h2593f22.ps, width=11cm} }\vspace{2mm} \par {... ...} }\vspace{2mm} \par {\epsfig{file=h2593f24.ps, width=11cm} } \par \end{figure}$	Figure 22: Correlation of the average electromotive force versus the mean $B_\phi$ , averaged for the northern and southern hemisphere separately. The quantities were derived, from top to bottom, from the runs of Models Ia, V, and VI. The sign of the correlation gives an indication for the sign of the dynamo $\alpha$ -effect.
Open with DEXTER

In fact, Models V and VI are those which show saturated kinetic energies. The other significant difference to the indefinite Models II and III is the 10 times higher diffusivity, $\eta =0.01$ . Notwithstanding, an indication for a negative $\alpha^{\rm dyn}$ was already found in simulations of a local box cut out of the disk by Brandenburg et al. (1995). The same sign is found by Ziegler & Rüdiger (2000b) from long runs of local simulations, although a clear influence of resistivity on the evolution was found there. Lowest magnetic Reynolds numbers are ${\rm Re}_{\rm m}=1450$ for Model V in this Paper, calculated with the velocity difference between inner and outer radial boundary, just as Ziegler & Rüdiger did for the local box.

These considerations are somewhat limited, since a straight-forward regression line in Fig. 22 will have a significant offset from the origin of the graph; a vanishing $\langle B_\phi\rangle$ does apparently not coincide with a vanishing $(\langle\vec{u}\rangle\times \langle\vec{B}\rangle)_\phi$ . The effect of the other field components is considered small though; while the average hemispherical $\langle B_\phi\rangle$ are of the order of 100, the $\langle B_{r}\rangle$ is of order 10 and $\langle B_z\rangle$ of order unity. A considerable contribution is suspected from the current density through the $\eta_{ijk}$ -term in Eq. (17), and it appears to be quite natural that an offset in the correlation graphs in Fig. 22 is found. For this reason, we did not plot such regression lines in the figure.

Additionally, the sole consideration of the induction equation in the "classical'' dynamo theory leads to the omission of the Lorentz force, which is not only essential for the onset of the magneto-rotational instability but for the maintenance of the turbulence and the transport of angular momentum as well. It is thus obvious that the neglect of the back-reaction of magnetic fields on the flow need not lead to representative models of Keplerian disks.

5 Conclusions

Three-dimensional global simulations have shown, that the magnetic shear-flow instability is a fast mechanism to generate a turbulent flow in a Keplerian disk. We presented models with a computational domain covering the full azimuthal range, accounting for accretion of matter. The Shakura-Sunyaev parameter $\alpha _{\rm SS}$ increases rapidly during the first revolutions after switching on the magnetic field and reaches 10^-2-10^-1 at this stage of the computations, though it appears not yet fully saturated in several of the models. Maxwell stresses exceed Reynolds stresses almost entirely. The latter undergo strong fluctuations and are partly negative.

Indications for a dynamo action in the disk are found, the corresponding dynamo- $\alpha^{\rm dyn}$ tending to be positive north of the equatorial plane and negative south of the equatorial plane. The generated magnetic fields may maintain the turbulence even if the external field ceases. The strong excitation of low-order azimuthal modes in the magnetic field is another promising fact for dynamo action with respect to the Cowling theorem.

5.1 The dynamo of Model V

For various reasons, a representative example for an accretion disk dynamo is provided by Model V (and nearly as well by the similar Model VI with higher inflow velocity). The time series is sufficiently long covering more than 18 revolution periods. The angular momentum transport reaches a saturated level of high efficiency during the last four orbital periods. The magnetic diffusivity is roughly two orders of magnitude higher than the numerical diffusivity. The striking dominance of one sign of the averaged EMF in Fig. 21 indicates a physically evolved, relevant state of the system.

The results connected with the outcomes of mean-field dynamo theory include the following facts: (i) A negative average toroidal magnetic field is found for the northern hemisphere, a positive on the southern one. The averaged EMF is negative in both hemisphere indicating a positive $\alpha^{\rm dyn}$ -effect on the northern side (negative on the southern side). (ii) The correlation plot of average toroidal field and EMF gives the same picture. (iii) Negative kinetic and current helicities on the northern hemisphere (positive on the southern one) are also consistent with a positive $\alpha^{\rm dyn}$ . (iv) The dipolar structure of the solution contrasts with the results from mean-field dynamo models of disks which yield quadrupolar solutions for positive $\alpha^{\rm dyn}$ . We have to add though that the disks of our Paper are not thin. Mean-field simulations by Covas et al. (1999) show a transition from quadrupolar to dipolar fields when the opening angle of the disk is enlarged, rather representing a torus than a disk.

The extension of the run-times of such simulations exceeding the order of diffusion times promises further results about dynamo action in accretion disks. The connection with mean-field concepts are most interesting as well as the implications for jet-launch models.

Global accretion disk simulations of magneto-rotational instability

1 Introduction