Issue |
A&A
Volume 498, Number 1, April IV 2009
|
|
---|---|---|
Page(s) | 1 - 12 | |
Section | Astrophysical processes | |
DOI | https://doi.org/10.1051/0004-6361/200811577 | |
Published online | 11 March 2009 |
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=R0, 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




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
![]() |
(3) |
where we have defined the vortex aspect-ratio



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
ui0 = S Aij xj, | (6) |
defining
![]() |
(7) |
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

![]() |
(9) | ||
![]() |
(10) |
with







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.
![]() |
Figure 1:
Vortex solution streamlines in a sheared flow with |
Open with DEXTER |
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

![]() |
(15) |
which leads to:
where k0,







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



![]() |
(18) |
where f is periodic with period





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, vz=0 and (17) is reduced to
![]() |
(21) | ||
![]() |
(22) |
This system describes horizontal epicyclic oscillations with frequency
having defined the rotation number as



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
(ux0 , uy0). The linearized equations of motion are
where

![]() |
(26) |
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



Replacing


![]() |
(29) |
there will be an instability. In fact for the problem on hand this is precisely equivalent to the above condition that

![]() |
(30) |
then both



or setting qy=-vy
Thus if




or equivalently
Thus it follows that vxqy grows faster than exponentially with growth rate

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.
![]() |
Figure 2: Growth rate (in shear time) of the 3D instability deduced numerically from the Floquet problem (17) in the non rotating case. This result is identical to the original elliptical instability described by Bayly (1986) and allow us to test our numerical approach. The colour scale represent the logarithm of the growth rate. |
Open with DEXTER |
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.
![]() |
Figure 3:
Growth rate of the 3D instability deduced numerically from the Floquet problem (17) in the Keplerian case for large |
Open with DEXTER |
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 kx and ky 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






![]() |
(36) |
where


![]() |
(37) |
where






![]() |
Figure 4:
Maximum growth rate in the Keplerian case (in shear time). We remark the low- |
Open with DEXTER |
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=z0, 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

These equations associated with the solution (16) lead to a stratified Floquet problem. As previously, this problem can be solved by diagonalizing a

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.
![]() |
Figure 5:
Growth rate of the 3D instability deduced numerically from the Floquet problem (41) in the Keplerian case with N/S=1.0. One observe once again the vertical modes for
|
Open with DEXTER |
![]() |
Figure 6:
Maximum growth rate in the Keplerian case (in shear time) for various stratification frequencies. As the stratification gets stronger, the growth rate decreases for weak vortices ( |
Open with DEXTER |
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
,
Ly=4 and Lz=2,
and Lx being adjusted according to the vortex aspect-ratio
.
For
,
we set Lx=10 whereas for
and
we set Lx=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.
![]() |
Figure 7:
Vortex streamlines in a sheared flow with |
Open with DEXTER |
![]() |
Figure 8:
Time evolution of
|
Open with DEXTER |
![]() |
Figure 9:
3D snapshots showing the evolution of a |
Open with DEXTER |
![]() |
Figure 10:
2D snapshots of the vertical velocity field in the (x,y) plane for |
Open with DEXTER |
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
(Lx= 3.46, Ly= 1.73, Lz= 2.31). The numbers of grid points used were (Nx= 196, Ny= 98, Nz= 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.
![]() |
Figure 11:
The root mean square vertical velocity obtained from the
compressible simulations is plotted as a function of time
for |
Open with DEXTER |
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.
![]() |
Figure 12:
Streamlines in the (x,y) plane for |
Open with DEXTER |
![]() |
Figure 13:
2D snapshots of the vertical velocity field in the (x,y) plane
for |
Open with DEXTER |
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
![]() |
(43) |
where



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.
![]() |
Figure 14:
Time evolution of
|
Open with DEXTER |
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
![]() |
Figure 1:
Vortex solution streamlines in a sheared flow with |
Open with DEXTER | |
In the text |
![]() |
Figure 2: Growth rate (in shear time) of the 3D instability deduced numerically from the Floquet problem (17) in the non rotating case. This result is identical to the original elliptical instability described by Bayly (1986) and allow us to test our numerical approach. The colour scale represent the logarithm of the growth rate. |
Open with DEXTER | |
In the text |
![]() |
Figure 3:
Growth rate of the 3D instability deduced numerically from the Floquet problem (17) in the Keplerian case for large |
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Maximum growth rate in the Keplerian case (in shear time). We remark the low- |
Open with DEXTER | |
In the text |
![]() |
Figure 5:
Growth rate of the 3D instability deduced numerically from the Floquet problem (41) in the Keplerian case with N/S=1.0. One observe once again the vertical modes for
|
Open with DEXTER | |
In the text |
![]() |
Figure 6:
Maximum growth rate in the Keplerian case (in shear time) for various stratification frequencies. As the stratification gets stronger, the growth rate decreases for weak vortices ( |
Open with DEXTER | |
In the text |
![]() |
Figure 7:
Vortex streamlines in a sheared flow with |
Open with DEXTER | |
In the text |
![]() |
Figure 8:
Time evolution of
|
Open with DEXTER | |
In the text |
![]() |
Figure 9:
3D snapshots showing the evolution of a |
Open with DEXTER | |
In the text |
![]() |
Figure 10:
2D snapshots of the vertical velocity field in the (x,y) plane for |
Open with DEXTER | |
In the text |
![]() |
Figure 11:
The root mean square vertical velocity obtained from the
compressible simulations is plotted as a function of time
for |
Open with DEXTER | |
In the text |
![]() |
Figure 12:
Streamlines in the (x,y) plane for |
Open with DEXTER | |
In the text |
![]() |
Figure 13:
2D snapshots of the vertical velocity field in the (x,y) plane
for |
Open with DEXTER | |
In the text |
![]() |
Figure 14:
Time evolution of
|
Open with DEXTER | |
In the text |
Copyright ESO 2009
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.