text:
# Numerical simulations of the Accretion-Ejection Instability in magnetised accretion disks

## 1 Introduction

Magnetic fields in accretion disks have been seen over the last decade to strongly effect the structural and dynamical properties of accretion disks as compared to the well defined standard steady-state, thin disk models (Shakura & Sunyaev ShakuraS73 (1973), Pringle Pringle81 (1981), Frank, King & Raine FrankKR92 (1992)). Since the early 1990s, numerical simulations of accretion disks have increasingly shown huge deviations from the smoothly varying, symmetrical (azimuthally and vertically), time-independent accretion disk models.

Although the magnetic field has been assumed to be of importance in generating instabilities within the disk (leading to turbulent transport), this has essentially been modeled as turbulent viscosity in the classical prescription $`\nu =\alpha c_sh`$ where $`\nu `$ is the viscosity, $`c_s`$ the sound speed, $`h`$ the disk height and $`\alpha `$ a dimensionless parameter. This then allows for the outward radial advection of angular momentum and inward radial transport of matter. However, this alone is insufficient to explain the time-variability of accretion disks such as relatively long timescale (hours to days) dwarf novae outburst of cataclysmic variables (e.g. Cannizzo & Mattei CannizzoM92 (1992)) or more rapid (of the order of seconds) QPOs (e.g. Swank et al SwankCMT97 (1997), Munoet al. Muno99 , Sobczak et al. Sobczaketal00 (2000)). The time variability of these systems has always been shown to be associated with the dynamics of the accretion disk itself - a feature that is not included within standard hydrodynamical treatments.

Since the discovery of the importance of the magneto-rotational (MRI) or Balbus-Hawley (BH) instability (Velikhov Velikhov59 (1959), Chandrasekhar Chandra60 (1960),Chandra61 (1961)) in the context of accretion disks by Balbus & Hawley (BalbusH91 (1991)), and its numerical demonstration by Hawley & Balbus (HawleyB91 (1991)), a number of local 3d simulations have been performed (e.g. Brandenburg et al. BNST95 (1995), Hawley et al. HawleyGB95 (1995), Stone et al. StoneHGB96 (1996), Hawley et al. HawleyGB96 (1996), Hawley 2000b ). These all show that turbulent transport is indeed generated through a local instability. The resulting transport, if it is measured in terms of the $`\alpha `$ prescription, is in the range expected from observations; however the resulting $`\alpha `$ parameter is strongly varying, meaning that the turbulence does not obey the simple prescriptions of the $`\alpha `$ model. Other features, such as time variability (Hawley & Krolik HawleyK00 (2000)), a vertically varying $`\alpha `$ (Brandenburg et al. 1996b ) and vertical structural asymmetries (Caunt Caunt98 (1998), Miller & Stone MillerS00 (2000)), imply that the magnetic field strongly affects the overall dynamics of an accretion disk (and in particular its vertical structure and its connection to the corona).

Jets, associated with accretion disks in Active Galactic Nuclei, X-ray binaries and Young Stellar Objects are certainly magnetised (Blandford & Payne BlandfordP82 (1982)) and the magnetic field at the disk surface is required to drive the matter across the slow magnetosonic point, after which centrifugal forces take over (Lovelace et al. LovelaceBC91 (1991)). The formation of jets is also observed to be connected to the accretion within the disk. In fact one the main attractions of MHD jet models is that they show that the jet is a very efficient way to extract angular momentum from the disk, in sharp contrast with $`\alpha `$-disk models where the angular momentum is evacuated outward radially by viscous stresses. A number of instabilities other than the MRI have been proposed for accretion disks (e.g. Parker Parker79 (1979), Papaloizou & Pringle PapP84 (1984), Spruit et al. SpruitSP95 (1995)); however the recent discovery of the accretion-ejection instability (AEI) (Tagger et al. TaggerHSP90 (1990), TP99) provides a possible link between radial and vertical angular momentum transport.

As shown by TP99, the instability occurs close to the inner edge of an accretion disk for which the quantity $`\mathrm{\Omega }\mathrm{\Sigma }/B_0^2`$ increases radially where $`\mathrm{\Omega }`$ is the angular velocity, $`\mathrm{\Sigma }`$ the surface density and $`B_0`$ the magnetic field for which the plasma beta (the ratio of gas to magnetic pressures) is of a moderate value (around unity). The instability appears as a low azimuthal wavenumber spiral wave, made unstable through the interaction with a Rossby vortex (associated with the gradient of vorticity) which it generates at its corotation radius. This also causes the emission of Alfvén waves from corotation vertically along magnetic field lines, as their footpoints are twisted. Hence, in the context of jet formation, the AEI naturally contains the highly desirable property of linking vertical propagation directly to the disk dynamics. Also, unlike the local MRI which is quenched for moderate fields, the AEI exists in situations for which the field pressure is comparable to the gas pressure, as may occur in the inner region of an accretion disk as magnetic field is advected inwards with matter and accretes towards the central object.

With the view of the low-frequency QPOs observed in X-ray binaries, in particular those hosting a black hole, the instability is also a good theoretical candidate (Varnière et al., VarniereRTD00 (2000); Rodriguez et al., RodriguezVTD00 (2000)). One can thus imagine a scenario in which, under initial low field strength situations the MRI initially acts, advecting the field radially inward by turbulent transport. If at least part of the vertical magnetic flux is advected with the gas the field in the inner disk region builds up, until a certain point is reached where the higher magnetic pressure quenches the MRI and the AEI sets in. This ‘magnetic flood scenario’ (Tagger Tagger99 (1999)) might explain the $``$ 30 minutes cycles of the micro-quasar GRS 1915+105 (Mirabel et al. Mirabel98 (1998)).

The inner region of accretion disks is increasingly becoming the focus of 3d numerical experiments. As computational power increases, global models (still limited by resolution and radial extent) are now being produced as well as the well-established local shearing-box models. Armitage (Armitage98 (1998)) describes a global model that omits vertical stratification and is limited in both vertical and radial extent, but illustrates effectively the global development of the MRI. More recently Hawley (2000a ) and Machida et al. (MachidaHM00 (2000)) describe simulations of accretion disks starting from toroidal configurations (i.e. a departure from Keplerian thin disks); at the time of preparation Hawley & Krolik (HawleyK00 (2000)) and Hawley (2000b ) have published numerical simulations of full accretion disk models.

We note here, and will discuss in our conclusions, that in these 3d simulations spiral waves are seen to occur close to the inner boundary of the disk as a result of magnetic stresses; this is very similar to the features expected from the AEI, and indeed very similar to the results presented here.

In a very different approach, Stehle & Spruit StehleS (2000) have presented a 2D ($`r,\varphi `$) model of a disk. The disk is considered as infinitely thin and embedded in vacuum. Their model is used to demonstrate the existence of the interchange instability, predicted to occur in strongly magnetised disks (Spruit & Taam SpruitT90 (1990), Spruit et al. SpruitSP95 (1995)), and associated with the radial gradient of the magnetic field. Unexpectedly, they find that the disk is also unstable for lower field strength than predicted.

