THE ROLE OF THE EQUATION OF STATE IN RESISTIVE RELATIVISTIC MAGNETOHYDRODYNAMICS

Magnetic fields play an important role in determining the evolution of the matter in many astrophysical objects. In highly conducting plasma, the magnetic field can be amplified by gas contraction or shear motion. Even when the magnetic field is weak initially, the magnetic field can grow rapidly and influence the gas dynamics of the system. This is particularly important for the high-energy astrophysical phenomena related to strongly magnetized relativistic plasmas associated with objects such as active galactic nuclei (AGNs; e.g., Urry & Padovani 1995), relativistic jets (e.g., Mirabel & Rodríguez 1999; Blandford 2002), pulsar winds (e.g., Gaensler & Slane 2006; Kirk et al. 2009), gamma-ray bursts (GRBs; Zhang & Mészáros 2004; Piran 2005; Mészáros 2006), and magnetars (e.g., Woods & Thompson 2006; Mereghetti 2008).

The ideal magnetohydrodynamic (MHD) approximation is a good description of the global properties and dynamics of such systems well into their nonlinear regimes. In this limit, the electrical resistivity η = 1/σ vanishes (infinite electrical conductivity). In this framework, many multidimensional ideal relativistic MHD (RMHD) codes have been developed to investigate relativistic astrophysical phenomena including fully nonlinear regimes (e.g., Komissarov 1999, 2001; Koide et al. 1999; Koldoba et al. 2002; Del Zanna et al. 2003, 2007; Leismann et al. 2005; Gammie et al. 2003; De Villiers & Hawley 2003; Anninos et al. 2005; Duez et al. 2005; Shibata & Sekiguchi 2005; Antón et al. 2006; Mignone & Bodo 2006; Mizuno et al. 2006; Neilsen et al. 2006; Giacomazzo & Rezzolla 2007; Farris et al. 2008; Mignone et al. 2009; Beckwith & Stone 2011; Inoue et al. 2011). The ideal MHD limit provides a convenient form for solving the equations of RMHD and is also an excellent approximation for many relativistic astrophysical phenomena. However, in extreme cases such as binary mergers (the merger of two neutron stars or of a neutron star with a black hole; e.g., Anderson et al. 2008; Shibata & Taniguchi 2011; Faber & Rasio 2012) or the central engine of long GRBs (collapsar; e.g., MacFadyen & Woosley 1999) the electrical conductivity can be small, and regions of high resistivity may appear.

Quite often numerical simulations using ideal RMHD exhibit violent magnetic reconnection. The magnetic reconnection observed in ideal RMHD simulations is due to purely numerical resistivity, occurs as a result of truncation errors, and hence fully depends on details of the numerical scheme and resolution. Magnetic reconnection is one of the most important phenomena in astrophysics. It is highly dynamic, and it converts magnetic energy into fluid energy. The magnetic reconnection process has been invoked to explain flaring events (e.g., Lyutikov 2006; Giannios et al. 2009) and magnetic annihilation (Coroniti 1990; Lyubarsky & Kirk 2001) in relativistic plasmas. Therefore, numerical codes solving the resistive RMHD (RRMHD) equations and that allow control of magnetic reconnection according to a physical model of resistivity are highly desirable.

Numerical simulation using the ideal RMHD equations is considerably easier than using the RRMHD equations because the equations become mixed hyperbolic with stiff relaxation terms. The pioneering work on RRMHD done by Komissarov (2007) solved the numerical flux by using the Harten–Lax–van-Leer (HLL) approximate Riemann solver and by using Strang's splitting technique for the stiff relaxation terms. More recently, Palenzuela et al. (2009) have proposed a numerical method that solves the stiff relaxation terms in the equations by an implicit–explicit (IMEX) Runge–Kutta method, and Dionysopoulou et al. (2012) have extended the work of Palenzuela et al. (2009) to three dimensions. A different approach has been taken by Dumbser & Zanotti (2009) who have applied the high-order P_NP_M scheme to solving the RRMHD equations, and also Takamoto & Inoue (2011) who have used the method of characteristics to solve the Maxwell equations accurately. Even more recently, a 3+1 resistive general relativistic MHD (GRMHD) code using mean-field dynamo closure has been developed by Bucciantini & Del Zanna (2013).

Plasma in the relativistic regime can have three major characteristics: the system has (1) relativistic fluid velocity (kinetic energy much greater than rest mass energy), (2) relativistic temperature (internal energy much greater than rest mass energy), or (3) relativistic Alfvén speed (magnetic energy much greater than rest mass energy). The second characteristic of relativistic temperature brings us to the issue of the equation of state (EoS) of the plasma. The EoS most commonly used in RMHD simulations is designed for plasmas with constant specific heat ratio (the so-called ideal EoS). However, this ideal EoS is valid only for plasmas with either ultrarelativistic temperature or non-relativistic temperature. The theory of relativistic perfect gases (Synge 1957) has shown that the specific heat ratio cannot be constant if consistency with kinetic theory is required. However, the exact EoS involves modified Bessel functions and is too complicated to be efficiently implemented in numerical codes. To get around this problem, Mignone et al. (2005) introduced an approximate EoS given by a simple analytical formulation in the context of relativistic non-magnetized flows. This approximate EoS was applied in the context of RMHD by Mignone & McKinney (2007). A different EoS approximation than that proposed by Mignone et al. (2005) has been proposed by Ryu et al. (2006). Clearly, a determination of the effects of a difference in the choice of the approximate EoS is important to further advances in RRMHD.

In this paper, we present the development of a new RRMHD simulation code including different realistic EoS approximations such as those proposed by Mignone et al. (2005) or by Ryu et al. (2006). This new RRMHD code is based on the ideal RMHD code RAISHIN (Mizuno et al. 2006, 2011) which uses a Godunov-type scheme to solve the conservation equations of ideal RMHD. In particular, we apply this new code to the role of the EoS in the RRMHD regime. We describe the basic equations of RRMHD in Section 2, three different EoSs are investigated in Section 3, and the numerical methods are described in Section 4. The various numerical tests in one dimension and multi dimensions are presented in Section 5. In Section 6 we conclude.

We have considered n^α to be the time-like translational killing vector field in a flat (Minkowski) space-time, so n^α = (− 1, 0, 0, 0), where we use Greek letters that take values from 0 to 3 for the indices of four-dimensional (4D) space-time tensors, while Roman letters take values from 1 to 3 for the indices of three-dimensional (3D) spatial tensors. We use the speed of light c = 1 and Lorentz–Heaviside notation for electromagnetic quantities, so that all $\sqrt{4} \pi$ factors disappear.

