Minimally implicit Runge-Kutta methods for Resistive Relativistic MHD

The Relativistic Resistive Magnetohydrodynamic (RRMHD) equations are a hyperbolic system of partial differential equations used to describe the dynamics of relativistic magnetized fluids with a finite conductivity. Close to the ideal magnetohydrodynamic regime, the source term proportional to the conductivity becomes potentially stiff and cannot be handled with standard explicit time integration methods. We propose a new class of methods to deal with the stiffness fo the system, which we name Minimally Implicit Runge-Kutta methods. These methods avoid the development of numerical instabilities without increasing the computational costs in comparison with explicit methods, need no iterative extra loop in order to recover the primitive (physical) variables, the analytical inversion of the implicit operator is trivial and the several stages can actually be viewed as stages of explicit Runge-Kutta methods with an effective time-step. We test these methods with two different one-dimensional test beds in varied conductivity regimes, and show that our second-order schemes satisfy the theoretical expectations.


Introduction
The Relativistic Resistive Magnetohydrodynamic (RRMHD) equations are a hyperbolic system of partial differential equations used to describe the fluid dynamics in the presence of magnetic fields with a finite (but potentially large) conductivity, and when the velocities involved are close to speed of light, c, or the fluid specific energies become comparable or larger than c 2 . In Astrophysical contexts, all these ingredients can be found in scenarios such as compact stars or relativistic jets.
The RRMHD system of equations can be written as a set of evolution equations corresponding to the conservation of rest-mass, momentum, energy, magnetic flux and electric charge. In addition to these conserved variables, Komissarov [1] has extended the system by two extra scalar potentials which control the evolution of the solenoidal constraint for the magnetic field and the divergence of the electric field following the ideas of Dedner et al. [2].
The physical or primitive variables in RRMHD are where B i (i = 1, 2, 3) are the components of the magnetic field, E i are the components of the electric field, q is the electric charge, ρ is the rest-mass density, is the specific internal energy and v i are the components of the flow velocity field. The conserved or numerical variables of this system are U = (Φ, B i , Ψ, E i , q, D, τ, P i ), (2) where D = ρW is the mass-density, W = 1/ √ 1 − v 2 is the Lorentz factor, τ = ρhW 2 − p − ρW + (E 2 + B 2 )/2, h is the specific enthalpy, p is the pressure, E 2 = E i E i , B 2 = B i B i , P i = ρhW 2 v i + (E × B) i and Φ and Ψ are two scalar potentials controlling the evolution of the constrains ∇ · B = 0 and ∇ · E = q, respectively. We note that the scalar potentials can be regarded also as part of the set primitive variables.
For practical purposes, the system of RRMHD equations can be split into two sets of equations, ones containing stiff source terms and the rest of them: where Y = (Φ, B i , Ψ, q, D, τ, P i ),S j E and S Y contain the corresponding fluxes and non-stiff source terms of the evolution equations, and σ is the conductivity of the system. In the ideal limit, σ → ∞ and the electric field can be expressed in terms of the rest of the primitive variables, so that the corresponding evolution equation for the electric field is therefore not included in the system of evolved equations. In the case of a finite and high-conductivity, the term proportional to σ becomes potentially a stiff source term.
The numerical evolution of the equations with non-stiff source termsS j E and S Y can be carried out with a TVD explicit Runge-Kutta (RK) method. In the high-conductivity regime, an alternative method has to be used instead to deal with the stiffness of the source term in the evolution equation of E. One possibility is the use of a semi-analytic integration of the stiff source terms (and/or the scalar fields), together with a Strang-splitting scheme [1], but then the order of convergence of the numerical scheme can reduce to first-order due to the fact that the kernel of the relaxation operator is non-trivial and corresponds to a singular matrix in the linear case.
A very interesting option to numerically evolve this hyperbolic system with stiff source terms is the use of implicit-explicit (IMEX) RK methods [3,4,5,6], as in the case of [7]. The problems arising in this case are, on one hand, the degradation to first or second-order for a range of values of σ, and, on the other hand, the unclear update of the variables when the primitive variables cannot be explicitly coined in terms of the conserved ones. The second point steams from the fact that the Lorentz factor and the velocity components appear in the stiff source terms. The updated values of these quantities are needed to update the electric field but, in general, the whole set of conserved variables (including the electric field) are needed to recover the velocity components. In [7] this fact was overcome by the addition of an iterative loop over the electric field, which is computationally expensive, and there is no guarantee of convergence to physically acceptable values.
We present here an alternative class of methods to the previous ones, which we name Minimally Implicit RK (MIRK) methods, to avoid numerical instabilities due to the presence of stiff terms without increasing the computational costs in comparison with explicit methods. The MIRK methods reduce to the optimal TVD explicit RK methods of Shu & Osher [8] for theS j E and S Y operators, and implicitly evolve the stiff source terms as we explain in next section; the proposed strategy needs no iterative loop over the electric field, the analytical inversion of the implicit operator is trivial and the several stages can actually be viewed as stages from explicit RK methods with an effective time-step.

