A Godunov-type scheme for relativistic magnetohydrodynamics

S. S. Komissarov 0 1 0 Astrospace Centre, Lebedev Physical Institute , Leninsky Prospect 53, Moscow B-333, 177924 Russia 1 Department of Applied Mathematics, The University of Leeds , Leeds LS2 9JT A B S T R A C T The basic properties of relativistic magnetohydrodynamics, as a hyperbolic system of quasilinear conservation laws, are discussed. These are then used to develop a multidimensional Godunov-type numerical scheme that enforces the magnetic flux conservation. This scheme is based on linear Riemann solvers and has second-order accuracy in smooth regions. The results of thorough test calculations demonstrate that the scheme is robust and can cope with truly ultrarelativistic problems. I N T R O D U C T I O N It is now well recognized that compact objects such as neutron stars and black holes can drive a wide variety of relativistic flows. Accretion flows on to supermassive black holes are generally regarded as being the main power source in active galactic nuclei (Zeldovich 1964; Begelman, Blandford & Rees 1984). Radio observations also show that the energy that is released frequently appears in the form of highly collimated relativistic jets (e.g. Bridle 1996). These jets are sources of synchrotron emission, which indicates the presence of significant magnetic fields. Various models for the production of these relativistic jets have been proposed (e.g. Wiita 1991). Although it is extremely difficult to test these models in observations, the most plausible ones include strong magnetic fields as a key ingredient. Accretion on to neutron stars and black holes of stellar mass is a basic feature of most models of compact galactic X-ray sources (Lewin, van Paradijs & van den Heuvel 1995). For a long time the precessing jets of SS433 (Milgrom 1979) were the only reliable example of relativistic outflows from such objects. The recent discovery of galactic superluminal radio sources shows that such outflows may also be quite common (Mirabel & Rodrguez 1994; Tingay et al. 1995). Highly relativistic blast waves are believed to be responsible for the gamma-ray bursts (Fenimore et al. 1992; Meszaros & Rees 1992). These developments in astrophysics explain the growing interest in relativistic gas dynamics and magnetohydrodynamics during the last couple of decades. Numerous attempts have been made to develop numerical schemes that would allow us to study such flows (e.g. Wilson 1972; Centrella & Wilson 1984; Dubal 1991; van Putten 1993, 1995). These early schemes could handle only weakly or moderately relativistic flows with Lorentz factors not much higher than 2. They were very useful in some applications, but also revealed certain problems in numerical relativity. One of these difficulties relates to the use of artificial viscosity required to stabilize shocks (Norman & Winkler 1986). Indeed, in relativistic gas dynamics dissipative processes not only contribute to the fluxes of conserved quantities but also to their space densities (e.g. Landau & Lifshitz 1959). If one ignores these modifications to the conserved variables then significant errors can arise in the values of primitive variables (velocity, density, etc.) computed from the conservative ones. On the other hand, the proper dissipative equations are very complicated and include mixed space and time derivatives. This seems to be the main reason for the failure of the schemes that require significant artificial viscosity. Recently, Koide, Nishikawa & Mutel (1996) and Nishikawa et al. (1997) have presented the results of multidimensional numerical simulations of a relativistic jet with an initial Lorentz factor as high as u0 4:56 using a scheme of LaxWendroff type. This is certainly an achievement. Unfortunately, the paper describing their scheme and the test simulations has not yet been published, which makes it difficult to comment on these results. A similar effect probably explains why schemes that rely heavily upon smoothing operators to stabilize shocks ( e.g. van Putten 1993, 1995) cannot handle problems with high Lorentz factors. Independent smoothing of the energy density and momentum density, which have very close magnitudes for ultrarelativistic velocities, is bound to introduce significant errors in the subsequently computed primitive variables and may even produce a case where no solution for primitive variables exists. This suggests (see also the discussion in Norman & Winkler 1986) the use of Godunov-type shock-capturing schemes, because they do not require large artificial viscosity or smoothing operators in order to handle shocks. Recent developments in this direction have been very successful (Eulderink & Melemma 1994; Font et al. 1994; Duncan & Hughes 1994; Mart & Muller 1996; Falle & Komissarov 1996). The test simulations presented in Mart & Muller (1996) and Falle & Komissarov (1996) included flows with Lorentz factor up to 200, and the simulations of relativistic jets presented in Mart et al. (1997) and Komissarov & Falle (1998) are as good as the equivalent classical simulations. In fact, it is the amount of artificial dissipation introduced that matters and we have to emphasize that a small amount of artificial dissipation is useful even in shock-capturing schemes (Falle & Komissarov 1996, see also Sections 4.5 and 5.4). The next step is the development of the shock-capturing scheme for relativistic magnetohydrodynamics (RMHD), since in many cases magnetic fields are believed to be dynamically important. Here, we describe an attempt to develop such a scheme. B A S I C E Q U AT I O N S Evolution equations of RMHD For a uniform chemical composition, the equations of RMHD can be written in the form of the covariant conservation laws (Dixon 1978; Anile 1989): Tab w b2uaub is the energy-momentum tensor. w, p and ua u0; u are the fluid enthalpy, pressure and four-velocity and gab is the metric tensor. Fab baub bbua where Fgd is the electromagnetic field tensor and habgd is the Levi Civita alternating tensor. In the fluid frame ba 0; B where B is the the usual three-vector of magnetic field. In an arbitrary frame the three-vectors of magnetic and electric fields are related to ba by B Fi0 bu0 ub0; b0 Bu; b B b0u=u0: As in classical magnetohydrodynamics (MHD), the second Maxwell equation is only used to compute the electric current: Jb aFba: 10 The only difference is that Jb now depends on the time derivatives of the electric field (the displacement current). The space component of (2) is the same as the classical induction equation: B = b B 0; 11 t where b is the three-velocity of plasma. Note that we are using units such that the factor 4p, the (constant) magnetic permeability and the speed of light do not appear in the equations (the speed of light is unity). We also assume that greek indices run from 0 to 3, latin ones run from 1 to 3 and the signature of spacetime is . Constraints F00 ; 0, which means that the time component of (2) is not in fact an evolution (...truncated)