3D Vortices in Stratified, Rotating, Shearing Protoplanetary Disks

1. Numerical Simulations of 3D Vortices in Stratified, Rotating, Shearing Protoplanetary Disks April 8, 2005 - I PAM Workshop I: Astrophysical Fluid Dynamics Philip Marcus – UC Berkeley Joe Barranco – KITP UCSB Xylar Asay-Davis – UC Berkeley Sushil Shetty – UC Berkeley

2. Observations of Protoplanetary Disks Mass 0.01 – 0.1 Msun Diameter ≈ 100 – 1000 AU Age ≤ 10 million years

3. Fluid Dynamics (along with radiation) Determines Transport Properties Angular momentum Dust Grains Migration of Planetesimals and Planets No observations of turbulence or fluid structures (yet)

4. Big PictureWhat Does a PPD Look Like? • Is it laminar or turbulent – is it both? • Is the mid-plane filled with vortices? • Are vortices long-lived or transient? • Is there an energy cascade? • Is it a 2D environment? • Is a Keplerian disk stable? • Can energy be extracted from the mean shear? • Are there pure hydro mechanisms at work?

5. Eddy Viscosity neddy * Replace the nonlinear, and difficult-to calculate advective term – (V¢r)V with a ficticious, linear, easy-to-calculate diffusion r¢neddyr V. * Set neddy = a cs H0 (Shakura & Sunyaev 1973) • < 1 because turbulent eddies are probably subsonic and not larger than a scale height H0 in extent. * Parameterize angular mom. transfer and/or rate of mass accretion in terms ofa. Origin of turbulence? Shear instabilities, convection, MHD instabilities … etc.

6. Does a work for Angular Momentum Transport (or anything else)? • Perhaps with heat transfer in non-rotating, (r¢ V), non-shearing flows such as thermal convection, e.g., Prandlt mixing-length theory. • Runs into problems when used for transporting vectors quantities and when there is competition among “advectively conserved” quantities. e.g., Spherical Couette Flow • Do not violate Fick’s Law. Angular momentum angular must be transported radially outward from a forming protostar; yet the disk’s ang. mom. increases with radius.

7. Formation of Planetesimals • 2 Competing Theories • Binary Agglomeration: Sticking vs. Disruption? • Gravitational Instability: Settling vs. Turbulence? Toomre criterion: vd < πGΣd/ΩK ≈ 10 cm/s mm grains  km planetesimals

8. Transport Implied by “Hot Jupiters” • Formed farther out in the disk where it was cooler, and then migrated to present location. • Type I Migration: Small protoplanet raises tides in the disk, which exert torques on the planet. • Type II Migration: Large protoplanet opens gap in the disk; both gap and protoplanet migrate inward on slow viscous timescale.

9. Vortices in Protoplanetary Disks? Anticyclonic shear leads to anticyclonic vortices

10. Vortices in Protoplanetary Disks? • Recipe for vortices: • Rapid rotation • Intense shear • Strong stratification

11. Assumptions for PPD Flow (being, or about to be, dropped) Disk is Cold: • Pressure small, so is radial pressure gradient • Gas is un-ionized and therefore not coupled to magnetic field • Sound speed cs¿ Keplerian velocity VK Cooling modeled as ¶ T /¶ t = L – T/tcool where is fixed.

12. Base Flow in Protoplanetary Disk Near balance between gravity and centrifugal force: Anticyclonic Shear No hydrostatic balance in the radial direction!

13. Base Flow in Protoplanetary Disk(and pre-computational bias) Vertical hydrostatic balance: Cool, thin disk:

14. Vortices in Protoplanetary Disks • Shear will tear a vortex apart unless: • Vortex rotates in same sense as shear. In PPD, vortices must be ANTICYCLONES. • Strength of the vortex is at least of the same order as the strength of the shear. • ωV ~ σk ~ Ωk • Ro ≡ ωz/2Ωk ~ 1

15. Vortices in Protoplanetary Disks • Velocity across vortex likely to be subsonic; otherwise shocks would rapidly dissipate kinetic energy of vortex: • V ~ σk L < cs • V ~ Ωk L < cs • ε ≡ V/cs~ (Ωk/cs)L < 1 • But from hydrostatic balance: cs ~ ΩkH • ε ≡ V/cs~ L/H < 1

16. PPD vs. GRS Wave speeds: cg /cs ¼ 1

17. Equations of Motion • Momentum: With Coriolis and Buoyancy • Divergence: Anelastic Approximation • Temperature: With Pressure-Volume Work

18. Two-Dimensional Approx. is not correct and misleading Too easy to make vortices due to limited freedom, conservation of w and inverse cascade of energy

19. Computational Method • Cartesian Domain: (r,f,z)  (x,y,z) • Valid when d≡ H0 /r0 << 1 and Dr/r0 << 1 • 3D Spectral Method • Horizontal basis functions: Fourier-Fourier basis for shearing box; Chebyshev otherwise • Vertical basis functions: Chebyshev polynomials for truncated domain OR Cotangent mapped Chebyshev functions for infinite domain • Parallelizes and scales

20. Spectral Methods Are Not L • Limited to Cartesian boxes • Limited to Fourier series • Limited to linear problems • Impractical without fast transforms

