Go to previous page Go up Go to next page

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 F μνuν must become zero, which results in a vanishing electric-field four-vector in the fluid rest frame, eμ = 0. This case corresponds to the ideal MHD condition. Under this assumption, the electric field measured by the Eulerian observer has components
0 i 0ijk E = 0, E = − αη vjBk, (47 )
and Maxwell’s equations ∇ ∗F μν = 0 ν reduce to the divergence-free condition plus the induction equation for the evolution of the magnetic field:
√ -- ∂-(-γBi-) -1---∂- √ -- i -1---∂- √ -- i j j i ∂xi = 0, √ γ∂x0 ( γB ) = √ γ∂xj { γ[α &tidle;vB − α&tidle;v B ]}. (48 )
In the induction equation we use the notation &tidle;vi = vi − βi∕ α.

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, T μν = Tμν + Tμν Fluid EM, where T μν Fluid is given by Equation (5View Equation) for a perfect fluid. On the other hand μν TEM can be obtained from the Faraday tensor as follows:

μν μλ ν 1- μν λδ T EM = F Fλ − 4 g F F λδ, (49 )
which, in ideal MHD, can be rewritten as
( 1 ) TEμνM = uμuν + --gμν b2 − bμbν, (50 ) 2
where b2 = bνbν. The total stress-energy tensor is thus given by
μν ∗ μ ν ∗ μν μ ν T = ρh u u + p g − b b , (51 )
with the definitions p∗ = p + b2∕2 and h ∗ = h + b2∕ ρ.

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 [429Jump To The Next Citation Point] and more recently by [86Jump To The Next Citation Point] and [20Jump To The Next Citation Point], 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 [86Jump To The Next Citation Point], D, Sμ and E, coincide with those reported in Equation (23View Equation), save for the definition of the total four-momentum Sμ, which includes the norm of the magnetic field four-vector bμ, namely S μ = (ρh + b2)W uμ, with W = αu0. With this choice, the “conservation” equations of mass density, momentum and energy are again written as a set of advection equations similar to Equations (20View Equation) – (22View Equation), except for the inclusion of the magnetic field terms, namely:

∂D-- -1---∂-- √ -- i ∂t + √ γ-∂xi (D γV ) = 0, (52 ) ∂-(S − αb bt) + √1--∂--√ γ(S V i − αb bi) ∂t j j γ∂xi j j ( 2 ) ( ) + α -∂-- p + b- + 1- S-αSβ − αb βbα -∂--gαβ = 0, (53 ) ∂xj 2 2 S0 ∂xj ∂E 1 ∂ √-- ∂W p ∂ √ -- ----+ √-----i( γEV i) + p----+ √-----i( γW Vi) = 0, (54 ) ∂t γ∂x ∂t γ ∂x

where i V 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 [86Jump To The Next Citation Point] such a viscosity prescription coincides with that used in the unmagnetized case [171Jump To The Next Citation Point], 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 √ --- ℬi = 4 π√ γW (bi − V ib0), the induction equation, which reads:

∂ ℬi ∂ ----− ---j(Viℬj − ℬiV j) = 0, (55 ) ∂t ∂x
where ℬi must be subject to the divergence-free constraint, j ∂-ℬ- = 0 ∂xj.

A similar internal-energy formulation was presented in [20Jump To The Next Citation Point]. There are, however, important subtle differences with respect to the approach of [86Jump To The Next Citation Point] in the particular form of the evolution equations, due to the different definition of the Lorentz (relativistic boost) factor employed by [20Jump To The Next Citation Point], namely √ --- 0 √ -- 0 W = − gu = α γu. In addition, there are a couple of points particular to the formulation of [20Jump To The Next Citation Point], which are worth mentioning. First, in the momentum-evolution equation, the term -∂-√ --- i ∂xi( − gbjb ) is cast into projected and divergence components √--- √ --- bi∂∂xi( − gbj) + − gbj ∂∂xibi, which has significant numerical advantages for achieving stable formulations of sheared Alfvén waves. Second, the induction equation is rewritten in the form

