The formation and stability of vortices in protoplanetary disks

We review vortex formation through instabilities associated with structural features in protoplanetary disks. We then describe the form of the generic parametric or elliptical instability associated with and, localized on, non circular streamlines. The effect of this type of instability on the potential accumulation of solid particles in the centre of such vortices has yet to be determined.


Introduction
Planets are believed to form in circumstellar disks of gas and dust circulating around around young stars.Planet formation requires that solid material accumulates from sub-micron sized interstellar grains producing planetesimals which in turn amalgamate to form planets. Large-scale, long-lived vortices may play an important role in this process, since gas drag can result in them collecting solid particles in their centres (Barge & Sommeria 1995) potentially speeding up the accumulation process (Lyra et al. 2009).
Such vortices require a driving process to maintain them against known potential instabilities and viscous decay.Possible processes involve instabilities associated with features in the protoplanetary disk such as rings and gaps and also baroclinic effects.We discuss these here as well as the ideas behind possible parametric or elliptical instabilities associated with non circular streamlines.

Instability of rings and edges
Rings and edges may form in protoplanetary disks as a result of a variety of effects.Structures of this type are expected at the boundary between so called live and dead zones.The inner live zone has magnetohydrodynamic turbulence and large effective viscosity causing it to be cleared quickly and thus be of low density, while matter piles up in the relatively inert (apart from active surface layers) dead zone that has relatively small effective viscosity.The transition region is expected to be about a scale height in width and to form an edge like structure with a mound (for recent simulations of vortex production in this context see eg.Meheut et al. 2012, Lyra & Mac Low 2012).
The other situation that is likely to produce edge structures is when a planet produces a gap in protoplanetary disk as a result of the gravitational perturbations it exerts.This effect was studied by de ratio 0.05.All these cases show instability to vortex formation on account of there being a maximum associated with a localised mound of inverse vortensity, or surface density per unit vorticity.
Nonlinear 2D simulations indicate that, in the larger planet case, small vortices merge resulting in the eventual formation of a single large vortex.Lin & Papaloizou (2011) studied the effect of selfgravity, showing that this reduced the scale of the original instability and hindered vortex merging, until being replaced by a different kind of edge instability when the effect of self-gravity becomes large enough (see Fig. 1).

Unstable normal modes in rings and edges
When considering disk stability one looks at linear normal modes with angular variation through a factor exp(imφ), where m is the azimuthal mode number and φ the azimuthal angle.It turns out that potentially unstable normal modes with low azimuthal mode number m r/H, the inverse aspect ratio (ratio of radius to scale height), are effectively confined to the edge region, whereas those associated with higher m radiate density waves.In the former case, the governing equation for a 2D isothermal disk can be taken to be Here W = Σ /Σ is the fractional perturbation to the surface density, Σ, c s is the isothermal speed and η = (rΣ) −1 d(r 2 Ω)/dr is the specific vorticity or vortensity, with Ω being the local angular velocity.
The quantity σ = (Ω − ω p ), where mω p is the mode frequency.When this is real there is a corotation resonance where σ vanishes.At marginal stability, for a given m, corotation coincides with the minimum in vortensity so that (1) is non singular (Papaloizou & Pringle 1985).For smaller m instability is expected. 01002-p.2

Similarity to Kelvin-Helmholtz instability
The above analysis of this instability is similar to that for the Kelvin-Helmholtz instability of an incompressible one dimensional flow for which the velocity in the Cartesian x direction is u(z).The stream function perturbation, f, satisfies the inviscid Orr-Sommerfeld equation (Drazin & Reid 1981): where k is the magnitude of the wavenumber in the direction of the unperturbed velocity and kc is the mode frequency.The parallels between (2) and ( 1) are clear.The role of the vortensity gradient is replaced by the vorticity gradient in the former.The instability can then be seen as the result of interaction or mixing of material on both sides of corotation with the same vorticity/vortensity.This process naturally taps the free energy available in the shear flow.
Early work on disk problems assumed a barotropic equation of state and was mostly 2D though limited analysis indicated that results would be preserved in 3D.Li et al. (2000) considered the 2D case allowing for a non barotropic equation of state calling the instability a 'Rossby wave instability'.The essential 2D nature of the instability has been confirmed by recent simulations of Lin (2012)a and Meheut et al. (2012) in the barotropic case.However, a linear analysis by Lin (2012)b indicates 3D effects may be more significant in the non barotropic case.