21. Spectral Methods • Beat 2nd order f.d. by ~ 4p per spatial dimension for 1% accuracy, (4p)2 for 0.1% • Often have diagonalizable elliptic operators • Non-dissipative & dispersive – must explicitly put in • Have derivative operators that commute • ¶2 / ¶ x2 = (¶ / ¶ x) (¶ / ¶ x) • Should not be used with discontinuities

22. Tests of the Codes • Energy, momentum, enstrophy and p.v. balance • Linear eigenmodes • Agrees with 2D solutions – Taylor Columns • Agrees with 3D asymptotics (equilibria) • Different codes (mapped, embedded and truncated) agree with each other • Gave unexpected results for which the algorithms were not “tuned”

23. Sliding box: shear s =(-3/2)W t=0 t=Lx/sd r d r Lx

24. Open Boundaries or 1 Domain • Mapping - 1 < z < 1 ! 0 < q < p • Arbitrary truncation to finite domain impose arbitrary boundary conditions • Embedding -Ldomain < -Lphysical < z < Lphysical < Ldomain

25. Mappings q q z = L cos(q) Chebyshev z = L cot(q) Others: Matsushima & Marcus JCP 1999 Rational Legendre for cylindrical coordinate at origin and harmonic at 1

26. Embedding Lose 1/3 of points in buffer regions Lembedding Buffer region Lphysical Domain of physical interest Lphysical Buffer region -Lembedding

27. No Free Lunch • Mappings lose 1/3 of the points for |z| > L advantage is (maybe) realistic b.c. • Truncations waste 1/3 of the points due to clustering at the artificial boundary • Embeddings lose 1/3 of the points at Lphysical < |z| < Lembedding

28. Internal Gravity Waves • Breaking: hri/ exp{-z2/2H02} (Gaussian) • Energy flux ~ hri V3 /2 • Huge energy source (locally makes slug flow of Keplerian, differential velocity) • Source of perturbations is oscillating vortices • Fills disk with gravity waves • Waves are neutrally stable • Gravity / z, stratification / z ! Brunt-Vaisalla frequency / z • For finite z domain (-L < z < L) there are wall-modes / exp{+z2/2H02}

29. Gravity Wave Damping 1 • ¶ Vx / ¶ t = L • ¶ Vy / ¶ t = L • ¶ Vz / ¶ t =L - ¶ h / ¶ z +g T/hTi • ¶ T / ¶ t =L- (d hSi / dz) (hTi/cp) Vz - T/tcool ¶2 Vz / ¶ t2 = - (d hSi / dz) (g/cp) Vz

30. Gravity Wave Damping 1 • ¶ Vx / ¶ t = L • ¶ Vy / ¶ t = L • ¶ Vz / ¶ t =L - ¶ h / ¶ z +g T/hTi – a1 ¶ T/¶ t • ¶ T / ¶ t =L- (d hSi / dz) (hTi/cp) Vz– a2 ¶ Vz/¶ t - T/tcool ¶2 Vz / ¶ t2 = - (d hSi / dz) (g/cp) Vz/(1 - a1a2) + [a1(d hSi / dz) (hTi/cp) - a2(g/hTi)]¶ Vz/¶ t

31. What does 2D mean?Columns or pancakes?

32. Tall Columnar Vortex

33. Vortex in the Midplane of PPD Ã  r = 2H0! Ã r0 = 4H0! Blue = Anticyclonic vorticity, Red = Cyclonic vorticity Ro = 0.3125

34. -z plane at r=r0 à  z = 8H0 ! Ro = 0.3125 Blue = Anticyclonic vorticity Red = Cyclonic vorticity à r0 = 4H0!

39. r-z plane at f=0 à  z = 8H0 ! Ro = 0.3125 Blue = Anticyclonic vorticity Red = Cyclonic vorticity à r = 2H0!

43. z = 1.5 H0 z = 1.0 H0 z = 0 H0

44. Spontaneous Formation of Off-Midplane Vortices

45. 3D Vortex in PPD

46.  v? High Pressure Coriolis Force z Equilibrium in Horizontal • Horizontal momentum equation: • For Ro · 1, Geostophic balance between gradient of pressure and the Coriolis force. • Anticyclones have high pressure centers . • For Ro À 1, low pressure centers.

47.  v? vz z Role of vz • In sub-adiabatic flow: rising cools the fluid while sinking warms it. • This in turn creates cold, heavy top lids and warm, buoyant bottom lids. • This balances the vertical pressure force (and has horizontal temperature gradients in accord with the thermal wind equation. • Numerical calculations show that after lids are created, vz! 0. • Magnitude of vz is set by dissipation rate, the faster of tcool or “advective cooling”.

48. Radial Cascade of Vortices • Internal gravity waves or small velocity perturbations with vz, or that create vz at outer edge of disk

49. Radial Cascade of Vortices • Anticlonic bands embedded in like-signed shear flow are unstable and break up into stable vorices • Cyclonic bands embedded in the opposite-signed Keplerian shear are stable • Energy provided by “step-functioning” the linearly stable background shear flow Marcus JFM1990 Coughlin & Marcus PRL1997

50. Radial Cascade of Vortices Create Radial bands of wz • ¶Vx / ¶ t = L • ¶ Vy / ¶ t = L • D wz / D t = (baroclinic) - (2W + s + wz) (vz /H) • ¶ Vz / ¶ t = L • ¶ T / ¶ t = L }