The total energy momentum tensor T^αβ describing a perfect fluid coupled to an electromagnetic field is defined as

$$\begin{equation} T^{\alpha \beta }= T^{\alpha \beta }_{{\rm fluid}} + T^{\alpha \beta }_{{\rm EM}}. \end{equation} \tag{ 1 }$$

The first term is due to matter:

$$\begin{equation} T^{\alpha \beta }_{{\rm fluid}}= \rho hu^{\alpha }u^{\beta }+pg^{\alpha \beta }, \end{equation} \tag{ 2 }$$

where u^α is the fluid four velocity, while h (= 1 + + p/ρ), ρ, p, and are the enthalpy, the proper rest mass density, the gas pressure, and the specific internal energy as measured in the fluid rest frame. The second term comes from the electromagnetic field:

$$\begin{equation} T^{\alpha \beta }_{{\rm EM}} = F^{\alpha \mu } F_{\mu }^{\beta } - {1 \over 4} \left(F^{\mu \nu }F_{\mu \nu }\right) g^{\alpha \beta }, \end{equation} \tag{ 3 }$$

where F^αβ and its dual *F^αβ are the Maxwell and Faraday tensors of the electromagnetic field given by

$$\begin{eqnarray} && F^{\alpha \beta } = n^{\alpha } E^{\beta } - n^{\beta } E^{\alpha } + n_{\nu } e^{\nu \alpha \beta \mu } B_{\mu }, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} && ^*F^{\alpha \beta } = n^{\alpha } B^{\beta } - n^{\beta } B^{\alpha } + n_{\nu } e^{\nu \alpha \beta \mu } E_{\mu }. \end{eqnarray} \tag{ 5 }$$

E^α and B^α are the electric and magnetic fields as measured by an observer moving along any time-like vector n^α, while $e^{\alpha \beta \mu \nu } = \sqrt{\vphantom{A^A}\smash{\hbox{$-g$}}} \epsilon _{\alpha \beta \mu \nu }$ is the Levi-Civita alternating tensor of space-time and _αβμν is the 4D Levi-Civita symbol.

In the global inertial frame with time-independent coordinate grid, the full system of Euler and Maxwell's equations is

$$\begin{eqnarray} && \partial _{t} D + \nabla \cdot D{\boldsymbol v} =0, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} && \partial _{t} {\boldsymbol m} + \nabla \cdot {\boldsymbol \Pi } = 0, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} && \partial _{t} \tau + \nabla \cdot {\boldsymbol Y} =0, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} && \partial _{t} {\boldsymbol E} - \nabla \times {\boldsymbol B} + \nabla \Psi = {-}{\boldsymbol j}, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} && \partial _{t} {\boldsymbol B} + \nabla \times {\boldsymbol E} + \nabla \Phi = {\boldsymbol 0}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} && \partial _{t} \Psi + \nabla \cdot {\boldsymbol E} = q - \kappa \Psi , \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} && \partial _{t} \Phi + \nabla \cdot {\boldsymbol B} = -\kappa \Phi , \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} && \partial _{t} q + \nabla \cdot {\boldsymbol j} = 0, \end{eqnarray} \tag{ 13 }$$

where j is the spatial current vector, q is the charge density, κ is the damping rate parameter, and the conserved variables

$$\begin{eqnarray} D &=& \rho \gamma , \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} {\boldsymbol m} &=& \rho h \gamma ^{2} {\boldsymbol v} + {\boldsymbol E} \times {\boldsymbol B}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \tau &=& \rho h \gamma ^{2} - p + {1 \over 2}(E^{2} + B^{2}) \end{eqnarray} \tag{ 16 }$$

express the relativistic mass density, the momentum density, and the total energy density. Here, v is the velocity measured by an inertial observer and $\gamma \equiv 1/\sqrt{\vphantom{A^A}\smash{\hbox{$1-v^{2}$}}}$ is the Lorentz factor. The energy flux density and the momentum flux density can then be given by

$$\begin{eqnarray} {\boldsymbol Y} &=& \rho h \gamma ^{2} {\boldsymbol v} + {\boldsymbol E} \times {\boldsymbol B}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} {\boldsymbol \Pi } &=& - {\boldsymbol E}\!{\boldsymbol E} - {\boldsymbol B}\!{\boldsymbol B} + \rho h \gamma ^{2} {\boldsymbol v}{\boldsymbol v} + \left[ {1 \over 2} (E^{2}+B^{2}) + p \right] {\boldsymbol g}.\quad \end{eqnarray} \tag{ 18 }$$

An EoS is needed to close the system, so we have adopted a variable EoS (e.g., Mignone et al. 2005; Mignone & McKinney 2007; Ryu et al. 2006). The details of the variable EoS are explained in next section.

Equations (9)–(12) evolve the augmented Maxwell's equations which contain two additional fields Ψ and Φ to control the system dynamics. In this approach, the two scalar fields Ψ and Φ indicate deviations of the divergence of the electric and magnetic fields from the values prescribed by Maxwell's equations, propagate at the speed of light, and decay exponentially over a timescale ∼1/κ when the damping rate parameter κ > 0. Following previous studies (Komissarov 2007; Palenzuela et al. 2009; Dumbser & Zanotti 2009), we have adopted the so-called hyperbolic divergence-cleaning approach used in the context of ideal MHD (Dedner et al. 2002).

The system of Equations (6)–(13) is closed by means of Ohm's law. Ohm's law for relativistic plasmas can be very complicated (e.g., Lichnerowicz 1967; Ardavan 1976; Blackman & Field 1994; Gedalin 1996; Melatos & Melrose 1996; Punsley 2001; Meier 2004). In this paper, we consider only the simplest kind of relativistic Ohm's law which assumes an isotropic plasma resistivity (e.g., Komissarov 2007; Palenzuela et al. 2009; Zenitani et al. 2010; Takamoto & Inoue 2011; Takahashi et al. 2011). In covariant form, the four vector of the electric current is obtained from

$$\begin{equation} I^{\alpha } = \sigma F^{\alpha \beta } u_{\beta } + q_{0} u^{\alpha }, \end{equation} \tag{ 19 }$$

where σ = 1/η is the conductivity, η is the resistivity, and q₀ = −I_αu^α is the electric charge density as measured in the fluid flame (Lichnerowicz 1967; Blackman & Field 1994). In a special relativistic inertial frame, we find