Our model is very similar to theirs, with numerical characteristics optimised according to our previous knowledge of the AEI. We also consider an infinitely thin disk threaded by a moderate vertical magnetic field: this is justified by the properties of the AEI, which is essentially constant vertically across the disk. The disk is also embedded in vacuum but, in order to separate different physical effects, we consider in the present paper only configurations where initially the equilibrium magnetic field is due to external currents. As discussed below we use a logarithmic radial grid: this allows us to get a better resolution in the inner region of the grid, where it is necessary, and on the other hand to model a much larger radial extent of the disk and to get rid of unwanted effects of the boundary condition at the outer radius of the simulation. We believe that such effects may explain the stronger instability observed in their work, compared with ours. This numerical setup (with, as discussed below, a magnetic potential given by a Poisson equation) allows us to use well-established methods developed in the description of self-gravitating disks.

In spite of these differences, our results are quite similar to the ones they obtain at moderate field amplitude and we believe that, although the initial conditions differ and the numerical methods are independent, similar physics are at the heart of the two simulations.

This model allows for the development of the spiral wave and the vortex – at the heart of the instability. The paper is organised as follows: Section 2 describes the basic model including the fluid equations and implementation of the magnetic field, Section 3 provides details of the numerical methods used, in Section 4 we present our results and in Section 5 our conclusions.

## 2 The model

We present, in cylindrical coordinates, a global two-dimensional model of the inner region of an accretion disk around a central object. The model is 2D in the sense that we solve only for the horizontal ($`r,\varphi `$) components of the velocity, and that all perturbed quantities are assumed independent of $`z`$ within the disk. We assume that the disk is threaded by a poloidal magnetic field which is symmetric about the midplane, hence purely vertical at $`z=0`$. As described later, assuming that the disk is of a finite vertical height and embedded in vacuum, the magnetic field outside the disk (and in particular at its surfaces) can be derived from a magnetic potential, which in turn, is calculated from the distribution of magnetic field within the disk (this is analogous to calculating the gravitational potential in a self-gravitating disk). Currents within the disk can be derived from the jump in the horizontal component of the field across the disk. Hence magnetic stresses can be included in an otherwise purely hydrodynamic model.

The variables we solve for are the radial and perturbed azimuthal velocity, $`u_r`$ and $`u_\varphi `$, the surface density, $`\mathrm{\Sigma }`$ (which comes from the vertical integration of the density, $`\rho `$) and the vertical component of the magnetic field, $`B_z`$. The magnetic potential is solved separately as $`\mathrm{\Phi }_m`$. The “perturbed” component of the azimuthal velocity is defined as its deviation from keplerian rotation (so that it contains both the true perturbation and the equilibrium contribution of the pressure gradient). By substituting for the full azimuthal velocity $`U_\varphi =u_K+u_\varphi `$, where $`u_K`$ is the Keplerian velocity, the axisymmetric component of $`U_\varphi `$ can be enforced more accurately and is less prone to numerical roundoff errors (encountered from large azimuthal velocities and gravitational forces) which can lead to anomalous radial motion.

The instability can be successfully described using a 2d model because the vertical scale height of the perturbations is greater than the disk thickness (TP99), unlike the Magneto-Rotational Instability (MRI) for which the vertical wavelength is critical, so that 3d local models are required to study the full nature of the instability (e.g. Brandenburg et al. BNST95 (1995), Hawley et al. HawleyGB96 (1996), Stone et al. StoneHGB96 (1996)). Hence our simulations will not show the full turbulent transport of angular momentum expected within accretion disks, but will however allow us to isolate the effects of the AEI alone in their inner region, providing a first approach of the physics of the AEI, and valuable details for future 3d models.

### 2.1 Fluid equations

The model solves for the momentum, induction and continuity equations integrated vertically across the disk, in the hypothesis that (as expected for the AEI) the vertical scale height of the perturbations is larger than the disk thickness. The momentum equations for $`r`$ and $`\varphi `$ are given by

$`{\displaystyle \frac{u_r}{t}}=u_r{\displaystyle \frac{u_r}{r}}{\displaystyle \frac{U_\varphi }{r}}{\displaystyle \frac{u_r}{\varphi }}+{\displaystyle \frac{u_\varphi ^2}{r}}+{\displaystyle \frac{2u_\varphi u_K}{r}}`$
$`{\displaystyle \frac{1}{\mathrm{\Sigma }}}{\displaystyle \frac{P}{r}}+{\displaystyle \frac{B_z\mathrm{\Delta }B_r}{4\pi \mathrm{\Sigma }}},`$ (1)
$`{\displaystyle \frac{u_\varphi }{t}}=u_r{\displaystyle \frac{U_\varphi }{r}}{\displaystyle \frac{U_\varphi }{r}}{\displaystyle \frac{U_\varphi }{\varphi }}+{\displaystyle \frac{u_rU_\varphi }{r}}`$
$`{\displaystyle \frac{1}{\mathrm{\Sigma }r}}{\displaystyle \frac{P}{\varphi }}+{\displaystyle \frac{B_z\mathrm{\Delta }B_\varphi }{4\pi \mathrm{\Sigma }}},`$ (2)
where $`\mathrm{\Sigma }`$ is the surface density, $`U_\varphi =u_K+u_\varphi `$ and $`u_K=\sqrt{GM/r}`$ is the Keplerian velocity. Since pressure stresses are relatively unimportant, we assume for simplicity an isothermal equation of state so that the pressure gradient only changes as a result of advection of fluid density. The vertically integrated pressure, $`P`$, is given by

$$P=c_s^2\mathrm{\Sigma },$$
(3)
where $`c_s`$ is the sound speed. The terms $`\mathrm{\Delta }B_r`$ and $`\mathrm{\Delta }B_\varphi `$ represent the currents associated with the jumps in radial and azimuthal magnetic field across the surfaces of the disk. Assuming symmetry conditions apply then $`\mathrm{\Delta }B_r=B_r^+B_r^{}=2B_r^+`$ where the labels + and - note the values at the upper and lower disk surfaces, and similarly $`\mathrm{\Delta }B_\varphi =2B_\varphi ^+`$.

The continuity equation, integrated vertically, gives

$$\frac{\mathrm{\Sigma }}{t}=\frac{1}{r}\frac{}{r}(\mathrm{\Sigma }ru_r)\frac{1}{r^2}\frac{}{\varphi }(\mathrm{\Sigma }rU_\varphi ).$$
(4)
Similarly the vertical component of the induction equation takes the form:

$$\frac{B_z}{t}=\frac{1}{r}\frac{}{r}(B_zru_r)\frac{1}{r^2}\frac{}{\varphi }(B_zrU_\varphi ).$$
(5)
The identical form of the continuity and induction equations results from flux freezing in the ideal MHD limit, i.e. from the fact that both surface density and vertical magnetic flux are conserved in each tube threading the disk. It will be used below in the calculation of the magnetic potential.

### 2.2 Magnetic potential

Equations 4 and 5 show that the horizontal advection of the vertical magnetic field has an identical form to the advection of density in the limits of a vertically averaged model. Another identical relationship can be derived from assuming that the disk lies in vacuum. With this condition the absence of currents outside the disk allows one to describe the perturbed magnetic field outside the disk as the gradient of a magnetic potential, $`\mathrm{\Phi }_M`$. Following the derivation in TP99, we describe the field outside the disk as:

$$B=\mathrm{sign}(z)\mathrm{\Phi }_M,$$
(6)
so that both $`B_z`$ and $`\mathrm{\Phi }_M`$ are even in $`z`$, while $`B_r`$ and $`B_\varphi `$ are odd. In vacuum the divergence-free field condition gives

$$^2\mathrm{\Phi }_M=0.$$
(7)
On the other hand symmetry gives:

