Efficient solver of relativistic hydrodynamics with implicit Runge-Kutta method

We propose a new method to solve the relativistic hydrodynamic equations based on an implicit Runge-Kutta method with an optimized fixed-point iterative solver. In the case of ideal hydrodynamics, the accuracy and computational cost of our new method are compared with those of explicit ones for the (1+1)-dimensional Riemann problem, as well as the (2+1)-dimensional Gubser flow and event-by-event initial conditions for heavy-ion collisions generated by TrENTo. We demonstrate that the solver converges with only one iteration in most cases, and as a result, the implicit method requires a smaller computational cost than the explicit one at the same accuracy in these cases.


Introduction
Relativistic hydrodynamics is a versatile tool for describing various long-range phenomena from astrophysics to nuclear physics at high energies.One of the most important applications is relativistic heavy-ion collisions (HIC).These collision experiments of heavy nuclei at ultrarelativistic energies are currently operated at the Large Hadron Collider (LHC) at CERN, the Relativistic Heavy Ion Collider (RHIC) at BNL, the Super Proton Synchrotron (SPS) at CERN, and the Heavy Ion Synchrotron SIS18 at GSI, and several other experiments are planned at large facilities all over the world.In the hot and dense reaction zone, the strongly coupled quark-gluon plasma (QGP) is created [1,2,3], which is well described by a perfect fluid.To investigate the properties of the QGP created in heavy-ion collisions, we need to perform the simulation of heavyion collision reactions based on the dynamical models and compare the results with high-statistics experimental data [4,5,6,7,8,9,10].In modern modeling of heavy-ion collisions, a relativistic viscous hydrodynamic code [11,12,13,14,15,16,17,18,19,20,21,22,23] forms the core of the dynamical description of the reactions together with an event-by-event initialstate model and a final-state hadronic cascade.
To carry out those simulations, large computing resources and high performances of numerical codes are required.While in nature the evolution of a heavy-ion collision takes several 10 fm from the initial state to the kinetic freezeout of final particles, a numerical simulation of a single event can take up to several hours depending on the computational devices.In addition, a large number of collision events are required, particularly, for the ongoing interest in global Bayesian analyses [24,25,26,27,28,29] and higher-order fluctuation observables, such as the skewness and the kurtosis of the event-by-event distributions of conserved charges [30,31,32,33].To properly extract the properties of the QGP, both large computational resources and high-performance codes are important.In the present study, we focus on the latter: the potential improvement of the numerical solver.
To efficiently solve relativistic hydrodynamics for heavy-ion collisions while preserving the needed accuracy, various space-discretization schemes have been implemented, such as SHASTA (sharp and smooth transport algorithm) [11,12,14,34,35], KT (Kurganov-Tadmor) [13,20,22,23,36], HLLE (Harten-Lax-van Leer-Einfeldt) [37,16,17,38,19], and other Gudunovtype schemes [18,21].In this study, we instead explore the time discretization.The time integration methods are classified into two, explicit and implicit methods.In an explicit method, the state of the next time step is obtained by straightforward evaluations of given expressions.In an implicit method, we need to solve implicit equations to obtain the next state.To the best of the authors' knowledge, the time integration within all of the existing codes for HIC is performed by explicit methods.The application of the implicit methods to the hydrodynamics in HIC is an unexplored territory.The goal of our work is to seek the possibility of developing an implicit Runge-Kutta (RK) method that is stable, accurate, and sufficiently efficient.
The computational cost of implicit methods is generally considered much higher than explicit ones because implicit methods typically require numerically solving equations, which would be one of the reasons that they have not been applied in the hydrodynamics for HIC.However, it is only true for a single-step time integration.There is another factor that affects the overall computational cost to obtain the solution at the final time, i.e., the number of time steps required to achieve sufficient accuracy and stability.Here, implicit methods have an important feature called the A-stability: accumulated numerical errors are bounded for any positive time step, ∆t, for linear ordinary differential equations.One can, therefore, hope that stable calculations can be performed using an implicit method at a larger ∆t even for non-linear equations leading to a smaller computational cost.
Stability is also important in the implementation of computational hydrodynamics.Various technical procedures for stabilization, such as flux limiters and monotonicity preservation, need to be developed using non-trivial assumptions and careful tests for each.
Choosing a stable implicit method may reduce such an effort for the stabilization, which would help deal with more complex systems such as fluctuating hydrodynamics [39,19,40,41,42,43], non-equilibrium chiral fluid dynamics [44,45], and magnetohydrodynamics [16,46].More specifically, these complex systems are anticipated to be stiff systems, where explicit methods require impractically small ∆t due to highfrequency modes.Genuine implicit methods can stably solve such stiff equation systems with a moderate time step ∆t.
Motivated by these advantages of implicit methods, we implement ideal relativistic hydrodynamics with an implicit time integration and compare its performance with a conventional explicit method quantitatively.To this end, we employ the Kurganov-Tadmor scheme [47] for the space discretization, which is used for example in Refs.[13,22].This spatial scheme is independent of the time discretization and thus allows us to systematically compare different time-discretization schemes.We will especially use the one-stage Gauss-Legendre method as an implicit scheme and Heun's method as an explicit scheme.Both are of second-order accuracy.
The most numerically demanding part of our implicit methods is the solver for the Runge-Kutta nonlinear equations.In our code, we use a simple fixed-point method for the solver but with several improvements.For example, we use the solution found in the previous time step as the initial guess of the fixed-point solver in the next step.We also check the convergence of iterations cell by cell and skip updating the cells in which the iteration converged.We show that these optimizations enable an efficient implicit scheme while preserving the accuracy and the simplicity of implementation.
We consider three initial conditions for benchmarks: the Riemann problem [48], the Gubser flow [49], and the event-by-event T R ENTo initial conditions for heavyion collisions [50].We show that contrary to the typical expectations, our implicit method is usually more efficient and accurate than the explicit one even with the same ∆t for the cases of these initial conditions.This paper is organized as follows.In Sec. 2, we introduce the hydrodynamic equations and their analytic solutions that will be used in later sections.In Sec. 3, after a brief review of the implicit and explicit RK methods, we discuss our space and time discretization schemes.The numerical setup and analyses are given in Sec. 4. We then discuss the numerical results in Sec. 5.The final section is devoted to discussions and a summary.

