Relativistic resistive magneto-hydrodynamics code for high-energy heavy-ion collisions

We construct a relativistic resistive magneto-hydrodynamic (RRMHD) numerical simulation code for high-energy heavy-ion collisions as a first designed code in the Milne coordinates. We split the system of differential equations into two parts, a non-stiff and a stiff part. For the non-stiff part, we evaluate the numerical flux using HLL approximated Riemann solver and execute the time integration by the second-order of Runge–Kutta algorithm. For the stiff part, which appears in Ampere’s law, we integrate the equations using semi-analytic solutions of the electric field. We employ the generalized Lagrange multiplier method to ensure the divergence-free constraint for the magnetic field and Gauss’s law. We confirm that our code reproduces well the results of standard RRMHD tests in the Cartesian coordinates. In the Milne coordinates, the code with high conductivity is validated against relativistic ideal MHD tests. We also verify the semi-analytic solutions of the accelerating longitudinal expansion of relativistic resistive magneto-hydrodynamics in high-energy heavy-ion collisions in comparison with our numerical result. Our numerical code reproduces these solutions.


Introduction
Relativistic hydrodynamics has been widely used for description of collective dynamics in various phenomena from nuclear physics to astrophysics. The high-energy heavy-ion collision experiment is one of the active fields of application of relativistic hydrodynamics. a e-mail: knakamura@hken.phys.nagoya-u.ac.jp (corresponding author) At Relativistic Heavy Ion Collider (RHIC), production of the strongly interacting quark-gluon plasma (QGP) was succeeded, which was achieved by measurement of key observables and theoretical interpretation to them [1][2][3][4]. In particular, at that time, the strong elliptic flow was successfully explained by hydrodynamic models. Along with other phenomenological analyses, this fact reached the conclusion that the QGP created at RHIC is not weakly interacting gas but strongly interacting plasma. The analysis of high-energy heavy-ion collisions based on the relativistic viscous hydrodynamics shed light on not only understanding dynamics of space-time evolution after collisions, but also QGP bulk properties such as temperature dependence of shear and bulk viscosities and a diffusion constant. Usually, the QGP bulk properties are discussed in comparison between experimental data and hydrodynamic model calculation. Intensive computation of Bayesian analysis is also performed for extracting detailed information of thermodynamic properties of QCD [5,6].
During the relativistic collision of positively charged heavy nuclei, the highest intense electromagnetic fields in our universe are produced, e.g. |eB| ∼ 10 15 T [7]. The highenergy heavy-ion collision is possible to address the property of interaction between deconfined strongly interacting matter and electromagnetic fields. The effect of electrical conductivity of the medium on the evolution of electromagnetic fields has been discussed with the semi-analytic solutions of electromagnetic fields in Refs. [8,9]. The contribution of electromagnetic fields to the charge-dependent anisotropic flow has been estimated from the semi-analytic solution of simplified relativistic hydrodynamic equations with electromagnetic fields in Refs. [10][11][12]. Furthermore, relativistic hydrodynamic simulation with background electromagnetic fields has been investigated [13,14].
The fully realistic analysis of high-energy heavy-ion collisions would require one to solve Maxwell equations together with relativistic hydrodynamic equations. It is called a full relativistic magneto-hydrodynamic (RMHD) framework, which describes the dynamics of the plasma coupled with electromagnetic fields. Such an analysis with infinite electrical conductivity, relativistic ideal MHD, has been carried out [15,16]. In the relativistic ideal MHD, the electric field measured in the fluid comoving frame vanishes. In order to consider the charge distribution of hadrons in high-energy heavy-ion collisions, the simulation based on the RMHD framework with finite electrical conductivity is needed. Therefore, we construct a model for high-energy heavy-ion collisions based on relativistic resistive magnetohydrodynamics (RRMHD).
The RRMHD simulation has been performed in astrophysical situations such as black hole accretion disks, jets, and neutron star magnetospheres. The dissipation associated with Ohmic conduction current affects the topology of the magnetic field and liberates the magnetic energy, which is known as magnetic reconnection discussed in astrophysical applications [17][18][19][20][21][22].
In astrophysical community, there are two known difficulties for constructing the RRMHD code. First one is that a time step size t becomes too much short for calculation of highly conducting medium with explicit time integration. Then it is time consuming to follow a long time evolution of dynamics. Several numerical schemes are proposed to solve this problem such as adopting a semi-analytic solution [23] and an implicit-explicit (IMEX) Runge-Kutta algorithm [24,25]. Then, the time step size is determined by the advection time of the fluid or the wave rather than the diffusion time. Another difficulty is found in keeping the constraints of the Maxwell equations, e.g., Gauss's law and the divergence-free constraint. Ampere's law and Faraday's law describe the time evolution of electric and magnetic fields, respectively. These equations ensure that constraints hold if they are satisfied at the initial state. These conditions are, however, not preserved numerically in a not well-designed scheme for multi-dimensional problems. The violation of the constraints leads to the growth of spurious oscillation and the simulation is unexpectedly terminated. In the framework of RRMHD, the generalized Lagrange multiplier (GLM) method is proposed to guarantee these conditions [23,26]. In this scheme, the numerical errors relating to the violation of the constraints are advected and dissipated. The GLM scheme is commonly used also in MHD [27] and RMHD [28] and shows the robustness of the numerical code. As another solution of this problem, a constrained transport (CT) method [29][30][31][32][33] is also developed. This approach takes the magnetic field as a staggered representation whereby the different components are situated on the cell surface which is normal to the corresponding components.
In this paper, we construct a RRMHD code by adopting the semi-analytic method for time integration and GLM method for preserving the conditions in the Milne coordinates for the purpose of understanding high-energy heavy-ion collisions. These methods are well-known in the astrophysical community, but we adopted them in the Milne coordinates to understand the high-energy heavy-ion collisions for the first time. We check the correctness and robustness of our code by performing several numerical tests which are common test problems for relativistic ideal MHD and RRMHD. We also propose the test problem of longitudinal expansion with an acceleration of relativistic resistive magnetohydrodynamics. The semi-analytic solution of this problem is shown by mimicking the electromagnetic configuration and fluid velocity in high-energy heavy-ion collisions [34]. We apply this solution as a test problem of RRMHD in the Milne coordinates.
This paper is organized as follows. In Sect. 2, we briefly review the formulation of RRMHD. We explain the numerical algorithm of the RRMHD systems in our simulation code in Sect. 3. We show numerical benchmark tests as a verification of our RRMHD simulation code in Sect. 4. A summary is given at the end in Sect. 5. Unless otherwise specified, we use natural unitsh = c = 0 = μ 0 = 1, where 0 and μ 0 are the electric permittivity and the magnetic permeability in vacuum, respectively. Throughout the paper, the components of the four tensors are indicated with Greek indices, whereas three vectors are denoted as boldface symbols.

