New Numerical Scheme for Resistive RMHD

A New numerical scheme for resistive relativistic MHD using method of characteristics

Makoto Takamoto Theoretical Astrophysics Group, Department of Physics, Kyoto University Tsuyoshi Inoue Division of Theoretical Astronomy, National Astronomical Observatory of Japan,
Abstract

We present a new numerical method of special relativistic resistive magnetohydrodynamics with scalar resistivity that can treat a range of phenomena, from nonrelativistic to relativistic (shock, contact discontinuity, and Alfvén wave). The present scheme calculates the numerical flux of fluid by using an approximate Riemann solver, and electromagnetic field by using the method of characteristics. Since this scheme uses appropriate characteristic velocities, it is capable of accurately solving problems that cannot be approximated as ideal magnetohydrodynamics and whose characteristic velocity is much lower than light velocity. The numerical results show that our scheme can solve the above problems as well as nearly ideal MHD problems. Our new scheme is particularly well suited to systems with initially weak magnetic field, and mixed phenomena of relativistic and non-relativistic velocity; for example, MRI in accretion disk, and super Alfvénic turbulence.

Subject headings:
plasma, relativistic resistive MHD, methods: numerical

1. Introduction

The magnetohydrodynamics (MHD) approximation has some interesting properties, for example, the flux freezing and magnetic pressure; the former can be used for the collimation of the jet, and the latter for the acceleration of the plasma. Thus, the magnetic field is considered an essential ingredient for many astrophysical phenomena. In particular, many observations indicate that most of the high energy phenomena in astrophysics are related to the strongly magnetized relativistic plasma around some compact objects, for example, AGN (Antonucci, 1993; Urry & Padovani, 1995), relativistic jet (Blandford & Konigl, 1979; Mirabel & Rodríguez, 1999), pulsar wind (Rees & Gunn, 1974; Camus et al., 2009), gamma-ray bursts (Woosley, 1993; Piran, 2004), and so on. Since it is extremely difficult to solve the relativistic MHD (RMHD) equations analytically, the theoretical investigations in fully nonlinear regimes are mainly based on the numerical simulations (McKinney & Gammie, 2004; Inoue et al., 2010). Most of these studies approximates the plasma as the ideal RMHD fluid. One reason for this is that the ideal RMHD is an excellent approximation of high energy phenomena for ordinary parameters. However, when one considers extreme phenomena, such as the neutron star mergers, or the central engines of GRB, the electrical conductivity can be small, and highly resistive regions may appear. In addition, when one considers the magnetic reconnection, the resistivity plays an essential role in this phenomenon. Magnetic reconnection is one of the most important phenomena, since it is highly dynamic, and it changes magnetic field energy into fluid energy  (Zweibel & Yamada, 2009; Zenitani et al., 2009; Li et al., 2007). Though numerical results of ideal RMHD exhibit magnetic reconnection, this originates in the purely numerical resistivity, and this is unphysical. For this reason, using resistive RMHD is important for the understanding of reconnection and related phenomena.

In order to consider Ohmic dissipation, one only has to take into account an additional term in the induction equation of non-relativistic MHD. However, similar to other non-relativistic dissipation, this induction equation is parabolic and it is well-known that this equation is acausal. As a result, if one takes into account Ohmic dissipation in a relativistic MHD in a similar way, the equation inevitably includes unphysical exponential growing modes, and unstable for small perturbations similar to other dissipation (Hiscock & Lindblom, 1983, 1985). This unphysical divergence results from the fact that one neglects the time derivative of the electric field in the induction equation with Ohmic dissipation. For this reason, when one takes into account the Ohmic dissipation, one has to consider the time evolution of the electric field, that is, one has to deal with the relativistic electromagnetic hydrodynamic equation. This equation is a telegrapher equation, and satisfies the causality.

In this paper, we present a new numerical scheme for the resistive RMHD. There are several examples of pioneering work for resistive RMHD, for example, Komissarov (2007, hereafter K07) proposed numerical method that solves hyperbolic fluxes by using the Harten-Lax-van Leer (HLL) prescription, and damping of the electric field by Ohmic dissipation that is very stiff by using Strang-splitting techniques; Palenzuela et al. (2009, hereafter P09) proposed a numerical method that solves hyperbolic fluxes by Local Lax-Friedrichs approximate Riemann solver, and the stiff part by using implicit-explicit (IMEX) Runge Kutta methods. However, these methods use light velocity as the characteristic velocity, and their numerical solutions are diffusive when one considers problems whose characteristic velocity is much lower than light velocity. This indicates that their numerical solutions are diffusive in many important high plasma dynamics, and also their solutions become highly diffusive when the characteristic velocity of phenomena is much lower than light velocity. In particular, when one solves the dynamics of the accretion disk around a black hole with a relativistic jet, one has to use relativistic resistive MHD code that can solve both highly relativistic and non-relativistic dynamics with resistivity for the following three reasons: (1) the saturation of the magnetorotational instability (MRI) depends on the resistivity; (2) the dynamics of an accretion disk are not ordinarily relativistic, especially, the dynamics of the MRI is sub-Alfvénic; (3) the dynamics of the jet are highly relativistic. For these reasons, previous schemes are diffusive in such phenomena, and we need more accurate numerical schemes. We are developing a new numerical scheme capable of accurately solving problems whose characteristic velocity is quite different from light velocity. In this scheme, we obtain numerical flux of fluid by using sound velocity as the characteristic velocity, and numerical flux of electromagnetic field by using appropriate characteristic velocities of RMHD. This enables us to obtain accurate numerical results when we consider problems whose characteristic velocity is much lower than light velocity. In addition, P09 pointed out that the Strang-splitting method used in the Komissarov method is unstable when applied to discontinuous flows with large conductivities. However, we find that this problem is not related to the Strang-splitting method, but the evolution of electric field during the primitive recovery, that is introduced in the method by P09. By considering this procedure, we can apply the Strang-splitting method to discontinuous flows with large conductivities.