Relativistic hydrodynamic
In this section, we summarize hydrodynamic equations and a few known analytical solutions that will be used in Sec. 5 for numerical tests.

Hydrodynamics equations
Hydrodynamic equations are given by the conservation law.We focus on the conservation of the energymomentum tensor T µν .In the Cartesian coordinate system with the flat metric g µν = diag(1, −1, −1, −1), it is written as The scheme described in this paper can be naturally extended to the general case including further conserved currents.
In this paper we consider the ideal hydrodynamics where the energy-momentum tensor is written in terms of thermodynamic quantities as where u µ is the four-velocity with u µ u µ = 1, ϵ the energy density, P the pressure, and ∆ µν = g µν − u µ u ν the projector orthogonal to u µ .Equation (2) contains five undetermined variables (u i , ϵ, P) with i = x, y, z.Since Eq. (1) has four equations, one needs another constraint to close the equations.This constraint is provided by the equation of state (EOS), which is the pressure as a function of the energy density P(ϵ).Equations ( 1) and ( 2) together with the EOS define the ideal hydrodynamics.
In the Cartesian coordinates, we choose the conserved densities T tt , T tx , T ty , and T tz as the state variables to be solved numerically.This is a requirement from the space discretization scheme as we will see later.In this treatment, we need an explicit way to calculate T iµ as a function of T tµ that satisfies Eq. ( 2).We first obtain ϵ and u µ from T tν and then calculate T iν using Eq.(2).To numerically determine ϵ and u µ satisfying Eq. ( 2), it is convenient to solve the following equations for the energy density ϵ and the velocity v iteratively: where Then, u µ is given by In high-energy heavy-ion collisions, the space-time evolution of the produced medium is well described by the Bjorken assumption of longitudinal boostinvariance [51].In this case, it is convenient to employ the Milne coordinates or written inversely, where z is the longitudinal direction parallel to the beam axis while the transverse coordinates, x and y, remain unchanged.The energy-momentum conservation in the Milne coordinate system reads where the last term is a geometric source term coming from the non-flat coordinate system with the metric To solve the hydrodynamic equations in the Milne coordinates, it is useful to rewrite Eq. (7) as with T µν = τT µν and to treat T µν as the state variables.
To describe the dynamics in the mid-rapidity region of a high-energy collision, it is also useful to reduce the spatial dimension to two by assuming u η = 0, T τη = 0, and the boost invariance (i.e., ∂ η = 0).It should be noted that in this case T ηη = P(ϵ)/τ 2 does not vanish even though T xη = T yη = 0.

Analytical solutions
In this subsection, we show two known analytical solutions of the hydrodynamic equations, to which we will compare our numerical results in later sections.

Riemann problem
Let us consider a system having translational invariance in the y and z directions in the Cartesian coordinates.We consider a Riemann problem [48] in the remaining x direction given by the initial condition, with a discontinuity at x = 0.The analytical solution is given for a particular choice of an EOS, P(ϵ) = c 2 s ϵ.Here, c s is the speed of sound.For t > 0, the solution ϵ(x, t) = ε( x) with x = x/t is given piecewise [48]: where The constants ϵ plateau and v shock = −v r are obtained by solving the following equations for (ϵ plateau , v l , v r ): where We note that Eqs. ( 16) and (17) represent the Taub equations [52,53].The solution ε( x) is illustrated in Fig. 1.The discontinuity at x = v shock is called the shock front, and the curve in the range −c s ≤ x < −v r-p is called the rarefaction wave.
In heavy-ion collisions, the created matter is surrounded by the vacuum, where ϵ(t, x) = 0, and the medium expands into it.Such a situation can be mimicked in the Riemann problem by setting ϵ min = 0, where the solution is simplified to

