A&A 490, 501-514 (2008)
DOI: 10.1051/0004-6361:200810385

High accretion rates in magnetised Keplerian discs mediated by a Parker instability driven dynamo

A. Johansen1,2 - Y. Levin1


1 - Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
2 - Work partially done at Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

Received 13 June 2008 / Accepted 22 August 2008

Abstract
Hydromagnetic stresses in accretion discs have been the subject of intense theoretical research over the past one and a half decades. Most of the disc simulations have assumed a small initial magnetic field and studied the turbulence that arises from the magnetorotational instability. However, gaseous discs in galactic nuclei and in some binary systems are likely to have significant initial magnetisation. Motivated by this, we performed ideal magnetohydrodynamic simulations of strongly magnetised, vertically stratified discs in a Keplerian potential. Our initial equilibrium configuration, which has an azimuthal magnetic field in equipartion with thermal pressure, is unstable to the Parker instability. This leads to the expelling of magnetic field arcs, anchored in the midplane of the disc, to around five scale heights from the midplane. Transition to turbulence happens primarily through magnetorotational instability in the resulting vertical fields, although magnetorotational shear instability in the unperturbed azimuthal field plays a significant role as well, especially in the midplane where buoyancy is weak. High magnetic and hydrodynamical stresses arise, yielding an effective $\alpha $-value of around 0.1 in our highest resolution run. Azimuthal magnetic field expelled by magnetic buoyancy from the disc is continuously replenished by the stretching of a radial field created as gas parcels slide in the linear gravity field along inclined magnetic field lines. This dynamo process, where the bending of field lines by the Parker instability leads to re-creation of the azimuthal field, implies that highly magnetised discs are astrophysically viable and that they have high accretion rates.

Key words: accretion, accretion disks - Galaxy: center - instabilities - magnetohydrodynamics (MHD) - turbulence

1 Introduction

Since the seminal work of Balbus & Hawley (1991) it has been widely recognised that hydromagnetic stresses play a central role in the dynamics of Keplerian gaseous discs. Much of the subsequent theoretical work has been devoted to the study of the magnetic dynamo driven by a combination of magnetorotational instability (MRI), rotation, and Keplerian shear (Brandenburg et al. 1995; Hawley et al. 1996, see Brandenburg & Subramanian 2005, for an excellent review of astrophysical dynamo theory). Typically, numerical simulations of the MRI-driven dynamo begin with an initial zero net flux magnetic field with an associated pressure which is a small fraction of the thermal pressure. The idea is that the MRI-driven turbulence should increase the characteristic coherence length of the magnetic field, and grow its mean strength to a significant fraction of the equipartition value. The simulations, however, have had at best a mixed success in explaining high and persistent inflow rates observed in astrophysical accretion discs. Firstly, a number of recent numerical experiments have shown that the effectiveness of the dynamo depends in a critical way on the value magnetic Prandtl number, defined as the ratio of the collisional kinematic viscosity to the magnetic diffusivity (Lesur & Longaretti 2007; Fromang et al. 2007, following suggestive non-disc simulations of Schekochihin et al. 2005 and a careful study of previously published MRI-turbulence results by Pessah et al. 2007). The upshot of this work is that when the Prandtl number is significantly less than one, as is expected in most accretion discs (Balbus & Henri 2008), the dynamo seems to fail[*]. Secondly, in the zero net flux simulations where all parameters are chosen so that the dynamo works, the measured value of the Shakura-Sunyaev parameter $\alpha $ is typically of order 10-3 and at most 10-2 (e.g. Fromang & Nelson 2006; Brandenburg et al. 1995; Johansen & Klahr 2005; Sano et al. 2004). This is 1-2 orders of magnitude smaller than what is required to explain the high accretion rates in dwarf novae systems (King et al. 2007).

The weak initial magnetic fields assumed in almost all disc-MRI simulations may be unrealistically small for many astrophysical discs. We give here two examples:

1.
Extended accretion discs around central supermassive black holes (such as the ones traced by megamasers in nearby galaxies; see e.g. Vlemmings et al. 2007; Greenhill 2007) are fed by molecular material in the interstellar medium. Molecular clouds are known to have large scale superthermal magnetic fields, and thus the initial magnetic fields in AGN discs are likely to be comparable to or larger than the thermal equipartition values. The $\sim$2 pc molecular circumnuclear disc in our galactic centre is permeated by large scale equipartition magnetic fields (Hildebrand et al. 1990; Wardle & Konigl 1990) which are certain to play an important dynamical role in its subsequent evolution.

2.
Accretion discs in close binary systems, which are fed by a Roche-lobe overflow from a tidally distorted low mass star, may be initially magnetised. The magnetisation $1/\beta$ of the gas as it departs from the donor star can be estimated as

 \begin{displaymath}% {1\over \beta}\sim {B_\star^2 \over 8\pi \rho c_{\rm s}^2}, \end{displaymath} (1)

where $B_\star$ (measured in Gauss) is the stellar magnetic field at the surface, $c_{\rm s}$ is the speed of sound at the stellar surface, and $\rho$ is the density of gas as it becomes detached from the star. The density $\rho$ can be estimated as

 \begin{displaymath}% \rho\sim {\dot{M}\over 4\pi R_\star^2 c_{\rm s}}, \end{displaymath} (2)