This paper is organized as follows. In Section 2, the basic equations of resistive RMHD are presented. In Section 3, we present the numerical method. Results of numerical test problems previously presented are shown in Section 4. In Section 5, we present results of numerical test problems that cannot be solved accurately by previous codes.

2. Basic Equations

Throughout this paper, we use the units

(1)

In Cartesian coordinates, the Minkowski metric tensor is given by

(2)

Variables indicated by Greek letters take values from to , and those indicated by Roman letters take values from to .

2.1. The Maxwell equations

The covariant Maxwell equations can be written as

(3)
(4)

where is the Maxwell tensor, the Faraday tensor, and the four-vector of electric current.

If we consider highly ionized plasma, the electric and magnetic susceptibilities can be neglected. Then, one has

(5)

where

(6)

is the Levi-Civita alternating tensor of space-time, and is the four-dimensional Levi-Civita symbol.

We introduce a future-directed unit timelike vector normal to a spacelike hypersurface . Using , we can decompose the Maxwell tensor into following forms:

(7)

Similarly, the current four-vector can be decomposed into:

(8)

where is the charge density observed in the rest frame of , and the conduction current satisfying . In the following, we consider only Minkowski spacetime, so .

By using the decomposition of the Maxwell tensor Eq. (7) and the current four-vector (8), the Maxwell equations can be split into the familiar set

(9)
(10)
(11)
(12)

From Maxwell equations, we can derive the electric charge conservation law

(13)

2.2. The hydrodynamic equations

The relativistic hydrodynamic equations can be obtained from the conservation of mass, momentum, and energy

(14)
(15)

where is the mass density current and the energy-momentum tensor defined respectively as

(16)
(17)

where

(18)
(19)

Here is the specific enthalpy, is the proper rest mass density, is the thermodynamic pressure, and is the specific internal energy.

The evolution equation of a relativistic resistive MHD is

(20)

where , , is the density, momentum density, total energy density. In the laboratory frame, , , are given by

(21)
(22)
(23)

where is the fluid three-velocity, is the Lorentz factor, and numerical fluxes are

(24)
(25)
(26)

This is the most common form of perfect fluid equations for the numerical hydrodynamics.

2.3. Ohm’s law

The system of Eqs. (9) - (12), (20) is closed by means of Ohm’s law. Although there are various forms of Ohm’s law, we consider only the simplest kind of relativistic Ohm’s law that accounts only for the plasma resistivity, and that assumes that it is isotropic similar to previous studies K07 and P09. In the covariant form, it is given by

(27)

where is the conductivity, is the resistivity, and is the electric charge density as measured in the fluid frame.

As the Maxwell equations and fluid equations, we can decompose Eq. (27) into 3 + 1 form, and then the space component of Eq. (27) is given by

(28)

In the fluid rest frame, Eq. (28) becomes

(29)

The ideal MHD limit of Ohm’s law can be obtained in the limit of infinite conductivity (). In this limit, Eq. (28) reduces to

(30)

Splitting this equation into the components that are normal and parallel to the velocity vector, it becomes

(31)
(32)

From these equations, we can obtain the usual result

(33)

3. Numerical Method

In this section, we present our new numerical scheme for the resistive RMHD. Since the pioneering studies of resistive RMHD K07 and P09 use light velocity as the characteristic velocity, their solution becomes highly diffusive when characteristic velocity is much lower than light velocity. In our new scheme, we obtain numerical flux of fluid by using sound velocity as the characteristic velocity, and numerical flux of the electromagnetic field by using Alfvén velocity as the characteristic velocity. This enables us to obtain accurate numerical results even when characteristic velocity is much lower than light velocity. In the following sections, we consider the one-dimensional case. The extension to the multi-dimensional scheme using the constrained transport method  (Evans & Hawley, 1988; Stone & Norman, 1992). will be shown in our next paper.

3.1. Strang Splitting method