Gubser flow
Another exact solution is the Gubser flow [49], which is a solution for the conformal EOS, P(ϵ) = ϵ/3, characterized by the boost invariance and the azimuthal symmetry (i.e., the rotational symmetry in the x-y plane).The solution is given in the Milne coordinates as with where q and ϵ 0 are free parameters.

Discretization of hydrodynamic equations
To solve the hydrodynamic equations numerically, we must discretize both the spatial and temporal coordinates.In this section, we first give a brief review of the Runge-Kutta (RK) method with a focus on the differences between implicit and explicit methods in Sec.3.1.The space and time discretization schemes applied in our study are then discussed in Secs.3.2 and 3.3, respectively.

Implict and explicit Runge-Kutta methods
The RK method is a common way to numerically solve a system of ordinary differential equations (ODEs) of the general form, where ⃗ y(t) represents the full set of state variables.The RK method solves the variables ⃗ y(t end ) at the end time t = t end given an initial condition ⃗ y(t init ) at the initial time t = t init by discretizing the time into steps of length ∆t and iteratively updating ⃗ y(t) for every time step.For one time step from t to t + ∆t, the state variable of the next time ⃗ y(t + ∆t) is constructed from ⃗ y(t) as where S is the number of stages of the RK method, and the coefficients a nm , b n , c n for n, m = 1, • • • , S are determined so that the method has desired properties.One specific set of the coefficients defines a particular RK method.The table that summarizes these coefficients is called a Butcher table, whose examples are shown in Table 1.
In a RK method, the numerical error of ⃗ y(t end ) from the exact solution ⃗ y * (t end ) is suppressed in the limit ∆t → 0 as where C(t) is specific to the ODE and is independent of ∆t.The largest integer p satisfying Eq. ( 30) is called the order of accuracy that is specific to each RK method.
The RK methods are classified into two categories: explicit and implicit ones.A RK method is said to be explicit when a nm = 0 for m ≥ n.In this case, ⃗ k (n) in Eq. ( 29) depend only on ⃗ y(t) and ⃗ k (m) for m < n, so ⃗ k (n) are sequentially obtained from n = 1 to S by simply substituting the known values.The other RK methods are said to be implicit.For implicit methods, one must solve the set of equations Eq. ( 29) for ⃗ k (n) (t), which is typically non-linear and practically carried out by a numerical solver.The implicit methods are thus usually considered to be computationally more expensive.
Despite this disadvantage, the implicit methods have a useful property, i.e., the implicit methods are always linearly stable (the A-stability) as far as the exact solution is not divergent [54].This means that a larger ∆t can be adopted for stiff equations in implicit methods while keeping numerical error finite.Another property of the implicit methods is that they can reach a given order of accuracy in fewer stages than explicit ones.More specifically, explicit methods require as many stages as the order up to p = 4 and more stages than the order for p > 4. On the other hand, the maximum accuracy for a given S is p = 2S in implicit methods [55].

Space discretization
In computational hydrodynamics, we discretize the spatial coordinates to represent the state of the fluid using a finite number of degrees.The spatial derivatives are then approximated in terms of discretized fields so that the partial differential equations become a set of ODEs, which can be solved by a RK method.This discretization has to be carefully chosen to suppress numerical diffusion but at the same time to avoid spurious oscillations.
In this study, we employ for the space discretization the Kurganov-Tadmor (KT) scheme with the MUSCL (monotonic upstream-centered scheme for conservation laws) reconstruction [47].This method has been applied to relativistic hydrodynamics for HIC in Refs.[13,22].The KT method has second-order accuracy for space discretization and is designed to capture discontinuities accurately and suppress spurious oscillations.Another important property of the KT scheme for our study is its independence of time discretization.This property allows us to combine this scheme with arbitrary time integration methods and compare them.
To discretize the space in the Cartesian coordinates, we assume a structured grid where the space is discretized into cells by an equal mesh size ∆x along each coordinate.In the KT scheme, one chooses the set of the conserved densities of the cells, {[T tν ] j }, as the state variables, where the subscript j of the square bracket specifies a cell.We spatially discretize the continuity equation with using the KT operators KT i .In our case, the KT operator in the x direction, KT x , is written as where H j+1/2 is the flux at the boundary,   3) and ( 4).The maximum local speed of propagation is given by where J αβ (T tν ) = ∂T xα /∂T tβ | T tν is the 4 × 4 Jacobian matrix, and ρ(M) is the spectral radius of a matrix M 1 .
In the original study of the KT scheme [47], the quantities on the cell boundary Q ± j+1/2 are calculated from the cell values Q j using the MUSCL reconstruction as The flux limiter in Eq. ( 38) guarantees the nonoscillatory property of Q with the parameter θ restricted to 1 ≤ θ ≤ 2. We choose θ = 1.1, which gives nonoscillatory results with a limited numerical diffusion.
In relativistic hydrodynamics, the reconstruction could be simply applied by choosing Q = T tν to directly obtain the boundary values [T tν ] ± j+1/2 in Eq. ( 35) from [T tν ] j .However, to enhance the numerical stability of our implementation, we instead choose Q = (m, T ti ) for the spatial reconstruction, where m is the mass density 1 An example of the analytic formula of ρ(J αβ (T tν )) can be found in Ref. [13].
of the fluid cell defined by More explicitly, we obtain [T tν ] ± j+1/2 as follows.First, a set of conserved densities [T tν ] j obtained in each intermediate numerical step can fail in reproducing the physical requirement m 2 ≥ 0. To avoid this situation, whenever a set of [T tν ] j is calculated, we replace where a small regulator s is chosen to be 10 −15 for the present study.Second, the MUSCL reconstruction would guarantee the non-oscillatory property of the reconstructed variables Q.We confirmed this property when Eqs. (37) and (38) are applied to Q = T tν .However, in this procedure, the non-oscillatory property of m 2 is not guaranteed, and thus the requirement m 2 > 0 can be broken.As a practical matter, we are faced with unphysical finite oscillations in ϵ and u µ obtained from T tν of this reconstruction.To suppress these oscillations, we instead apply the reconstruction by the choice Q = (m, T ti ) as in the following procedure: 1. Starting from [T tν ] j , construct the local mass density [m] j using Eq. ( 40).
We confirmed that this prescription significantly reduces the oscillating behavior of ϵ and u µ in our numerical tests.
The KT operators in the other directions, KT y and KT z , are defined similarly.The right-hand side of Eq. ( 32) combined with the MUSCL reconstruction (37)-( 39) references two cells away from the current cell in every direction.This property plays an important role in optimizing the implicit solver discussed below.
In the Milne coordinates, we solve hydrodynamics equations in the form of Eqs. ( 8) and ( 9) with T tν as the state variables.In this case, we use the same prescription as above with the following modifications: To obtain the energy density and velocity, Eqs. ( 3) and ( 4) are applied after transforming T τν to the local orthonormal frame, and the result is transformed back to the Milne basis.In the spatial reconstruction, we use m = τm in place of the mass density m.The geometric source term P(ϵ) in Eq. ( 8) is added to the right-hand side of Eq. (32).