Relativistic resistive magneto-hydrodynamics
The conservation laws for the number density current N μ and for the total energy-momentum tensor of the plasma T μν in the dynamics of the whole system are written by, where ∇ μ is the covariant derivative. The electromagnetic fields follow Maxwell equations, where F μν is a Faraday tensor and F μν = 1 2 μνρσ F ρσ is its dual tensor, with μνρσ = (−g) −1/2 [μνρσ ], g = det(g μν ) and [μνρσ ] is a completely anti-symmetric tensor. Here g μν is a metric tensor. If the magnetization and polarization effects are ignored, the energy-momentum tensor of the electromagnetic fields is written by, and this tensor follows ∇ μ T μν f = J μ F μν , from Maxwell equations. The total energy-momentum tensor is the sum of the contribution of matter and electromagnetic fields T μν = T μν m + T μν f . The conservation law of the total system shown by Eq. (2) gives, For the ideal fluid, the energy-momentum tensor and the charged current of fluids are written by, where u μ (u μ u μ = −1) is a single fluid four-velocity, ρ is the number density, e = T μν m u μ u ν is energy density and p = 1 3 μν T μν m is pressure of the fluid. We have introduced the projection tensor μν = g μν +u μ u ν . The Faraday tensor and its dual tensor can be rewritten as, where, are the electric and magnetic fields measured in the comoving frame of the fluid. We introduce the electric field E i and the magnetic field B i , which are related to e μ and b μ by where γ is the Lorentz factor of the fluid and v i is the fluid three-velocity. The projections of ∇ μ T μν m = −J μ F μν along the perpendicular and parallel directions with respect to u μ are given by, where D = u μ ∇ μ and = ∇ μ u μ . Equations (15) and (16) correspond to the equation of motion and energy equation, respectively. They are auxiliary equations that are not solved in our numerical code. However, they are useful for the discussion of Bjorken flow in resistive medium, shown in Sects. 4.6 and 4.7.
Since the system of equations Eqs. (1)-(4) is closed by Ohm's law, we adopt the simplest form of it shown in Ref. [35]. In the covariant form, Ohm's law is written by, where σ is electrical conductivity and q = −J μ u μ is electric charge density of the fluid in the comoving frame. Maxwell equations lead to the charge conservation law, When we take the ideal limit (σ → ∞) of Ohm's law, Eq. (17) is reduced to, 3 Numerical procedure We now represent governing equations of the RRMHD in a suitable form for numerical calculation. Also, we show a numerical scheme to solve RRMHD equations which are appropriate for studying high-energy heavy-ion collisions.

