On the local stability of vortices in differentially rotating discs
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 2 October 2018 (MN L A TEX style file v2.2)
On the local stability of vortices in di ff erentially rotating discs A. D. Railton (cid:63) and J. C. B. Papaloizou † Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Centre for Mathematical Sciences,Wilberforce Road, Cambridge, CB3 0WA, UK
ABSTRACT
In order to circumvent the loss of solid material through radial drift towards the central star,the trapping of dust inside persistent vortices in protoplanetary discs has often been suggestedas a process that can eventually lead to planetesimal formation. Although a few special caseshave been discussed, exhaustive studies of possible quasi-steady configurations available fordust-laden vortices and their stability have yet to be undertaken, thus their viability or oth-erwise as locations for the gravitational instability to take hold and seed planet formation isunclear. In this paper we generalise and extend the well known Kida solution to obtain a seriesof steady state solutions with varying vorticity and dust density distributions in their cores, inthe limit of perfectly coupled dust and gas. We then present a local stability analysis of theseconfigurations, considering perturbations localised on streamlines. Typical parametric insta-bilities found have growthrates of 0 . Ω P , where Ω P is the angular velocity at the centre ofthe vortex. Models with density excess can exhibit many narrow parametric instability bandswhile those with a concentrated vorticity source display internal shear which significantlya ff ects their stability. However, the existence of these parametric instabilities may not neces-sarily prevent the possibility of dust accumulation in vortices. Key words: planetary systems: formation — planetary systems: protoplanetary discs
Studies of planet formation have been ongoing since the for-mulation of the nebula hypothesis for the formation and earlyevolution of the Solar System (Swedenborg 1734; Kant 1755;Weizs¨acker 1944). Centrifugal forces for the most part balancethe gravitational attraction of the central star and a protoplanetary(PP) disc that forms along with it. The disc contains dust grains0 . − µ m in size which undergo coagulation (eg. Safronov 1969;Dominik & Tielens 1997) through the action of electrostatic ratherthan gravitational forces, the latter being expected to dominateduring the later stages of planet formation (eg. Safronov 1969;Lissauer 1993; Papaloizou & Terquem 2006).The notion of planetesimal formation, through the sticking togetherof dust grains in PP discs through two-body collisions, was first de-veloped by Chamberlin (eg. Chamberlin 1900). However, bodiesabove about a meter in size have very poor sticking qualities (Benz2000) so that collisions between them will result in fragmentationor bouncing rather than growth. In addition, particles in a typicalPP disc, with mid plane pressure decreasing monotonically withradius, experience a headwind in the azimuthal direction which (cid:63) E-mail: [email protected] † E-mail: [email protected] causes them to lose angular momentum and drift radially towardsthe central star. As a consequence, metre-sized bodies may spiralinto the star on timescales as short as a hundred years (eg. Weiden-schilling 1977; Papaloizou & Terquem 2006). This indicates thatplanetesimals must be formed within the rapid radial drift time ofthese bodies. These two di ffi culties for planetesimal formation con-stitute the metre-size barrier.However, in a disc for which the mid plane pressure does not de-crease monotonically with radius, aerodynamic e ff ects can concen-trate solids in allowing planetesimals to form. As a consequenceof the fact that particles tend to drift in the direction of the pres-sure gradient, Whipple (1972) showed that a pressure maximumlocated in an axisymmetric ring is a very e ff ective particle trap. Ifsu ffi cient concentration can occur, planetesimal formation can thenbe assisted by gravitational instability (Safronov 1969; Goldreich& Ward 1973). In the context of the above scenario, MRI simula-tions indicate the possibility of zonal flows that produce long-livedaxisymmetric pressure bumps (Johansen, Youdin & Klahr 2009;Fromang & Stone 2009).Isolated pressure maxima can exist in the centre of anticyclonic vor-tices. Barge & Sommeria (1995) showed that such vortices couldalso be natural localised particle traps, proposing that they could besites of planetesimal formation. A number of authors have shownthat anticyclonic vortices can form coherent and long-lived struc-tures in PP discs (Bracco et al. 1999; Chavanis 2000; Barranco c (cid:13) a r X i v : . [ a s t r o - ph . E P ] O c t Railton and Papaloizou & Marcus 2005). In addition there are diverse means of generat-ing these vortices in discs, namely through the Baroclinic Instabil-ity (Petersen, Julien & Stewart 2007; Lesur & Papaloizou 2010),the instability at the interface between MRI active and dead zones(Lovelace et al. 1999; M´eheut et al. 2010), and the edge instabilityassociated with gaps in the disc produced by existing planets (Lin& Papaloizou 2011).However, it is well known that such vortices are prone to the socalled elliptical instability, which can be regarded as a local para-metric instability associated with periodic motion on streamlines(see Lesur & Papaloizou 2009, and references therein). However,this has only been analysed in full detail for the special case ofa Kida vortex (see Kida 1981) with no dust present. These solu-tions apply to a local patch of the disc that can be represented usingthe well known shearing box formalism (Goldreich & Lynden-Bell1965). However, as discussed in this paper, a large variety of vortexconfigurations can be constructed with and without dust concen-trations. Their vorticity and density profiles have a degree of arbi-trariness when no frictional or di ff usive processes operate, althoughthey would be expected to be determined by the form of these whenthey do.Since the existence of instabilities in such vortices could be a threatto their survival or dust attracting capability, a comprehensive sta-bility analysis is desirable. However, up to now only velocity pro-files that give a constant period of circulation such as in the Kidavortex has been considered in theoretical developments. In this con-text Chang & Oishi (2010) considered the 2D (independent of thevertical direction) stability of a dust laden vortex assuming such aprofile but did not consider the issue of whether the profiles adoptedprovided either a steady state or matched onto a suitable disc back-ground. They assumed a separability that applies strictly to the uni-form density Kida solution but not more general cases. It is thepurpose of this paper to further consider the structure and stabilityof vortices in a protoplanetary disc background.For simplicity we consider vortices with radial length scale lessthan the disc scale height for which the dust stopping time is veryshort. In this case an incompressible fluid model with frozen-indensity distribution can be adopted. We consider vortices allow-ing for both non-uniform vorticity and density distributions in theircores. We consider local stability but with an emphasis of keepingthe analysis as general as possible so that it is not necessary to haveparticular velocity profiles and we can address some of the issuesmentioned above. Apart from recapturing the existing results forthe Kida vortex we are also able to consider vortices with moregeneral vorticity and density distributions and make an assessmentof instabilities on the dust accumulation process.The plan of the paper is as follows: In Section 2 we give the basicequations governing the fluid model that we use. We go on to derivea partial di ff erential equation for the stream function in Section 3.In Section 3.5 we adapt the Kida solution to apply to the situationwhen the vortex has a high density core. This polytropic solution isapplicable to a Keplerian background when the aspect ratio is 7. Forother values the background has a superposed pressure extremum.In Section 4 we formulate the stability analysis governing local per-turbations to incompressible vortical flows. Perturbations localisedon streamlines are considered. In an Eulerian description, these canbe associated with a time-independent wavenumber that can leadto exponentially growing modes or in the generic case, where theperiod of circulation in the vortex is not constant a wavenumberwith a magnitude that ultimately increases linearly with time. Theformer class of modes is shown to give rise to the known insta-bilites of the Kida vortex. We generalise the instability seen there that is associated with a central saddle point in the pressure dis-tribution to more general cases and indicate how parametric insta-bility bands appear for streamlines close to the centre. Modes withwavenumbers whose magnitude ultimately increases linearly withtime cannot have amplitudes that increase exponentially with timeindefinitely, but they may undergo temporary amplification. In thespecial case of a vortex with constant period, as was consideredby Chang & Oishi (2010), the corresponding amplitudes may growexponentially with time.We go on to describe our numerical procedures and results in Sec-tion 5. We give results for steady state vortices with non-uniformvorticity sources in their cores, with and without increased centraldensity on account of a dust component, and discuss their stabil-ity. We augment the discussion of stability by considering di ff erentmodes for the analytic polytropic solution and also a simple pointvortex model which should be a limiting case applicable to stream-lines distant from the core. These models are found to behave ac-cording to expectation from the other numerical models. Often nar-row instability bands are found that incoming dust would have toencounter.In Section 6 we summarise and discuss our results, providing ar-guments why instabilities of the type found here may not preventsignificant dust accumulation in vortices with large aspect ratio. We begin by considering a single fluid model of the dust and gascirculating in a protoplanetary accretion disc. We consider unmag-netized regions of the disc such as dead zones and so neglectLorentz forces. The basic equations for the fluid are those of conti-nuity and momentum conservation. In a frame rotating with angularvelocity Ω P = Ω P ˆ k , with ˆ k being the unit vector in the fixed direc-tion of rotation (here called the vertical direction) and Ω P being themagnitude of the angular velocity, these take the form ∂ρ∂ t + ∇· ( ρ v ) = (1)and ρ (cid:32) ∂ v ∂ t + v ·∇ v + Ω P × v (cid:33) = −∇ P − ρ ∇ Φ . (2)Here, P is the pressure, ρ is the density, Φ is sum of the gravita-tional potential due to the central mass M ∗ , Φ gr = − GM ∗ / | r | andthe centrifugal potential, Φ rot = − Ω P | r × ˆk | / , with r being the po-sition vector measured from the central star. The fluid velocity is v .(For a full list of symbols see Table 1.) It is expected that the evolu-tion of the dust particle distribution can be modelled as a pressure-less fluid, which has a frictional interaction with the gas as long asthe dimensionless parameter Ω P τ s (cid:28) τ s → ∇ · v = ∂ρ∂ t + v · ∇ ρ = . (4) c (cid:13) , 000–000 tability of vortices Table 1.
Table of parameters, variables and symbolsSymbol Definition Ω P Magnitude of angular velocity Ω P ˆk Unit vector in z direction ρ, ρ (cid:48) Density and density perturbation v ≡ ( v , v , v ) Velocity v (cid:48) ≡ ( v (cid:48) x , v (cid:48) y , v (cid:48) z ) Velocity perturbation P , P (cid:48) Pressure and pressure perturbation
Φ = Φ gr + Φ rot Sum of gravitational and centrifugal potentials r Position vector measured from the cental star of mass M , r = | r | G Gravitational constant τ s Dust stopping time c s = H Ω P Sound speed with H being the disc scale height ψ, ψ , ψ Stream functions for the general flow, the backgroundflow and for the superposed vortex, ψ = ψ + ψ F ( ψ ) Arbitrary function that appears in Poisson form of mo-mentum equation, see equation (7) A ( ψ ) Bernoulli source term in equation (12) B ( ψ ) Density source term in equation (12) α, β Power-law indices in expressions for A ( ψ ) , B ( ψ ), re-spectively A , B Scaling factors in expressions for A ( ψ ) , B ( ψ ), respec-tively χ Aspect ratio of vortex patch ψ b The value of ψ evaluated on vortex boundary S Shear of background flow ω t Total vorticity in vortex patch ω v Vorticity imposed on background to produce vortexpatch b Parameter determining magnitude of central density ex-cess (polytropic model) n Power-law index in expression for density profile (poly-tropic model) β P Scaling factor in expression for pressure (polytropicmodel) ξ Lagrangian displacement S A Phase function in WKBJ ansatz (local analysis) λ Large parameter in WKBJ ansatz (local analysis) k ≡ ( k x , k y , k z ) Wave vector θ, k , ¯ t Constant angle, wavenumber scaling parameter & con-stant of integration occurring in equation (43) for k σ Eigenfrequency E Symmetric matrix on the RHS of equation (40) S ⊥ Part of S A that is a function of x and y only C K Amplitude factor and constant of integration in equa-tion (53) for particle trajectories in Kida vortices S K , φ Constant factor and angle in expression for S ⊥ applica-ble to a Kida vortex η Scaling amplitude for the Eulerian velocity perturbationin vertical stability analysis τ, W Scaled time and η in vertical stability analysis Z The quantity ∝ D ρ (cid:48) / Dt used in the vertical stabilityanalysis which satisfies a Hill equation q Constant in Hill equation for Z ω m Measure of total imposed vorticity for numerical vortexmodels ρ m Measure of total imposed mass excess for numericalvortex models γ Growth rate of instability γ Parameter taking asymptotic form γ → γ t / P = π/ω, ˜ P kida Period to circulate around a streamline in the generaland Kida cases
The density is thus conserved for fluid elements corresponding to afrozen-in dust distribution. Dissipative processes would cause thisdistribution to evolve slowly in time. However, in this paper forsimplicity we shall assume any assosiated time scale is much longerthan evolutionary tine scales of interest, such as those associatedwith dynamical instabilities. Thus we adopt equations (2), (3) and(4) throughout.
In a steady state, the equation of motion (2) reduces to v · ∇ v + Ω P × v = − ∇ P ρ − ∇ Φ . (5)In order to consider local steady state solutions within a Keple-rian disc in detail, we adopt a local shearing box with origin cen-tred on a point of interest and rotating with its Keplerian angularvelocity (see Goldreich & Lynden-Bell 1965; Regev & Umurhan2008). This specifies Ω P . A local Cartesian coordinate system isadopted with the x -axis in the radial direction, the y -axis in thedirection of shear and the z -axis normal to the disc mid-plane.For a general vector a we adopt the equivalent representations a ≡ ( a x , a y , a z ) ≡ ( a , a , a ) . To within an arbitrary constant andup to order x , the combined centrifugal and gravitational potential Φ = − Ω P (3 x − z ) /
2. The length scale associated with each dimen-sion of the box can be taken to be the vertical scale height which, inthe thin disc approximation, is assumed to be very much less thanthe local radius or distance to the central mass. z We look for solutions of equations (3)-(5) for which the fluid statevariables are independent of z and have v z ≡ v =
0. In orderto do this the z dependence of Φ is ignored. In order to satisfythe condition ∇ · v = , we adopt a stream function ψ, such that v = ( ∂ψ/∂ y , − ∂ψ/∂ x , v = (0 , − Ω P x / ,
0) and thus ψ = ψ = Ω P x / z , we remark that theymay apply to horizontal planes of an isothermal disc for which hy-drostatic equilibrium holds in the z direction (see also Lesur & Pa-paloizou 2009). In that case ρ ∝ exp( − ( Ω P z / (2 c s )) , where c s is theconstant local isothermal sound speed and in the thin disc approxi-mation the vertical scale height H = c s / Ω P (cid:28) r is implicit (see e.g.Pringle 1981). It is readily seen that a factor ∝ exp ( − Ω P z / (2 c s ))may be applied to the two dimensional solutions for ρ and P ob-tained from (3)-(5). Then when the z dependence is restored to Φ , hydrostatic equilibrium will hold in the z direction. Note thatthis feature depends on there being no x dependence in the aboveexpression for ρ , which in turn depends on the adoption of thequadratic potential Φ = − Ω P (3 x − z ) / , which is valid in thethin disk limit to within a correction of order H / r . The character-istic velocity associated with the box being c s in the thin disc limiton dimensional grounds, the characteristic velocity associated withthis correction is then expected to be of order c s H / r . This may beassumed to be small for thin enough discs even for vortices withsubsonic velocities. We remark that the above description of localsolutions in vertical hydrostatic equilibrium, with negligible ver-tical flows in the thin disc limit, has been found numerically to beapplicable to vortices generated by the Rossby wave instability (see c (cid:13) , 000–000 Railton and Papaloizou eg. Lin 2012). Note too that in the limit of zero stopping time con-sidered here, the dust is frozen into the fluid so that vertical settlingdoes not occur.
For a steady state, equation (4) becomes v · ∇ ρ = . For a two di-mensional flow this implies that the density is constant on stream-lines and is thus a function of ψ alone. Accordingly we may write, ρ = ρ ( ψ ) , where ρ ( ψ ) is an arbitrary function of ψ . This cannotbe determined if τ s = τ s = Ω P + ω ) × v = − P ρ ∇ ρ − ∇ (cid:32) P ρ + Φ + | v | (cid:33) , (6)where ω z = ˆ k · ∇ × v = −∇ ψ is the z component of vorticity, ω , asobserved in the rotating frame. Expressing quantities in terms of ψ, equation (6) becomes (cid:32) −∇ ψ + Ω P + P ρ d ρ d ψ (cid:33) ∇ ψ = −∇ (cid:32) P ρ + Φ + |∇ ψ | (cid:33) ≡ −∇ F (7)As both sides of the above have to be proportional to ∇ ψ, it followsthat F is a function of ψ alone, or F = F ( ψ ) . In the absence ofdi ff usive processes, this arbitrary function, the derivative of whichin the absence of a density gradient represents a conserved vorticity,also has to be input externally.Equation (7) can thus be written as a second order partial di ff eren-tial equation for the stream function in the form ∇ ψ = dFd ψ + Ω P + P ρ d ρ d ψ . (8)We remark that once F ( ψ ) is specified, the pressure is expressed interms of the stream function through the relation P ρ = − Φ − |∇ ψ | + F ( ψ ) . (9)Solutions of (8) corresponding to local vortices with central dustconcentrations may be sought once the arbitrary functions ρ and F , and appropriate boundary conditions are specified. In this con-text we note that, after making an appropriate adjustment to F ,equation (7) is invariant to adding an arbitrary constant to P . Forconvenience, when constructing steady state vortices numerically,we shall choose this to make the pressure zero in the limit whenthe flow becomes a pure Keplerian background flow with no addedvortex. F ( ψ ) and ρ ( ψ ) in practice We begin by separating out the solution corresponding to an undis-turbed Keplerian background flow for which ψ = ψ . To do this wewrite ψ = ψ + ψ = Ω P x + ψ , (10)where ψ corresponds to the superposed vortex. In addition we set F = − Ω P ψ + F ( ψ ) , (11) where F vanishes for the background flow. Equation (8) thenyields ∇ ψ = dF d ψ + P ρ d ρ d ψ = A ( ψ ) + P ρ B ( ψ ) . (12)Here we denote A ( ψ ) as the Bernoulli source and B ( ψ ) as the den-sity source, these both being regarded as sources of vorticity.As in the limit τ s = , these functions are invariants that have to bespecified. Accordingly we specify A ( ψ ) and B ( ψ ) so as to enable alarge class of steady state solutions with varying vorticity and den-sity profiles to be considered. The Bernoulli and density sourcesare superposed on horizontal planes on which there is a uniformbackground Keplerian flow with density ρ . They are non zero onlyon streamlines that circulate interior to a bounding streamline, withthe location of the point where it crosses the y axis being specified.The arbitrary unit of length is chosen so that this point is at (0 , , the ignorable z coordinate being from now on suppressed. The con-figurations are symmetric with respect to reflections in both the x and y axes. The unit of time is chosen so that Ω p =
1, while thearbitrary unit of mass is then chosen so that ρ = . In order toperform calculations we adopt power law functions of the form A ( ψ ) = A | ψ − ψ b | α (13) ρ ( ψ ) − ρ = B | ψ − ψ b | β . (14)Here ψ b denotes the value of ψ evaluated on the vortex core bound-ary which intersects (0 , . The functions A ( ψ ) and ρ ( ψ ) − ρ are setto be zero on streamlines exterior to those with with ψ = ψ b . Theconstants A and B are chosen to scale the total vorticity and rela-tive mass excesses associated with the Bernoulli and mass sourcesrespectively and α and β are constant indices. Note that A > α = β = B = A and B and the indices α and β . The Kida vortex provides a well-known 2D steady state analyticsolution. The core is an elliptical patch with constant vorticity ω t = − S + ω v . Here − S is the vorticity associated with the back-ground flow as seen in the rotating frame, with S = Ω P / ω v is the vorticity imposed on the back-ground. Outside the core, the vorticity is that of the backgroundflow. This corresponds to a Bernoulli source inside the core givenby A ( ψ ) = − ω v in the above notation.Kida vortices are characterised by their aspect ratio χ = a / b (cid:62) a and b are the semi-major and semi-minor axes of the el-liptical streamlines within the core. One can find steady solutionswhen the semi-major axis of the vortex is aligned with the back-ground shear and we have ω v = − S ( χ + χ ( χ − . (15)The stream function in the vortex core is given by ψ = ψ + ψ = S χ − (cid:32) χ x + y χ (cid:33) , (16)This is found by looking for a solution with elliptical streamlines c (cid:13) , 000–000 tability of vortices -1-0.5 0 0.5 1-1 -0.5 0 0.5 1 y x (a) 0 0.5 1 1.5 2 2.5 3 -1-0.5 0 0.5 1-1 -0.5 0 0.5 1 y x (b) -0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 -1-0.5 0 0.5 1-1 -0.5 0 0.5 1 y x (c) -0.07-0.06-0.05-0.04-0.03-0.02-0.01 0 Figure 1.
The pressure distributions inside Kida vortices with a range of aspect ratios, χ . In (a), the aspect ratio χ = / χ = χ = ω m of 3 . , .
125 and0 .
09 respectively. All vortices with χ > which is connected to an exterior solution of Laplace’s equationunder the condition of continuity of ψ and ∇ ψ on the vortex bound-ary (see Lesur & Papaloizou 2009, and references therein for moredetails ). From equation (9), for a fixed value of z , the pressure isthen given to within an arbitrary constant by P ρ = Ω P x − S χ − (cid:32) x χ + y χ (cid:33) + S ψ ( χ + χ ( χ − − Ω P ψ, (17)where the last two terms constitute F ( ψ ). For a Keplerian disk withuniform background S = Ω P /
2. In that case there are a rangeof pressure profiles associated with di ff erent aspect ratios as canbe seen in Figure 1. This is particularly relevant when consideringdusty gases as particles tend to drift in the direction of the pressuregradient towards pressure maxima. We remark that it is possible to consider di ff erent values of S (i.e.a non-Keplerian background flow) while retaining the potential Φ that is appropriate for a Keplerian disc. This requires the pressuregradient to be non zero in the background flow and as a conse-quence enables us to consider situations where the vortex is centredon a background where there is a pressure extremum. We commentthat this is of interest as dust is expected to accumulates at the cen-tre of a ring where there is a pressure maximum (Whipple 1972)and in addition, the Rossby wave instability can result in vorticesforming at such locations (see e.g. M´eheut et al. 2010, 2012a).For example if we set S = (cid:115) χ − χ + Ω P , (18)then we have to within a constant in the vortex core that P ρ = − S ψχ ( χ − + F ( ψ ) . (19)This makes P /ρ a function of ψ alone which will be a linear func-tion of ψ provided that F ( ψ ) is. This turns out to be useful for con-structing models with non-uniform density. However, when such asolution is matched to an exterior solution (e.g. Lesur & Papaloizou2009) the background flow will correspond to one with v y = − S x .From (18) this corresponds to the Keplerian case strictly only when χ =
7. For other values of S , consideration of the exterior Kida so-lution implies that there is an implied background pressure struc-ture which corresponds to a background pressure maximum at the coorbital radius for χ > χ < S given by (18) inside the vortex core wherethe vorticity source will be uniform. To obtain this solution we set ρ = (1 − b ( ψ − ψ b ) /ψ b ) n inside the vortex, where ψ b is the streamfunction on the core boundary and the background density is takento be unity. The quantities b and n are constants determining theprofile and magnitude of the density excess above the background.At the vortex centre ρ = (1 + b ) n while at the boundary ρ = P = β P ψ b ρ + / n / (( n + b ), where β P is a constant determinedsuch that the equilibrium conditions apply. These conditions areobtained from (8) and (19) which require that F is a linear functionof ψ . They give S ( χ + χ ( χ − = dF ( ψ ) d ψ + Ω P − n β P n + β P ψ b ( n + b (cid:32) − b ( ψ − ψ b ) ψ b (cid:33) = − S ψχ ( χ − + F ( ψ ) . (20)Together they imply that β P = Ω p − χ (cid:115) χ − , (21)which determines β P . The parameter n > b can then be chosen to scale the density excess abovethe background in the centre of the vortex provided χ > . In thispaper we have limited consideration to the case n = . We now consider the stability of steady state flows of the type in-troduced above. We find it useful to consider both the Eulerian andLagrangian formulation of the linear stability problem as they arefound to be convenient for di ff erent purposes. Following the La-grangian approach developed by Lynden-Bell & Ostriker (1967),we introduce the Lagrangian variation ∆ such that the change to astate variable Q as seen following a fluid element is ∆ Q . The La- c (cid:13) , 000–000 Railton and Papaloizou grangian displacement is given by ∆ r = ξ and we have ∆ v = D ξ Dt , (22)where D / Dt denotes the convective derivative for the unperturbedflow. Thus DDt ≡ ∂∂ t + v · ∇ . The Eulerian variation, Q (cid:48) , the change in Q as seen in a fixed coor-dinate system, is given by Q (cid:48) = ∆ Q − ξ ·∇ Q . Taking the Lagrangianvariation of the equation of motion (2), we obtain D ξ Dt + Ω P × D ξ Dt = ∆ F , (23)with ∆ F = − ∇ P (cid:48) ρ + ρ (cid:48) ρ ∇ P − ξ · ∇ (cid:32) ∇ P ρ + ∇ Φ (cid:33) . (24)The Lagrangian variation in the density is zero, thus ∆ ρ = ρ (cid:48) + ξ · ∇ ρ = . (25)Equations (23) and (25) together with the incompressibility con-dition ∇ · ξ =
0, lead to system of equations for the horizontalcomponents of ξ that is fourth order in time (see below). Whilethe above Lagrangian formulation is convenient for some aspectssuch as the analytic discussion of saddle point instability in Section4.4, the Eulerian formulation presented below is found to be moreconvenient in other contexts. The corresponding equations in terms of the Eulerian variations are D v (cid:48) Dt + Ω P × v (cid:48) + v (cid:48) · ∇ v = − ρ ∇ P (cid:48) + ρ (cid:48) ρ ∇ P and D ρ (cid:48) Dt = − v (cid:48) · ∇ ρ, (26)which is a system that is third order in time. It is lower order thanthe Lagrangian system on account of trivial solutions correspond-ing to a relabelling of fluid elements being present in the latter case(see Friedman & Schutz 1978). Equations expressed in terms ofEulerian variations were found to be simpler to use when analysingvertical stability in Section 4.9 and for the same reason were solvednumerically when considering vortex stability in Section 5.2. We consider perturbations that are localized on streamlines whichhave short wavelengths in the directions perpendicular to the unper-turbed velocity, but can have a long wavelength in the direction ofthe unperturbed velocity. The latter is a natural outcome of shear-ing motions. To do this we begin by assuming that any perturbationquantity takes the form ∆ Q = ∆ Q exp ( i λ S A ) . (27)Here we adopt a WKBJ ansatz with the phase function S A left ar-bitrary for the time being and the constant λ taken to be a largeparameter. For a discussion of the approach followed here in a va-riety of contexts, see Lifschitz & Hameiri (1991), Sipp & Jacquin(2000) and Papaloizou (2005). The e ff ective wavenumber k = λ ∇ S A (28) then has a large magnitude. The amplitude factor ∆ Q is the WKBJenvelope. Due to the rapid variation of the complex phase λ S A , onecan perform a WKBJ analysis, such that the state variables ∆ Q are expanded in inverse powers of λ. The lowest order term in ξ is constant, while the lowest order term in P (cid:48) is ∝ λ − . To lowestorder (23) gives DS A Dt = . (29)When working to the next order, only the variation of the rapidlyvarying phase S A needs to be considered when taking spatial deriva-tives, apart from when considering expressions involving the oper-ator D / Dt ≡ ∂/∂ t + v · ∇ as this annihilates S A . Accordingly thecontribution v · ∇ ( ∆ Q ) must be retained. Noting the above, we cansubstitute perturbations of the form (27) into the governing equa-tions and remove the factor exp ( i λ S A ) , thus obtaining equations forthe lowest order contribution to the quantities ∆ Q alone. For easeof notation we drop the subscript 0 from now on.Following this procedure (23) gives D ξ Dt + Ω P × D ξ Dt = − i k P (cid:48) ρ − ξ · ∇ ρρ ∇ P + ξ · ∇ ( v · ∇ v + Ω P × v ) , (30)where we recall that the lowest order term for P (cid:48) ∝ λ − . In additionto this, the incompressibility condition gives k · ξ = . (31)Using (28) and (29) we can also find an equation for the evolutionof k in the form D k Dt = − k j ∇ v j . (32)Equations (30), (31) and (32) give a complete system for the evo-lution of ξ and k , after the elimination of P (cid:48) , as an initial valueproblem. Because the evolution consists of advection of data alongstreamlines, it is possible to consider disturbances localized on in-dividual streamlines (Papaloizou 2005). Localization amplitudesare una ff ected by the evolution considered. In general one couldstart with an arbitrary initial S A and then k would depend on time. A relatively simple class of solutions for k can be obtained by set-ting S A to be independent of time and a function of quantities con-served on unperturbed streamlines so that v · ∇ S A = k is fixed for all time andwe only have to solve for ξ . For the simple case of a two dimen-sional vortex with initial state-independent of z , we may take S A = g ( ψ ) + k z z /λ. (33)Here ψ is the unperturbed stream function, g is an arbitrary functionand k z is the constant vertical wavenumber. We assume for now that k z (cid:44)
0, appropriate to the physically realistic case where perturba-tions are localized in z .We further remark that although the above form of k is not themost general solution of (32), apart from when the velocity is lin-ear in the coordinates as for the Kida vortex, other solutions aresuch that the magnitude of the wavenumber ultimately increaseslinearly with time. In that situation we expect that although theremay be temporary amplification, the system may not ultimatelyshow growth of linear perturbations exponentially with time.. Thissituation is already well known in the context of the shearing box(see Goldreich & Lynden-Bell 1965). c (cid:13) , 000–000 tability of vortices For now we continue the discussion adopting the time-independentform of k derived from S A given by (33) and return to the discus-sion of more general k and the special nature of the Kida vortex inSection 4.6 below.We may use the vertical component of (30) together with (31) toeliminate P (cid:48) and ξ z and thus obtain a pair of equations for ( ξ x , ξ y ) ≡ ( ξ , ξ ). These can be written in the form (cid:32) δ ij + k i k j k z (cid:33) D ξ j Dt + (cid:32) (cid:15) i j Ω P + k i k z Dk j Dt (cid:33) D ξ j Dt + k i ξ j k z D k j Dt = H i (34)where H = − ξ · ∇ ρρ ∇ P + ξ · ∇ ( v · ∇ v + Ω P × v ) . (35)We remark that the neglect of vertical stratification in this calcula-tion can be generally justified if it is assumed that k z / ( k x + k y ) islarge, otherwise the modes can be assumed to be localized in thevicinity of the midplane where the vertical startification is least.Note that as only horizontal components are considered, the sum-mation is for j = , D k j Dt = ∂∂ x q (cid:34) k µ (cid:32) v µ ∂ v q ∂ x j − v q ∂ v µ ∂ x j (cid:33) (cid:35) . (36)We remark that although there is an arbitrary function g in the defi-nition of k used in this section, because the derivatives in (34) cor-respond to advecting around streamlines and dg / d ψ is constant onstreamlines, the latter quantity e ff ectively behaves as a multiplica-tive constant merely scaling the magnitude of the wavenumber. The equivalent of (34) written in terms of the Eulerian variations( v (cid:48) x , v (cid:48) y ) ≡ ( v (cid:48) , v (cid:48) ) obtained from (26) is (cid:32) δ ij + k i k j k z (cid:33) Dv (cid:48) j Dt + (cid:15) i j Ω P v (cid:48) j + v (cid:48) j ∂ v i ∂ x j + k i v (cid:48) j k z Dk j Dt = ρ (cid:48) ρ ∂ P ∂ x i . (37)In addition, the Eulerian density variation satisfies D ρ (cid:48) Dt = − v (cid:48) · ∇ ρ. (38) The analyses in the previous sections reduce the stability problemto solving an initial value problem of integrating a set of simulta-neous ordinary di ff erential equations around streamlines. The inde-pendent variable measures the location on a streamline.As the unperturbed motion on a streamline is periodic, these equa-tions have periodic coe ffi cients through their dependence on k andthe gradient of ρ etc. Thus Floquet theory may be applied (eg. Whit-taker & Watson 1996). According to this, if some internal modewith a natural oscillation frequency dependent on k is describedby (34), unstable bands of exponential growth are expected as k is varied to allow resonances of this frequency with the frequencyof motion around the streamline. To see how this can come aboutwe shall specialise to the case when | k z | (cid:29) (cid:113) k x + k y which cor-responds to the so called horizontal instability, as in this limit themotion occurs in uncoupled horizontal planes. Taking the limit k z → ∞ , the horizontal components of (34) yielda pair of equations for the horizontal components of displacement,while the z component and the condition ∇ · ξ = ξ z →
0. Wethus have D ξ Dt + Ω P × D ξ Dt = − ξ · ∇ ρρ ∇ P + ξ · ∇ ( v · ∇ v + Ω P × v ) , (39)where we now have ξ = ( ξ x , ξ y ) . Using equation (5), this can bewritten in the equivalent form D ξ Dt + Ω P × D ξ Dt = − (cid:32) ξ ρ · ∇ (cid:33) ∇ P − ( ξ · ∇ ) ∇ Φ (40)We remark that (40) becomes an equation with constant coe ffi cientsfor ξ when ρ is constant and P and Φ are quadratic in x and y . Asthis is always the case arbitrarily close to the centre of any regularvortex where there is a stagnation point, there are generic conse-quences. In the horizontal limit we solve (39) by setting ξ = ξ exp i σ t ,where ξ is a constant vector, and finding an algebraic equationfor σ . In doing this we find it convenient to define the symmetricmatrix E by writing the right hand side of (40) as ( E ( ξ x , ξ y ) T ) T . Theequation for σ is readily found to be given by σ − σ (4 Ω P − Tr ( E )) + det( E ) = , (41)with Tr and det denoting the trace and determinant respectively.A su ffi cient condition for instability, or at least one complex root for σ , is that det E <
0. Note that in the limit approaching the vortexcentre, the elements of E are given by E i , j = − ρ ∂ ∂ x i x j (cid:32) P − ρ Ω P x (cid:33) (42)evaluated in the limit | r | →
0. This condition is equivalent to P − ρ Ω P x / P has a saddle point that appears as a maximum alongthe x -coordinate line and a minimum along the y -coordinate line,as can be seen to occur directly from the analytic solution for Kidavortices with 3 / < χ < , and for the vortex illustrated in the bot-tom panels of Figure 2. Saddle points of this type are genericallyassociated with instability in all cases, independently of density orvorticity profile. We recall that (40) applies in the limit k z → ∞ and that it becomesan equation with constant coe ffi cients when ρ is constant and P and Φ are quadratic in the coordinates. This is always the case for anystreamline in the core of a Kida vortex so moving away from thecentre has no e ff ect. However, in more general cases P , ρ and hence E will be represented by a power series in x and y . Thereforefor a fluid element, E will be periodic with period one half of thatassociated with circulating around the streamline. Thus equation(40) will have coe ffi cients that are periodic in time and parametricinstability becomes possible.Close to the vortex centre the time-dependence can be treated asa perturbation and parametric instability derived analytically, al-beit in terms of unknown coe ffi cients. This is done in detail in Pa-paloizou (2005), where a procedure is followed that can be applied c (cid:13) , 000–000 Railton and Papaloizou − ∇ ψ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 16.75Model˜ P kida (a) { α, β, ρ m , ω m } = { , , , . } , the Kida vortex where χ = − ∇ ψ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 12.55Model˜ P kida (b) { . , , , . } , with χ = . − ∇ ψ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 10.88Model˜ P kida (c) { . , , , . } , with χ = . Figure 2.
Vortices with ω m = .
09 and density equal to the background value. From left to right we have: (i) the vorticity distribution, (ii) the ψ distribution,(iii) the pressure distribution and (iv) a plot of the period P around a streamline against the value of the positive y -coordinate where it intersects the y axis.Increasing α from 0 to 4 results in a vortex that is stronger (i.e. sheared less by the background flow) and accordingly has a smaller aspect ratio. Note that allthese vortices, except the one illustrated in the bottom panels, have a pressure maximum at the vortex centre. The pressure distribution in the latter vortex hasa saddle point. There is also significant shear, as indicated by the variation in the period P to circulate around a streamline. directly to the problem considered here, so for the sake of brevitywe refer the reader to Appendix B of that paper.We note that parametric instability is first expected to occur whenthe epicyclic oscillation period is equal to the period to circulatearound the streamline. Higher order bands are expected to be gen-erated when the ratio of epicyclic oscillation period to circulationperiod is 1 / , / , . . . For a vortex with a core like the Kida vortex,these resonances occur when χ = . , .
89 and 7 .
32 respectively(Lesur & Papaloizou 2009).
The form of wavenumber described in Section 4.2.1 di ff ers fromthat adopted in the stability analysis of the Kida vortex core givenby Lesur & Papaloizou (2009). However, because the governingequations (37), a with ρ (cid:48) =
0, are independent of the choice of theorigin of time, they turn out to be equivalent. Results are illustratedin section 5.2.Equations (30) - (32) of Section 4.2 apply in this case with thephase function being given by λ S A = k ( t ) · r where the wave vector k is given by k ( t ) = k (cid:16) χ sin( θ ) sin (cid:0) φ ( t ) (cid:1) , sin( θ ) cos (cid:0) φ ( t ) (cid:1) , cos( θ ) (cid:17) , (43)where k and θ are constants and φ ( t ) = Ω P / (2( χ − t − t ) , t being constant. It is important to note that (43) only works when thevelocity components are linear functions of the coordinates. Notethat in spite of this restriction Chang & Oishi (2010) used it inthe problem of the stability of a vortex with a density gradient, forwhich this condition would not be expected to be self-consistentlysatisfied. Thus an assessment of the situation that occurs when thevelocity in the vortex is not a linear function of the coordinatesshould be carried out. S A for an arbitrary two dimensionalincompressible vortical flow The general form of S A is obtained from the solution of equation(29) in the form DS A Dt = . (44) c (cid:13) , 000–000 tability of vortices ρ ( ψ ) − ρ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 16.49Model˜ P kida (a) ( α, β, ρ m ) = { , . , . } , with χ = . ρ ( ψ ) − ρ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 16.19Model˜ P kida (b) { , . , . } , with χ = . ρ ( ψ ) − ρ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 16.01Model˜ P kida (c) { , . , . } , with χ = . Figure 3.
Vortices with varying density excess for which ω m = . α = β =
1. From left to right we have: (i) the density distribution, (ii) the ψ distribution, (iii) the pressure distribution and (iv) a plot of the period P around a streamline against the value of the positive y -coordinate where it intersectsthe y axis. For the case when the background flow is independent of z , we canwrite S A = S ⊥ ( x , y ) + k z z /λ , where k z is the component of k in the z direction and S ⊥ satisfies ∂ S ⊥ ∂ t + ∂ψ∂ y ∂ S ⊥ ∂ x − ∂ψ∂ x ∂ S ⊥ ∂ y = . (45)The general solution of (45) requires that S ⊥ be a function only ofquantities that are invariant of the particle trajectories obtained bysolving dxdt = ∂ψ∂ ydydt = − ∂ψ∂ x . (46)The solutions for x and y define orbits or streamlines that are peri-odic in time and on which ψ is constant. The period is 2 π/ω, where ω in general will be a function of ψ . Quantities such as x and y canbe expressed as a Fourier series in the form x = ∞ (cid:88) n = −∞ x n ( ψ ) exp (i n φ ) , (47) y = ∞ (cid:88) n = −∞ y n ( ψ ) exp (i n φ ) , (48) where φ = ω ( t − t ) and t is a constant on an orbit that can be takento be the time at which x passes through its maximum value.The general solution of equation (45), which states that S ⊥ is a con-stant on an orbit defining a streamline, is that S ⊥ is an arbitraryfunction of ψ and t . As the orbits are periodic, this function shouldbe periodic in t with period 2 π/ω . Accordingly S ⊥ can also bewritten as a Fourier series in the form S ⊥ = ∞ (cid:88) n = −∞ C n ( ψ ) exp ( − i n ω t ) = ∞ (cid:88) n = −∞ C n ( ψ ) exp (i n ( φ − ω t )) . (49)We may now find k = λ ∇ S A leading to k x = λ (cid:32) ∂ψ∂ x ∂ S ⊥ ∂ψ + ω ∂ y ∂ψ ∂ S ⊥ ∂φ (cid:33) k y = λ (cid:32) ∂ψ∂ y ∂ S ⊥ ∂ψ − ω ∂ x ∂ψ ∂ S ⊥ ∂φ (cid:33) . (50)In obtaining the above we note that quantities are either expressedas functions of ( x , y ) ≡ r or ( φ, ψ ) as independent variables. Trans-forming between these representation is facilitated by noting that v = ω∂ r /∂φ and that the Jacobian ∂ ( φ, ψ ) /∂ ( x , y ) is equal to ω . Weremark that when only terms with n = ∂ S ⊥ /∂φ = c (cid:13) , 000–000 Railton and Papaloizou ρ ( ψ ) − ρ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 13.88Model˜ P kida (a) ( α, β, ρ m , ω m ) = { . , , . , . } , with χ = . ρ ( ψ ) − ρ y -1-0.500.51-1 -0.5 0 0.5 1 y x -0.100.10.20.30.40.50.60.70.8 ψ -1-0.500.51-1 -0.5 0 0.5 1 y x P P y ˜ P kida = 13.55Model˜ P kida (b) { . , , . , . } , with χ = . Figure 4.
Vortices with a non-uniform Bernoulli vorticity source , with ω m = .
09 and non-zero density enhancement parameter ρ m . From left to right wehave: (i) the density distribution, (ii) the ψ distribution, (iii) the pressure distribution and (iv) a plot of the period P around a streamline against the value of thepositive y -coordinate where it intersects the y axis. The vorticity profiles in these vortices are non uniform, resulting in significant variation of the period forcirculating around internal streamlines and hence significant internal shear. On the other hand if terms with n (cid:44) d ω/ d ψ (cid:44) , thewavenumber is expected to depend on time as well as on x and y .To emphasise this point we rewrite (50) in the form k x = λ (cid:32) ∂ψ∂ x ∂ S ⊥ ∂ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:32) ω ∂ y ∂ψ − d ω d ψ t ∂ψ∂ x (cid:33) ∂ S ⊥ ∂φ (cid:33) k y = λ (cid:32) ∂ψ∂ y ∂ S ⊥ ∂ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − (cid:32) ω ∂ x ∂ψ + d ω d ψ t ∂ψ∂ y (cid:33) ∂ S ⊥ ∂φ (cid:33) . (51)Here | denotes that a derivative is to be taken ignoring the ψ de-pendence of ω , which is now taken into account by the terms witha factor t . In the limit t → ∞ we have k x + k y ∼ λ (cid:32) d ω d ψ (cid:33) (cid:32) ∂ S ⊥ ∂φ (cid:33) | v | t . (52)The right hand side of the above is the product of | v | t and a factorthat is constant on a streamline indicating that the magnitude of thewavenumber increases to arbitrarily large values at all points on it. = ± for a Kida vortex For a Kida vortex in a Keplerian background, ω is constant so thatterms ∝ t in (51) are absent. We also note that only terms with n = ± x = C K ψ / cos( φ ) , y = − C K χψ / sin( φ ) , (53)with the amplitude factor C K being given by C K = (4( χ − / (3 Ω P χ )) / . (54)We now adopt S ⊥ = − S K ψ / cos( φ + φ − ω t ) , where S K and φ = ω t − π/ ω = Ω P / (2( χ − k = λ ∇ S A .With the help of (51), we find that k is indeed given by (43) pro-vided that we identify k z = k cos( θ ) and λω C K S K / = k sin( θ ).We confirm that although this time-dependent wavenumber was de-rived from terms with n = ± n =
0, one obtains equations governing the stabilityof a Kida vortex that are independent of which is chosen. This isbecause, for the Kida vortex, the equations are invariant to a shiftin the origin of time on a streamline and thus independent of t . Thismeans that we may specify ( k x , k y ) ∝ ∇ ψ in either case.However, it is important to note that in the generic case for which ω is not the same on di ff erent streamlines, one must adopt the time-independent form if modes growing exponentially with time in theusually expected manner are to be obtained. If a wavenumber in-creases linearly with time only temporary exponential growth isexpected (eg. Goldreich & Lynden-Bell 1965), with perturbationsultimately subject to at most power law growth with time thus re-quiring a nonlinear analysis to determine the outcome. This can beshown to be the case for the systems considered here (see below).Accordingly we extend the linear stability analysis for the Kidavortex to more general cases by adopting the time-independentwavenumber ( n = ω is constant.We comment that the situation here is analogous to the one that oc-curs for a di ff erentially rotating disc for which fluid elements orbiton circles. The time-independent wavenumber modes here corre-spond to axisymmetric modes there. The modes with wavenumbersthat increase with time correspond to non axisymmetric modes inthe disc case. c (cid:13) , 000–000 tability of vortices We now discuss vertical stability for which k z = k z / ( k x + k y ) → t → ∞ for choices of wavenum-ber that ultimately increase linearly with time. Accordingly, in thislimit the discussion of vertical stability given here should apply. Weadopt the Eulerian formulation leading to the linear equations (26).The linearized incompressibility condition gives k · v (cid:48) =
0. Thus wemay set v (cid:48) x = η k y , and v (cid:48) y = − η k x for some scalar η . The linearizedform of the condition that ρ is fixed on fluid elements gives D ρ (cid:48) Dt = − v (cid:48) · ∇ ρ, or equivalently D ρ (cid:48) Dt = η k · v d ρ d ψ . (55)Eliminating P (cid:48) from the x and y components of the first of the equa-tions (26) and making use of (32) we obtain a relation between η and ρ (cid:48) in the form D ( η | k | ) Dt = − ρ (cid:48) ρ ( k × ∇ P ) · ˆ k . (56)Equations (55) and (56) provide a pair of first order ordinary di ff er-ential equations for which the integration is taken around stream-lines. We note that in the general case from (51) it is seen that k isnot a periodic function of time and so they do not lead to a Floquetproblem. To discuss stability further set η = W / | k | and a new scaled timevariable τ , defined through d τ = | k | / dt , to get the pair of equations D ρ (cid:48) D τ = W k · v | k | / d ρ d ψ , (57) DWD τ = − ρ (cid:48) ρ ( k × ∇ P ) · ˆ k | k | / . (58)In this form, provided ∂ S ⊥ /∂φ (cid:44) , the coe ffi cients of W in (57) and ρ (cid:48) in (58) tend to zero for large τ in such a way that there can be noexponentially growing solutions that apply at large times (althoughweaker growth could occur). Note in this context that although k increases linearly with time, equation (51) implies that k · v remainsbounded. We now consider the situation for a vortex which is assumed tohave stream function given by (16) and non-constant ρ = ρ ( ψ ),namely the polytropic model considered in Section 3.5. In this casethe pressure P = P ( ψ ) and d ω ( ψ ) / d ψ =
0. On account of the lat-ter relation we can adopt the solution given by (43) for θ = π/ , which corresponds to k z =
0, without encountering problems ofthe wavenumber increasing linearly with time. The coordinates ona streamline can be specified as indicated in Section 4.8.1. Thenequations (55) and (56) can be combined to give a second orderordinary di ff erential equation for ρ (cid:48) of the form DDt (cid:34)(cid:16) χ + + ( χ −
1) cos 2( ω t − φ ) (cid:17) D ρ (cid:48) Dt (cid:35) = − Ω P χ sin ( ω t − φ ) ψ ( χ − ρ dPd ψ d ρ d ψ ρ (cid:48) , (59)Note that the terms on the right hand side of (59) multiplying ρ (cid:48) areconstant on a streamline. Thus if we set Z = (cid:16) χ + + ( χ −
1) cos 2( ω t − φ ) (cid:17) D ρ (cid:48) Dt , (60) we find that Z satisfies a form of Hill’s equation that can be writtenin the form D Z Dt = − Ω P χ sin ( ω t − φ ) ψ ( χ − (cid:18) χ + + ( χ −
1) cos 2( ω t − φ ) (cid:19) ρ dPd ψ d ρ d ψ Z . (61)This equation can be interpreted as describing the evolution of agravity wave with a time-dependent wavenumber. It can rewrittenin the form D Z Dt = − q + χ + ( χ −
1) cos 2( ω t − φ ) Z , (62)where q is a constant that can be scaled up to a maximum value q max = Ω P χ ( χ − − ψρ − ( dP / d ψ )( d ρ/ d ψ ) by adjusting the valueof sin ( ω t − φ ) through the specification of ω t − φ . The quan-tity q max can be interpreted as the square of a buoyancy frequency.The possibility of parametric instability is expected when this fre-quency is large enough to be comparable to the vortex frequency ω . Solutions of this are discussed below. We perform calculations to obtain steady state solutions for ψ forvortices on a 512 × x and y occupying the interval [ − , B = ∇ ψ = A ( ψ ) . (63)We start by defining the Bernoulli source term A to be a non-zeroconstant inside the unit circle passing through (0 ,
1) and zero out-side this. We then solve equation (63) by utilising the 2D Green’sFunction, obtaining ψ ( r ) = π (cid:34) log | r − r (cid:48) | A ( ψ ) d r (cid:48) (64)We then add in the streamfunction for the background shear, ψ = Ω P x , to get the total streamfunction ψ . Setting the bound-ing streamline of the vortex core to pass through (0 , , we thenrecalculate the boundary, rescale the Bernoulli source function asindicated below, and apply the Green’s Function again to obtain anupdated solution for ψ . We define πω m = (cid:82) A ( ψ ) dS , where the integral is taken over thetotal area inside the vortex core. This quantity is the magnitudeof the circulation around the vortex contributed by the non back-ground source. In order to rescale the Bernoulli source, we recal-culated A , adopting the required functional form, after every iter-ation. We then renormalized it, so as to ensure that πω m remainedfixed. This procedure was iterated until there was no discernible c (cid:13) , 000–000 Railton and Papaloizou χ θ log γ (a) Kida vortex χ θ log γ (b) Polytropic model Figure 5.
The stability of the Kida vortex (a) and polytropic model with uniform density in a non Keplerian background (b). χ θ log γ (a) { α, β, ρ m } = { . , , } , with y max = . χ θ log γ (b) { α, β, ρ m } = { . . , } , with y max = . Figure 6.
The stability of motion in vortices with α = .
25, on di ff erent streamlines. The parameter ω m is varied to produce the range of aspect ratios shownfor this streamline. di ff erence in successive streamfunction and pressure contours. Forvortices with B = , this required generally about 50 iterations.Comparison with the analytic solution for the Kida vortex case (seeSection 3.4) enabled our numerical procedure to be checked (seebelow). In order to obtain solutions with increased density in the centralparts of the vortex core we have also found solutions of equation(12) with both A (cid:44) , and B (cid:44)
0. For these cases the pressuredistribution has to be calculated at each iteration. This is found fromequation (9) which yields P ρ = Ω P x − |∇ ψ | − Ω P ψ + (cid:90) ψψ b A ( ψ (cid:48) ) d ψ (cid:48) . (65) Our numerical models are characterised by input values of α and β that are required to implement equations (13) and (14) respectively.It is convenient to scale the mass per unit length added to the vortexby imposing a fixed value of ρ m = (cid:82) [ ρ ( ψ ) − ρ ] dS and, as before πω m = (cid:82) A ( ψ ) dS , which is a measure of the additional imposedcirculation. In this case the constants (13) and (14) are rescaledafter every iteration to ensure the required fixed values of ρ m and πω m are maintained.Solutions are thus characterised by the parameters { α, β, ρ m , ω m } ,with the last two quantities expressed in our dimensionless units.Convergence for models with constant density was straightforward,although resolution issues arise for a fixed grid when α becomes toolarge on account of peaking of the vorticity profile in the centre ofthe vortex.For models with increased density inside the vortex core, conver-gence was more di ffi cult and required a starting model close to the c (cid:13) , 000–000 tability of vortices final one. In order to deal with this we began by imposing onlysmall increases to either of, or both, β and ρ m from their values ap-propriate to an existing solution. These changes moved them in thedirection of our target parameters. This procedure was especiallynecessary for the cores of vortices with non-uniform vorticity pro-files. Having done this, a new form of ψ was obtained as aboveby iterating the Green’s function solution N (cid:62)
10 times, at whichpoint additional small increments of the order of 1% were made to β and ρ m and the process repeated until the target values were at-tained. After that the Green’s function solution could be iterated toconvergence. The procedure we used was to solve the system of equations (37 )together with equation (38) for the Eulerian perturbations as an ini-tial value problem. As the integration is around steady state stream-lines, state quantities need to be specified on them. A particularstreamline has to be specified by identifying an initial location eg.(0 , y max ) . The location on it as the integration proceeds is then spec-ified by solving the equations
DxDt = ∂ψ∂ yDyDt = − ∂ψ∂ x . (66)The growth rate of any instability present was obtained by solvingthe additional equation D γ Dt =
12 ln | v (cid:48) | . (67)For a system with growth rate γ, we expect that ultimately γ → γ t /
2. Accordingly we determined γ by making a parabolic fit to γ at large times. We found that integrations running for 1000 cir-culations around streamlines could detect growth rates down to ∼ . Ω P .In order to get a global view of stability we need to consider mod-els over a large range of aspect ratio and k z . For cases with time-independent wavenumber, the latter quantity is specified throughthe ratio (cid:113) k x + k y / | k z | ≡ tan θ evaluated at (0 , y max ) on the streamline. We remark that this quantity is ∝ | v | on the streamline and sois minimised at the chosen location. Regions of instability often re-quire high resolution in this phase space in order to be adequatelymonitored. Typically for a model with specified α, β and ρ m , werequired a 300 ×
300 grid in the ( χ, θ ) plane (see results below).In order to perform these time-consuming calculations, we usedthe fact that all required quantities on a given streamline can bedetermined in terms of ψ . To facilitate this we made two dimen-sional least square polynomial fits of up to sixth order to ψ insidethe vortex core for a grid of calculated models. We note that thisshould give a good representation for streamlines close enough tothe centre and the onset of parametric instability (see section 4.5).In practice we found that even fourth order polynomial fits gaveresults not significantly di ff erent from those presented below. We show the streamlines for vortices with ω m = .
09 and densityequal to the background value obtained using our numerical pro-cedure in Figure 2. These include a Kida vortex with χ = α = α =
2. In Figure 2 the vorticity distribution, the ψ contours and the pressure distribution are shown.In addition, the period ˜ P required to circulate around a streamlineis plotted against the value of the positive y -coordinate y max , whereit intersects the y -axis. The values of ˜ P kida indicated are calculatedanalytically using the formula˜ P kida = (cid:73) ds |∇ ψ | = π χ − , (68)which is independent of the streamline chosen (the Kida solutionhas constant period throughout the vortex patch). The results forthe Kida vortex are in good agreement with analytic expectation.Note that the periods we plot were obtained by locating coordinateextrema from the results of numerical integrations of fluid particlesmoving around streamlines stored on a coarse temporal grid. Thisresults in some low-level jitter at a relative level ∼ − . This isalso a measure of the departure from a constant value of the periodobtained numerically in the Kida vortex case.Note too that a significant non-constant value of the period to circu-late around a streamline indicates the presence of shear inside thevortex. The latter becomes more noticeable as α increases. Apartfrom in the case of the vortex illustrated in the bottom panels, thereis a pressure maximum at the centre and in the latter case there is asaddle point. We find that the transition from central saddle point tomaximum occurs at smaller aspect ratios as α increases and the vor-ticity distribution becomes more centrally concentrated. This shiftsthe central saddle point instability described above to central valuesof χ <
4. We recall that cases for which the pressure distributionhas a maximum at the centre of the vortex are of interest as they areexpected to attract dust (Whipple 1972; Cardoso et al. 1996). In ad-dition, streamlines for vortices with centrally concentrated vorticitysources tend to become pinched towards the y -axis as compared tothe Kida case as one moves out from the centre. This leads to arelative increase in the circulation period.Vortices with varying density excess [for which ω m = . α = β =
1] are similarly illustrated in Figure 3. In these cases be-cause the Bernoulli source is uniform, the density distribution isshown. For these models the central density ranges between 1 . . ∼ . Ω P . However, this could be significant for linear perturba-tions (see below).Finally vortices with a non-uniform Bernoulli vorticity source, with α = . ω m = .
09, and non-zero density enhancement parameter ρ m are illustrated in Figure 4. The vorticity profiles in these vorticesare non-uniform, resulting in significant variation of the period andhence significant internal shear. The stability of both Kida vortices and the polytropic model withuniform density in a non-Keplerian background (see Section 3.5),determined using our procedures as viewed in the ( χ, θ ) plane areillustrated in Figures 5(a) and (b) respectively. The results for theKida vortices are in good agreement with previous work (see Lesur& Papaloizou 2009). The results for the polytropic model are quali-tatively similar but with the region of strong saddle point instabilityshifted to smaller values of χ < ff ects of introducing a variable vorticity pro-file in the core, we illustrate the stability of vortices with α = . c (cid:13) , 000–000 Railton and Papaloizou χ θ log γ (a) { α, β, ρ m } = { . , , } χ θ log γ (b) { α, β, ρ m } = { . , , } χ θ log γ (c) { α, β, ρ m } = { . , , } χ θ log γ (d) Point vortex Figure 7.
As in Figure 6 but for di ff erent values of α , with y max = .
85 as standard. In addition the stability of the point vortex model is shown. In this casethe streamline aspect ratio is varied by changing the strength of the point vortex. in the ( χ, θ ) plane in Figure 6. The stability of motion on stream-lines with y max = . y max = .
85 (Figures 6(a) and (b) respec-tively) is shown. The parameter ω m is varied to produce the rangeof aspect ratios shown for this streamline. As indicated in Section4.5, moving out from the vortex centre we expect parametric insta-bility to occur for χ ∼ .
65. This is visible for the streamline with y max = . χ extends towards the χ -axis at smaller values of χ . Inaddition the region associated with saddle point instability shifts tosmaller values of χ .In Figure 7 we show similar results to those presented in Figure 6for streamlines with y max = .
85 and α = { . , . , . } . As α increases, the instability band originating from χ = .
65 widensand moves to smaller values of χ , while the small χ region asso-ciated with saddle point instability eventually disappears. Severaladditional instability bands appear at larger values of χ . We alsoconsidered a vortex point source model. This adopts the stream-function ψ = K ln | r | + Ω x / α . The aspect ratio of the streamline for a fixed y max = . K . The stabilityproperties are obtained using the same procedures as for the othermodels. Results are given in Figure 7. The behaviour is very simi-lar to that seen for the case with α =
2, with the ( χ, θ ) plane fillingup with instability bands. Growthrates at larger values of χ (cid:38) ∼ . Ω P . We comment that although these models do not haveany density excess, because they are relevant to streamlines out-side a high density core they are of generic significance for vorticesaccumulating dust in their cores.The stability of vortices with α = β = ρ m areillustrated in Figure 8. Similarly, results for vortices with α = { . , . } and non-uniform density excess are illustrated in Fig-ure 9. For these calculations, χ was varied by changing ω m while ρ m was chosen such that a fixed mass per unit length was addedto the vortex with specified α and β for all χ . As the vortex coreboundaries all pass through (0 , χ .The central density for the vortices illustrated in Figures 8 and 9 c (cid:13) , 000–000 tability of vortices χ θ log γ (a) { α, β, ρ m } = { , , . } , streamline passing through (0 . , . χ θ log γ (b) { α, β, ρ m } = { , , } , streamline passing through (0 , . χ θ log γ (c) { α, β, ρ m } = { , , } , streamline passing through (0 , . χ θ log γ (d) { α, β, ρ m } = { , , } , streamline passing through (0 , . Figure 8.
The stability vortices with α = , β = ρ m . Results are shown for streamlines with y max = { . , . , . } . The parameter ω m is variedto produce the range of aspect ratios shown. are shown as a function of χ in Figure 10. Central densities rangefrom 3 to 12 times the background level at large values of χ .In Figure 8 we can see the e ff ect of introducing a small density ex-cess on the stability of a Kida vortex. For the case with ρ m = . y max = . .
85. In theformer case where departures from a quadratic stream function aresmaller, a parametric instability band is seen emanating from theexpected location χ = .
65. This appears to attempt to connectwith the band originating from large χ and θ in the Kida vortexcase, leaving a region that is very weakly unstable, or possibly sta-ble, between them. Bands with such regions seem to be a commonfeature in these calculations. They require time consuming calcu-lations at high resolution to locate them. In our case we found itimpractical to resolve instabilities with growth rates < . Ω P .When the streamline with y max = .
85 is considered, the band orig-inating from χ ∼ .
65 has broadened to produce an unstable regionwith growth rate ∼ . Ω P for 4 < χ <
5. Two additional narrowbands appear at larger χ with characteristic growth rates ∼ . Ω P .When the ensemble of streamlines is considered for this model, it is di ffi cult to find any χ for which there is stability for all θ . Theparameter ω m is varied to produce the range of aspect ratios shown.The lower panels of Figure 8 illustrate the stability of a model with ρ m = .
0. This higher density model shows qualitatively similarbehaviour but with the appearance of a larger number of weak in-stability bands and, for y max = . , an instability band for χ < α = .
25 and α = . . , . χ , ∼ . Ω P beingsimilar to the concentrated Bernoulli source cases without densityexcesses.In Figure 11 we show the stability of the polytropic models with n = χ, θ ) plane with central densities 3 and 5 times thebackground level. The results are for streamlines close to the coreboundary with y max = .
95. We recall that for these models, vor-tices associated with values of χ (cid:44)
7, have to be considered to be c (cid:13) , 000–000 Railton and Papaloizou χ θ log γ (a) { α, β, ρ m } = { . , , . } χ θ log γ (b) { α, β, ρ m } = { . , , . } Figure 9.
The stability of vortices with α = . . , . immersed in a non-Keplerian background (see Section 3.5 ). Thepresence of instability bands will be noted, three in the lower cen-tral density case and four in the higher central density case. In thisrespect the results are qualitatively similar to models with centraldensity excess immersed in a Keplerian background for all χ . Thecharacteristic growth rates for χ ∼ ∼ . Ω P .Finally, we have studied the parametric instability that can occurwhen the vertical stability of the polytropic model is considered(see Section 4.9.2) because that has no internal shear. This connectsto a Keplerian background for χ =
7. Accordingly we consider thatcase. The growth rate is plotted as a function of q / Ω P in Figure 12.It will be seen that there is instability for q (cid:38) . Ω P and the growthrate reaches a maximum of ∼ . Ω P for q = . Ω P . This is aboutan order of magnitude larger than the growth rates illustrated in Fig-ure 11 for modes with time-independent wavenumber at this valueof χ . Thus modes with k z = q for the polytropic modelsoccurs on the core boundary and is q max = n b Ω P /
4. Thus stabilityis indicated for su ffi ciently small n b < ∼ / In this paper we have generalised the steady state, two dimensionalKida solution in the shearing sheet to allow for general vorticityprofiles and dust density distributions. The length scales are as-sumed to be small enough that the gas flow can be approximated asbeing incompressible. The solutions apply in the limit of very smalldust stopping time so that the dust density distribution is frozen-insuch that the system can be modelled as an incompressible fluidwith variable density. As this can be specified arbitrarily, a widerange of solutions is in principle possible. In particular, we gen-eralised the well known Kida solution so that it can incorporatea central dust density enhancement (see the polytropic models ofSection 3.5). However, the background disc is Keplerian only foran aspect ratio χ =
7. Other values of χ correspond to backgrounds χ ρ Figure 10.
The density excess above the background measured in units ofthe background density for the vortices with { α, β, ρ m } = { , , } (solidline), { α, β, ρ m } = { , , . } (dotted line), { α, β, ρ m } = { . , , . } (dashedline) and { α, β, ρ m } = { . , , . } (dashed-dot line). with density maxima or minima. Other solutions obtained numeri-cally were discussed in Section 5.3.These solutions, although two dimensional, can apply to horizontalplanes in an isothermal disc while being in hydrostatic equilibriumin the vertical direction. It is important to note that the Kida solutionis special in that the velocity inside the core is a linear functionof coordinates that results in in a constant period for circulatingaround streamlines, resulting in zero internal shear. This is not thecase outside the core or inside it for general vorticity distributions,with significant consequences for a discussion of stability.We have generalised the stability analysis applied to the core of c (cid:13) , 000–000 tability of vortices χ θ log γ (a) Streamline at (0 , . χ θ log γ (b) Streamline at (0 , . Figure 11.
The stability of the polytropic models in a non-Keplerian background. Cases with maximum density equal to three times the background value a)and six times the background value b). q ρ Figure 12.
The growth rate of the parametric instability of the polytropicmodel as a function of q expressed in units of Ω P for k z = the Kida vortex so that it can be extended to more general mod-els in Section 4. The core of a Kida vortex is special because thestability problem becomes separable for a specific choice of time-dependent wavenumber. This does not apply in general. However,we found that it is possible to look for modes localised on indi-vidual streamlines and that, from an Eulerian viewpoint, these canhave either a time-independent wavenumber or a wavenumber thatdepends on time. When the circulation period is not constant onstreamlines, the latter form ultimately increases linearly with timeand so cannot be associated with conventional exponentially grow-ing linear modes (see Section 4.8), a situation familiar in shearingbox calculations.The modes with time-independent wavenumber can be exponen- tially growing and turn out to yield those found in the Kida vortexcase. For only these can the analysis can be extended to apply tomore general vortices. We remark that when comparing the situ-ation to that of perturbations of an axisymmetric disc, the time-independent modes correspond to axisymmetric ( m =
0) modes,while the time-dependent wavenumber modes correspond to m (cid:44) χ < ∼ . Ω P . This is seen for the limitingcase of a point vortex and we can infer that this is a generic is-sue when considering dust spiralling inwards from the outer disctowards the vortex core. Models with a density excess can showmany narrow parametric instability bands though those with flattervorticity profiles show less strongly growing modes with growthrates ∼ . Ω P .We also considered the stability of the dust laden polytropic modelwith χ = k z =
0. These are e ff ectively those considered by Chang &Oishi (2010). This is possible as in this special case the vortex hasno internal shear. Indeed we found that modes could occur withgrowth rates, γ, up to 0 . Ω P . However, these modes would be af-fected by shear if an attempt is made to generalise them to othercases. Even cases with weak shear shown in Figure 3 correspond toshear ∼ . Ω P . For modes with putative growth rates ∼ . Ω P wemight accordingly expect temporary growth for 10 growth timesor a temporary amplification factor ∼ . We believe that a non-linear analysis is required to resolve the outcome in such a case.Note too as remarked in section 4.2.1 these modes are also likely to c (cid:13) , 000–000 Railton and Papaloizou be a ff ected by vertical stratification unless streamlines very close tothe midplane are considered.However, the results presented here taken together imply that dustparticles attracted from the outer disc to a vortex core with high as-pect ratio χ may well encounter parametric instabilities with char-acteristic growth rates of a few × − Ω P up to a maximum of0 . Ω P . This is the case even outside any high density core and soit is important to assess potential consequences for dust accumu-lation in vortices. The instabilities are parametric and local so thatthey can be inhibited by either dissipative e ff ects or e ff ects that dis-rupt the periodicity of the circulating motion. As the magnitude ofdissipative e ff ects is uncertain and for example dependent on theassumed location in a protoplanetary disc, we shall consider thesecond type of e ff ect.In this context we note that the accumulation of dust in vorticesmay occur rapidly (eg. Lyra et al. 2009; M´eheut et al. 2012b) suchthat the dust particle motion departs significantly from being peri-odic. This may be the case for the gas motion also (M´eheut et al.2012b). For these reasons, parametric instability may not have beenseen in simulations up to now. However, if we suppose parametricinstability is present, it is likely to lead to some form of low levelturbulence, This was indicated by work of both Lesur & Papaloizou(2010) and Lyra & Klahr (2011) who found that such instabilitiesdo not have a strongly disruptive e ff ect on large aspect ratio vor-tices produced by the subcritical baroclinic instability (SBI). Thusalthough parametric instability may act to cause a vortex to ulti-mately decay, it may be maintained if there is some mechanism togenerate it, such as the SBI or Rossby wave instability. For Kidavortices, strong, exponentially growing and potentially fast vortex-destroying instability only occurs on account of the saddle pointinstability which for which 3 / (cid:46) χ (cid:46) ff usion (eg. Lyra & Lin 2013). To make a crude estimate this wenote that the inflow rate for small particles driven by the pressuregradient is | v | ∼ |∇ P | /ρτ s (eg. Papaloizou & Terquem 2006).Supposing the vortex has a length scale L in the minor axis direc-tion, we estimate P ∼ ρ Ω P L and | v | ∼ Ω P L τ s . As the unstablemodes are local, the wavelength should be (cid:28) L . For the purposeof making crude estimates we adopt π/ | k | = L /
10. An estimateof the associated di ff usion coe ffi cient based on dimensional scal-ing is D = γ/ | k | . Balancing pressure driven inflow against dif-fusion we obtain |∇ ρ | / | ρ | = / L ρ ∼ | v | / D ∼ (10 π Ω P ) τ s ./ ( γ L ).Hence L ρ ∼ f L , where f ∼ ( γ/ Ω P ) / (100 π Ω P τ s ). Significant dustconcentrations become possible once f <
1. For γ = . Ω P , thisbecomes equivalent to Ω P τ s > ∼ − . Note that the characteristicinflow time is then 1 / ( Ω P τ s ). Although estimates of the above typeare highly uncertain they are made to indicate that the existence ofparametric instabilities of the type considered here does not neces-sarily prevent the possibility of dust accumulation in vortices. ACKNOWLEDGEMENTS
The authors would like to thank the anonymous reviewer for theirhelpful and insightful comments which lead to numerous improve-ments in the manuscript. ADR acknowledges support by STFC.
REFERENCES
Barge P., Sommeria J., 1995, A&A, 295, L1 Barranco J. A., Marcus P. S., 2005, ApJ, 623, 1157Benz W., 2000, Space Sci. Rev., 92, 279Bracco A., Chavanis P. H., Provenzale A., Spiegel E. A., 1999,Phys. Fluids, 11, 2280Cardoso O., Gluckmann B., Parcollet O., Tabeling P., 1996, Phys.Fluids, 8, 209Chamberlin T. C., 1900, J. Geo., 8, 58Chang P., Oishi J. S., 2010, ApJ, 721, 1593Chavanis P. H., 2000, A&A, 356, 1089Dominik C., Tielens A. G. G. M., 1997, ApJ, 480, 647Friedman J. L., Schutz B. F., 1978, ApJ, 221, 937Fromang S., Stone J. M., 2009, A&A, 507, 19Garaud P., Lin D. N. C., 2004, ApJ, 608, 1050Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125Goldreich P., Ward W. R., 1973, ApJ, 183, 1051Johansen A., Youdin A., Klahr H., 2009, ApJ, 697, 1269Kant I., 1755, Allgemeine Naturgeschichte und Theorie des Him-melsKida S., 1981, J. Phys. Soc. Jpn., 50, 3517Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1Lesur G., Papaloizou J. C. B., 2010, A&A, 513, 60Lifschitz A., Hameiri E., 1991, Phys. Fluids A, 3, 2644Lin M.-K., 2012, ApJ, 754, 21Lin M.-K., Papaloizou J. C. B., 2011, MNRAS, 415, 1426Lissauer J. J., 1993, ARA&A, 31, 129Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ,513, 805Lynden-Bell D., Ostriker J. P., 1967, MNRAS, 136, 293Lyra W., Johansen A., Zsom A., Klahr H., Piskunov N., 2009,A&A, 497, 869Lyra W., Klahr H., 2011, A&A, 527, 138Lyra W., Lin M.-K., 2013, ApJ, 775, 17M´eheut H., Casse F., Varniere P., Tagger M., 2010, A&A, 516, 31M´eheut H., Keppens R., Casse F., Benz W., 2012a, A&A, 542, 9M´eheut H., Meliani Z., Varniere P., Benz W., 2012b, A&A, 545,134M´eheut H., Varniere P., Casse F., Tagger M., 2010, in SF2A-2010:Proceedings of the Annual meeting of the French Society of As-tronomy and Astrophysics, Boissier S., Heydari-Malayeri M.,Samadi R., Valls-Gabaud D., eds., p. 265Papaloizou J. C. B., 2005, A&A, 432, 743Papaloizou J. C. B., Terquem C., 2006, Rep. Prog. Phys., 69, 119Petersen M. R., Julien K., Stewart G. R., 2007, ApJ, 658, 1236Pringle J. E., 1981, ARA&A, 19, 137Regev O., Umurhan O. M., 2008, A&A, 481, 21Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka., Vol. 26Sipp D., Jacquin L., 2000, Phys. Fluids, 12, 1740Swedenborg E., 1734, The Principia. Swedenborg InstituteWeidenschilling S. J., 1977, MNRAS, 180, 57Weizs¨acker C. F. V., 1944, Z. Astrophys., 22, 319Whipple F. L., 1972, in From Plasma to Planet, Elvius A., ed., p.211Whittaker E., Watson G., 1996, A Course of Modern Analysis,Cambridge Mathematical Library. Cambridge University Press c (cid:13)000