Time discretization
The time evolution of the differential equations (32), which are obtained by the space discretization of Eq. ( 31), can be solved by RK methods introduced in Sec.3.1, where ⃗ y(t) and ⃗ h(⃗ y) in Eq. ( 29) correspond to the state variables {[T tν ] j } and the right-hand side of Eq. (32), respectively.In this study, we solve Eq. ( 32) using the implicit and explicit RK methods and demonstrate that the former can be more efficient for various cases.
We have tested various implicit and explicit RK methods with different stages and orders.They include the Gauss-Legendre family for the implicit methods, and Heun's and the Euler methods for the explicit.For the purpose of our study, we choose Heun's method as the representative of explicit methods, as it is commonly used in hydrodynamic codes for HIC [13,22].This method has two stages (S = 2) and second order of accuracy (p = 2), and its Butcher table is given in the left column of Table 1.For the implicit method we choose the one-stage Gauss-Legendre method [56,57,55], which has the same order of accuracy as Heun's method.The Butcher table of this method is given in the right column of Table 1.In Sec. 5, we compare the numerical results of these methods.In the following, we refer to Heun's and one-stage Gauss-Legendre methods as the explicit and implicit methods, respectively.
Another conventional explicit method at the same order of accuracy as Heun's method is the midpoint method, whose Butcher table is given in the middle column of Table 1.In Appendix A, we show that the performance of this method is almost the same as Heun's method in the Riemann problem.
Figure 2: Schematic picture of the iterative process K (0) , K (1) , • • • in the fixed-point method toward the solution K * for the single-variable case.