The relativistic resistive MHD is hyperbolic-relaxation equations. In previous work K07 and P09, they assume that characteristic velocity is the speed of light. Thus, their schemes are highly diffusive when the characteristic velocity is lower than light velocity. For this reason, we apply the Strang splitting method (Strang, 1968) and solve the basic equations by using each appropriate characteristic velocity.

First, we split fluid equations Eq. (20) as follows:

(34)

where

(35)
(36)
(37)
(38)
(39)

The flux of the fluid component can be calculated by using the Riemann solver; the flux of the electromagnetic component can be calculated by the method of characteristics.

Next, we consider the Maxwell equations Eqs. (9) - (12). Eqs. (9) and (10) are not evolution equations but constraint equations, and we treat them separately from evolution equations. The evolution equations of and are Eqs. (11) and (12). By using Ohm’s law Eq. (28), Eq. (11) reduces to

(40)

The source term of this equation includes evolving variables , so this equation is a hyperbolic equation with stiff relaxation terms that requires special care to capture the dynamics in a stable and accurate manner. Thus, we split the charge current into two parts similar to K07

(41)
(42)

Then, we split Eq. (40) into two parts

(43)
(44)

Eq. (43) is non-stiff equations, and Eq. (44) is stiff equations.

As a result, the evolution part of the Maxwell equations can be rewritten as

(45)
(46)
(47)

In component form, Eqs. (45) and (46) reduce to

(48)
(49)
(50)
(51)
(52)
(53)

We solve Eqs. (49), (50), (52), and (53) using method of characteristics (MOC), which will be shown in Sec. 3.2. Eq. (51) is solved using the Runge-Kutta method. The numerical scheme for the stiff equation Eq. (47) will be shown in Sec. 3.3.

3.2. Method of characteristics

The method of characteristics can be used to solve the initial value problems of advective and hyperbolic equations. As is well known, the Maxwell equations are hyperbolic, so we can solve the Maxwell equations accurately by using this method.

The Maxwell equations for the transverse fields are Eqs. (49), (50), (52), and (53). By adding and subtracting these equations, for , , and , we obtain

(54)
(55)

where is the characteristic velocity, and this is equal to the speed of light in ordinal Maxwell equations.

The transverse fields are recovered from by

(56)
(57)

The left-hand side of Eq. (54) is the total derivative for an observer moving at velocity .

Let us consider conservative discretizations of Eqs. (50) and (52):

(58)
(59)

where superscript means the time-step, and subscript means the coordinate of cell center. Using Eqs. (56) and (57), we can obtain the numerical flux of Eqs. (58) and (59). (See Fig. 1.). The same procedure can be done for time advance of , .

Figure 1.— A schematic drawing of Eulerian-like characteristics when one uses piecewise linear interpolation. is the characteristic velocity. On the left is the subsonic case, and on the right is the supersonic case. These figures show that half time-step transverse electromagnetic field and are determined by the fields at the base of two characteristics.

The characteristic velocity of the Maxwell equations in vacuum is light velocity. However, since we consider the electromagnetic hydrodynamics equations, appropriate characteristic velocity has to be used for them. Also, because we consider resistive systems, the characteristic velocity varies with the conductivity and the scale of wave modes. For example, as shown in Appendix. the transverse electromagnetic hydrodynamic waves propagate with the light velocity when is large, where is the wave number, and they propagate with the Alfvén velocity when is smaller than a critical value depending on , and . Because of the finite resistivity, the frequency of the transverse waves has an imaginary part (damping rate), which is a increasing function of . In this scheme, we use the light velocity as the characteristic velocity when is smaller than the critical value; when is larger than the critical value, we use appropriate magnetohydrodynamic characteristic velocities. The critical value is determined so that the transverse waves whose phase velocities are light velocity are dissipated within the numerical integration timestep . A detailed procedure to judge whether we use the light velocity or magnetohydrodynamic characteristic velocities is given in Appendix.

In addition to the numerical flux of the Maxwell equations, the characteristic velocity is also required to construct the Maxwell stress terms and the Poynting flux term. When the characteristic velocity obtained from the analysis of the transverse waves is the light velocity, we use the light velocity as the characteristic velocity for them; when the transverse wave characteristic velocity is the Alfvén velocity, we use the characteristic velocities for them as follows. Note that if the following characteristic velocities are not used, numerical integration becomes unstable, and unstable numerical oscillation occurs.

