On Some Recent Developments in Numerical Methods for Relativistic MHD

1. On Some Recent Developments in Numerical Methods for Relativistic MHD as seen by an astrophysicist with some experience in computer simulations Serguei Komissarov School of Mathematics University of Leeds UK

2. Recent reviews in Living Reviews in Relativity (www.livingreviews.org): (i) Marti & Muller, 2003, “Numerical HD in Special Relativity (ii) Font, 2003, “Numerical HD in General Relativity”;

3. Optimistic plan of the talk • Conservation laws and hyperbolic waves. • Non-conservative (orthodox) and conservative (main stream) schools. • Causal and central numerical fluxes in conservative schemes. • Going higher order and adaptive. • Going multi-dimensional. • Keeping B divergence free. • 7. Going General Relativistic. • 8. Stiffness of magnetically-dominated MHD. • 9. Intermediate (trans-Alfvenic) shocks.

4. II. CONSERVATION LAWSAND HYPERBOLIC WAVES • Single conservation law U - conserved quantity, F- flux of U, S - source of U • System of conservation laws

5. 1D system of conservation laws with no source terms In many cases F is known as only an implicit function of U, f(U,F)=0 . In relativistic MHD the conversion of U into F involves solving a system of complex nonlinear algebraic equations numerically; computationally expensive ! • Non-conservative form of conservation laws Usually there exist auxiliary (primitive) variables, P, such that U and F are simple explicit functions of P. where

6. Continuous hyperbolic waves - Jacobean matrix Eigenvalue problem: • wavespeed of k-th mode - transported information Fast, Slow, Alfven, and Entropy modes in MHD

7. Shock waves - shock equations s – shock speed • Hyperbolic shocks: • As UraUl one has • (i) (Ur -Ul) ar ; • (ii) s alk. • e.g. Fast, Slow, Alfven, and • Entropy discontinuities in MHD continuous hyperbolic wave There exist other, non-hyperbolic shock solutions ! hyperbolic shock

8. III. NON-CONSERVATIVEAND CONSERVATIVE SCHEMES (a) Non-conservative school (orthodox) • Wilson (1972)a De Villiers & Hawley (2003), Anninos et al.(2005); • Finite-difference version of • Artificial viscosity (physically motivated dissipation) is utilised to construct • stable schemes; • (i) poor representation of shocks; • (ii) only low Lorentz factors (g<3); • “Why Ultra-Relativistic Numerical Hydrodynamics is Difficult” • by Norman & Winkler(1986); • Anninos et al.,(2005) : Go conservative!

9. (b) Conservative school where - exchange by the same amount of U between the neighbouring cells

10. IV. CAUSALAND CENTRAL NUMERICAL FLUXES (a) Causal (upwind) fluxes Utilize exact or approximate solutions for the evolution of the initial discontinuity at the cells interfaces (Riemann problems) to evaluate fluxes. Initial discontinuity Its resolution Implemented in the Relativistic MHD schemes by Komissarov (1999,2002,2004); Anton et al.(2005).

11. Linear Riemann Solver due to Roe (1980,1981) linearization at t = tn Riemann problem: Wave strengths: • a system of linear equations • for the wave strengths, a(k) Constant flux through the interface x = xi+1/2 : transported information wave strengths wave speeds

12. (b) Non-causal (central) fluxes ? Why not to try something simpler, like Well, this leads to instability. Why not to dump it with indiscriminate diffusion?! The modified equation where artificial diffusion This leads to the following numerical flux , where L is the highest wavespeed on the grid. Very high diffusion! • Lax:

13. Kurganov-Tadmor (KT): where l(k) are the local wavespeeds (Local Lax flux) • Harten, Lax & van Leer (HLL) : artificial diffusion where This makes some use of causality: i+1

14. Implemented in the Relativistic MHD schemes by: * HLL: Del Zanna & Bucciantini (2003), Gammie et al. (2003); Duez et al.(2005), Anton et al. (2005). * KT: Anton et al. (2005), Anninos et al.(2005) + Koide et al.(1996,1999) The central schemes are claimed to be as good as the causal ones ! Are they really?

15. 1D test simulations • (i) Stationary fast shock LRS – linear Riemann solver HLL

16. (ii) Stationary tangent discontinuity LRS HLL

17. (iii) Stationary slow shock LRS HLL

18. (iv). Fast moving slow shock LRS HLL