Fixed-point solver with local optimization
In implicit RK methods, ⃗ k (n) in Eq. ( 29) are obtained by solving the system of implicit equations.In this study we adopt the fixed-point method as a numerical solver for this procedure.To illustrate this method, we first rewrite Eq. ( 29) in a compact form with The fixed-point method solves Eq. ( 42) by starting from an initial guess ⃗ K (0) and iteratively updating the guess as A schematic illustration of this procedure for a single-variable case is shown in Fig. 2.
A necessary condition for the convergence of the fixed-point method is that the spectral radius of the Ja- at the solution ⃗ K = ⃗ K * satisfies R < 1.Moreover, the convergence is faster for smaller R. From Eqs. (43) and (44) one finds that the Jacobian matrix ∂ ⃗ F/∂ ⃗ K, and hence R, is proportional to ∆t.Therefore, the convergence is guaranteed for a sufficiently small ∆t and is faster for a smaller ∆t.However, the convergence is not guaranteed in general.In fact, we will see later that our method can be aborted for ∆t/∆x ≳ 0.4 due to a failure of the fixed-point solver.The convergence in these cases can be established with additional prescriptions as will be discussed in the forthcoming publications [58].
In the present study, we simply neglect these cases since our method works stably in the range of ∆t that is practically important.
The speed of convergence in the fixed-point method also depends on the choice of the initial guess ⃗ K (0) strongly.In this study, for the initial guess we use the solution of ⃗ K in the previous RK step.This choice comes from an expectation that ⃗ K hardly changes in a small time interval, and hence ⃗ K at the previous RK time step is a good guess for the next one.As we will see in Sec. 5, this prescription works quite well.For the first time step of RK at t = t init , we use ⃗ K (0) = (0, • • • , 0) for the initial guess.
To improve the efficiency, we utilize the locality of the hydrodynamic equations, i.e., the right-hand side of Eq. ( 32) depends only on cells up to two neighbors in each direction.By taking this advantage, in our implementation, we check the convergence cell by cell using the criterion with a pre-determined parameter e, where [ ⃗ K] j represents the set of variables in ⃗ K at cell j, and ⟨T tt ⟩ is the spatial average of T tt .Note that [ ⃗ K] j includes 4S variables corresponding to T tν for all stages.
Since ⟨T tt ⟩/∆t has the same unit as [ ⃗ K (l) ] j , e is dimensionless.This parameter should be properly tuned since smaller e leads to better accuracy but degrades the efficiency, and vice versa.In our numerical analysis in Sec. 5 we choose e = 10 −3 and e = 2 × 10 −4 for one and two space dimensions, respectively.
In Eq. ( 46), we require that the threshold scales as ∆t p−1 .Since the error of [ ⃗ K (l+1) ] j is parametrically of order ∆t higher than that of [ ⃗ K (l) ] j , the error of [ ⃗ K (l+1) ] j is of order ∆t p when the condition (46) is met.This order of [ ⃗ K (l+1) ] j is in accordance with the accuracy of RK in Eq. (30) 2 .In Sec. 5 we show that the expected scaling (30) is nicely obtained using the criterion (46).
More specifically, we update all the cells for the first iteration.For the succeeding iterations, we update cell j only when any of the following conditions are met: 1. Eq. ( 46) was not satisfied at cell j in the last update.2. Any of the cells adjacent to cell j was updated in the previous iteration and did not satisfy Eq. ( 46).
We update all the cells that need updating at once.We proceed with the iterations until Eq. ( 46) is reached for all the cells.We note that even once the condition ( 46) is satisfied at a cell, the condition can be violated later due to updates of neighboring cells.As we will show in the next section, this optimization drastically reduces the computational cost.In fact, only a single update is enough for most cells.We refer to this prescription as the local optimization in the following.
In the above procedure, we check only the most adjacent cells.Although the right-hand side of Eq. ( 32) also depends on the second neighboring cells, we have checked that their inclusion hardly changes the quality of numerical results despite an extra computational cost.We also notice that an update of a cell does not necessarily lead to the update of its neighboring cells in the next iteration.
Strictly speaking, the local optimization violates the energy-momentum conservation.This is because the flux at a boundary between two cells can be inconsistent when one cell is updated but the other is not.The numerical error from this violation, however, is controlled within the accuracy as we will show in Sec. 5.

Boundary conditions and vacuum
Since the numerical evolution of hydrodynamics is performed on a finite grid, one has to introduce boundary conditions.In our space scheme the evaluation of [KT i ] j references two neighboring cells.Therefore, the variables on these cells must be specified when they are outside the numerical grid.In this study, we employ the same boundary condition as that used in Ref. [17].In this boundary condition, when cells j ± 1 are outside the numerical grid, we use the variables on cell j for those of cell j ± 1.If cells j ± 2 are outside the numerical grid, the variables on cells j ± 1 are substituted there.Since this boundary condition roughly approximates an outflow without reflection, it is sufficient for an application to heavy-ion collisions when the grid is large enough that the boundary is essentially kept in the vacuum.
In the numerical analysis in Sec. 5, we deal with the vacuum region with ϵ(x j ) = 0, where x j represents the spatial coordinates of cell j.However, the treatment of ϵ(x j ) = 0 is problematic since it leads to an undefined velocity through Eq. ( 4).To avoid this situation, we force the energy-density ϵ(x j ) to be always at least ϵ vac = 10 −100 fm −4 in Eq. (3).

Dimensions
In this study, we compare the implicit and explicit RK methods in (1+1)-and (2+1)-dimensional numerical simulations.For the (1+1)-dimensional case, we deal with the Riemann problem discussed in Sec.2.2.1 with finite and zero minimum energy density ϵ min .We also perform (2+1)-dimensional simulations with the boost invariance for the Gubser flow and event-by-event initial conditions for heavy-ion collisions generated by T R ENTo [50].
In all cases, we consider a system of length 40 fm in each direction.To investigate the dependence on the mesh size ∆x, we perform numerical simulations with ∆x = 0.2 and 0.4 fm where the spatial extent is divided into 200 and 100 cells, respectively.The time step ∆t is fixed in a single simulation, and we perform simulations with different ∆t down to ∆t = 0.1×2 −5 ∆x to investigate the accuracy of the numerical results.