where $\dot{M}$ is the average accretion rate of the companion, and $R_\star$is the radius of the star. We get

 \begin{displaymath}% {1\over \beta}\sim 0.3 \left({\dot{M}\over 10^{-9}~M_{\odot... ...right)^2 \left({R_\star\over 10^{11}~\hbox{cm}}\right)^2\cdot \end{displaymath} (3)

The streaming and shearing motion after the gas detaches from the star can further amplify the magnetic field. Clearly, there is a realistic parameter range where the initial magnetisation is high.
Motivated by these astrophysical considerations, in this paper we perform numerical experiments on gas discs which contain initially strong magnetic fields, with magnetic pressure comparable to that of the gas. Specifically, we initialise the azimuthal field so that it is subject to the Parker instability (PI, Parker 1966) in the vertical stratification.

In terms of initial conditions our physical set up is close to Machida et al. (2000). However, Machida et al. (2000) focused on the formation and evolution of a disc corona, and the spirit or their numerical experiment was different from ours. They have simulated the whole circular disc, and have introduced reflecting boundary conditions at the midplane of the disc, which suppresses field anchoring in the midplane. By contrast, our purpose is to understand the long-term dynamics of the disc, possible field confinement, and the dynamo processes related to strong magnetic fields. Therefore, our shearing-box simulations focus on a small part of the disc and study it with high numerical resolution. We simulate the fluid both below and above the midplane, and thus have no artificial boundary conditions at the midplane of the disc.

We find that both the short term and, more importantly, the long term behaviour of initially strongly magnetised discs is radically different from that of their weakly magnetised counterparts. We observe the following three-step dynamics: (a) the Parker instability expels azimuthal field in huge arcs, creating vertical field which becomes the seed for a strong magnetorotational instability; (b) matter sliding down inclined field lines stretches the azimuthal magnetic field and creates a vertically dependent large scale mean radial field; and (c) the Keplerian shear recreates azimuthal field from the stretching of the radial field. The latter step closes the dynamo cycle, in much the same way as was sketched by Tout & Pringle (1992) a decade and a half ago, although non-axisymmetric magnetorotational instability in the azimuthal field also plays an important role in creating accretion stresses in our simulations (Foglizzo & Tagger 1995; Balbus & Hawley 1992; Terquem & Papaloizou 1996). The azimuthal field remains strongly concentrated towards the disc midplane; this is in contrast with the simulations of Parker instability which do not include strong Keplerian shear in the radial direction (Kim et al. 1998). We show that the dynamo is robust and stable over at least tens of orbital periods, and that the accretion torque increases if we use a finer grid. We observe $\alpha $-values as high as 0.1 in our highest resolution run.

The plan of the paper is as follows. In Sect. 2 we detail the mechanics of our numerical experiments, and in the following section (Sect. 3) we test our computer code by comparing the results of Parker instability simulations without Keplerian shear to the extensive literature that exists on the PI under rigid rotation. In Sect. 4 we perform simulations with the Keplerian shear included, and present our main results. We devote the next section (Sect. 5) to analysing the confinement of azimuthal flux to the disc by a dynamo process. In Sect. 6 we conclude with the discussion of the astrophysical implications and possible future improvements of our work.

   
2 Description of the numerical experiment

In this section we describe the numerical method we use to solve the equations of ideal magnetohydrodynamics, the set up of an initial condition that is unstable to the Parker instability, and technical issues such as dissipation terms and boundary conditions. We use the Pencil Code[*] to solve the relevant partial differential equations.

2.1 Dynamical equations

We consider a local corotating patch of a Keplerian accretion disc. The coordinate axes are oriented such that x points radially away from the central gravity source, y points along the main orbital flow, while z points vertically out of the disc parallel to the Keplerian rotation vector  $\vec{\Omega}$. In the shearing sheet approximation the equation of motion, for the velocity field ${\vec u}$ relative to the Keplerian flow, is

 
$\displaystyle % \frac{\partial{\vec u}}{\partial t} + ({\vec u} \cdot {\vec \nabla}) {\vec u} + u_y^{(0)} \frac{\partial{\vec u}}{\partial y}$ = $\displaystyle 2 \varOmega u_y {\vec e}_x - (2-q) \varOmega u_x {\vec e}_y - \varOmega^2 z {\vec e}_z$  
    $\displaystyle + \frac{1}{\rho^{*}} {\vec J} \times {\vec B} - \frac{1}{\rho} {\vec \nabla}P + {\vec f}_\nu({\vec u}).$ (4)

Here $q=\partial\ln\varOmega/\partial\ln r$ is the shear parameter, with q=0 for rigid rotation and q=3/2 for Keplerian rotation. The left hand side includes advection both by the velocity field ${\vec u}$ itself and by the linearised Keplerian flow $u_y^{(0)}=-q \varOmega x$. The first two terms on the right hand side represent the Coriolis force in the x- and y-directions, modified in the y-component by the radial advection of the Keplerian flow, $\dot{u}_y=-u_x \partial u_y^{(0)}/\partial x$. The third term is the linearised vertical component of the gravity of the central object, while the Lorentz and pressure gradient forces appear in the usual way. The high order numerical scheme of the Pencil Code has very little numerical dissipation from time-stepping the advection term (Brandenburg 2003), so we add explicit viscosity through the term  ${\vec f}_\nu({\vec u})$, described in detail in Sect. 2.2.

2.1.1 Alfvén limiter

In this paper we consider stratified models with 4-16 orders of magnitude of range in densities. In low-density regions the Alfvén speed

 \begin{displaymath}% v_{\rm A} = \frac{B}{\sqrt{\mu_0 \rho}} \end{displaymath} (5)

may become very high when strong magnetic fields are transported away from regions near the midplane to form a low density corona, and we must reduce $v_{\rm A}$ artificially in order to avoid very low time-steps in the explicit numerical scheme of the Pencil Code[*]. We have done so by replacing the density $\rho$ in the Lorentz force term of Eq. (4) by a modified density $\rho^{*}$ defined as

 \begin{displaymath}% \frac{1}{\rho^*} = \frac{1}{\rho} \left[1+\left(\frac{v_{\rm A}^2}{M^2}\right)^n \right] ^{-1/n}\cdot \end{displaymath} (6)

Here M is a limiting value for the Alfvén speed and n is an index that smooths the transition from the regime where the Alfvén speed is not affected ( $v_{\rm A}<M$) to the region where the limiter applies ( $v_{\rm A}>M$). This yields an effective Alfvén speed $v_{\rm A}^* \approx v_{\rm A}$ in the limit $v_{\rm A}\ll M$, while in the limit of high Alfvén speeds the effective speed is

 \begin{displaymath}% v_{\rm A}^* = \frac{B}{\sqrt{\mu_0 \rho^*}} \rightarrow \fr... ...ac{M}{v_{\rm A}} = M \quad \textrm{for} \quad v_{\rm A} \gg M. \end{displaymath} (7)

We used a smoothing index of n=5. We are grateful to Heinemann for implementing this form of the Alfvén limiter in the Pencil Code (Heinemann et al. 2007). A similar Alfvén speed limiter was used by Miller & Stone (2000). The Alfvén limiter in the form presented here does not conserve momentum, but we have experimented with $M^2/c_{\rm s}^2=10$ and $M^2/c_{\rm s}^2=100$ and found no qualitative differences in our results. Thus we apply $M^2/c_{\rm s}^2=100$ to all 2D simulations, and $M^2/c_{\rm s}^2=10$ to the 3D simulations to get a longer time-step in those computationally expensive runs.

2.1.2 Induction equation

The magnetic vector potential ${\vec A}$ is evolved through the induction equation

 \begin{displaymath}% \frac{\partial{\vec A}}{\partial t} + u_y^{(0)} \frac{\part... ...vec B} + q \varOmega A_y {\vec e}_x + {\vec f}_\eta({\vec A}). \end{displaymath} (8)

The explicit resistivity term ${\vec f}_\eta({\vec A})$ is discussed in Sect. 2.2. Working with the vector potential ${\vec A}$ rather than the magnetic field ${\vec B}$ has the advantage that ${\vec B}={\vec \nabla}\times{\vec A}$ is kept solenoidal ( ${\vec \nabla}\cdot{\vec B}=0$) by mathematical construction. The stretching of the radial component of the magnetic field by the Keplerian shear is represented by the second term on the right hand side of Eq. (8) (Brandenburg et al. 1995). The current density ${\vec J}$, needed for the Lorentz force in Eq. (4), is calculated from Ampère's law ${\vec J}=\mu_0^{-1}{\vec \nabla}\times{\vec B}$. We set the vacuum permeability $\mu_0=1$.

2.1.3 Continuity equation and equation of state

The mass density is evolved in its logarithmic form

 \begin{displaymath}% \frac{\partial\ln \rho}{\partial t} + {\vec u}\cdot{\vec \n... ...tial y} = - {\vec \nabla}\cdot {\vec u} + f_{\rm D}(\ln \rho). \end{displaymath} (9)

The last term on the right hand side is an explicit diffusion term (see Sect. 2.2). It is advantageous to evolve the logarithmic mass density when a large dynamical range in density is considered. Also, in the initial magnetohydrostatic equilibrium, $\ln \rho$ varies parabolically with z, and the Pencil Code's finite difference approach to spatial derivatives yields a perfect derivative of parabolic curves (Brandenburg 2003).

We close the dynamical equation system by applying an isothermal equation of state with $P=c_{\rm s}^2 \rho$. The sound speed $c_{\rm s}$ is assumed constant.

  
2.2 Dissipation terms

The viscosity term ${\vec f}_\nu$ of Eq. (4) appears in full generality as

 
$\displaystyle % {\vec f}_\nu = \nu_3 \left[\nabla^6{\vec u}+({\vec S}^{(3)}\cdo... ...bla}\ln\rho)\right] + {\vec \nabla}\nu_{\rm shock}({\vec \nabla}\cdot{\vec u}).$     (10)

The viscosity includes both hyperviscosity with coefficient $\nu_3$ and shock viscosity with coefficient  $\nu_{\rm shock}$. A simplified third order rate-of-strain tensor  ${\vec S}^{(3)}$ is here defined as

\begin{displaymath}% S_{ij}^{(3)} = \frac{\partial^5 u_i}{\partial x_j^5}\cdot \end{displaymath} (11)

The shock viscosity coefficient is obtained by taking the negative part of the divergence of ${\vec u}$, then taking the maximum over three grid cells in each direction, and finally smoothing over three grid cells, to obtain

 \begin{displaymath}% \nu_{\rm shock}=c_{\rm shock} \langle {\rm max}[-{\vec \nabla}\cdot{\vec u}]_+ \rangle (\delta x)^2. \end{displaymath} (12)

Here $c_{\rm shock}$ is a parameter of order unity (Haugen et al. 2004).

For the resistivity term ${\vec f}_\eta$ in Eq. (8) we include both hyperresistivity and shock resistivity,

$\displaystyle % {\vec f}_\eta({\vec A}) = \eta_3 \nabla^6 {\vec A} + \eta_{\rm shock} \nabla^2{\vec A}+{\vec \nabla}\eta_{\rm shock}{\vec \nabla}\cdot{\vec A}.$     (13)

This form is the resistivity conserves all components of the magnetic field ${\vec B}$. We set $\eta_{\rm shock}=\nu_{\rm shock}$.

Mass diffusion consists of shock diffusion,

$\displaystyle % {\vec f}_D(\ln\rho) = D_{\rm shock} [\nabla^2\ln\rho+({\vec \nabla}\ln\rho)^2] + {\vec \nabla} D_{\rm shock}\cdot{\vec \nabla}\ln\rho,$     (14)

where $D_{\rm shock}=\nu_{\rm shock}$ is defined in Eq. (12). We also upwind the finite differencing of the advection term in the continuity equation (Dobler et al. 2006). Mass diffusion and upwinding are necessary to damp out spurious small scale modes that are left behind due to the dispersion error in the finite differencing of the advection term in the continuity equation.

We set shock diffusivity coefficients to $c_{\rm shock}=\eta_{\rm shock}=D_{\rm shock}=6.0$ in all simulations to dissipate energy in shocks forming far away from the midplane of the disc. In run S3D_256_Lz18, which is extended vertically compared to the other 3-D simulations, we had to double $c_{\rm shock}$ to dissipate enough energy in the regions above six scale heights from the midplane.

We test the dependence of our results on dissipation parameters by keeping the mesh Reynolds number approximately constant with increasing resolution, i.e. scaling hyperdissipation parameters by $(\delta x)^5$ and shock dissipation parameters by $(\delta x)^2$. Higher resolution simulations thus probe the effect of decreasing the dissipation coefficients.

Table 1: Simulation parameters.

   
2.3 Initial conditions

We assume that the initial magnetic field is toroidal, a reasonable assumption since Keplerian shear efficiently generates coherent toroidal field from a small radial component. Furthermore, we tune the initial field so that (a) magnetic pressure is the constant fraction  $\beta^{-1}$ of the gas pressure, and (b) vertical hydrostatic equilibrium is enforced, i.e.

 \begin{displaymath}% 0 = - \rho(z) \varOmega^2 z - {\partial\over \partial z}\le... ...rOmega^2 z - (1+\beta^{-1}) \frac{\partial P}{\partial z}\cdot \end{displaymath} (15)

Mathematically, the initial density, magnetic field, and vector potential are given by


 
  $\displaystyle \rho(z) = \rho_0 \exp[-z^2/(2 H_\beta^2)],$ (16)
  $\displaystyle \frac{B_y^2}{2\mu_0} = \beta^{-1} c_{\rm s}^2 \rho_0 \exp[-z^2/(2 H_\beta^2)],$ (17)
  $\displaystyle A_x(z) = \sqrt{2 \pi \mu_0 \beta^{-1} c_{\rm s}^2 \rho_0} H_\beta {\rm erf}[z/(2 H_\beta)].$ (18)

Here $H_\beta=\sqrt{1+\beta^{-1}} c_{\rm s}/\varOmega$ is the gas scale height, and $\rho_0$ is the initial density in the midplane. The thermal pressure alone would give rise to a scale height of $H=H_\infty=c_{\rm s}/\varOmega$. We will use this more familiar scale height as our unit of length, setting H=1throughout the paper.

Since we shall consider many (six to twelve) scale heights above and below the midplane, we run into the problem that the finite differencing of the magnetic pressure gradient underestimates its actual value and leads to spurious accelerations far away from the midplane. This is not a problem for the equilibrium density stratification, since the logarithmic density varies as a parabola with height over the midplane, and any parabolic shape has a perfect numerical derivative in the symmetric finite difference scheme of the Pencil Code. In order to give magnetic stratification equal opportunity as the density, and to obtain a perfect initial magnetohydrostatic equilibrium, we measure magnetic field and current density relative to their initial values, as explained below.

We split the field and the current into two components, as follows. For the Lorentz force in Eq. (4) and the induction term in Eq. (8) we replace


    $\displaystyle {\vec B} \rightarrow {\vec B}^{(0)} + B_y^{(1)}(z) {\vec e}_y,$ (19)
    $\displaystyle {\vec J} \rightarrow {\vec J}^{(0)} + J_x^{(1)}(z) {\vec e}_x.$ (20)

The initial magnetic field $B_y^{\rm (0)}(z)$ is given by Eq. (17), while the initial current density  $J_x^{\rm (0)}(z)$ is calculated as
                             Jx(0)(z) = $\displaystyle -\frac{1}{\mu_0}\frac{\partial B_y}{\partial z}$  
  = $\displaystyle \sqrt{2\mu_0^{-1}\beta^{-1} c_{\rm s}^2 \rho_0} [z/(2 H_\beta^2)] \exp[-z^2/(4 H_\beta^2)].$ (21)

The effect of this splitting is that the initial, mean field is not exposed to any (numerical or explicit) resistivity. Since we wish for resistive effects to go to zero with increasing resolution[*], the splitting does not introduce any extra spurious effects. We have checked that the Parker instability develops similarly in models with and without this field splitting, and have found that the models without splitting indeed converge towards the models with splitting as the resolution increases.


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{0385fi1a.eps}\hspace*{2mm}... ....eps}\hspace*{2mm} \includegraphics[width=8.4cm,clip]{0385fi1d.eps}\end{figure} Figure 1: Evolution of the Parker instability in 2D rigid rotation (run R2D_256). Overlaid on the density are magnetic field streamlines (white lines) and velocity field vectors (white arrows, averaged over 8 grid points, the longest arrows represent approximately four times the sound speed). The initial stratification is unstable to magnetic buoyancy, and magnetic field arcs begin to rise from the midplane. The arcs merge to form longer arcs, and eventually the system settles down into a new equilibrium state with two superarcs and four dense pockets of matter in the midplane.
Open with DEXTER

