EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 516 (June-July 2010) A&A, 516 (2010) A31 Full HTML
Free Access
Issue
A&A
Volume 516, June-July 2010
Article Number A31
Number of page(s) 9
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201014000
Published online 23 June 2010
A&A 516, A31 (2010)

Rossby wave instability and three-dimensional vortices in accretion disks

H. Meheut1 - F. Casse1 - P. Varniere1 - M. Tagger2,1

1 - AstroParticule et Cosmologie (APC), Université Paris Diderot, 10 rue A. Domon et L. Duquet, 75205 Paris Cedex 13, France
2 - Laboratoire de Physique et Chimie de l'Environnement et de l'Espace (Université d'Orléans/CNRS), Orléans, France

Received 6 January 2010 / Accepted 31 March 2010

Abstract
Context. The formation of vortices in accretion disks is of high interest in various astrophysical contexts, in particular for planet formation or in the disks of compact objects. But despite numerous attempts it has thus far not been possible to produce strong vortices in fully three-dimensional simulations of disks.
Aims. The aim of this paper is to present the first 3D simulation of a strong vortex, established across the vertically stratified structure of a disk by the Rossby wave instability.
Methods. Using the ersatile advection code (VAC), we set up a fully 3D cylindrical stratified disk potentially prone to the Rossby wave instability.
Results. The simulation confirms the basic expectations obtained from previous 2D analytic and numerical works. The simulation exhibits a strong vortex that grows rapidly and saturates at a finite amplitude. On the other hand the third dimension shows unexpected additional behaviors that could be of strong importance in the astrophysical roles that these vortices can play.

Key words: Accretion, accretion disks - protoplanetary disks - hydrodynamics - instabilities - methods: numerical

1 Introduction

Vortices have been actively searched for in accretion disk theory and numerical simulations because of their multiple astrophysical interests. In particular for protoplanetary disks, vortices have been invoked (Barge & Sommeria 1995) to shorten time scale for the growth of grains to planetesimals to a level that be feasible. The streaming of the gas inside a vortex could speed up this growth by concentrating dust grains in its center, by what is sometimes called the ``tea leaf'' effect, whereby tea leaves accumulate in the center of a cup when it is stirred. This mechanism has been studied numerically (Lyra et al. 2008,2009; Johansen et al. 2004) but only for the 2D case of an infinitely thin disk. We note that an alternative mechanism for trapping the dust grains based on high pressure region was proposed by Johansen et al. (2007).

On the other hand, despite multiple attempts to generate realistic vortices in fully cylindrical 3D numerical simulations of disks, even the most advanced ones (Barranco & Marcus 2005) have only found off-midplane vortices. This may be explained (Tagger 2001) by the fact that due to both differential rotation and differential vorticity (which is absent in the shearing box model used in many simulations) vortices are very rapidly sheared away and could survive only by being continuously re-generated.

We present the first 3D numerical simulations of the Rossby wave instability (RWI), which is precisely a mechanism to do this: in certain conditions this instability generates Rossby vortices that grow exponentially with time. The existence of this instability in protoplanetary disks was first mentioned in Lovelace et al. (1999), and it has been suggested later (Varnière & Tagger 2006) that at the edges of the ``dead zone'' of protoplanetary disks, where the gas is not ionized and thus should not be turbulent, the conditions are such that the RWI can grow; it should thus both help accretion to proceed across the dead zone and help planets to grow there. The simulations we present confirm the basic properties of the instability, but they show new features in its 3D structure, whose effects on accretion and the growth of planetesimals will need to be understood.

The paper is organized as follows: after recalling the basic properties of the RWI we will present the physical setup, then the numerical one. Section 4 will present the results, and Sect. 5 a discussion and conclusion.

2 The Rossby wave instability

The RWI has been found and discussed in different contexts of differentially rotating disks. A discussion of this history has been given by Varnière & Tagger (2006), and we will not repeat it here. It can exist in galactic disks (Sellwood & Kahn 1991; Lovelace & Hohlfeld 1978) as well as in accretion disks (Lovelace et al. 1999; Papaloizou & Pringle 1985; see also Li et al. 2001,2000). It can be seen as the form that the Kelvin-Helmoltz instability takes in differentially rotating disks, and it has a similar instability criterion: an extremum in the quantity $\mathcal{L}$ related to the vorticity of the equilibrium flow. In an unmagnetized thin disk this quantity can be written as (Li et al. 2000)

\begin{displaymath}\mathcal{L}=\frac{\Sigma \Omega}{2\kappa ^2}\frac{p}{\Sigma ^... ...gma}{2(\vec \nabla \times \vec v)_z}\frac{p}{\Sigma ^ \gamma}, \end{displaymath} (1)

where $\Sigma$ is the surface density, $\vec v$ is the velocity of the fluid, $\gamma$ the adiabatic index, $\Omega$ the rotation frequency and $\kappa^2=4\Omega+2\Omega \Omega '$ the epicyclic frequency squared (so that $\kappa^2/2\Omega$ is the vorticity). Here the prime notes a radial derivative.

\begin{figure} \par\includegraphics[width=8.8cm,clip]{14000f1.eps} \end{figure} Figure 1:

Schematic view of the Rossby wave instability. Rossby waves are generated in the region of the density extremum, and density waves are emitted away from this region due to differential vorticity, which couples these two waves. The Rossby waves have their corotation radius at the density maximum.

Open with DEXTER

Two possibilities that can result in an extremum in $\mathcal{L}$ have been investigated: the first one occurs near the marginally stable orbit around a compact object, where relativistic effects create a maximum of $\kappa^2/2\Omega$. The growth rate of the RWI is strongly increased by a poloidal magnetic field threading the disk (Tagger & Varnière 2006), but decreased by a toroidal one (Yu & Li 2009). This has been proposed as an explanation for the high-frequency quasi-periodic oscillation (QPO) of microquasars (Tagger & Varnière 2006; Tsang & Lai 2009; Lai & Tsang 2009).

The second possibility is to have an extremum of the surface density. The model of Varnière & Tagger (2006) for protoplanetary disks relies on the assumption that extrema of $\Sigma$ should occur at the edges of the dead zone of these disks, so that the RWI should be unstable there. Tagger & Melia (2006) and Falanga et al. (2007) have also used the RWI in an explanation for the quasi-periodicity that may have been observed during the flares of SgrA*.

The RWI has been studied both analytically (Lovelace et al. 1999; Li et al. 2000) and numerically (Li et al. 2001; Varnière & Tagger 2006). It is formed by Rossby waves trapped in the extremum of $\mathcal{L}$, as shown in Fig. 1. In the $\beta$-plane approximation differential rotation is neglected and the dispersion relation of Rossby waves is given by

\begin{displaymath}\nonumber \omega=\frac{m\alpha}{k^2+m^2}, \end{displaymath}  

where k is the adimensional radial wavenumber, m the azimuthal wavenumber, and $\alpha$ the vorticity gradient. In an accretion disk, differential vorticity and differential rotation couple them to spiral waves, emitted on both sides of the extremum region. The wave dispersion relation calculated from the equations obtained by Lovelace et al. (1999) and Tagger (2001) is

\begin{displaymath}\nonumber \tilde\omega=\frac{c_{\rm s}^2m\alpha}{c_{\rm s}^2(k^2+m^2)+r^2\kappa^2}, \end{displaymath}  

where $\tilde\omega=\omega-m\Omega$, $c_{\rm s}$ is the sound speed and r the radius.

Rossby waves have their corotation radius (where their phase velocity equals the rotation velocity of the gas) at that extremum. The standing wave pattern they form appears as a vortex located in the region of the extremum, and grows exponentially. The waves have a positive energy flux and angular momentum beyond that radius and a negative energy flux inside it, so that (as will be checked in the simulation) the pattern can grow as it causes the gas inside the corotation to lose angular momentum and accrete, while the gas beyond the corotation gains that angular momentum and moves outward.

As a result the instability tends to destroy the extremum of $\mathcal{L}$, as seen for instance in Tagger & Melia (2006). In the context of the mechanism proposed by Varnière & Tagger (2006) for protoplanetary disks, we expect the instability to saturate at the amplitude where this is balanced by the continuous regeneration of the density extremum by the gas accreting from the outer disk and piling up at the edge of the dead zone (and reversely at the inner edge of the dead zone).

However, these studies have all been done in the approximation of an infinitely thin disk because of the complexity of a full 3D analytical study, or of the numerical resources needed for 3D simulations. Such a study is however highly desirable, both to consider the full complexity of the gas (and grains) flow, and to validate the thin disk approximation: this approximation was introduced (and its limits defined) by Goldreich & Lynden-Bell (1965b,a), but this does not apply as such to Rossby perturbations. We will see below that indeed 3D effects add significant and potentially important qualitative differences.

3 Gaseous accretion disk

3.1 Hydrodynamic equations

We work in cylindrical coordinates $(r, \phi, z)$ with the 3D Euler equation

$\displaystyle \partial_{\rm t} \rho +\vec \nabla \cdot (\vec v \rho)$ = 0 (2)
$\displaystyle \partial_{\rm t} (\rho \vec v) + \nabla \cdot (\vec v \rho \vec v)+\vec \nabla p$ = $\displaystyle -\rho \vec \nabla \Phi_{\rm G},$             (3)

where $\rho$ is the mass density of the fluid, and $\vec v$ its velocity, and p the pressure. $\Phi_{\rm G}={-GM_*}/{(r^2+z^2)^{1/2}}$ is the gravity potential of the central object with G the gravitational constant and M* the mass of the central object. We consider a barotropic flow, i.e. the entropy S is constant in the entire system

\begin{displaymath}\nonumber p=S\rho^\gamma, \end{displaymath}  

with the adiabatic index $\gamma=5/3$. The sound speed is given by $c_{\rm s}^2=\gamma p/\rho=S\gamma \rho^{\gamma-1}$ and the temperature by $T=p/\rho=S\rho ^{\gamma -1}$.

\begin{figure} \par\includegraphics[width=9cm,clip]{14000f2.eps} \end{figure} Figure 2:

Isocontours of the initial density, with a bump at $r_{\rm B}=3r_{\rm i}$. The density is normalized to the density in the midplane at the inner edge.

Open with DEXTER

\begin{figure} \par\includegraphics[width=9cm,clip]{14000f3.eps} \end{figure} Figure 3:

The initial azimuthal velocity of the midplane gas (continuous line) is compared to Keplerian rotation (dots). A zoom on the bump region is shown in the upper right corner.

Open with DEXTER

\begin{figure} \par\includegraphics[width=9cm]{14000f4.eps} \end{figure} Figure 4:

Plot of the logarithm of the square of the epicyclic frequency $\kappa ^2$ (averaged over z) as a function of the radius.

Open with DEXTER

3.2 Initial conditions

We choose as initial equilibrium a density profile decreasing radially as r-1/2, to which is added a density bump

\begin{displaymath}\frac{\rho^{\rm ini}(r,z=0)}{\rho_{\rm m}}=\left(1+\chi \exp{... ... i})^2}{2\sigma ^2}}\right)\frac{1}{\sqrt{r/r_{\rm i}}}\cdot \end{displaymath} (4)