$$\frac{}{z}\mathrm{\Phi }_M(z=0^+)=\frac{}{z}\mathrm{\Phi }_M(z=0^{})=B_z^D,$$
(8)
where $`B_z^D`$ is the perturbed value of $`B_z`$ in the disk. Thus we can write throughout space:

$$^2\mathrm{\Phi }_M=2B_z^D\delta (z),$$
(9)
where $`\delta (z)`$ is the Dirac function. The form of this is identical to the standard Poisson equation for the gravitational potential of a self-gravitating disk:

$$\mathrm{\Delta }\mathrm{\Phi }=4\pi G\mathrm{\Sigma }\delta (z).$$
(10)
Hence, from the similarity between the Poisson equations on one hand, and on the other hand the similarity between the continuity equation (4) and the induction equation (5), we can solve for the magnetic potential using the classical Poisson kernel and more generally the methods devised in the simulation of self-gravitating disks (Binney & Tremaine BinneyT87 (1987)), as described in Section 3.3. The integration of the potential from individual field lines within the disk yields an expression for the total magnetic potential as:

$$\mathrm{\Phi }_M(r,\varphi )=\frac{1}{2\pi }_{r_1}^{r_2}_0^{2\pi }\frac{B_z^D(r^{},\theta )r^{}dr^{}d\theta }{\sqrt{r^2+r_{}^{}{}_{}{}^{2}2rr^{}\mathrm{cos}\theta }}.$$
(11)
This similarity explains that, as found by Tagger et al. (TaggerHSP90 (1990)), the waves we will describe are very similar to the spiral waves of self-gravitating disks. An important difference, however, is the negative sign in Equation 9, meaning that the magnetic field plays the role of a negative self-gravity: its presence “stiffens” the plasma against compression, whereas self-gravity helps it.

### 2.3 Initial Conditions

For the models presented here we assume, in the context of the accretion disk of a black hole in an X-ray binary, that the central object has a mass of $`M=10M_{}`$ and the disk extends from between $`r=1,000\mathrm{km}`$ to $`r=50,000\mathrm{km}`$ from the origin at $`r=0`$ with a constant aspect ratio of $`ϵ=h/r=0.1`$ where $`h`$ is the disk thickness. The only critical parameters, however, are $`ϵ`$ and the ratio $`\beta `$ of thermal to magnetic pressure, so that similar results would apply in disks of similar aspect ratio but very different dimensions, such as the accretion disks of young stellar objects or Active Galactic Nuclei.

The sound speed within the disk is calculated directly from the thin-disk approximation (e.g. Frank, King & Raine FrankKR92 (1992)) such that $`c_s=h\mathrm{\Omega }_K`$ where $`\mathrm{\Omega }_K`$ is the Keplerian angular velocity; hence $`c_s`$ varies radially as $`r^{1/2}`$. The disk starts from a steady-state velocity profile to which we add small scale random fluctuations with, unless otherwise stated, a maximum Mach number of 0.1 (hence the absolute magnitude of the initial fluctuations also decreases radially). The initial steady-state azimuthal velocity is set from the equilibrium condition, such that the centrifugal force counteracts radial forces acting on the fluid which include both gravity and gas pressure. Initially there are no currents within the disk and the magnetic field is of a purely external origin, being purely vertical with a radial gradient. This results in magnetic stresses which are initially vanishing, but build up with the magnetic potential during the simulation. The exact profiles of $`B_z`$, the initial external field, and $`\mathrm{\Sigma }`$, the surface density, are varied to examine the criterion of instability. Linear theory (TP99) lets us expect the disk to be unstable for initial profiles such that the quantity $`\mathrm{\Omega }\mathrm{\Sigma }/B_z^2`$ increases radially. This criterion is investigated in Section 4.3. The strength of the external field is also a variable parameter. The plasma beta, using the vertically integrated quantities, is defined as

$$\beta =\frac{8\pi P}{B^2}=\frac{8\pi \mathrm{\Sigma }c_s\mathrm{\Omega }_K}{B^2},$$
(12)
taking into account standard thin disk approximations. We vary the field strength using different values of $`\beta `$ defined at the inner edge of the disk. The instability occurs for moderately strong magnetic fields and is theoretically most active for $`\beta 1`$. We therefore perform tests on $`\beta `$ in the range $`0.1\beta 10`$.

### 2.4 Non-dimensionalisation

We choose typical length, time and density scales to non-dimensionalise the quantities. We take \[l\]=$`10^8\mathrm{cm}=1000\mathrm{km}`$, \[t\]=$`1000s`$ and $`[\rho ]=10^8\mathrm{g}\mathrm{cm}^3`$. These have the added advantage that velocities are returned as $`[v]=1\mathrm{km}\mathrm{s}^1`$. The disk dimensions are then $`r=\{1,50\}`$, the surface density, we take to be $`\mathrm{\Sigma }=10[\rho ][ł]=1\mathrm{kg}\mathrm{cm}^2`$ and the inner orbital period is then $`T_{r=1}=1.7\times 10^4[t]=0.17\mathrm{s}`$.

## 3 The numerical methods

Equations 1 to 5 are solved on a 2d cylindrical ($`r,\varphi `$) grid. A conservative scheme of the type described by Stone & Norman (StoneN92 (1992)) is used to advance the fluid properties in time. The scheme uses a staggered mesh such that $`\mathrm{\Sigma }`$ and $`B_z`$ are cell centered quantities and $`u_r`$ and $`u_\varphi `$ are edge centered. This is advantageous in the fact that derivatives, required for a particular variable, are typically calculated at the desired position and with increased accuracy. The grid itself is evenly spaced in azimuth and logarithmically in radius. This has two advantages for this particular model: the region of interest (close to the inner boundary) has greater resolution and waves propagating to the outer boundary are damped via inherent numerical dissipation, which guarantees the outer boundary condition for these simulations (radiation condition, i.e. no information coming in from the outer boundary). This is particularly important here in reference to the related Papaloizou-Pringle instability (Papaloizou and Pringle, 1984), which we should recover in the limit of vanishing magnetic field ($`\beta \mathrm{}`$). In this limit it is well known that a reflecting outer boundary can result in a much more violent instability (Narayan et al., 1987). We believe that this might explain the apparently more unstable numerical results of Stehle and Spruit (StehleS (2000)).

Averaging of quantities onto the staggered mesh takes into account the non-uniform nature of the grid. The resolution is typically taken to be 256x128 cells in $`r`$ and $`\varphi `$ although lower resolution tests have been made for comparison. The high radial resolution is used to capture the fine details of the instability. The azimuthal resolution is arguably excessive since the instability generates low azimuthal wavenumber spiral waves. However we found inaccuracies, or even numerical instabilities, associated with the regularisation of the logarithmic singularity of the Poisson kernel described by Binney & Tremaine (BinneyT87 (1987)) when the azimuthal and radial grid sizes are very different. Tests with modified smoothing schemes have been made with lower azimuthal resolution and indicate that the results presented here are not dependent upon the scheme, or higher resolution. Further details of the Poisson kernel are given in Section 3.3.

### 3.1 Evolution of the fluid

The equations are solved in a conservative form such that angular momentum, mass and total magnetic flux are conserved to numerical precision throughout the simulation. The method follows that of Stone & Norman (StoneN92 (1992)) such that the terms are split into source and transport terms: the source terms being centrifugal, gravitational, pressure and magnetic terms in the momentum equation and the transport terms being the advection of matter.