Error estimate
The results obtained by computational hydrodynamics have discretization errors.Numerical simulations are also affected by floating-point errors in general.In the present study, we evaluate these errors using two quantitative measures.First, we compare the energy density obtained by the numerical analysis, ϵ num (x j ), with that of the exact solution ϵ exact (x j ) the continuum case using the measure with Here, we take the maximum of |e 1 | and |e 2 | for the denominator of Eq. ( 48), instead of either one of them, to avoid a situation where Eq. ( 47) becomes extremely large when one of them is vanishingly small.As a result, the absolute value of Eq. ( 48) is bounded by 1 for e 1 > 0 and e 2 > 0. As for the second measure, we choose a numerical result obtained at small ∆t, denoted by ϵ ref (x j ), as a reference.We then compare ϵ num (x j ) with the reference as For the reference ϵ ref (x j ), we use the result by the explicit method with the smallest ∆t = 0.1 × 2 −5 ∆x.
To compare numerical errors of different time discretization methods, ∆ ref (x j ) is more convenient than ∆ exact (x j ).This is because a finite space-discretization error remains in ∆ exact (x j ) even in the ∆t → 0 limit, while one can discuss genuinely the time-discretization error using ∆ ref (x j ).
In Sec. 5, we inspect the magnitude of numerical error using the maximum and spatial average of |∆ ref (x j )|.However, we found that |∆ ref (x j )| tends to be close to unity near the vacuum region where ϵ exact (x j ) is vanishingly small.In this region ϵ num (x j ) and ϵ ref (x j ) have non-vanishing values due to the numerical diffusion, i.e., the unphysical propagation from the matter region, which is inevitable in computational hydrodynamics.Nevertheless, ϵ num (x j ) is still negligibly small and the large |∆ ref (x j )| should not be taken into account as an essential error.
To avoid this apparent error in the vacuum region, in our error analysis we exclude the cells where ϵ num (x j ) is smaller than a threshold value ϵ thr , where ϵ thr is determined so that the sum of the energy in the excluded cells is less than 10 −6 times the total energy.We also exclude cells within 1 fm from the boundaries to remove possible boundary effects.We then evaluate the maximum and the average of |∆ ref (x j )| in the remaining cells.

Estimate of computational cost
To compare the performance of different numerical methods, we also need a quantitative measure of the computational cost.However, the actual computational cost depends on the implementation, and its appropriate comparison is difficult in general.In this study, we focus on the fact that the most computationally demanding part of our analysis is the evaluation of the local KT operator [h KT ] j in Eq. (32).This means that the number of evaluations N KT serves as an approximate measure of the computational cost.We thus use the number of evaluations per cell, for the measure of the computational cost, where N cell stands for the number of total cells.For our explicit method with S = 2, [h KT ] j is evaluated twice per cell in a single time step.In the implicit method, the cost of a single time step depends on the number of iterations of the fixed-point solver, where one local iteration needs one evaluation of [h KT ] j with S = 1.Moreover, due to the local optimization, n KT is not necessarily an integer.
We have checked that n KT indeed reflects the real computational times well with appropriate optimizations.

1+1d Riemann problem with shock propagation
Now, let us discuss the numerical results.We begin in this subsection to deal with the Riemann problem, whose exact solution is presented in Sec.2.2.1.This is a (1+1)-dimensional problem in flat spacetime with initial conditions given by Eqs.(10) and (11).We numerically solve the hydrodynamic evolution from the initial time t init = 0 to t end = 15 fm.We set c 2 s = 1/3 and ϵ max = 10 fm −4 and perform our analysis for the two values of ϵ min = 1 and 0 fm −4 .The latter is important for the application to heavy-ion collisions where the matter expands into the vacuum, which typically causes numerical problems.

Numerical solutions
In the upper panel of Fig. 3, we first show the numerical result of ϵ(x j ) at t = t end for ϵ min = 1 fm −4 as a function of x/t.The red (blue) curves show the results obtained by the explicit (implicit) method with ∆t = 0.1 × 2 −5 ∆x while the dotted and dashed lines represent the results for ∆x = 0.4 and 0.2 fm, respectively.The black solid line shows the analytic solution ϵ exact (x j ).In the lower panel, we also show ∆ exact (x j ).
One sees that the results of the explicit and implicit methods degenerate well at this ∆t.This suggests that the error from the time discretization is well suppressed in both methods.By comparing the results of ∆x = 0.4 and 0.2 fm, one also finds that the numerical result approaches the analytic solution as ∆x becomes smaller.Nevertheless, a clear deviation is still observed at the smaller ∆x.This is significant around the rarefaction front and shock at x/t = −c s and v shock , respectively, which are shown in the enlarged subpanels of the upper panel.In particular, a large deviation ∆ exact ≃ 0.5 survives even for the small ∆x at x/t = v shock , which is a consequence of the discontinuity of ϵ exact (x j ) at this point.
In Fig. 4, the numerical result for ϵ min = 0 is shown in the same manner as Fig. 3.The figure shows that ∆ exact (x j ) becomes almost unity in the vacuum region x/t > 1.As explained in Sec.4.2, this result comes from vanishing ϵ exact (x j ) in the vacuum region.
In the following analysis, we use ∆ ref (x j ) defined in Eq. ( 49) in place of ∆ exact (x j ) to analyze the timediscretization error exclusively.