The density is normalized by $\rho_{\rm m}$, the midplane density at the inner boundary of the simulation $r_{\rm i}$, and $r_{\rm B}$ is the position of the bump. For the parameters $\chi$ and $\sigma$, which control the amplitude and the width of the density bump, we take the values $\chi=0.4$ and $\sigma=0.1$ respectively, which gives a bump similar to that obtained by Varnière & Tagger (2006) in their simulation of the dead zone. These values allow the instability to grow relatively fast, decreasing the numerical load, while the effect of the pressure gradient on the equilibrium rotation velocity leaves us a margin with respect to the Rayleigh criterion $\kappa^{2}>0$ (see Fig. 4).

The bump is centered at $r_{\rm B}=3r_{\rm i}$, far enough from the inner edge of the disk to avoid a strong effect of the boundary condition there. The vertical density profile is initially at hydrostatic equilibrium, giving an aspect ratio of the order of 10-1,

$\displaystyle \rho^{\rm ini}(r,z)\!=\!\rho(r,z\!=\!0)\left[1\!-\!\frac{\gamma-1... ...{\rm i}}{r}-\frac{r_{\rm i}}{\sqrt{r^2+z^2}}\right)\right]^{\frac{1}{\gamma-1}}$   (5)

with and $S=10^{-3}r_{\rm i}^2\Omega_{\rm K}(r_{\rm i})^2\rho_{\rm m}^{1-\gamma}$. Finally, we initially use a ``floor'' density $\rho_{\min}=10^{-2}\rho_{\rm m}$ to avoid getting too low densities in the corona. Figure 2 shows the resulting isodensity contours in a vertical cut of the disk.