Then, the necessary procedures are as follows:

  1. For the numerical flux of Maxwell equation Eqs. (58) and (59), we use the Alfvén wave velocity in laboratory frame because the information of transverse electromagnetic fields are transmitted by the Alfvén wave 111 In the case of small limit, the Alfvén velocity in laboratory frame becomes . We find that when this is also small, numerical oscillations occur and numerical integration becomes occasionally unstable. We can prevent this purious oscillations, if we use the characteristic velocity as when . In this case, our scheme for solving the Maxwell equations becomes equivalent to the HLLE scheme. Note that introduction of this modification does not change any results presented in this paper. . In relativistic MHD, the Alfvén velocity in laboratory frame can be obtained by solving (Anile, 1990)

    (60)

    where , , , and is the covariant magnetic field defined as

    (61)
  2. For the Maxwell tension terms in Eq. (38), we use the Alfvén wave velocity in fluid comoving frame because magnetic tension force is originated by the Alfvén wave. In relativistic MHD, the Alfvén velocity in fluid comoving frame is given by

    (62)
  3. For the Poynting flux of energy equation Eq. (39) and the Maxwell pressure terms in Eq. (38), we use the fast magnetosonic wave velocity in laboratory frame because the magnetic pressure originates in the magnetosonic wave. In relativistic MHD, fast magnetosonic wave velocity in the laboratory frame can be obtained by solving

    (63)

    Eq. (63) is a quartic equation, and in an ordinary one has to use the Newton-Raphson method or the quartic formula for obtaining solutions. However, since our scheme splits the fluid part and the electromagnetic part, the sound velocity can be set equal to zero. Then, the characteristic equation Eq. (63) reduces to

    (64)

    By using the quadratic formula, one can obtain solutions of above equation:

    (65)

To sum up, we only have to substitute the appropriate characteristic velocities , and into in Eq. (54), and calculate the electromagnetic field at half time step. Then, the numerical fluxes of electromagnetic hydrodynamics equations are given by

(66)
(67)

where means that they are calculated by using the Alfvén velocity in comoving frame, and by using the fast magnetosonic wave velocity in laboratory frame. For the numerical flux of the Maxwell equation, one has to use the Alfvén velocity in laboratory frame for the calculation.

3.3. Stiff part

As explained Sec. 3.1, Eq. (44) contains stiff terms. Following the previous work K07, we split the equation into components normal and parallel to the velocity vector.

(68)
(69)

Since we use the Strang splitting method, the right-hand side of the above equations can be considered constant other than the electric field . As a result, these equations can be solved analytically

(70)
(71)

where and suffix 0 indicates the initial component. If we use the explicit integrator, the stiff equation has to be solved in very small time steps . However, since Eqs. (70) and (71) are formal solutions, we can avoid the stability constraints of the time step. In the context of ambipolar diffusion in partially ionized plasma, a similar numerical technique using the piecewise formal solution of stiff part is known to be useful scheme  (Inoue et al., 2007; Inoue & Inutsuka, 2008, 2009).

3.4. Constraint Equations

It is well known that Eqs. (9) and (10) are constraints on the Cauchy surface. Though Maxwell equations ensure that these constraints are preserved at all times, straightforward numerical integration of Maxwell equations does not preserve these properties because of the accumulated numerical error. This causes corruption of numerical results, and results in a crash in the end. For this reason, there are a number of numerical techniques for avoiding this problem. We have implemented hyperbolic divergence cleaning for the electric field. The main idea of the hyperbolic divergence cleaning is that one defines new variable as the deviation from constraint equations, and arranges a system of equations to decay or carry the deviation out of the computational domain by high speed waves. For the magnetic field, if one sets constant, the constraint equation can be satisfied in one-dimensional case. In the multi-dimensional case, we can implement constrained transport method (Evans & Hawley, 1988; Stone & Norman, 1992). The detailed implementation will be presented in our next paper.

For hyperbolic divergence cleaning, we modify Eqs. (9) and (51)

(72)
(73)

where is a new dynamic variable and a positive constant. Clearly, when we set , we can recover standard Maxwell equation Eq. (9). From these equations, we can obtain the telegrapher equation for

(74)

Thus, propagates at the speed of light and decays exponentially over a timescale .

Similar to Eq. (44), Eq. (72) contains stiff source terms. Thus, we split the equation into a stiff part and non-stiff part

(75)
(76)

The analytical solution of Eq. (76) is

(77)

where is the initial value of .

3.5. Primitive recovery

In order to compute numerical flux (35), (36), (37), (38), and (39), the primitive variables have to be recovered from the conserved variables . In conserved variables, and can be obtained by evolving the Maxwell equations. However, as pointed out by P09, it is more stable to perform evolution of stiff part Eqs. (70) and (71) during this primitive recovery process when is large, i.e., ideal MHD approximation is valid. This is because when we consider MHD approximation, the electric field is equal to ; however, in general, primitive recovered does not satisfy this relation. In what follows we explain the primitive recovery procedure following P09.

  1. Set an initial guess for the velocity by using previous time step value Then, evolve electric field using Eqs. (70) and (71).

  2. Subtract Poynting flux and electromagnetic energy density from conserved variables, and new variables can be defined as follows:

    (78)
    (79)

    Then, variables are the ideal relativistic fluid conserved variables, and can be recovered by using the ordinary procedures.

  3. Replace the initial guess for the velocity with the obtained velocity , and repeat the steps 1 - 3 until the primitive variables converge.

3.6. Algorithm

In this section, we provide the detailed numerical algorithm.

In Cartesian coordinates, the relativistic resistive MHD equations written in conservative fashion are

(80)

where

(81)
(82)
(83)
(84)
(85)
(86)
(87)
(88)