$$\begin{equation} {\boldsymbol j} = \sigma \gamma [ {\boldsymbol E} + {\boldsymbol v} \times {\boldsymbol B} - ({\boldsymbol E} \cdot {\boldsymbol v}){\boldsymbol v}] + q {\boldsymbol v}. \end{equation} \tag{ 20 }$$

In the fluid rest frame, this equation becomes

$$\begin{equation} {\boldsymbol j} = \sigma {\boldsymbol E}. \end{equation} \tag{ 21 }$$

The ideal MHD limit of Ohm's law is given by the limit of infinite conductivity (σ → ∞). In this limit, Equation (20) reduces to

$$\begin{equation} {\boldsymbol E}+{\boldsymbol v} \times {\boldsymbol B} - ({\boldsymbol E} \cdot {\boldsymbol v}) {\boldsymbol v} = 0. \end{equation} \tag{ 22 }$$

Splitting this equation into components normal and parallel to the velocity vector gives

$$\begin{eqnarray} && \partial _{t} {\boldsymbol E}_{\parallel } + \sigma \gamma [{\boldsymbol E}_{\parallel } - ({\boldsymbol E} \cdot {\boldsymbol v}) {\boldsymbol v}] =0, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} && \partial _{t} {\boldsymbol E}_{\perp } + \sigma \gamma [{\boldsymbol E}_{\perp } + {\boldsymbol v} \times {\boldsymbol B}] = 0. \end{eqnarray} \tag{ 24 }$$

From these equations, one obtains the well-known ideal MHD condition

$$\begin{equation} {\boldsymbol E} = - {\boldsymbol v} \times {\boldsymbol B}. \end{equation} \tag{ 25 }$$

In this limit, the electric field is orthogonal to both magnetic and velocity fields.

An EoS relating thermodynamic quantities is required to close the system of Equations (6)–(13). In general, an EoS is written as

$$\begin{equation} h \equiv h(p, \rho), \end{equation} \tag{ 26 }$$

and general forms for the polytropic index n and the sound speed c_s are given by

$$\begin{equation} n=\rho {\partial h \over \partial \rho } -1,\quad \ c_{s}^{2} = - {\rho \over nh} {\partial h \over \partial \rho }. \end{equation} \tag{ 27 }$$

The most commonly used EoS, a constant Γ law (ideal) EoS, is given by

$$\begin{equation} h= 1+ {\Gamma \over \Gamma -1} \Theta , \end{equation} \tag{ 28 }$$

where Γ is the constant specific heat ratio and Θ = p/ρ is the temperature. The sound speed is calculated from

$$\begin{equation} c_{s}^{2} = {\Gamma \Theta \over h}. \end{equation} \tag{ 29 }$$

The constant Γ law EoS may be applied correctly to a plasma with non-relativistic temperature where Γ = 5/3 or to a plasma with an ultrarelativistic temperature where Γ = 4/3. However, in the high-temperature limit, i.e., Θ → ∞ with Γ > 4/3, the sound speed exceeds relativistic limit ($c_{s} >1/\sqrt{3}$). Moreover, a constant Γ law EoS is not consistent with relativistic kinetic theory, the so-called Taub's fundamental inequality, which requires the specific enthalpy to satisfy:

$$\begin{equation} (h-\Theta)(h-4\Theta) \ge 1. \end{equation} \tag{ 30 }$$

This rules out a constant Γ law EoS with Γ > 4/3, if applied to 0 < Θ < ∞.

The theory of single-component perfect gases in the relativistic regime shows that the specific enthalpy is a function of the temperature Θ = p/ρ only and has the form (Synge 1957)

$$\begin{equation} h = {K_{3} (1/\Theta) \over K_{2} (1/\Theta)}, \end{equation} \tag{ 31 }$$

where K₂ and K₃ are the second- and third-order modified Bessel functions of the second kind, respectively. Using an equivalent Γ_eq = (h − 1)/(h − 1 − Θ) in the non-relativistic temperature limit (Θ → 0) yields Γ_eq → 5/3, and in the ultrarelativistic temperature limit (Θ → ∞) yields Γ_eq → 4/3 (see Figure 1(a)). However, this EoS requires extra computational costs because the thermodynamics of the fluid is expressed in terms of the modified Bessel functions (Falle & Komissarov 1996).

Zoom In Zoom Out Reset image size

Recently, Mignone et al. (2005) proposed an EoS, the so-called TM EoS, that follows Equation (31) well. The TM EoS, which was first introduced by Mathews (1971), is given by

$$\begin{equation} p={\rho \epsilon (\rho \epsilon + 2 \rho) \over 3 (\rho \epsilon + \rho)} \quad{\rm or }\; h={5 \over 2} \Theta + \sqrt{{9 \over 4} \Theta ^{2} +1}, \end{equation} \tag{ 32 }$$

and the sound speed is calculated from

$$\begin{equation} c_{s}^{2} = {\Theta \over 3h} {5h - 8 \Theta \over h-\Theta }. \end{equation} \tag{ 33 }$$

The TM EoS corresponds to the lower bound of Taub's fundamental inequality, i.e., (h − Θ)(h − 4Θ) = 1, and produces the correct asymptotic values for Γ_eq.

Ryu et al. (2006) proposed an EoS which is a simpler algebraic function of Θ, hereafter referred to as the RC EoS, that satisfies Taub's inequality for all Θ. The RC EoS is given by

$$\begin{equation} {p \over \epsilon - p} = {3p + 2 \rho \over 9 p + 3 \rho } {\rm or } h=2{6 \Theta ^{2} + 4 \Theta + 1 \over 3 \Theta +2}, \end{equation} \tag{ 34 }$$

and the sound speed is calculated from

$$\begin{equation} c_{s}^{2} = {\Theta (3\Theta +2)(18 \Theta ^{2} + 24 \Theta + 5) \over 3 (6 \Theta ^{2} + 4 \Theta +1)(9 \Theta ^{2}+ 12 \Theta +2)}. \end{equation} \tag{ 35 }$$

Figure 1 shows the equivalent Γ, the specific enthalpy, and the sound speed as a function of Θ for TM EoS, RC EoS, Synge's EoS, as well as a constant Γ law EoS with Γ = 5/3 and 4/3. The specific enthalpy and the sound speed computed using the TM EoS and the RC EoS are well matched to Synge's EoS and cannot be distinguished on the plots. The approximations to Synge's EoS such as the TM EoS and RC EoS are hereafter referred to as approximate EoSs.