The initial velocity field is purely toroidal, $v_r^{\rm ini}=v_z^{\rm ini}=0$ and $v_{\phi}^{\rm ini}$ was chosen for the disk to be in radial equilibrium in a Newtonian potential

\begin{displaymath}v_{\phi}^{\rm ini}=\sqrt{ \frac{ r^2GM_* r_{\rm i}}{(r^2+z^2)^{3/2}}+\frac{r\rho_{\rm m}}{\rho r_{\rm i}}\nabla_r p}. \end{displaymath} (6)

Figure 3 shows the deviation of the disk azimuthal velocity from a pure Keplerian rotation. The pressure term in the last equation induces only a weak variation.

4 Numerical setup

4.1 Numerical code and scheme

The numerical simulations presented here use the versatile advection code (VAC) developped by Tóth (1996). In the version we use the code solves the 3D hydrodynamics equations for an isentropic flow. We use VAC with the total variation diminishing monotonic upstream-centered scheme for conservation laws (TVD-MUSCL), a Roe Riemann solver, a Hancock predictor step and a Woodward limiter (Colella & Woodward 1984). The TVD-MUSCL scheme detailed in Tóth & Odstrcil (1996) is one of the less dissipative schemes included in the VAC code, which is useful to observe the full development of the instability and properly characterize its saturation.

4.2 Grid and numerical resolution

The grid is cylindrical ($r,\phi,z$) with a resolution of $154\times68\times68$, including two ghosts cells at each boundary to impose the boundary conditions. The grid is non-uniform with a streching factor of 20 in the radial direction, and 10 in the vertical one.

This achieves a higher resolution in the regions of physical or numerical interest. The mesh extends from ( $r_{\rm i},0,0$) to ( $6r_{\rm i}$, $r_{\rm i}$, $2\pi$) where r and z have been normalized by the disk inner edge radius $r_{\rm i}$. The radial and vertical extensions are large enough for the different waves of the RWI to develop without unwanted effects from the boundary conditions.

The minimum value of the mass density is fixed to $10^{-6}\rho_{\rm m}$. With the initial conditions presented before this configuration gives a bump with approximatively 60 cells radially and 40 vertically.

4.3 Boundary conditions

The boundary conditions are imposed using the two ghost cells added at each boundary of the mesh, allowing us to compute the derivatives in all the cells. Because the simulation is restricted to a small part of the disk, we have chosen continuous boundary conditions in the radial direction, meaning that the values of the physical quantities are copied from the first cell of the computational domain in the ghost zones. This boundary condition is thus partially transparent. In the azimuthal direction the boundary condition is periodic. Only the upper part of the disk is simulated because the vortices we are interested in are even in z. In Sect. 5.2 we discuss a test run without this condition, which confirms the validity of this choice.

4.4 Disk equilibrium and stability

As mentioned in Sect. 3.2 we start the simulation with a low but finite density in the corona. On the other hand, we needed to start from a true disk equilibrium to avoid a rapid relaxation that would dominate the evolution of the system. The main difficulty was the constant density we chose to start from in the corona and the transition to the power law vertical profile in the disk; this introduced a jump in the vertical derivative of the density. We thus applied a smoothing over $n_{\rm smooth}$ vertical grid points (with $n_{\rm smooth}=4$ here) to this density profile to enforce a smooth transition. All this means that we are not at exact hydrostatic equilibrium in the gravity field. A different but related difficulty arose from the discretization of the grid and the use of slope limiters, which introduced residual spurious forces for each cell, but a zero mean value on the whole grid. Although they are small (at the level of roundoff error), they are systematic (i.e. always act in the same direction) and become important in low-density regions, especially at the disk-corona transition.