The source terms are calculated initially using standard second order finite difference operators on the staggered mesh (modified for a logarithmically spaced grid) followed by the transport terms. The transport terms are calculated using a standard second order Van Leer (vanLeer79 (1979)) upwind scheme.

We also invoke the FARGO scheme (Masset Masset00 (2000)) to allow for increased timesteps, typically limited by the underlying large velocity Keplerian motion. The Courant-Friedrichs-Lewy condition on the maximal timestep for numerical stability is thus reduced from $`dt\stackrel{<}{}dx/v_\theta `$ to $`dt\stackrel{<}{}dx/v_g`$, where $`dx`$ is the grid size and $`v_g`$ the group velocity of waves which can propagate in the simulation box. Typically $`v_g`$ will be of the order of the fast magnetosonic speed, i.e. comparable with the sound speed if $`\beta 1`$. The allowed timestep is thus increased by a factor $`v_\theta /v_gϵ^1`$, since $`v_\theta `$ is dominated by the keplerian velocity.

In fact $`v_g`$ can be estimated more precisely: from the dispersion relation given by TP99, the maximum group velocity can be approximated as

$$v_g\frac{v_A^2}{c_s},$$
(13)
where $`v_A=B_z^2ϵr/8\pi \mathrm{\Sigma }`$ is the Alfvén velocity. The magnetic contribution to the timestep is therefore of the order of

$$dt=\mathrm{min}\left(\mathrm{\Delta }x\frac{c_s}{v_A^2}\right),$$
(14)
taken over the whole of the numerical domain.

### 3.2 Boundary conditions

The effect of boundary conditions on the behavior and development of the AEI was discussed in detail in TP99.

#### 3.2.1 Inner boundary

Concerning the inner boundary, they showed that any reflecting condition (e.g. $`u_r=0`$, $`u_\theta =0`$, etc.) results in the existence of modes with the same instability criterion, but different frequency. A detailed discussion of the physics of the inner boundary of the disk is impossible: indeed, in the example of the accretion disks of black holes in X-ray binaries, we still do not know what exists between the disk and the black hole horizon; one could consider an ADAF, or the magnetic structure threading the black hole, associated with the Blandford-Zjanek (1977) mechanism. In both cases the detailed physics of the inner disk boundary is unknown. Hawley and Krolik (HawleyK00 (2000)) and Hawley(2000b ) have recently used a pseudo-newtonian gravitational potential to model the flow of the gas at the Last Stable Orbit near a black hole. This removes the need for an artificial boundary condition at the inner disk edge. However the low-frequency QPO which we consider as a manifestation of the AEI in X-ray binaries is observed mostly when the disk stays away from the Last Stable Orbit (Varnière et al., VarniereRTD00 (2000); Rodriguez et al., RodriguezVTD00 (2000)). Here for simplicity we adopt the commonly used condition $`u_r=0`$.

A disadvantage of this condition is that, since magnetic flux is advected with the gas, it piles up in the inner few grid-zones. However, we find that overall, the effective change in the field strength is acceptable for our current models, where we find that the plasma $`\beta `$ changes by 70% over the course of 200 orbits in the most extreme case. We feel that, although this is not the most desirable effect, it is acceptable for the purposes of this paper. Future models will allow for the magnetic flux to spread in the inner “hole”, between the disk and the central object. This will limit the build-up of the magnetic field (and its gradient) near the inner disk boundary.

#### 3.2.2 Outer boundary

At the outer boundary, reflecting conditions are often used for simplicity but are non-physical, unless one can argue for the existence of a sharp outer disk radius. For instance in the case of an X-ray binary, a sharp outer boundary may exist, associated with the physics of tidal effects and the gas flow from the companion. We have already mentioned that, as shown by Narayan et al. (1987), a reflecting outer boundary may make the instability much more violent. However in realistic conditions the ratio of the outer to inner disk radius is very large, and this effect should not act. Here the logarithmic radial grid allows us to extend the simulation grid to very large radius, and there the grid size also becomes large. On the other hand for modes of interest the dispersion relation (Tagger et al., TaggerHSP90 (1990)) shows that the radial wavelength becomes small at large radius. Numerically, when it becomes comparable with the grid size, the waves are heavily damped (this is a natural consequence of the use of low order derivatives and advection). This ensures that waves generated in the region of physical interest, at small or moderate radius, will be damped before they reach the outer boundary, and are thus not reflected. This effectively implements the correct boundary condition that no information should flow from infinity. A further advantage of the large radial interval is also that, for a double check, we can push the outer boundary far enough that waves cannot travel to the outer boundary within the limited duration of the simulations.

### 3.3 Calculation of the magnetic potential

As discussed in Section 2.2, the magnetic potential can be solved analogously to the gravitational potential in a self-gravitating disk. For speed of calculation we use the FFT method discussed by Binney & Tremaine (BinneyT87 (1987)). Rather than using the perturbed density, the perturbed vertical magnetic field is used as the source of the potential, and the Poisson kernel is modified to account for the change in coefficients (and sign) between Equations 9 and 10. As in Binney & Tremaine, we perform a change of variables to reflect the logarithmic grid and simplify the calculation. It can thus be shown that after transforming the magnetic field and potential such that

$`V_M(u,\varphi )`$ $`=`$ $`\mathrm{exp}(u/2)\mathrm{\Phi }_M[r(u),\varphi ]`$ (15)
$`b_z(u,\varphi )`$ $`=`$ $`\mathrm{exp}(3u/2)B_z^D(r(u),\varphi )`$ (16)
where $`V_M`$ is the reduced potential, $`b_z`$ the reduced perturbed magnetic field and $`u`$ the logarithm of the radius, $`r`$, Equation 11 can be rewritten as

$$V_M(u,\varphi )=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}_0^{2\pi }K(uu^{},\varphi \varphi ^{})b_z(u^{},\varphi ^{})𝑑u^{}𝑑\varphi ^{}$$
(17)
with the Poisson kernel defined as

$$K(uu^{},\varphi \varphi ^{})=\frac{1}{\sqrt{2[\mathrm{cosh}(uu^{})\mathrm{cos}(\varphi \varphi ^{})]}}.$$
(18)
This kernel has a logarithmic singularity where $`uu^{}=0`$ and $`\varphi \varphi ^{}=0`$. For this case Binney & Tremaine (BinneyT87 (1987)) give an approximate smoothing:

$$K(0,0)=\frac{1}{\pi }\left[\frac{1}{\mathrm{\Delta }\varphi }\mathrm{arcsinh}\left(\frac{\mathrm{\Delta }\varphi }{\mathrm{\Delta }u}\right)+\frac{1}{\mathrm{\Delta }u}\mathrm{arcsinh}\left(\frac{\mathrm{\Delta }u}{\mathrm{\Delta }\varphi }\right)\right]$$
(19)
where $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }\varphi `$ are mesh sizes. However this form is only valid for small and comparable values of $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }\varphi `$. Here the physics of the AEI lets us expect (as confirmed by the results below) only low-$`m`$ instabilities (where $`m`$ is the azimuthal wavenumber), and for speed of computation we often wish to use much fewer grid cells in $`\varphi `$ than in $`u`$. In that case the approximation (19) gives the unphysical property that $`K(0,0)<K(0,1)`$, which in turn generates numerical instabilities. In this case we use an alternate method, often used in self-gravitating disks, modeling the finite thickness of the disk (which does physically regularise the Poisson kernel). In order to do this we replace, in the computation of the kernel, the denominator $`|rr^{}|`$ by $`\sqrt{(}(rr^{})^2+h^2`$. This makes the kernel regular at $`r=r^{}`$; it removes the need for a separate evaluation of $`K(0,0)`$, which can be written as