We perturb the initial Keplerian velocity field by Gaussian noise of amplitude $\delta u \sim 10^{-3}$. This perturbation seeds the linear Parker instability, and the noise also contains leading waves that can be amplified by non-axisymmetric magnetorotational instability in the azimuthal magnetic field in simulations with Keplerian shear.

2.4 Boundary conditions

We set the usual shearing sheet boundary conditions in the x and y directions (shear-periodic in x and periodic in y). At the upper and lower boundaries we impose free-slip conditions, $\partial u_x/\partial z=\partial u_y/\partial z=u_z=0$, for the velocity field ${\vec u}$ and perfect conductor, $A_x=A_y=\partial A_z/\partial z=0$, for the magnetic vector potential ${\vec A}$. This effectively allows the radial and azimuthal components of the magnetic field to evolve freely at the upper and lower boundary planes, while the vertical component of the magnetic field is forced to be zero. Globally the average value of all components of the magnetic field are then conserved. Perfect conductor boundary conditions do not allow azimuthal flux to flow out of the box, and therefore there is a concern that the presence of the boundaries will artificially enhance the field in the disc. To address this concern, we have made sure that our results have converged with increasing vertical extent of the box; see Sect. 5 and Figs. 7 and 12.

For the mass density $\ln \rho$ we set a symmetric boundary condition with $\partial\ln\rho/\partial z=0$ at the upper and lower boundaries. Although this precludes the gas pressure gradient from participating in magnetohydrostatic equilibrium there, the free-slip boundary condition for the velocity field, with zero vertical velocity component, means that this does not cause any unwanted accelerations.

