# Efficient solver of relativistic hydrodynamics with implicit Runge-Kutta method

Nathan Touroux, Masakiyo Kitazawa, Koichi Murase, Marlene Nahrgang

► To cite this version:

Nathan Touroux, Masakiyo Kitazawa, Koichi Murase, Marlene Nahrgang. Efficient solver of relativistic hydrodynamics with implicit Runge-Kutta method. Progress of Theoretical and Experimental Physics [PTEP], 2024, 2024 (6), pp.063D01. 10.1093/ptep/ptae058. hal-04152119

HAL Id: hal-04152119

https://hal.science/hal-04152119v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Solver of Relativistic Hydrodynamics with an Implicit Runge–Kutta Method

Nathan Touroux ID<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, Masakiyo Kitazawa ID<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Koichi Murase ID<fnref n="3" />, and Marlene Nahrgang ID<fnref n="1" />

Received February 2, 2024; Revised April 5, 2024; Accepted April 11, 2024; Published April 23, 2024

We propose a new method to solve the relativistic hydrodynamic equations based on implicit Runge–Kutta methods with a locally optimized fixed-point iterative solver. For numerical demonstration, we implement our idea for ideal hydrodynamics using the one-stage Gauss–Legendre method as an implicit method. The accuracy and computational cost of our new method are compared with those of explicit ones for the (1+1)D Riemann problem, as well as the (2+1)D Gubser flow and event-by-event initial conditions for heavy-ion collisions generated by TRENT0. 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, while it may not converge with an unrealistically large $\Delta t$. By showing a relationship between the one-stage Gauss–Legendre method with the iterative solver and the two-step Adams–Bashforth method, we argue that our method benefits from both the stability of the former and the efficiency of the latter.

Subject Index D31

## 1. 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–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 a simulation of heavy-ion collision reactions based on the dynamical models and compare the results with high-statistics experimental data [4–10]. In modern modeling of heavy-ion collisions, a relativistic viscous hydrodynamic code [11–23] forms the core of the dynamical description of the reactions together with an event-by-event initial-state model and a final-state hadronic cascade.

<footnote n="1">¹SUBATECH UMR 6457 (IMT Atlantique, Université de Nantes, IN2P3/CNRS), 4 rue Alfred Kastler, 44307 Nantes, France</footnote>

<footnote n="2">²Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan</footnote>

<footnote n="3">³Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="4">⁴J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies, KEK, Tokai, Ibaraki 319-1106, Japan</footnote>

<footnote>*E-mail: nathan.touroux@subatech.in2p3.fr</footnote>
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–29] and higher-order fluctuation observables, such as the skewness and the kurtosis of the event-by-event distributions of conserved charges [30–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 necessary 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) [16, 17, 19, 37, 38], and other Gudunov-type schemes [18, 21]. In this study, we instead explore the time discretization. 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 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 could be one of the reasons that they have not been applied in the hydrodynamics for HIC. However, this 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, $\Delta t$, for linear ordinary differential equations. One can, therefore, hope that stable calculations can be performed using an implicit method at a larger $\Delta t$ even for nonlinear 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 nontrivial 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 [19, 39–43], nonequilibrium 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 $\Delta t$ due to high-frequency modes. Genuine implicit methods can stably solve such stiff equation systems with a moderate time step $\Delta 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 in e.g. 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 TRENTO initial conditions for heavy-ion 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 $\Delta t$ for the cases of these initial conditions. We also find that the fixed-point iteration may fail to converge with large $\Delta t$, where the Courant-Friedrichs-Lewy (CFL) condition [51,52] is unsatisfied. Nevertheless, as long as a solution is obtained, its stability is ensured.

We argue that the better performance of our method comes from the fact that it can be interpreted as a combination of the one-stage Gauss-Legendre (GL1) and the two-step Adams-Bashforth (AB2) methods. We show that our method with a single iteration can be related to AB2, which has the accuracy order $p = 2$ but less stability, while it converges to GL1 after iterations. Multiple iterations are needed only when AB2 diverges from GL1. In this sense, our method makes use of GL1 and AB2 depending on the local stiffness, thereby exploiting the advantages in stability and efficiency of both methods in an optimal way. Because of this property, our method is expected to be advantageous even when applied to more complex systems of equations such as viscous hydrodynamics.

This paper is organized as follows. In Section 2, we review the hydrodynamic equations and their analytical solutions that will be used in later sections. In Section 3, we next review the implicit and explicit RK methods and the spatial discretization scheme in the present study with technical details. In Section 4, we propose new improvements to implicit RK methods and discuss the accuracy, efficiency, and stability of the improved implicit method. The numerical setup and analyses are given in Section 5. We then discuss the numerical results in Section 6. The final section is devoted to discussions and a summary.

## 2. Relativistic hydrodynamics

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

### 2.1. Hydrodynamics equations

Hydrodynamic equations are given by the conservation law. We focus on the conservation of the energy-momentum tensor $T^{\mu\nu}$:

$$ \begin{aligned} \partial_{,\mu} T^{\mu\nu} &= \partial_{\mu} T^{\mu\nu} + \Gamma^{\mu}_{\alpha\mu} T^{\alpha\nu} + \Gamma^{\nu}_{\alpha\mu} T^{\mu\alpha} \\ &= 0, \end{aligned} \quad (1) $$
where $\partial_\mu$ denotes the covariant derivative, and $\Gamma^{\nu}_{\mu\alpha}$ the associated Christoffel symbols. For the general case including further conserved currents, the scheme described in this paper can be naturally extended.

In this paper we consider the ideal hydrodynamics where the energy–momentum tensor is written in terms of thermodynamic quantities as

$$T^{\mu\nu} = \epsilon u^{\mu} u^{\nu} - P \Delta^{\mu\nu}, \qquad (2)$$

where $u^\mu$ is the four-velocity with $u^\mu u_\mu = 1$, $\epsilon$ the energy density, $P$ the pressure, and $\Delta_{\mu\nu} = g_{\mu\nu} - u_\mu u_\nu$ the projector orthogonal to $u^\mu$. Equation (2) contains five undetermined variables ($u^i$, $\epsilon$, $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(\epsilon)$. Equations (1) and (2) together with the EOS define the ideal hydrodynamics.

In the Cartesian coordinates, we choose the *conserved densities*,<fnref n="1" /> $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\mu}$ as a function of $T^{i\mu}$ that satisfies Eq. (2). We first obtain $\epsilon$ and $u^\mu$ from $T^{\nu\mu}$ and then calculate $T^{\nu\mu}$ using Eq. (2). To numerically determine $\epsilon$ and $u^\mu$ satisfying Eq. (2), it is convenient to solve the following equations [17] for the energy density $\epsilon$ ($> 0$) and the velocity $v$ ($\in [0, 1]$) iteratively:

$$\epsilon = T^{tt} - vK, \qquad v = \frac{K}{T^{tt} + P(\epsilon)}, \qquad (3)$$

where $K = \sqrt{\sum_{i=1}^3 (T^{ti})^2}$. These equations require the condition $T^{ti} \ge K$, which is ensured by the prescription described in Section 3.2. Then, $u^\mu$ is given by

$$u^i = \frac{T^{ti}}{\sqrt{(T^{tt} + P(\epsilon))(\epsilon + P(\epsilon))}}, \qquad u^i = \sqrt{1 + (u^i)^2}. \qquad (4)$$

In high-energy heavy-ion collisions, the space-time evolution of the produced medium is well described by the Bjorken assumption of longitudinal boost invariance [53]. In this case, it is convenient to employ the Milne coordinates

$$\tau = \sqrt{t^2 - z^2}, \qquad \eta = \tanh^{-1} \frac{z}{t}, \qquad (5)$$

or, written inversely,

$$t = \tau \cosh \eta, \qquad z = \tau \sinh \eta, \qquad (6)$$

where $z$ is the longitudinal direction parallel to the beam axis while the transverse coordinates, $x$ and $y$, remain unchanged. In this coordinate system, the metric is $g_{\mu\nu} = \text{diag}(1, -1, -1, -\tau^2)$, and the nonzero Christoffel symbols are $\Gamma^\eta_{\tau\eta} = \Gamma^\eta_{\eta\tau} = 1/\tau$ and $\Gamma^\tau_{\eta\eta} = \tau$.

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 the *boost invariance* in the Milne coordinate system where $\partial_\eta = 0,^2 u^\eta = 0,^3$ and $T^\tau\eta = 0$. In this case, the conservation equation (1) is given

<footnote n="1">¹In this paper, for convenience, we mean by the *conserved densities* the densities of conserved charges even though the densities in each cell are not conserved separately.</footnote>