The electric field and magnetic field are evolved by the Maxwell equations. If the Ohmic dissipation is considered, the Maxwell equations have stiff and non-stiff part. The non-stiff part is

(89)
(90)
(91)

where . The Maxwell equations are consistent with the equation of charge conservation. However, numerical errors in general destroy the conservation law in a way similar to the constraint equations. Thus, the above equation contains the equation of charge conservation.

As explained in Sec. 3.3 and 3.4, the stiff-part is evolved by using the formal solution

(92)
(93)
(94)

Using the above system equations, the second-order numerical algorithm is given as follows.

  1. Advance the Stiff-part equations over by using the formal solutions Eqs. (92) - (94).

  2. Advance the non-stiff part of Maxwell equations Eqs. (89) and (91) over by using method of characteristics as explained in Sec. 3.2, and calculate numerical flux (38) and (39). On the other hand, numerical flux (35) - (37) can be calculated by using approximate Riemann solver  (Martí & Müller, 1994; Martí, 1996; Banyuls et al., 1997; Aloy et al., 1999; Pons et al., 2000; Font et al., 2000; Del Zanna & Bucciantini, 2002; Martí Müller, 2003; Mignone & Bodo, 2005; Mignone et al., 2005). In this paper, we use the HLLC solver.

  3. Advance conserved variables D, , e over the half time-step by using Eqs. (34) - (39). Then, calculate primitive variables of half time step by primitive recovery explained in Sec. 3.5. In our scheme, electric field has to be evolved by using formal solution (92) and (93) during primitive recovery. Primitive variables obtained through this procedure are used for the calculation of the numerical flux at .

  4. Again, advance initial stiff variables over by using the formal solutions Eqs. (92) - (94).

  5. Calculate temporal second-order numerical flux (35) - (39) by using primitive variables obtained through the procedure 3. Then, advance conserved variables over by Eq. (34), and electric field and magnetic field by the Maxwell equations of stiff-part Eq. (91).

  6. Calculate primitive variables by a primitive recovery process. During this process, the electric field is advanced over by using formal solutions (92) and (93).

For the spatial second-order, we use the MUSCL scheme by Van Leer explained in Appendix.

Note that if we evolve electric field in integration of stiff equations or primitive recovery procedure, we have to evolve other primitive variables. This is because conserved variables are not changed during those procedures, and this means that the change of electric field affects all other primitive variables.

4. Test simulations

In this section, several one-dimensional test simulations given in previous studies K07 and P09 are presented. For the numerical flux of fluid, we use the HLLC solver (Mignone & Bodo, 2005). We use an ideal equation of state with , and Courant number, .

4.1. Large amplitude CP Alfvén waves

This test consists of the propagation of a large amplitude circularly-polarized Alfvén waves along a uniform background field . The analytical exact solution of this problem is given by Del Zanna et al. (2007) (Del Zanna et al., 2007); and this problem is used as the ideal-MHD limit test problem by P09. We use the same condition as P09.

(95)
(96)

where , , is the wave number, and is the amplitude of the wave. The special relativistic Alfvén speed is given by

(97)

For the initial data parameters, we have used , and . Using these parameters, the Alfvén velocity is . For the boundary condition, the periodic one is used. In addition, we use a high uniform conductivity following P09, since this is the exact solution of ideal relativistic MHD.

Fig. 2 is results of our new code at (one Alfvén crossing time) for three different resolution cases with . The computational domain is . This result indicates that our new code reproduces ideal relativistic MHD solutions when the conductivity is high.

Figure 2.— The results of large amplitude circularly-polarized Alfvén wave test with large conductivity . This test is carried out for three different grid points .

In these test problems, we cannot achieve full second-order accuracy. The left-hand side of Fig. 3 is the norm errors of the tangential magnetic field of this test problem. This figure shows that our numerical result is nearly -order convergence. We estimate this is because our scheme uses many operator splittings, and the time accuracy of our scheme worsens. Note that this problem is one of the most difficult to solve in relativistic resistive MHD, since this is the limit of large conductivity . The right-hand side of Fig. 3 is the norm errors of the of the next test problem. Since the conductivity is moderate value in that test problem, our new scheme achieves second-order convergence.

Figure 3.— norm errors of the tangential magnetic field under different grid resolution for the second-order schemes using the new scheme. The left-hand side is the result of Large amplitude CP Alfvén waves, and the right-hand side is the result of the self-similar current sheet.

4.2. Self-similar current sheet

This problem is used as the test problem of highly resistive cases in K07 and P09. In this test, it is assumed that the magnetic pressure is much smaller than gas pressure, so that the background fluid is not influenced by the evolution of the magnetic field. We assume the magnetic field has only tangential component , and changes its sign within this current sheet. Since we are interested only in the evolution of the magnetic field, the background fluid is set initially in equilibrium, . In addition, we assume that the conductivity is high, and the diffusion timescale is much longer than the light propagating timescale. Although the resistive relativistic MHD equation is hyperbolic, this assumption allows us to neglect the displacement currents at least in the rest frame. As the result, the evolution equation is reduced to