$$K(\mathrm{\Delta }u,\mathrm{\Delta }\varphi )=\frac{1}{\sqrt{2[\mathrm{cosh}(\mathrm{\Delta }u)\mathrm{cos}(\mathrm{\Delta }\varphi )+ϵ^2\mathrm{exp}(\mathrm{\Delta }u)]}},$$
(20)
where $`\mathrm{\Delta }u=uu^{}`$, $`\mathrm{\Delta }\varphi =\varphi \varphi ^{}`$ and $`ϵ`$ is the aspect ratio of the disk. For the purposes of the magnetic potential alone, we have experimented with modifying $`ϵ`$ to effectively alter the smoothing length (however keeping the disk height identical) and find that it makes only a very slight difference to the extent of the disk that is affected strongly by the instability. The overall development of the instability is unchanged. This was to be expected since the smoothing affects only short-range interactions, while the instability has large wavelengths except at large radii where we wish it to be damped by the limited resolution of the grid. For the purposes of this paper, in showing the general characteristics of the instability, we generally use the Binney & Tremaine kernel; however we include later an example of a low azimuthal resolution simulation using the modified kernel (20).

The power of the FFT method of calculation is that the kernel (and its Fourier transform) need only be calculated once at the beginning of the simulation. At the beginning of each time step the Fourier transform of $`b_z`$ is calculated. The inverse Fourier transform of the convolution of this and the Fourier transform of $`K`$ is then performed to generate $`V_M`$ which in turn gives $`\mathrm{\Phi }_M`$.

The radial and azimuthal magnetic field at the disk surface can then be calculated from

$$𝐁=\mathrm{sign}(z)\mathrm{\Phi }_M$$
(21)
yielding $`B_r`$ and $`B_\varphi `$. Their jump across the disk (coming from the sign of $`z`$ in Equation 21) gives the current in the disk, and thus the Lorentz force in the source terms of Equations 1 and 2.

The initial setup of the field is somewhat unphysical since it assumes that the field is initially totally of external origin, and we have no currents in the disk. The current density in the azimuthal direction is given by

$$j_\varphi =\frac{B_r}{z}\frac{B_z}{r}=0,$$
(22)
therefore, for a magnetic field which is purely vertical at the disk, and decreases radially one has

$$\frac{B_r}{z}<0.$$
(23)
This implies that the field lines must bend inward towards the central object; in a physically realistic situation the opposite is usually considered, giving the familiar “hourglass” shape commonly used in disk-jet models. However, for the purposes of this paper, we prefer to use this setup which allows us to identify the effect associated with the field gradient, and leaves out the action of field curvature and equilibrium currents. Preliminary linear calculations have shown us that an “hourglass” configuration is slightly more unstable than the one we use here. The code is perfectly able to describe such initial configurations, with currents both in the central hole and within the disk. These will be considered in future publications.

## 4 Results

Our main goal is to characterize the instability, its dependence on the properties of the disk, and its influence on the global evolution and accretion rate. We will thus first vary different parameters affecting the instability, and its strength. We will then discuss in selected cases the detailed properties of the instability, and finally evaluate its efficiency for accretion.

### 4.1 Effect of the field strength

Let us first consider the influence of the magnetic field pressure in the disk. For this we choose equilibrium magnetic field and density profiles which, according to the linear theory, should give a rather strong growth rate: the instability criterion (TP99) is that

$$\frac{}{r}\mathrm{ln}\frac{\mathrm{\Omega }\mathrm{\Sigma }}{B_0^2}>0.$$
(24)
We thus take an equilibrium magnetic field varying as

$$B_0(r)=B_0(\mathrm{r}_{\mathrm{min}})\left(\frac{\mathrm{r}_{\mathrm{min}}}{r}\right)^{5/4},$$
(25)
which is the profile used by Stehle and Spruit (StehleS (2000)). This should be unstable if the surface density profile is flatter than $`\mathrm{\Sigma }r^1`$. For this first series of runs we choose a somewhat unrealistic flat density profile, and vary the parameter $`\beta `$. We know that the magneto-rotational instability exists only for $`\beta >1`$, whereas the AEI is most unstable for $`\beta 1`$.

Figure 1 shows the evolution of the rms velocity within the disk for values of $`\beta =0.1,0.5`$, $`1.0`$, $`5.0`$ and $`10.0`$. As is shown, the perturbed velocity grows rapidly for $`\beta =0.5`$ and $`\beta =1.0`$, exceeding the initial random noise after approximately 10 orbits. The growth is approximately exponential during this time. We will show in Sections 4.5 and 4.6 that the instability appears as a spiral mode, in agreement with the theory. For the other values of $`\beta `$ we see very little action in comparison. The $`\beta =0.1`$ and $`\beta =10.0`$ cases are virtually stable with the $`\beta =5.0`$ being somewhat more active but failing to achieve the same amplitude as the most unstable cases. This run is still interesting to observe as the growth of the instability is slower and the evolution of the spiral wave is somewhat cleaner, as presented later.

Also shown in the figure is the rms magnetic field. This remains relatively constant over time for all cases, with only moderate increases in $`\beta `$ occuring for the unstable cases (i.e. $`\beta =0.5`$, 1.0 and 5.0) as shown in Figure 2. The maximum change in $`\beta `$ occurs for the $`\beta =1.0`$ case. We see it decreasing to approximately $`\beta =0.6`$ after the initial development of the instability as the field and matter are accreted to the inner boundary. After this time we see some oscillatory features and finally a very gradual decrease to a value of $`\beta 0.3`$ at around 200 orbits. Hence the plasma beta decreases to approximately 30% its initial value for the most active case.

In the subsequent sections we examine more closely the most active values of plasma beta (namely $`\beta =0.5`$, $`\beta =1.0`$ and $`\beta =5.0`$). As seen from the growth and development of the perturbed velocity, the cases $`\beta =0.5`$ and $`\beta =1.0`$ are similar and this holds true for many other characteristics in the structural evolution. Similar features hold also for the $`\beta =5.0`$ case with development being somewhat slower.

### 4.2 Effect of the spatial resolution

As discussed in Section 3.3, the Poisson kernel given by Binney & Tremaine (BinneyT87 (1987)), and used in the computation of the magnetic potential, is only valid for high azimuthal resolution (comparable with the radial one). For this reason the simulations presented in the majority of this paper are for $`256\times 128`$ grid zones in $`u`$ and $`\varphi `$ respectively. The high radial resolution is needed to obtain a good description of the AEI, which depends critically on a resonance localised at the corotation of the modes (where their angular phase velocity equals the rotation velocity of the gas). On the other hand the instability is predicted to occur only for very low azimuthal wavenumber $`m`$, so that a high angular resolution would be necessary only to describe in details the dissipation of wave energy at small scales, as the non-linear evolution leads to wave steepening an possibly to the formation of a spiral shock. This goes far beyond the purpose of the present paper and would require much more elaborate simulations.