2.5 Simulation parameters

We perform a series of simulations of the evolution of Parker and magnetorotational instabilities, usually starting from an equipartition ($\beta=1$) accretion disc in magnetohydrostatic equilibrium (as described in Sect. 2.3). Simulation parameters are written in Table 1. The main runs are the shearing sheet simulations S3D_*, representing varying resolution (S3D_256, S3D_512), higher initial ratio of gas to magnetic pressure (S3D_256_b3), increased vertical extent of the box (S3D_256_Lz18) and an exploration of the effect of a moderate amount of net vertical field (S3D_256_Bz0.03). For testing purposes we also run a series of 2D and 3D rigid rotation simulations (R2D_256, R2D_512, R2D_256_Lz24, R3D_256), which we present in the next section.

   
3 Testing the code: evolution of the Parker instability under rigid rotation

In this section we study the evolution of the ``pure'' Parker instability in 2D and 3D simulations without shear (Parker 1966; Kim et al. 1997). This allows us to test our code against the previously published results, and to review previous work done on the evolution of the Parker instability in systems with no systematic shear.

Basu et al. (1997) performed 2D inertial frame simulations of the non-linear evolution of the Parker instability. They found that primarily antisymmetric midplane crossing modes are excited by the random noise initial condition. The instability develops quickest in the tenuous gas high over the midplane. As the magnetic field lines rise in big arcs, matter streams unobstructed down the lines that are no longer horizontal. Dense pockets of gas (approximately a factor two times more dense than the average state) form in the midplane, anchoring the rising magnetic field lines there. The final state of the Parker instability, in which the tension in the bent magnetic field lines balances the buoyancy, has a lower energy than the initial equilibrium state (Mouschovias 1974) and is therefore preferred.

In Fig. 1 we show the evolution of the Parker instability obtained with the Pencil Code under rigid rotation in the 2D (y,z) plane (run R2D_256). Magnetic field lines initially rise in arcs of wavelength $\lambda_{\rm PI}\approx 7 H$. As matter streams sideways along the rising field lines (see arrows in Fig. 1) dense clumps of gas appear in the midplane after 10 orbits. Eventually the magnetic arcs merge and form two superarcs with four dense gas clumps in the midplane.

In Fig. 2 we show the maximum gas density in the box as a function of time. The clumps in the midplane are overdense by around factor 1.6-2.2 compared to the initial state, and this contrast is relatively well converged already at 256 $\times $ 128 grid points (as is the non-linear evolution of the Parker instability in general). Our simulations are in good qualitative agreement with the Basu et al. (1997) simulations.

3.1 Evolution of the azimuthal field

In Fig. 3 we show the mean azimuthal magnetic field as a function of height over the midplane in the saturated state of the 2D Parker instability under rigid rotation. The initial azimuthal field is redistributed by the Parker instability so that some of the midplane flux ends above four scale heights from the midplane. This scale is similar to the typical azimuthal wavelength of the Parker instability, so the final scale of vertical variation of the field likely comes about as the magnetic tension in the rising field lines grows to match the buoyant rise of the field lines. The non-linear transition from the initially stratified state to a state with magnetic arcs anchored in the midplane leads to the formation of an approximately force-free corona extending from $\sim$4 scale heights from the midplane. In the corona, the magnetic tension and pressure terms are large (compared to their difference) and oppositely directed, and current flows parallel to the magnetic field lines. Increasing the vertical extent of the box from 12 to 24 scale heights leads to a transport of some magnetic flux to regions with |z|>6 H, at the cost of a 20-30% decrease in the magnetic flux in the region |z|<6 H. The overall shape of the azimuthal magnetic field distribution is nevertheless relatively unchanged.


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0385fig2.eps} \end{figure} Figure 2: The maximum gas density as a function of time for 2D simulations of the Parker instability (under rigid rotation). Increasing resolution leads to no significant change of the maximum density, showing that the Parker instability operates in a robust way already at 256 $\times $ 128 grid points.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0385fig3.eps} \end{figure} Figure 3: Mean azimuthal magnetic field, averaged over the y-direction, as a function of height z over the midplane in 2D simulations with rigid rotation. The initial Gaussian magnetic pressure profile is lifted by the Parker instability, reducing the azimuthal flux by a factor of approximately 1.5 in the midplane. The flux is transported to above four scale heights from the midplane, and an approximately force free corona develops at these high altitudes. Increasing the vertical extent of the computational domain from 12 to 24 scale heights causes the flux to be transported slightly further away from the midplane, but the overall shape of the magnetic field distribution stays approximately the same.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0385fig4.eps} \end{figure} Figure 4: Mean azimuthal magnetic field versus height over the midplane for the Parker instability in 3D rigid rotation. The azimuthal field quickly spreads evenly over the entire vertical extent of the box, in contrast to the 2D case (shown in Fig. 3) which makes a non-linear transition to a state that still has azimuthal flux differences.
Open with DEXTER

3.2 Rigid rotation in 3D

Two-dimensional models of the Parker instability allow only the undular mode to develop, but not the related interchange mode (Matsumoto et al. 1993), and hence the field has no possibility to escape from the disc. In three-dimensional rigidly rotating disc simulations, the field does escape. In Fig. 4 we show the evolution of the azimuthal field in a 3D model including rotation (but still no shear). Now the azimuthal field quickly redistributes itself evenly with height over the midplane, leading to a state where vertical gravity is balanced only by thermal pressure. Allowing the mixed undular/interchange instability to take place leads to a genuinely turbulent state (Kim et al. 1998), which mixes the azimuthal field vertically and eventually spreads it uniformly throughout the entire extent of the box. In the next section we will show, however, that Keplerian shear completely suppresses the vertical spreading of the azimuthal field, which remains concentrated towards the disc midplane due to a dynamo coupling between Bx and By.


  \begin{figure} \par\mbox{\includegraphics[width=1.5cm]{0385fi5a.eps} \includegra... ...1.5cm]{0385fi5g.eps} \includegraphics[width=8.8cm]{0385fi5h.eps} } \end{figure} Figure 5: The vertical magnetic field at the sides of the simulation box for run S3D_512. The initial fastest growing wavelength of the Parker instability is around 5 scale heights in the azimuthal direction and 0.5 scale heights in the radial direction. The magnetorotational instability eventually sets in as the vertical field grows in amplitude, and the box turns turbulent with a strongly fluctuating vertical field.