Error and computational cost
Next, let us compare the numerical error and the cost of the explicit and implicit methods.In Figs. 5 and 6, we plot the maximum and average of the error defined by Eq. ( 49) at t = t end for ϵ min = 1 and 0 fm −4 , respectively, for ∆x = 0.4 fm (upper plot) and 0.2 fm (lower plot).The horizontal axis is the computational cost n KT in Eq. (50).The red squares and blue circles show the explicit and implicit results, respectively.The large point in each line shows the result for ∆t/∆x = 0.1.The distance between the points along one line is obtained from varying ∆t by a factor of two.The numerical simulation fails when the time step ∆t is too large, so only successful results are shown in the figures.In the explicit method, the cost is n KT = S (t end − t init )/∆t = 30 fm/∆t.We observe several notable features in Figs. 5 and 6.First, all the results are approximately proportional to n −2 KT .This behavior is in accordance with the order of accuracy p = 2 of our explicit and implicit methods.Second, the computational cost of the implicit method is significantly smaller than the explicit method for a fixed |∆ ref | in both maximum and average numerical errors; The former is about 2-3 times as efficient as the latter.This means that with the same accuracy, the implicit method is faster than the explicit one.This is a remarkable outcome of the present study since it is widely considered that implicit methods are slower than explicit methods.The figures also show that this advantage of the implicit method is obtained through the reduction both in the numerical error and computational cost for a given ∆t.
To better understand the origin of these results, in Fig. 7 we show the number of iterations performed by the fixed-point solver in the implicit method on the t-x plane.The number of iterations is represented by different colors for the numerical simulations with ∆t/∆x = 0.1 and 0.1 × 2 −3 fm (left and right), and ∆x = 0.4 and 0.2 fm (top and bottom).The dashed and dash-dotted lines show the locations of the shock and the rarefaction front at which ϵ exact (x) and ∂ϵ exact (x)/∂x, respectively, have a discontinuity.The figure shows that for almost all cells one iteration is enough to achieve the convergence.This result is reasonable since the hydrodynamic variables would change slowly when they are smooth.As explained in Sec.3.3, we use the solution in the previous time step for the initial guess in the fixed-point solver, and it is already a good estimate for the next time step.Moreover, more than one iteration is required only around the shock and rarefaction, where the number of iterations is mostly two except for a few cells with three or four iterations.As a result, the average number of the [h KT ] j evaluations for an update is smaller than two in the implicit method.This makes the implicit method advantageous compared with the explicit one where two evaluations are required for every update.
In Figs. 7 and 8 one finds that the regions of two iterations have similar outlines in the left and right panels.This behavior would be attributed to fine structures that arise due to the space discretization.
As already discussed, in Figs. 5 and 6 we only show the successful results.We found that the explicit method fails for ∆t/∆x ≥ 1.6.This result is consistent with the Courant-Friedrichs-Lewy (CFL) condition [59,60] for stability, ∆t/∆x < C CFL where C CFL = O(1) depends on the detailed scheme.On the other hand, failures of the implicit method occur for ∆t/∆x ≥ 0.8.These failures happen when the fixed-point solver does not con-   verge.Despite A-stability of the implicit method itself, the fixed-point solver is not ensured to converge for large ∆t.

Gubser flow
Next, we perform the same comparison for the Gubser flow as discussed in Sec.2.2.2.We solve transverse (2+1)-dimensional hydrodynamics assuming the boost invariance.In this case, the solution contains neither discontinuity nor vacuum region.We employ ϵ 0 = 1 fm −4 and q = 1 for the parameters in Eq. ( 20) and set τ init = 1 fm for the initial proper time.We then solve the hydrodynamic equations numerically until the end time τ end = 10 fm.
In the top panels of Figs. 9 and 10, we first show the energy density at several values of τ for ∆x = 0.4 and 0.2 fm, respectively, obtained by the implicit method with ∆t = 0.1 × 2 −5 ∆x fm.Here we do not show the result of the explicit method since its behavior is almost the same as anticipated from the discussion in Sec.5.1.The panels show that the numerical solution has noncircular behaviors, while the exact solution, Eq. ( 20), has azimuthal symmetry.Correspondingly, ∆ exact (x j ) shown in the middle panels has oscillating behaviors.These errors are attributed to the space discretization, as we have checked that their structure hardly changes for small ∆t.By comparing Figs. 9 and 10, one also finds that the deviation from the exact solution is suppressed as ∆x becomes smaller.
Shown in the lower panels of Figs. 9 and 10 are the number of iterations in the fixed-point solver for each τ.These panels show that two iterations are needed around the region where ϵ(x j ) change rapidly, but only one iteration is enough for the other wide area.This result is consistent with the one in Sec.5.1.
In Fig. 11, we compare the maximum and the average of |∆ ref (x j )| in the explicit and implicit methods as functions of n KT .The meanings of the lines and symbols are the same as before.The figure shows that the implicit method is advantageous again in reducing the computational cost even in this case.

