### 3.1 Numerical approaches

Maxwell’s equations can be simplified if the fluid is a perfect conductor. In this case is infinite and, to keep the current finite, the term must become zero, which results in a vanishing electric-field four-vector in the fluid rest frame, . This case corresponds to the ideal MHD condition. Under this assumption, the electric field measured by the Eulerian observer has components
and Maxwell’s equations reduce to the divergence-free condition plus the induction equation for the evolution of the magnetic field:
In the induction equation we use the notation .

For a fluid endowed with a magnetic field, the stress-energy tensor is the sum of that of the fluid and that of the electromagnetic field, , where is given by Equation (5) for a perfect fluid. On the other hand can be obtained from the Faraday tensor as follows:

which, in ideal MHD, can be rewritten as
where . The total stress-energy tensor is thus given by
with the definitions and .

#### 3.1.1 Nonconservative formulations

As happens for the equations of general-relativistic hydrodynamics discussed previously, available approaches to cast the corresponding magnetohydrodynamics equations in forms suitable for numerical work also follow to a great extent the basic distinction between conservative and nonconservative formulations. The latter, which are covered in this section, were first laid out by Wilson in the mid 1970s, who pioneered a new emerging field in numerical astrophysics, just as he had done a few years earlier for numerical Eulerian hydrodynamics. The original approach was followed by [429] and more recently by [86] and [20], who adopted an internal energy formulation with artificial viscosity for shock capturing and a method of characteristics for handling the magnetic fields.

The dynamic variables used in the formulation of [86], and , coincide with those reported in Equation (23), save for the definition of the total four-momentum , which includes the norm of the magnetic field four-vector , namely , with . With this choice, the “conservation” equations of mass density, momentum and energy are again written as a set of advection equations similar to Equations (20) – (22), except for the inclusion of the magnetic field terms, namely:

where is again the transport velocity defined in Section 2.1.2. As in the purely hydrodynamic case, the time component of the momentum is obtained from the corresponding normalization condition together with the algebraic relations between the magnetic-field components. Similarly, an artificial viscosity must be added for the numerical solution of the last two equations to increase the entropy of the fluid at shocks. In [86] such a viscosity prescription coincides with that used in the unmagnetized case [171], except for the modification of the definition of the inertia to include the magnetic energy.

In addition, a supplementary evolution equation is needed for the constrained transport magnetic field , the induction equation, which reads:

where must be subject to the divergence-free constraint, .

A similar internal-energy formulation was presented in [20]. There are, however, important subtle differences with respect to the approach of [86] in the particular form of the evolution equations, due to the different definition of the Lorentz (relativistic boost) factor employed by [20], namely . In addition, there are a couple of points particular to the formulation of [20], which are worth mentioning. First, in the momentum-evolution equation, the term is cast into projected and divergence components , which has significant numerical advantages for achieving stable formulations of sheared Alfvén waves. Second, the induction equation is rewritten in the form

where is a coefficient related to the largest characteristic speed in the flow multiplying the divergence cleanser function. Such a coefficient is included for purely numerical reasons in order to reduce discretization errors in time, either propagating any nonzero divergence off the computational domain or damping it (see [35] for tests of divergence-cleaning methods in classical MHD). This approach is adopted by [20] in order to avoid the implementation of the constrained transport method in their code (see Section 4.3) which, due to the need of an additional grid staggered with respect to the base grid, would much complicate the adaptive-mesh refinement (AMR) module of the code. Similar strategies are adopted in the AMR GRMHD code of [286].

#### 3.1.2 Conservative formulations

We turn now to the discussion of conservative formulations for the GRMHD equations. Very recently, a number of groups have proposed various such conservative approaches to formulating these equations, all aimed at solving them using methods specifically designed for hyperbolic systems of conservation (or balance) laws [149202011083662428626591]. The existing approaches all share the key ingredient of casting the GRMHD equations in first-order flux-conservative form, yet the choice of variables (state-vector, fluxes, etc) differs (albeit slightly) among most of them. Therefore, for the sake of clarity, we start here by discussing one such conservative formulation first, namely that of [24] (which coincides with the one obtained by [286] and is also employed in the codes of [265151]), leaving to the second part of this section the analysis of similarities and differences present in the additional formulations.

