Recent advances in Astrophysical MHD

1. Recent advances in Astrophysical MHD Jim Stone Department of Astrophysical Sciences & PACM Princeton University, USA Recent collaborators: Tom Gardiner (Cray Research) John Hawley (UVa) Peter Teuben (UMd)

2. Eagle nebula (M16) Numerical methods for MHD are crucial for understanding the dynamics of astrophysical plasmas.

3. Modern schemes solve the equations of ideal MHD in conservative form CD Many schemes are possible: central schemes WENO schemes TVD schemes We have adopted finite-volume techniques using Godunov’s method.

4. rEV Bx By Basic Algorithm: Discretization Scalars and velocity at cell centers Magnetic field at cell faces Cell-centered quantities volume-averaged Face centered quantities area-averaged Area averaging is the natural discretization for the magnetic field.

5. Key element of Godunov method is Riemann solver • For pure hydrodynamics of ideal gases, exact/efficient nonlinear • Riemann solvers are possible. • In MHD, nonlinear Riemann solvers are complex because: • There are 3 wave families in MHD – 7 characteristics • In some circumstances, 2 of the 3 waves can be degenerate (e.g. VAlfven = Vslow ) Equations of MHD are not strictly hyperbolic (Brio & Wu, Zachary & Colella) Thus, in practice, MHD Godunov schemes use approximate and/or linearized Riemann solvers.

6. Many possible approximations are possible: • Roe’s method – keeps all 7 characteristics, but treats each as a simple wave. • Harten-Lax-van Leer (HLLE) method – keeps only largest and smallest characteristics, averages intermediate states in-between. • HLLD method – Adds entropy and Alfven wave back into HLLE method, giving four intermediate states. Good resolution of all waves Requires characteristic decomposition in conserved variables Expensive and difficult to add new physics Does not preserve positivity Very simple and efficient Guarantees positivity Very diffusive for contact discontinuities Reasonably simple and efficient Guarantees positivity Better resolution of contact discontinuities

7. D 2 F D . D Keeping B = 0 • Do nothing. Assume errors remain small and bound. • Evolve B using vector potential defined through B= • Remove solenoidal part of B using “flux-cleaning”. That is, set B  B – where + B = 0 • Use Powell’s “8-wave solver” (adds div(B) source terms) • Evolve integral form of induction equation so as to conserve magnetic flux (constrained transport). Requires taking second difference numerically to compute Lorentz force fD Requires solving elliptic PDE every timestep – expensive May smooth discontinuities in B Gives wrong jump conditions for some shock problems Requires staggered grid for B (although see Toth 2000)

8. The CT Algorithm • Finite Volume / Godunov algorithm gives E-field at face centers. • “CT Algorithm” defines E-field at grid cell corners. • Arithmetic averaging: 2D plane-parallel flow does not reduce to equivalent 1D problem • Algorithms which reconstruct E-field at corner are superior Gardiner & Stone 2005

9. Simple advection tests demonstrate differences Field Loop Advection (b = 106): MUSCL - Hancock Arithmetic average Gardiner & Stone 2005 (Balsara & Spicer 1999)

10. Which Multidimensional Algorithm? CornerTransportUpwind [Colella 1991] (12 R-solves) • Optimally Stable, CFL < 1 • Complex & Expensive forMHD... CTU (6 R-solves) • Stable for CFL < 1/2 • Relatively Simple... MUSCL-Hancock • Stable for CFL < 1/2 • Very Simple, but diffusive...

11. Verification: Linear Wave Convergence(2N x N x N) Grid

12. Validation: Hydro RT instability 2562 x 512 grid Random perturbations Isosurface and slices of density

13. experiment Dimonte et al, Phys. Fluids (2004) have used time evolution of rising “bubbles” and falling “spikes” from experiments to validate hydro codes: Asymptotic slope too small in ALL codes by about factor 2 Probably because of mixing at grid scale Comparison with expts. using miscible fluids much better

14. Application: 3D MHD RT instability 2562 x 512 grid Random perturbations Isosurface and slices of density B = (Bx, 0, 0) lcrit = Lx/2

15. Codes are publicly available • Latest project is funded by NSF ITR; source code public. • Code, documentation, and training material posted on web. • 1D, 2D, and 3D versions are/will be available. Download a copy from www.astro.princeton.edu/~jstone/athena.html • Current status: • 1D version publicly available • 2D version publicly available • 3D version will be released in ~1 month

16. EXAMPLE APPLICATION: MHD of Accretion disks e.g., mass transfer in a close binary

17. W a r-3/2 L a r1/2 If accreting plasma has any angular momentum, it will form a rotationally supported disk Profiles of specific angular momentum (L) and orbital frequency (W) for Keplerian disk: • Thereafter, accretion can only occur if angular momentum is transported outwards. • microscopic viscosity too small • anomalous (turbulent?) viscosity required • MHD turbulence driven by magnetorotational instability (MRI) dominates

18. Start from a vertical field with zero net flux: Bz=B0sin(2px) Sustained turbulence not possible in 2D – dissipation rate after saturation is sensitive to numerical dissipation: Code Test Animation of angular velocity fluctuations: dVy=Vy+1.5W0x CTU with 3rd order reconstruction, 2562 grid, bmin=4000, orbits 2-10

19. Magnetic Energy Evolution Numerical dissipation is ~1.5 times smaller with CTU & 3rd order reconstruction than Zeus.

20. 3D MRI Animation of angular velocity fluctuations: dVy=Vy+1.5W0x Initial Field Geometry is Uniform By CTU with 3rd order reconstruction, 128 x 256 x 128 Grid bmin= 100, orbits 4-20 In 3D, sustainedturbulence

21. Stress & Energy for <Bz>  0 No qualitative difference with ZEUS results (Hawley, Gammie, & Balbus 1995)

22. Vertical structure of stratified disks Now using static nested-grids to refine midplane of thin disks with cooling in shearing box Next step towards nested-grid global models Density Angular momentum fluctuations

23. Future Extensions to Algorithm • Curvilinear coordinates • Nested and adaptive grids (already implemented) • full-transport radiation hydrodynamics (w. Sekora) • non-ideal MHD (w. Lemaster) • special relativistic MHD (w. MacFayden) Future Applications Star formation Global models of accretion disks