(98)

This equation has exact solution

(99)
(100)

where erf is the error function. Following K07 and P09, we set the initial condition at with , , , and . The computational domain is , and the number of grid points is . Fig. 4 is the numerical result at . This figure shows that our scheme can solve a highly resistive problem accurately. The convergence rate is consistent with the second-order spatial and temporal discretization.

Figure 4.— The result of self-similar current sheet test comparing the exact solution. The solid line is the exact solution, and the dotted line is the numerical result with conductivity .

4.3. The propagation of Alfvénic transverse waves with Ohmic dissipation

In order to confirm the capability of our method for the relativistic resistive MHD, we perform the test calculation of the propagation of Alfvénic transverse waves with Ohmic dissipation, and compare the results with the exact dispersion relation Eq. (A10) obtained in Appendix.

As explained in Appendix, the resistive relativistic magnetohydrodynamic equation contains transverse wave modes that become the light wave in large region, and become the Alfvén wave in small region. To demonstrate the propagation of transverse waves, we set the initial condition by eigenfunctions of the mode obtained from Eqs. (A6) - (A9)

(101)
(102)
(103)
(104)

where and , and is the solution of the dispersion relation Eq. (A10). We set the same parameters in Appendix.

(105)

Since the enthalpy includes the information of the equation of state, one can take any value of . In this calculation, we set and . The computational domain covers the region where the periodic boundary condition is imposed, and the number of grid points is .

The propagation speed of the numerical waves can be determined by tracing the position where is maximum. We measure the propagation speed and evaluate Re based on the time when the maximum of reaches again, i.e. one-wave crossing period. The damping rate Im is measured by using where is the maximum of after the one-wave crossing time.

In Figs. 5, we plot the real and imaginary part of against . The solid line is the exact dispersion relation obtained in Appendix. We have performed the calculation in the cases of . These figures show that our new numerical code can reproduce the propagation of Alfvénic transverse waves accurately for any value of the conductivity .

Figure 5.— The result of the propagation of Alfvénic transverse waves with Ohmic dissipation test problem. The solid line is the exact dispersion relation, and the dots are the numerical solutions. The test calculations are performed in the cases of .

4.4. Shock-tube problem

For the shock tube test problem, we consider the simple MHD version of the Brio and Wu test as P09. The initial left and right states are given by

(106)
(107)

All the other fields are set to .

Fig. 6 is the numerical results at that change grid points . The computational domain covers the region . We also plot an ideal RMHD solution by the solid line computed by a publicly available code developed by Giacomazzo and Rezzolla (Giacomazzo & Rezzolla, 2006). The conductivity is uniform with . The solution of this Riemann problem contains a rarefaction moving to the left, a shock moving to the right, and a tangential discontinuity between them. Fig. 6 shows that our numerical solution of the resistive MHD can reproduce the profile of an ideal MHD shock tube problem using high conductivity . In addition, our numerical solution captures contact discontinuity as sharp as P09.

Fig. 7 is the numerical results of the same problem that changes the conductivity . We also plot the ideal RMHD solution by the solid line. The number of grid points is . This result shows that our numerical solution reproduces nearly the same results as P09.

P09 reports that Strang’s splitting method becomes unstable for moderately high values of the conductivity for this shock tube problem, and one has to use the implicit method. However, this is not related to whether one uses Strang’s splitting or implicit method, but to the revision of the electric field during the iteration of the primitive recovery (H. R. Takahashi 2010, private communication). Our scheme uses Strang’s splitting, but can solve this shock tube problem stably even when , if we revise the electric field during the primitive recovery as explained in Sec. 3.5.

Figure 6.— The numerical results of the Riemann shock tube test problem for three different grid points . We use the conductivity . The solid line is the ideal solution.
Figure 7.— The numerical results of the Riemann shock tube test problem for different conductivity cases: . The number of grid points is . The solid line is the ideal solution.

5. Test Simulations for fluid dominated case

The previous studies K07 and P09 use light velocity for the characteristic velocity. Thus, their numerical solutions become highly diffusive when one considers problems whose sound velocity or Alfvén velocity is much lower than light velocity. In this section, we perform test problems in such cases, and compare the results of the HLL code with that of our code.

5.1. Shock tube test problem

In this section, we compute a high plasma shock tube problem, and compare results of our code with those of the HLL code. The initial left and right states are given by

(108)
(109)

All the other fields are set to .

Figs. 8 are the numerical results of our code and the HLL one, being compared with ideal solutions at . The number of grid points is .

Figure 8.— The numerical result of shock tube problem for fluid energy dominated case comparing with that of HLL and ideal solution. On the left is the density profile, and on the right is the profile of the tangential magnetic field . The number of grid points is .