Metric
We split the space-time into 3 + 1 components by spacelike hypersurface defined as the iso-surfaces of a scalar time function t and assume a metric of the form, Since we consider only the Cartesian and the Milne coordinates in this paper, the metric tensor is diagonal. We will use the Cartesian coordinates (t, x, y, z) in Sects. 4.1, 4.2, 4.3, 4.4 and 4.5.1. The Milne coordinates (τ, x, y, η s ) will be adopted in Sects. 4.5.2, 4.6 and 4.7, where τ := √ t 2 − z 2 is the longitudinal proper-time and η s := 1 2 ln t+z t−z is the space rapidity. We note that in the Cartesian coordinates, the three metric g i j = diag{1, 1, 1}, with √ −g = 1, whereas in the Milne coordinates, g i j = diag{1, 1, τ 2 }, with √ −g = τ . In both coordinates, ∂ j g ik = 0 ( j = 1, 2, 3), source terms of space-components in conservation equations vanish. However, in the Milne coordinates, where g 33 = τ 2 , the source terms for the energy conservation equation contain a nonzero term proportional to 1 2 ∂ 0 g 33 = τ .

Constraint equations
Maxwell equations contain two constraint equations, Though Maxwell equations ensure that these constraints are satisfied at all time steps, a simple integration of Maxwell equations in numerical simulation does not preserve these conditions because of numerical error. The accumulation of numerical errors leads to unphysical oscillation and numerical simulation crashes in the end for multi-dimensional problems. For this reason, a number of numerical techniques for avoiding this problem are investigated. In this model, we adopt the GLM method to guarantee these conditions [23,26,28]. The main idea is that one introduces two variables, ψ and φ as the deviation from constraints. One manages a system of equations to decay or carry the deviation ψ and φ out of the computation domain by relatively high-speed waves.
For the GLM method, we extend the Eqs. (3) and (4) to, where ψ and φ are new variables and κ is a positive constant. In the Cartesian coordinates, we obtain the telegraph equations for ∇ · E − q and ∇ · B, Consequently, ∇ · E − q and ∇ · B propagate at the speed of light and decay exponentially over a timescale 1/κ. In the form of the metric Eq. (20), a time-like normal vector is reduced to n μ = (−1, 0). The modified divergence-free equation of B and the Faraday's law become, Also, the modified Gauss's law and Ampere's law are obtained as, The conservation law of electric charge is written by,

Basic equations
Let us rewrite the equations of motion Eqs. (1), (2), (23), (24) and (31) in a conservative form which is appropriate for numerical integration, where U, F i , S e and S s are the sets of conservative variables, fluxes, source terms which are explicitly solved, and source terms of the stiff part, respectively. These variables contain the following components, where the number density D, the total momentum i , the stress tensor T i j , the total energy density ε and the conduction electric current J i c are given by, where the electromagnetic energy density p em is defined as . We adopt the operator splitting method for time integration. Then Eq. (32) can be divided into two equations as, Equation (39) can be integrated in time with an explicit manner, while Eq. (40) becomes stiff for large σ and κ. The one-dimensional discretization of Eq. (39) can be written by, where x is the grid spacing and f is the numerical flux. The subscript i denotes the grid point, x = i x, and the superscript n shows the number of time steps, t = n t. We need to reconstruct the primitive variables on the cell surface from those of the cell center to evaluate the numerical flux. We adopt the second-order accurate scheme for this reconstruction [36] given by, Then the numerical flux is computed by using the Harten-Lax-van Leer (HLL) method [37] given by, where λ + and λ − represent the maximum and minimum wave speeds of the system, respectively. We evaluate them by the speed of light for simplicity. The superscript L (R) denotes that the variables are calculated from the primitive variables reconstructed by these at grid on the left (right) side. Then, we numerically integrate Eq. (39) using numerical fluxes.

Stiff part
In this subsection, we show how to solve the stiff part Eq. (40) based on Ref. [23], which contains Ampere's law and equations of φ and ψ. For Ampere's law, we split the stiff relaxation equation Eq. (40) into the perpendicular and the parallel directions with respect to v i by the operator splitting method and we redefine where E i ⊥ and E i are the perpendicular and the parallel directions of the electric field with respect to v i . Here, v i and γ are assumed to be a constant during a small time step t. The solutions of the initial value problem for these relaxation equations are given by, where E * i ⊥ = − i jk v j B k and suffix 0 represents the initial value of E i .
In this scheme, the fluid velocity and Lorentz factor are assumed to be constant within the time step t. Thus our numerical scheme reduces to a first-order accuracy in time.
The formal solution of Ampere's law should be adopted in a small time step because the velocity and Lorentz factor are not the constant during t. Hence, the numerical accuracy worsens if the velocity drastically changes with time. We note that the higher order accuracy scheme is proposed by the previous studies, e.g., Ref. [25].