<footnote n="2">²Here, when we apply $\partial_\eta$ to tensors, we assume that the local basis of the tangent space is chosen to be the natural basis of the Milne coordinate system. The relation $\partial_\eta = 0$ does not hold in mixed basis cases such as $\partial_\eta u^i$ and $\partial_\eta u^z$.</footnote>

<footnote n="3">³More specifically, the flow velocity can be written in the Milne coordinates as $(u^\tau, u^x, u^y, u^\eta) = (u^\tau(\tau, x, y), u^x(\tau, x, y), u^y(\tau, x, y), 0)$. This corresponds to $(u^i, u^x, u^y, u^z) = (u^\tau(\tau, x, y)\cosh \eta, u^x(\tau, x, y), u^y(\tau, x, y), u^z(\tau, x, y)\sinh \eta)$ in the Cartesian coordinate basis.</footnote>
by

$$ \partial_{\tau} T^{\tau\nu} = -\partial_{x} T^{x\nu} - \partial_{y} T^{y\nu} - \frac{[\epsilon + P(\epsilon)]u^{\tau}u^{\nu}}{\tau}, \quad (\nu = \tau, x, y). $$

It is practically useful to rewrite Eq. (7) as

$$ \partial_{\tau} \tilde{T}^{\tau\tau} = -\partial_{x} \tilde{T}^{x\tau} - \partial_{y} \tilde{T}^{y\tau} - P(\epsilon), \quad (7) $$

$$ \partial_{\tau} \tilde{T}^{\tau l} = -\partial_{x} \tilde{T}^{x l} - \partial_{y} \tilde{T}^{y l}, \quad (l = x, y), \quad (8) $$

with $\tilde{T}^{\mu\nu} = \tau T^{\mu\nu}$ and to treat $\tilde{T}^{\mu\nu}$ as the state variables.

## 2.2. 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.

### 2.2.1. 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,

$$ \epsilon(t = 0, x) = \begin{cases} \epsilon_{\max}, & x < 0, \\ \epsilon_{\min}, & x \ge 0, \end{cases} \quad (9) $$

$$ u^l(t = 0, x) = 1, \quad u^x(t = 0, x) = 0, \quad (10) $$

with a discontinuity at $x = 0$. The analytical solution is given for a particular choice of an EOS, $P(\epsilon) = c_s^2\epsilon$. Here, $c_s$ is the speed of sound. For $t > 0$, the solution $\epsilon(x, t) = \tilde{\epsilon}(\tilde{x})$ with $\tilde{x} = x/t$ is given piecewise [48]:

$$ \tilde{\epsilon}(\tilde{x}) = \begin{cases} \epsilon_{\max}, & \tilde{x} < -c_s, \\ \epsilon_{\max} \left( \frac{(1-c_s)(1-\tilde{x})}{(1+c_s)(1+\tilde{x})} \right)^{\frac{1+c_s^2}{c_s}}, & -c_s \le \tilde{x} < -v_{\text{r-p}}, \\ \epsilon_{\text{plateau}}, & -v_{\text{r-p}} \le \tilde{x} < v_{\text{shock}}, \\ \epsilon_{\min}, & v_{\text{shock}} < \tilde{x}, \end{cases} \quad (11) $$

$$ u^x(\tilde{x}) = \frac{v(\tilde{\epsilon}(\tilde{x}))}{\sqrt{1 - v(\tilde{\epsilon}(\tilde{x}))^2}}, \quad u^l(\tilde{x}) = \sqrt{1 + (u^x(\tilde{x}))^2}, \quad (12) $$

where

$$ v_{\text{r-p}} = \frac{v(\epsilon_{\text{plateau}}) - c_s}{1 - v(\epsilon_{\text{plateau}})c_s}, \quad (13) $$

$$ v(\epsilon) = \frac{1 - (\epsilon / \epsilon_{\max})^{2c_s/(1+c_s^2)}}{1 + (\epsilon / \epsilon_{\max})^{2c_s/(1+c_s^2)}}, \quad (14) $$

The constants $\epsilon_{\text{plateau}}$ and $v_{\text{shock}} = -v_r$ are obtained by solving the following equations for $(\epsilon_{\text{plateau}}, v_l, v_r)$:

$$ \epsilon_{\text{plateau}}\gamma(v_l)v_l = \epsilon_{\min}\gamma(v_r)v_r, \quad (15) $$

$$ 4\epsilon_{\text{plateau}}\gamma(v_l)v_l^2 + \epsilon_{\text{plateau}} = \epsilon_{\min}\gamma(v_r)v_r^2 + \epsilon_{\min}, \quad (16) $$

$$ v_l = \frac{v(\epsilon_{\text{plateau}}) + v_r}{1 + v(\epsilon_{\text{plateau}})v_r}, \quad (17) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Analytical solution to the Riemann problem for $t > 0$ with $\tilde{x} = x/t$.</figcaption></figure>

where $\gamma(v) = 1/\sqrt{1-v^2}$. We note that Eqs. (15) and (16) represent the Taub equations [54,55]. The solution $\tilde{\epsilon}(\tilde{x})$ is illustrated in Fig. 1. The discontinuity at $\tilde{x} = v_{\text{shock}}$ is called the shock front, and the curve in the range $-c_s \le \tilde{x} < -v_{\text{r-p}}$ is called the rarefaction wave.

In heavy-ion collisions, the created matter is surrounded by the vacuum, where $\epsilon(t, x) = 0$, and the medium expands into it. Such a situation can be mimicked in the Riemann problem by setting $\epsilon_{\min} = 0$, where the solution is simplified to

$$ \tilde{\epsilon}(\tilde{x}) = \begin{cases} \epsilon_{\max}, & \tilde{x} < -c_s, \\ \epsilon_{\max} \left( \frac{(1-c_s)(1-x)}{(1+c_s)(1+x)} \right)^{\frac{1+c_s^2}{2c_s}}, & -c_s \le \tilde{x} < 1, \\ 0, & 1 \le \tilde{x}. \end{cases} \qquad (18) $$

### 2.2.2. Gubser flow.
Another exact solution is the Gubser flow [49], which is a solution for the conformal EOS, $P(\epsilon) = \epsilon/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

$$ \epsilon(\tau, r) = \frac{\epsilon_0(2q)^{8/3}}{\tau^{4/3} (1 + 2q^2(\tau^2 + r^2) + q^4(\tau^2 - r^2)^2)^{4/3}}, \qquad (19) $$

$$ u^{\tau}(\tau, r) = \cosh k(\tau, r), \qquad (20) $$

$$ u^{x}(\tau, r) = \frac{x}{r} \sinh k(\tau, r), \qquad (21) $$

$$ u^{y}(\tau, r) = \frac{y}{r} \sinh k(\tau, r), \qquad (22) $$

$$ u^{\eta}(\tau, r) = 0, \qquad (23) $$

with

$$ k(\tau, r) = \tanh^{-1} \frac{2q^2 \tau r}{1 + q^2 \tau^2 + q^2 r^2}, \qquad (24) $$

$$ r = \sqrt{x^2 + y^2}, \qquad (25) $$

where $q$ and $\epsilon_0$ are free parameters.

## 3. 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
<table><caption>**Table 1.** Properties and the Butcher tableaux of Heun's, the midpoint, and the GL1 methods.</caption><thead><tr><th></th><th colspan="2">Heun</th><th>Midpoint</th><th colspan="2">Gauss-Legendre 1 (GL1)</th></tr><tr><th>Type</th><th colspan="2">Explicit</th><th>Explicit</th><th colspan="2">Implicit</th></tr></thead><tbody><tr><th>Stage <i>S</i></th><td colspan="2">2</td><td>2</td><td colspan="2">1</td></tr><tr><th>Order <i>p</i></th><td colspan="2">2</td><td>2</td><td colspan="2">2</td></tr><tr><th>c<sub>n</sub></th><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>a<sub>nm</sub></th><td>1</td><td>1</td><td>0</td><td>0.5</td><td>0.5</td></tr><tr><th>b<sub>m</sub></th><td colspan="2">0.5</td><td>0.5</td><td>0</td><td>1</td></tr><tr><th></th><td colspan="2"></td><td></td><td>0.5</td><td>0.5</td></tr><tr><th></th><td colspan="2"></td><td></td><td colspan="2">1</td></tr></tbody></table>

with a focus on the differences between implicit and explicit methods in Section 3.1. The space discretization scheme and the boundary condition applied in our study are then discussed in Sections 3.2 and 3.3.

### 3.1. Implicit 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

$$ \partial_t \vec{y}(t) = \vec{h}(t, \vec{y}(t)), \qquad (26) $$