These figures show that HLL solver becomes more diffusive than our code. In addition, Figs. 8 show that the density profile of the shock heated region somewhat overshoots that of the ideal solution, and tangential magnetic field slightly undershoots that of the ideal solution. These results show that when the plasma is high, the HLL solver becomes highly diffusive and does not reproduce the correct value of the shock heated region. In contrast, our numerical results reproduce ideal solutions very well even for high problems.

5.2. The propagation of contact discontinuity

In this section, we calculate the propagation of a contact discontinuity, and study the accuracy of capturing the contact discontinuity for various advection velocities. When one uses the HLL code by Komissarov, the numerical results can be expected to be diffusive for the case of very slow advection velocity, since the HLL code uses light velocity for the characteristic velocity. In contrast, our new code uses sound velocity for the fluid characteristic velocity, and the numerical results will be more accurate for any advection velocity.

We consider the propagation of contact discontinuity of magnetohydrodynamics. The initial condition is

(110)
(111)

All the other fields are set to .

Since we want to consider the ideal fluid case, we consider high conductivity . We use an equation of state with , and the computational domain covers the region with grid points. The CFL number is , and the integration is carried out until fluid crossing time. For the boundary condition, the periodic one is used. For the advection velocity, we use the following velocities:

(112)
Figure 9.— The numerical results of the propagation of contact discontinuity of RMHD by using our new scheme for different advection velocity: . The number of grid points is .
Figure 10.— The numerical results of the propagation of contact discontinuity of RMHD by using the HLL code for different advection velocity: . The number of grid points is .

On the left of Figs. 9 are the numerical results of the density profile calculated by using our new code, and on the left of Figs. 10 are the numerical results of the density profile calculated by using the HLL code, The solid lines are the ideal solution. This figure shows that the numerical results of density profiles by HLL code of is nearly equal to that of our new code. However, the numerical results by HLL code become more diffusive than by our code as the advection velocity becomes small; in contrast, the accuracy of the numerical results by our code is nearly independent of the advection velocity. The right hand side of Figs. 9 are the numerical results of tangential magnetic field calculated by using our new code, and the right hand side of Figs. 10 are the numerical results of tangential magnetic field calculated by using the HLL code. Similar to the density profile, the numerical results by using the HLL code become more diffusive than by using our code as the advection velocity becomes small; in contrast, the accuracy of the numerical results by our code is nearly independent of the advection velocity.

In conclusion, the HLL code is not capable of accurately solving problems whose advection velocity is smaller than light velocity, since the HLL code uses light velocity for the characteristic velocity. The diffusive result of HLL code is always problematic for any discontinuity when the propagation velocity is much smaller than light velocity. In contrast, since our code uses appropriate characteristic velocities, the numerical dissipation does not depend on the characteristic velocity. For this reason, our new code can solve any advection velocity problems accurately, especially for the problems including discontinuities.

5.3. The propagation of small amplitude Alfvén wave

In this section, we consider the propagation of small amplitude Alfvén waves in high plasma. The integration is performed for different resolutions, and we compare the numerical results of Komissarov’s HLL code and our code. For the application to the numerical simulation of MRI, the integration is performed for a small number of grid points: for one wavelength of the Alfvén wave; this corresponds to the number of grid points for resolving the wavelength of maximum growth rate of MRI.

For the initial condition, we consider

(113)
(114)
(115)

In this case, the Alfvén velocity and plasma beta are given by

(116)

where the Alfvén velocity and the plasma beta is defined as

(117)
(118)

Since the initial magnetic field is very weak for most of the MRI phenomenon, a weak magnetic field is considered. In order to consider the ideal fluid case, we set a high conductivity . We use an equation of state with . The computational domain covers the region . The CFL number is , and the integration is carried out until Alfvén wave crossing time. For the boundary condition, the periodic one is used.

Figure 11.— The numerical result of the propagation of a small slow Alfvén wave. On the left is the result of our new scheme, and on the right is that of the HLL code. The number of grid points is .

The numerical results of our code and HLL are presented in Figs. 11. Although the amplitude of both results falls because of the numerical diffusion, it can be seen that HLL results are more diffusive than our numerical results when the number of grid points is . When the number of grid points is , the numerical result of HLL code is a little more accurate than that of our code. This is because our new scheme uses an operator split for the accuracy, and the convergence rate is a little less than second order in time. However, from a practical point of view, it is impossible to cost 64 grid points for the wavelength of maximum growth rate of MRI in many cases, and still our new code can integrate the growth of magnetic field by MRI more accurately.

In conclusion, when one considers the high plasma, our code is more accurate than the HLL code because our code uses sound velocity and Alfvén velocity as the characteristic velocity. In particular, the above results show that our new method is useful for the application to the phenomena including MRI. This instability occurs in the system whose angular momentum changes as , and the amplitude of the perturbative Alfvén waves grows exponentially in over the duration of nearly one Kepler rotation. Since the above condition is satisfied in most of the differential rotating systems in gravity, MRI is one of the most important astrophysical phenomena. In order to reproduce this instability numerically, one has to resolve the wavelength of maximum growth rate. However, this is difficult for most problems, since this wavelength is proportional to the initial weak magnetic field. For this reason, in order to reproduce MRI numerically, one has to use numerical schemes that can integrate small amplitude Alfvén waves accurately by smaller number of grid points. Then, the results of test problems in this section show that our new numerical scheme can deal such problems more accurately than previous codes.