( ) ∂ℬi- -∂-- i j j∂V-i -∂-- ∂ℬi- ∂t + ∂xj (ℬ V ) = ℬ ∂xj + η∂xj ∂xi , (56 )
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 ∇ ⋅ B 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 [20Jump To The Next Citation Point] 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 [286Jump To The Next Citation Point].

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 [149Jump To The Next Citation Point20Jump To The Next Citation Point201Jump To The Next Citation Point108Jump To The Next Citation Point366Jump To The Next Citation Point24Jump To The Next Citation Point286Jump To The Next Citation Point265Jump To The Next Citation Point91Jump To The Next Citation Point]. 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 [24Jump To The Next Citation Point] (which coincides with the one obtained by [286Jump To The Next Citation Point] and is also employed in the codes of [265Jump To The Next Citation Point151Jump To The Next Citation Point]), leaving to the second part of this section the analysis of similarities and differences present in the additional formulations.

Following [24Jump To The Next Citation Point] 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 (32View Equation). The state vector and the vector of fluxes of the GRMHD system of equations read

U (w ) = (D, S ,τ,Bk ), (57 ) j Fi(w ) = (D &tidle;vi,Sj&tidle;vi + p∗δij − bjBi∕W, i ∗ i 0 i i k k i τ&tidle;v + p v − αb B ∕W, &tidle;v B − &tidle;v B ), (58 )
where, again, &tidle;vi = vi − βi∕ α. Here, the conserved quantities for the continuity, Euler, and energy equations are now respectively defined as D = ρW, Sj = ρh ∗W 2vj − αb0bj, and ∗ 2 ∗ 2 0 2 τ = ρh W − p − α (b ) − D. The corresponding vector of sources coincides with the one given by Equation (33View Equation) 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 (57View Equation) and (58View Equation) 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
i i 0 i b0 = W-B--vi, bi = B--+--αb-u-. (59 ) α W

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

∘--2- v = Bx-, (60 ) a ρ (| ┌│ (-------------------)2--------)| 2 1-{ 2 B2x +-B2y-+-B2z │∘ 2 B2x +-B2y-+-B2z 22} vf,s = 2 |(cs + ρ ± cs + ρ − 4vacs|) . (61 )

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 [15Jump To The Next Citation Point] (see also [16404]) studied the characteristic structure of these equations (eigenvalues, right and left-eigenvectors) in the (ten-dimensional) space of covariant variables (uμ,bμ,p,s), s 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 [24Jump To The Next Citation Point] discussed in this section, and along each coordinate direction, they are given by

i i λe = αv −√β-,------ (62 ) bi ± ρh + B2ui λa ± = -0---√--------2-0 . (63 ) b ± ρh + B u

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 λf±, and the two solutions in between the slow magnetosonic waves λs±. The seven waves can be ordered as follows

λf− ≤ λa − ≤ λs− ≤ λe ≤ λs+ ≤ λa+ ≤ λf+, (64 )
which defines the domain of dependence of the GRMHD system.

We note that, as a result of working with the unconstrained, 10 × 10 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, [199Jump To The Next Citation Point] and [198Jump To The Next Citation Point] eliminate the unphysical eigenvectors by enforcing the waves to preserve the values of the invariants uμuμ = − 1 and u μbμ = 0, as suggested by [15]. Correspondingly, in [34Jump To The Next Citation Point] 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 [199Jump To The Next Citation Point], 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, [24Jump To The Next Citation Point] derives a covariant characterization of such degenerate states, which can be checked in any reference frame. In addition, [23Jump To The Next Citation Point] 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 [340Jump To The Next Citation Point150Jump To The Next Citation Point].

3.1.2.1 Additional conservative approaches
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. [149Jump To The Next Citation Point] and Komissarov [201Jump To The Next Citation Point], 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.,  ∇ T μ = 0 μ ν. The continuity equation, however, remains identical to the expression used in [24Jump To The Next Citation Point]. Therefore, the vector of conserved variables in [201Jump To The Next Citation Point149Jump To The Next Citation Point] reads

√ --- 0 0 0 i U = − g(ρu ,T i,T0,B ). (65 )