Subcritical baroclinic instability
An alternative method for amplifying vortices has been indicated by Petersen et al. (2007) and Lesur & Papaloizou (2010).The conditions required are a disk unstable to the radial Schwarzschild (convection) criterion.together with non zero thermal diffusion or cooling.In addition, streamlines that are closed around the vortex centre are required from the outset, implying that a finite amplitude initial perturbation corresponding to a subcritical instability is required.The amplification works through the action of heating and cooling on material circulating around the vortex resulting in buoyancy forces that tend to accelerate the circulation or amplify the vortex (see Lesur & Papaloizou 2010).This effect is of course opposed by viscosity.
In the compressible shearing boxes considered by Lesur & Papaloizou (2010), vortices produce density waves associated with weak outward angular momentum transport .Values of the standard α parameter up to several ×10 −3 have been seen.However, shearing box simulations tend to evolve towards a state with one vortex per box resulting in the whole of space being filled with a web of interacting vortices (see Fig. 2) that is probably unrealistic.Proper global simulations are thus required to address the angular momentum transport issue.Having discussed the way vortices may be set up through the operation of various (mainly 2D) instabilities, we now discuss the issue of their stability.First we consider the setup of steady configurations and then go on to discuss methods for analysing stability.

Structure and stability of vortex configurations
The governing equation of motion in a frame rotating with angular velocity, Instabilities and Structures in Proto-Planetary Disks

EPJ Web of Conferences
Here the velocity is v, the density is ρ and the system is subject to a force per unit mass F which has components due to pressure, P, and gravity, with the fixed gravitational potential being Ψ.Self-gravity is neglected.For steady state configurations equation (3) reduces to For the domain of consideration, we adopt a local shearing box with origin centered on a point of interest and rotating with its Keplerian angular velocity, thus fixing Ω P (Goldreich & Lynden-Bell 1965).A Cartesian coordinate system is adopted with the y axis in the direction of shear and the z axis normal to the disk mid-plane.For a general vector a we adopt the representations a ≡ (a x , a y , a z ) ≡ (a 1 , a 2 , a 3 ).For an undisturbed background Keplerian flow, v = (0, −3Ω P x/2, 0).It is possible to consider very small particles and gas as a variable density incompressible flow.Because the stopping time for drift relative to the gas is zero, the particle distribution is frozen into the gas so that we have Dρ/Dt = 0, as well as ∇ • v = 0.For a steady 2D flow with stream function ψ, this implies that ρ = ρ(ψ) is a function of ψ alone.This would have to be a snapshot of a configuration with solid material slowly migrating towards the vortex centre.Such locations are potential sites of planet formation and accordingly the focus of much attention.

The Kida solution
A well known steady state 2D solution with no z dependence is a vortex, the core of which is an elliptical patch of constant vorticity ω t = −3Ω P /2+ω v , where −3Ω P /2 if the background flow vorticity as seen in the rotating frame and ω v is the vorticity of the vortex (Kida 1981) .Outside the core, the vorticity is that of the background flow.The aspect ratio, or ratio of major to minor axis is related to The components of v for a vortex centered on the origin are (3Ω P y/(2χ(χ − 1)), −3Ω P χx/(2(χ − 1)), 0) and accordingly are linear functions of the coordinates.Note that the vorticity of the background is positive in the inertial frame so that ω v being negative corresponds to an anticyclonic vortex.This solution is an important reference case and importantly it can be generalized to apply to arbitrary vorticity distributions and also to incorporate a variable density, but can then only be found numerically (see Railton & Papaloizou this volume).However, the latter step is essential for a selfconsistent stability analysis of configurations of interest.

The stability of general incompressible vortical flows allowing for density gradients
We now consider the stability of steady state flows of the type introduced above.A convenient formalism is the Lagrangian approach developed by Lynden-Bell & Ostriker (1967).Following these authors we introduce the Lagrangian variation ∆ such that the change to a state variable Q as seen following a fluid element is ∆Q.Thus the Lagrangian displacement is given by ∆r = ξ.We also have The Eulerian variation, Q is the change in Q as seen in a fixed coordinate system.It is given by Q = ∆Q − ξ • ∇Q.We now take the Lagrangian variation of the equation of motion (3), noting that the operators D/Dt and ∆ commute ( Lynden-Bell & Ostriker1967).Thus we obtain 01002-p.4 with The Lagrangian variation in the density is zero, so that we have Equations ( 6) and ( 8) together with the incompressibility condition ∇ • ξ = 0, lead to effectively four equations for P and the three components of ξ.