Minimally implicit Runge-Kutta methods
In order to derive the MIRK methods, by construction we impose that they must recover the optimal TVD explicit RK methods for theS j E and S Y operators. We also consider updated and non-updated values of the electric field in the evaluation of the stiff source term in Eq. (3), but only non-updated values of the velocity field and Lorentz factor. The later strategy is key to avoid employing extra iterative loops to recover the primitive variables from the conserved ones (see below). In order to simplify the notation, we defineσ = σ W and The proposed method can be written as an usual explicit, s−stages RK method with an effective time-step in the following way: where | n and | n+1 denote the values of variables at time levels t n and t n+1 , respectively,ã ij denotes the coefficients of the Butcher matrix associated to the explicit RK method, a ijl , b ij are the coefficients associated to the implicit part of the method, and ∆t is the time-step. Note that, differently from standard RK-IMEX methods, there is no final explicit substep. Indeed, Eqs. (9) and (10), simply mean that the updated values of the variables at t n+1 are the ones computed after the previous s−substeps of the MIRK algorithm. Equivalently, in the standard notation of Butcher tables (e.g., [7]), the coefficients ω andω are identically zero. We point out that the time update of the conserved variables is explicit, i.e., all the terms in the right-hand side of the equations are evaluated in time levels of the algorithm that have already been computed. Singularly, the updated values of the electric field depend upon already known values of it and of the primitive variables. This particularity of the MIRK schemes is very useful for the RRMHD system, since it avoids the need of performing a double iterative procedure to recover the primitive variables from the conserved ones. This is a clear advantage over the usage of RK-IMEX schemes to deal with the stiffness of the RRMHD system (e.g., [7]). It is also relevant to point out that, in order to write the implicit scheme in an explicit form, a much simpler matrix than in [7] has been inverted analytically, since the terms E l v l v j in the expression for S j E are treated explicitly in MIRK schemes.

First-order method
The first-order MIRK method corresponds to s = 1 in the previous expressions, with all coefficients zero except forã the absolute value of the determinant of the matrix that updates the values of the conserved variables from one timestep to the next one, |M ∆t |, is also bounded by 1. This necessary (not sufficient) condition restricts the range of the coefficients very efficiently and allows to find values for the coefficient(s) very close to the optimal ones. Under the hypothesis that |M ∆t | is bounded by 1 when the stiff source terms are neglected, it can be checked that the choice c 1 = 0 yields a determinant of M ∆t also bounded by 1 when the stiff source terms are included, for ∆t sufficient small. The final expression of the method for this optimal value of the coefficient can be written as: We note that the factor ∆t 1 + ∆tσ| n , can be regarded as an effective time step, that decreases in magnitude asσ increases. Indeed, in the ideal limit, Eq. (13)

