Numerical Relativistic Magnetohydrodynamics
Compact objects like black holes and neutron stars interacting with the relativistic plasma in the sorrounding region are believed to be responsible for many high-energy phenomena in astrophysics. Active galactic nuclei, relativistic jets from compact objects, gamma-ray bursts, extraction of rotational energy from the ergosphere of a Kerr black hole represent physical examples where the interplay between strong magnetic and strong gravitational fields is the crucial factor for understanding their dynamics. The theoretical framework in which these mechanisms are treated is that of general relativistic magnetohydrodynamics (GRMHD). Since relativistic magnetized flows are often associated with the formation of strong shocks at high Lorentz factors, the numerical solution of the GRMHD equations is very efficiently performed through shock-capturing Godunov-type methods, that are based on the conservative formulation of the equations.
We can see how the conservative form of the GRMHD equations are obtained by starting from the expression of the energy momentum tensor, which is given by the joint contribution of matter and of electromagnetic fields
The Faraday electromagnetic tensor F satisfies Maxwell's equations
While the four-current can be decomposed according to the fluid velocity u
As in classical MHD, the current must be a derived quantity, thus we still need to impose an extra condition on F, and this is represented by Ohm's law for a perfectly conducting fluid, that takes the form:
In other words, we are assuming that the freely moving charges can screen any electric field that may arise locally. Therefore, the electromagnetic energy-momentum tensor and the Faraday tensor become
The system of the GRMHD equations is then given by the continuity equation, the energy-momentum equation, and by the induction equation for the magnetic field, augmented by the divergence free condition.
However, if we want to perform numerical simulations, this covariant representation of the fundamental equations is not very useful, and an unambiguous separation between time and space coordinates (with respect to a some given observer) becomes necessary. This can be done by resorting to a consolidated approach in numerical relativity, namely the so called 3+1 formalism that imports into the relativistic language some basic intuitions of the Newtonian framework.
In the 3+1 formalism the space-time is foliated into non-intersecting space-like hyper-surfaces. The unit normal to defines the four-velocity of a special observer, referred to as the Eulerian observer:
The generic metric (with Β=shift vector) of the spacetime is then split into a temporal and a spatial part
Similarly, any tensor can be split into its temporal and spatial components
By replacing the decomposed tensors into the covariant form of the equations, the conservative form of the GRMHD equations becomes available
where
In compact form the GRMHD equations read
There are a number of peculiar features that distinguish these equations with respect to the purely hydrodynamic ones, when magnetic fields are absent. The eigenstructure is much more complex with 7 modes as in classical MHD (2 Alfvens, 4 magnetosonics, 1 entropy wave). As in classical MHD, there are degenerate states in which two or more speeds coincide. The curl nature of the induction equation and of the diverge free condition must be maintained in the numerical scheme. Discontinuous solutions satisfy jump relations just for the tangential component of the magnetic field, while the normal field component is continuous. Moreover, no closed analytic expressions between physical and conserved variables exist:
The problem of the conversion from conservative to primitive variables can be reduced to the solution of a 2x2 set of non-linear equations, to be solved with standard root-finding procedures.
where 
and 
depends on the Equation of State used
Once x and y have been computed, the primitive variables are easly found from the relations:
The strategy for obtaining a high order shock capturing numerical scheme for the GRMHD equations in conservation form proceeds along three mains steps. Firstly, the conserved variable U must be reconstructed to build the flux F at the interfaces between adjacent cells. Secondly, a proper approximate Riemann solver must be chosen for the computation of the fluxes. Finally, a high order computation of the flux is needed.
We can obtain high order reconstruction, by adopting routines that avoid decomposition into characteristic waves, and use non-oscillatory interpolation techniques. There are several algorithms available in the literature.
- ENO3 for third order original ENO method (Harten et al. 1997)
- CENO3 for the Convex-ENO scheme (Liu & Osher 1998)
- WENO5 for the Weighted-ENO fifth order scheme (Jiang & Shu 1996)
- MP5 for the Monotonicity Preserving scheme by Suresh & Huynh (1997)
WENO5: perform a weighted combination of 3-point substencils
MP5: Correct it with a non-linear filter based on limiting algorithms if spurious oscillations are found. Accuracy is preserved near smooth extrema and only degrades to first order accuracy at discontinuities.
A very convenient flux formula is provided by the HLL Riemann solver
The 2 fast and 2 slow magnetosonic waves are the roots of a 4-th order polinomial:
which admits an approximate reduction (corresponding to the degenerate case of normal propagation 
) to a quadratic equation for the computation of the 2 fast magnetosonics.
Finally, we can provide a high order approximation of the numerical flux function from the point values at the interface between adjacent cells
A standard Taylor expansion allows to obtain...
In order to guarantee the divergence free condition for the magnetic field, we use the Constraint Transport method:
In 2D Cartesian coordinates this means:
Some Numerical tests:
The exact solution of a typical Riemann problem can be compared with the numerical solution at a given time
A strong shock propagating into a magnetized medium (cylindrical explosion)
Initial conditions: at the center of the grid place a circle with 
while outside 
. There is a uniform magnetic field 
and zero velocity everywhere. Performed with MP5 at third order.
Radial accretion onto a Schwarzschild black hole. The numerical solution is compared to the stationary one.
Stationary and axisymmetric geometrically thick disc solution around a Kerr black hole
Essential bibliography:
Marti J.M., Mueller E., http://www.livingreviews.org/lrr-2003-7/
Font J.A., http://www.livingreviews.org/lrr-2003-4/
Anton et al., ApJ, 637, 296, 2006
Del Zanna L., Bucciantini N., Londrillo P., 2003, A&A 400, 397
Londrillo P., Del Zanna L., 2004, J. Comput. Phys., 195, 17
Del Zanna L., Zanotti O., Bucciantini, N., Londrillo, P., A&A, 473, 11, 2007
|Webmaster: Claudio Marchesini|