Open with DEXTER

The Parker instability in a shearing, rotating frame (which is the relevant frame for most astrophysical purposes) must necessarily be more complicated than in the rigidly rotating frame. The intrinsically non-axisymmetric instability is deformed by the radial shear, while the whole coordinate frame can tap energy from the gravitational potential through Maxwell and Reynolds stresses. Already Shu (1974) expanded the Parker instability to include the effect of shear. The conclusion of the analytical stability analysis is that no amount of shear can completely stabilise against the Parker instability, as the instability can develop at such small radial scales that growth happens faster than shear.

   
4 The Parker instability in Keplerian shear

Having gained insight into the evolution of the Parker instability in 2D and 3D models with rotation but no shear, we now turn to the combined effect of Parker and magnetorotational instabilities in a shearing sheet model of a Keplerian accretion disc.

   
4.1 Creation of vertical field

In Fig. 5 we show the evolution of the vertical component of the magnetic field at the sides of the box for run S3D_512. The following sequence of events is observed: first the (sheared) Parker instability grows unaffected by the MRI, but as the vertical magnetic fields gets a significant amplitude, the MRI sets in and leads to a turbulent state where both PI and MRI determine the dynamics. The Parker instability initially sets in high up in the atmosphere where the buoyancy is strongest, at the typical wavelength $\lambda_y=5 H$ and $\lambda_x=0.5 H$. We actually observe that the radial wavelength of the Parker instability decreases with increasing resolution as the reduced dissipation allows growth at higher wave numbers. Eventually the whole disc goes into a non-linear turbulent state with strongly fluctuating vertical fields and twisted magnetic field loops extending far above the midplane of the disc. The typical radial coherence scale in this state is much longer than in the initial linear growth phase (compare $t=2.0 T_{\rm orb}$ to $t=15.0 T_{\rm orb}$ in Fig. 5).

The initial azimuthal field is also unstable to a (transient) non-axisymmetric magnetorotational ``instability'' (Keppens et al. 2002; Foglizzo & Tagger 1995; Balbus & Hawley 1992; Terquem & Papaloizou 1996). Magnetorotational instability in the azimuthal field has highest growth rate for modes with asmall vertical wavelength[*]. The absence of high kz modes in the first panel of Fig. 5 indicates that non-axisymmetric MRI is not the dominant driver of the linear dynamics of the system. However, in Fig. 6 we show the vertical velocity field at the y=-Ly/2 plane for run S3D_512 during the linear growth phase. The top panel exhibits the usual signs of Parker instability high up in the atmosphere, but changing the contrast by two orders of magnitude (bottom panel) reveals modes of very short radial and vertical wavelength in the midplane of the disc, representing trailing waves that have already been amplified by magnetorotational instability in the azimuthal field. Transient amplification of swinging waves may also be important in determining the non-linear evolution of the system at later times (e.g. Lesur & Ogilvie 2008; Rincon et al. 2007).


  \begin{figure} \par\includegraphics[width=7cm,clip]{0385fig6.eps} \end{figure} Figure 6: The vertical velocity component at the y=-Ly/2 plane during the linear phase of run S3D_512 (at $t=2.4~T_{\rm orb}$). The top panel shows the high radial wavenumber structures typical for the 3D Parker instability. However, changing the contrast by two orders of magnitude ( bottom panel) reveals structures of high vertical wavenumber in the midplane of the disc, likely a result of magnetorotational instability in the azimuthal field.
Open with DEXTER

Although production of vertical field is not a main ingredient of azimuthal MRI, any unstable mode with a non-axisymmetric vertical velocity component can create significant vertical field by stretching the large scale azimuthal field (although highly variable in the vertical direction). The production of large scale vertical field by the Parker instability, on the other hand, is demonstrated in Fig. 7. Here we plot the mean azimuthal field and the root-mean-square of the vertical field as a function of height over the midplane for the mixed PI/MRI shearing sheet simulations. Significant vertical magnetic fields of magnitude Bz $\sim$ 0.2 (with associated pressure $\beta_z$ $\sim$ 50) arise due to magnetic buoyancy. There is a reversal of the azimuthal magnetic field at around four scale heights from the midplane. The correspondence between run S3D_256 and run S3D_256_Lz18 is extremely close, so the sign reversal is not due to the presence of the vertical boundary. Increasing resolution leads to some increase in vertical field strength in the midplane, but the field reversal still takes place at the same height over the midplane. Although the non-linear state of the combined PI and MRI is time-dependent and fluctuating, there is a clear saturation and vertical confinement of the mean azimuthal magnetic field (we return to the issue of field confinement in Sect. 5).


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0385fig7.eps} \end{figure} Figure 7: Confinement of the azimuthal field in shearing box simulations. The mean azimuthal field and the root-mean-square of the vertical magnetic field are plotted as a function of height over the midplane. The Parker instability creates strong vertical fields from the initial azimuthal flux (dotted line in top plot). Vertical loss of azimuthal flux is nevertheless completely suppressed outside of a field reversal occurring at approximately four scale heights from the midplane. Doubling the resolution (S3D_512) leads to an increased vertical field strength and a decreased azimuthal field in the midplane, likely due to a significant increase in turbulent motions at higher resolution (see Table 2).
Open with DEXTER

4.2 Stresses

The magnetorotational instability in turn feeds off both the vertical field lines created by the Parker instability and the azimuthal field component still present in the non-linear state. The resulting Maxwell stress $\langle B_x B_y \rangle$, averaged over time and radial and azimuthal directions, is shown as a function of height over the midplane in Fig. 8. We have divided the stress into contributions from the mean field, $\langle \langle B_x \rangle_{xy} \langle B_y \rangle_{xy} \rangle_t$, and fluctuating field, $\langle (B_x-\langle B_x \rangle_{xy}) (B_y-\langle B_y \rangle_{xy}) \rangle_{xyt}$. The occurrence of a large scale Bx is discussed in detail in Sect. 5. The radial field couples with the large scale azimuthal field to yield a large scale structure in the Maxwell stress as well, relatively well converged when increasing the resolution.

The fluctuating magnetic field is associated with significant stresses around the midplane, with a Maxwell stress in the range $0.03\ldots0.1$ in the regions within a couple of scale heights from the midplane. The stress from this fluctuating field increases significantly with increasing resolution, as the decreasing dissipation on small radial scales allows the Parker instability to create stronger vertical fields, while the decreased numerical and artificial dissipation lets the both vertical field and azimuthal field MRI develop faster.


  \begin{figure} \par\includegraphics[width=8.1cm,clip]{0385fig8.eps} \end{figure} Figure 8: The Maxwell stress $\langle B_x B_y \rangle$ as a function of height over the midplane. The top panel shows the stress from the mean field $\langle B_x \rangle \langle B_y \rangle$, while the bottom panel shows the contribution from the fluctuating field ${\vec B}'={\vec B}-\langle{\vec B}\rangle$. The magnetorotational instability feeds off both the azimuthal fields and the large scale vertical fields created by the Parker instability, to create significant stresses, up to $\langle B_x B_y \rangle\approx-0.1$ in the midplane of the high resolution run S3D_512.
Open with DEXTER