Local Analysis
We consider eventually perturbations that are localized on streamlines, that is they have short wavelengths in the directions perpendicular to the unperturbed velocity, but can have a long wavelength in the direction of the unperturbed velocity.The latter is a natural outcome of shearing motions.To do this we begin by assuming that any perturbation quantity takes the form Here we adopt a WKBJ ansatz with the phase function S left arbitrary for the time being and the constant λ taken to be a large parameter.For a discussion of the approach followed here in avariety of contexts, see Lifschitz & Hameiri (1991), Sipp & Jacquin (2000) and Papaloizou (2005).The effective wavenumber k = λ∇S (10) then has a large magnitude.The amplitude factor ∆Q 0 is the WKBJ envelope.Because of the rapid variation of the complex phase λS , one can perform a WKBJ analysis, such that the state variables, 01002-p.5 Instabilities and Structures in Proto-Planetary Disks EPJ Web of Conferences ∆Q 0 are expanded in inverse powers of λ.The lowest order term in ξ is constant, while the lowest order term in P is ∝ λ −1 .To lowest order (6) gives When working to the next order, on account of the rapidly varying phase, only the variation of this needs to be considered when taking spatial derivatives, apart from when considering expressions involving the operator, D/Dt ≡ ∂/∂t + v • ∇, because this annihilates S .Accordingly the contribution v • ∇(∆Q 0 ) must be retained.Noting the above, we can substitute perturbations of the form ( 9) into the governing equations and remove the factor exp(iλS ), thus obtaining equations for the lowest order contribution to the quantities ∆Q 0 alone.For ease of notation we drop the subscript 0 from now on.Following this procedure (6) gives where we recall that the lowest order term for P ∝ λ −1 .In addition to this the incompressibility condition gives Using (10) and (11) we can also find an equation for the evolution of k in the form Equations ( 12), ( 13) and ( 14) give a complete system for the evolution of ξ and k, after the elimination of P , as an initial value problem.Because the evolution consists of advection of data along streamlines, it is possible to consider disturbances localized on individual streamlines (Papaloizou 2005).Localization amplitudes are unaffected by the evolution considered.In general one could start with an arbitrary initial S and then k would depend on time.

Solutions for a time independent wavenumber
A relatively simple class of solutions for k can be obtained by setting S to be independent of time and a function of quantities conserved on unperturbed streamlines so that v • ∇S = 0 (Papaloizou 2005).Then k is fixed for all time and we only have to solve for ξ.For the simple case of a two dimensional vortex with initial state independent of z, we may take S = g(ψ) + k z z/λ.Here ψ is the unperturbed stream function, g is an arbitrary function and k z is the constant vertical wavenumber.Adopting this form of S and the corresponding k, assuming k z 0, which is appropriate to the physically realistic case where perturbations are localized in z, we may use the vertical component of (12) together with (13) to eliminate P and ξ z and thus obtain a pair of equation for (ξ x , ξ y ) ≡ (ξ 1 , ξ 2 ).These can be written in the form where 01002-p.6 Note that as only the horizontal components are considered, the summation is for j = 1, 2 only.In addition we readily find from ( 14) that We remark that although there is an arbitrary function g in the definition of k used in this section, because the derivatives in (15) correspond to advecting around streamlines and dg/dψ is constant on streamlines it effectively behaves as a multiplicative constant merely scaling the magnitude of the wavenumber.

Generic instability
The analyses in the previous sections reduces the stability problem to solving an initial value problem of integrating a set of simultaneous ordinary differential equations around streamlines.The independent variable measures the location on a streamline.As the unperturbed motion on streamline is periodic, these equations have periodic coefficients through their dependence on k and the gradient of ρ etc.Thus the Floquet theory may be applied (eg.Whittaker and Watson 1996).According to this, if some internal mode with a natural oscillation frequency dependent on k is described by (15) unstable bands of exponential growth are expected as k is varied to allow resonances of this frequency with the frequency of motion around the streamline.
For example if we assume v = (0, −3Ω P x/2, 0), and take k = (k x , 0, k z ) to be constant, corresponding to the usual shearing sheet, (15) has oscillation modes corresponding to inertia-gravity waves satisfying the dispersion relation , where the square of the buoyancy frequency frequency N 2 = ∇ρ • ∇P/ρ 2 .The background flow is for the most part only slightly modified for an elongated vortex so we might expect the analogues of such modes to be parametrically excited.That would also be expected in the situation where the straight background streamlines are perturbed to a new steady state, while maintaining periodicity in the y direction.Note that because we can choose the ratio k x /k z arbitrarily in these examples a continuous range of mode frequencies is accessible so that resonance conditions can always be satisfied.This makes parametric instability generic.