Thus simulations with low angular resolution, and increased radial resolution, would be enough for our present purposes. As discussed in Section 3.3, this leads to a numerical instability if we use the regularisation, given by Binney & Tremaine (BinneyT87 (1987)), of the logarithmic singularity of the Poisson kernel; we have made low azimuthal resolution runs with a different regularisation, corresponding to a smoothing by the finite vertical thickness of the disk. In this section we compare two simulations: one using the standard Binney & Tremaine kernel at $`256\times 128`$ and another with a smoothed kernel (assuming a disk aspect ratio $`ϵ=0.1`$) with a resolution of $`256\times 32`$. All other parameters are identical for the two runs, with $`\beta =1.0`$.

Figure 3 shows the different evolutions of the two cases. As seen, the initial evolution of the instability is virtually identical. A similar magnitude of $`u_{rms}`$ is achieved after approximated 20 orbits. Only the final, saturated state of the instability (which depends at least in part on the dissipation of wave energy at small scales) is slightly affected, both for the perturbed velocity and magnetic fields. This confirms that the smoothed kernel does give the correct physics we expect.

### 4.3 Effect of the radial density gradient

We now retain the same radial magnetic field profile, but change the density profile to check its effect on the growth of the instability. We have chosen not to try and reproduce the profiles used in the linear analysis of TP99: these profiles, imposed by the numerical method used in TP99 to solve the linear system, have a sharp and localised gradient at the corotation radius, in order to limit the range of density and/or magnetic field in the disk. They would thus be smeared rapidly by the non-linear evolution of the instability, precluding the more direct comparison of theory with numerical simulation one could expect. Work in progress with a different scheme for the solution of the linear problem should relax this condition, and allow a direct comparison of the theory with the non-linear simulations.

To study the stability criterion we chose two extra cases: one with a surface density varying as $`\mathrm{\Sigma }r^{3/2}`$ and a steeper one with $`\mathrm{\Sigma }r^2`$. These should be increasingly stable. The first case is identical to that chosen by Stehle & Spruit. According to the criterion (24), only the flat density profile should be subject to the AEI.

Figure 4 displays the growth of the instability for these three different gradients of surface density, along with a run of different resolution. The initial random velocities vary slightly since the deviation from Keplerian rotation (thus the effect of the pressure gradient on the equilibrium rotation curve) is included in the rms values.

This shows that, as expected, a radial surface density gradient stabilizes the modes. For $`\mathrm{\Sigma }r^2`$ the instability has totally disappeared. For $`\mathrm{\Sigma }r^{3/2}`$ it is still present, though weaker than with a flat density profile. Thus numerically the instability threshold seems to be steeper (in terms of the density profile) than given by the criterion (24). This can be attributed to the fact that, as discussed by TP99, there are actually two different contributions to the growth rate: the corotation resonance gives the criterion (24). However a second mechanism (actually the first discovered by Tagger et al., TaggerHSP90 (1990)) corresponds, as for galactic spirals, to the emission of an outgoing spiral wave, emitted beyond the corotation radius. This outgoing spiral is very visible in the plots, figures 5-7. This mechanism does not depend on the radial profiles (we refer to TP99 for a more complete discussion). Although we expect it to be weaker than the corotation resonance, it shifts the instability threshold leaving the $`\mathrm{\Sigma }r^{3/2}`$ case weakly unstable. The presence of the outgoing spiral will be discussed later in more details. In order to verify the instability of the $`\mathrm{\Sigma }r^{3/2}`$ profile, we have performed a very high ($`512\times 512`$) resolution run, giving results similar to the lower resolution one (though the growth rate is slightly lower, as shown in Table 1).

### 4.4 Growth rates

The plots of rms velocity (Figures 1, 3 and 4), show that in general the instability grows exponentially until it reaches a saturated state. We give in Table 1 the growth rates in terms of $`\mathrm{\Omega }_0`$ \- the angular velocity at the inner radius - of all the runs presented here covering the different resolutions, $`\beta `$ and radial density gradients. The growth rate indeed appears to be strongest for the flat density gradient with $`\beta =1.0`$, and to almost totally disappear when $`\beta `$ is one order of magnitude higher or lower. Similarly, as the density gradient gets steeper, the growth rate decreases and the case with $`\mathrm{\Sigma }r^2`$ is totally stable.

These values are comparable with the theoretically predicted ones, as given in TP99. As explained in the previous section, a detailed comparison is difficult, since TP99 use quite unrealistic radial profiles, in order to limit the difficulties associated with their solution method.

### 4.5 Evolution of the spiral wave

The evolution of the instability into a well developed spiral wave is best seen for the $`\beta =5.0`$ case. Figure 5 shows the evolution from the initial random velocity field at intervals of 20 orbits. The plots are of the radial velocity in the inner region of the disk.

By $`t=60`$ orbits (at the inner edge of the disk), an $`m=3`$ (3-armed) spiral is clearly evident. This then develops over time to lower wavenumbers, eventually becoming a tightly wound $`m=1`$ spiral after about 200 orbits. The disk goes through some intermediate stages for which a mixture of wavenumbers co-exist. This is discussed in more detail in the following section.

The development of the spiral from higher to lower wavenumbers is again in agreement with the theory: the most unstable azimuthal wavenumber decreases from $`m=`$ a few to $`m=1`$ as $`\beta `$ decreases. This decrease of $`\beta `$, shown in Figure 2, occurs as magnetic flux is advected inward with the gas to the inner region of the disk.

All the runs show similar behaviour to the $`\beta =5.0`$ case. The spiral always develops from high to small wavenumbers, eventually becoming $`m=1`$ in all cases.

### 4.6 Time analysis of the development of the instability

We further this analysis of the evolution of the instability by performing a spectral analysis of the development of its structure with time.

Figure 6 shows a contour plot of the radial velocity, measured on the axis $`\varphi =0`$, as a function of radius and time, for the $`\beta =1`$ run, with lighter regions indicating higher velocities. This provides some essential details on the nature of the instability. It becomes visible after $`t20`$ orbits. Higher velocities are seen to occur close to the inner radius of the disk. In the innermost region ($`r\stackrel{<}{}2`$), the contours shown are horizontal, whereas farther out they are oblique with a positive slope: the latter represent a wave travelling outward, while the former indicate a standing pattern. This is the structure expected for the AEI which (like self-gravity driven galactic spirals) distinguishes two regions in the disk, separated by the corotation radius (where the angular phase velocity of the wave equals the rotation frequency of the gas): inward from corotation, the instability establishes a standing pattern formed of an outgoing and an ingoing spiral. Beyond corotation, it emits a trailing spiral wave travelling outward.

This all occurs very close to the centre of the disk. Virtually no effects are observed for $`r`$ much greater than $`20`$. This is also a good indication that the outer boundary will have little effect on the results. Virtually no waves travel to the outer boundary and hence reflections and subsequent interactions can be neglected. Any significant signal returning from the outer boundary would appear in this plot as an oblique pattern with a negative slope.

Figure 7 shows in more details the outgoing character of the wave in the outer region. Figure 8 shows the vorticity in the flow. A sharp feature is seen at $`r1.5`$. This again is expected from the AEI, which grows by coupling the spiral pattern to a Rossby vortex it generates at its corotation radius. The sharp vorticity feature is a manifestation of this vortex. This will play a crucial role in future studies: the instability mechanism lies in the fact that the spiral inward from corotation extracts energy and angular momentum from the flow (thus causing accretion), and transfers them either to the Rossby vortex or to the outgoing spiral beyond corotation. The vortex is usually more efficient (as discussed in Section 4.3). It can be considered as a twisting of the footpoints of magnetic field lines in the disk. Thus, if we took into account a small but finite density above the disk, this twisting would propagate as Alfvén waves, carrying vertically toward the corona (where they could power a jet or an outflow) the energy and angular momentum extracted from the disk. This is the process which inspired the name of the Accretion-Ejection Instability.