where $\vec{y}(t)$ represents the full set of state variables. The RK method solves the variables $\vec{y}(t_{\text{end}})$ at the end time $t = t_{\text{end}}$ given an initial condition $\vec{y}(t_{\text{init}})$ at the initial time $t = t_{\text{init}}$ by discretizing the time into steps of length $\Delta t$ and iteratively updating $\vec{y}(t)$ for every time step. For one time step from $t$ to $t + \Delta t$, the state variable of the next time $\vec{y}(t + \Delta t)$ is constructed from $\vec{y}(t)$ as

$$ \vec{y}(t + \Delta t) = \vec{y}(t) + \Delta t \sum_{n=1}^{S} b_n \vec{k}_{(n)}, \qquad (27) $$

$$ \vec{k}_{(n)} = \vec{h} \left( t + c_n \Delta t, \vec{y}(t) + \Delta t \sum_{m=1}^{S} a_{nm} \vec{k}_{(m)} \right), \qquad (28) $$

where $S$ is the number of stages of the RK method, and the coefficients $a_{nm}$, $b_n$, $c_n$ for $n, m = 1, \dots, S$ are determined so that the method has the desired properties. One specific set of the coefficients defines a particular RK method. The table that summarizes these coefficients is called a Butcher tableau, examples of which are shown in Table 1.

In a RK method, the numerical error of $\vec{y}(t_{\text{end}})$ from the exact solution $\vec{y}^*(t_{\text{end}})$ is suppressed in the limit $\Delta t \to 0$ as

$$ \|\vec{y}(t_{\text{end}}) - \vec{y}^*(t_{\text{end}})\| \le C(t_{\text{end}}) \Delta t^p, \qquad (29) $$

where $C(t)$ is specific to the ODE and is independent of $\Delta t$. The largest integer $p$ satisfying Eq. (29) 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 \ge n$. In this case, $\vec{k}_{(n)}$ in Eq. (28) depend only on $\vec{y}(t)$ and $\vec{k}_{(m)}$ for $m < n$, so $\vec{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 (28) for $\vec{k}_{(n)}(t)$, which is typically nonlinear 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; they are always linearly stable (A-stability) as long as the exact solution is not divergent [56]. This means that a larger $\Delta t$ can be adopted for stiff equations in implicit methods while keeping the 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 [57].

## 3.2. 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 $\Delta x$ along each coordinate. In the KT scheme, one chooses the set of the conserved densities of the cells, $\{[T^{tv}]_j\}$, as the state variables, where the subscript $j$ of the square bracket specifies a cell. We spatially discretize the continuity equation

$$ \partial_t T^{tv} = -\partial_x T^{xv} - \partial_y T^{yy} - \partial_z T^{zv}, \quad (v = t, x, y, z), \quad (30) $$

as

$$ \partial_t [T^{tv}]_j = [h_{KT}(T^{tv})]_j, \quad (31) $$

with

$$ [h_{KT}(T^{tv})]_j = -[KT_x(T^{tv})]_j - [KT_y(T^{tv})]_j - [KT_z(T^{tv})]_j, \quad (32) $$

using the KT operators $KT_i$. In our case, the KT operator in the $x$ direction, $KT_x$, is written as

$$ [KT_x(T^{tv})]_j = \frac{H_{j+1/2} - H_{j-1/2}}{\Delta x}, \quad (33) $$

$$ H_{j+1/2} = \frac{T^{xv}([T^{tv}]_{j+1/2}^{+}) + T^{xv}([T^{tv}]_{j+1/2}^{-})}{2} + \frac{a_{j+1/2}}{2} ([T^{tv}]_{j+1/2}^{+} - [T^{tv}]_{j+1/2}^{-}), \quad (34) $$

where $H_{j+1/2}$ is the flux at the boundary, $H_{j-1/2} = H_{(j-1)/2+1/2}$ with $j \pm n$ representing the $\pm n$th neighboring cell in the $x$ direction from $j$, and $[T^{tv}]_{j+1/2}^{\pm}$ are values reconstructed on the cell boundary that will be specified below. The function $T^{xv}([T^{tv}]_{j+1/2}^{\pm})$ calculates the flux as a function of the reconstructed $[T^{tv}]_{j+1/2}^{\pm}$ using Eqs. (3) and (4). In this study, we use the Newton
method to solve Eq. (3) for v. The maximum local speed of propagation is given by

$$a_{j+1/2} = \max \left\{ \rho(J_{\alpha\beta}([T^{tv}]_{j+1/2}^{+})), \rho(J_{\alpha\beta}([T^{tv}]_{j+1/2}^{-})) \right\}, \quad (35)$$

where $J_{\alpha\beta}(T^{tv}) = \partial T^{x\alpha}/\partial T^{t\beta}|_{T^{tv}}$ is a 4 × 4 Jacobian matrix, and $\rho(M)$ is the spectral radius of a matrix $M$.<fnref n="4" />

In the original study of the KT scheme [47], the quantities on the cell boundary $Q_{j+1/2}^{\pm}$ are calculated from the cell values $Q_j$ using the MUSCL reconstruction as

$$Q_{j+1/2}^{\pm} = Q_{j+1/2 \pm 1/2} \mp \frac{\Delta x}{2} [\hat{\partial}_x Q]_{j+1/2 \pm 1/2}, \quad (36)$$

$$[\hat{\partial}_x Q]_j = \minmod\left(\minmod\left(\theta \frac{Q_j - Q_{j-1}}{\Delta x}, \frac{Q_{j+1} - Q_{j-1}}{2\Delta x}\right), \theta \frac{Q_{j+1} - Q_j}{\Delta x}\right), \quad (37)$$

$$\minmod(a, b) = \frac{\operatorname{sgn}(a) + \operatorname{sgn}(b)}{2} \min(|a|, |b|). \quad (38)$$

The flux limiter in Eq. (37) guarantees the nonoscillatory property of $Q$ with the parameter $\theta$ restricted to $1 \le \theta \le 2$. Choosing $\theta < 1$ would produce excessive numerical diffusion, whereas choosing $\theta > 2$ would produce unphysical oscillation. We choose $\theta = 1.1$ following Ref. [13], which gives nonoscillatory results with a limited numerical diffusion. We verified that changing the value of $\theta$ hardly changes the tendency of our results for accuracy and efficiency in Section 6 unless $\theta$ is too close to two.

In relativistic hydrodynamics, the reconstruction could be simply applied by choosing $Q = T^{tv}$ to directly obtain the boundary values $[T^{tv}]_{j+1/2}^{\pm}$ in Eq. (34) from $[T^{tv}]_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 of the fluid cell defined by

$$m^2 = g_{\mu\nu} T^{\mu\nu} T^{\nu\mu} = g_{tt} (T^{tt})^2 - K^2. \quad (39)$$

More explicitly, we obtain $[T^{tv}]_{j+1/2}^{\pm}$ as follows.

First, a set of conserved densities $[T^{tv}]_j$ obtained in each intermediate numerical step can fail in reproducing the physical requirement $m^2 \ge 0$. To avoid this situation, whenever a set of $[T^{tv}]_j$ is calculated, we replace $T^{ti}$ with

$$\max(T^{tt}, (1+s)K), \quad (40)$$

where a small regulator $s$ is chosen to be $10^{-15}$ for the present study.

Second, the MUSCL reconstruction would guarantee the nonoscillatory property of the reconstructed variables $Q$. We confirmed this property when Eqs. (36) and (37) are applied to $Q = T^{tv}$. However, in this procedure, the nonoscillatory 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 $\epsilon$ and $u^{\mu}$ obtained from $T^{tv}$ 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^{tv}]_j$, construct the local mass density $[m]_j$ using Eq. (39).

(2) Reconstruct $[m]_{j+1/2}^{\pm}$ and $[T^{ti}]_{j+1/2}^{\pm}$ from $[m]_j$ and $[T^{ti}]_j$, respectively, using Eq. (36).

(3) Determine $[T^{ti}]_{j+1/2}^{\pm}$ from $[m]_{j+1/2}^{\pm}$ and $[T^{ti}]_{j+1/2}^{\pm}$ by inverting Eq. (39).

We confirmed that this prescription significantly reduces the oscillating behavior of $\epsilon$ and $u^{\mu}$ in our numerical tests.

<footnote n="4">⁴An example of the analytic formula of $\rho(J_{\alpha\beta}(T^{tv}))$ can be found in Ref. [13].</footnote>
The KT operators in the other directions, KT<sub>y</sub> and KT<sub>z</sub>, are defined similarly. The right-hand side of Eq. (31) combined with the MUSCL reconstruction (36)–(38) 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. (7) and (8)
with $\tilde{T}^{\tau\nu}$ 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^{\tau\nu}$ to the local orthonormal frame, and the result is transformed back to
the Milne basis. In the spatial reconstruction, we use $\tilde{m} = \tau m$ in place of the mass density $m$.
The geometric source term $P(\epsilon)$ in Eq. (7) is added to the right-hand side of Eq. (31).