Primitive recovery
In our numerical code, we need to reconstruct primitive variables {ρ, p, v i } from evolved conservative variables {D, ε, i } in order to calculate the numerical flux at each time step. To reconstruct the primitive variables, we introduce new variables, where ε and are the energy density and the momentum of fluid, respectively. The set of the variables {D, ε , i } is the relativistic ideal fluid conservative variables. Therefore, we adopt the ordinary primitive reconstruction method of relativistic hydrodynamic simulation for these variables [38]. We obtain a one-dimensional equation about the gas pressure, where, The p and ρ dependence of e is determined by an equation of state (EoS). This equation is numerically solved by the Newton-Raphson algorithm. The other primitive variables are reconstructed as, We consider the ideal gas EoS p = ( − 1)(e − mρ) with the mass of conserved species m = 1 in this paper. We note that we take the non-zero ρ in Sects. 4.1, 4.2, 4.3, 4.4, and 4.5.1. In Sects. 4.5.2, 4.6, and 4.7, we set ρ = 0.

Numerical algorithm
Our numerical simulation code is based on finite difference schemes. The advection equation Eq. (39) including the source term S e is solved by explicit time integration with the second-order TVD Runge-Kutta algorithm [39]. The primitive variables are interpolated from the cell center to the cell surface by using the second-order accurate scheme represented by Eqs. (42) and (43) [36]. The numerical flux is calculated by the HLL flux shown in Eq. (44) [37]. The stiff part Eq. (40) is integrated by the analytic solutions explained in Sect. 3.4. We adopt the GLM method to guarantee these conditions [23,26,28] as the divergence cleaning. The primitive recovery is executed by solving Eq. (53) by the Newton-Raphson algorithm explained in Sect. 3.5 [38]. We summarize our numerical algorithm as follows: -the values of the primitive variables are reconstructed from the cell center to the cell surface by using the second-order accurate scheme [36] shown in Eqs. (42) and (43).

Test problems
In this section, numerical tests are performed for the verification of our numerical simulation code. We will use the Cartesian coordinates (t, x, y, z) in Sects. 4.1, 4.2, 4.3, 4.4 and 4.5.1. The Milne coordinates (τ, x, y, η s ) will be adopted in Sects. 4.5.2, 4.6 and 4.7. We take the CFL constant as C CFL = 0.1 for all test problems. In all tests, φ and ψ are initialized as zero at all computational grids. The artificial constant κ for two constraints is determined by the relation α := hκ/c h where h = min( x, y, z) and the maximum speed of system c h = 1. We set κ = 5.5 in all multidimensional tests [23], which keeps α ≤ 1.

Shock tube test problem
In order to test the shock-capturing properties of our numerical simulation code, we consider the simple MHD version of the Brio-Wu test in Ref. [25]. The initial left and right states are given by, Fig. 1 We display the magnetic field B y at t = 0.4 in the Brio-Wu type shock tube test problem and all the other variables are initially set to 0. Figure 1 shows the results at t = 0.4 for σ = 0, 10, 10 2 , 10 3 , and 10 4 . The simulation box is bounded by x ∈ [0, 1.0] and the number of grid points is 400. For σ = 10 4 , the left going rarefaction and right going shock wave, and a tangential discontinuity between them appear. Our numerical solution with σ = 10 4 reproduces the solution of relativistic ideal MHD simulations [40]. In addition, the solutions with low conductivity are similar to the results of the other RRMHD numerical simulations [25,41]. Especially the electromagnetic field does not interact with gas for σ = 0, so that the electromagnetic waves propagate with speed of light. The wave fronts should be located at 0.1 and 0.9 for left and right going waves, respectively. Our results are consistent with these analytic solutions, although wave fronts slightly have a smooth profile due to the numerical diffusion.

Large amplitude circularly polarized Alfvén wave
We perform a test for the propagation of a large amplitude circularly polarized Alfven waves along a uniform background magnetic field B 0 . The analytical solution of relativistic ideal MHD is proposed in Ref. [42], which is given by, where B x = B 0 , v x = 0, k = 2π is the wave number, and η A is the amplitude of the wave. The square of the special relativistic Alfven wave's speed v 2 A is given by, where h = (e + p)/ρ is the specific enthalpy. We set the initial parameters as ρ = p = η A = 1, and B 0 = 1.1547. We use an ideal gas EoS, p = ( − 1)(e − mρ) with = 2 in this test. From these parameters, the Alfven wave's speed is estimated to v A = 1/2.