Second-order method
The second-order MIRK method corresponds to s = 2 in the general expression of Sect. 2.
Following the same prescription as in the first-order method, i.e., imposing that the method recovers the optimal TVD second-order explicit RK method when stiff source terms are neglected, all coefficients are zero except for (1 − c 1 ) = 0 = (c 1 /2 − c 2 ). If the values of the conserved variables (E i , Y)| (n) are bounded at the time step n, this method guarantees that the values (E i , Y)| (n+1) will also be bounded, even whenσ 1. In order to save RAM memory and improve on computational efficiency, the second-order MIRK scheme can be rewritten as: Written in this form, the second-order MIRK scheme does not need any extra memory storage as compared with an optimal explicit second-order RK TVD scheme. After linearization of the non-stiff source terms, we focus again on the condition |M ∆t | < 1. We find two possible sets of values for the coefficients. If we choose c 2 = (1−c 1 ) 2 2c 1 , we have to restrict to c 1 < 0 or 0.67 < c 1 < 0.75. We find by numerical experimentation (see Sect that only the region c 1 < 0 gives stable results, with the optimal value for c 1 located around −0.1. If we consider the case c 1 = 0, we have to restrict to −1 < c 2 < −0.5. After numerical experimentation we find that the optimal value for c 2 is located around −0.97. Values of both coefficients below zero are favored, and all the options have to be verified numerically.

Numerical results
The MIRK methods described in the previous sections have been implemented in the MRGENESIS code [9,10,11,12]. In order to compute the numerical fluxes we employ a Local Lax-Friedrichs approximate solver [1], in which the states left and right to each interface are obtained by a second-order monotone central (MC) reconstruction of primitive variables. We have also tested the new methods in combination with a more diffusive minmod reconstruction, but the results are affected by a larger diffusion close to discontinuities or big gradients (as expected).
In order to test the newly proposed MIRK schemes, we have considered two different onedimensional tests, each of which proves a different regime of conductivities and, hence, different degrees of stiffness of the RRMHD system. In all the tests considered in the next sections, we assume that the equation of state is that of an ideal gas with a fixed adiabatic index γ.

Self-similar current sheet
Komissarov [1] proposed a test to prove the resistive regime in which σ 100, i.e., relatively far away from the ideal limit. This is the so-called self-similar current sheet, where we have a magnetic field B = (0, B(x, t), 0). The magnetic pressure is much smaller than the gas pressure and B(x, 0) changes sign within a thin current layer. Except for the magnetic field, the rest of the primitive variables are set to be constant in the whole domain: p = 50, ρ = 1, γ = 4/3, E = 0, v = 0 and σ = 100. This test has an analytic (self-similar) solution for t > 0 of the form where B 0 = 1, x 0 = 1.5, the initial time is set to t = 1, and the error function is defined in the usual way: We employ 200 computational cells to cover the domain [0, 3]. We let the system to evolve until t = 10 and compare our results with the analytic solution in Fig. 1. In this case we have used the second order MIRK scheme with c 1 = −0.1 and c 2 = −0.97383794 ( Fig. 1; left) and also with c 1 = 0 and c 2 = −0.85 ( Fig. 1; right). As can be observed, there is basically no difference between the numerical and the analytical solutions. The results do not qualitatively change if the MIRK coefficients are slightly changed with respect to their canonical values. We note that the CFL employed in this test is 0.7.

Circularly polarized Alfvén waves
A second one-dimensional test that proves the ideal limit assuming a very large value for the conductivity is considered by Palenzuela et al. [7]. This test was originally proposed by Del Zanna et al. [13] and has an exact solution in ideal Relativistic Magnetohydrodynamics. The solution describes the propagation of large amplitude Alfvén waves along a uniform background magnetic field with strength B 0 , in a domain with periodic boundary conditions. In the ideal limit, the velocity along the x−direction is zero (v x = 0), and the exact solution reads: where B x = B 0 , k is the wave vector, η A is the wave amplitude, and v A is the relativistic Alfvén speed, that in this case is Though there is no exact solution in the RRMHD regime for circularly polarized Alfvén waves, in the limit of very large conductivity, Eqs. (20) and (21) are approximately satisfied. Thus we compute the evolution from the initial state and until the train of waves completes one period.
In the ideal limit, the solution should be identical to that in the initial state and we expect that if σ 1, also the resistive solution shall be quite similar to the initial one. For this test we fix γ = 2, k = 2π, B 2 0 = 4/3, p = ρ = η A = 1, which leads to a period T = 2, and an Alfvén speed v A = 1/2. Figure 2 shows the comparison of the analytic versus the numerical solutions with different resolutions. It is evident that our numerical solution lies on top of the analytic one, even for very moderate resolutions. As in the previous test, we have computed the solution with different MIRK coefficients. We find that, in this case, much more demanding from the numerical point of view (because of the large value of the conductivity), the fact that we do have only approximately optimal MIRK coefficients (see Sect. 2) has a clear impact on the stability of the method. While the set of coefficients (c 1 , c 2 ) = −(0.1, −0.97383794) yields unstable results for CFL values larger than approximately 0.4 ( Fig. 2; left), the set of parameters (c 1 , c 2 ) = (0, −0.85) is still stable for CFL=0.45. In any case, the range of allowed values of the CFL yielding stable results depends on the magnitude of σ. From an extensive numerical experimentation we conclude that for σ 5 × 10 6 , CFL values up to about 0.7 yield stable results for (c 1 , c 2 ) = (0, −0.85) and small variations of these parameters. For σ 5 × 10 6 progressively smaller values of the CFL coefficient must be used to maintain the stability of the numerical solution.

Order of convergence
We have also performed a number of convergence tests employing as reference the one dimensional problems at hand, varying the value of the second-order MIRK coefficients, c 1 and c 2 . In Fig. 3 we display the L 1 -norm errors for the B y component of the magnetic field as a function of the resolution. The precise definition of the error for a given test including N x uniform cells reads: where L is the domain size, and B y (ref),i is a numerical solution computed at sufficiently high resolution to be considered as a reference solution. In our case, B y (ref),i is computed with N x = 3200.
Our second-order MIRK scheme converges to a bit more than second-order for the self-similar current sheet test for a relatively large range of coefficients c 1 and c 2 (Fig. 3). The order of convergence is slightly smaller than 2 (p = 1.83) for c 1 = −0.1 and c 2 = (1 − c 1 ) 2 /(2c 1 ) = −6.05 and a bit larger than 2 (p = 2.19) for c 1 = 0 and c 2 = −0.9. Thus, in the limit of moderate conductivity the method performs as expected.
In Fig. 4 we display the L 1 -norm errors for the B y component of the magnetic field as a function of the resolution for the test of circularly polarized Alfvén waves. We observe that at low resolution and for very high conductivity values, σ = 10 6 , the second-order schemes have an order of convergence closer to first-order than to second-order. However, as the resolution increases, we find that the order of the method evolves towards the nominal second-order of the algorithm. We also find by numerical experimentation that the best values of the second-order MIRK coefficients lie in the range −0.1 ≤ c 1 ≤ 0, and −1 < c 2 < −0.85. We believe that the difference between the theoretical and numerical order of convergence are triggered by the very large conductivity, which not only increases the stiffness of the system of equations, but also introduces severe looses of accuracy as a result of the very different orders of magnitude among different variables in arithmetic operations (e.g., to compute the electric currents).

Conclusions
We propose here a new kind of RK methods, named minimally implicit RK methods, to numerically evolve some systems of hyperbolic equations with stiff source terms. As a first target for the new MIRK methods, we have considered the numerical evolution of the RRMHD equations as a motivation of the development of the methods. However, we foresee applications to other systems of hyperbolic equations with stiff source terms. At the moment, we are working on applying MIRK methods to coupled system of General Relativistic Radiation Hydrodynamics.
In this system, the neutrino transport equations become very stiff in the opaque regime (when neutrinos are trapped in the matter; see, e.g., [14,15,16]). MIRK methods recover the optimal TVD explicit RK algorithm when the stiff source terms are neglected. They have a very simple analytical inversion of the implicit operators included to deal with the stiff source terms, and their structure can be viewed as an explicit RK method with an effective time-step.
Currently, one dimensional numerical test problems have been considered, but we are working on their application to multidimensional scenarios. Our aim is to improve the current treatment of the resistivity employed in ideal MHD simulations (where the source of such resistivity is not physical but numerical). In this line, we aim to improve the treatment of resistivity in some local numerical simulations of the Kevin-Helmholtz instability in the context of neutron star mergers [17], as well as in global models (e.g., [18]). Likewise, we aim to incorporate physical resistive dissipation to the problem of internal shocks in blazars [19,20,21,22] and in the shocks associated to afterglows [23,24,25].
MIRK methods present some advantages in comparison with the well-known RK IMEX methods: CPU time saving, RAM memory saving and accuracy; these properties have to be explored in more detail. We also seek to compare our methods with Strang-splitting ones.