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 Godunovs
method.
4
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
• Roes 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
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 Powells 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 lt 1
• Complex Expensive for MHD...
• CTU (6 R-solves)
• Stable for CFL lt 1/2
• Relatively Simple...
• MUSCL-Hancock
• Stable for CFL lt 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
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/jst
one/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
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
BzB0sin(2px)
Sustained turbulence not possible in 2D
dissipation rate after saturation is sensitive to
numerical dissipation Code Test
Animation of angular velocity fluctuations
dVyVy1.5W0x
CTU with 3rd order reconstruction, 2562 grid,
bmin4000, 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
dVyVy1.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,
sustained turbulence
21
Stress Energy for ltBzgt ? 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