---

On the stability of elliptical vortices in accretion discs

G. Lesur - J. C. B. Papaloizou

Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK

Received 23 December 2008 / Accepted 3 March 2009

Abstract
Context. The existence of large-scale and long-lived 2D vortices in accretion discs has been debated for more than a decade. They appear spontaneously in several 2D disc simulations and they are known to accelerate planetesimal formation through a dust trapping process. In some cases, these vortices may even lead to an efficient way to transport angular momentum in protoplanetary discs when MHD instabilities are inoperative. However, the issue of the stability of these structures to the imposition of 3D disturbances is still not fully understood, and it casts doubts on their long term survival
Aims. We present new results on the 3D stability of elliptical vortices embedded in accretion discs, based on a linear analysis and several non-linear simulations.
Methods. We introduce a simple steady 2D vortex model which is a non-linear solution of the equations of motion, and we show that its core is made of elliptical streamlines. We then derive the linearised equations governing the 3D perturbations in the core of this vortex, and we show that they can be reduced to a Floquet problem. We solve this problem numerically in the astrophysical regime, including a simplified model to take into account vertical stratification effects. We present several analytical limits for which the mechanism responsible for instability can be explained. Finally, we compare the results of the linear analysis to some high resolution numerical simulations obtained with spectral and finite difference methods. A discussion is provided, emphasising the astrophysical consequences of our findings for the dynamics of vortices.
Results. We show that most anticyclonic vortices are unstable due to a resonance between the turnover time and the local epicyclic oscillation period. A small linearly stable domain is found for vortex cores with an aspect-ratio of around 5. However, our simulations show that it is only the vortex core that is stable, with the instability still appearing on the vortex boundary. In addition, we find numerically that results obtained under the assumption of incompressibility are not affected by the introduction of a moderate compressibility. Finally, we show that a strong vertical stratification does not create any additional stable domain of aspect ratio, but it significantly reduces growth rates for relatively weak (and therefore elongated) vortices.
Conclusions. Elliptical vortices are always unstable, whatever the horizontal or vertical aspect-ratio is. The instability can however be weak and is often found at small scales, making it difficult to detect in low-order finite-difference simulations.

Key words: accretion, accretion disks - instabilities - hydrodynamics

1 Introduction

The existence of 2D long-lived vortices in accretion discs was first proposed by von Weizsäcker (1944) in an outmoded model of planet formation. This idea was revived by Barge & Sommeria (1995) to accelerate planetesimals formation by a dust trapping process. This kind of vortex is often observed in 2D simulations of discs (see e.g. Bodo et al. 2007; Godon & Livio 1999; Umurhan & Regev 2004; Johnson & Gammie 2005b), since 2D turbulence is known to generate an inverse cascade of energy leading to large 2D vortices (Onsager 1949). Vortices may also be generated by 2D instabilities such as Rossby wave instabilities (Lovelace et al. 1999) or baroclinic instabilities (Klahr & Bodenheimer 2003; Petersen et al. 2007), although the latter is still a matter of ongoing debate (Johnson & Gammie 2005a,2006). These vortices may play at least two important roles regarding accretion disc dynamics. First, they could lead to an efficient angular momentum transport process in regions in which the magneto-rotational instability (Balbus & Hawley 1998) doesn't operate, such as in dead zones (Gammie 1996). Second, they are a very efficient way to accelerate the planetesimal formation process in protoplanetary discs (Johansen et al. 2004). However, the stability of these vortices when small 3D disturbances are imposed is largely unknown.

In the astrophysics community, this issue has been investigated mainly numerically. Shen et al. (2006) examined the formation of 2D vortices starting from 2D turbulence in fully compressible simulations. According to their results, a small 3D noise added to their initialy 2D configuration destroys the coherent vortices in a few orbits, relaxing the flow to its laminar state. Barranco & Marcus (2005) also computed the evolution of 3D vortices using an anelastic code incorporating vertical stratification. As Shen et al. (2006), they found that midplane vortices were destroyed by 3D perturbations. However, they also showed that off-midplane vortices could survive for several hundreds of orbits, leading to the possibility of a stabilizing effect due to the stratification