3.3. Boundary conditions and vacuum

Since the numerical evolution of hydrodynamics is performed on a finite grid, one has to intro-
duce boundary conditions. In our space scheme the evaluation of [KT<sub>j</sub>] references two neigh-
boring 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 out-
flow 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 Section 6, we deal with the vacuum region with $\epsilon(\mathbf{x}_j) = 0$, where $\mathbf{x}_j$ represents the spatial coordinates of cell $j$. However, the treatment of $\epsilon(\mathbf{x}_j) = 0$ is problematic since it leads to an undefined velocity through Eq. (4). To avoid this situation, we force the energy density $\epsilon(\mathbf{x}_j)$ to be always at least $\epsilon_{\text{vac}} = 10^{-100} \text{ fm}^{-4}$ in Eq. (3).

4. Locally optimized fixed-point implicit solver

The time evolution of the differential equations (31), which are obtained by the space discretiza-
tion of Eq. (30), can be solved by RK methods introduced in Section 3.1, where $\vec{y}(t)$ and $\vec{h}(\vec{y})$
in Eq. (28) correspond to the state variables $\{[T^{\tau\nu}]_j\}$ and the right-hand side of Eq. (31), re-
spectively. For implicit RK methods, we need to numerically solve implicit equations for sub-
stantial problems. However, implicit methods usually cannot be efficient without modifications
compared to explicit ones.

To overcome this excessive computational cost of implicit methods, we propose an improved
solver for implicit RK methods in this study. After defining our method in Section 4.1, we
discuss the convergence condition of our method and its connection to the stiffness and stability
in Section 4.2. In Section 4.3, we finally discuss the accuracy, efficiency, and stability for a
specific case of the implicit RK method called the one-stage Gauss-Legendre (GL1) method
by relating it to the two-step Adams-Bashforth (AB2) method.

4.1. *Fixed-point solver with local optimization*

In implicit RK methods, $\vec{k}_{(n)}$ in Eq. (28) 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
<figure><img src="image_3.png" /><figcaption>Fig. 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.</figcaption></figure>

illustrate this method, we first rewrite Eq. (28) in a compact form,

$$ \vec{K} = \vec{F}(\vec{K}), \qquad (41) $$

with $\vec{K} = (\vec{k}_{(1)}, \dots, \vec{k}_{(S)})$ and

$$ \vec{F}(\vec{K}) = (\vec{f}_{(1)}(\vec{K}), \dots, \vec{f}_{(S)}(\vec{K})), \qquad (42) $$

$$ \vec{f}_{(n)}(\vec{K}) = \vec{h} \left( t + c_n \Delta t, \vec{y}(t) + \Delta t \sum_m a_{nm} \vec{k}_{(m)} \right). \qquad (43) $$

The fixed-point method solves Eq. (41) by starting from an initial guess $\vec{K}^{(0)}$ and iteratively updating the guess as $\vec{K}^{(l+1)} = \vec{F}(\vec{K}^{(l)})$ for $l = 0, 1, 2, \dots$ until the convergence is reached. A schematic illustration of this procedure for a single-variable case is shown in Fig. 2.

The convergence in the fixed-point method strongly depends on the choice of the initial guess $\vec{K}^{(0)}$. The initial guess should not be too far from the true solution. In this study, for the initial guess, we use the solution of $\vec{K}$ in the previous RK step. This choice comes from an expectation that $\vec{K}$ only changes at the order of $O(\Delta t)$, and hence $\vec{K}$ at the previous RK time step is a good guess for the next one. As we will see below, this choice of $\vec{K}^{(0)}$ plays a crucial role in realizing good efficiency in our method. For the first time step of RK at $t = t_{init}$, we use $\vec{K}^{(0)} = (0, \dots, 0)$ for the initial guess.

To improve the efficiency, we utilize the locality of discretized partial differential equations; i.e. a finite number of cells is referenced in an evaluation of the time derivative for one cell. In our case, the right-hand side of Eq. (31) 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

$$ \|[\vec{K}^{(l+1)}]_j - [\vec{K}^{(l)}]_j\| = \|[\vec{F}(\vec{K}^{(l)})]_j - [\vec{K}^{(l)}]_j\| < e^{\frac{\langle T^{tt} \rangle}{\Delta t}} \left(\frac{\Delta t}{\Delta x}\right)^p, \qquad (44) $$
with a predetermined parameter *e*, where $[\vec{K}]_j$ represents the set of variables in $\vec{K}$ at cell *j*, and $\langle T^{il} \rangle$ is the spatial average of $T^{il}$ in the entire numerical grid. Note that $[\vec{K}]_j$ includes 4S variables corresponding to $T^{iv}$ for all stages.

Since $\langle T^{il} \rangle / \Delta t$ has the same unit as $[\vec{K}^{(l)}]_j$, *e* is dimensionless. This parameter should be properly tuned for each calculation setup since smaller *e* leads to better accuracy but degrades the efficiency, and vice versa. In our numerical analysis in Section 6 we choose $e = 10^{-3}$ and $e = 2 \times 10^{-4}$ for one and two space dimensions, respectively.

In Eq. (44), we require that the threshold scales as $\Delta t^{p-1}$. This $\Delta t$ dependence is necessary to preserve the accuracy order *p* of the numerical results. Since the error of $[\vec{K}^{(l+1)}]_j$ is parametrically of order $\Delta t$ higher than that of $[\vec{K}^{(l)}]_j$, the error of $[\vec{K}^{(l+1)}]_j$ is of order $\Delta t^p$ when the condition (44) is met. This order of $[\vec{K}^{(l+1)}]_j$ is in accordance with the accuracy of RK in Eq. (29).<fnref n="5" /> On the right-hand side of Eq. (44), $\Delta x$ is used to make *p* dependence dimensionless. In Section 6 we show that the expected scaling (29) is nicely obtained using the criterion (44).

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) Equation (44) was not satisfied at cell *j* in the last update.

(2) Any of the cells adjacent to cell *j* were updated in the previous iteration and did not satisfy Eq. (44).

We update all the cells that need updating at once. We proceed with the iterations until Eq. (44) is reached for all the cells. We note that even once the condition (44) 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 *local optimization* in the following.

In the above procedure, we check only the most adjacent cells. Although the right-hand side of Eq. (31) also depends on the second neighboring cells, we have checked that their inclusion hardly changes the quality of the numerical results despite an extra computational cost.<fnref n="6" /> 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 Section 6.

## 4.2. Convergence condition and stiffness

As is common in iterative methods, the fixed-point method does not necessarily converge to the true solution. A necessary condition for the convergence of the fixed-point method is written

<footnote n="5">⁵The numerical error of one RK step is of order $\Delta t^{p+1}$. Also, notice $\Delta t$ on the right-hand side of Eq. (27).</footnote>

<footnote n="6">⁶The success of this prescription may rely on the fact that Eq. (31) only depends on at most second neighboring cells. When one employs a space scheme that depends on more distant cells, this procedure may need to be modified to keep the quality of the numerical results.</footnote>
in terms of the spectral radius R of the Jacobian matrix of $\vec{F}$ at the solution $\vec{K} = \vec{K}^*$:

$$R = \rho \left( \frac{\partial \vec{F}}{\partial \vec{K}}(\vec{K}^*) \right) < 1. \qquad (45)$$

Moreover, the convergence is faster for smaller R. From Eqs. (42) and (43) one finds that the Jacobian matrix $\partial \vec{F}/\partial \vec{K}$, and hence R, is proportional to $\Delta t$. Therefore, the convergence is guaranteed for a sufficiently small $\Delta t$ and is faster for a smaller $\Delta t$. However, the convergence is not guaranteed in general. In fact, we will see later that our method can be aborted for $\Delta t/\Delta x \gtrsim 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 forthcoming publications [N. Touroux, M. Kitazawa, K. Murase, M. Nahrgang (in preparation)]. In the present study, we simply neglect these cases since our method works stably in the range of $\Delta t$ that is practically important.

The convergence condition (45) is related to the stiffness of the system of equations. In fact, Eq. (45) is rewritten as

$$R = \rho \left( \Delta t \left[ a_{nm} \left. \frac{\partial \vec{h}}{\partial \vec{y}} \right|_{(t+c_n \Delta t, \vec{y}+\Delta t \sum_{l=1}^{S} a_{nl} \vec{k}_{(l)}^*)} \right]_{nm} \right) \\ = \Delta t \rho(a_{nm}) \rho \left( \frac{\partial \vec{h}(t, \vec{y})}{\partial \vec{y}} \right) + \mathcal{O}((\Delta t)^2), \qquad (46)$$