The stability analysis of a Kida vortex core
The viewpoint described in section 3.2.1 differs from that found in the stability analysis of the Kida vortex core given by Lesur & Papaloizou (2009), the results of which are illustrated in Fig. 2.However, equations ( 12) -( 14) of section 3.2 apply in this case.They are satisfied by the phase function λ S = k(t) • x where the wave vector k is given by k(t) = k 0 χ sin(θ) sin φ(t) , sin(θ) cos φ(t) , cos(θ) , where k 0 and θ are constants.(18) with φ(t) = 3Ω P /(2(χ − 1))(t − t 0 ), t 0 being constant.This is not of the form (9) which is specifically taylored to disturbances with short wavelengths in the cross stream direction and is used in section 3.2.1.Note that the latter form can be applied generally whereas (18) only works for a linear shear.However, both types of solution satisfy the equations governing local stability in latter case.This raises the issue as to whether representations are complete.Nonetheless it has been shown that if exponential growth is found for any solution of the local problem, instability is guaranteed (see eg. Lifschitz & Hameiri 1991, Papaloizou 2005).Thus the approach is more useful for showing instability rather than stability.Note that in spite of the restriction on the wavenumber form given by ( 18), Chang & Oishi (2010) used it in the problem of the stability of a vortex with a density gradient, for which the shear would not be expected to be self-consistently linear, and obtained a heavy core instability.

01002-p.7
Instabilities and Structures in Proto-Planetary Disks EPJ Web of Conferences

Discussion
The generation of vortices through instabilities at gap edges induced by orbiting planets, or at an interface between a live and dead zone, or through amplification through the subcritical baroclinic instability has been shown to be robust.The centres of such vortices have been proposed be the sites of concentration of solid particles that migrate there through the action of gas drag.Their survival and evolution however, depend on a balance between formation mechanisms and the action of parametric or elliptical instabilities.While the vortices may survive or be regenerated in some form (eg. Lesur & Papaloizou 2010), the effect on solid accumulation is less clear.In particular caution should be exerted with respect to the analysis of Chang & Oishi (2010) who assumed (18) for a vortex with a density gradient for which the shear is no longer expected to be linear in a self-consistent equilibrium (see Railton & Papaloizou in this volume).Note too that when finite stopping times are taken into account, instabilities dependent on density gradients may not operate if particles migrate through a wavelength within a putative growth time.The instabilities we described, being local, may work on small length scales that are difficult to resolve in 3D simulations.This underlines the need for a rigorous stability analysis for a self-consistent background state which we hope will be undertaken in the near future.

Figure 1 .
Figure 1.Vortices produced through the instability of the edge of a gap produced by a Saturn-mass planet in a disc with aspect ratio 0.05.The values of the relative surface density change from the background disk profile are indicated at the number of orbital periods indicated after initiation, taken from Lin & Papaloizou (2011).The units are such that the circular orbit is at radius r = 5.The minimum Toomre parameter in the background disc is from left to right 1.5, 2 and 8 respectively so that self-gravity decreases in importance in this progression (see Lin & Papaloizou 2011 for details).

Figure 2 .
Figure 2. Left: Trailing density waves excited by vortices amplified by the baroclinic instability formed by sticking 9 shearing boxes together consistently with the boundary conditions, taken from Lesur & Papaloizou (2010) which may be consulted for more details .The vortices occupy the low density (dark but not darkest) regions.Right: Growth rate of the instability of a Kida vortex with Keplerian background as a function of the aspect ratio, χ, and θ, with the ratio of vertical to horizontal wavenumber being cot θ, expressed as a multiple of 3Ω P /2, taken from Lesur & Papaloizou (2009).Note the three instability bands including a weak one for large χ.