A well-known and challenging feature of the system of Equations (6)–(13) is that they have source terms for the evolution of the electric field that become stiff in the high-conductivity (low-resistivity) limit. Following Komissarov (2007), we will use the Strang splitting technique (Strang 1968).

The system of Equations (6)–(13) can be written as a single-phase vector equation:

$$\begin{equation} {\partial U(P) \over \partial t } + {\partial F^{m} (P) \over \partial x^{m}} = S(P), \end{equation} \tag{ 36 }$$

where

$$U = \left(\begin{array}{@{}c@{}}D \\ m^{i} \\ \tau \\ E^{i} \\ B^{i} \\ \Psi \\ \Phi \\ q \end{array} \right), \ P=\left(\begin{array}{@{}c@{}}\rho \\ v^{i} \\ p \\ E^{i} \\ B^{i} \\ \Psi \\ \Phi \\ q \end{array} \right), \ S=\left(\begin{array}{@{}c@{}}0 \\ 0^{i} \\ 0 \\ -j^{i} \\ 0^{i} \\ q-\kappa \Psi \\ -\kappa \Phi \\ 0 \end{array} \right)$$

are the vectors of conserved variables, primitive variables, and sources, respectively, and

$$F^{m} = \left(\begin{array}{@{}c@{}}\rho \gamma v^{m} \\ \Pi ^{im} \\ Y^{m} \\ -e^{{\rm imk}} B_{k} + \Psi g^{im} \\ e^{{\rm imk}} E_{k} + \Phi g^{im} \\ E^{m} \\ B^{m} \\ j^{m} \end{array} \right)$$

is the vector of numerical fluxes, where e^ijk is the Levi–Civita alternating tensor of space.

The source term can be split into two parts:

$$S_{a}(P)=\left(\begin{array}{@{}c@{}}0 \\ 0^{i} \\ 0 \\ -qv^{i} \\ 0^{i} \\ q \\ 0 \\ 0 \end{array} \right) {\rm and} \ S_{b}(P)=\left(\begin{array}{@{}c@{}}0 \\ 0^{i} \\ 0 \\ -j_{c}^{i} \\ 0^{i} \\ -\kappa \Psi \\ -\kappa \Phi \\ 0 \end{array} \right),$$

where

$$\begin{equation} {\boldsymbol j}_{c} = \sigma \gamma [{\boldsymbol E} + {\boldsymbol v} \times {\boldsymbol B} - ({\boldsymbol E} \cdot {\boldsymbol v}) {\boldsymbol v}] \end{equation} \tag{ 37 }$$

is the conductivity current. The source term S_b is a stiff relaxation term that requires special care to capture the dynamics in a stable and accurate manner. In the Strang time-step splitting technique, first the solution is advanced using the stiff part of the equations

$$\begin{equation} {\partial U(P) \over \partial t} = S_{b}(P) \end{equation} \tag{ 38 }$$

over the half-time step, Δt/2. Second, advance of the non-stiff part of the equations is made via second-order accurate numerical integration of

$$\begin{equation} {\partial U(P) \over \partial t} + {\partial F^{m} (P) \over \partial x^{m}}= S_{a}(P) \end{equation} \tag{ 39 }$$

over the full time step. Third, again the solution is advanced by the stiff part of the equations over the half-time step.

Time advance of the non-stiff part of the equations is given by

$$\begin{equation} U_{n+1} = U_{n} + \Delta t \sum ^{n_d}_{m=1} {F_{m-1/2, n+1/2} - F_{m+1/2, n+1/2} \over \Delta x^{m}} + \Delta t S_{a, n+1/2}, \end{equation} \tag{ 40 }$$

where U_n represents the cell-centered conserved variables at t = t_n, U_{n + 1} represents the cell-centered conserved variables at t = t_n + Δt, S_{a, n + 1/2} represents the cell-centered non-stiff source term at t = t_n + Δt/2, F_{m + 1/2, n + 1/2} is the numerical flux though the right-hand side cell interface, F_{m − 1/2, n + 1/2} is the numerical flux though the left-hand side cell interface in the direction of x^m at time t = t_n + Δt/2, Δx^m is the cell size in this direction, and n_d is the number of spatial dimensions. The non-stiff sources and numerical fluxes at the half-time step are calculated from

$$\begin{equation} U_{n+1/2} = U_{n} + {\Delta t \over 2} \sum ^{n_d}_{m=1} {F_{m-1/2, n} - F_{m+1/2, n} \over \Delta x^{m}} + {\Delta t \over 2} S_{a, n}. \end{equation} \tag{ 41 }$$

The numerical fluxes at the cell-interface F_{m + 1/2, n} are calculated using the simplified HLL approximate Riemann solver (Harten et al. 1983; Komissarov 2007) where the maximum characteristic speed of the system in each direction equals the speed of light. The left-hand and right-hand states of each cell interface using the HLL approximate Riemann solver are computed from various reconstruction schemes as in our ideal RMHD code (Mizuno et al. 2006, 2011). In this paper, we use a piecewise linear method (PLM) reconstruction such as the minmod slope-limited linear interpolation scheme or the Monotonized Central (MC) slope-limited linear interpolation scheme as these are the simplest reconstruction schemes that capture a shock sharply. The resulting scheme is second-order accurate in time and space.

Following the work of Komissarov (2007), the split evolution equations of the electric field, Equations (23) and (24), can be solved analytically:

$$\begin{eqnarray} {\boldsymbol E}_{\parallel } &=& {\boldsymbol E}_{\parallel }^{0} \exp \left[ - {\sigma \over \gamma } t \right], \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} {\boldsymbol E}_{\perp } &=& {\boldsymbol E}_{\perp }^{*} + ({\boldsymbol E}_{\perp }^{0} - {\boldsymbol E}_{\perp }^{*}) \exp [-\sigma \gamma t], \end{eqnarray} \tag{ 43 }$$

where E*_⊥ = −v × B and suffix 0 indicates the initial component. The stiff part of the equations related to the two scalar fields Ψ and Φ are also solved analytically with solution

$$\begin{eqnarray} \Psi &=& \Psi _{0} \exp [-\kappa t], \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \Phi &=& \Phi _{0} \exp [-\kappa t], \end{eqnarray} \tag{ 45 }$$

where Ψ₀ and Φ₀ are the initial values of Ψ and Φ.

In order to evolve this system of equations, the numerical fluxes F^m must be computed at each time step. These fluxes depend on the primitive variables P, which must be recovered from the evolved conserved variables U. In conserved variables, E and B can be calculated by evolving Maxwell's equations. However, it is more stable to evolve the stiff part of Equations (42)–(43) during the primitive recovery process when σ becomes large (Palenzuela et al. 2009). The primitive recovery procedure adapted to an approximate EoS such as the TM EoS and the RC EoS follows that used by Palenzuela et al. (2009).