In order to show the vortex more clearly, we have plotted in Figure 9 the flow pattern in the inner region of the disk. A turning point in the flow is seen close to the inner boundary in the top-right quadrant of the plot. This vortex is more clearly shown in the flow pattern of Figure 10. For clarity, both figures show the deviations from the underlying Keplerian flow. A flow pattern around a singular point is seen in both plots. This point is at the corotation radius.

Figure 11 shows details of the spatial structure of the waves during four intervals of time, for the $`\beta =5.0`$ run, which shows most clearly the spectral evolution. For each time interval it shows two plots: the left one is a contour plot of the radial velocity, analogous to the one shown in Figure 6. The right plot shows contours of the Fourier transform of the left one, giving as a function of radius, the frequencies at which perturbations evolve during the given time interval. Frequencies of $`\omega /\mathrm{\Omega }_02.2,1.5\mathrm{and}.7`$ appear succesively, corresponding to the $`m=3,2\mathrm{and}1`$ spirals seen in Figure 5. These frequencies correspond to corotation radii which increase ($`\omega /m`$ decreases) as $`m`$ decreases, as expected from the theory.

We note that results recently released by Hawley & Krolik (HawleyK00 (2000)) and Hawley (2000b ) show similar behaviour close to the inner edge of an accretion disk. Their numerical setup is much more complex than ours, and describes a 3D MHD model of the inner region of an accretion disk around a black hole with a pseudo-Newtonian potential. Although their goal is to study the magneto-rotational instability in a realistic 3D geometry, they note that after some time a significant part of the radial transport is done by waves. They show, in the inner region of the disk, plots very similar to our Figure 6. A detailed comparison is difficult, since this appears after an important evolution of the disk and magnetic field geometry. Furthermore, the relatively poor radial resolution imposed by the full 3D simulation would probably rule out a proper description of the AEI. However the similarity between their results and our Figure 6 (including the horizontal and oblique patterns discussed previously) is striking. We also note that, when this appears, the magnetic field has strongly evolved from their initial setup, as some vertical magnetic flux has been advected inward with the gas. In the inner region of the disk, they thus have an essentially vertical field of rather high amplitude (although the value of $`\beta `$ is not given). We thus consider that the wave patterns they obtain could very well be a manifestation of the AEI, rather than waves generated non-linearly by the magneto-rotational turbulence.

### 4.7 Accretion rates

From the simulations we are able to estimate the efficiency of the instability in accreting matter towards the central object. At a given radius the accretion rate is calculated as:

$$\dot{M}(r)=_0^{2\pi }ru_r(r,\theta )\mathrm{\Sigma }(r,\varphi )𝑑\varphi .$$
(26)
However, this value is highly variable. Plotting against $`r`$ at one point in time shows very little of the overall effect of the instability. We therefore show values that are averaged over space or time, in order to obtain a better understanding of the overall accretion within the disk. All figures shown are for the $`\beta =1.0`$ run.

Figure 12 shows the accretion rate against time, averaged between $`r=1`$ and $`r=10`$. This region is chosen since the instability has little effect at larger radii. This shows that even though the accretion rate is highly intermittent, it is predominantly positive. Three bursts of accretion are seen before the first 100 orbits, at $`t30`$, $`t50`$ and $`t75`$ orbits. These burst are also seen in the plots of rms velocity in Figure 1.

However, the physics of the instability lets us expect a more complex behavior: spiral density waves carry a negative flux of energy and angular momentum inside their corotation radius, and a positive one beyond corotation: this means that a wave inside corotation grows by extracting energy and angular momentum from the disk, and transfering them either to the Rossby vortex, at the corotation radius, or to the outgoing spiral, beyond it. Accordingly, we expect to see (at least during the initial phase when the instability grows following its linear behavior) accretion in the inner part of the disk, and outward motion (negative $`\dot{M}`$) beyond the corotation radius, which is typically of the order of $`1.52`$. This is exactly what we see in Figure 13a, which shows the accretion flux in the inner region of the disk, averaged between times 10 and 25 and for the same run with $`\beta =1`$. However the result is still quite noisy, especially at large radius: this turns out to come from our initialization, where we have introduced strong ($`.1`$) random fluctuations of the radial velocity, from which the instability grows. By the time these fluctuations have dissipated or traveled away, they have caused significant radial flows. In order to check this we have performed a new simulation with the initial rms noise reduced to $`.001`$: Figure 13b shows that the noise has strongly decreased at large radius, while the inner region is qualitatively similar (note that absolute values of $`\dot{M}`$ cannot be compared, since the instability now takes longer to grow from the initial conditions to a given amplitude). Finally we run a new case where the initial random perturbation now depends on $`r`$ as a Gaussian centered around $`r=5`$: Figure 13c shows the resultant accretion flow, which has the same structure as in Figure 13b, but where the fluctuations at large radius have totally disappeared. Figure 13d is a zoom on the inner disk region, showing that as expected $`\dot{M}`$ is positive inside a radius $`r1.7`$, and negative beyond, going to zero at large radius.

At larger times, however, accretion becomes positive beyond the corotation radius: this is shown by Figure 14, where $`\dot{M}`$ is averaged over the whole simulation. It is interesting to see that, in the region between corotation and $`r=10`$, most of this positive accretion occurs by bursts. These (already noted in Figures 1 and 12) are shown on Figure 15, where $`\dot{M}`$ is shown as a function of $`r`$ and $`t`$. The bursts are seen as strong streaks propagating outward from the corotation region. We conclude that this is due to non-linear effects (e.g. shocks, or global re-arrangements of the current profile) not included in the linear theory, but responsible for the saturation of the instability after a finite time. We defer their full characterization to future work where more realistic magnetic field and current profiles will be used.

It is customary to measure the efficiency of accretion in terms of a turbulent viscosity, parametrised with the parameter $`\alpha _{ss}`$ of Shakura & Sunyaev (ShakuraS73 (1973)). This allows a comparison of the instability mechanism with standard $`\alpha `$-disk model. We thus convert the accretion rate to an equivalent value of $`\alpha _{ss}`$, although we emphasize that the physics of the AEI (which forms large-scale, quasi-stationary patterns where the accretion energy is carred away by waves) is very different from the assumptions underlying the Shakura-Sunyaev model, i.e. small-scale turbulence leading to local deposition of the accretion energy. Here, as explained in previous sections and in more details in TP99, the instability mechanism relies essentially on the excitation of the Rossby vortex by the spiral wave. This occurs through the long-range action of the magnetic potential $`\mathrm{\Phi }_M`$. Thus the accretion energy is not deposited locally, but transported at large distances by the magnetic stresses figuring in equations (1-2). The present estimate can thus be taken only as a convenient measure of the efficiency of the instability to cause turbulent transport. Furthermore, a reliable assesment of the accretion rate due to the AEI would require taking into account the vertical emission of Alfvén waves (i.e. include the effect of a finite density corona above the disk), and effects associated with the finite thickness of the disk.