where $\rho(a_{nm})$ is the spectral radius of the Butcher tableau of a given RK method. Here, $\rho(\partial \vec{h}/\partial \vec{y})$ characterizes the stiffness of the system since the largest absolute eigenvalue $\max_{\lambda} |\lambda| = \rho(\partial \vec{h}/\partial \vec{y})$ (i.e. the inverse of the shortest timescale) gives the most severe impact on the stability. Equation (46) shows that R is proportional to the stiffness $\rho(\partial \vec{h}/\partial \vec{y})$.

This means that the convergence becomes slow in stiff systems. Since we locally truncate the iteration number of the fixed-point method using Eq. (44), the iteration number becomes automatically larger in the stiff spatial region where short-timescale phenomena are happening. In short, our local optimization automatically detects the stiffness of local regions.

### 4.3. Accuracy, efficiency, and stability

In later sections, we particularly consider the one-stage Gauss-Legendre (GL1) method [57,58,59] as an implicit method of order p = 2, whose Butcher tableau is given in the right-hand column of Table 1. In this subsection, we focus on the locally optimized GL1 and discuss its accuracy, efficiency, and stability.

We first show that the method applied to GL1 keeps the accuracy order p = 2 even when the iteration number is fixed to one. The single-iteration GL1 is explicitly written down as

$$\vec{y}(t + \Delta t) = \vec{y}(t) + \Delta t \vec{k}^{(1)}, \qquad (47)$$

$$\vec{k}^{(1)} = \vec{h} \left( t + \frac{\Delta t}{2}, \vec{y}(t) + \frac{\Delta t}{2} \vec{k}^{(0)} \right), \qquad (48)$$

where the initial guess $\vec{k}^{(0)}$ of the iteration is chosen to be the solution in the previous time step. By inverting Eq. (47) of the previous time step, the initial guess can be expressed as

$$\vec{k}^{(0)} = \frac{\vec{y}(t) - \vec{y}(t - \Delta t)}{\Delta t}. \qquad (49)$$
Plugging Eqs. (48) and (49) into Eq. (47), we obtain

$$ \vec{y}(t + \Delta t) = \vec{y}(t) + \Delta t \vec{h} \left( t + \frac{\Delta t}{2}, \frac{3\vec{y}(t) - \vec{y}(t - \Delta t)}{2} \right). \quad (50) $$

The Taylor expansion of this expression matches the analytic expansion up to the second order in $\Delta t$:

$$ \vec{y}(t + \Delta t) = \vec{y}(t) + \Delta t \vec{h}(t, \vec{y}(t)) + \frac{(\Delta t)^2}{2} [\partial_t + \vec{h}(t, \vec{y}(t)) \cdot \partial_{\vec{y}}] \vec{h}(t, \vec{y}(t)) + \mathcal{O}((\Delta t)^3). \quad (51) $$

Thus, the single-iteration fixed-point method for GL1 already has the desired accuracy of order $p = 2$.

Next, we see that the single-iteration GL1 can be related to AB2. By regarding

$$ \begin{aligned} \vec{z} \left( t + \frac{\Delta t}{2} \right) &= \frac{3\vec{y}(t) - \vec{y}(t - \Delta t)}{2} \\ &= \vec{y}(t - \Delta t) + \frac{3}{2} \Delta t \vec{h} \left( t - \frac{\Delta t}{2}, \vec{z} \left( t - \frac{\Delta t}{2} \right) \right). \end{aligned} \quad (52) $$

as state variables, Eq. (50) can be shown to be equivalent to the AB2 method:

$$ \vec{z}(t + \Delta t) = \vec{z}(t) + \frac{3}{2} \Delta t \vec{h}(t, \vec{z}(t)) - \frac{1}{2} \Delta t \vec{h}(t - \Delta t, \vec{z}(t - \Delta t)). \quad (53) $$

Here, the initial guess $\vec{k}^{(0)}$ plays an important role in the equivalence of the single-iteration GL1 (50) to AB2 (53). In the single-iteration GL1, the initial guess carries the information of the previous step $\vec{y}(t - \Delta t)$, which is similar to the linear multistep methods such as AB2 using the states of previous steps.

Both the single-iteration GL1 and AB2 are second-order methods with one-stage evaluation of $\vec{h}$, so they are more efficient than the second-order explicit RK with two stages. However, the problem with those one-stage methods is the stability in the stiff equations. Those methods effectively use an extrapolation using the previous step $\vec{y}(t - \Delta t)$ as observed in Eqs. (50) and (53) and are known to be less stable than the explicit RKs. This is a common problem with linear multistep methods such as the Adams-Bashforth method. In our method, this problem is automatically remedied by the iterations. As already discussed, the number of fixed-point iterations increases in the stiff region so that the method approaches the original GL1, which is A-stable. Thus, we can view our method as automatic switching between the two second-order methods, the less stable but efficient AB2 and the stable GL1, based on the local stiffness. In this way, our method benefits from both the efficiency of AB2 and the stability of GL1.

As we will observe in the numerical tests in the later sections, the stiff part appears near the vacuum or the region with large gradients. Such regions typically appear only in limited space-time volumes in practical applications, so we can expect the computational cost of the locally optimized implicit method to be close to the cost of AB2. Since AB2 with a one-stage evaluation of $\vec{h}(t, \vec{y}(t))$ has half the cost of the second-order explicit RK methods with two stages, we may expect a better performance of the locally optimized implicit method as long as the total size of the stiff regions is not too large.

The application of the idea of the locally optimized implicit method is not limited to relativistic hydrodynamics. For example, the idea can be applied to the general (discretized) field equations with local interactions. In principle, this idea can be applied to any ODEs with the locality, where a degree of freedom directly interacts with a fixed number of other degrees regardless of the size of the entire system. Thus, we expect broader applications of this idea in various fields.
Finally, we comment on the future extension of our analysis to the higher-order implicit RK methods. In this subsection, we considered the specific case of GL1 utilizing its particular form. It should be in principle possible to design the proper choice of the initial guess $\vec{K}^{(0)}$ for the general implicit RK methods by taking linear combinations of the solutions $\vec{K}$ in the previous time steps so that the accuracy order of the implicit method is preserved even with a single iteration. These analyses are left for future study.

# 5. Numerical setup

## 5.1. Runge-Kutta methods

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 ones. For the demonstration 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 tableau is given in the left-hand column of Table 1. For the implicit method, we choose the one-stage Gauss-Legendre method [57,58,59], which has the same order of accuracy as Heun's method. In Section 6, we compare the numerical results of these methods. In the following, we refer to Heun's and the 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 tableau 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.

## 5.2. Dimensions

In this study, we compare the implicit and explicit RK methods in (1+1)- and (2+1)D numerical simulations. For the (1+1)D case, we deal with the Riemann problem discussed in Section 2.2.1 with finite and zero minimum energy density $\epsilon_{min}$. We also perform (2+1)D simulations with the boost invariance for the Gubser flow and event-by-event initial conditions for heavy-ion collisions generated by TRENT0 [50].

In all cases, we consider a system of length 40 fm in each direction. To investigate the dependence on the mesh size $\Delta x$, we perform numerical simulations with $\Delta x = 0.2$ and $0.4$ fm where the spatial extent is divided into 200 and 100 cells, respectively. The mesh size $\Delta x = 0.2$ fm is a typical value used in applications to HIC [13,17,37]. The value $\Delta x = 0.4$ fm is employed for comparison. The time step $\Delta t$ is fixed in a single simulation, and we perform simulations with different $\Delta t$ down to $\Delta t = 0.1 \times 2^{-5}\Delta x$, which is determined so that our results in Section 6 contain enough points to investigate the accuracy order $p$ through Eq. (29). We also checked that our results behave consistently even at $\Delta t = 0.1 \times 2^{-6}\Delta x$ and $0.1 \times 2^{-7}\Delta x$ in several cases.

## 5.3. 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, $\epsilon_{num}(\mathbf{x}_j)$, with that of the exact solution $\epsilon_{exact}(\mathbf{x}_j)$ in the continuum case
using the measure

$$ \Delta_{\text{exact}}(\mathbf{x}_j) = D(\epsilon_{\text{num}}(\mathbf{x}_j), \epsilon_{\text{exact}}(\mathbf{x}_j)), \qquad (54) $$

with

$$ D(\epsilon_1, \epsilon_2) = \frac{\epsilon_1 - \epsilon_2}{\max(\epsilon_1, \epsilon_2)}. \qquad (55) $$

Here, we take the maximum of $\epsilon_1$ and $\epsilon_2$ for the denominator of Eq. (55), instead of either one of them, to avoid a situation where Eq. (54) becomes extremely large when one of them is vanishingly small. As for the second measure, we choose a numerical result obtained at small $\Delta t$, denoted by $\epsilon_{\text{ref}}(\mathbf{x}_j)$, as a reference. We then compare $\epsilon_{\text{num}}(\mathbf{x}_j)$ with the reference as

