Mon. Not. R. Astron. Soc. 000, 1–14 (2007) Printed 11 February 2013 (MN LATEX style file v2.2)
Equation of State in Relativistic Magnetohydrodynamics:
variable versus constant adiabatic index
A. Mignone1,2 ⋆ and Jonathan C. McKinney3⋆
1INAF Osservatorio Astronomico di Torino, 10025 Pino Torinese, Italy
2Dipartimento di Fisica Generale dell’Università, Via Pietro Giuria 1, I-10125 Torino, Italy
3Institute for Theory and Computation, Center for Astrophysics, Harvard University, 60 Garden St., Cambridge, MA, 02138
Accepted 2007 April 12. Received 2007 April 12; in original form 2007 January 25
ABSTRACT
The role of the equation of state for a perfectly conducting, relativistic magnetized
fluid is the main subject of this work. The ideal constant Γ-law equation of state,
commonly adopted in a wide range of astrophysical applications, is compared with a
more realistic equation of state that better approximates the single-specie relativistic
gas. The paper focus on three different topics. First, the influence of a more realis-
tic equation of state on the propagation of fast magneto-sonic shocks is investigated.
This calls into question the validity of the constant Γ-law equation of state in problems
where the temperature of the gas substantially changes across hydromagnetic waves.
Second, we present a new inversion scheme to recover primitive variables (such as
rest-mass density and pressure) from conservative ones that allows for a general equa-
tion of state and avoids catastrophic numerical cancellations in the non-relativistic
and ultrarelativistic limits. Finally, selected numerical tests of astrophysical relevance
(including magnetized accretion flows around Kerr black holes) are compared using
different equations of state. Our main conclusion is that the choice of a realistic equa-
tion of state can considerably bear upon the solution when transitions from cold to
hot gas (or viceversa) are present. Under these circumstances, a polytropic equation
of state can significantly endanger the solution.
Key words: equation of state - relativity - hydrodynamics shock waves - methods:
numerical - MHD
1 INTRODUCTION
Recent developments in numerical hydrodynamics have
made a breach in the understanding of astrophysical phe-
nomena commonly associated with relativistic magnetized
plasmas. Existence of such flows has nowadays been largely
witnessed by observations indicating superluminal motion in
radio loud active galactic nuclei and galactic binary systems,
as well as highly energetic events occurring in proximity of
X-ray binaries and super-massive black holes. Strong evi-
dence suggests that the two scenarios may be closely related
and that the production of relativistic collimated jets results
from magneto-centrifugal mechanisms taking place in the in-
ner regions of rapidly spinning accretion disks (Meier et al.
2001).
Due to the high degree of nonlinearity present in the
equations of relativistic magnetohydrodynamics (RMHD
henceforth), analytical models are often of limited appli-
cability, relying on simplified assumptions of time inde-
⋆ E-mail:mignone@to.astro.it(AM);jmckinney@cfa.harvard.edu(JCM)
pendence and/or spatial symmetries. For this reason, they
are frequently superseded by numerical models that appeal
to a consolidated theory based on finite difference meth-
ods and Godunov-type schemes. The propagation of rel-
ativistic supersonic jets without magnetic field has been
studied, for instance, in the pioneering work of van Putten
(1993); Duncan & Hughes (1994) and, subsequently, by
Mart́ı et al. (1997); Hardee et al. (1998); Aloy et al. (1999);
Mizuta et al. (2004) and references therein. Similar inves-
tigations in presence of poloidal and toroidal magnetic
fields have been carried on by Nishikawa et al. (1997);
Koide (1997); Komissarov (1999) and more recently by
Leismann et al. (2005); Mignone et al. (2005).
The majority of analytical and numerical models, in-
cluding the aforementioned studies, makes extensive use of
the polytropic equation of state (EoS henceforth), for which
the specific heat ratio is constant and equal to 5/3 (for a
cold gas) or to 4/3 (for a hot gas). However, the theory of
relativistic perfect gases (Synge 1957) teaches that, in the
limit of negligible free path, the ratio of specific heats can-
not be held constant if consistency with the kinetic theory
http://arxiv.org/abs/0704.1679v1
2 A. Mignone and J.C. McKinney
is to be required. This was shown in an even earlier work by
Taub (1948), where a fundamental inequality relating spe-
cific enthalpy and temperature was proved to hold.
Although these results have been known for many
decades, only few investigators seem to have faced this im-
portant aspect. Duncan et al. (1996) suggested, in the con-
text of extragalactic jets, the importance of self-consistently
computing a variable adiabatic index rather than using a
constant one. This may be advisable, for example, when
the dynamics is regulated by multiple interactions of shock
waves, leading to the formation of shock-heated regions in
an initially cold gas. Lately, Scheck et al. (2002) addressed
similar issues by investigating the long term evolution of
jets with an arbitrary mixture of electrons, protons and
electron-positron pairs. Similarly, Meliani et al. (2004) con-
sidered thermally accelerated outflows in proximity of com-
pact objects by adopting a variable effective polytropic index
to account for transitions from non-relativistic to relativis-
tic temperatures. Similar considerations pertain to models
of Gamma Ray Burst (GRB) engines including accretion
discs, which have an EoS that must account for a combi-
nation of protons, neutrons, electrons, positrons, and neu-
trinos, etc. and must include the effects of electron degen-
eracy, neutronization, photodisintegration, optical depth of
neutrinos, etc. (Popham et al. 1999; Di Matteo et al. 2002;
Kohri & Mineshige 2002; Kohri et al. 2005). However, for
the disk that is mostly photodisintegrated and optically
thin to neutrinos, a decent approximation of such EoS is
a variable Γ-law with Γ = 5/3 when the temperature is be-
low mec
2/kb and Γ = 4/3 when above mec
2/kb due to the
production of positrons at high temperatures that gives a
relativistic plasma (Broderick, McKinney, Kohri in prep.).
Thus, the variable EoS considered here may be a reasonable
approximation of GRB disks once photodisintegration has
generated mostly free nuclei.
The additional complexity introduced by more elabo-
rate EoS comes at the price of extra computational cost since
the EoS is frequently used in the process of obtaining numer-
ical solutions, see for example, Falle & Komissarov (1996).
Indeed, for the Synge gas, the correct EoS does not have a
simple analytical expression and the thermodynamics of the
fluid becomes entirely formulated in terms of the modified
Bessel functions.
Recently Mignone et al. (2005a, MPB henceforth) in-
troduced, in the context of relativistic non-magnetized flows,
an approximate EoS that differs only by a few percent from
the theoretical one. The advantage of this approximate EoS,
earlier adopted by Mathews (1971), is its simple analytical
representation. A slightly better approximation, based on an
analytical expression, was presented by Ryu et al. (2006).
In the present work we wish to discuss the role of the
EoS in RMHD, with a particular emphasis to the one pro-
posed by MPB, properly generalized to the context of rel-
ativistic magnetized flows. Of course, it is still a matter of
debate the extent to which equilibrium thermodynamic prin-
ciples can be correctly prescribed when significant deviations
from the single-fluid ideal approximation may hold (e.g.,
non-thermal particle distributions, gas composition, cosmic
ray acceleration and losses, anisotropy, and so forth). Nev-
ertheless, as the next step in a logical course of action, we
will restrict our attention to a single aspect - namely the use
of a constant polytropic versus a variable one - and we will
ignore the influence of such non-ideal effects (albeit poten-
tially important) on the EoS.
In §2, we present the relevant equations and discuss the
properties of the new EoS versus the more restrictive con-
stant Γ-law EoS. In §3, we consider the propagation of fast
magneto-sonic shock waves and solve the jump conditions
across the front using different EoS. As we shall see, this
calls into question the validity of the constant Γ-law EoS
in problems where the temperature of the gas substantially
changes across hydromagnetic waves. In §4, we present nu-
merical simulations of astrophysical relevance such as blast
waves, axisymmetric jets, and magnetized accretion disks
around Kerr black holes. A short survey of some existing
models is conducted using different EoS’s in order to deter-
mine if significant interesting deviations arise. These results
should be treated as a guide to some possible avenues of
research rather than as the definitive result on any individ-
ual topic. Results are summarized in §5. In the Appendix,
we present a description of the primitive variable inversion
scheme.
2 RELATIVISTIC MHD EQUATIONS
In this section we present the equations of motion for rel-
ativistic MHD, discuss the validity of the ideal gas EoS as
applied to a perfect gas, and review an alternative EoS that
properly models perfect gases in both the hot (relativistic)
and cold (non-relativistic) regimes.
2.1 Equations of Motion
Our starting point are the relativistic MHD equations in
conservative form:
+ ∇ ·
vv − bb + Ipt
vB − Bv
= 0 , (1)
together with the divergence-free constraint ∇·B = 0, where
v is the velocity, γ is the Lorentz factor, wt ≡ (ρh+p+b2) is
the relativistic total (gas+magnetic) enthalpy, pt = p+ b
is the total (gas+magnetic) fluid pressure, B is the lab-frame
field, and the field in the fluid frame is given by
= γ{v · B, B
(v · B)}, (2)
with an energy density of
|b|2 = |B|
+ (v · B)2. (3)
Units are chosen such that the speed of light is equal to one.
Notice that the fluxes entering in the induction equation
are the components of the electric field that, in the infinite
conductivity approximation, become
Ω = −v × B . (4)
The non-magnetic case is recovered by letting B → 0 in the
previous expressions.
Equation of state in RMHD 3
Figure 1. Equivalent Γ (top left), specific enthalpy (top right),
sound speed (bottom left) and specific internal energy (bottom
right) as functions of temperature Θ = p/ρ. Different lines corre-
spond to the various EoS mentioned the text: the ideal Γ = 5/3-
law (dotted line), ideal Γ = 4/3-law (dashed line), TM EoS (solid
line). For clarity the Synge-gas (dashed-dotted line) has been
plotted only in the top left panel, where the “unphysical region”
marks the area where Taub’s inequality is not fulfilled.
The conservative variables are, respectively, the labora-
tory density D, the three components of momentum mk and
magnetic field Bk and the total energy density E:
D = ργ , (5)
mk = (Dhγ + |B|2)vk − (v · B)Bk , (6)
E = Dhγ − p + |B|
|v|2|B|2 − (v · B)2
, (7)
The specific enthalpy h and internal energy ǫ of the gas
are related by
h = 1 + ǫ +
, (8)
and an additional equation of state relating two thermody-
namical variables (e.g. ρ and ǫ) must be specified for proper
closure. This is the subject of the next section.
Equations (5)–(7) are routinely used in numerical codes
to recover conservative variables from primitive ones (e.g.,
ρ, v, p and B). The inverse relations cannot be cast in
closed form and require the solution of one or more non-
linear equations. Noble et al. (2006) review several methods
of inversion for the constant Γ-law, for which ρǫ = p/(Γ−1).
We present, in Appendix A, the details of a new inversion
procedure suitable for a more general EoS.
2.2 Equation of State
Proper closure to the conservation law (1) is required in
order to solve the equations. This is achieved by specifying
an EoS relating thermodynamic quantities. The theory of
relativistic perfect gases shows that the specific enthalpy is
a function of the temperature Θ = p/ρ alone and it takes
the form (Synge 1957)
K3(1/Θ)
K2(1/Θ)
, (9)
where K2 and K3 are, respectively, the order 2 and 3 modi-
fied Bessel functions of the second kind. Equation (9) holds
for a gas composed of material particles with the same mass
and in the limit of small free path when compared to the
sound wavelength.
Direct use of Eq. (9) in numerical codes, however,
results in time-consuming algorithms and alternative ap-
proaches are usually sought. The most widely used and pop-
ular one relies on the choice of the constant Γ-law EoS
h = 1 +
Γ − 1
Θ , (10)
where Γ is the constant specific heat ratio. However, Taub
(1948) showed that consistency with the relativistic kinetic
theory requires the specific enthalpy h to satisfy
(h − Θ) (h − 4Θ) > 1 , (11)
known as Taub’s fundamental inequality. Clearly the con-
stant Γ-law EoS does not fulfill (11) for an arbitrary choice
of Γ, while (9) certainly does. This is better understood in
terms of an equivalent Γeq, conveniently defined as
Γeq =
h − 1
h − 1 − Θ
, (12)
and plotted in the top left panel of Fig. 1 for different EoS.
In the limit of low and high temperatures, the physically
admissible region is delimited, respectively, by Γeq 6 5/3
(for Θ → 0) and Γeq 6 4/3 (for Θ → ∞). Indeed, Taub’s
inequality is always fulfilled when Γ 6 4/3 while it cannot
be satisfied for Γ > 5/3 for any positive value of the tem-
perature.
In a recent paper, Mignone et al. (2005a) showed that
if the equal sign is taken in Eq. (11), an equation with the
correct limiting values may be derived. The resulting EoS
(TM henceforth), previously introduced by Mathews (1971),
can be solved for the enthalpy, yielding
Θ2 + 1 , (13)
or, using ρh = ρ + ρǫ + p in (11) with the equal sign,
ρǫ (ρǫ + 2ρ)
3 (ρǫ + ρ)
ǫ + 2
ǫ + 1
. (14)
Direct evaluation of Γeq using (13) shows that the TM EoS
differs by less than 4% from the theoretical value given by
the relativistic perfect gas EoS (9). The proposed EoS be-
haves closely to the Γ = 4/3 law in the limit of high tem-
peratures, whereas reduces to the Γ = 5/3 law in the cold
gas limit. For intermediate temperatures, thermodynamical
quantities (such as specific internal energy, enthalpy and
sound speed) smoothly vary between the two limiting cases,
as illustrated in Fig. 1. In this respect, Eq. (13) greatly im-
proves over the constant Γ-law EoS and, at the same time,
offers ease of implementation over Eq. (9). Since thermody-
namics is frequently invoked during the numerical solution
of (1), it is expected that direct implementation of Eq. (13)
4 A. Mignone and J.C. McKinney
in numerical codes will result in faster and more efficient
algorithms.
Thermodynamical quantities such as sound speed and
entropy are computed from the 2nd law of thermodynamics,
− d log p , (15)
where S is the entropy. From the definition of the sound
speed,
, (16)
and using de = hdρ (at constant S), one finds the useful
expression
ḣ − 1
Γ-law EoS ,
5h − 8Θ
h − Θ TM EoS .
where we set ḣ = dh/dΘ. In a similar way, direct integration
of (15) yields S = k log σ with
Γ-law EoS ,
(h − Θ) TM EoS .
with h given by (13).
3 PROPAGATION OF FAST
MAGNETO-SONIC SHOCKS
Motivated by the previous results, we now investigate the
role of the EoS on the propagation of magneto-sonic shock
waves. To this end, we proceed by constructing a one-
parameter family of shock waves with different velocities,
traveling in the positive x direction. States ahead and be-
hind the front are labeled with U 0 and U 1, respectively, and
are related by the jump conditions
vs [U ] = [F (U )] , (19)
where vs is the shock speed and [q] = q1 − q0 is the jump
across the wave for any quantity q. The set of jump condi-
tions (19) may be reduced (Lichnerowicz 1976) to the fol-
lowing five positive-definite scalar invariants
[J ] = 0 , (20)
[hη] = 0 , (21)
[H] =
= 0 , (22)
p + b2/2
[h/ρ]
= 0 , (23)
+ 2H [p] + 2
= 0 , (24)
where
J = ργγs(vs − vx) , (25)
is the mass flux across the shock, and
η = −J
(v · B) + γs
Bx . (26)
Figure 2. Compression ratio (top panels), internal energy (mid-
dle panels) and downstream Mach number (bottom panels) as
functions of the shock four-velocity γsvs. The profiles give the so-
lution to the shock equation for the non magnetic case. Plots on
the left have zero tangential velocity ahead of the front, whereas
plots on right are initialized with vy0 = 0.99. Axis spacing is
logarithmic. Solid, dashed and dotted lines correspond to the so-
lutions obtained with the TM EoS and the Γ = 4/3 and Γ = 5/3
laws, respectively.
Here γs denotes the Lorentz factor of the shock. Fast or
slow magneto-sonic shocks may be discriminated through
the condition α0 > α1 > 0 (for the formers) or α1 < α0 < 0
(for the latters), where α = h/ρ −H.
We consider a pre-shock state characterized by a cold
(p0 = 10
−4) gas with density ρ = 1. Without loss of gen-
erality, we choose a frame of reference where the pre-shock
velocity normal to the front vanishes, i.e., vx0 = 0. Notice
that, for a given shock speed, J2 can be computed from the
pre-shock state and thus one has to solve only Eqns. (21)–
(24).
3.1 Purely Hydrodynamical Shocks
In the limit of vanishing magnetic field, only Eqns. (23) and
(24) need to be solved. Since J2 is given, the problem sim-
Equation of state in RMHD 5
plifies to the 2 × 2 nonlinear system of equations
[h/ρ]
= 0 , (27)
= 0 . (28)
We solve the previous equations starting from vs = 0.2,
for which we were able to provide a sufficiently close guess to
the downstream state. Once the p1 and ρ1 have been found,
we repeat the process by slowly increasing the shock velocity
vs and using the previously converged solution as the initial
guess for the new value of vs.
Fig. 2 shows the compression ratio, post-shock inter-
nal energy ǫ1 and Mach number v1/cs1 as functions of the
shock four velocity vsγs. For weakly relativistic shock speeds
and vanishing tangential velocities (left panels), density and
pressure jumps approach the classical (i.e. non relativistic)
strong shock limit at γsvs ≈ 0.1, with the density ratio be-
ing 4 or 7 depending on the value of Γ (5/3 or 4/3, respec-
tively). The post-shock temperature keeps non-relativistic
values (Θ ≪ 1) and the TM EoS behaves closely to the
Γ = 5/3 case, as expected.
With increasing shock velocity, the compression ratio
does not saturate to a limiting value (as in the classical
case) but keeps growing at approximately the same rate
for the constant Γ-law EoS cases, and more rapidly for the
TM EoS. This can be better understood by solving the
jump conditions in a frame of reference moving with the
shocked material and then transforming back to our origi-
nal system. Since thermodynamics quantities are invariant
one finds that, in the limit h1 ≫ h0 ≈ 1, the internal energy
becomes ǫ1 = γ1 − 1 and the compression ratio takes the
asymptotic value
= γ1 +
γ1 + 1
Γ − 1
, (29)
when the ideal EoS is adopted. Since γ1 can take arbitrarily
large values, the downstream density keeps growing indefi-
nitely. At the same time, internal energy behind the shock
rises faster than the rest-mass energy, eventually leading to
a thermodynamically relativistic configuration. In absence
of tangential velocities (left panels in Fig. 2), this transi-
tion starts at moderately high shock velocities (γsvs & 1)
and culminates when the shocked gas heats up to relativis-
tic temperatures (Θ ∼ 1 ÷ 10) for γsvs & 10. In this regime
the TM EoS departs from the Γ = 5/3 case and merges on
the Γ = 4/3 curve. For very large shock speeds, the Mach
number tends to the asymptotic value (Γ−1)−1/2, regardless
of the frame of reference.
Inclusion of tangential velocities (right panels in Fig. 2)
leads to an increased mass flux (J2 ∝ γ20) and, consequently,
to higher post-shock pressure and density values. Still, since
pressure grows faster than density, temperature in the post-
shock flow strains to relativistic values even for slower shock
velocities and the TM EoS tends to the Γ = 4/3 case at even
smaller shock velocities (γsvs & 2).
Generally speaking, at a given shock velocity, density
and pressure in the shocked gas attain higher values for lower
Γeq. Downstream temperature, on the other hand, follows
the opposite trend being higher as Γeq → 5/3 and lower
when Γeq → 4/3.
Figure 3. Compression ratio (top), downstream plasma β (mid-
dle) and magnetic field strength (bottom) as function of the shock
four-velocity γsvs with vanishing tangential component of the
velocity. The magnetic field makes an angle π/6 (left) and π/2
(right) with the shock normal. The meaning of the different lines
is the same as in Fig. 2.
3.2 Magnetized Shocks
In presence of magnetic fields, we solve the 3 × 3 nonlin-
ear system given by Eqns. (22), (23) and (24), and di-
rectly replace η1 = η0h0/h1 with the aid of Eq. (21).
The magnetic field introduces three additional parameters,
namely, the thermal to magnetic pressure ratio (β ≡ 2p/b2)
and the orientation of the magnetic field with respect to
the shock front and to the tangential velocity. This is ex-
pressed by the angles αx and αy such that Bx = |B| cos αx,
By = |B| sin αx cos αy , Bz = |B| sin αx sin αy . We restrict
our attention to the case of a strongly magnetized pre-shock
flow with β0 ≡ 2p0/b20 = 10−2.
Fig. 3 shows the density, plasma β and magnetic pres-
sure ratios versus shock velocity for αx = π/6 (left panels)
and αx = π/2 (perpendicular shock, right panels). Since
there is no tangential velocity, the solution depends on one
angle only (αx) and the choice of αy is irrelevant. For small
shock velocities (γsvs . 0.4), the front is magnetically driven
with density and pressure jumps attaining lower values than
the non-magnetized counterpart. A similar behavior is found
in classical MHD (Jeffrey & Taniuti 1964). Density and
6 A. Mignone and J.C. McKinney
Figure 4. Density ratio (top), downstream plasma β (middle)
and magnetic field strength (bottom) as function of γsvs when
the tangential component of the upstream velocity is vt = 0.99.
The magnetic field and the shock normal form an angle π/6. The
tangential components of magnetic field and velocity are aligned
(left) and orthogonal (right).Different lines have the same mean-
ing as in Fig. 2.
magnetic compression ratios across the shock reach the clas-
sical values around γsvs ≈ 1 (rather than γsvs ≈ 0.1 as in
the non-magnetic case) and increase afterwards. The mag-
netic pressure ratio grows faster for the perpendicular shock,
whereas internal energy and density show little dependence
on the orientation angle αx. As expected, the TM EoS mim-
ics the constant Γ = 5/3 case at small shock velocities. At
γsvs . 0.46, the plasma β exceeds unity and the shock starts
to be pressure-dominated. In other words, thermal pressure
eventually overwhelms the Lorentz force and the shock be-
comes pressure-driven for velocities of the order of vs ≈ 0.42.
When γsvs & 1, the internal energy begins to become com-
parable to the rest mass energy (c2) and the behavior of the
TM EoS detaches from the Γ = 5/3 curve and slowly joins
the Γ = 4/3 case. The full transition happens in the limit of
strongly relativistic shock speeds, γsvs . 10.
Inclusion of transverse velocities in the right state af-
fects the solution in a way similar to the non-magnetic case.
Relativistic effects play a role already at small velocities
because of the increased inertia of the pre-shock state in-
troduced by the upstream Lorentz factor. For αx = π/6
Figure 5. Density contrast (top), plasma β (middle) and mag-
netic field strength (bottom) for vt = 0.99. The magnetic field is
purely transverse and aligned with the tangential component of
velocity on the left, while it is orthogonal on the right. Different
lines have the same meaning as in Fig. 2.
(Fig. 4), the compression ratio does not drop to small val-
ues and keeps growing becoming even larger (. 400) than
the previous case when vt = 0. The same behavior is re-
flected on the growth of magnetic pressure that, in addi-
tion, shows more dependence on the relative orientation of
the velocity and magnetic field projections in the plane of
the front. When αy = π/2, indeed, magnetic pressure at-
tains very large values (b2/b20 . 10
4, bottom right panel in
Fig. 4). Consequently, this is reflected in a decreased post-
shock plasma β. For the TM EoS, the post-shock properties
of the flow begin to resemble the Γ = 4/3 behavior at lower
shock velocities than before, γsvs ≈ 2 ÷ 3. Similar consid-
erations may be done for the case of a perpendicular shock
(αx = π/2, see Fig. 5), although the plasma β saturates
to larger values thus indicating larger post-shock pressures.
Again, the maximum increase in magnetic pressure occurs
when the velocity and magnetic field are perpendicular.
4 NUMERICAL SIMULATIONS
With the exception of very simple flow configurations, the
solution of the RMHD fluid equations must be carried out
Equation of state in RMHD 7
Figure 6. Solution to the mildly relativistic blast wave (problem
1) at t = 0.4. From left to right, the different profiles give den-
sity, thermal pressure, total pressure (top panels), the three com-
ponents of velocity (middle panel) and magnetic fields (bottom
panels). Computations with the TM EoS and constant Γ = 5/3
EoS are shown using solid and dotted lines, respectively.
numerically. This allows an investigation of highly nonlinear
regimes and complex interactions between multiple waves.
We present some examples of astrophysical relevance, such
as the propagation of one dimensional blast waves, the prop-
agation of axisymmetric jets, and the evolution of magne-
tized accretion disks around Kerr black holes. Our goal is to
outline the qualitative effects of varying the EoS for some in-
teresting astrophysical problems rather than giving detailed
results on any individual topic.
Direct numerical integration of Eq. (1) has been
achieved using the PLUTO code (Mignone et al. 2007) in
§4.1, §4.2 and HARM (Gammie et al. 2003) in §4.3. The new
primitive variable inversion scheme presented in Appendix
A has been implemented in both codes and the results pre-
sented in §4.1 were used for code validation. The novel in-
version scheme offers the advantage of being suitable for a
more general EoS and avoiding catastrophic cancellation in
the non-relativistic and ultrarelativistic limits.
4.1 Relativistic Blast Waves
A shock tube consists of a sharp discontinuity separat-
ing two constant states. In what follows we will be con-
sidering the one dimensional interval [0, 1] with a discon-
tinuity placed at x = 0.5. For the first test problem,
states to the left and to the right of the discontinuity are
given by (ρ, p,By , Bz)L = (1, 30, 6, 6) for the left state and
(ρ, p, By , Bz)R = (1, 1, 0.7, 0.7) for the right state. This re-
sults in a mildly relativistic configuration yielding a max-
imum Lorentz factor of 1.3 6 γ 6 1.4. The second test
Figure 7. Solution to the strong relativistic blast wave (problem
2) at t = 0.4. From left to right, the different profiles give den-
sity, thermal pressure, total pressure (top panels), the three com-
ponents of velocity (middle panel) and magnetic fields (bottom
panels). Computations with the TM EoS and constant Γ = 5/3
EoS are shown using solid and dotted lines, respectively.
consists of a left state given by (ρ, p,By , Bz)L = (1, 10
3, 7, 7)
and a right state (ρ, p, By, Bz)R = (1, 0.1, 0.7, 0.7). This con-
figuration involves the propagation of a stronger blast wave
yielding a more relativistic configuration (3 6 γ 6 3.5). For
both states, we use a base grid with 800 zones and 6 levels
of refinement (equiv. resolution = 800 · 26) and evolve the
solution up to t = 0.4.
Computations carried with the ideal EoS with Γ = 5/3
and the TM EoS are shown in Fig. 6 and Fig. 7 for the
first and second shock tube, respectively. From left to right,
the wave pattern is comprised of a fast and slow rarefac-
tions, a contact discontinuity and a slow and a fast shocks.
No rotational discontinuity is observed. Compared to the
Γ = 5/3 case, one can see that the results obtained with
the TM EoS show considerable differences. Indeed, waves
propagate at rather smaller velocities and this is evident at
the head and the tail points of the left-going magneto-sonic
rarefaction waves. From a simple analogy with the hydrody-
namic counterpart, in fact, we know that these points prop-
agate increasingly faster with higher sound speed. Since the
sound speed ratio of the TM and Γ = 5/3 is always less
than one (see, for instance, the bottom left panel in Fig. 1),
one may reasonably predict slower propagation speed for the
Riemann fans when the TM EoS is used. Furthermore, this is
confirmed by computations carried with Γ = 4/3 that shows
even slower velocities. Similar conclusions can be drawn for
the shock velocities. The reason is that the opening of the
Riemann fan of the TM equation state is smaller than the
Γ = 5/3 case, because the latter always over-estimates the
sound speed. The higher density peak behind the slow shock
8 A. Mignone and J.C. McKinney
Figure 8. Jet velocity as a function of the Mach number for
different values of the initial density contrast η. The beam Lorentz
factor is the same for all plots, γb = 10. Solid, dashed and dotted
lines correspond to the solutions obtained with the TM EoS and
the Γ = 4/3 and Γ = 5/3 laws, respectively.
follows from the previous considerations and the conserva-
tion of mass across the front.
4.2 Propagation of Relativistic Jets
Relativistic, pressure-matched jets are usually set up by
injecting a supersonic cylindrical beam with radius rb
into a uniform static ambient medium (see, for instance,
Mart́ı et al. 1997). The dynamical and morphological prop-
erties of the jet and its interaction with the surrounding are
most commonly investigated by adopting a three parameter
set: the beam Lorentz factor γb, Mach number Mb = vb/cs
and the beam to ambient density ratio η = ρb/ρm. The
presence of a constant poloidal magnetic field introduces a
fourth parameter βb = 2pb/b
2, which specifies the thermal
to magnetic pressure ratio.
4.2.1 One Dimensional Models
The propagation of the jet itself takes place at the velocity
Vj , defined as the speed of the working surface that sepa-
rates shocked ambient fluid from the beam material. A one-
dimensional estimate of Vj (for vanishing magnetic fields)
can be derived from momentum flux balance in the frame of
the working surface (Mart́ı et al. 1997). This yields
ηhb/hm
1 + γb
ηhb/hm
, (30)
where hb and hm are the specific enthalpies of the beam and
the ambient medium, respectively. For given γb and den-
sity contrast η, Eq. (30) may be regarded as a function of
Figure 9. Computed results for the non magnetized jet at t = 90
for the ideal EoS (Γ = 5/3 and Γ = 4/3, top and middle panels)
and the TM EoS (bottom panel), respectively. The lower and
upper half of each panels shows the gray-scale map of density
and internal energy in logarithmic scale.
Figure 10. Position of the working surface as a function of time
for Γ = 5/3 (circles), Γ = 4/3 (stars) and the TM EoS (dia-
monds). Solid, dotted and dashed lines gives the one-dimensional
expectation.
Equation of state in RMHD 9
Figure 11. Density and magnetic field for the magnetized jet at
t = 80 (first and second panels from top) and at t = 126 (third
and fourth panels). Computations were carried with 40 zones per
beam radius with the TM EoS.
the Mach number alone that uniquely specifies the pres-
sure pb through the definitions of the sound speed, Eq. (17).
For the constant Γ-law EoS the inversion is straightforward,
whereas for the TM EoS one finds, using the substitution
Θ = 2/3 sinh x,
pb = η
1 − t2m
, (31)
where tm satisfies the negative branch of the quadratic equa-
15 − 6M
24 − 10M
+ 9 = 0 , (32)
with t = tanh x. In Fig. 8 we show the jet velocity for in-
creasing Mach numbers (or equivalently, decreasing sound
speeds) and different density ratios η = 10−5, 10−3, 10−1, 10.
The Lorentz beam factor is γb = 10. Prominent discrepan-
cies between the selected EoS arise at low Mach numbers,
where the relative variations of the jet speed between the
constant Γ and the TM EoS’s can be more than 50%. This
regime corresponds to the case of a hot jet (Θ ≈ 10 in the
η = 10−3 case) propagating into a cold (Θ ≈ 10−3) medium,
for which neither the Γ = 4/3 nor the Γ = 5/3 approxima-
tion can properly characterize both fluids.
4.2.2 Two Dimensional Models
Of course, Eq. (30) is strictly valid for one-dimensional flows
and the question remains as to whether similar conclusions
can be drawn in more than one dimension. To this end we
investigate, through numerical simulations, the propagation
of relativistic jets in cylindrical axisymmetric coordinates
(r, z). We consider two models corresponding to different
sets of parameters and adopt the same computational do-
main [0, 12] × [0, 50] (in units of jet radius) with the beam
being injected at the inlet region (r 6 1, z = 0). Jets are in
pressure equilibrium with the environment.
In the first model, the density ratio, beam Lorentz fac-
tor and Mach number are given, respectively, by η = 10−3,
γb = 10 and Mb = 1.77. Magnetic fields are absent. Inte-
gration are carried at the resolution of 20 zones per beam
radius using the relativistic Godunov scheme described in
MPB. Computed results showing density and internal en-
ergy maps at t = 90 are given in Fig. 9 for Γ = 5/3, Γ = 4/3
and the TM EoS. The three different cases differ in several
morphological aspects, the most prominent one being the
position of the leading bow shock, z ≈ 18 when Γ = 5/3,
z ≈ 48 for Γ = 4/3 and z ≈ 33 for the TM EoS. Smaller
values of Γ lead to larger beam internal energies and there-
fore to an increased momentum flux, in agreement with the
one dimensional estimate (30). This favors higher propaga-
tion velocities and it is better quantified in Fig. 10 where
the position of the working surface is plotted as a function
of time and compared with the one dimensional estimate.
For the cold jet (Γ = 5/3), the Mach shock exhibits a larger
cross section and is located farther behind the bow shock
when compared to the other two models. As a result, the
jet velocity further decreases promoting the formation of a
thicker cocoon. On the contrary, the hot jet (Γ = 4/3) prop-
agates at the highest velocity and the cocoon has a more
elongated shape. The beam propagates almost undisturbed
and cross-shocks are weak. Close to is termination point,
the beam widens and the jet slows down with hot shocked
gas being pushed into the surrounding cocoon at a higher
rate. Integration with the TM EoS reveals morphological
and dynamical properties more similar to the Γ = 4/3 case,
although the jet is ≈ 40% slower. At t = 90 the beam does
not seem to decelerate and its speed remains closer to the
one-dimensional expectation. The cocoon develops a thin-
ner structure with a more elongated conical shape and cross
shocks form in the beam closer to the Mach disk.
In the second case, we compare models C2-pol-1 and
B1-pol-1 of Leismann et al. (2005) (corresponding to an
ideal gas with Γ = 5/3 and Γ = 4/3, respectively) with
the TM EoS adopting the same numerical scheme. For this
model, η = 10−2, vb = 0.99, Mb = 6 and the ambient
medium is threaded by a constant vertical magnetic field,
2pb. Fig. 11 shows the results at t = 80 and
t = 126, corresponding to the final integration times shown
in Leismann et al. (2005) for the selected values of Γ. For the
sake of conciseness, integration pertaining to the TM EoS
only are shown and the reader is reminded to the original
work by Leismann et al. (2005) for a comprehensive descrip-
tion. Compared to ideal EoS cases, the jet shown here pos-
sesses morphological and dynamical properties intermediate
between the hot (Γ = 4/3) and the cold (Γ = 5/3) cases. As
expected, the jet propagates slower than in model B1-pol-1
(hot jet), but faster than the cold one (C2-pol-1). The head
of the jet tends to form a hammer-like structure (although
less prominent than the cold case) towards the end of the
integration, i.e., for t & 100, but the cone remains more con-
fined at previous times. Consistently with model C2-pol-1,
the beam develops a series of weak cross shocks and outgo-
ing waves triggered by the interaction of the flow with bent
magnetic field lines. Although the magnetic field inhibits the
formation of eddies, turbulent behavior is still observed in
cocoon, where interior cavities with low magnetic fields are
10 A. Mignone and J.C. McKinney
Figure 12. Magnetized accretion flow around a Kerr black hole
for the ideal Γ-law EoS with Γ = 4/3. Shows the logarithm of the
rest-mass density in colour from high (red) to low (blue) values.
The magnetic field has been overlayed. This model demonstrates
more vigorous turbulence and a thicker corona that leads to a
more confined magnetized jet near the poles.
Figure 13. As in figure 12 but for Γ = 5/3. Compared to the
Γ = 4/3 model, there is less vigorous turbulence and the corona
is more sharply defined.
formed. In this respect, the jet seems to share more features
with the cold case.
Figure 14. As in figure 12 but for the TM EoS. This EoS leads
to turbulence that is less vigorous than in the Γ = 4/3 model but
more vigorous than in the Γ = 5/3 model. Qualitatively the TM
EoS leads to an accretion disk that behaves somewhere between
the behavior of the Γ = 4/3 and Γ = 5/3 models.
4.3 Magnetized Accretion near Kerr Black Holes
In this section we study time-dependent GRMHD numerical
models of black hole accretion in order to determine the ef-
fect of the EoS on the behavior of the accretion disk, corona,
and jet. We study three models similar to the models stud-
ied by McKinney & Gammie (2004) for a Kerr black hole
with a/M ≈ 0.94 and a disk with a scale height (H) to ra-
dius (R) ratio of H/R ∼ 0.3. The constant Γ-law EoS with
Γ = {4/3, 5/3} and the TM EoS are used. The initial torus
solution is in hydrostatic equilibrium for the Γ-law EoS, but
we use the Γ = 5/3 EoS as an initial condition for the TM
EoS. Using the Γ = 4/3 EoS as an initial condition for the
TM EoS did not affect the final quasi-stationary behavior
of the flow. The simplest question to ask is which value of Γ
will result in a solution most similar to the TM EoS model’s
solution.
More advanced questions involve how the structure
of the accretion flow depends on the EoS. The previ-
ous results of this paper indicate that the corona above
the disk seen in the simulations (De Villiers et al. 2003;
McKinney & Gammie 2004) will be most sensitive to the
EoS since this region can involve both non-relativistic and
relativistic temperatures. The corona is directly involved
is the production of a turbulent, magnetized, thermal disk
wind (McKinney & Narayan 2006a,b), so the disk wind is
also expected to depend on the EoS. The disk inflow near
the black hole has a magnetic pressure comparable to the gas
pressure (McKinney & Gammie 2004), so the EoS may play
a role here and affect the flux of mass, energy, and angular
momentum into the black hole. The magnetized jet asso-
ciated with the Blandford & Znajek solution seen in simu-
lations (McKinney & Gammie 2004; McKinney 2006) is not
expected to depend directly on the EoS, but may depend in-
Equation of state in RMHD 11
directly through the confining action of the corona. Finally,
the type of field geometries observed in simulations that
thread the disk and corona (Hirose et al. 2004; McKinney
2005) might depend on the EoS through the effect of the
stiffness (larger Γ leads to harder EoSs) of the EoS on the
turbulent diffusion of magnetic fields.
Figs. 12, 13 and 14 show a snapshot of the accretion
disk, corona, and jet at t ∼ 1000GM/c3 . Overall the re-
sults are quite comparable, as could be predicted since the
Γ = {4/3, 5/3} models studied in McKinney & Gammie
(2004) were quite similar. For all models, the field geometries
allowed are similar to that found in McKinney (2005). The
accretion rate of mass, specific energy, and specific angular
momentum are similar for all models, so the EoS appears to
have only a small effect on the flow through the disk near
the black hole.
The most pronounced effect is that the soft EoS (Γ =
4/3) model develops more vigorous turbulence due to the
non-linear behavior of the magneto-rotational instability
(MRI) than either the Γ = 5/3 or TM EoSs. This causes
the coronae in the Γ = 4/3 model to be slightly thicker and
to slightly more strongly confine the magnetized jet resulting
in a slight decrease in the opening angle of the magnetized
jet at large radii. Also, the Γ = 4/3 model develops a fast
magnetized jet at slightly smaller radii than the other mod-
els. An important consequence is that the jet opening angle
at large radii might depend sensitively on the EoS of the ma-
terial in the accretion disc corona. This should be studied in
future work.
5 CONCLUSIONS
The role of the EoS in relativistic magnetohydrodynamics
has been investigated both analytically and numerically. The
equation of state previously introduced by Mignone et al.
(2005a) (for non magnetized flows) has been extended to
the case where magnetic fields are present. The proposed
equation of state closely approximates the single-specie per-
fect relativistic gas, but it offers a much simpler analyti-
cal representation. In the limit of very large or very small
temperatures, for instance, the equivalent specific heat ratio
reduces, respectively, to the 4/3 or 5/3 limits.
The propagation of fast magneto-sonic shock waves has
been investigated by comparing the constant Γ laws to the
new equation of state. Although for small shock veloci-
ties the shock dynamics is well described by the cold gas
limit, dynamical and thermodynamical quantities (such as
the compression ratio, internal energy, magnetization and so
forth) substantially change across the wave front at moder-
ately or highly relativistic speeds. Eventually, for increasing
shock velocities, flow quantities in the downstream region
smoothly vary from the cold (Γ = 5/3) to the hot (Γ = 4/3)
regimes.
We numerically studied the effect of the EoS on shocks,
blast waves, the propagation of relativistic jets, and magne-
tized accretion flows around Kerr black holes. Our results
should serve as a useful guide for future more specific stud-
ies of each topic. For these numerical studies, we formu-
lated the inversion from conservative quantities to primitive
quantities that allows a general EoS and avoids catastrophic
numerical cancellation in the non-relativistic and ultrarela-
tivistic limits. The analytical and numerical models confirm
the general result that large temperature gradients cannot
be properly described by a polytropic EoS with constant
specific heat ratio. Indeed, when compared to a more re-
alistic EoS, for which the polytropic index is a function of
the temperature, considerable dynamical differences arises.
This has been repeatedly shown in presence of strong dis-
continuities, such shocks, across which the internal energy
can change by several order of magnitude.
We also showed that the turbulent behavior of magne-
tized accretion flows around Kerr black holes depends on the
EoS. The Γ = 4/3 EoS leads to more vigorous turbulence
than the Γ = 5/3 or TM EoSs. This affects the thickness of
the corona that confines the magnetized jet. Any study of
turbulence within the accretion disk, the subsequent genera-
tion of heat in the coronae, and the opening and acceleration
of the jet (especially at large radii where the cumulative dif-
ferences due to the EoS in the disc are largest) should use
an accurate EoS. The effect of the EoS on the jet opening
angle and Lorentz factor at large radii is a topic of future
study.
The proposed equation state holds in the limit where
effects due to radiation pressure, electron degeneracies and
neutrino physics can be neglected. It also omits potentially
crucial physical aspects related to kinetic processes (such
as suprathermal particle distributions, cosmic rays), plasma
composition, turbulence effects at the sub-grid levels, etc.
These are very likely to alter the equation of state by ef-
fectively changing the adiabatic index computed on merely
thermodynamic arguments. Future efforts should properly
address additional physical issues and consider more gen-
eral equations of state.
ACKNOWLEDGMENTS
We are grateful to our referee, P. Hughes, for his worthy
considerations and comments that led to the final form of
this paper. JCM was supported by a Harvard CfA Institute
for Theory and Computation fellowship. AM would like to
thank S. Massaglia and G. Bodo for useful discussions on
the jet propagation and morphology.
REFERENCES
Aloy, M. A., Ibáñez, J. M. , Mart́ı, J. M. , Gómez, J.-L.,
Müller, E. 1999, ApJL, 523, L125
Aloy, M. A., Ibáñez, J. M., Mart́ı, J. M., Müller, E. 1999,
ApJS, 122, 151
Anile, M., & Pennisi, S. 1987, Ann. Inst. Henri Poincaré,
46, 127
Anile, A. M. 1989, Relativistic Fluids and Magneto-fluids
(Cambridge: Cambridge University Press), 55
Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984,
Reviews of Modern Physics, 56, 255
Bernstein, J.P., & Hughers, P.A. 2006, astro-ph/0606012
Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
Einfeldt, B., Munz, C.D., Roe, P.L., and Sjögreen, B. 1991,
J. Comput. Phys., 92, 273
Del Zanna, L., Bucciantini, N., & Londrillo, P. 2003, As-
tronomy & Astrophysics, 400, 397 (dZBL)
http://arxiv.org/abs/astro-ph/0606012
12 A. Mignone and J.C. McKinney
De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ,
599, 1238
Di Matteo, T., Perna, R., & Narayan, R. 2002, ApJ, 579,
Duncan, G. C., & Hughes, P. A. 1994, ApJL, 436, L119
Duncan, C., Hughes, P., & Opperman, J. 1996, ASP
Conf. Ser. 100: Energy Transport in Radio Galaxies and
Quasars, 100, 143
Falle, S. A. E. G., & Komissarov, S. S. 1996, mnras, 278,
Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ,
589, 444
Giacomazzo, B., & Rezzolla, L. 2005, J. Fluid Mech., xxx
Hardee, P. E., Rosen, A., Hughes, P. A., & Duncan, G. C.
1998, ApJ, 500, 599
Harten, A., Lax, P.D., and van Leer, B. 1983, SIAM Re-
view, 25(1):35,61
Hirose S., Krolik J. H., De Villiers J.-P., Hawley J. F., 2004,
ApJ, 606, 1083
Jeffrey A., Taniuti T., 1964, Non-linear wave propagation.
Academic Press, New York
Kohri, K., & Mineshige, S. 2002, ApJ, 577, 311
Kohri, K., Narayan, R., & Piran, T. 2005, ApJ, 629, 341
Koide, S. 1997, ApJ, 478, 66
Komissarov, S. S. 1997, Phys. Lett. A, 232, 435
Komissarov, S. S. 1999, mnras, 308, 1069
Leismann, T., Antón, L., Aloy, M. A., Müller, E., Mart́ı,
J. M., Miralles, J. A., & Ibáñez, J. M. 2005, Astronomy
& Astrophysics, 436, 503
Lichnerowicz, A. 1976, Journal of Mathematical Physics,
17, 2135
Mart́ı, J. M. & Müller, E. 2003, Living Reviews in Relativ-
ity, 6, 7
Mart́ı, J. M. A., Müller, E., Font, J. A., Ibáñez, J. M. A.,
& Marquina, A. 1997, ApJ, 479, 151
Mathews, W. G. 1971, ApJ, 165, 147
Meier, D. L., Koide, S., & Uchida, Y. 2001, Science, 291,
Meliani, Z., Sauty, C., Tsinganos, K., & Vlahakis, N. 2004,
Astronomy & Astrophysics, 425, 773
McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977
McKinney, J. C. 2005, ApJL, 630, L5
McKinney, J. C. 2006, MNRAS, 368, 1561
McKinney J. C., Narayan R., 2006a, MNRAS, in press
(astro-ph/0607575)
McKinney J. C., Narayan R., 2006b, MNRAS, in press
(astro-ph/0607576)
Mignone, A., Plewa, T., and Bodo, G. 2005, ApJS, 160,
Mignone, A., Massaglia, S., & Bodo, G. 2005, Space Science
Reviews, 121, 21
Mignone, A., & Bodo, G. 2006, MNRAS, 368, 1040
Mignone, A., Bodo, G., Massaglia, S., Matsakos, T.,
Tesileanu, O., Zanni, C., and A. Ferrari 2006, accepted
for publication on ApJ.
Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Grav-
itation, San Francisco: W.H. Freeman and Co., 1973
Mizuta, A., Yamada, S., & Takabe, H. 2004, ApJ, 606, 804
Nishikawa, K.-I., Koide, S., Sakai, J.-I., Christodoulou,
D. M., Sol, H., & Mutel, R. L. 1997, ApJL, 483, L45
Noble, S. C., Gammie, C. F., McKinney, J. C., & Del
Zanna, L. 2006, ApJ, 641, 626
Popham, R., Woosley, S. E., & Fryer, C. 1999, pJ, 518, 356
Ryu, D., Chattopadhyay, I., & Choi, E. 2006, ApJS, 166,
Scheck, L., Aloy, M. A., Mart́ı, J. M., Gómez, J. L., Müller,
E. 2002, MNRAS, 331, 615
Synge, J. L. 1957, The relativistic Gas, North-Holland Pub-
lishing Company
Taub, A. H. 1948, Physical Review, 74, 328
Tchekhovskoy, A, McKinney, J. C.,& Narayan, R. 2006,
MNRAS, submitted
Toro, E. F. 1997, Riemann Solvers and Numerical Methods
for Fluid Dynamics, Springer-Verlag, Berlin
van Putten, M. H. P. M. 1993, ApJL, 408, L21
APPENDIX A: PRIMITIVE VARIABLE
INVERSION SCHEME
We outline a new primitive variable inversion scheme that
is used to convert the evolved conserved quantities into so-
called primitive quantities that are necessary to obtain the
fluxes used for the evolution. This scheme allows a gen-
eral EoS by only requiring specification of thermodynamical
quantities and it also avoids catastrophic cancellation in the
non-relativistic and ultrarelativistic limits. Large Lorentz
factors (up to 106) may not be uncommon in some astro-
physical contexts (e.g. Gamma-Ray-Burst) and ordinary in-
version methods can lead to severe numerical problems such
as effectively dividing by zero and subtractive cancellation,
see, for instance, Bernstein & Hughes (2006).
First, we note that the general relativistic conserva-
tive quantities can be written more like special relativis-
tic quantities by choosing a special frame in which to mea-
sure all quantities. A useful frame is the zero angular mo-
mentum (ZAMO) observer in an axisymmetric space-time.
See Noble et al. (2006) for details. From their expressions,
it is useful to note that catastrophic cancellations for non-
relativistic velocities can be avoided by replacing γ − 1 in
any expression with (uαu
α)/(γ + 1), where here uα is the
relative 4-velocity in the ZAMO frame. From here on the
expressions are in the ZAMO frame and appear similar to
the same expressions in special relativity.
A1 Inversion Procedure
Numerical integration of the conservation law (1) proceeds
by evolving the conservative state vector U = (D, m, B, E)
in time. Computation of the fluxes, however, requires veloc-
ity and pressure to be recovered from U by inverting Eqns.
(5)–(7), a rather time consuming and challenging task. For
the constant-Γ law, a recent work by Noble et al. (2006)
examines several methods of inversion. In this section we
discuss how to modify the equations of motion, interme-
diate calculations, and the inversion from conservative to
primitive quantities so that the RMHD method 1) permits
a general EoS; and 2) avoids catastrophic cancellations in
the non-relativistic and ultrarelativistic limits.
Our starting relations are the total energy density (7),
E = W − p + 1 + |v|
|B|2 − S
, (A1)
and the square modulus of Eq. (6),
http://arxiv.org/abs/astro-ph/0607575
http://arxiv.org/abs/astro-ph/0607576
Equation of state in RMHD 13
|m|2 =
W + |B|2
)2 |v|2 − S
2W + |B|2
, (A2)
where S ≡ m · B and W = Dhγ. Note that in order for
this expression to be accurate in the non-relativistic limit,
one should analytically cancel any appearance of E in this
expression. Eq. (A2) can be inverted to express the square
of the velocity in terms of the only unknown W :
|v|2 = S
2(2W + |B|2) + |m|2W 2
(W + |B|2)2W 2
. (A3)
After inserting (A3) into (A1) one has:
E = W − p + |B|
|B|2|m|2 − S2
2(|B|2 + W )2 . (A4)
In order to avoid numerical errors in the non-relativistic
limit one must modify the equations of motion and sev-
eral intermediate calculations. One solves the conservation
equations with the mass density subtracted from the en-
ergy by defining a new conserved quantity (E′ = E − D)
and similarly for the energy flux. In addition, operations
based upon γ can lead to catastrophic cancellations since
the residual γ − 1 is often requested and is dominant in the
non-relativistic limit. A more natural quantity to consider
is |v|2 or γ2|v|2. Also, in the ultrarelativistic limit calcula-
tions based upon γ(|v|2) have catastrophic cancellation er-
rors when |v| → 1. This can be avoided by 1) using instead
|u|2 ≡ γ2|v|2 and 2) introducing the quantities E′ = E −D
and W ′ = W − D, with W ′ properly rewritten as
D|u|2
1 + γ
to avoid machine accuracy problems in the nonrelativistic
limit, where χ ≡ ρǫ+p. Thus our relevant equations become:
′ − p + |B|
|B|2|m|2 − S2
2(|B|2 + W ′ + D)2 , (A6)
|m|2 = (W + |B|2)2 |u|
1 + |u|2 −
2W + |B|2
, (A7)
where W = W ′ + D.
Equations (A6) and (A7) may be inverted to find W ′,
p and |u|2. A one dimensional inversion scheme is derived
by regarding Eq. (A6) as a single nonlinear equation in the
only unknown W ′ and using Eq. (A7) to express |u|2 as a
function of W ′. Using Newton’s iterative scheme as our root
finder, one needs to compute the derivative
= 1 − dp
|B|2|m|2 − S2
(|B|2 + W ′ + D)3 . (A8)
The explicit form of dp/dW ′ depends on the particular EoS
being used. While prior methods in principle allow for a
general EoS, one has to re-derive many quantities that in-
volve kinematical expressions. This can be avoided by split-
ting the kinematical and thermodynamical quantities. This
also allows one to write the expressions so that there is no
catastrophic cancellations in the non-relativistic or ultrarel-
ativistic limits. Assuming that p = p(χ, ρ), we achieve this
by applying the chain rule to the pressure derivative:
. (A9)
Partial derivatives involving purely thermodynamical quan-
tities must now be supplied by the EoS routines. Derivatives
with respect to W ′, on the other hand, involve purely kine-
matical terms and do not depend on the choice of the EoS.
Relevant expressions needed in our computations are given
in the Appendix.
Once W ′ has been determined to some accuracy, the
inversion process is completed by computing the velocities
from an inversion of equation (6) to obtain
W + |B|2
, (A10)
One then computes χ from an inversion of equation (A5) to
obtain
− D|u|
(1 + γ)γ2
, (A11)
from which p or ρǫ can be obtained for any given EoS. The
rest mass density is obtained from
, (A12)
and the magnetic field is trivially inverted.
In summary, we have formulated an inversion scheme
that 1) allows a general EoS without re-deriving kinematical
expressions; and 2) avoids catastrophic cancellation in the
non-relativistic and ultrarelativistic limits. This inversion in-
volves solving a single non-linear equation using, e.g., a one-
dimensional Newton’s method. A similar two-dimensional
method can be easily written with the same properties, and
such a method may be more robust in some cases since the
one-dimensional version described here involves more com-
plicated non-linear expressions.
One can show analytically that the inversion is accurate
in the ultrarelativistic limit as long as γ.ǫ
machine for γ and
p/(ργ2)&ǫmachine for pressure, where ǫmachine ≈ 2.2 × 10−16
for double precision. The method used by Noble et al. (2006)
requires γ.ǫ
machine/10 due to the repeated use of the ex-
pression γ = 1/
1 − v2 in the inversion. Note that we
use γ =
1 + |u|2 that has no catastrophic cancellation.
The fundamental limit on accuracy is due to evolving en-
ergy and momentum separately such that the expression
E − |m| appears in the inversion. Only a method that
evolves this quantity directly (e.g. for one-dimensional prob-
lems one can evolve the energy with momentum subtracted)
can reach higher Lorentz factors. An example test problem
is the ultrarelativistic Noh test in Aloy et al. (1999) with
p = 7.633 × 10−6, Γ = 4/3, 1 − v = 10−11 (i.e. γ = 223607)
This test has p/(ργ2) ≈ 1.6 × 10−16, which is just below
double precision and so the pressure is barely resolved in
the pre-shock region. The post-shock region is insensitive to
the pre-shock pressure and so is evolved accurately up to
γ ≈ 6 × 107. These facts are have been also confirmed nu-
merically using this inversion within HARM. Using the same
error measures as in Aloy et al. (1999) we can evolve their
test problem with an even higher Lorentz factor of γ = 107
and obtain similar errors of .0.1%.
A2 Kinematical and Thermodynamical
Expressions
The kinematical terms required in equation (A9) may be
easily found from the definition of W ′,
′ ≡ Dhγ − D = D(γ − 1) + χγ2 , (A13)
14 A. Mignone and J.C. McKinney
by straightforward differentiation. This yields
(D + 2γχ)
d|v|2
, (A14)
d(1/γ)
= −Dγ
d|v|2
, (A15)
where
d|v|2
= − 2
3W (W + |B|2) + |B|4
+ |m|2W 3
(W + |B|2)3
, (A16)
is computed by differentiating (A3) with respect to W (note
that d/dW ′ ≡ d/dW ). Equation (A14) does not depend on
the knowledge of the EoS.
Thermodynamical quantities such as ∂p/∂χ, on the
other hand, do require the explicit form of the EoS. For
the ideal gas EoS one simply has
p(χ, ρ) =
Γ − 1
χ , (A17)
where χ = ρǫ+ p. By taking the partial derivatives of (A17)
with respect to χ (keeping ρ constant) and ρ (keeping χ
constant) one has
Γ − 1
= 0 . (A18)
For the TM EoS, one can more conveniently rewrite
(14) as
3p(ρ + χ − p) = (χ − p)(χ + 2ρ − p) , (A19)
which, upon differentiation with respect to χ (keeping ρ con-
stant) yields
2χ + 2ρ − 5p
5ρ + 5χ − 8p
. (A20)
Similarly, by taking the derivative with respect to ρ at con-
stant χ gives
2χ − 5p
5ρ + 5χ − 8p
. (A21)
In order to use the above expressions and avoid catas-
trophic cancellation in the non-relativistic limit, one must
solve for the gas pressure as functions of only ρ and χ and
then write the pressure that explicitly avoids catastrophic
cancellation as {χ, p} → 0. One obtains:
p(χ, ρ) =
2χ(χ + 2ρ)
5(χ + ρ) +
9χ2 + 18ρχ + 25ρ2
. (A22)
Also, for setting the initial conditions it is useful to be able
to convert from a given pressure to the internal energy by
using
ρǫ(ρ, p) =
9p2 + 4ρ2
, (A23)
which also avoids catastrophic cancellation in the non-
relativistic limit.
A3 Newton-Raphson Scheme
Equation (A6) may be solved using a Newton-Raphson it-
erative scheme, where the (k + 1)-th approximation to the
W ′ is computed as
′(k+1)
′(k) − f(W
df(W ′)/dW ′
W ′=W ′(k)
, (A24)
where
) = W
′ − E′ − p + |B|
|B|2|m|2 − S2
2(|B|2 + W ′ + D)2
, (A25)
and df(W ′)/dW ′ ≡ dE′/dW ′ is given by Eq. (A8).
The iteration process terminates when the residual
∣W ′(k+1)/W ′(k) − 1
∣ falls below some specified tolerance.
We remind the reader that, in order to start the iter-
ation process given by (A24), a suitable initial guess must
be provided. We address this problem by initializing, at the
beginning of the cycle, W ′(0) = W̃+ − D, where W̃+ is the
positive root of
P(W, 1) = 0 , (A26)
and P(W, |v|) is the quadratic function
P(W, |v|) = |m|2−|v|2W 2+(2W+|B|2)(2W+|B|2−2E) .(A27)
This choice guarantees positivity of pressure, as it can be
proven using the relation
P(W, |v|)
2(2W + |B|2) , (A28)
which follows upon eliminating the (S/W )2 term in Eq. (A2)
with the aid of Eq. (A1). Seeing that P(W, |v|) is a con-
vex quadratic function, the condition p > 0 is equivalent to
the requirement that the solution W must lie outside the
interval [W−, W+], where P(W±, |v|) = 0. However, since
P(W, |v|) > P(W, 1), it must follow that W̃+ > W+ and
thus W̃+ lies outside the specified interval. We tacitly as-
sume that the roots are always real, a condition that is al-
ways met in practice.
	Introduction
	Relativistic MHD Equations
	Equations of Motion
	Equation of State
	Propagation of Fast Magneto-sonic Shocks
	Purely Hydrodynamical Shocks
	Magnetized Shocks
	Numerical Simulations
	Relativistic Blast Waves
	Propagation of Relativistic Jets
	Magnetized Accretion near Kerr Black Holes
	Conclusions
	Primitive Variable Inversion Scheme
	Inversion Procedure
	Kinematical and Thermodynamical Expressions
	Newton-Raphson Scheme

  The role of the equation of state for a perfectly conducting, relativistic