One-dimensional propagation
For one-dimensional propagation, the computational domain is x ∈ [−0.5, 0.5]. We use a periodic boundary condition. Since the solution given by Eqs. (61) and (62) is obtained by solving special relativistic ideal MHD equations, we take a sufficiently large conductivity σ = 10 6 in this test problem. Figure  The results with high conductivity σ mean that our simulation code can handle the ideal MHD limit. Also, there is no dependence of the number of grid points, which suggests that even 100 grid points are enough for the description of the analytical solution.
For this problem, we cannot achieve full second-order accuracy. Figure 3 shows the L 1 (u, x −1 ) norm errors of the tangential magnetic field B y as a function of the inverse of grid-cell size x −1 , where u EXACT (x i ) is an exact solution at x = x i . The behavior of L 1 norm shows that our numerical simulation is nearly 1.6-order convergence. The main reason for it is that we execute the second-order Runge-Kutta algorithm with many operator splittings. It makes the time accuracy of our scheme worsen [41]. This is the most difficulty to solve in RRMHD since this is the limit of large electrical conductivity.

Two-dimensional propagation
We consider the two-dimensional propagation of large amplitude circularly polarized Alfven waves to confirm the numerical convergence in the multi-dimensional test. The initial vector field such as velocity and the magnetic field is rotated in the x − y plane from one-dimensional propagation. The wave vector k = (k x , k y , k z ) in the three-dimensional Cartesian coordinates is determined by the rotation angle α, where k z is set to be zero since we consider two-dimensional propagation. Then, the vector A is rotated to We set k x = 2π and the computational domain x ∈ [−0.5, 0.5] × y ∈ [−0.5/ tan α, 0.5/ tan α]. The rotation angle is set to be α = tan −1 2. The number of grid points in y-direction is taken to be N y = N x /2 where N x is the number of grid points in x-direction. We take the periodic boundary condition in each direction. The artificial constant κ for the GLM method is chosen to be κ = 5.5 in each case. Figure 4 shows the z-component of magnetic field at t = 2.2 and y = 0. The green long dashed-dotted, red dashed, and black solid lines represent results with σ = 10 2 and N x = {100, 200, 400}, respectively. In the low conductivity case, there is no analytical solution. To see the convergence rate, we compare the reference case (N x = 400) with the other low-resolution cases. We calculate the convergence rate using the equation, We confirm that the convergence rate is achieved with 1.3order convergence. This convergence rate is the same as that of the one-dimensional propagation test with N = {100, 200, 400}.

Self-similar current sheet
Next, we consider the evolution of the self-similar current sheet proposed in Refs. [23,25]. In this test problem, we assume that the magnetic pressure is much smaller than the pressure of fluid. The magnetic field has only a tangential component B = (0, B(x, t), 0) and B(x, t) changes the sign across the current sheet. The initial pressure of background fluid is set to be uniform, p = const. We take the high conductivity σ , so that the diffusion timescale is much longer than the light crossing timescale. In this assumption, the evolution equation of the magnetic field is reduced to, The analytic solution of this equation is given by, where erf is the error function. We set the initial condition at t = 1 with p = 50, ρ = 1, E = v = 0, and σ = 100. The computational domain is x ∈ [−1.5, 1.5], and the number of grid points is N = 200. Figure 5 shows the numerical result at t = 9. The black long dashed-dotted line represents the initial profile of B y at t = 1. The blue solid and red dotted