In this section, one-dimensional (1D) and two-dimensional (2D) tests are presented. Three 1D tests have been used to validate our new RRMHD code in different regimes. A 1D shock-tube test and three 2D tests have been used to investigate the effect of different EoSs. The damping coefficient of the hyperbolic divergence cleaning is set to κ = 1. The magnetic field is divergence free and charge is preserved at the truncation error level.

5.1. One-dimensional Tests

5.1.1. Large Amplitude CP Alfvén Wave Test

This test consists of the propagation of a large-amplitude circularly polarized Alfvén wave along a uniform background magnetic field B₀ in a domain with periodic boundary conditions. The exact solution is given by Del Zanna et al. (2007) in the ideal MHD limit, and was used as an ideal MHD limit test problem by Palenzuela et al. (2009) and Takamoto & Inoue (2011). Here we use conditions similar to previous studies with

$$\begin{eqnarray} (B_{y}, B_{z}) &=& \zeta _{A} B_{0} (\cos [k(x-v_{A}t], \sin [k(x-v_{A})t]), \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} (v_{y}, v_{z}) &=& - {v_{A} \over B_{0}} (B_{y}, B_{z}), \end{eqnarray} \tag{ 47 }$$

where B_x = B₀, v_x = 0, k is the wave number, and ζ_A is the amplitude of the wave. The special relativistic Alfvén speed v_A is given by

$$\begin{equation} v_{A}^{2} = {2B_{0}^{2} \over h+ B_{0}^{2}\left(1+ \zeta _{A}^{2}\right)} \left(1+ \sqrt{ 1- \left({2 \zeta _{A} B_{0}^{2} \over h + B_{0}^{2} \left(1 + \zeta _{A}^{2}\right)}\right) } \right)^{-1}. \end{equation} \tag{ 48 }$$

For this test, we use initial parameters ρ = p = 1 and B₀ = 0.46188, the Alfvén velocity v_A = 0.25c, and we adopt a constant gamma law EoS with Γ = 2. Following Palenzuela et al. (2009), we use a high uniform conductivity σ = 10⁵ with three different resolutions of 50, 100, and 200 cells covering the computational domain x ∈ [ − 0.5, 0.5]. Figure 2 shows the numerical results at t = 4 (one Alfvén crossing time) for the three different resolutions. This result shows that the new RRMHD code reproduces ideal RMHD solutions when the conductivity σ is high. The L₁ norm errors of the magnetic field component B_y in this test are shown in Figure 3. The numerical result is slightly shallower than second-order convergence. This is likely caused by the periodic boundary.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1.2. One-dimensional Self-similar Current Sheet Test

This test problem has been used for moderate resistivity cases (e.g., Komissarov 2007; Palenzuela et al. 2009; Takamoto & Inoue 2011). In this test problem, the magnetic pressure is much smaller than the gas pressure everywhere. The magnetic field configuration is given by B = [0, B_y(x, t), 0], where B_y(x, t) changes sign within a thin current layer of thickness Δl. An initial solution is provided in equilibrium with p = constant. The evolution of this thin current layer is a slow diffusive expansion due to the resistivity and described by the diffusion equation

$$\begin{equation} \partial _{t} B_{y} - {1 \over \sigma } \partial ^{2}_{x} B_{y} =0. \end{equation} \tag{ 49 }$$

As the thickness of the layer becomes much larger than Δl the expansion becomes self-similar with

$$\begin{equation} B_{y}(x,t) = B_{0}\ {\rm erf} \left( \frac{\sqrt{\sigma}} {2 \chi } \right), \end{equation} \tag{ 50 }$$

where $\chi = \sqrt{t/x}$ and erf is the error function. This analytic result can be used for testing the moderate resistivity regime.

In the test problem, we have chosen an initial solution at t = 1 with p = 50, ρ = 1, E = v = 0, and σ = 100 (η = 1/σ = 0.01). A constant gamma law EoS with Γ = 2 is used. The computational domain is uniform with 200 cells in [ − 1.5, 1.5]. The numerical simulation is evolved up to t = 10 and then the numerical solution is compared in Figure 4 with the analytical solution. The numerical and analytical solutions cannot be distinguished on the plot. This indicates that the moderate resistivity regime is well described by the code.

Zoom In Zoom Out Reset image size

5.1.3. One-dimensional Shock-tube Tests

As a first shock-tube test in the RRMHD regime, we consider a simple MHD version of the Brio and Wu test as in Palenzuela et al. (2009) and Takamoto & Inoue (2011). The initial left and right states are separated at x = 0.5 and are given by

$$\begin{eqnarray} \left(\rho ^{L}, p^{L}, B^{L}_{y}\right) &=& (1.0, 1.0, 0.5) \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} \left(\rho ^{R}, p^{R}, B^{R}_{y}\right) &=& (0.125, 0.1, -0.5). \end{eqnarray} \tag{ 52 }$$

All other fields are set to 0. We use a constant Γ law EoS with Γ = 2. The computational domain covers the region x ∈ [ − 0.5, 0.5] with 200 cells.

Figure 5 shows the numerical results at t = 0.4 for conductivities σ = 0, 10, 10², 10³, and 10⁵. The exact solution to the ideal RMHD Riemann problem was found by Giacomazzo & Rezzolla (2006). When B_x = 0, the solution contains only two fast waves, a left-moving rarefaction wave, and a right-moving shock with a tangential discontinuity between them. The results show that the solution smoothly changes from a wave-like solution for σ = 0 to the ideal-MHD solution for high conductivity σ = 10⁵. Note that for σ = 0 the solution describes a discontinuity propagating at the speed of light corresponding to Maxwell's equations in vacuum. This result is nearly the same as test results for other codes (Palenzuela et al. 2009; Takamoto & Inoue 2011).

Zoom In Zoom Out Reset image size

Palenzuela et al. (2009) reported that Strang's splitting technique became unstable for moderately high conductivity in this shock-tube test, and they suggested using an implicit method. However, Takamoto & Inoue (2011) found that this instability is not related to Strang's splitting technique but instead to the calculation of the electric field during the primitive recovery procedure. In this new RRMHD code, the shock-tube test is solved stably using Strang's splitting technique even when σ ≃ 10⁶.