magnetized fluid is the main subject of this work. The ideal constant
$\Gamma$-law equation of state, commonly adopted in a wide range of
astrophysical applications, is compared with a more realistic equation of state
that better approximates the single-specie relativistic gas. The paper focus on
three different topics. First, the influence of a more realistic equation of
state on the propagation of fast magneto-sonic shocks is investigated. This
calls into question the validity of the constant $\Gamma$-law equation of state
in problems where the temperature of the gas substantially changes across
hydromagnetic waves. Second, we present a new inversion scheme to recover
primitive variables (such as rest-mass density and pressure) from conservative
ones that allows for a general equation of state and avoids catastrophic
numerical cancellations in the non-relativistic and ultrarelativistic limits.
Finally, selected numerical tests of astrophysical relevance (including
magnetized accretion flows around Kerr black holes) are compared using
different equations of state. Our main conclusion is that the choice of a
realistic equation of state can considerably bear upon the solution when
transitions from cold to hot gas (or viceversa) are present. Under these
circumstances, a polytropic equation of state can significantly endanger the
solution.

Introduction
	Relativistic MHD Equations
	Equations of Motion
	Equation of State
	Propagation of Fast Magneto-sonic Shocks
	Purely Hydrodynamical Shocks
	Magnetized Shocks
	Numerical Simulations
	Relativistic Blast Waves
	Propagation of Relativistic Jets
	Magnetized Accretion near Kerr Black Holes
	Conclusions
	Primitive Variable Inversion Scheme
	Inversion Procedure
	Kinematical and Thermodynamical Expressions
	Newton-Raphson Scheme