Cylindrical explosion
The symmetric explosions are useful standard tests for MHD codes even though there are no exact solutions because of the existence of the shock waves in all possible angles [23]. This problem including a strong multidimensional shock wave is useful to check the robustness of the code. We take the same condition as that in Ref. [23]. The computational domain is x ∈ [−6.0, 6.0] × y ∈ [−6.0, 6.0] in the two-dimensional Cartesian coordinates. The number of grid points is 200×200 with the uniform grid. The initial cylinder radius is taken to r = 1 centered at the origin, where r = x 2 + y 2 . The pressure and number density of fluid are set to p = 1 and ρ = 0.01 for r ≤ 0.8 and exponentially decrease with increasing radius for 0.8 ≤ r < 1.0. The fluid in the exterior of the initial cylinder has p = ρ = 0.001 for r ≥ 1.0. We take the uniform magnetic field, B = (0.1, 0.0, 0.0), and the velocity v = 0 at the initial state. The plasma resistivity and decay constant for variables ψ and φ are set to be σ = 55 and κ = 5.5. This calculation with σ = 55 is the case of the highest possible value of σ before our code crashes. This reason is that the spatial resolution h = 0.06 is slightly larger than τ σ = 0.018. Figures 6 and 7 show the two dimensional profiles of p and B x at t = 4, respectively. We can see that strong shock waves are formed and the initial uniform magnetic fields are deformed due to the cylindrical explosion. Such a deformation of the magnetic field may lead to the violation of constraints ∇ · B = 0 and ∇ · E = q. If these constraints are violated, the calculation would crash due to the growth of unphysical spurious oscillation. We can avoid these problems by adopting the GLM method. Figure 8 represents the two dimensional profile of Lorentz factor γ at t = 4. The  profile of γ is not isotropic at this time since the magnetic force decelerates gas in y-direction. Figure 9 stands for the two-dimensional profile of the 10 3 × φ at t = 4. The magnitude of φ is very small, |φ| ∼ 10 −3 . The front of φ moves to the outside of the computational domain faster than the expansion of fluids does. We have confirmed that the ψ is zero within a round error in all computational domains. The results show the robustness of our code in multidimensional problems. We also note that our results are consistent with those of other groups (e.g., Ref. [23]). We get the same results on x−z and y−z planes.

Rotor test
The rotor test is important for a calibration of resistive as well as ideal MHD numerical multidimensional codes [43,44]. We perform the resistive rotor test both in the Cartesian and in the Milne coordinates.

Cartesian coordinates
The Cartesian computational domain is x ∈ [−0.5, 0.5] × y ∈ [−0.5, 0.5] with 300 equidistant grid points in each direction. Inside r ≤ 0.1 centered at the origin, the fluid rotates with constant angular velocity = 8.5 with uniform number density ρ = 10 at the initial state. Outside this region (r ≥ 0.1), the medium is static and uniform (ρ = 1). Both the pressure ( p = 1) and the magnetic field B = (1.0, 0, 0) are uniform in the whole region. The initial electric field is given by the ideal condition, −v × B. The adiabatic index is = 4/3. Figures 10 and 11 show snapshots of the gas pressure p and the electric field component E z at t = 0.3 with electrical conductivity σ = 10 6 . These results agree with those of other simulation codes [43][44][45]. Again, the GLM method does work in our code, so that the growth of unphysical oscillation suppresses. The same results are reproduced by taking x−z and y−z planes.

Milne coordinates
The rotor test problem for the relativistic ideal MHD in the Milne coordinates is proposed in Ref. [15]. The computational domain in the Milne coordinates is x ∈ [−0.5, 0.5] × y ∈ [−0.5, 0.5] with 400 equidistant grid points in each Fig. 10 The pressure of the fluid at t = 0.3 is displayed in the 2D resistive rotor test in the Cartesian coordinates Fig. 11 The electric field component E z at t = 0.3 is shown in the 2D resistive rotor test in the Cartesian coordinates direction. We adopt the free boundary condition in both directions. In this test, the EoS is assumed to be the ultrarelativistic ideal gas EoS, p = e/3, and the number density is set to ρ = 0 at all grid points. We set the initial coordinate time to t = 1.0. Instead of imposing high density around the origin, we set a high pressure region p = 5 inside r = 0.1. The initial transverse velocity (v x , v y ) is given by angular speed = 9.7. Outside this region, the pressure is p = 1. The initial longitudinal velocity is set to v η = 0, which is assumed to be longitudinal Bjorken expansion amount to v z = z/t. The initial magnetic field is set to B = (2.0, 0, 0) and the electric field is taken to be −v × B. Figure 12a and b show the results of p and p em of our simulation code at t = 1.4 with σ = 10 3 , respectively. In the Milne coordinates, we observe the decay of the pressure of the fluid and magnetic pressure, which occurs in the whole computational domain because of the longitudinal expansion of the system. An asymmetric shaped compression wave by the magnetic field is observed in our results because of the higher initial pressure inside the radius r < 0.1 and the rotation of the cylinder. Our result with σ = 10 3 has a similar configuration of the result in ECHO-QGP simulation [15]. In the high conductive case, our results capture the features of relativistic ideal MHD. Figure 13a and b represent the results of p and p em of our simulation code at t = 1.4 with σ = 10, respectively. We can see that the anisotropy becomes weaker for a lower conductivity. It indicates that the fluid is weakly coupled with electromagnetic fields. The fluid vorticity does not affect the dynamics of electromagnetic fields. These features can be captured in our code. This problem is an extension of the one-dimensional boost invariant flow proposed by Bjorken [46] to the ideal MHD [47]. We consider the relativistic boost invariant flow in zdirection. In magnetized Bjorken flow, we consider that the fluid follows the ultra-relativistic ideal gas EoS, p = e/3. The pressure of fluid is assumed to be uniform in space. The fluid velocity is given by v z = z/t, which is obtained by assuming a longitudinal boost invariant expansion. The corresponding four-velocity in the Cartesian coordinates is given by u μ = (cosh η s , 0, 0, sinh η s ). In the Milne coordinates, the four-velocity becomes simply u μ = (1, 0, 0, 0). The comoving derivative and the expansion rate become D = ∂ τ and = 1/τ , respectively. The magnetic field is given by b μ = (0, b x , b y , 0). In this assumption, from the energy conservation equation Eq. (2), one derives [47], where b 2 = b μ b μ . From the energy equation Eq. (16), we obtain, for the ideal MHD. The evolution equation of the magnetic field is obtained as, The analytic solutions of these equations for the ultrarelativistic ideal gas EoS, p = e/3, are given by, where τ 0 = 0.5 fm is initial time, e 0 = 10 GeV/fm 3 is initial energy density, and b 0 = 1.0 GeV 1/2 /fm 3/2 is initial magnetic field. Figures 14 and 15 show the time evolution of energy density and magnetic field strength with σ = 10 6 , respectively. The solid curves show analytic solutions, while the dashed curves denote the numerical results. These results are in good agreement with the analytic solutions.