We thus chose to slightly modify the gravitational potential, so that the initial pressure profile given in Sect. 3.2 is really at equilibrium. We applied the following procedure: a first time-step was done to calculate the difference between the physical and numerical equilibrium. This allowed us to compute the source term $\rho \vec {\nabla}\Phi_{\rm G}$ in the Euler equations that enforced numerical equilibrium: this source term should be equal to the gravity term if we were in exact equilibrium, and the difference arises from these spurious forces. We enforced equilibrium by subsequently considering this constant modified source term. We will a posteriori check that this does not substantially affect our result. Moreover we checked that with the whole simulation box inside the disk (without corona) this modified source term is not needed, but wave reflection at the upper boundary makes it impossible to properly study the physics of the instability.

After several time-steps we added low level perturbations to these initial conditions to provide seeds for the unstable modes that can develop if the disk is unstable. We fond it more convenient to do this with perturbations with low azimuthal wavenumber, because (as confirmed in the simulation) they are the only ones expected to be unstable. We thus applied an initial radial velocity perturbation

$\displaystyle v_{r{\rm p}}=v_r+\epsilon \sin\big(2\pi(r-1.2)/.8\big)[\sin(\phi)+\sin(2\phi) +\sin(3\phi)+\sin(5\phi)]$   (7)

with the amplitude of the perturbation $\epsilon$ of the order of 10-7.

\begin{figure} \par\includegraphics[width=8.8cm,clip]{14000f5.eps} \end{figure} Figure 5:

Plot of the normalized density in the midplane of the disk at $t=176/\Omega _{\rm K}(r_{\rm i})$. One can identify the position of the crescent-shaped vortex at the bump radius ( $r=3r_{\rm i}$) and the spiral arms on both side of the bump region.

Open with DEXTER

5 Results

In this section we describe the results obtained with this simulation. We first describe the general properties of the instability, which allow us to identify it and to study its structure and evolution. We then focus on the velocity stream to understand the flow pattern induced by the instability in the disk.

\begin{figure} \par\includegraphics[width=7.65cm,clip]{14000f6.eps}\hspace*{1mm}... ...phics[width=7.65cm,clip]{14000f9.eps}\vspace*{-1mm}\vspace*{-2mm} \end{figure} Figure 6:

The compressional and rotational parts of the total flow ( left) and its non-axisymmetric part ( right) at $t=176/\Omega _{\rm K}(r_{\rm i})$, showing the two waves that form the instability. Top: the amplitude of $(\vec \nabla \times \vec v)_z$ traces the vorticity, i.e. the Rossby vortex centered at the extremum of density. Bottom: the amplitude of $\vec \nabla \cdot\vec v$traces the compressional waves (spiral density waves) that develop away from that region. In these figures, the whole disk is represented with the azimuthal coordinate on the vertical axis and the radius on the horizontal one. Note that although the top-right panel shows both a cyclonic and an anticyclonic vortex (as expected for an m=1 mode), once combined with the bulk flow only the anticyclonic vortex appears in the top-left panel.

Open with DEXTER

5.1 General study of the instability

Because the instability was only studied in 2D, we will first base the discussion on what we obtain in the midplane of the 3D simulation. As expected the simulation shows a coherent perturbed pattern forming in the region of the density bump after a few rotation times. Figure 5 shows the normalized density in the midplane. The non-axisymmetric pattern corresponds to the RWI, extending spiral arms on both sides of the bump.

As explained in Sect. 2, the structure of the instability is formed by Rossby waves on both sides of the extremum of $\mathcal{L}$ (here, the density bump at $r=3r_{\rm i}$), resulting in a standing vortex there, and spiral density waves emitted on both sides of that region. We can visualize this by plotting the curl and divergence of the flow, showing respectively its rotational (Rossby waves) and compressional (density waves) parts. Figure 6 shows that the expected patterns do indeed appear.

The different stages of the development of the instability can be seen Fig. 7, which presents the amplitude of the density perturbation on a logarithmic scale as a function of time. After a short stage where the non-pertinent perturbations decrease, the unstable mode enters the linear stage, i.e. the perturbation grows exponentially. We find a growth rate of $.39\Omega_{\rm K}(r=r_{\rm i})/\Omega_{\rm K}(r=3r_{\rm i})$. This is comparable with the growth rates obtained from 2D theory (Li et al. 2000) or numerical simulation (Tagger & Melia 2006). The different physical and numerical setups (including the equation of state and the boundary conditions) prevent a more detailed comparison. The local dispersion equation (Eq. (3)) does not allow us to estimate the growth rate of the global modes. This would anyway be possible only with a fully 3D linear computation, which would require an extensive effort. The perfect exponential growth we obtain in the present non-linear simulation over four decades in amplitude, its coherence with the 2D results, and its independence from initial conditions show that we do capture here the linear phase of the instability properly. This stage lasts about 10 Keplerian orbits at the bump radius, and finally the instability reaches saturation. The entire progression of the instability is thus captured in this simulation.