The Balsara Test 2 (Balsara 2001) is used as a second shock-tube test to investigate the effect of different EoSs in the RRMHD regime. In ideal RMHD, Mignone & McKinney (2007) have already performed this test to check the effect of different EoSs. In this test the initial left and right states are separated at x = 0.5 and are given by

$$\begin{eqnarray} \left(\rho ^{L}, p^{L}, B^{L}_{x}, B^{L}_{y}, B^{L}_{z}\right) &=& (1.0, 30.0, 5.0, 6.0, 6.0) \end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} \left(\rho ^{R}, p^{R}, B^{R}_{x}, B^{R}_{y}, B^{R}_{z}\right) &=& (1.0, 1.0, 5.0, 0.7, 0.7). \end{eqnarray} \tag{ 54 }$$

The computational domain covers the region x ∈ [ − 0.5, 0.5] with 800 cells. This test shows that a mildly relativistic blast wave propagates to the right with maximum Lorentz factor of 1.3 ⩽ γ ⩽ 1.4.

The numerical results at t = 0.4 for the constant Γ law EoS with Γ = 5/3, the TM EoS and the RC EoS using conductivities σ = 0, 10, 10², 10³ are shown in Figures 6, 7, and 8, respectively. The solutions show fast and slow rarefaction waves, a contact discontinuity, and slow and fast shocks from left to right. In these cases, no rotational discontinuity is seen. The results obtained from the TM EoS and RC EoS cases are considerably different from the constant Γ law EoS with Γ = 5/3. In the approximate EoS cases, the rarefaction waves and shocks propagate with smaller velocities. This is predicted from the lower sound speed in the approximate EoS cases relative to overestimated sound speed in the ideal EoS case with Γ = 5/3. Behind the slow shock, the approximate EoS cases have a higher peak density, which follows from the previous considerations. These properties are consistent with those in the ideal RMHD case (Mignone & McKinney 2007). On the other hand, the results obtained from TM EoS and RC EoS cases are very similar at our numerical resolution. This similarity reflects the similarity in the distributions of specific enthalpy (see Figure 1). Again, the results show a smooth change from a wave-like solution for σ = 0 toward an ideal MHD solution for the highest conductivity, σ = 10³, for all the different EoS cases. The differences between the approximate EoS cases and the constant Γ law case are larger at higher σ where the approximate EoS cases approach the ideal MHD case more slowly. Even for low conductivity, i.e., σ = 10, we clearly see a difference between the ideal EoS case with Γ = 5/3 and the approximate EoS cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. Two-dimensional Tests

5.2.1. The Cylindrical Explosion

We now consider tests involving shocks in multi-dimensions. First we choose a test involving a cylindrical blast wave expanding into an initially uniform magnetic field. This is a standard test for ideal RMHD codes even though there is no exact solution because this test will reveal subtle bugs and potential weaknesses in the numerical implementation. For this test, we use a Cartesian computational domain (x, y) ∈ [ − 6, 6] with 200 uniform cells in each direction. The initial explosion is initialized by setting the gas pressure and density to p = 1 and ρ = 0.01 within a cylinder of radius r < 0.8 centered on the origin. In an intermediate region 0.8 < r < 1.0, the pressure and density decrease exponentially to that of the ambient gas which has p = ρ = 0.001. The initial magnetic field is uniform in the x-direction with B = (0.05, 0, 0). The other quantities are set to zero (i.e., v = E = q = 0).

Figure 9 shows the magnetic field components B_x and B_y at t = 4 using the ideal EoS with Γ = 4/3 and the TM EoS. The ideal-MHD simulation is performed using a high conductivity of σ = 10⁵. The results are qualitatively similar to those reported in previous studies in ideal RMHD (Komissarov 1999; Leismann et al. 2005; Neilsen et al. 2006; Del Zanna et al. 2007; Mizuno et al. 2011) and RRMHD (Komissarov 2007; Palenzuela et al. 2009). The results obtained from the ideal EoS with Γ = 4/3 and the TM EoS are qualitatively very similar. This means that an ideal EoS with Γ = 4/3 satisfactorily captures all the shock properties.

Zoom In Zoom Out Reset image size

Figure 10 shows 1D profiles of the gas and magnetic pressure along the y-axis for conductivities of σ = 0, 10, 10², 10³, and 10⁵ using the ideal EoS with Γ = 4/3 and the TM EoS. At high conductivity σ ⩾ 10³, we do not see any significant difference. This means that high conductivity recovers the ideal-MHD solution. As the conductivity decreases, the maximum gas and magnetic pressure decrease. Of course there is no magnetic pressure increase for σ = 0.

Zoom In Zoom Out Reset image size

5.2.2. Kelvin–Helmholtz Instability Test

We present calculations of the linear and nonlinear growth of the 2D Kelvin–Helmholtz instability (KHI) to investigate the effect of conductivity and the EoS on the development of turbulence in the RRMHD regime.

Initial conditions for this test are taken from a combination of previous studies in ideal RMHD (Bucciantini & Del Zanna 2006; Mignone et al. 2009; Beckwith & Stone 2011). The shear velocity profile is given by

$$\begin{equation} v_{x}=\left\lbrace \begin{array}{@{}ll} v_{{\rm sh}} \tanh \left({y-0.5 \over a} \right) & {\mbox{if } y > 0.0} \\\ms -v_{{\rm sh}} \tanh \left({y+0.5 \over a} \right) & {\mbox{if } y \le 0.0.} \end{array} \right. \end{equation} \tag{ 55 }$$

Here, a = 0.01 is the characteristic thickness of the shear layer, and v_sh = 0.5 corresponds to a relative Lorentz factor of 2.29. The initial uniform pressure is p = 1.0. The density is initialized using the shear velocity profile, with ρ = 1.0 in regions with v_sh = 0.5 and ρ = 10⁻² in regions with v_sh = −0.5. The magnetic field components are given in terms of the poloidal and toroidal magnetization parameters μ_p and μ_t as

$$\begin{equation} (B_{x}, B_{y}, B_{z})=\left(\sqrt{2 \mu _{p} p}, 0, \sqrt{2 \mu _{t} p}\right), \end{equation} \tag{ 56 }$$

with μ_p = 0.01 and μ_t = 1.0. The instability is seeded by a single-mode perturbation of the form

$$\begin{equation} v_{y}=\left\lbrace \begin{array}{@{}ll}A_{0} v_{{\rm sh}} \sin (2 \pi x) \exp \left[ -\left({y-0.5 \over \alpha } \right)^{2} \right] & {\mbox{if } y > 0.0} \\\ms - A_{0} v_{{\rm sh}} \sin (2 \pi x) \exp \left[ -\left({y+0.5 \over \alpha } \right)^{2} \right] & {\mbox{if } y \le 0.0.} \end{array} \right. \end{equation} \tag{ 57 }$$