Accelerating longitudinal expansion
This test problem is presented in Ref. [34]. This is a resistive extension of the magnetized Bjorken flow [47]. The derivation of the semi-analytic solutions is shown in Appendix A in detail, but we here show the results. In this problem, we do not suppose the boost invariant flow. We assume that the fluid velocity is parallel to the longitudinal (z-) direction, and the fluid is uniform on the transverse (x-y) plane. Let where Y is the rapidity withγ The electric and magnetic four vectors are considered on the transverse plane and orthogonal to each other, In Ref. [34], electromagnetic fields are taken to be the following forms, The solutions for the rapidity and electromagnetic fields are given by, Here c(η s ) is an arbitrary function. We suppose, which is the same as that in Ref. [34]. Here α is an arbitrary constant. The arbitrary constant c 0 is taken to be 0.0018τ 0 , which determines the initial magnetic field strength on the laboratory frame in the Minkowski coordinates B y L (τ 0 , 0) = 0.0018 GeV 2 /e.
The energy density is determined by solving the following equations, where, Here, we take the ideal gas EoS, p = κe and κ = 1/3. We numerically solve these differential equations on a grid of points on the (τ, η s ) plane. First, we give the arbitrary function c(η s ), which is assumed to have a form given in Eq. (84). Then the rapidity, four-velocity, and electromagnetic fields are obtained from Eqs. Here, we perform the RRMHD simulation using the numerical solution of the above ODEs as the initial condition. Since the resistive effect is prominent with small α [34], we take the initial energy density e 0 = 1.0 GeV/fm 3 , α = 0.1, and the electrical conductivity σ = 0.023 fm −1 . The number of grid points in the computational domain η s ∈ [−3.0, 3.0] is N = 200. We adopt free boundary conditions at η s = ±3.0. The waves are, however, sometimes reflected at the boundaries and they affect the numerical results. To avoid this, we take the boundaries far away from the central region η s ∈ [−1.0, 1.0] and stop the calculation before the reflected waves reach the central region. Figures 16 and 17 show the results of our simulation code in comparison with the semi-analytic solution. Our results are in good agreement with the semi-analytic solutions. In Fig. 16a, the energy density decays and expands to the lon- gitudinal direction by the resistive effects. In Fig. 16b, the fluid velocity in the Milne coordinates has a finite value in the forward and backward rapidity regions and it decays with time. This feature comes from the resistive effect described in Eq. (81). Our results reproduce this behavior of the fluid. In Fig. 17a, the magnetic field component in the comoving frame b y as a function of η s is represented. The magnetic field decays with time by the longitudinal expansion similar to the magnetic Bjorken flow in Sect. 4.6. However, the magnetic field has a non-uniform profile in the η s direction by the acceleration induced from the resistive effect, which is different from the magnetic Bjorken flow. Figure 17b shows the electric field component measured in the coming frame e x as a function of η s . The electric field has a positive value in the backward rapidity region and it decreases with rapidity. The electric field changes its sign at η s = 0. This feature describes that the electric field is produced by the two colliding nuclei in high-energy heavy-ion collisions. Our results capture these features and its diffusion which is consistent with the semi-analytic solutions. Figure 18 represents the electrical conductivity dependence of u η s /u τ at τ = 1 fm. The blue, red, and black dotted lines show the semi-analytic solutions given Eqs. (80) and (81) in the cases of σ = 0.023, 0.033 and 0.043 fm −1 , respectively. All our numerical results are in good agreement with the semi-analytic solutions. The magnitude of u η s /u τ in forward and backward rapidity regions decreases with increasing σ . This means that the case of higher conductivity is close to the magnetized Bjorken flow explained in Sect. 4.6. In fact, if we take the limit of infinite electrical conductivity, Eq. (81) becomes zero, which is consistent with magnetized Bjorken flow. Figure 19 shows the electrical conductivity dependence of e x at τ = 1 fm. The blue, red, and black dotted lines show the semi-analytic solutions of the Eq. (82) in the cases of σ = 0.023, 0.033 and 0.043 fm −1 , respectively. The electrical conductivity dependence of e x in our results is in good . The magnitude of the electric field in forward and backward rapidity regions decreases with σ . This is consistent with the ideal limit of RRMHD. In the relativistic ideal MHD, the electric field measured in the comoving frame of fluid vanishes. Our numerical results capture this feature of ideal limit of RRMHD.