In Fig. 9 we show the measured mean Maxwell stress $\langle B_x B_y \rangle$ as a function of time. The run S256_3D_Bz0.03 has a moderate vertical field (Bz=0.03) imposed through the box[*]. With this set up the magnetorotational instability sets in before the Parker instability does, creating significant stresses already after a few orbits. Eventually however the Parker instability develops as usual, and the stresses reach saturation at a level that is only slightly higher (in absolute value) than for the zero net vertical flux run. Thus it seems that it is not very important whether the magnetorotational instability develops before the Parker instability or vice versa. This situation may nevertheless change when going to either weaker azimuthal fields or stronger vertical fields, in which case the turbulent diffusion created by the magnetorotational instability may reduce the midplane azimuthal flux quicker than the Parker instability can grow (but see Blaes et al. 2007, where the Parker instability arises in simulations with a net azimuthal field that is much weaker than ours). In the case of a strong, imposed vertical field, there is however already an in inexhaustible source of accretion stresses present in the disc without the need to invoke an additional mechanism based on azimuthal fields and Parker instability.

Table 2: Statistical flow properties in the saturated, turbulent state.

In Table 2 we show the measured statistical properties of the non-linear state of the combined Parker and magnetorotational instabilities. We divide into the midplane regions with |z|<2 (containing 85% of the gas mass) and atmosphere regions with |z|>2 (containing 15% of the gas mass). The midplane regions are characterised by strong magnetic fields and high accretion through magnetic stresses, whereas the atmosphere has stronger velocity fields and higher Reynolds stress, but weaker magnetic fields and Maxwell stresses. Decreasing the initial magnetic pressure by a factor of three (run S3D_256_b3) leads to an expected decrease in both magnetic energies and stresses.


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0385fig9.eps} \end{figure} Figure 9: The mean Maxwell stress as a function of time. The high resolution run S3D_512 has faster growth of the stress and higher saturation level. Imposing the gas to a moderate net vertical field (dash-dotted line) lets the vertical field MRI develop from the beginning, but eventually the stresses saturate at an only slightly higher level (in its absolute value) than in the zero net vertical flux case.
Open with DEXTER

The z-dependent Reynolds and Maxwell stresses can be translated into an average turbulent viscosity (following Brandenburg et al. 1995),


  
    $\displaystyle \langle \rho u_x u_y \rangle = \frac{3}{2} \varOmega \nu_{\rm t}^{\rm (kin)} \langle \rho \rangle,$ (22)
    $\displaystyle -\frac{1}{\mu_0} \langle B_x B_y \rangle = \frac{3}{2} \Omega \nu_{\rm t}^{\rm (mag)} \langle \rho \rangle.$ (23)

Here $\nu_{\rm t}^{\rm (kin)}$ and $\nu_{\rm t}^{\rm (mag)}$ are the turbulent viscosities due to the velocity field and the magnetic field, respectively. We can furthernormalise the turbulent viscosities by the sound speed $c_{\rm s}$and gas scale height H (Pringle 1981; Shakura & Sunyaev 1973),


  
    $\displaystyle \nu_{\rm t}^{\rm (kin)} = \alpha_{\rm t}^{\rm (kin)} c_{\rm s} H,$ (24)
    $\displaystyle \nu_{\rm t}^{\rm (mag)} = \alpha_{\rm t}^{\rm (mag)} c_{\rm s} H.$ (25)

For the gas scale height we use the expression $H=H_\beta=\sqrt{1+\beta^{-1}}c_{\rm s}/\Omega$, with $\beta\equiv P/P_{\rm mag}$ given by the initial value. The resulting $\alpha $-values are given in Table 3. The large scale vertical fields arising from the Parker instability, together with a non-axisymmetric instability in the azimuthal field, give rise to $\alpha_{\rm t}\approx0.05\ldots0.1$, resulting in very high accretion rates through $\dot{M}=3\pi\nu_{\rm t}\varSigma_{\rm gas}$ (Pringle 1981).

Machida et al. (2000) observed similarly high $\alpha $-values in their global simulation of a strongly magnetised accretion disc, and they attributed the Maxwell stresses to a non-axisymmetric magnetorotational instability of the azimuthal field. Kim et al. (2002) considered the competition between Parker and gravitational instabilities in the context of a galactic potential. Simulations including rotation and shear nevertheless did not show any sign of the magnetorotational instability and the authors note that the magnetorotational instability in the azimuthal field grows too slowly to show up in the simulation time-scale of 4-5 orbits. In this work we find that the onset of turbulence fits a two layer scenario where the atmosphere is dominated by Parker instability and MRI in the ensuing vertical fields, whereas azimuthal field MRI drives the linear growth within a couple of scale heights of the disc midplane.

4.3 Energy spectra

The kinetic and magnetic energy spectra are shown in Fig. 10 for runs S3D_256 and S3D_512. We have taken the power at scale ${\vec k}$, $P({\vec k})$, and summed over shells of constant wave number $k=\vert{\vec k}\vert$, excluding the anisotropic scale $k=2\pi/(24 H)$ that is only present in the y-direction. Figure 10 shows that kinetic energy dominates over magnetic energy by approximately a factor of two at most scales, except for the 1-2 largest scales which are dominated by the magnetic field.

Table 3: Turbulent viscosity coefficients and $\alpha $-values, based on Eqs. (22)-(25).


  \begin{figure} \par\includegraphics[width=8.1cm,clip]{0385fi10.eps} \end{figure} Figure 10: Kinetic and magnetic energy spectrum for runs S3D_256 and S3D_512. The power  $P({\vec k})$ has been averaged over snapshots between 15 and 20 orbits and summed over shells of constant wave number $k=\vert{\vec k}\vert$. Kinetic energy dominates over magnetic energy at most scales, except for the few largest scale of the box.
Open with DEXTER

   
5 Field confinement


  \begin{figure} \par\includegraphics[width=8.1cm,clip]{0385fi11.eps} \end{figure} Figure 11: Different terms in the azimuthal component of the induction equation averaged over x and y and over evenly spaced snapshots between $t=10~T_{\rm orb}$ and $t=30~T_{\rm orb}$, for run S3D_256. The azimuthal magnetic field is in equilibrium between turbulent transport, representing compression and advection due to both the mean and the fluctuating velocity field (turbulent resistivity), and the Keplerian stretching term that creates By out of Bx.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{0385fi12.eps} \end{figure} Figure 12: The radial component of the magnetic field, Bx, averaged over the radial and azimuthal directions, for runs with Lz=12 H ( top panel) and Lz=18H ( bottom panel). The vertical structure quickly develops a peak in the midplane and one additional peak on each side of the midplane, and this structure stays statistically unchanged for the duration of the simulation. This radial field creates azimuthal field, by stretching in the Keplerian shear, that balances the turbulent resistivity acting on the mean azimuthal field (see Fig. 11).
Open with DEXTER

An important issue related to our proposed path to accretion is whether the original azimuthal flux can stay confined in the disc or whether it will escape to infinity by the action of turbulent resistivity (i.e. magnetic buoyancy). The combination of Parker and interchange instabilities in the non-shearing frame will eventually redistribute the azimuthal magnetic field evenly over the entire box height (see Fig. 4). We emphasise that we do not see the same behaviour in the shearing sheet. The magnetic field and velocity field stay statistically constant for the entire duration of the simulations, with no sign of decay or gradual loss of azimuthal flux (Fig. 7).

   
5.1 Vertical structure of the field