\begin{figure} \par\includegraphics[width=9cm,clip]{14000f10.eps}\vspace*{2mm} \end{figure} Figure 7:

Instability growth: the maximum value of the non-axisymmetric part of the density is plotted as a function of time with a logarithmic scale for the vertical axis. We added here a second time axis (lower) that shows the number of Keplerian orbits at the bump radius to estimate the growth rate in the same units as Li et al. (2000). The equation of the continuous line is $-4.1+0.39t\Omega _{\rm K}(r_{\rm i}) / \Omega _{\rm K}(3r_{\rm i})$. The graph thus shows the linear stage and the saturation after about 10 Keplerian times at $r=3r_{\rm i}$.

Open with DEXTER

Figure 9 shows that the saturation is not due to the flattening of the initial density bump and thus of the extremum of $\mathcal{L}$, as in the 2D simulations of Tagger & Melia (2006), so that non-linearities must be suspected. In order to check this we ran a new simulation, taking as initial radial density profile the average one at the end of the present simulation (but again with only weak velocity perturbation, as at its startup), and observed the RWI growing again.

The azimuthal density profile at the bump radius is shown in Fig. 10. It reveals a departure from a pure sinusoid, corresponding to the presence of m > 1 modes, whose growth are shown in Fig. 8. Because they do not show the exponential growth of the m = 1 but do show modes that are not present in the initial perturbation, we analyze them as non-linearly generated harmonics of that mode. Figure 12 shows that their effect remains relatively weak, so that the saturation of the mode cannot be ascribed to them.

\begin{figure} \par\includegraphics[width=18cm,clip]{14000f11.eps} \end{figure} Figure 8:

The evolution of the amplitude of the ten first azimuthal modes in the same coordinates as Fig. 7. The amplitude of each mode is calculated with an azimuthal Fourier transform, $\rho _{\rm i}$ being the amplitude of the ith mode. The development of modes that were not present in the initial perturbation proves the existence of non-linearities.

Open with DEXTER

On the other hand, as shown in Fig. 13, small-scale structures appear in the vertical flow, and become stronger when the mode approaches saturation. They cascade to the smallest scale of the simulation, making us believe that dissipation in these structures is the main cause of the saturation we observe. Assessing their role will require both improved numerics, to deal with these structures, and physical understanding because one can expect that vertically propagating sound waves generated within the disk will be subject to wavebreaking when they reach the low-density, low-sound speed corona. This will be considered in future work.

\begin{figure} \par\includegraphics[width=9.cm]{14000f12.eps} \end{figure} Figure 9:

Comparison between the density integrated over z and $\phi $ ( left) and of the RWI criterion $\mathcal{L}$ ( right) at t=0 in solid line and $t=176/\Omega _{\rm K}(r_{\rm i})$ in dashed line. One can see that the extremum of $\mathcal{L}$ at $r=3r_{\rm i}$ has decreased but not disappeared when the instability saturates.

Open with DEXTER

The flattening of the bump (Fig. 9) is explained by the accretion rate, shown in Fig. 11: it is positive inward from the density extremum and negative beyond it. This is expected from a mode with its corotation at that radius which grows by exchanging angular momentum between the inner and outer region. As in 2D simulations, the transition between positive and negative accretion is so sharp that it is actually the best diagnostic of the corotation radius and thus of the mode frequency.

Finally, considering the possibility of other unstable modes, we show in Fig. 12 the time evolution of the azimuthal Fourier components, m=0 to 3. Although perturbations do appear at m=2 and 3, that they evolve together with the m=1 indicates that they are probably harmonics of that mode and not independently growing ones. The most linearly unstable m depends on the local conditions. We ran 2D simulations with the same density profile and found in that case a dominant m=2. We will see below that the 3D simulation does show new features in the structure of the mode, which seem to contribute to its linear growth. They might be at the origin of this difference.

5.2 Study of the velocity stream

In this section we consider the flow pattern generated by the instability. Because the stream in the $(r,\phi)$plane is dominated by Keplerian rotation, we consider the non-axisymmetric part of the velocity, i.e. we study the disk in a sheared frame rotating as the mean flow at each radius. Figure 13 does show in this plane the expected vortex structure, with both a cyclonic and anticyclonic feature, because the m=1 structure generates both. As seen in Fig. 5 however, when combined with the Keplerian flow only the anticyclonic component appears, forming a crescent-shaped feature, while the cyclonic one is at the x-point where the tips of the crescent join.

\begin{figure} \par\includegraphics[width=9.cm]{14000f13.eps} \end{figure} Figure 10:

The azimuthal surface density profile at $r=3r_{\rm i}$ at different times with the same normalization. The line is at $t=32/\Omega _{\rm K}(r_{\rm i})$, the departure from the sinusoid comes from the initial perturbations that is a sum of sinus. The dots is at $t=104/\Omega _{\rm K}(r_{\rm i})$, the profile is closer to a sinusoid, the instability is in the linear phase, the nonpertinent modes (m>1) are negligible. The dashed line is at $t=176/\Omega _{\rm K}(r_{\rm i})$, the departure from a sinusoid is believed to be a direct consequence of non-linearities.

Open with DEXTER

\begin{figure} \par\includegraphics[width=9cm,clip]{14000f14.eps} \end{figure} Figure 11:

The accretion rate defined as $- \int r\rho v_r {\rm d}z {\rm d}\phi / \!\!\int r\rho_{\rm m} c_{\rm s} {\rm d}z {\rm d}\phi$ as a function of radius at the end of the linear phase ( $t=176/\Omega _{\rm K}(r_{\rm i})$).

Open with DEXTER

On the other hand Fig. 13 also shows unexpected new features: first, the plot in the ($\phi,\ z$) plane shows an downdraft at the center of the anticyclonic vortex, and a updraft at the center of the cyclonic one. Also, plots in the ($r,\ z$) plane show that strong convection rolls form on both sides of the vortex, in the anticlockwise direction for the outer one and clockwise direction for the inner one. Streamlines at the outer edges of these horizontal and vertical vortices connect them smoothly in a manner reminiscent of orbits near separatrices in hamiltonian dynamics. The main streamlines are schematized in Fig. 14.

\begin{figure} \par\includegraphics[width=13.5cm,clip]{14000f15.eps} \end{figure} Figure 12:

Temporal evolutions of the first four azimuthal modes. A different color bar is used for m=0 which is approximately ten times stronger than the others.

Open with DEXTER
\begin{figure} \mbox{\includegraphics[width=8.4cm,clip]{14000f16.eps} \includeg... ...14000f18.eps} \includegraphics[width=8.4cm,clip]{14000f19.eps} } \end{figure} Figure 13:

Top: velocity streamlines in a frame rotating with the disk in the midplane of the disk (left) and in the vertical plane at $r=3r_{\rm i}$ (right). Bottom: velocity streamlines in a vertical frame at $\phi /2\pi =0.64$ (left) and 0.66 (right).

Open with DEXTER

We checked the robustness of these convection rolls by performing numerical simulations with the entire vertical structure of the disk and without assuming a vertical symmetry, all the other parameters kept unchanged. The vertical symmetry was broken by the initial perturbation added on the radial velocity, but resulted in similar flow patterns. Furthermore, when we applied only antisymmetric initial perturbations, we found that this delayed the development of the instability without changing its vertical structure. We conclude that the delay is due to the difficulty for antisymmetric initial perturbations to seed the unstable mode and that the convection rolls are part of the linear structure of the instability. To confirm this diagnostic, we also checked that the maximal vertical velocity involved in these convection rolls follows the linear growth of the instability as in Fig. 7.

This strong vertical structure of the instability and of the vortex pattern could not be expected from the 2D results. We believe that the vertical vortices are resonantly excited where the local Doppler-shifted frequency, $\omega-m\Omega(r)$, is equal to the mean vertical oscillation frequency. The relation between this and the Brünt-Väisälä frequency, depending on the equation of state of the gas, will be discussed elsewhere. This would thus be the equivalent, in a fluid disk, of the ``peanuts''  structure found in observations and numerical simulations of barred spiral galaxies (Athanassoula 2008; Combes et al. 1990; Combes & Sanders 1981). This is believed (although an explanation based on the firehose instability has also been proposed) to result from the resonant excitation of the vertical motion of stars where $\omega-m\Omega(r)$ equals their vertical frequency of motion.

6 Discussion and outlook

We presented here the first fully three-dimensional, cylindrical numerical simulation showing a strong and persistent vortex in a differentially rotating disk. Its growth results from the Rossby wave instability, which has been invoked in various astrophysical contexts like in protoplanetary disks or in the accretion disk of compact objects.

The simulation essentially confirms the main properties of the instability, which forms a vortex extending across the vertically stratified structure of the disk. The vortex is localized at the extremum of vortensity and very efficiently acts to destroy this extremum, by permitting an exchange of angular momentum between regions located on either side of it. In our simulation we found no sign of the elliptical instability (Kerswell 2002; Lesur & Papaloizou 2009) which has been claimed to destroy vortices. We note however that these works were done in the shearing sheet approximation, which neglects the equilibrium vorticity gradient and thus can consider local small-scale modes, whereas the RWI is a large-scale global instability. The same remark applies to the work of Latter & Balbus (2009) and Lesur & Papaloizou (2010).

The most important new result found in our 3D simulation is the occurrence of strong vertical convection rolls, excited on both sides of the horizontal vortex. Figure 13 shows that the vertical velocity in these rolls is comparable to the radial velocity involved in the Rossby vortex, while their growth shows that they are an inherent part of the structure of the mode, tightly connecting the horizontal and vertical circulations.

Future work should allow us to assess the importance of this convection in the astrophysical contexts where the RWI may act. In particular we expect them to play a very important role in the concentration of dust grains in protoplanetary disks, in addition to the mechanism of Barge & Sommeria (1995), which will accumulate them at the center of the horizontal vortex. Indeed Fig. 13 also shows that in the midplane of the disk, the flow lines rapidly spiral inward in the cyclonic vortex and outward in the anticyclonic one. This coincides with the updraft at the cyclonic vortex and downdraft at the anticyclonic one. We note that dust grains will thus be rapidly transported with the gas, at low z, toward the center of the cyclonic vortex, but that their weight should make it very difficult for them to be dragged in the updraft. The 3D structure of the vortex presented here will then have a direct impact on the accumulation of grains that will happen on an even lower timescale than in a 2D vortex. We also note however that counter to both intuition and the effect discussed by Barge & Sommeria (1995), this would tend to accumulate the grains at the cyclonic rather than anticyclonic vortex.

\begin{figure} \par\includegraphics[width=8.8cm,clip]{fig20.ps} \end{figure} Figure 14:

Schematic view of the circulation, showing the horizontal and vertical vortices in a sheared frame rotating at the local mean angular velocity.

Open with DEXTER

In order to explore this mechanism, which could be of primordial importance for the growth of planetesimals and planet formation, work should proceed in two directions: performing simulations using varied density profiles, including that expected at the edge of the dead zone of protoplanetary disks (Varnière & Tagger 2006), and following the motion and growth of dust grains in the resulting vertical flows. An analytical description of the vertical structure observed in the simulation is also needed.

Another direction for future work concerns the initial goal of this work, which was to analyze the 3D structure of MHD instabilities in magnetized disks.

In accretion disks threaded by a vertical magnetic field, it has been shown that an accretion-ejection instability (AEI) can occur (Tagger & Pellat 1999) and explain the low-frequency quasi-periodic oscillation of X-ray binaries (Varnière et al. 2002; Rodriguez et al. 2002). The instability grows by coupling magnetically-driven spiral density waves and a Rossby vortex.It has also been shown (Varnière & Tagger 2002) that this vortex can re-emit vertically, as Alfvén waves propagating to the corona along magnetic field lines, a substantial fraction of the accretion energy and angular momentum; it was also suggested that this could be a source for accretion-driven winds and jets. In this context we can expect that the convection rolls observed here would be replaced by resonantly excited slow magnetosonic waves. Because the main action of these waves is to move gas along magnetic field lines we can expect them to provide mass-load to the resulting jet, much as in MHD steady-state models of jets (Ferreira & Pelletier 1995; Casse & Ferreira 2000) mass-loading occurs at the slow magnetosonic point and acceleration occurs in the region of the Alfvénic point. In order to study this prospect a fully 3D simulation that can handle both the geometry (a disk threaded by open poloidal magnetic field lines) and the magnitude (on the order of equipartition with the gas pressure) of the relevant magnetic field needs to be developed.

Acknowledgements
This work was granted access to the HPC resources of IDRIS under the allocation 2009-i2009042125 made by GENCI (Grand Equipement National de Calcul Intensif).

References

  1. Athanassoula, E. 2008, in IAU Symp. 245, ed. M. Bureau, E. Athanassoula, & B. Barbuy, 93 [Google Scholar]
  2. Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
  3. Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [NASA ADS] [CrossRef] [Google Scholar]
  4. Casse, F., & Ferreira, J. 2000, A&A, 353, 1115 [NASA ADS] [Google Scholar]
  5. Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  6. Combes, F., & Sanders, R. H. 1981, A&A, 96, 164 [NASA ADS] [Google Scholar]
  7. Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82 [NASA ADS] [Google Scholar]
  8. Falanga, M., Melia, F., Tagger, M., Goldwurm, A., & Bélanger, G. 2007, ApJ, 662, L15 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807 [NASA ADS] [Google Scholar]
  10. Goldreich, P., & Lynden-Bell, D. 1965a, MNRAS, 130, 97 [NASA ADS] [CrossRef] [Google Scholar]
  11. Goldreich, P., & Lynden-Bell, D. 1965b, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
  12. Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  14. Kerswell, R. R. 2002, Annual Review of Fluid Mechanics, 34, 83 [Google Scholar]
  15. Lai, D., & Tsang, D. 2009, MNRAS, 393, 979 [NASA ADS] [CrossRef] [Google Scholar]
  16. Latter, H. N., & Balbus, S. A. 2009, MNRAS, 399, 1058 [NASA ADS] [CrossRef] [Google Scholar]
  17. Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874 [NASA ADS] [CrossRef] [Google Scholar]
  20. Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] [Google Scholar]
  21. Lovelace, R. V. E., & Hohlfeld, R. G. 1978, ApJ, 221, 51 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 491, L41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009, A&A, 497, 869 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Papaloizou, J. C. B., & Pringle, J. E. 1985, MNRAS, 213, 799 [NASA ADS] [Google Scholar]
  26. Rodriguez, J., Varnière, P., Tagger, M., & Durouchoux, P. 2002, A&A, 387, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Sellwood, J. A., & Kahn, F. D. 1991, MNRAS, 250, 278 [NASA ADS] [Google Scholar]
  28. Tagger, M. 2001, A&A, 380, 750 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Tagger, M., & Melia, F. 2006, ApJ, 636, L33 [NASA ADS] [CrossRef] [Google Scholar]
  30. Tagger, M., & Pellat, R. 1999, A&A, 349, 1003 [NASA ADS] [Google Scholar]
  31. Tagger, M., & Varnière, P. 2006, ApJ, 652, 1457 [NASA ADS] [CrossRef] [Google Scholar]
  32. Tóth, G. 1996, Astrophys. Lett. Comm., 34, 245 [Google Scholar]
  33. Tóth, G., & Odstrcil, D. 1996, J. Comput. Phys., 128, 82 [Google Scholar]
  34. Tsang, D., & Lai, D. 2009, MNRAS, 400, 470 [NASA ADS] [CrossRef] [Google Scholar]
  35. Varnière, P., & Tagger, M. 2002, A&A, 394, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Varnière, P., & Tagger, M. 2006, A&A, 446, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Varnière, P., Rodriguez, J., & Tagger, M. 2002, A&A, 387, 497 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Yu, C., & Li, H. 2009, ApJ, 702, 75 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{14000f1.eps} \end{figure} Figure 1:

Schematic view of the Rossby wave instability. Rossby waves are generated in the region of the density extremum, and density waves are emitted away from this region due to differential vorticity, which couples these two waves. The Rossby waves have their corotation radius at the density maximum.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{14000f2.eps} \end{figure} Figure 2:

Isocontours of the initial density, with a bump at $r_{\rm B}=3r_{\rm i}$. The density is normalized to the density in the midplane at the inner edge.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{14000f3.eps} \end{figure} Figure 3:

The initial azimuthal velocity of the midplane gas (continuous line) is compared to Keplerian rotation (dots). A zoom on the bump region is shown in the upper right corner.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm]{14000f4.eps} \end{figure} Figure 4:

Plot of the logarithm of the square of the epicyclic frequency $\kappa ^2$ (averaged over z) as a function of the radius.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{14000f5.eps} \end{figure} Figure 5:

Plot of the normalized density in the midplane of the disk at $t=176/\Omega _{\rm K}(r_{\rm i})$. One can identify the position of the crescent-shaped vortex at the bump radius ( $r=3r_{\rm i}$) and the spiral arms on both side of the bump region.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=7.65cm,clip]{14000f6.eps}\hspace*{1mm}... ...phics[width=7.65cm,clip]{14000f9.eps}\vspace*{-1mm}\vspace*{-2mm} \end{figure} Figure 6:

The compressional and rotational parts of the total flow ( left) and its non-axisymmetric part ( right) at $t=176/\Omega _{\rm K}(r_{\rm i})$, showing the two waves that form the instability. Top: the amplitude of $(\vec \nabla \times \vec v)_z$ traces the vorticity, i.e. the Rossby vortex centered at the extremum of density. Bottom: the amplitude of $\vec \nabla \cdot\vec v$traces the compressional waves (spiral density waves) that develop away from that region. In these figures, the whole disk is represented with the azimuthal coordinate on the vertical axis and the radius on the horizontal one. Note that although the top-right panel shows both a cyclonic and an anticyclonic vortex (as expected for an m=1 mode), once combined with the bulk flow only the anticyclonic vortex appears in the top-left panel.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{14000f10.eps}\vspace*{2mm} \end{figure} Figure 7:

Instability growth: the maximum value of the non-axisymmetric part of the density is plotted as a function of time with a logarithmic scale for the vertical axis. We added here a second time axis (lower) that shows the number of Keplerian orbits at the bump radius to estimate the growth rate in the same units as Li et al. (2000). The equation of the continuous line is $-4.1+0.39t\Omega _{\rm K}(r_{\rm i}) / \Omega _{\rm K}(3r_{\rm i})$. The graph thus shows the linear stage and the saturation after about 10 Keplerian times at $r=3r_{\rm i}$.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=18cm,clip]{14000f11.eps} \end{figure} Figure 8:

The evolution of the amplitude of the ten first azimuthal modes in the same coordinates as Fig. 7. The amplitude of each mode is calculated with an azimuthal Fourier transform, $\rho _{\rm i}$ being the amplitude of the ith mode. The development of modes that were not present in the initial perturbation proves the existence of non-linearities.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9.cm]{14000f12.eps} \end{figure} Figure 9:

Comparison between the density integrated over z and $\phi $ ( left) and of the RWI criterion $\mathcal{L}$ ( right) at t=0 in solid line and $t=176/\Omega _{\rm K}(r_{\rm i})$ in dashed line. One can see that the extremum of $\mathcal{L}$ at $r=3r_{\rm i}$ has decreased but not disappeared when the instability saturates.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9.cm]{14000f13.eps} \end{figure} Figure 10:

The azimuthal surface density profile at $r=3r_{\rm i}$ at different times with the same normalization. The line is at $t=32/\Omega _{\rm K}(r_{\rm i})$, the departure from the sinusoid comes from the initial perturbations that is a sum of sinus. The dots is at $t=104/\Omega _{\rm K}(r_{\rm i})$, the profile is closer to a sinusoid, the instability is in the linear phase, the nonpertinent modes (m>1) are negligible. The dashed line is at $t=176/\Omega _{\rm K}(r_{\rm i})$, the departure from a sinusoid is believed to be a direct consequence of non-linearities.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{14000f14.eps} \end{figure} Figure 11:

The accretion rate defined as $- \int r\rho v_r {\rm d}z {\rm d}\phi / \!\!\int r\rho_{\rm m} c_{\rm s} {\rm d}z {\rm d}\phi$ as a function of radius at the end of the linear phase ( $t=176/\Omega _{\rm K}(r_{\rm i})$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=13.5cm,clip]{14000f15.eps} \end{figure} Figure 12:

Temporal evolutions of the first four azimuthal modes. A different color bar is used for m=0 which is approximately ten times stronger than the others.

Open with DEXTER
In the text

  \begin{figure} \mbox{\includegraphics[width=8.4cm,clip]{14000f16.eps} \includeg... ...14000f18.eps} \includegraphics[width=8.4cm,clip]{14000f19.eps} } \end{figure} Figure 13:

Top: velocity streamlines in a frame rotating with the disk in the midplane of the disk (left) and in the vertical plane at $r=3r_{\rm i}$ (right). Bottom: velocity streamlines in a vertical frame at $\phi /2\pi =0.64$ (left) and 0.66 (right).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig20.ps} \end{figure} Figure 14:

Schematic view of the circulation, showing the horizontal and vertical vortices in a sheared frame rotating at the local mean angular velocity.

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.