In these three test problems, we consider extremely high density cases in order to distinguish differences easily. However, this can always happen when the magnetic field is weak. As a result, if one considers problems including an initially weak magnetic field like magnetic rotational instability (MRI) in the accretion disk, our code can produce more accurate results.

6. Conclusion

In this paper, we have presented a new numerical scheme of resistive RMHD for one-dimensional case which can solve matter dominated problems more accurately than the existing numerical method. Since this new scheme uses different characteristic velocity for obtaining the numerical flux of fluid and electromagnetic field, one can solve accurately and stably problems whose characteristic velocity is much lower than that of light.

When one considers relativistic problems, one has to solve stiff equations for electric fields. In general, it is difficult to deal with stiff equations, and special methods have been presented; for example, K07 uses the Strang’s splitting method, and P09 use the implicit method. P09 report that Strang’s splitting method is incapable of solving problems that include discontinuity, such as shock. However, we find that this is not related to the method for the stiff equations, and one can solve problems including shock if one evolves the electric field during the primitive recovery; we use Strang’s splitting method, and the solver is well behaved for shock tube problems. The results of other test problems show that our new scheme is capable of accurately solving both highly resistive problems and nearly ideal MHD ones. In addition, it has been shown that our code can solve low characteristic velocity problems more accurately than the HLL code.

The problems of high density and high plasma appear when one considers magnetorotational instability (MRI) in the accretion disk with a relativistic jet, for example. In this case, one has to use relativistic resistive MHD code that can solve both highly relativistic and non-relativistic dynamics with resistivity for the following three reasons: (1) the saturation of the magnetorotational instability (MRI) depends on the resistivity; (2) the dynamics of an accretion disk are not ordinarily relativistic, especially, the dynamics of the MRI is sub-Alfvénic; (3) the dynamics of the jet are highly relativistic. Our new scheme can solve such problems accurately even when the initial magnetic field is very weak.

We present multi-dimensional extension of our scheme in our next paper.

We would like to thank Bruno Giacomazzo for providing the code computing the exact solution of the Riemann problem in ideal MHD. Numerical computations were in part carried out on Cray XT4 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. This work is supported by Grant-in-aids from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, No. 223369 (T. I.).

Appendix A The dispersion relation of the relativistic electromagnetic fluid

As explained in Sec. 3.2, we solve the evolution of electric and magnetic field by the method of characteristics. For the characteristic velocity, we use the appropriate MHD characteristic velocity when is large, that is, ideal MHD approximation is valid. However, when the conductivity is not so large, we have to replace the characteristic velocity with the speed of light. In this section, we discuss when to switch the characteristic velocity from appropriate MHD characteristic velocity to light velocity. In the following, we calculate the linear perturbation of the relativistic electromagnetic equation in order to obtain characteristic velocity.

The relativistic electromagnetic fluid equations are given by

(A1)
(A2)
(A3)
(A4)
(A5)

To obtain the dispersion relation, we start by expanding physical variables around an unperturbed state in the following frame:

  • The fluid is at rest:

  • The x-coordinate is parallel to the :

  • The magnetic field is in the x-direction:

  • charge neutrality:

Since we only want to judge when to switch characteristic velocity, we consider propagation of the transverse waves along the magnetic field. When one uses this procedure during the numerical simulation, one only has to calculate of the simulation data, and substitute its square root into the above . This is because is a scalar, and becomes the square of the magnetic field in the fluid comoving frame because of the assumption of the charge neutrality. Since the magnetic field appears only in the form of in the following procedure, one can neglect the sign of magnetic field. In the following, we consider only the characteristic velocity of transverse waves.

In the above condition, the Alfvén mode is included in the z-component of the velocity and magnetic field , and decouples from other variables. For this reason, we consider only variables related to and .

We replace the current vector in Eq. (A1) with Eq. (A3). Then, the perturbed equations are

(A6)
(A7)
(A8)
(A9)

From these equations, the following dispersion relation is obtained:

(A10)

Eq. (A10) is the biquadratic equation with respect to , and has the formula of radicals. However, the analytical formula is very complex and hard to analyze, and is not suitable for obtaining the characteristic velocity. As shown below, the transverse wave becomes an Alfvén wave in the long wavelength regime, and the light wave in the short wavelength regime shown in Fig. 12. Fig. 12 shows that the damping rate is a monotonically increasing function of . For this reason, we establish the following criterion for the characteristic velocity: when all light modes damp during one time step we use appropriate MHD characteristic velocity for the method of characteristics; when some light modes do not damp during one time step , we use light velocity for the method of characteristics. We will discuss this method in detail in the following.

First, we substitute into Eq. (A10), and divide the dispersion relation into a real part and imaginary part. Then, the real part is