In Fig. 11 we dissect the equilibrium state of By by averaging different terms from the induction equation over x and y and over evenly spaced snapshots between $t=10~T_{\rm orb}$ and $t=30~T_{\rm orb}$. There is an almost perfect equilibrium between turbulent transport of magnetic fields ( $\dot{B}_y=-{\vec \nabla}\cdot({\vec u} B_y)$), representing compression and advection due to both the mean part of the velocity field and the fluctuating part (``turbulent resistivity'') and azimuthal field created by the Keplerian shear from the radial field [ $\dot{B}_y=-(3/2)\varOmega B_x$]. This Keplerian stretching allows for a closure of the entire process by which strong accretion occurs from initially purely azimuthal fields. The regular magnetic stretching term  $({\vec B}\cdot{\vec \nabla})u_y$(which does not include stretching by the Keplerian shear because we measure velocities relative to the main shear) is generally of the same sign as the turbulent transport term and thus acts oppositely of the Keplerian shear term.


  \begin{figure} \par\includegraphics[width=8.1cm,clip]{0385fi13.eps} \end{figure} Figure 13: The power contribution to the radially and azimuthally averaged radial field, as a function of time measured in orbits. Power is provided almost exclusively by the magnetic stretching term ${\vec B}\cdot{\vec \nabla}u_x$, while advection and compression are both sink terms.
Open with DEXTER

The vertical and temporal dependence of the radial magnetic field Bx is shown in Fig. 12 as a function of height over the midplane z and time t measured in orbits. After around 5 orbits a clear structure evolves out of the flow, with a peak of Bx in the midplane of the disc, followed by a valley and an additional peak at each side of the midplane. This structure stays statistically unchanged for the entire duration of the run. A similar steady state structure of the radial field was observed by Hanasz et al. (2002) in rigid rotation simulations of the galactic Parker instability. In the simulations of Hanasz et al. (2002) Bx would be either positive or negative in the midplane, depending on the radial wave number of the initial perturbation. In the shearing sheet the radial wave number of the PI is forced by the shear to be quite high, and our observation of a positive radial field in the midplane agrees with the high radial wave number simulations of Hanasz et al. (2002). Because we include shear in our simulations we additionally see the constant creation of a z-dependent By, by the Keplerian shear, which balances out the turbulent transport and prevents the azimuthal field from spreading evenly over the box height.

In Fig. 13 we show the power contribution to the large scale Bx as a function of time. We have first averaged Bx and the advection, compression and stretching terms of the radial component of the induction equation over the x- and y-directions. The power contribution of each term is extracted by multiplying the time derivative $\dot{B}_x$ with Bx and averaging over z. Power is provided almost exclusively by the magnetic stretching term $\dot{B}_x={\vec B}\cdot{\vec \nabla}u_x$, whereas both advection and compression extract energy from Bx at all times (except for a few peaks to around zero power). The magnetic stretching term provides a direct coupling between Bx and By, indicating that Bx is created from Byin a dynamo process.


  \begin{figure} \par\includegraphics[width=16.7cm,clip]{0385fi14.eps} \end{figure} Figure 14: Derivation of the dynamo $\alpha _{\rm dyn}$ from a single snapshot ( top panels) and from the average of several snapshots ( bottom panels). The (negative) integral of the magnetic stretching term $-\int \overline{{\vec B}\cdot{\vec \nabla}u_x}{\rm d}z$ is shown together with the azimuthal magnetic field  $\overline {B}_y$, both as functions of height over the midplane, in the left panels. The electromotoric force contribution from the stretching term is of opposite sign on each side of the midplane. The right panels show the correlation between $\overline {B}_y$ and $-\int \overline{{\vec B}\cdot{\vec \nabla}u_x}{\rm d}z$, with orange (grey) denoting points above the midplane and blue (dark) denoting points below. There is a very clear correlation (anticorrelation) above (below) the midplane, indicating a dynamo $\alpha $ that is positive above the midplane and negative below, and of order $\alpha_{\rm dyn}\approx0.023~c_{\rm s}$.
Open with DEXTER

5.2 Dynamo

The creation of a systematic radial field component can be understood from gas that plunges down field lines that have been bent by the Parker instability (Parker 1992; Hanasz et al. 2004; Rozyczka et al. 1995; Hanasz & Lesch 1993). As gas slides azimuthally along the field lines, the Coriolis force causes a counter clockwise rotation around dense clumps and clockwise rotation around underdense regions, twisting magnetic field lines in such a way that the perturbed electromotoric force points parallel to the (unperturbed) azimuthal field line. Radial field is subsequently created from the stretching of the perturbed field lines, in much the same way as imagined for the canonical $\alpha $-mechanism (Parker 1955; Moffatt 1978).

To explore this scenario in more quantitative terms we follow Moffatt (1978) and expand the velocity and magnetic fields in constant and fluctuating parts, ${\vec u}=\overline{{\vec U}}+{\vec u}'$, ${\vec B}=\overline{{\vec B}}+{\vec b}'$, leading to the following evolution equation for the mean magnetic field

 \begin{displaymath}% \dot{\overline{{\vec B}}} = {\vec \nabla}\times(\overline{... ... dyn}\overline{{\vec B}} - \eta_{\rm t} \overline{{\vec J}}). \end{displaymath} (26)

Here overlines denote azimuthal and radial (and possible temporal) averaging, while $\alpha _{\rm dyn}$ is a parameter that describes the proportionality between mean field  $\overline{{\vec B}}$ and fluctuating electromotoric force  $\overline{{\vec u}'\times{\vec b}'}$. The parameter  $\eta_{\rm t}$ represents turbulent resistivity, proportional to the current density $\overline{{\vec J}}={\vec \nabla}\times\overline{{\vec B}}$ of the mean field. Writing out the radial component of Eq. (26) we get