From standard thin disk approximations of steady accretion disks, one has

$`\nu \mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \frac{\dot{M}}{3\pi }}\left[1\left({\displaystyle \frac{r_{}}{r}}\right)^{1/2}\right],`$ (27)
$``$ $`{\displaystyle \frac{\dot{M}}{3\pi }},`$
(Frank, King & Raine FrankKR92 (1992)) where $`\nu `$ is the viscosity, $`\mathrm{\Sigma }`$ the surface density and $`r_{}`$ is the radius of the central object. The standard prescription for the viscosity is given by

$$\nu =\alpha _{ss}c_sh.$$
(28)
Combining Equations 27 and 28 along with another thin disk assumption, that $`c_s=h\mathrm{\Omega }`$, yields

$$\alpha _{ss}(r,t)=\frac{\dot{M}(r,t)}{3\pi ϵu_\varphi (r,t)r\mathrm{\Sigma }(r,t)},$$
(29)
where $`ϵ`$ is the disk aspect ratio. The values of $`u_\varphi (r,t)`$ and $`\mathrm{\Sigma }(r,t)`$ represent azimuthally averaged quantities resulting in an $`\alpha _{ss}`$ that depends on radius and time.

Figure 16 shows the effective $`\alpha _{ss}`$ against time when averaged between $`r=1`$ and $`r=10`$, for the $`\beta =1`$ run with Gaussian initial perturbations. It is negative while the wave is in its linear growth phase, but becomes positive in the non-linear stage. It is again strongly intermittent, showing the same bursts as in previous figures.

It is useful to compare values of $`\alpha _{ss}`$ with values derived from numerical models of turbulent accretion disks, for which the use of the $`\alpha `$ parameter is more applicable. The use of shearing-box models (e.g. Brandenburg et al. BNST95 (1995), Brandenburg et al. 1996a , Stone et al. StoneHGB96 (1996), Hawley et al. HawleyGB96 (1996)) over recent years have given insights to the turbulent transport of energy via calculations of Reynolds and Maxwell stresses as well as other methods (such as from the overall increase in thermal energy). Various values of $`\alpha `$ calculated from different quantities by Brandenburg et al. (1996a ) indicate that their $`\alpha _{ss}`$ is typically of the order of 0.001 to 0.01 with maxima generally around 0.01 to 0.02. These are comparable to those of Stone et al. (StoneHGB96 (1996)) and and Hawley et al. (HawleyGB96 (1996)).

With the AEI simulation in 2d for $`\beta =1.0`$ we see similar values of the viscosity parameter, $`\alpha _{ss}`$. The averages over time and space yield values typically between 0.001 and 0.006. Maximum values at a given position and time are around 0.06. These are still somewhat smaller than the recent values given by Hawley & Krolik (HawleyK00 (2000)) for their 3d model of the inner accretion disk and also that of Armitage (Armitage98 (1998)).

The strength of the instability, depending on $`\beta `$ and other parameters, affects the value of $`\alpha _{ss}`$. Table 2 shows the maximum values of $`\alpha _{ss}`$ from each of the runs.

Again we see that the instability is strongest for moderate values of $`\beta `$, i.e. around $`\beta 0.5`$. Indeed, we see that values of $`\alpha `$ are strongly peaked for $`\alpha =0.5`$ and are more comparable to values desirable for models of cataclysmic variables (Cannizzo Cannizzo93 (1993)). Recent global simulations of the MRI (Armitage Armitage98 (1998), Hawley 2000a ) indicate that values of $`\alpha `$ derived from global numerical simulations are generally higher than local shearing-box simulations (e.g. Brandenburg et al. 1996a , Stone et al. StoneHGB96 (1996), Miller & Stone MillerS00 (2000)) and hence the possibility exists of cataclysmic variables occuring as a result of a combination of both the AEI and the MRI. We must also emphasize that the physics discussed here concerns only the innermost part of the disk. Presumably, if accretion proceeds faster at larger radii, the resulting inflow of matter and magnetic flux will strongly affect the evolution of the AEI and its consequences on inward transport.

## 5 Conclusions

In this paper we have presented a two-dimensional model of the inner region of an accretion disk threaded by an initially vertical magnetic field, surrounded by vacuum. Starting from an equilibrium situation with additional random small scale fluctuations, an instability is seen to occur on a timescale of a few tens of orbits. The instability appears as a low azimuthal wavenumber spiral wave, amplified mainly by (and feeding) a Rossby vortex. These results are comparable (with a slightly different physical setup and a more adapted numerical one) with the ones obtained recently by Stehle and Spruit (StehleS (2000)).

The instability is seen to occur most strongly for a plasma beta of $`\beta 1`$ (i.e. equipartition between the gas thermal pressure and magnetic pressure in the disk), with its amplitude and growth rate rapidly decreasing for both stronger and weaker fields. Also, the growth rates (and the existence of the instability) depend on the radial gradients of equilibrium quantities.

As the spiral wave develops and affects the disk structure (in particular decreasing its $`\beta `$ value) it is seen to progress from larger to smaller azimuthal wavenumbers until it reaches $`m=1`$, i.e. a one-armed spiral.

These properties (timescales, development, characteristics, stability criterion) allow us to identify the instability as the Accretion-Ejection Instability, and globally to confirm qualitatively and quantitatively the behaviours predicted theoretically by TP99.

The instability is seen to generate accretion in the inner region of the disk. The amplitude of the accretion – calculated through an equivalent value of the standard $`\alpha `$ parameter of viscous transport (Frank, King & Raine FrankKR92 (1992)) – is comparable, when compared to other numerical simulations of accretion disks (e.g. Brandenburg et al. BNST95 (1995), Stone et al. StoneHGB96 (1996)), with the one generated at lower magnetic pressure by the magneto-rotational instability (MRI) of Balbus & Hawley (BalbusH91 (1991)) .

Furthermore, we have noted that waves observed in more recent publications of global, 3D simulations of the inner region of accretion disks (Hawley & Krolik HawleyK00 (2000), Hawley 2000b ) show properties (development of spiral waves with a characteristic propagation pattern) which are very similar to the ones described here. These waves, which seem to cause a significant part of the accretion, develop once the global magnetic structure has evolved in the inner region of the disk to become quite similar to the one used here. It is thus quite possible that the AEI is also present in these simulations, once the Magneto-Rotational Instability has caused the disk to evolve to a favorable magnetic configuration (in terms of geometry and strength) in its inner region.

Future work needs to extend the current simulations to 3D, in two steps: the first one will be to consider a disk which for simplicity is still infinitely thin, but embedded in a corona of small but finite density: this will allow us to describe the vertical ejection, as Alfvén waves, of the energy and angular momentum extracted from the disk and presently stored in a Rossby vortex. It will thus be possible to consider how they can ultimately be deposited in the corona, where they could energize a jet or an outflow.

A second step will be to turn to fully 3D simulations. This leads to numerical constraints which are quite different from the ones encountered in the simulation of magneto-rotational turbulence, since the physics of the AEI is quite different (e.g. it needs very few grid points within the disk in the vertical direction, but a large number in the radial one). This will allow us to address the acceleration of the gas, across a slow magnetosonic point, to finally form a jet.

###### Acknowledgements.

The authors gratefully acknowledge the help of F. Masset in providing valuable support in the development of many of the numerical methods described in this paper. They also acknowledge many discussions with R. Stehle and H. Spruit on their comparable simulations.