Here, A₀ = 0.1 is the perturbation amplitude and α = 0.1 is the characteristic length scale over which the perturbation amplitude decreases exponentially. The computational domain covers x ∈ [ − 0.5, 0.5], y ∈ [ − 1, 1] with 256 × 512 cells.

Figure 11 shows the perturbation amplitude (Δv_y ≡ (v_{y, max} − v_{y, min})/2) and volume-averaged poloidal magnetic field (B_pol = $\sqrt{\vphantom{A^A}\smash{\hbox{$B_{x}^{2}+B_{y}^{2}$}}}$) as a function of time for conductivities of σ = 0, 10, 10², 10³, and 10⁵ using an ideal EoS with Γ = 4/3 and using the TM EoS. All cases show an initial linear growth phase. Except for σ = 0, the ideal EoS and the TM EoS have almost the same growth rate and the maximum amplitude is reached at t ∼ 2. The maximum amplitude indicates the transition from the linear to the nonlinear phase. The σ = 0 cases exhibit a lower growth rate and later transition to the nonlinear phase than the higher conductivity cases. Thus, differences in the conductivity and EoS do not affect the growth of KHI, except for σ = 0.

Zoom In Zoom Out Reset image size

Poloidal field amplification via stretching due to the main vortex developed by KHI follows the growth of KHI (see Figure 12). In high-conductivity cases, poloidal field amplification is very large, an increase by almost one order of magnitude. Saturation in the poloidal field amplitude occurs later than the transition from the linear to the nonlinear KHI growth phase. This means that magnetic field amplification via stretching continues even after KHI is fully developed. Poloidal field amplification is weaker and saturation occurs earlier when the conductivity is low. Larger poloidal field amplification occurs for the TM EoS case than for the ideal EoS case. Therefore, we find that magnetic field amplification via stretching due to the main vortex developed by KHI is strongly affected by the conductivity and the EoS.

Zoom In Zoom Out Reset image size

Figures 12 and 13 show the time evolution of the density and the poloidal to toroidal field ratio ($B_{{\rm pol}}/B_{{\rm tor}}=\sqrt{\vphantom{A^A}\smash{\hbox{$B^{2}_{x} + B^{2}_{y}$}}}/B_{z}$) for high conductivity, σ = 10⁵, using the ideal EoS with Γ = 4/3 (Figure 12) and using the TM EoS (Figure 13). Both cases show formation of a main vortex by growth of KHI in the linear growth phase. In the ideal EoS case, a secondary vortex appears, although not fully developed. However, development of a secondary vortex is not found in the TM EoS case. Beckwith & Stone (2011) found that a secondary vortex did not appear even in very high resolution simulations using the HLL approximate Reimann solver in ideal RMHD. In our simulations, we also used the HLL approximate Riemann solver to calculate the numerical flux, but do find a secondary vortex in the ideal EoS case. The difference is likely the result of the reconstruction scheme used here and the different reconstruction scheme used by Beckwith & Stone. In the nonlinear phase, the main vortex is distorted and stretched. The magnetic field is strongly amplified by shear motion in the vortex in the linear phase and by stretching in the nonlinear phase. As the mixing layer grows the field lines are bunched into a filamentary-like stretched structure. In the TM EoS case, the vortex becomes strongly elongated along the flow direction. The structure created in the nonlinear phase is very different in the ideal EoS and the TM EoS cases.

Zoom In Zoom Out Reset image size

The field amplification structure for different conductivities from σ = 0 to 10³ is shown in Figure 14. As seen in the time evolution of the averaged poloidal field shown in Figure 11, the magnetic field amplification is weaker when the conductivity is low. Field amplification is a result of fluid motion in the vortex. In the high-conductivity case, the magnetic field follows the fluid motion, like ideal MHD, and is strongly twisted. When the conductivity declines, the magnetic field is no longer strongly coupled to the fluid motion. Therefore, the magnetic field is not strongly twisted. In the case using the TM EoS, we see the same trend for different conductivity and do not show the result here.

Zoom In Zoom Out Reset image size

5.2.3. Relativistic Magnetic Reconnection Test

The final test involves relativistic magnetic reconnection. Pioneering work on relativistic magnetic reconnection using an RRMHD code and 2D simulations was performed by Watanabe & Yokoyama (2006) who considered Petschek-type reconnection. Zenitani et al. (2010) also have studied details of Petschek-type reconnection via 2D RRMHD simulations. Zanotti & Dumbser (2011) have investigated the dependence of Petschek-type relativistic magnetic reconnection by performing 2D and 3D simulations over a broad range of conductivities and magnetizations. Takahashi et al. (2011) have studied Sweet–Parker-type relativistic magnetic reconnection using 2D RRMHD simulations. In this test, we present simulations of Petschek-type reconnection and investigate the effect of the EoS.

We use initial conditions similar to that used in previous work (Watanabe & Yokoyama 2006; Zenitani et al. 2010; Zanotti & Dumbser 2011). The density and gas pressure are given by

$$\begin{eqnarray} \rho &=& \rho _{b} + \mu _{m} \cosh ^{-2} (y) \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} p &=& p_{b} + \mu _{m} \cosh ^{-2} (y), \end{eqnarray} \tag{ 59 }$$

where ρ_b = p_b = 0.1 are the uniform density and gas pressure outside the current sheet, and μ_m = B²₀/(2γ²₀) = 1.0 is the magnetization parameter. The velocity field is initially zero, hence γ₀ = 1. The magnetic field changes orientation across the current sheet according to

$$\begin{equation} B_{x} = B_{0} \tanh (y), \end{equation} \tag{ 60 }$$

where B₀ is calculated from the magnetization parameter. The current distribution is given by

$$\begin{equation} j_{z} = B_{0} \cosh ^{-2} (y). \end{equation} \tag{ 61 }$$

Over the whole computational domain there is a small background uniform resistivity η_b = 1/σ_b = 10⁻³, except within a circle of radius r_η = 0.8 which defines a region of anomalous resistivity with amplitude η₀ = 1.0. The resistivity can be written as

$$\begin{equation} \eta = \left\lbrace \begin{array}{@{}ll} \eta _{b} + \eta _{0}[2 (r/r_{\eta })^{3} - 3(r/r_{\eta })^{2} + 1] & {\mbox{for } r \le r_{\eta },} \\ \ms \eta _{b} & {\mbox{for } r > r_{\eta },} \end{array} \right. \end{equation} \tag{ 62 }$$