19. IV. GOING HIGHER ORDERAND ADAPTIVE • Fully causal fluxes provide better numerical representation of stationary • and slow moving shocks/discontinuities ( see also Mignone & Bodo, 2005) • However for fast moving moving shocks/discontinuities they give similar • results to central fluxes. (Lucas-Serrano et al. 2004, Anton et al. 2005). • How to improve the representation of shocks moving rapidly • across the grid ? • Use adaptive grids to increase resolution near shocks. • Falle & Komissarov (1996), Anninos et al. (2005); • (ii) Use sub-cell resolution to reduce numerical diffusion.

20. The nature of numerical diffusion

21. second order scheme first order scheme

22. first order accurate - piece-wise constant reconstruction; • second order accurate - piece-wise linear reconstruction; • Komissarov (1999,2002,2004), Gammie et al.(2003), • Anton et al. (2005). • third order accurate - piece-wise parabolic reconstruction; • Del Zanna & Bucciantini (2003), Duez et al. (2005). • In astrophysical simulations Del Zanna & Buucciantini are • forced to reduce their scheme to second order (oscillations at shocks) !? • THERE IS THE OPTIMUM?

23. V. GOING MULTI-DIMENSIONAL. F Ui,j F F F

24. VI. KEEPINGB DIVERGENCE FREE. • Differential equations • the evolution equation can • keep B divergence free ! • Difference equations may not have such a nice property. What do we do about this ? (i) Absolutely nothing. Treat the induction equation as all other conservation laws ( Koide et al. 1996,1999). Such schemes crash all too often! magnetic monopoles with charge density -“magnetostatic force”

25. (ii) Toth’s constrained transport. Use the “modified flux” F that is such a linear combination of normal fluxes at neighbouring interfaces that the “corner- -centred” numerical representation of divB is kept invariant during integration. Implemented in Gammie et al.(2003), and Duez et al.(2005)

26. (iii)Constrained Transport of Evans & Hawley. Use staggered grid (with B defined at the cell interfaces) and evolve magnetic fluxes through the cell interfaces using the electric field evaluated at the cell edges. This keeps the following “cell-centred” numerical representation of divB invariant Implemented in Komissarov (1999,2002,2004), de Villiers & Hawley (2003), Del Zanna et al.(2003), and Anton et al.(2005)

27. (iv) Diffusive cleaning Integrate this modified induction equation (not a conservation law ) - diffusion of div B Implemented in Anninos et al (2005)

28. (v) Telegraph cleaning by Dedner et al.(2002) Introduce new scalar variable, Y, additional evolution equation (for Y), and modify the induction equation as follows: conservation laws • the “telegraph equation” for div B

29. VII. GOING GENERAL RELATIVISTIC • GRMHD equations can also be written as conservation laws; - covariant continuity equation - continuity equation in partial derivatives • determinant of the • metric tensor • conservative form of the continuity equation. • t=x0/c and t=const defines a space-like • hyper-surface of space-time (absolute space) U F i Utilization of central fluxes is straightforward.

30. Riemann problems can be solved in the frame of the local Fiducial Observer (FIDO) • using Special Relativistic Riemann solvers. A FIDO is at rest in the absolute space • but generally is not at rest relative to the coordinate grid (Papadopoulos & Font 1998) - abg-representation of the metric form. Vector b is the grid velocity in FIDO’s frame Here we have got a Riemann problem with moving interface Implemented in Komissarov (2001,2004), Anton et al. (2005) FIDO’s frame coordinate grid

31. VIII. “STIFFNESS”OF MAGNETICALLY-DOMINATED MHD Magnetohydrodynamics Magnetodynamics This has 4 independent components.This has only 2 !

32. What to do if such magnetically-dominated regions do develop ? • Solve the equations of Magnetodynamics (e.g. Komissarov 2001,2004); • “Pump” new plasma in order to avoid running into the “danger zone”.

33. VIII. INTERMEDIATE SHOCKS This numerical solution of the relativistic Brio & Wu test problem is corrupted by the presence of non-physical compound wave which involves a non-evolutionary intermediate shock. Such shocks are known to pop up in non-relativistic MHD simulations. Brio & Wu (1988); Falle & Komissarov (2001); de Sterk & Poedts (2001); Torrilhon & Balsara (2004); Almost nothing is known about the relativistic intermediate shocks. How to avoid them? Use very high resolution. Torrilhon & Balsara (2004) fast rarefaction compound wave slow rarefaction intermediate shock fast rarefaction slow shock

34. Yes! This is finally over ! Thank you!

35. Stationary Contact LRS HLL

36. Slowly Moving Contact LRS HLL

37. Stationary Current Sheet LRS HLL

38. Fast Rarefaction Wave (LRS) 1st order; no diffusion. 1st order; LLF-type diffusion rarefaction shock 1st order 2nd order