It is often assumed that these vortices are unstable because of the elliptical instability. The elliptical instability is a parametric instability appearing when a multiple of the vortex turnover frequency matches an inertial wave frequency, leading to a positive resonance. It is observed when the backgound flow follows closed streamlines, and being localized on individual streamlines is a local instability (in particular it doesn't need to involve the vortex boundaries). This instability was first found numerically by Pierrehumbert (1986) and described using Craik & Criminale (1986) solutions by Bayly (1986) for pure elliptical flows. The rotating case was studied by Craik (1989), who showed that anticyclonic elliptical flows can be stable for some rotation rates. Interested readers may consult Kerswell (2002) for a more extensive discussion of the elliptical instability and its development in fluid mechanics.

In the present paper, we investigate the elliptical instability in the context of accretion disc vortices. We first present a steady 2D vortex model, which is a non-linear solution of the local disc equations. We then present the linearised equations governing 3D perturbations inside the vortex. A criterion for the instability is derived from these equations and a physical understanding of the mechanism responsible for the instability is provided. We briefly extend these results to a simplified stratified case, and we compare our findings to fully non-linear simulations of accretion discs vortices. Finally, we provide a discussion and a comparison with previous work.

2 Local model of an elliptical vortex

2.1 Shearing-sheet model

In the following, we will assume a local model for the accretion disc, following the shearing-sheet approximation. The reader may consult Hawley et al. (1995), Balbus (2003) and Regev & Umurhan (2008) for an extensive discussion of the properties and limitations of this model. As a simplification, we will assume the flow is incompressible, consistently with the small shearing box model (Regev & Umurhan 2008). The shearing box equations are found by considering a Cartesian box centred at r=R₀, rotating with the disc at angular velocity . We define and for consistency with the standard notation for plane Couette flows (e.g. Drazin & Reid 1981). Note that this definition differs from the standard notation used in shearing boxes (Hawley et al. 1995) with , and . In this rotating frame, one obtains the following set of governing equations

In these equations, we have defined the mean shear , which is set to for a Keplerian disc. The generalised pressure is calculated solving a Poisson equation with the incompressibility condition. One can check easely that the velocity field is a steady solution of these equations.

2.2 The Kida (1981) solution

We want to study the stability of a steady elliptical vortex embedded in the sheared flow described in the previous subsection. Since we are looking for 2D solutions in the (x,y) plane, we can omit the Coriolis force as it will only change the pressure distribution. We define this vortex by an elliptical patch of constant vorticity (the ``core'') where -S is the background flow vorticity and is the vorticity of the vortex itself. Outside of this core, the vorticity is assumed to be , extending to infinity. According to Kida (1981) (Eq. (2.9)), such a vortex is steady if the semi-major axis is aligned with x and if its vorticity satisfies

where we have defined the vortex aspect-ratio , a and b being respectively the vortex semi-major and semi-minor axis. One deduces from this result that only anticyclonic vortices ( and having opposite sign) are stable in Keplerian flows. Since this solution is steady, no streamline goes through the core boundaries, and the streamlines inside the core have to be elliptical, with the same aspect-ratio as the vortex core.

Thanks to this property, we can write the velocity field in the vortex core, assuming it's centered on x=y=0, as

This solution can be written in the simpler form

defining

2.3 Explicit solution

A complete solution for the velocity field can be found defining the streamfunction so that and . This streamfunction satisfies:

One solves these equations in elliptical coordinates, defining  by

with . In these coordinates, the core boundary is found at with . Requiring that  and are continuous at , one finds the following expressions for :

where the subscripts i and o stand for inside and outside the core.

As expected, the inner solution reproduces the elliptical flow described in the previous subsection (Eqs. (4)-(5)). In the outer solution, one recognises the background shear (3rd term), which dominates for large . Moreover, this solution exhibits a linear term in  corresponding to the long-range perturbation of the vortex. In cartesian coordinates, this linear dependance translates into a logarithmic tail for , showing that the vortex presence can be felt far from the core. As we will see in Sect. 5 below, this property leads to numerical artefacts when one tries to fit this solution in a finite-size box. For completeness, we show in Fig. 1 the streamlines obtained from solution (11)-(12). As expected, the streamlines inside the core are ellipses of constant aspect ratio . Outside the core, one still finds closed streamlines but with a much more elongated structure.

2.4 Perturbation equations

In the following, we concentrate on the evolution of perturbations in the vortex core, assuming the latter is infinite. This corresponds to a limit in which the perturbations are small compared to the typical horizontal size of the core. To study the evolution of 3D perturbations in a 2D elliptical flow, we write the total velocity field as where is supposed to have an infinitely small amplitude. Following Kelvin (1880) and Craik & Criminale (1986), we use a time explicit Fourier decomposition for the perturbation . Plugging this solution in (1) leads at first order to:

where we have included the tidal term of (1) in . To satisfy (13) for any x_k, one has to solve:

which leads to:

where k₀, and t₀ are integration constants and is the turnover angle. Therefore, the perturbations will have a rotating wavevector with a turnover time equal to the vortex turnover time . Because of the incompressibility condition, will correspond to horizontal motion perturbations whereas will imply vertical ones. We also note that the rotating wavevector will involve smaller structures in the y direction. According to this result, the horizontal aspect ratio of the perturbation wavevector is equal to that of the vortex ().

We can then simplify (13) eliminating the pressure, which leads to the final system

where and . Plugging solution (16) in this equation leads to a Floquet problem for v, as already pointed out by Bayly (1986). Interestingly, the Floquet problem doesn't depend on the norm of the wavevector but just on its direction. Therefore, this problem is scale-independant, at least in the inviscid limit. The solution to this problem is known to be a superposition of Floquet modes written

where f is periodic with period . To determine the Floquet exponents , we compute the eigenvalues of the matrix where  satisfies the generalised equation

with the initial condition

A numerical approach is required to solve this problem in most cases. However, some limits can be understood using analytical approaches, as we will see in the next section.

3 The elliptical instability

3.1 Horizontal instability

It is possible to derive an analytical criterion for the instability in the limit . In this limit, v_z=0 and (17) is reduced to

This system describes horizontal epicyclic oscillations with frequency

having defined the rotation number as . In the limit where the vortex is infinitely weak, we find the classical epicyclic frequency , which is equal to the Keplerian frequency in discs. This epicyclic frequency is imaginary when

If a columnar vortex is in this regime, its horizontal layers will tend to drift exponentially in the (x,y) plane, leading to the destruction of the vertical coherence of the structure. In accretion disc vortices, this regime is found for .

3.2 A generalization

We note that the above discussion may be generalized to apply to a wider class of steady flows (u_x⁰ , u_y⁰). The linearized equations of motion are

where is the pressure perturbation. The operator

denotes the convective derivative. We now assume perturbations are local in z so that we may assume that the z dependence is through a factor where k_zis very large in the sense that the length scale k_z^-1 is smaller than any other in the problem. From the incompressibility condition , and the z component of (25) we conclude that . Thus it may be dropped from the x and y components of (25) which then give the pair

Replacing by with the understanding that the time derivative is to be taken on a fixed streamline, we see that in general the system (27)-(28) leads to a second order ordinary differential equation with periodic coefficients, the period being the time to circulate round the chosen streamline. But note that for the problem on hand the coefficients are constant because the unperturbed velocity is a linear function of the coordinates. In more general cases one has a Floquet problem of the type described above for every streamline indicating that unstable modes are indeed localized on streamlines. In order to connect with (24) we note that one can prove that if everywhere on a chosen streamline

there will be an instability. In fact for the problem on hand this is precisely equivalent to the above condition that which direcrly leads to (24). To show this we first note that if

then both and must be either positive or negative. Without loss of generality we assume both to be positive (the case when they are both negative may be recovered by setting ). Then from Eqs. (27)-(28) we may obtain

or setting q_y=-v_y

Thus if and are the minimum values of and respectively, we have

or equivalently

Thus it follows that v_xq_y grows faster than exponentially with growth rate which means that there is a Floquet exponent corresponding to exponential growth.

3.3 Numerical results

In the general case ( and ), the Floquet problem (17) can't be solved analytically. We therefore use the following numerical approach. We first evolve the M matrix following Eq. (19) for a given  and  with a 4th order Runge-Kutta-Fehlberg algorithm, keeping the error down to 10^-10. The eigenvalues are then computed from  using the GNU scientific library routines. This procedure is repeated for 1000 values of  and 1000 values of  to get a 2D representation of the regions for which the instability is found.

To test our numerical results, we have first reproduced the results of Bayly (1986) for a vortex in a non rotating sheared flow. One find in this case that vortices are unconditionally unstable to fully 3D disturbances (Fig. 2). Moreover, in the limit , only modes with are found to be unstable, consistently with Bayly's findings. Finally, we note that the growth rate for the instability weakens as we elongate the vortex. This finding, although apparently opposite to the classical result for the elliptical instability, is due to the fact that our growth rate is measured in inverse shear time and not in inverse vortex turnover time.

We have used our numerical method to solve the Floquet problem in the Keplerian case (Fig. 3). In this case, we observe essentially a two bands structure. The first band (low- band) is located between  and . It is associated with a short growth time (typically of the order of one shear time), and it includes the horizontal instability described in Sect. 3.1 for . In the limit (infinitely strong vortex), we find the same result as in the non rotating case (instability for ), which was to be expected since the vortex turnover time is much smaller than the rotation time in this limit. The second band (high- band) is found for and is much weaker since it involves growth rate of the order of 10^-2 S. This band gets wider and weaker as we go to large aspect ratios (and weak vortices), with a growth rate maximum found for We have tried to get a larger resolution in the region -6 and it seems that this band reaches for , but the growth rates involved are very small. Note that another narrower band is observed in the regime of high and smaller in Fig. 3. As we will see in the following, this band is due to higher order resonance and is therefore weaker than the bands described previously.

To get a simpler representation of the elliptical instability in the Keplerian case, we have plotted the maximum growth rate as a function of  in Fig. 4. We observe the same basic features as on the colormaps with 2 bands. We also observe that the high- band decays as gets larger. This property is to be expected since the limit corresponds to a pure shear flow, which is known to be linearly stable. However, for a finite (but large) , vortices are always unstable, contrary to claims in the literature (Lithwick 2009).

Therefore, it seems that no elliptical instability is found for . However, we recall that up to now we have only considered perturbations interior to the vortex core. But as we will see later, instabilities may also appear outside the vortex core.

3.4 Physical interpretation

It may be noted that the system (17) can be written as a single Hill's equation (see for example Waleffe 1990). In this case, one finds that the periodic excitation frequency is twice the vortex frequency , since the time-dependent wave numbers k_x and k_y always appear quadratically. In the limit for which the excitation is small, one would expect a series of instability bands. By analogy with Mathieu's equation, one should find such a band when the local epicyclic frequency  is an harmonic of the vortex frequency:

with . In the Keplerian case, one would therefore expect instability bands arising from for , , ... However, according to Fig. 3, the first unstable band (n=1) doesn't appear. This peculiar property can be understood writing the epicyclic frequency as

where is the absolute vorticity of the background vortex. According to this expression, the n=1 resonance condition is equivalent to the constraint . Interestingly, the absolute vorticity also appears explicitly in the linearized vertical vorticity equation:

where is the vertical vorticity of the perturbation. According to this equation, is a conserved quantity if the background absolute vorticity is zero. We conclude from this argument that the n=1 resonance condition cannot lead to an instability, since  could not be conserved in that case. Therefore, the first unstable band appears for n=2 ( ), which is consistent with our numerical result. As stated before, higher order bands (n>2) are also observed in the computation (Fig. 3), but they are much weaker, and one can't check the origin of these bands on the axis with enough precision to compare with Eq. (35).

4 Stratified case

4.1 Model and equations

Barranco & Marcus (2005) showed using anelastic numerical simulations that off-midplane vortices were able to survive for several hundred orbits whereas midplane vortices were destroyed rapidly. This suggests that stratification might be a way to stabilise vortices and suppress the elliptical instability. To study this problem in a simple configuration, we have considered the case of a 3D vortex centered at z=z₀, i.e. above the midplane. We then follow the evolution of perturbations with inside the vortex. In first approximation, the evolution of these perturbations can be described using the Boussinesq approximation which can be written

where is the potential temperature and N is the Brunt-Väisälä frequency. In the following, we will assume a stable vertical stratification, i.e. a real Brunt-Väisälä frequency. The vortex solution found previously is still valid in the Boussinesq approximation. Therefore, we can follow the same procedure as the one used in the non stratified case to get

These equations associated with the solution (16) lead to a stratified Floquet problem. As previously, this problem can be solved by diagonalizing a evolution matrix M, following the same procedure as in the non stratified case. The numerical method used to solve these equations is the same as in the previous section, except we have added an extra dimension in the solver to handle the evolution of the potential temperature.

4.2 Results

We plot the growth rate of the instability as a function of and  for a moderately stratified vortex N=S in Fig. 5. We note that the influence of the stratification is strong in the domain in which the vortex is weak (), and significant differences are observed when comparing with Fig. 3. First, we remark that the stable region observed in the non stratified case ( ) is not present in the N=S stratified case. This result is due to the presence of high frequency buoyancy modes which can match the turnover frequency when the inertial modes cannot. Moreover, the high- band has been replaced by a series of weaker bands. These bands can be understood as reminiscences of the resonance condition (35), since in the limit the effect of stratification disappears.

When we compute the maximum growth rate as a function of for various stratification frequencies (Fig. 6), we note that the growth rate at large decreases significantly with stronger stratification. In this sense, the stratification tends to weaken the elliptical instability, but does not suppress it.

5 Non-linear simulations

5.1 Resolving the elliptical instability

As pointed out in the previous section, the elliptical instability is always present for highly elongated vortices (), although it is quite weak. We have also found that no elliptical instability was observed for when stratification was omitted. To check these results, one has to perform non-linear simulations in which this instability is resolved. A condition to get this instability in a numerical simulation may be found assuming we have a 2D vortex embedded in a disc with a height H, the vortex core being larger than H. In the limit , we find the elliptical instability for (Fig. 3). According to (16), unstable modes will have in this case and . Moreover, since 3D perturbations have to fit vertically inside the disc. If we assume we need n points to resolve one wavelength (n>2), then the longest wavelength modes unstable to the elliptical instability will be resolved if the x resolution is n points per scale height and the y resolution is points per scale height. This last condition is not often satisfied in 3D global simulations, and may explain why long-lived and stable vortices are sometimes observed. Note that the value n may be quite high when using low order finite difference methods (one might expect ) since one wants the (numerical) dissipation rate for these wavelengths to be smaller than the growth rate, which is itself small compared to the shear frequency. Therefore, very accurate (e.g. spectral) numerical methods are preferable to study this kind of weak instability.

To check for resolution artefacts, the non stratified incompressible simulations and the compressible simulations carried out with NIRVANA below have been computed with two different resolutions (the ``high'' resolution being twice the ``low'' resolution in each direction). No significant difference was found between high resolution and low resolution runs. In particular, the same growth rates were obtained, with the same localisation properties for the instabilities. However, in low resolution runs carried out with NIRVANA that were seeded with truncation errors, the unstable scales were initially a factor 2 larger, but with the growth rates and later nonlinear behaviour being very similar. This can be understood as a consequence of the scale free property of the linear instability that holds once scales are sufficiently small but still larger than the dissipation scale, here expected to be the grid scale. Noting that conclusions would be unaltered if the low resolution runs were adopted, only high resolution results will be presented in this section.

5.2 Numerical solution for an isolated vortex

To check the existence of the elliptical instability in accretion disc vortices, one needs a solution for an isolated vortex. In Sect. 2, we have derived an exact solution for a vortex in an infinite box. Numerically however, one has to work in a finite size domain, using essentially the shearing-sheet boundary conditions (Hawley et al. 1995). As pointed out previously, the infinite solution has a slowly decaying tail which can lead to artefacts when one starts with this solution in a finite-size shearing box. When one tries this procedure, the 2D solution rapidly develops instabilities near the boundaries, especially along the x axis, leading to large amplitude oscillations. To prevent this, we have chosen to solve (8) using fast Fourier transforms (FFT), which are a natural choice for our boundary conditions. A solution obtained by this method in a box is shown in Fig. 7. The solution found by this procedure is very similar to the infinite solution, except near where an X point is observed. This difference probably explains the unsteadiness observed when one uses the infinite solution numerically. We find that periodic solutions derived using the FFT procedure are better suited for numerical simulations. Although these solutions are not exactly steady (mainly because of time varying boundary conditions), they keep the same shape and the same amplitude for several turnover times when one uses high precision numerical methods, which is enough for our purpose.

As pointed out previously, this vortex model is not a proper infinite elliptical flow, and the approximation used in our linear analysis may not be valid. Nevertheless, the core of these elliptical vortices is made of elliptical streamlines. Since our linear analysis is essentially a local analysis of stability on individual streamlines, one expects our results to apply to the core of such vortices. Note however that other parametric instabilities may appear outside of the vortex core because closed streamlines still exist in this region (see Fig. 7).

We have computed the evolution of several isolated elliptical vortices in a shearing box with dimensions using an incompressible spectral method (see Lesur & Longaretti 2005, for a complete description of the numerical method). We have used points and an eighth order hyperviscosity instead of the classical viscosity to prevent the dissipation of the weaker unstable modes. The vortices are initialised in the center of the box setting a vorticity patch with dimension ( ). We always set , L_y=4 and L_z=2, and L_x being adjusted according to the vortex aspect-ratio . For , we set L_x=10 whereas for and we set L_x=20. For each run, we initialise the velocity field solving (8) in Fourier space, and we add a small (10^-7) 3D white noise to the 2D solution.

5.3 Non stratified case

To follow the evolution of 3d motions, we have computed , where denotes a volume average procedure. Time histories of this quantity are given in Fig. 8 for , and . As expected from our linear analysis, we find an exponential growth for both and . Since the growth rate for is very fast, the instability saturates at and the vortex structure is destroyed, as shown in Fig. 9. In the exponential regime, we find a growth rate for and  for , in agreement with theoretical values ( ; ). We also find an instability for with which was not expected in our linear analysis ( elliptical streamlines are stable according to Fig. 4).

To check the localisation of the instability, we present (x,y) slices of the vertical velocity for , and in Fig. 10. We find that unstable modes are localised inside the vortex for and , as expected from our linear analysis. However, for , the instability appears outside of the vortex core, in a region in which one find closed non-elliptical streamlines (see Fig. 1). In this particular case, no vertical motions are found in the vortex core, which is consistent with our analysis.

The existence of this outer instability can be understood quite simply following the physical description used for the elliptical instability (Sect. 3.4). We know that these instabilities arise when the epicyclic frequency is an harmonic of a closed streamline frequency multiplied by (the streamline frequency being the inverse of the time needed by a fluid particle to cover the whole streamline and get back to its starting point). Moreover, as one moves away from the vortex core, the streamline frequency tends to 0 since the flow tends towards a pure shear flow, whereas the epicyclic frequency tends to the Keplerian frequency. Therefore, one will always be able to find a resonance outside of the vortex core, and consequently a parametric instability.

6 Effects of compressibility

We have also investigated the effects of compressibility. To do this we performed three dimensional hydrodynamic simulations using NIRVANA (Ziegler & Yorke 1997). NIRVANA has been used frequently in the past to study various problems involving MHD turbulence in the shearing box (Fromang & Papaloizou 2006; Papaloizou et al. 2004). Because of its diffusive character, it is best suited to vortices with small  that are not too elongated. For this reason we limit consideration to We consider isothermal shearing boxes with vertical gravity. In units such that the full radial width of the elliptical vorticity patch is unity, the box dimensions were (L_x= 3.46, L_y= 1.73, L_z= 2.31). The numbers of grid points used were (N_x= 196, N_y= 98, N_z= 128). A vertical gravitational acceleration was applied over 80% of the vertical computational domain, being set to zero in two domains extending from the vertical boundaries, each being of an extent equal to 10% of the whole vertical domain. The boundary conditions of periodicity in shearing coordinates could then be applied. We comment that the use of periodic boundary conditions in the vertical direction, with small buffer zones near the boundaries that allow these conditions to be applied consistently, has been found to be effective in vertically stratified shearing box simulations of MHD turbulence (e.g. Stone et al. 1996). This approach allows for regular treatment of the boundary without significant artefacts being produced in the mid-plane regions of the flow, as would be expected for local instabilities of the type considered here. To test this as well as the effects of expanding the computational domain in general, we have performed a test simulation with the compuational domain doubled in size in each direction for the case with the largest value of considered below. The dimensions of the elliptical vorticity patch were not changed. The growth rate and nonlinear development of the instability were indeed essentially the same in the mid-plane regions, which because of higher inertia carried most of the energy, as in the smaller domain simulations presented below. The instability in the larger domain case readily spreads to the more extensive regions outside the original vorticity patch.

We considered three values of the isothermal sound speed . Adopting the shearing time as the unit of time, these were such that were 1.3, 0.65, and 0.13. Here the velocity being unity in our dimensionless units is the maximum velocity difference across the vorticity patch due to the background shear.

The incompressible two dimensional flow field calculated in Sect. 5.2 was applied on horizontal planes. As above this was calculated by solving the vorticity equation with the imposition of periodic boundary conditions. To make this possible, one must subtract out the mean vorticity in the patch so that the net surface integral over the (x,y) plane is zero, which results in some unsteadiness as described above. However, in the compressible case, additional unsteadiness occurs because an associated unbalanced density change is generated. This was calculated and included in the set up by noting that for a strictly steady state horizontal flow of the type considered here, , with being the gravitational potential, should be a function only of the stream function . This is readily determined from the imposed steady flow from the condition

where is taken to be the same function of  as for the infinite vortex patch. The latter is of course an approximation as we modified the infinite patch by imposing periodic boundary conditions. Nonetheless its use to calculate the initial density on horizontal plane resulted in initial states that became increasingly more steady as increased.

In each of the simulations the two dimensional vortex structure was found to remain for about 40 shearing times before undergoing a linear instability which began with a very short wavelength in the vertical direction. The root mean square (weighted with density) vertical velocity is plotted in Fig. 11. Initial values increase with because the initially imposed flow is further from being stationary.

Streamlines in the horizontal midplane are shown for at t= 37, and at t=47 in Fig. 12. Although linear instability has begun, the initially imposed flow is found to be well preserved for the case with which most closely resembles the incompressible case. In this case the maximum growth rate in the linear phase of 0.75S is in good agreement with the linear prediction obtained from Eq. (23).

However, when the flow shows significantly greater time variability. The vortex is squeezed in the radial direction and trailing features begin to form outside the initial vortex core which turn into density waves which, however, are not well represented in our small computational domain. Global vortical structures survive until they are eventually destroyed by the instability Snapshots of the vertical velocity field in the horizontal midplane are shown after the onset of linear instability in Fig. 13 for the three simulations. The 3D perturbations are found to be localized inside the original elliptical boundary of the vorticity patch just as in the incompressible case.

6.1 Stratified case

To study the vertically stratified case, we have implemented the Boussinesq Eqs. (38) in our spectral code. Compared to (38), we have added an eighth order hyper diffusivity to the velocity and potential temperature equations. We followed the same procedure as in the non stratified case to set up the vortex structure. We have run the simulations with a moderate stratification (N=S), for which one expects a departure from the non stratified case when  (Fig. 5). The time history of , for various vortex aspect-ratios, is presented in Fig. 14.

For we find and for . No instability is observed in our simulations for . These measured growth rates are always below the expected growth rate from our linear analysis ( , and ) which might be due to excessive numerical dissipation in the stratified case, despite the high resolution and the hyper viscosity prescription used. For both and , the instability appears to be localised inside the vortex core, as expected. We can conclude from our results that the elliptical instability can be detected in stratified simulations, but resolution and numerical dissipation must be carefully controlled to get accurate results.

7 Discussion

We have described a 3D instability appearing in elliptical vortices (Kida vortices) embedded in accretion discs, known as the elliptical instability. We have shown analytically, for the case of no stratification that this elliptical instability is always found in vortex cores, except for a narrow range of vortex aspect-ratio ( ). When the vortex is weak (elongated vortex), the instability involves small radial wavelengths compared to azimuthal and vertical ones, the ratio being of the order of . Moreover, the inclusion of a stable stratification tends to weaken the instability for large  but it does not suppress it.

Numerical simulations have essentially confirmed our linear analysis. We have found the predicted growth rate in the case where the instability was expected inside the core. Furthermore, instabilities outside of the vortex core were observed when the core was linearly stable, in agreement with our physical interpretation. The inclusion of compressibility in the simulations didn't change the behaviour of the instability, and analytical results were recovered. We also studied numerically the effect of vertical Boussinesq stratification. Again, analytical results were recovered, except for weak vortices () for which stability was observed, probably because of numerical diffusion.

These results tend to indicate that elliptical instabilities, or more generally parametric instabilities, are always present in accretion disc vortices. However, they often involve small radial wavelengths (compared to the disc thickness) and small growth rates (compared to the orbital time scale). Put together, these properties make elliptical instabilities hard to capture numerically because of numerical diffusion at the grid scale. For example, highly accurate spectral methods with high-order hyper-diffusivity were required to marginally resolve the instability in a  vortex.

In this work, we have not addressed the saturation of this instability. This is however rather intentional, as we think that the non-linear outcome should be studied together with the production mechanism for these vortices. In fact, the evolution of these structures will be dictated ultimately by the competition between the production mechanism and the saturated 3D instability. In this context, studying the saturation of this instability on its own will be of little interest, since it will be modified by the vortex production process. Note too that the time-scale needed by the saturated instability to destroy the vortex structure is unknown. In particular, there is no reason to think that this time-scale is related to the growth rate in the linear stage. For example, a weakly unstable vortex could be destroyed in one orbit once saturation is reached. Therefore, no firm conclusion can be drawn from this work concerning the evolution of the vortex structure itself. However, in the absence of a sustaining production process, eventually one would expect the destruction of the vortex structure when the perturbation amplitude reaches the background flow amplitude as was indeed observed for our simulation.

This work has also several limitations. In particular, we have assumed an isolated 2D vortex as a starting point. However, one might want to study 3D vortices including a complicated vertical structure, with e.g. a mean vertical flow in the core (similar to cyclones on Earth for example). This kind of structure was suggested by Barranco & Marcus (2005) for off-midplane vortices, and there is no doubt that stability properties will be modified in this case. Moreover, we have mostly restricted our study to elliptical streamlines, since these are tractable analytically. Although some extensions to more general cases were shown, we have not provided any formal proof of instability in the most general case. However, the instability mechanism exhibited in the elliptical case appears to be quite general, and we think it's unlikely that a steady vortex could be stable everywhere to 3D perturbations.

Finally, these results point out the necessity of carrying high-resolution 3D simulations of the mechanisms suggested in the literature for vortex formation. This includes for example Rossby wave instabilities (Li et al. 2001), baroclinic instabilities (Klahr & Bodenheimer 2003; Petersen et al. 2007) or gaps due to giant planets (de Val-Borro et al. 2007). In any case, the presence of 3D parametric instabilities will certainly modify the global outcome of these processes, with probably new physical effects.

Acknowledgements

The authors thank Jeremy Goodman, Gordon Ogilvie, Orkan Umurhan and François Rincon for useful discussions. The simulations presented in this paper were performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England. GL acknowledges support by STFC.

References

• Balbus, S. A. 2003, ARA&A, 41, 555 [NASA ADS] [CrossRef] (In the text)
• Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] (In the text)
• Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] (In the text)
• Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [NASA ADS] [CrossRef] (In the text)
• Bayly, B. J. 1986, Phys. Rev. Lett., 57, 2160 [NASA ADS] [CrossRef] (In the text)
• Bodo, G., Tevzadze, A., Chagelishvili, G., et al. 2007, A&A, 475, 51 [NASA ADS] [CrossRef] [EDP Sciences]
• Craik, A. D. D. 1989, J. Fluid Mech., 198, 275 [NASA ADS] [CrossRef] (In the text)
• Craik, A. D. D., & Criminale, W. O. 1986, Proc. R. Soc. Lond. A, 406, 13 [NASA ADS] [CrossRef] (In the text)
• de Val-Borro, M., Artymowicz, P., D'Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Drazin, P., & Reid, W. 1981, Hydrodynamic stability (Cambridge Univ. Press) (In the text)
• Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751 [NASA ADS] [CrossRef] [EDP Sciences]
• Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
• Godon, P., & Livio, M. 1999, ApJ, 523, 350 [NASA ADS] [CrossRef]
• Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] (In the text)
• Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Johnson, B. M., & Gammie, C. F. 2005a, ApJ, 626, 978 [NASA ADS] [CrossRef]
• Johnson, B. M., & Gammie, C. F. 2005b, ApJ, 635, 149 [NASA ADS] [CrossRef]
• Johnson, B. M., & Gammie, C. F. 2006, ApJ, 636, 63 [NASA ADS] [CrossRef]
• Kelvin, L. 1880, Phil. Mag., 10, 155 (In the text)
• Kerswell, R. R. 2002, Ann. Rev. Fluid Mech., 34, 83 [NASA ADS] [CrossRef] (In the text)
• Kida, S. 1981, Phys. Soc. Japan J., 50, 3517 [NASA ADS] [CrossRef] (In the text)
• Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef]
• Lesur, G., & Longaretti, P.-Y. 2005, A&A, 444, 25 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874 [NASA ADS] [CrossRef] (In the text)
• Lithwick, Y. 2009, ApJ, 693, 85 [NASA ADS] [CrossRef] (In the text)
• Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] (In the text)
• Onsager, L. 1949, Nuov. Cim. Supp., 6, 279 [CrossRef] (In the text)
• Papaloizou, J. C. B., Nelson, R. P., & Snellgrove, M. D. 2004, MNRAS, 350, 829 [NASA ADS] [CrossRef]
• Petersen, M. R., Julien, K., & Stewart, G. R. 2007, ApJ, 658, 1236 [NASA ADS] [CrossRef]
• Pierrehumbert, R. T. 1986, Phys. Rev. Lett., 57, 2157 [NASA ADS] [CrossRef] (In the text)
• Regev, O., & Umurhan, O. M. 2008, A&A, 481, 21 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Shen, Y., Stone, J. M., & Gardiner, T. A. 2006, ApJ, 653, 513 [NASA ADS] [CrossRef] (In the text)
• Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656 [NASA ADS] [CrossRef] (In the text)
• Umurhan, O. M., & Regev, O. 2004, A&A, 427, 855 [NASA ADS] [CrossRef] [EDP Sciences]
• von Weizsäcker, C. F. 1944, Z. Astrophys., 22, 319 (In the text)
• Waleffe, F. 1990, Phys. Fluids A, 2, 76 [NASA ADS] [CrossRef] (In the text)
• Ziegler, U., & Yorke, H. W. 1997, Comp. Phys. Comm., 101, 54 [NASA ADS] [CrossRef] (In the text)

All Figures

---

Copyright ESO 2009