where $r = \sqrt{\vphantom{A^A}\smash{\hbox{$x^{2} + y^{2}$}}}$. The electric field is calculated from the resistivity distribution as

$$\begin{equation} E_{z} = \eta j_{z}. \end{equation} \tag{ 63 }$$

The computational domain is x ∈ [ − 50, 50], y ∈ [ − 20, 20] with 2000 × 800 cells, and outflow boundary conditions are used in both directions.

Figure 15 shows the density, the x-component of the four-velocity, γv_x, and the out-of-plane current j_z at t = 100 using an ideal EoS with Γ = 4/3 and the TM EoS. In this figure, we confirm the essential features of Petschek-type relativistic reconnection reported in previous work (Watanabe & Yokoyama 2006; Zenitani et al. 2010; Zanotti & Dumbser 2011). In both the ideal and TM EoS cases, we see similar time evolution and morphology. After an initial adjustment stage of t ⩽ 10, the reconnection process starts around the point at (x, y) = (0, 0) triggered by anomalous resistivity. The magnetic field shows a typical X-type topology. As a result of reconnection, the magnetic energy is converted into both thermal and kinetic energy. Two magnetic islands (the so-called plasmoids) which correspond to the high-density region in Figure 15 move in opposite directions (the figure shows only half of the simulation region and only one magnetic island) and are accelerated along the direction of the magnetic field. A fast reconnection jet is formed inside a narrow nozzle within a pair of slow shocks (Petschek slow shock). The reconnection jet collides with a plasmoid in front of the current sheet further downstream. The plasmoid is surrounded by strong currents (see Figures 15(e) and (f)), which also correspond to the slow shocks (Ugai 1995). These slow shocks surround the plasmoid and are connected to the Petschek slow shocks. In the TM EoS cases, the plasmoid has a faster speed than in the ideal EoS case with Γ = 4/3, and the plasmoid in the TM EoS case propagates further.

Zoom In Zoom Out Reset image size

The time evolution of the maximum outflow velocity (v_{x, max}), the volume-averaged magnetic energy (B²), and the normalized reconnection rate is shown in Figure 16. The outflow gradually accelerates and nearly saturates by t ∼ 60 with v_x ∼ 0.8c. The figure clearly shows that the outflow speed in the TM EoS case is slightly faster than in the ideal EoS case. The magnetic energy gradually decreases with time. Dissipated magnetic energy results in an increase to both the thermal and kinetic energy. Increase in the kinetic energy accompanies acceleration of the outflow (plasmoid). The normalized reconnection rate is defined as $\mathcal {R}=E_{z}^{*}/v_{A, {\rm in}} B_{{\rm in}} \sim v_{{\rm in}}/v_{{\rm out}}$, where E*_z is the electric field at the reconnection point and the upstream properties with subscript ``in" are evaluated at (x, y) = (0, 3) (Zenitani et al. 2010).¹ The reconnection rate saturates at about t = 50 in both cases. However, the TM EoS case has a larger reconnection rate than the ideal EoS case ($\mathcal {R} \sim 0.16$ in the ideal EoS case with Γ = 4/3 and $\mathcal {R} \sim 0.17$ in the TM EoS case). Therefore, the different EoSs lead to a quantitative difference in relativistic magnetic reconnection.

Zoom In Zoom Out Reset image size

The role of the EoS in RRMHD using a newly developed RRMHD code has been investigated. A number of numerical tests in 1D and multi-dimensions have been performed to check the robustness and accuracy of the new code. All of the tests show the effectiveness of the new code in situations involving both small and large uniform conductivities.

The 1D tests of the propagation of a large-amplitude circularly polarized Alfvén wave show that the new RRMHD code reproduces ideal RMHD solutions when the conductivity σ is high and that the code has second-order accuracy. The 1D self-similar current sheet tests indicate that the analytical solution in the moderate conductivity regime is well described by the new code. In a simple MHD version of the Brio and Wu 1D shock-tube test, the code is stable using Strang's splitting technique even when the conductivity is high (σ ∼ 10⁶). In the 2D cylindrical explosion tests, at the high conductivity σ ⩾ 10³, the results recover the solution from ideal RMHD. The results of the KHI test show that the growth rate of KHI is independent of the conductivity, except for very low conductivity (σ ≃ 0). However, magnetic field amplification via stretching of the main vortex developed by KHI strongly depends on the conductivity. The effect of conductivity on magnetic field amplification via KHI in 3D is a topic for future study.

EoSs proposed by Mignone et al. (2005) and Ryu et al. (2006), which closely approximate Synge's single-component perfect relativistic gas EoS, have been incorporated in the code. In the limit of non-relativistic and ultrarelativistic temperatures, the equivalent specific heat ratio associated with the EoSs that approximate Synge's EoS appropriately changes from the 5/3 to the 4/3 limits.

The numerical tests studied the effect of the EoS on shocks, blast waves, the KHI, and relativistic magnetic reconnection. The results provide a useful guide for future more specific studies of each topic. The tests confirm the general result that large temperature gradients cannot be properly described by an ideal EoS with a constant specific heat ratio. The results using a more realistic EoS, which we have studied here, show considerable dynamical differences. The 1D shock-tube tests (Balsara Test2) show that the results obtained from the TM EoS and RC EoS cases are considerably different from the constant γ-law EoS case with γ = 5/3. In the 2D KHI tests, the nonlinear behavior depended on the EoS, even though the growth rate of the KHI was almost the same. In reconnection tests, the approximate EoS cases resulted in a faster reconnection outflow speed and a larger reconnection rate than the ideal EoS case. We conclude that any studies of shocks, instabilities, and relativistic magnetic reconnection should use a realistic approximation to Synge's EoS.

Y.M. thanks H. Takahashi, S. Zenitani, K.-I. Nishikawa, and P. E. Hardee for useful comments and discussion. This work has been supported by NSF awards AST-0908010 and AST-0908040, NASA award NNX09AD16G, and the Taiwan National Science Council under the grant NSC 100-2112-M-007-022-MY3. Y.M. also acknowledges partial support from the NAOJ Visiting Scholar Program (Short-term). The simulations were performed on the Columbia Supercomputer at the NAS Division of the NASA Ames Research Center, the SR16000 at YITP in Kyoto University, and Nautilus and Kraken at the National Institute for Computational Sciences in the XSEDE project supported by National Science Foundation.