\begin{displaymath}% \dot{\overline{B}}_x = -\frac{\partial}{\partial z}(\overli... ...eta_{\rm t} \frac{\partial^2\overline{B}_x}{\partial z^2}\cdot \end{displaymath} (27)

The first term is due to advection of the large scale radial field by the large scale vertical velocity component. Figure 13 showed that the source of magnetic energy in the x-component is the magnetic stretching term. We may thus identify the positive contribution to $\dot{\overline{B}}_x$ with the properly averaged stretching term,

\begin{displaymath}% \dot{\overline{B}}_x = \overline{{\vec B}\cdot{\vec \nabla}u_x}. \end{displaymath} (28)

We can subsequently estimate $\alpha _{\rm dyn}$ from

\begin{displaymath}% \alpha_{\rm dyn} \overline{B}_y = -\int \overline{{\vec B}\cdot{\vec \nabla}u_x} {\rm d}z \end{displaymath} (29)

by integrating the magnetic stretching term over z. This is of course only a crude approximation of $\alpha _{\rm dyn}$ that ignores many of the complications of analysing the evolution of the mean field component (see e.g. Brandenburg et al. 2008; Brandenburg 2001), but this method gives a good order of magnitude estimate of the efficiency of creating radial field from the large scale azimuthal field.

In the left panels of Fig. 14 we plot the integral $-\int \overline{{\vec B}\cdot{\vec \nabla}u_x}{\rm d}z$ as a function of height over the midplane and compare it to the azimuthal magnetic field  $\overline {B}_y$. The electromotoric force contribution of the magnetic stretching term is of opposite sign at each side of the midplane. The correlation between $\overline {B}_y$ and $-\int \overline{{\vec B}\cdot{\vec \nabla}u_x}{\rm d}z$ is shown in the right panels of Fig. 14, with orange (grey) symbols indicating points above the midplane and blue (dark) symbols indicating points below the midplane. Anticorrelation below the midplane and correlation above indicates a positive  $\alpha _{\rm dyn}$ above the midplane and a negative  $\alpha _{\rm dyn}$ below the midplane, of the order $\alpha_{\rm dyn}\approx0.023~c_{\rm s}$ for run S3D_256. The higher resolution run S3D_512 gives $\alpha_{\rm dyn}\approx0.032~c_{\rm s}$. The fact that the helicity dynamo increases in efficiency when going to higher resolution, even though both the collisional hyper and shock resistivity coefficients are reduced, may indicate that the outlined dynamo is a fast dynamo, although future simulations applying resistivity on physical rather than on numerical grounds will be needed to corroborate this point (see e.g. Hanasz et al. 2002).

Considering the creation of radial field by the Parker instability in a galactic environment, Hanasz & Lesch (1993) predict $\alpha\simeq0.4\langle v_\phi \rangle$, where the ``cyclonic velocity'' $\langle v_\phi \rangle$ of Parker (1979) is approximately $0.1~c_{\rm s}$ in our simulations. The resulting dynamo coefficient $\alpha_{\rm dyn}\simeq0.04~c_{\rm s}$ is quite similar to what we find here from integrating the magnetic stretching term. In the theoretical framework of Moffatt (1978), on the other hand, the value of $\alpha _{\rm dyn}$ should be comparable to the correlation time of the turbulence $\tau_{\rm corr}$ times the mean helicity, $\alpha_{\rm dyn}\sim-(1/3)\tau_{\rm corr}\overline{({\vec \nabla}\times {\vec u}')\cdot{\vec u}'}$, at least in the limit of vanishing collisional resistivity. The helicity in our simulations is positive below the midplane and negative above the midplane, with an amplitude of $\vert({\vec \nabla}\times{\vec u}')\cdot{\vec u}'\vert\sim1$. Thus the inferred sign of $\alpha _{\rm dyn}$ fits well with the analytical theory, but the absolute value of $\alpha _{\rm dyn}$ is at least ten times smaller than the expectation based on a correlation time of order unity. This discrepancy may be simply due to the fact that kinematic dynamo theory is not applicable to our simulations, because the Lorentz force plays a significant role in determining the evolution of the velocity field and the magnetic field.

   
6 Summary and conclusions

In this paper we consider the evolution of strongly magnetised Keplerian accretion discs. Our numerical experiments show that the hydromagnetic state of the gas flow is very different from what is seen in zero net flux simulations. The Parker instability leads to huge magnetic arcs rising several scale heights from the disc midplane, and the magnetorotational instability in turn feeds off the vertical fields and creates a highly turbulent flow, an interaction that was predicted analytically by Tout & Pringle (1992). Although the flow is stochastic and time fluctuating, we have identified an underlying dynamo process that couples the vertically dependent mean radial and azimuthal magnetic field components. As gas slides down inclined field lines, it obtains a helical motion due to Coriolis forces, and thus the azimuthal field lines are twisted in such a way as to create a mean electromotoric force in the direction of the unperturbed field line - a configuration prone to create radial field. In turn the large scale radial field regenerates the azimuthal field through stretching in the Keplerian shear. Although Parker instability dominates the linear growth phase, we have found evidence for magnerotational instability in the azimuthal field as well. In the midplane of the disc, where the buoyancy is weak, azimuthal MRI drives the initial evolution towards turbulence (Foglizzo & Tagger 1995,1994; Terquem & Papaloizou 1996). These two related instabilities, magnetorotational instability in the vertical fields created by the Parker instability and magnetorotational swing instability in the azimuthal fields, both rely on azimuthal flux confinement and can coexist in the linear as well as in the non-linear state of transmagnetic accretion discs.

Such a path to accretion, based on the interaction of Parker and magnetorotational instabilities, has at least two appealing traits. First of all that the vertical fields that feed the magnetorotational instability are created in a transparent way by the Parker instability. Zero magnetic flux models must most likely rely on a small scale dynamo in order to create vertical fields, and there is mounting evidence that such a dynamo would not operate in the bulk part of accretion discs where the magnetic Prandtl number is much lower than unity (Schekochihin et al. 2005; Fromang et al. 2007). The second appealing result of our model is that the Maxwell and Reynolds stresses are significant ( $\alpha\approx0.1$). Such high accretion stresses could solve the problem that observed accretion rates are often one or two orders of magnitude higher than the accretion rates obtained in zero net flux MRI simulations (King et al. 2007).

The regeneration of azimuthal field by the shearing of an appropriate radial field was seen in all our simulations that included Keplerian shear. As magnetohydrostatic equilibrium is compromised by the Parker instability, gas streams down along inclined field lines. Coriolis force diverts the gas to the right, and a radial magnetic field is created as the azimuthal field is subjected to shear-regions typically the size of the Parker instability. Eventually magnetic reconnection leads to a coherent large scale radial magnetic field. This dynamo was predicted by Parker (1992) and subsequently observed in the rigid rotation simulations of Hanasz et al. (2002). To our knowledge we are the first to point out the relevance of Parker's fast galactic dynamo to accretion discs and how it closes the accretion loop by replenishing the azimuthal field that is lost by magnetic buoyancy.

The fine-tuned initial conditions with a purely azimuthal magnetic field and a constant ratio of magnetic to thermal pressure may be questioned. However our experiments with a combined azimuthal and vertical field shows that the Parker instability is robust even if the azimuthal field coexists with a moderately strong net vertical field, and that the additional vertical field component may indeed increase the accretion rate further.

Our results may also be relevant for star formation in the galactic centre. Although there is currently no coherent accretion disc structure, the population of young, massive stars in a disc-like structure close to the galactic centre points towards the brief existence of an accretion disc some million years ago (Levin & Beloborodov 2003; Nayakshin et al. 2007; Alexander et al. 2008; Milosavljevic & Loeb 2004). The disc was likely to be initially strongly magnetised, as indicated by the current high magnetisation of the circumnuclear molecular ring. Hence the type of MHD processes studied in this paper may be of central significance for the disc dynamics. The presented simulations of transmagnetic ( $\beta\sim1$) discs argues that discs dominated by magnetic pressure $\beta \ll 1$ are astrophysically viable. The existence of such discs was conjectured by Pariev et al. (2003) and their limitations and observational consequences were explored by Begelman & Pringle (2007). A number of problems in accretion disc theory are alleviated by the presence of super-equipartition magnetic fields, among them is the long-standing issue of self-gravity and fragmentation of AGN discs (Goodman 2003).

Future research into strongly magnetised accretion discs should also focus on the self-consistent modelling of the magnetisation of the material that feeds accretion discs in various environments and on the evolution of magnetised disc coronae, in light of effects such as reconnection, shearing of foot points, Ohmic heating and radiative cooling that take place there.

Acknowledgements
Computer simulations were performed at the PIA cluster of the Max Planck Institute for Astronomy and at the SARA Computing and Networking Service in Amsterdam. We would like to thank Micha\l Hanasz, Steve Balbus, and Andrew King for inspiring discussions. Axel Brandenburg, Jim Pringle and Richard Alexander are thanked for commenting an early version of the manuscript. We are grateful to Gordon Ogilvie and to the referee, Thierry Foglizzo, for pointing out the relevance of the azimuthal field MRI for our simulations.

References

 
Copyright ESO 2008