The choice of 0 T μ 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 x μ; see also [309] for a similar approach in the purely hydrodynamic case). This is not the situation in the formulation of [24Jump To The Next Citation Point] in which the energy-momentum equations are obtained, in the spirit of the 3+1 splitting, from projections of the conservation law, namely ν μ γi ∇ μT ν = 0 and ν μ n ∇ μTν = 0.

Regarding the induction equation, both [149Jump To The Next Citation Point] and [201Jump To The Next Citation Point] use the following expression

√ --- ∂( − gBi) ∂ √ --- i j √--- i j -----------+ ---j( − gV B − − gB V ) = 0, (66 ) ∂t ∂x
which differs from the expression used by [24Jump To The Next Citation Point] both from the use of i &tidle;v and the presence of √ -- γ and √ --- − g accompanying the time and spatial derivatives, respectively, in the latter case.

On the other hand, in the approach followed by Shibata and Sekiguchi [366Jump To The Next Citation Point], the evolution equations for the energy-momentum equations are also derived as in [24Jump To The Next Citation Point], 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 e6φ, φ being a spacetime conformal factor, motivated by the formulation used in the numerical code for solving Einstein’s equations, namely BSSN [363Jump To The Next Citation Point39Jump To The Next Citation Point]. (We note that such factors are also present in the formulation used by Shibata in the purely hydrodynamic case; see, e.g., [355Jump To The Next Citation Point].) Another important difference with [24Jump To The Next Citation Point] is the choice of the transport velocity i i 0 v = u ∕u throughout (instead of &tidle;vi) 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 ℬi = e6φ(W bi − αbtui), with 0 W = αu.

Correspondingly, in the conservative formulation of Anninos et al. [20Jump To The Next Citation Point] 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 [20Jump To The Next Citation Point] is different to those of [24Jump To The Next Citation Point] and [366Jump To The Next Citation Point], 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 [149Jump To The Next Citation Point] and [201Jump To The Next Citation Point], as the indices of the energy-momentum tensor are both chosen to be contravariant (instead of 0 T 0 and 0 T i). More precisely, the corresponding conserved variables are j √ --- 0j S = − gT and √ --- 00 ɛ = − gT. 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 4π 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

j ( i) ∂-ℬ- + -∂--(ℬjV i − ℬiV j) = η-∂- ∂-ℬ- , (67 ) ∂t ∂xi ∂xj ∂xi
with i √---√ -- i i0 ℬ = 4 π γW (b − V b ) and η the same coefficient as in Equation (56View Equation).

Next, the formulation of Duez et al. [108Jump To The Next Citation Point] 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 [20Jump To The Next Citation Point], 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 [149Jump To The Next Citation Point201Jump To The Next Citation Point], i.e., direct components of the energy-momentum tensor, namely √ --- − gT0i. However, for the remaining part of the energy-momentum conservation law, which leads to the energy equation, Duez et al. [108Jump To The Next Citation Point] do not follow [149Jump To The Next Citation Point201Jump To The Next Citation Point] and, instead of using √ −-gT 0 0, they choose τ, the same variable as [24Jump To The Next Citation Point], also subtracting the continuity equation, namely √ -- μν τ = γn μnνT − ρ ⋆. Finally, the induction equation is written as in Equation (67View Equation), but without the divergence cleanser coefficient (η = 0).

Finally, in the formulation of Del Zanna et al. [91Jump To The Next Citation Point] 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. [24Jump To The Next Citation Point] the same definition for the three-velocity is used, along with the same variable for the relativistic density D (densitized with a factor γ1∕2). For the Euler equation, components of the stress-energy tensor 0 Tj are used, as in the approaches of [149Jump To The Next Citation Point201Jump To The Next Citation Point108Jump To The Next Citation Point], 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 [91Jump To The Next Citation Point] take √ −-gT μνnν. The continuity equation is not subtracted from the energy equation. As far as the induction equation is concerned, the same conservative expression used by [149Jump To The Next Citation Point201Jump To The Next Citation Point] is employed.


  Go to previous page Go up Go to next page