Following [24] the conservation equations for the energy-momentum tensor, together with the continuity equation and the equation for the evolution of the magnetic field, can be written as a first-order, flux-conservative, hyperbolic system equivalent to Equation (32). The state vector and the vector of fluxes of the GRMHD system of equations read

where, again, . Here, the conserved quantities for the continuity, Euler, and energy equations are now respectively defined as , , and . The corresponding vector of sources coincides with the one given by Equation (33) save for the use of the complete (fluid plus electromagnetic field) stress-energy tensor (the magnetic-field evolution equation is source-free). We note that expressions (57) and (58) contain components of the magnetic-field four-vector in the co-moving frame and of the magnetic-field three-vector measured by the Eulerian observer. The two are related by

The hyperbolic structure of the flux-vector Jacobian matrices corresponding to Equations (57)-(58), the knowledge of which is essential for building up a numerical code based on upwind HRSC schemes, is analyzed in [24] (see also [23]). In the classical MHD case, the wave structure has been analyzed by [61], which shows that there are seven physical waves: two Alfvén waves (with eigenvalues , and being the fluid and Alfvén speeds, respectively), two fast and two slow magnetosonic waves (, ), and one entropy wave (). The expressions for the Alfvén and magnetosonic speeds read

We note in passing the numerical analysis of the hyperbolic structure of the classical ideal MHD system carried out by [400], where nonuniform convergence to the true solution was found for a number of finite volume methods in coarse grid simulations.

The wave structure of the above GRMHD hyperbolic system is, understandably, more involved than that corresponding to the GRHD case. As mentioned before, Anile [15] (see also [16404]) studied the characteristic structure of these equations (eigenvalues, right and left-eigenvectors) in the (ten-dimensional) space of covariant variables , being the specific entropy. It was found that only the entropic waves and the Alfvén waves can be explicitly written in closed form. For the formulation of [24] discussed in this section, and along each coordinate direction, they are given by

On the other hand, the (fast and slow) magnetosonic waves are given by the numerical solution (through a root-finding procedure) of a quartic characteristic equation. Among the magnetosonic waves, the two solutions with maximum and minimum speeds are called fast magnetosonic waves , and the two solutions in between the slow magnetosonic waves . The seven waves can be ordered as follows

which defines the domain of dependence of the GRMHD system.

We note that, as a result of working with the unconstrained, system of equations, there appear superfluous eigenvalues associated with unphysical waves, which need to be removed (the entropy waves as well as the Alfvén waves appear as double roots). Any attempt to develop a numerical procedure based on the wave structure of the RMHD equations must, hence, remove these nonphysical waves (and the corresponding eigenvectors) from the wave decomposition. In the particular case of special-relativistic MHD, [199] and [198] eliminate the unphysical eigenvectors by enforcing the waves to preserve the values of the invariants and , as suggested by [15]. Correspondingly, in [34] the physical eigenvectors are selected by comparing with the equivalent expressions in the nonrelativistic limit.

We also note that, as in the classical case, the relativistic MHD equations have degenerate states in which two or more wavespeeds coincide. This breaks the strict hyperbolicity of the system and may affect numerical schemes for its solution. These degeneracies are analyzed in [199], which finds that, with respect to the fluid rest frame, the degeneracies in both classical and relativistic MHD are the same; either the slow and Alfvén waves have the same speed as the entropy wave when propagating perpendicularly to the magnetic field (Degeneracy I), or the slow or the fast wave (or both) have the same speed as the Alfvén wave when propagating in a direction aligned with the magnetic field (Degeneracy II). On the other hand, [24] derives a covariant characterization of such degenerate states, which can be checked in any reference frame. In addition, [23] also works out a single set of right and left eigenvectors, which are regular and span a complete basis in any physical state, including degenerate states. Such renormalization procedure can be seen as a relativistic generalization of the work performed in [61] in classical MHD.

Finally, it is worth pointing out that major advances in the physical understanding of the wave structure of the GRMHD equations have been possible in recent years owing to the remarkable derivation of the exact solution of the Riemann problem in special relativistic MHD [340150].