$$ \Delta_{\text{ref}}(\mathbf{x}_j) = D(\epsilon_{\text{num}}(\mathbf{x}_j), \epsilon_{\text{ref}}(\mathbf{x}_j)). \qquad (56) $$

For the reference $\epsilon_{\text{ref}}(\mathbf{x}_j)$, we use the result by the explicit method with the smallest $\Delta t = 0.1 \times 2^{-5}\Delta x$.

To compare numerical errors of different time discretization methods, $\Delta_{\text{ref}}(\mathbf{x}_j)$ is more convenient than $\Delta_{\text{exact}}(\mathbf{x}_j)$. This is because a finite space discretization error remains in $\Delta_{\text{exact}}(\mathbf{x}_j)$ even in the $\Delta t \to 0$ limit, while one can discuss genuinely the time discretization error using $\Delta_{\text{ref}}(\mathbf{x}_j)$.

In Section 6, we inspect the magnitude of the numerical error using the maximum and spatial average of $|\Delta_{\text{ref}}(\mathbf{x}_j)|$. However, we found that $|\Delta_{\text{ref}}(\mathbf{x}_j)|$ tends to be close to unity near the vacuum region where $\epsilon_{\text{exact}}(\mathbf{x}_j)$ is vanishingly small. In this region $\epsilon_{\text{num}}(\mathbf{x}_j)$ and $\epsilon_{\text{ref}}(\mathbf{x}_j)$ have nonvanishing values due to numerical diffusion, i.e. the unphysical propagation from the matter region, which is inevitable in computational hydrodynamics. Nevertheless, $\epsilon_{\text{num}}(\mathbf{x}_j)$ is still negligibly small and the large $|\Delta_{\text{ref}}(\mathbf{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 $\epsilon_{\text{num}}(\mathbf{x}_j)$ is smaller than a threshold value $\epsilon_{\text{thr}}$, where $\epsilon_{\text{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 of the boundaries to remove possible boundary effects. We then evaluate the maximum and the average of $|\Delta_{\text{ref}}(\mathbf{x}_j)|$ in the remaining cells.

## 5.4. 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_{\text{KT}}]_j$ in Eq. (31). This means that the number of evaluations $N_{\text{KT}}$ serves as an approximate measure of the computational cost. We thus use the number of evaluations per cell,

$$ n_{\text{KT}} = \frac{N_{\text{KT}}}{N_{\text{cell}}}, \qquad (57) $$

for the measure of the computational cost, where $N_{\text{cell}}$ stands for the number of total cells. For our explicit method with $S = 2$, $[h_{\text{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_{\text{KT}}]_j$ with $S = 1$. Moreover, due to the local optimization, $n_{\text{KT}}$ is not necessarily an integer. We have checked that $n_{\text{KT}}$ indeed reflects the real computational times well with appropriate optimizations.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Energy density $\epsilon(\mathbf{x}_j)$ (upper panel) and its error compared with the exact solution $\Delta_{\text{exact}}(\mathbf{x}_j)$ (lower panel) for the Riemann problem at $t_{\text{end}} = 15$ fm with $\epsilon_{\min} = 1 \text{ fm}^{-4}$ and $\Delta t = 0.1 \times 2^{-5} \Delta x$.</figcaption></figure>

# 6. Numerical results

## 6.1. 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 Section 2.2.1. This is a (1+1)D problem in flat space-time with initial conditions given by Eqs. (9) and (10). We numerically solve the hydrodynamic evolution from the initial time $t_{\text{init}} = 0$ to $t_{\text{end}} = 15$ fm. We set $c_s^2 = 1/3$ and $\epsilon_{\max} = 10 \text{ fm}^{-4}$ and perform our analysis for the two values of $\epsilon_{\min} = 1$ and $0 \text{ 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.

### 6.1.1. *Numerical solutions.*
In the upper panel of Fig. 3, we first show the numerical result of $\epsilon(\mathbf{x}_j)$ at $t = t_{\text{end}}$ for $\epsilon_{\min} = 1 \text{ fm}^{-4}$ as a function of $x/t$. The red (blue) curves show the results obtained by the explicit (implicit) method with $\Delta t = 0.1 \times 2^{-5} \Delta x$ while the dotted and dashed lines represent the results for $\Delta x = 0.4$ and $0.2$ fm, respectively. The black solid line shows the analytic solution $\epsilon_{\text{exact}}(\mathbf{x}_j)$. In the lower panel, we also show $\Delta_{\text{exact}}(\mathbf{x}_j)$.

One sees that the results of the explicit and implicit methods degenerate well at this $\Delta t$. This suggests that the error from the time discretization is well suppressed in both methods. By comparing the results of $\Delta x = 0.4$ and $0.2$ fm, one also finds that the numerical result approaches the analytic solution as $\Delta x$ becomes smaller. Nevertheless, a clear deviation is still observed at the smaller $\Delta x$. This is significant around the rarefaction front and shock at $x/t = -c_s$ and $v_{\text{shock}}$, respectively, which are shown in the enlarged subpanels of the upper panel. In particular, a large deviation $\Delta_{\text{exact}} \simeq 0.5$ survives even for the small $\Delta x$ at $x/t = v_{\text{shock}}$, which is a consequence of the discontinuity of $\epsilon_{\text{exact}}(\mathbf{x}_j)$ at this point.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Same as Fig. 3 but for $\epsilon_{\min} = 0$.</figcaption></figure>

In Fig. 4, the numerical result for $\epsilon_{\min} = 0$ is shown in the same manner as Fig. 3. The figure shows that $\Delta_{\text{exact}}(\mathbf{x}_j)$ becomes almost unity in the vacuum region $x/t > 1$. As explained in Section 5.3, this result comes from vanishing $\epsilon_{\text{exact}}(\mathbf{x}_j)$ in the vacuum region.

In the following analysis, we use $\Delta_{\text{ref}}(\mathbf{x}_j)$ defined in Eq. (56) in place of $\Delta_{\text{exact}}(\mathbf{x}_j)$ to analyze the time discretization error exclusively.

6.1.2. *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. (56) at $t = t_{\text{end}}$ for $\epsilon_{\min} = 1$ and $0 \text{ fm}^{-4}$, respectively, for $\Delta x = 0.4 \text{ fm}$ (upper plot) and $0.2 \text{ fm}$ (lower plot). The horizontal axis is the computational cost $n_{\text{KT}}$ in Eq. (57). The red squares and blue circles show the explicit and implicit results, respectively. The large point in each line shows the result for $\Delta t/\Delta x = 0.1$. The distance between the points along one line is obtained from varying $\Delta t$ by a factor of two. The numerical simulation fails when the time step $\Delta t$ is too large, so only successful results are shown in the figures. In the explicit method, the cost is $n_{\text{KT}} = S(t_{\text{end}} - t_{\text{init}})/\Delta t = 30 \text{ fm}/\Delta t$.

We observe several notable features in Figs. 5 and 6. First, all the results are approximately proportional to $n_{\text{KT}}^{-2}$. 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 $|\Delta_{\text{ref}}|$ in both maximum and average numerical errors; the former is about 2–3 times as efficient as the latter.<fnref n="7" /> This means that with the same accuracy, the implicit method is faster than the explicit one. This is a remarkable

<footnote n="7">⁷One might want to compare the values of $|\Delta_{\text{ref}}|$ at the same $\Delta t$ for different $\Delta x$. However, $|\Delta_{\text{ref}}|$ represents the difference from $\epsilon_{\text{ref}}(x)$, not $\epsilon_{\text{exact}}(x)$, so it does not contain the spatial discretization error. Since $\epsilon_{\text{ref}}(x)$ is different for each $\Delta x$, their values at different $\Delta x$ are not comparable.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Maximum and average of the numerical error $|\Delta_{\text{ref}}|$ defined in Eq. (56) as functions of the computational cost $n_{\text{KT}}$ for the Riemann problem with $\epsilon_{\min} = 1 \text{ 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 $\Delta x = 0.4$ and $0.2$ fm, respectively. The large points show the results at $\Delta t/\Delta x = 0.1$, while the other small points are obtained by varying $\Delta t$ by a factor of 2. Note that the absolute magnitudes of errors with different $\Delta x$ are not comparable because the reference $\epsilon_{\text{ref}}(\mathbf{x}_j)$ depends on $\Delta x$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same result as Fig. 5, but for $\epsilon_{\min} = 0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Number of iterations required for the fixed-point solver at each space-time point for the Riemann problem with $\epsilon_{\min} = 1 \text{ fm}^{-4}$. The dashed and dash-dotted lines show the positions of the shock and rarefaction front, respectively. Four results for $(\Delta t/\Delta x, \Delta x) = (0.1, 0.4)$, $(0.1 \times 2^{-3}, 0.4)$, $(0.1, 0.2)$, and $(0.1 \times 2^{-3}, 0.2)$ are shown in the top-left, top-right, bottom-left, and bottom-right panels.</figcaption></figure>

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 $\Delta 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 $\Delta t/\Delta x = 0.1$ and $0.1 \times 2^{-3} \text{ fm}$ (left and right), and $\Delta x = 0.4$ and $0.2 \text{ fm}$ (top and bottom). The dashed and dash-dotted lines show the locations of the shock and the rarefaction front at which $\epsilon_{\text{exact}}(x)$ and $\partial\epsilon_{\text{exact}}(x)/\partial 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 Section 4, 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_{\text{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 could 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 $\Delta t/\Delta x \ge 1.6$. This result is consistent with the CFL condition [51,52] for stability, $\Delta t/\Delta x < C_{\text{CFL}}$ where $C_{\text{CFL}} = O(1)$ depends on the detailed scheme. On the other hand, failures of the implicit method can occur for $\Delta t/\Delta x \gtrsim 0.4$. These failures happen when the fixed-point solver does not converge. Despite the A-stability of the implicit method itself, the fixed-point solver is not guaranteed to converge for large $\Delta t$.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Same result as Fig. 7, but for $\epsilon_{min} = 0$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Energy density $\epsilon(x_j)$ (upper panels), error $\Delta_{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 $\tau$ for $\Delta t = 0.1 \times 2^{-5}\Delta x$ and $\Delta x = 0.4$ fm.</figcaption></figure>

## 6.2. Gubser flow

Next, we perform the same comparison for the Gubser flow as discussed in Section 2.2.2. We solve transverse (2+1)D hydrodynamics assuming the boost invariance. In this case, the solution contains neither discontinuity nor vacuum region. We employ $\epsilon_0 = 1 \text{ fm}^{-4}$ and $q = 1$ for the parameters in Eq. (19) and set $\tau_{init} = 1 \text{ fm}$ for the initial proper time. We then solve the hydrodynamic equations numerically until the end time $\tau_{end} = 10 \text{ fm}$.

In the top panels of Figs. 9 and 10, we first show the energy density at several values of $\tau$ for $\Delta x = 0.4$ and $0.2$ fm, respectively, obtained by the implicit method with $\Delta t = 0.1 \times 2^{-5}\Delta x$ fm. Here we do not show the result of the explicit method since its behavior is almost the same as
<figure><img src="image_11.png" /><figcaption>Fig. 10. Same results as Fig. 9, but with $\Delta x = 0.2$ fm.</figcaption></figure>

anticipated from the discussion in Section 6.1. The panels show that the numerical solution has noncircular behaviors, while the exact solution (19) has azimuthal symmetry. Correspondingly, $\Delta_{\text{exact}}(\mathbf{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 $\Delta t$. By comparing Figs. 9 and 10, one also finds that the deviation from the exact solution is suppressed as $\Delta 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 $\tau$. These panels show that two iterations are needed around the region where $\epsilon(\mathbf{x}_j)$ changes rapidly, but only one iteration is enough for the other wide area. This result is consistent with the one in Section 6.1.

In Fig. 11, we compare the maximum and the average of $|\Delta_{\text{ref}}(\mathbf{x}_j)|$ in the explicit and implicit methods as functions of $n_{\text{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.

### 6.3. 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 initial conditions using TRENTO [50]. We perform (2+1)D event-by-event simulations for lead-lead collisions with impact parameter $b = 3$ fm and collision energy $\sqrt{s_{NN}} = 2.76$ TeV. We use $\tau_{\text{init}} = 0.48$ fm as the initial time of the hydrodynamics. For the other parameters of TRENTO, we take the same values as in Refs. [24,25]. For the EOS, we use the parametrized form [20,60] of the Wuppertal-Budapest lattice EOS [61]. We then solve the hydrodynamic evolution until $\tau_{\text{end}} = 10$ fm.

In the upper panels of Figs. 12 and 13, we show the transverse profile of energy density $\epsilon(\mathbf{x}_j)$ for one event in the x-y plane obtained with the finest $\Delta t = 0.1 \times 2^{-5} \Delta x$ in the implicit method for $\Delta x = 0.4$ and 0.2 fm, respectively. The cells satisfying $\epsilon(\mathbf{x}_j) < \epsilon_{\text{thr}}$ are shown by the white
<figure><img src="image_12.png" /><figcaption>Fig. 11. Same as Fig. 5, but for the analysis of the Gubser flow in Figs. 9 and 10.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Energy density $\epsilon(x_j)$ (upper panels) and the number of iterations of the fixed-point solver in each cell (bottom panels) for a TRENTO initial condition (lead-lead collision with $b = 3$ fm and $\sqrt{s_{NN}} = 2.76$ TeV) at several values of proper time $\tau$ for $\Delta t = 0.1 \times 2^{-5}\Delta x$ and $\Delta x = 0.4$ fm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Same as Fig. 12, but for $\Delta x = 0.2$ fm.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Same as Fig. 5, but for 10 initial conditions generated by TRENT0.</figcaption></figure>

color. As discussed in Section 5.3, 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 $|\Delta_{\text{ref}}|$ as functions of $n_{\text{KT}}$ for 10 initial conditions generated randomly by TRENT0. 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.

## 7. Conclusion

In this study we have 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 Kurganov-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. We also showed that our locally optimized implicit method can be viewed as an adaptive mixing of two second-order methods: it corresponds to the two-step Adams-Bashforth method when a single iteration is performed in smooth regions, while it approaches the one-stage Gauss-Legendre method as the number of iterations increases in stiffer regions. As a result, it benefits from the efficiency of the two-step Adams-Bashforth method and the stability of the one-stage Gauss-Legendre method.

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 $\Delta t/\Delta x$. The advantage of the former mainly comes from the fact that the fixed-point solver requires fewer 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 fewer 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)D viscous hydrodynamics 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 $\Delta t/\Delta x$. A solution to improve the convergence will be discussed in a future publication. We have already started an implementation of viscous hydrodynamics in (3+1) dimensions with the implicit method. As will be reported in a forthcoming publication, our preliminary analyses show that our method is advantageous even in this case. We also confirmed the efficiency of our method in solving the diffusion equation and the Klein-Gordon equation. These results indicate the advantage of our method in a broader class of problems.

### Acknowledgments

The authors thank Marcus Bluhm and Gregoire Pihan for fruitful discussions in the initial phase of this work.

### Funding

N.T. is supported by MEXT (Ministry of Education, Culture, Sports, Science and Technology). This work was supported by JSPS (Japan Society for the Promotion of Science) KAKENHI (Grant Nos. JP19H05598, JP20H01903, JP22K03619, JP23H04507, and JP23K13102) and the TYL-FJPPL (Toshiko Yuasa Laboratory - France Japan Particle Physics Laboratory) program of IN2P3-CNRS (National Institute of Nuclear and Particle Physics - National Centre for Scientific Research) and KEK (High Energy Accelerator Research Organization). M.N. and N.T. acknowledge the support of the program "Etoiles montantes en Pays de la Loire 2017".

## Appendix A. Comparison with the midpoint method

In Section 6, we compare the computational efficiencies of Heun's and the GL1 methods as representatives of explicit and implicit methods, respectively. In this appendix, we employ the midpoint method as another explicit method and compare its efficiency with Heun's and the GL1 methods. The Butcher tableau of the midpoint method is shown in the middle column of Table 1. As shown there, this method has the same order of accuracy $p = 2$ as Heun's and the GL1 methods.

In Figs. A1 and A2, we plot the maximum and the average of the numerical error $|\Delta_{\text{ref}}|$ as functions of $n_{\text{KT}}$ for the (1+1)D Riemann problem with $\epsilon_{\min} = 1$ and $0 \text{ fm}^{-4}$, respectively.
<figure><img src="image_16.png" /><figcaption>**Fig. A1.** Maximum and average of the numerical error $|\Delta_{\text{ref}}|$ as functions of $n_{\text{KT}}$ for the Riemann problem with $\epsilon_{\min} = 1 \text{ fm}^{-4}$. The result of the midpoint method is shown by the green triangle symbols. The results for Heun's and the GL1 methods are the same as in Fig. 5.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. A2.** Same result as Fig. A1, but for $\epsilon_{\min} = 0 \text{ fm}^{-4}$.</figcaption></figure>
The results show that the performance of the midpoint method is almost the same as Heun's method.

## References

[1] U. W. Heinz and M. Jacob, arXiv:nucl-th/0002042 (Preprint) [Search in SPIRE].

[2] M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30 (2005) [arXiv:nucl-th/0405013] [Search in SPIRE].

[3] B. Muller and J. L. Nagle, Ann. Rev. Nucl. Part. Sci. 56, 93 (2006) [arXiv:nucl-th/0602029] [Search in SPIRE].

[4] D. A. Teaney, Quark-gluon plasma 4, 207-266 (2010) [arXiv:0905.2433] [Search in SPIRE].

[5] B. V. Jacak and B. Muller, Science 337, 310 (2012).

[6] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013) [arXiv:1301.2826] [Search in SPIRE].

[7] R. Derradi de Souza, T. Koide, and T. Kodama, Prog. Part. Nucl. Phys. 86, 35 (2016) [arXiv:1506.03863] [Search in SPIRE].

[8] P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Phys. Rept. 621, 76 (2016) [arXiv:1510.00442] [Search in SPIRE].

[9] P. Romatschke and U. Romatschke, Relativistic Fluid Dynamics In and Out of Equilibrium (Cambridge University Press, Cambridge, UK, 2019), Cambridge Monographs on Mathematical Physics [arXiv:1712.05815] [Search in SPIRE].

[10] W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018) [arXiv:1802.04801] [Search in SPIRE].

[11] H. Song and U. W. Heinz, Phys. Rev. C 77, 064901 (2008) [arXiv:0712.3715] [Search in SPIRE].

[12] A. K. Chaudhuri, arXiv:0801.3180 (Preprint) [Search in SPIRE].

[13] B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010) [arXiv:1004.1408] [Search in SPIRE].

[14] V. Roy and A. K. Chaudhuri, Phys. Rev. C 85, 024909 (2012); 85, 049902 (2012) [erratum] [arXiv:1109.1630] [Search in SPIRE].

[15] P. Bozek, Phys. Rev. C 85, 034901 (2012) [arXiv:1110.6742] [Search in SPIRE].

[16] L. Del Zanna, V. Chandra, G. Inghirami, V. Rolando, A. Beraudo, A. De Pace, G. Pagliara, A. Drago, and F. Becattini, Eur. Phys. J. C 73, 2524 (2013) [arXiv:1305.7052] [Search in SPIRE].

[17] I. Karpenko, P. Huovinen, and M. Bleicher, Comput. Phys. Commun. 185, 3016 (2014) [arXiv:1312.4160] [Search in SPIRE].

[18] Y. Akamatsu, S.-I. Inutsuka, C. Nonaka, and M. Takamoto, J. Comput. Phys. 256, 34 (2014) [arXiv:1302.1665] [Search in SPIRE].

[19] K. Murase, Causal hydrodynamic fluctuations and their effects on high-energy nuclear collisions, Ph.D. Thesis, Tokyo University (2015).

[20] D. Bazow, U. W. Heinz, and M. Strickland, Comput. Phys. Commun. 225, 92 (2018) [arXiv:1608.06577] [Search in SPIRE].

[21] K. Okamoto, Y. Akamatsu, and C. Nonaka, Eur. Phys. J. C 76, 579 (2016) [arXiv:1607.03630] [Search in SPIRE].

[22] L.-G. Pang, H. Petersen, and X.-N. Wang, Phys. Rev. C 97, 064918 (2018) [arXiv:1802.04449] [Search in SPIRE].

[23] L. Du and U. Heinz, Comput. Phys. Commun. 251, 107090 (2020) [arXiv:1906.11181] [Search in SPIRE].

[24] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 101, 024911 (2020).

[25] J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nat. Phys. 15, 1113 (2019).

[26] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. Lett. 126, 202301 (2021) [arXiv:2010.15130] [Search in SPIRE].

[27] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. C 103, 054909 (2021) [arXiv:2010.15134] [Search in SPIRE].

[28] J. E. Parkkila, A. Onnerstad, and D. J. Kim, Phys. Rev. C 104, 054904 (2021) [arXiv:2106.05019] [Search in SPIRE].

[29] G. Nijs and W. van der Schee, Phys. Rev. C 106, 044903 (2022) [arXiv:2110.13153] [Search in SPIRE].
[30] M. Asakawa and M. Kitazawa, Prog. Part. Nucl. Phys. **90**, 299 (2016) [arXiv:1512.05038] [Search inSPIRE].

[31] M. Bluhm et al. Nucl. Phys. A **1003**, 122016 (2020) [arXiv:2001.08831] [Search inSPIRE].

[32] J. Adam et al. Phys. Rev. Lett. **126**, 092301 (2021) [arXiv:2001.02852] [Search inSPIRE].

[33] J. Adamczewski-Musch et al. Phys. Rev. C **102**, 024914 (2020) [arXiv:2002.08701] [Search inSPIRE].

[34] H. Petersen, J. Steinheimer, G. Burau, M. Bleicher, and H. Stöcker, Phys. Rev. C **78**, 044901 (2008) [arXiv:0806.1695] [Search inSPIRE].

[35] H. Holopainen, H. Niemi, and K. J. Eskola, Phys. Rev. C **83**, 034901 (2011) [arXiv:1007.0368] [Search inSPIRE].

[36] M. McNelis, D. Bazow, and U. Heinz, Comput. Phys. Commun. **267**, 108077 (2021) [arXiv:2101.02827] [Search inSPIRE], doi:10.1016/j.cpc.2021.108077.

[37] T. Hirano, K. Tsuda, and K. Kajimoto, arXiv:nucl-th/0011087 (Preprint) [Search inSPIRE].

[38] Y. Tachibana, Hydrodynamic response to jet propagation in quark-gluon plasma, Ph.D. Thesis, Tokyo University (2014).

[39] J. I. Kapusta, B. Muller, and M. Stephanov, Phys. Rev. C **85**, 054906 (2012) [arXiv:1112.6405] [Search inSPIRE].

[40] K. Murase and T. Hirano, Nucl. Phys. A **956**, 276 (2016) [arXiv:1601.02260] [Search inSPIRE].

[41] M. Singh, C. Shen, S. McDonald, S. Jeon, and C. Gale, Nucl. Phys. A **982**, 319 (2019) [arXiv:1807.05451] [Search inSPIRE].

[42] M. Bluhm, M. Nahrgang, T. Schäfer, and S. A. Bass, EPJ Web Conf. **171**, 16004 (2018) [arXiv:1804.03493] [Search inSPIRE].

[43] A. Sakai, K. Murase, and T. Hirano, Phys. Rev. C **102**, 064903 (2020) [arXiv:2003.13496] [Search inSPIRE].

[44] M. Nahrgang, S. Leupold, C. Herold, and M. Bleicher, Phys. Rev. C **84**, 024912 (2011) [arXiv:1105.0622] [Search inSPIRE].

[45] C. Herold, M. Nahrgang, Y. Yan, and C. Kobdaj, Phys. Rev. C **93**, 021902 (2016) [arXiv:1601.04839] [Search inSPIRE].

[46] K. Nakamura, T. Miyoshi, C. Nonaka, and H. R. Takahashi, Phys. Rev. C **107**, 014901 (2023) [arXiv:2209.00323] [Search inSPIRE].

[47] A. Kurganov and E. Tadmor, J. Comput. Phys. **160**, 241 (2000).

[48] F. D. Lora-Clavijo, J. P. Cruz-Perez, F. S. Guzman, and J. A. Gonzalez, Rev. Mex. Fis. E **59**, 28 (2013) [arXiv:1303.3999] [Search inSPIRE].

[49] S. S. Gubser, Phys. Rev. D **82**, 085027 (2010).

[50] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C **92**, 011901 (2015).

[51] R. Courant, K. Friedrichs, and H. Lewy, Math. Ann. **100**, 32 (1928).

[52] R. Courant, K. Friedrichs, and H. Lewy, IBM J. Res. Dev. **11**, 215 (1967).

[53] J. D. Bjorken, Phys. Rev. D **27**, 140 (1983).

[54] A. H. Taub, Phys. Rev. **74**, 328 (1948).

[55] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon Press, Oxford, UK, 1987).

[56] E. Hairer and G. Wanner, J. Comput. Appl. Math. **111**, 93 (1999).

[57] I. Arieh, A First Course in the Numerical Analysis of Differential Equations (Cambridge University Press, Cambridge, UK, 1996).

[58] J. Kuntzmann and Z. Angew Math. Mech. **41**, T28 (1961).

[59] J. C. Butcher, Math. Comput. **18**, 50 (1964).

[60] D. Bazow, GPU-VH (2017) (available at: https://github.com/bazow/gpu-vh, date last accessed May 17, 2024).

[61] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, J. High Energy Phys. **1011**, 077 (2010) [arXiv:1007.2580] [Search inSPIRE].