Summary and discussion
We constructed a new RRMHD simulation code in the Milne coordinates which is suitable for study on high-energy heavyion collisions. We split the system of RRMHD equations into two parts, a non-stiff and a stiff part. The primitive variables are interpolated from the cell center to the cell surface by using the second-order accurate scheme [36]. For the nonstiff part, we evaluated the numerical flux using the HLL approximated Riemann solver and explicitly integrated the equations in time by the second-order of Runge-Kutta algorithm [39]. For the stiff part appeared in Ampere's law, we executed time integration using semi-analytic solutions to avoid unexpected small time steps. Though Maxwell equations ensure that the divergence-free constraints are satisfied at all times, in numerical simulation, the integration of Maxwell equations in a not well-designed scheme does not preserve these conditions because of the numerical error. In order to avoid this problem, we employed the generalized Lagrange multiplier method to guarantee these conditions [23,26,28].
We checked the correctness of our algorithm from the comparison between numerical calculations and analytical solutions or the other RMHD simulations such as Brio-Wu type shock tubes, propagation of the large amplitude circularly polarized Alfven waves, self-similar current sheet, cylindrical explosion, resistive rotor, and Bjorken flow. In these test problems, our numerical solutions were in good agreement with analytic solutions or results of the other RRMHD simulations. Furthermore, we investigated the accelerating longitudinal expansion of relativistic resistive magneto-hydrodynamics in high-energy heavy-ion collisions in comparison with semi-analytic solutions [34]. Our numerical code reproduced these solutions. We conclude that our numerical simulations capture the characteristic features of dynamics in high-energy heavy-ion collisions.
We shall employ our RRMHD code to investigate experimental data at RHIC and the LHC and understand the detailed QGP bulk property such as electrical conductivity. As the first application, we have demonstrated calculation of the directed flow in RRMHD expansion in Au-Au and Cu-Au collision systems at RHIC [48]. We found that electrical conductivity prevents the growth of the directed flow in asymmetric systems. The Ohmic conduction current may affect the charge-dependent flow of hadrons. Recently, the charge-dependent directed flow of the heavy mesons and dileptons has been proposed as a signature of the electromagnetic fields [49,50]. It can extract the initial electromagnetic fields and the QGP bulk properties from experimental results.
The highly intense electromagnetic fields induce novel quantum phenomena such as the chiral magnetic effect (CME) [51,52] and the chiral magnetic wave (CMW) [53]. The CME and CMW signals in experimental data have been explored in iso-bar collisions, Zr-Zr and Ru-Ru collisions at RHIC [54,55]. However, there are no striking evidences of these phenomena in the iso-bar experiments. For analytical studies, the initial condition of electromagnetic fields with CME in high-energy heavy-ion collisions has been studied [56,57], but the evolution of these fields has not been well understood. The recent development of numerical simulation for the chiral RMHD has been investigated [58]. The comprehensive studies using RRMHD would shed light on these problems. This issue remains an important future work.
are uniform in the transverse plane. Let us parametrize the four-velocity in (1 + 1)D as follows: where Y is the rapidity and v z = tanh Y . In the Milne coordinates, u μ is rewritten by, Under this parameterization, the comoving derivative and the expansion rate are given by, (A.14) In Ref. [34], in order to solve these equations, we can take the following Ansatz: with the arbitrary parameters α and c 0 .