As mentioned previously, a number of conservative formulations for the GRMHD equations have been proposed recently. Most of them show some differences with respect to the formulation just discussed, and some similarities too. Let us now consider these issues.

First, in the formulations of Gammie et al. [149] and Komissarov [201], the choice of conserved variables corresponding to the momentum and energy equations follows directly from the conservation of energy momentum, which is written with the free index down, i.e.,  . The continuity equation, however, remains identical to the expression used in [24]. Therefore, the vector of conserved variables in [201149] reads

The choice of for the energy-momentum equations is not capricious, but rather motivated by the reduction on the number of source terms in the corresponding evolution equations for the case of spacetime metrics with symmetries (ignorable coordinates ; see also [309] for a similar approach in the purely hydrodynamic case). This is not the situation in the formulation of [24] in which the energy-momentum equations are obtained, in the spirit of the 3+1 splitting, from projections of the conservation law, namely and .

Regarding the induction equation, both [149] and [201] use the following expression

which differs from the expression used by [24] both from the use of and the presence of and accompanying the time and spatial derivatives, respectively, in the latter case.

On the other hand, in the approach followed by Shibata and Sekiguchi [366], the evolution equations for the energy-momentum equations are also derived as in [24], from projections of the conservation law. However, the choice of conserved variables is slightly different, mainly due to the presence in all equations of a common factor , being a spacetime conformal factor, motivated by the formulation used in the numerical code for solving Einstein’s equations, namely BSSN [36339]. (We note that such factors are also present in the formulation used by Shibata in the purely hydrodynamic case; see, e.g., [355].) Another important difference with [24] is the choice of the transport velocity throughout (instead of ) and the fact that in the final form of the energy equation the continuity equation has not been subtracted. Additionally, the induction equation also adopts a different form as it is written for a three-magnetic field given by , with .

Correspondingly, in the conservative formulation of Anninos et al. [20] the main differences appear in the energy-momentum equations, leaving the continuity equation unaltered with respect to the nonconservative formulation discussed in Section 3.1.1. Regarding the momentum and energy-conservation equations the choice of variables in [20] is different to those of [24] and [366], as they are components of the energy-momentum tensor and not projections of the conservation law. But the choice is also different to that in [149] and [201], as the indices of the energy-momentum tensor are both chosen to be contravariant (instead of and ). More precisely, the corresponding conserved variables are and . Another important point to mention regarding the set of equations of this conservative formulation (and also in their nonconservative internal-energy formulation) is the presence of factors of in some terms, which have not been absorbed by the definition of the electromagnetic tensor (as in most of the other approaches). Finally, regarding the induction equation, we note that it is also cast in fully conservative form

with and the same coefficient as in Equation (56).

Next, the formulation of Duez et al. [108] also shares some similarities and has some differences with all previous approaches discussed so far. The transport velocity is used throughout. The choice of conserved variable for the continuity equation, which the authors name , is the same as in [20], being the equation also cast as a source-free advection equation. As for the momentum equation the choice of conserved variable coincides with that of [149201], i.e., direct components of the energy-momentum tensor, namely . However, for the remaining part of the energy-momentum conservation law, which leads to the energy equation, Duez et al. [108] do not follow [149201] and, instead of using , they choose , the same variable as [24], also subtracting the continuity equation, namely . Finally, the induction equation is written as in Equation (67), but without the divergence cleanser coefficient ().

Finally, in the formulation of Del Zanna et al. [91] an effort is made in the choice of the state vector so that the GRMHD equations keep a somewhat close relation to classical MHD and to simplify the expression of the source terms. As in Antón et al. [24] the same definition for the three-velocity is used, along with the same variable for the relativistic density (densitized with a factor ). For the Euler equation, components of the stress-energy tensor are used, as in the approaches of [149201108], which reduces the complexity of the corresponding source term. The most important difference with respect to other formulations is in the choice of conserved variable for the energy equation, for which [91] take . The continuity equation is not subtracted from the energy equation. As far as the induction equation is concerned, the same conservative expression used by [149201] is employed.