Event-by-event initial condition for heavy-ion collisions
Finally, we apply the same analysis to realistic cases of the relativistic heavy-ion collisions by generating the initial conditions using T R ENTo [50].We perform the (2+1)-dimensional event-by-event simulations for leadlead collisions with impact parameter b = 3 fm and collision energy √ s NN = 2.76 TeV.We use τ init = 0.48 fm as the initial time of the hydrodynamics.For the other parameters of T R ENTo, we take the same values as in Refs.[24,25].For the EOS, We use the parametrized form [20,61] of the Wuppertal-Budapest lattice EOS [62].We then solve the hydrodynamic evolution until τ end = 10 fm.
In the upper panels of Figs. 12 and 13, we show the transverse profile of energy density ϵ(x j ) for one event in the x-y plane obtained with the finest ∆t = 0.1 × 2 −5 ∆x in the implicit method for ∆x = 0.4 and 0.2 fm, respectively.The cells satisfying ϵ(x j ) < ϵ thr are shown by the white color.As discussed in Sec.4.2, these cells are discarded in the error analysis.In the lower panels, the corresponding number of iterations in the fixed-point solver is shown.These results again show that the average number of iterations is smaller than two, which is a reasonable result in light of the arguments in the previous subsections.
In Fig. 14, we plot the maximum and the average of |∆ ref | as functions of n KT for ten initial conditions generated randomly by T R ENTo.The figure shows that the implicit method tends to carry out the analysis with a smaller computational cost than the explicit method for a given error, although its advantage is less significant compared with the results in previous sections.

Conclusion
In this study we proposed a new method that numerically solves relativistic hydrodynamics with an implicit Runge-Kutta (RK) time integrator and implemented an efficient solver for ideal hydrodynamics.The one-stage Gauss-Legendre method is employed as the implicit method while the Kurgamov-Tadmor (KT) scheme is used for the space discretization.To solve the implicit RK equations, we propose a simple fixed-point solver with local optimization.For the initial guess of the solver, we employ the solution in the previous RK step.
The performance of our method was compared with Heun's method, which is one of the conventional explicit methods, with the same space discretization for different initial conditions in one and two space dimensions.We found that the computational cost of the implicit method is always smaller than the explicit one to attain the same accuracy except for impractically large ∆t/∆x.The advantage of the former mainly comes from the fact that the fixed-point solver requires less than two iterations on average for convergence.In comparison, Heun's method as a two-stage RK method needs two substeps in each RK update.This makes the number of evaluations of the KT operator smaller in the implicit case.These results show that our implicit method can be efficient and useful for the application to HIC.This also suggests that the implicit RK method can be used for more general setups of computational hydrodynamics to attain efficient and simple implementations.Because implicit methods generally have an advantage in their numerical stability, they should lead to better implementation of the relativistic hydrodynamics with less technical procedures for stabilization, such as the limiters and regulations of variables.This advantage will be more effective in dealing with more complex equations such as fluctuating hydrodynamics and magnetohydrodynamics.
Although we limited our analysis to ideal hydrodynamics throughout this paper, it is important to extend the analysis to (3+1)-dimensional viscous hydrodynamics [58] and to incorporate conserved charges for practical applications to HIC.Our implicit method has a problem in the convergence of the fixed-point method at large ∆t/∆x.A solution to improve the convergence will be discussed in the forthcoming publication [58].

Figure 1 :
Figure 1: Analytical solution to the Riemann problem for t > 0 with x = x/t.

Figure 3 :
Figure 3: Energy density ϵ(x j ) (upper panel) and its error compared with the exact solution ∆ exact (x j ) (lower panel) for the Riemann problem at t end = 15 fm with ϵ min = 1 fm −4 and ∆t = 0.1 × 2 −5 ∆x.

Figure 5 :
Figure 5: Maximum and average of the numerical error |∆ ref | defined in Eq. (49) as functions of the computational cost n KT for the Riemann problem with ϵ min = 1 fm −4 .The maximum and average are shown by the solid and dashed lines, and the explicit and implicit results are shown by red square and blue circle symbols, respectively.The upper and lower panels show the results for ∆x = 0.4 and 0.2 fm, respectively.The large points show the results at ∆t/∆x = 0.1, while the other small points are obtained by varying ∆t by a factor of 2.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: Energy density ϵ(x j ) (upper panels), error ∆ exact (x j ) (middle panels) and the number of iterations of the fixed-point solver in each cell (lower panels) for the Gubser flow at several proper times τ for ∆t = 0.1 × 2 −5 ∆x and ∆x = 0.4 fm.

Table 1 :
Properties and the Butcher tables of Heun's, the midpoint, and GL1 methods.cell in the x direction from j, and [T tν ] ± j+1/2 are values reconstructed on the cell boundary that will be specified below.The function T xν ([T tν ] ± j+1/2 ) calculates the flux as a function of the reconstructed [T tν ] ± j+1/2 using Eqs.(