A new relativistic hydrodynamics code for high-energy heavy-ion collisions

We construct a new Godunov type relativistic hydrodynamics code in Milne coordinates, using a Riemann solver based on the two-shock approximation which is stable under the existence of large shock waves. We check the correctness of the numerical algorithm by comparing numerical calculations and analytical solutions in various problems, such as shock tubes, expansion of matter into the vacuum, the Landau-Khalatnikov solution, and propagation of fluctuations around Bjorken flow and Gubser flow. We investigate the energy and momentum conservation property of our code in a test problem of longitudinal hydrodynamic expansion with an initial condition for high-energy heavy-ion collisions. We also discuss numerical viscosity in the test problems of expansion of matter into the vacuum and conservation properties. Furthermore, we discuss how the numerical stability is affected by the source terms of relativistic numerical hydrodynamics in Milne coordinates.


Introduction
Relativistic hydrodynamics has been widely used for the description of macroscopic dynamics in various fields ranging from nuclear physics to astrophysics. The high-energy heavy-ion collision experiment is one of the active areas of relativistic hydrodynamics applications.
In 2005 at the relativistic heavy-ion collider (RHIC), the production of strongly interacting quark-gluon plasma (QGP) was achieved, which was supported not only by the experimental data but also theoretical analyses [1]. a e-mail: okamoto@hken.phys.nagoya-u.ac.jp Studies based on relativistic hydrodynamics have shown remarkable success in understanding various observables such as particle distributions, collective flows, particle correlations, and so on [2][3][4][5]. The strong elliptic flow at RHIC is a highlight of the success of hydrodynamic models and is one piece of evidence that the QGP is not a weakly interacting gas but a strongly interacting matter. Since then the hydrodynamic model has been one of the promising phenomenological models for the description of dynamics of hot and dense matter produced in the heavy-ion collisions.
The construction of relativistic viscous hydrodynamic model has been of practical importance and has formed a basis for the analyses of the heavy-ion collisions [6][7][8][9][10][11][12]. In the last decade, the hydrodynamic model itself has also been developed through the analyses of experimental data of heavy-ion collisions at RHIC and the large hadron collider (LHC). By comparing the hydrodynamic model calculations and the experimental observables such as particle distributions and collective flows, detailed bulk properties of QGP such as the QCD equation of state and its transport coefficients have been investigated. Also, physical QCD equation of state is now available by lattice QCD simulations at vanishing chemical potential [13,14] and is applied to the hydrodynamic model. This enables us to bridge the first principle lattice QCD simulations and the experimental data in the heavy-ion collisions. In spite of the success of hydrodynamic models in high-energy heavy-ion collisions, there are still several issues under discussion. Currently the hydrodynamic models often adopt Israel-Stewart theory [15] and a second-order viscous hydrodynamics from AdS/CFT correspondence [16] as their basic equations. However, we have not reached a conclusion on which relativistic viscous hydrodynamic equation is suitable for the description of relativistic heavy-ion collisions. This is because the extension from a relativistic ideal hydrodynamic equation to a viscous hydro-dynamic equation is not straightforward and several possible candidates exist. Also viscous second-order anisotropic hydrodynamics is proposed which reproduces the exact solution of the Boltzmann equation in the relaxation-time approximation [17]. Furthermore it remains an enormous challenge to understand why hydrodynamics can be applied to the dynamics shortly after a heavy-ion collision takes place. Conclusive understanding of the mechanism of thermalization and hydrodynamization on such a short time scale is still missing.
Here we emphasize that a numerical algorithm for solving the relativistic hydrodynamic equation is one of the important ingredients in developing the hydrodynamic models. Recent high statistical experimental data at RHIC and the LHC imposed a more rigorous numerical treatment on the hydrodynamical models. For example, at RHIC and the LHC, higher harmonic anisotropic flow, which is expressed by the higher Fourier coefficient of particle yields as a function of azimuthal angle, is reported [18][19][20][21][22]. The origin of the higher harmonics is considered to be event-by-event initial fluctuations in the particle distributions. When comparing with those high statistical data, reducing the numerical dissipation of the numerical algorithm for relativistic hydrodynamic equations should allow us an access to more precise value of transport coefficients of the QGP. Usually each algorithm has advantages or disadvantages in terms of coding, computational time, numerical precision, and stability. Up to now, unfortunately, only little attention has been paid to the numerical aspects in the hydrodynamic models for high-energy heavy-ion collisions.
Recently we developed a state-of-the-art numerical algorithm for solving the relativistic hydrodynamic equation with the QGP equation of state [23]. In the algorithm, we use a Riemann solver based on the two-shock approximation [25][26][27][28] which is stable under the existence of large shock waves [29]. The new numerical scheme is stable even with a small numerical viscosity and can reduce the numerical uncertainly when extracting the physical viscosity of the QGP from the experimental data. However, this algorithm in Ref. [23] is developed in Cartesian coordinates. Meanwhile, at the high-energy heavy-ion collisions such as RHIC and the LHC, the expansion in longitudinal direction is rapid compared with that in transverse direction. For the description of a space-time evolution of high-energy heavyion collisions, Milne coordinates are more suitable than Cartesian coordinates. Therefore we extend our algorithm of relativistic ideal hydrodynamics in Cartesian coordinates to that in Milne coordinates so that we can efficiently apply it to the analyses of high-energy heavy-ion collisions. The algorithm that we shall present here plays an important role in solving the relativistic viscous hydrodynamic equation numerically [23,24]. For the viscous hydrodynamics, we split the hydrodynamic equations into an ideal part and a viscous part. The ideal part can be solved by the Riemann solver for ideal hydrodynamics.
The present article is organized as follows. We begin in Sect. 2 by showing the basic equations for the hydrodynamic models in Milne coordinates. In Sect. 3 we explain the numerical algorithm; Riemann problem in Milne coordinates and our numerical scheme. Section 4 is devoted to several numerical tests, such as relativistic shock tubes and a comparison with analytic solutions which describe the dynamics of realistic high-energy heavy-ion collisions. In addition, we discuss the propagation of longitudinal fluctuations. In Sect. 5 we investigate the conservation property of our code. We end in Sect. 6 with our conclusions.

Relativistic hydrodynamics
Relativistic hydrodynamics is based on the conservation equations of net charge, energy, and momentum, where J µ is the baryon number current and T µν is the energy-momentum tensor. For the ideal fluid, the energy-momentum tensor and the baryon number current are given by where n is the baryon number density, e is the energy density, p is the pressure, u µ is the normalized four-velocity of the fluid, u µ u µ = 1 and g µν is the metric tensor.
In high-energy heavy-ion collisions at RHIC and the LHC, approximate invariance under the longitudinal Lorentz boost is observed in particle rapidity distributions around mid-rapidity [30][31][32][33][34][35][36]. In such situations, Milne coordinates are suitable for the description of space-time evolution of the hot and dense matter after the collisions. Milne coordinates η and τ are described by the rapidity η = tanh −1 (z/t) and the proper time τ = √ t 2 − z 2 with Cartesian coordinates. The coordinate transformation of the four-velocity between Milne coordinates and Cartesian coordinates is given by where the transverse components of the four-velocity, u x and u y are the same in both coordinates. The four-velocity in Milne coordinates is written by three-dimensional velocity as where w i = u i /u τ (i = x, y, η) is the three-dimensional velocity in Milne coordinates and W represents the Lorentz factor, The coordinate transformation of the three-dimensional velocity vector between Milne coordinates and Cartesian coordinates is given by where v i = u i /u t (i = x, y, z) is the three-dimensional velocity in Cartesian coordinates. In contrast to u x and u y , x and y components of the three-dimensional velocity in Milne coordinates are different from those in Cartesian coordinates. The metric tensor is given by g αβ = diag(1, −1, −1, −1/τ 2 ) and the nonzero components of the Christoffel symbols are In Milne coordinates, the charge conservation equation Eq. (1) and the equation of energy and momentum conservation Eq. (2) are written by There are geometric source terms in the right-hand side of Eqs. (12)- (15), which contain the effect from the coordinate expansion with τ. One can rewrite Eqs. (12)- (15) Here the effect from the coordinate expansion with τ is absorbed into the Jacobian τ in the derivative terms. There are the source terms in the right-hand side of Eqs. (18) and (19), which indicates that T ττ and T τη are not the conserved quantities. Instead of them, the conserved quantities are T τt and T τz . T τt and T τz are related to T ττ and T τη through the coordinate transformation between Milne and Cartesian coordinates, Using the conserved quantities, T τt and T τz , one can express the hydrodynamic equations in the conservative forms [37]

Riemann problem in Milne coordinates
The Riemann problem is an initial-value problem for the hydrodynamic equation. The initial condition is given by two arbitrary constant hydrodynamic states V V V L and V V V R separated by a discontinuity, where t 0 and z i stand for an initial time and a location of the discontinuity, respectively. The hydrodynamic state contains information on fluid variables, namely baryon number density, fluid velocities, and pressure. The initial discontinuity at z = z i decays into three nonlinear waves [38,39]. Two of them are shock waves and/or rarefaction waves. The other is a contact discontinuity moving with hydrodynamic flow. These waves evolve between the constant hydrodynamic states V V V L and V V V R with a constant velocity. The hydrodynamic states V V V L and V V V R do not change until characteristic information from the discontinuity arrives. Therefore, the hydrodynamic state outside the light cone of the discontinuity (t 0 , We can also define a Riemann problem in Milne coordinates. The initial condition of the Riemann problem in Milne coordinates is set to where τ 0 and η i are the initial proper time and the location of the discontinuity and they represent the same point as (t 0 , z i ) in Eq. (26) in Cartesian coordinates. Note that the components of V V V are not (n, w x , w y , w η , p) but the same as those in Eq. (26), (n, v x , v y , v z , p). The velocity fields in Cartesian coordinates u i /u t = (v x , v y , v z ) are constant in the rapidity direction. However, the velocity fields in Milne coordinates u i /u τ = (w x , w y , w η ) depend on rapidity. Now we show that the analytical solution for the Riemann problem in Milne coordinates is obtained from that in Cartesian coordinates by proper coordinate transformations, in which the key issue is to represent the hydrodynamic states as variables independent of η. Now we compare the two initial-value problems Eqs. (26) and (27). Without loss of generality, we can assume that the initial discontinuity represented by the solid circle is located at z i = η i = 0 as in Figs.1 (a) and (b). In Fig. 1 (a), the thick solid line stands for the time at which we define the initial condition of the Riemann problem in Cartesian coordinates, Eq. (26). In Milne coordinates, the initial condition of the Riemann problem is set on a hyperbola τ = τ 0 as in Fig. 1 (b). By comparison between Figs.1 (a) and (b), the hyperbola τ = τ 0 in Milne coordinates is located inside the constant hydrodynamic state V V V L or V V V R in Cartesian coordinates. This suggests that the initial condition of Eq. (27) is satisfied by the solution of the Riemann problem in Cartesian coordinates Eq. (26). In brief, the Riemann problems in both coordinates are identical and thus so are their solutions. A detailed explanation as regards this proof is given in Appendix A.

Numerical scheme
Assuming that the hydrodynamic variables in the x and y directions are constant, from Eqs. (23)-(25) we obtain where ν = t, x, y or z. If the transverse components of the four-velocity u x and u y are vanishing, Eq. (28) expresses the one-dimensional longitudinal expansion. In our numerical scheme, we utilize the Lagrange step [40] in which the grid-cell boundary itself moves together with the hydrodynamic flow during a time step from τ n to τ n+1 = τ n + ∆ τ. We discretize Eq. (28) by space-time integration in a grid cell based on the Lagrangian approach (Fig. 2), where α = τ or η. Here η i is the location of the ith grid-cell boundary at the proper time τ n , δ η i (τ ′ ) expresses a moving distance of the grid-cell boundary from τ = τ n to τ = τ n + τ ′ (0 ≤ τ ′ ≤ ∆ τ), and η i + δ η i (τ ′ ) indicates the location of the ith grid-cell boundary at the proper time τ. The center of the ith grid cell is located at η i + ∆ η/2 at τ = τ n , where ∆ η = η i − η i−1 . Using Gauss' theorem for integration of Eq. (29), we find the value of T τν of the ith grid cell at the next time step τ n+1 , where ∆ η lag i ≡ ∆ η + (δ η i (∆ τ) − δ η i−1 (∆ τ)) is the Lagrange grid-cell size at the proper time τ n+1 , C i is the trajectory of the Lagrange grid-cell boundary at η i and n α,i is the unit normal vector to C i (Fig. 2). The average values of the conserved quantities in the grid cell at the proper times τ n and τ n+1 are defined respectively, by The second term of Eq. (30) indicates the flux of the conserved quantities passing through the grid-cell boundary.
Using the analytical solution of the Riemann problem in Milne coordinates, we evaluate the ∆ η lag i and the flux term in Eq. (30). The grid-cell boundary corresponds to the initial discontinuity in the Riemann problem and moves together with the contact discontinuity in the solution of the Riemann problem. This indicates that the physical quantities on the trajectory δ η i (τ ′ ) are given by those on the contact discontinuity. We obtain the moving distance δ η i (∆ τ) of the grid-cell boundary In Eqs. (33) and (34), V z 0.i is the velocity of the grid-cell boundary at η i seen from an observer sitting at η = η i in Milne coordinates. In the construction of algorithm, we use the Lorentz boost transformation which is explained in the next paragraph. Here we show the explicit form of Eq. (30).
For detailed derivation of it please see Appendix B. Up to the third order in ∆ τ, the flux terms are given by where P i is the pressure on the cell boundary located at η i . From now we explain the numerical algorithm for solving the discretized equation Eq. (30). We express the hydrodynamic variables in Milne coordinates as W W W ≡ (n, w x , w y , w η , p). The first step is the interpolation procedure in which the left and right hydrodynamic states at the grid-cell boundary are determined from the reconstruction of the distribution of volume-averaged hydrodynamic variables W W W in a grid cell. If a linear interpolation method is used for reconstruction of the distribution of hydrodynamic variables, second-order accuracy is achieved (the MC limiter [40]). For third-order accuracy, we need to utilize a quadratic curve (the piecewise parabolic method (PPM) [28,41,42]) in reconstruction of the distribution of hydrodynamic states. For the test calculation in the next section, we use the PPM. Next, using the constructed left and right states W W W S,i in Milne coordinates, we prepare the initial condition of the Riemann problem in terms of V V V S,i in Cartesian coordinates. To obtain V V V S,i , we move a grid cell to η = 0 with the Lorentz boost transformation and perform the coordinate transformation to V V V . Here W W W is invariant under the Lorentz boost transformation. Using these relations, v x,y = w x,y and v z = τw η at η = 0, we obtain In the second step, we solve the Riemann problem Eq. (26) with the initial condition Eq. (38). P i and V z 0,i in Eqs. (33)-(37) are determined by the analytical solution of the Riemann problem. 1 For solving the Riemann problem we employ the two-shock approximation [25][26][27][28].
In the approximation, we can avoid solving the ordinary differential equation for the rarefaction wave which takes a lot of computational time in the multidimensional problem.
The numerical procedure of this step depends on the equation of state. The Riemann solution with the QCD equation of state is described in Ref. [23]. After solving the Riemann problem, we perform the inverse Lorentz transformation to the original frame. Here ∆ η lag i is Lorentz invariant. The flux terms Eqs. (35)-(37) are defined on the original frame but written in terms of V z 0,i . In the third step, solving the discretized hydrodynamic equation Eqs. (30), (33)-(37) with the values of P i and V z 0,i , we obtain distribution of the conservative quantities T τν (ν = t, x, y, z) at the next time step using the Lagrange scheme. We remap the grids which move in the Lagrange step on the Eulerian coordinate [40]. 2 In the final step, we construct the primitive variables W W W from the conserved quantities [23]. We obtain T τα (α = τ, x, y, η) from T τ µ (µ = t, x, y, z) from the coordinate transformation. The numerical method for construction of W W W from T τα is the same as that in Ref. [23].

(3+1) dimensional systems
The one-dimensional code is easily extended to a multidimensional code by the Strang splitting method [43], that is to say, multidimensional hydrodynamic evolution is realized by successive one-dimensional hydrodynamic calculations. To avoid counting the expansion effect of coordinates more than once, we extract one-dimensional hydrodynamic equations from Eqs. (12)- (15), which have source terms in the right-hand side. Then we rewrite each one-dimensional hydrodynamic equation as that without the source terms.
To be explicit, we express Eqs. (13)-(15) as where the source term S α is given by Applying the dimensional splitting method to Eq. (39), we obtain the one-dimensional hydrodynamic equations, where P i and V x 0,i are determined by the solution of the Riemann problem whose initial condition is given by Eq. (38). The grid-cell size after the Lagrange step is given by Using the operator L k i which represents one-dimensional evolution in the i direction during the proper time k∆ τ, twodimensional expansion in (x, η) coordinates is given by L k x and L k η , Similarly the three-dimensional expansion in (x, y, η) coordinates is written by

The Courant-Friedrichs-Lewy condition
The Courant-Friedrichs-Lewy (CFL) condition helps us determine an appropriate time-step size in solving partial differential equations. The CFL condition is determined so that the numerical propagating speed on the grid is larger than the physical propagating speed, where ∆ η/∆ τ and w η c are the numerical signal velocity and the characteristic velocity in the rapidity direction, respectively. Note that the proper time τ is multiplied for dimensionless expression. This condition is important for the stability of the numerical calculations. If the CFL condition Eq. (50) is written as a function of the characteristic velocity in z direction v z c , the CFL condition in Milne coordinates has a rapidity dependence. At large rapidity the right-hand side of Eqs. (51) and (52) approaches 1, which is the same condition as that in the case of |v z c | = 1. Therefore the CFL condition in Milne coordinates is determined by Equation (53) indicates that the proper time of the denominator is larger, a possible time-step size ∆ τ is larger.
With these considerations, we define the Courant number as with an initial proper time τ 0 .

Boundary conditions
When we numerically solve partial differential equations and update the value of a cell, the values of its neighboring cells are necessary. For the cell on the boundary of hydrodynamic grid, we need to prepare additional cells which are called ghost cells and put appropriate information to them. In the next section we shall discuss several numerical tests using our code; shock tube problems in Milne coordinates, expansion of matter to the vacuum, Landau-Khalatnikov solution, fluctuations in longitudinal expansion, Gubser flow, and the conservation property. For the shock tube problem in Milne coordinates, we input an analytical solution as a boundary condition at the ghost cells. For the numerical test of expansion of matter to the vacuum, we use the physical values of the vacuum at the cell on the boundary. We employ the periodic boundary condition for the investigation of fluctuations in longitudinal expansion. For the Landau-Khalatnikov solution and the Gubser flow, we copy the values of the cells on the boundary onto those of the ghost cells.
In addition, we point out that in some cases we need a careful procedure at the boundary. For example, for the shock tube problem in Milne coordinates, we use the MC limiter at the boundary and ghost cells to reduce numerical errors which originate from inward flow at the boundary. In the expansion of matter to the vacuum, we observe that a numerical instability occurs at the discontinuity between matter and the vacuum. To stabilize the difficulty we employ the minmod limiter which is dissipative and smears out discontinuities compared to the MC limiter.

Numerical tests
We employ several problems to check correctness of the numerical algorithm in Milne coordinates. For onedimensional tests, we analyze with our code the Riemann problem and Landau-Khalatnikov solution [44,45] in order to verify our Riemann solver. The Landau-Khalatnikov solution is used for understanding the experimental data of the particle rapidity distributions in the high-energy heavy-ion collisions. Next we discuss propagation of fluctuations around Bjorken flow. We derive analytical solutions from linearized hydrodynamics and compare them to numerical calculations with our code. For multidimensional tests, we use the Gubser flow [46,47], which gives us a three-dimensional hydrodynamic expansion of hot and dense matter created after the high-energy heavy-ion collisions. In the test problems we use the ideal gas equation of state, e = 3p.

Shock tubes
Using the property that the initial-value problem of Eq. (27) in Milne coordinates is the same as the Riemann problem in Cartesian coordinates, we carry out the shock tube test in Milne coordinates. We solve the initial-value problem of Eq. (27) with our numerical scheme in Milne coordinates and compare the numerical results with the analytical solution of the Riemann problem Eq. (26). The analytical solution of the Riemann problem in Milne coordinates is obtained by the coordinate transformation from that in Cartesian coordinates [38,39] (see Appendix C).
We perform test calculations in two cases. In the first case we set the discontinuity at η = 0 and in the second case we put it at η = 1 at the initial time τ 0 = 1 fm. In the first test problem, the temperature in η > 0 (η < 0) is set to T L = 400 MeV (T R = 200 MeV) and in both regions the z component of the velocity v z is vanishing. In our numerical scheme, we directly calculate the time evolution of the rapidity component of the velocity w η , instead of v z . The transformation between v z and w η is given by Eq. (9), which suggests that w η is not vanishing even if v z = 0. The initial condition of w η is given by w η = 1 τ 0 −sinhη coshη which has an η dependence. We perform the numerical calculations with grid size ∆ η = 0.01 and time-step size ∆ τ = 0.1τ 0 ∆ η. Figure 3 shows the energy density distribution, the velocity w η , and the velocity v z , which is transformed by Eq. (9) from w η as a function of rapidity η together with the analytical solutions (solid lines). The rarefaction wave moves to the negative direction from η = 0 and the shock wave moves to the positive direction For the second numerical test, we put the initial discontinuity at η = 0. Other conditions, T L , T R and v z are the same as those in the first test calculation. In Cartesian coordinates the first and second numerical tests are essentially identical. In Milne coordinates, however, the second numerical problem is different from the first problem, because the w η depends on η differently from the first problem. The energy density distribution e, the velocity distribution w η , and the velocity distribution v z are shown in Fig. 4. Again, our numerical calculations show good agreement with the analytical solutions in the second numerical test.

Expansion of matter into the vacuum
As one of the specific problems of the Riemann problem, we consider the one-dimensional expansion of matter into the vacuum; a rarefaction wave appears at the discontinuity and expands between the matter and the vacuum [48]. This problem is useful for a realistic description of expansion of the QGP and hadronic matter into the vacuum in the high-energy heavy-ion collisions. We set the initial condition to p = 1000 fm −4 , w η = − sinhη coshη for |η| ≤ 1.5, where w η in |η| ≤ 1.5 corresponds to v z = 0 in Cartesian coordinates and in |η| > 1.5 lies the vacuum. Setting the vacuum for the boundary condition as in Eq. (55), one can avoid the matter from flowing into the system through the boundaries. In the shock tube problems in Milne coordinates in Sect.4.1.1, the matter comes in through the boundaries, which is a possible source of numerical error.
Here we discuss the importance of description of the hydrodynamic equation in the conservative form in developing numerical algorithm. For the investigation, we discretize the hydrodynamic equation with the source terms Eq. (43) and construct a code based on the same procedure explained in Sect. 3.2.1. In the code, T ττ and T τη are updated in the Lagrange step, instead of T τt and T τz . Then we compare the results of the code without the source terms, Eq. (28), and those with the source terms. In both numerical calculations the grid size ∆ η and the time-step size are set to 0.02 and 0.1τ 0 ∆ η, respectively. The hydrodynamic expansion starts at τ 0 = 1 fm.  Figure 5 shows the w η distribution as a function of η at τ = 1.1 fm. In 0 < η < 1.46, the matter is at rest in Cartesian coordinates, which corresponds to the negative flow in Milne coordinates. In 1.46 < η < 1.58, the rarefaction wave starts to expand; it propagates with the sound velocity inward the matter at rest in Cartesian coordinates and expands with the speed of light to the vacuum. In the rarefaction wave the steep velocity gradient is produced, where both of the codes with and without the source terms cannot reproduce the analytical result, though the code without the source terms is closer to the analytical solution. Figure 6 shows the numerical results of the energy density, w η and v z distributions as a function of η together with the analytical solution at later time τ = 4 fm. In |η| < 0.27 the velocity v z is vanishing in Cartesian coordinates, which indicates the negative w η in Milne coordinates. In 0.27 < |η| < 2.89 the rarefaction wave is spreading and at η = ±2.89 the boundary to the vacuum exits, which moves out to the vacuum at the speed of light. For the energy density distribution, both codes with and without the source terms reproduce the analytical solution, though near the boundary between the matter and the vacuum a small difference between them is observed. In both codes, the value of w η around the boundary is larger than that of the analytical solution near the boundary. The stronger flow near the boundary in the numerical solution causes the smaller energy density compared with the analytical solution. On the other hand, there is no difference between the behaviors of v z of both codes. The stability of numerical calculation is sensitive to the differences in w η of codes which are seen near the boundary between matter and the vacuum, because at the boundary pressure becomes zero and the Lorentz factor becomes infinity in the analytical solution.
To investigate the numerical accuracy of the codes with and without the source terms, we calculate the L1 norm which is defined by  where u is the energy density e or the rapidity component of the velocity w η , ∆ η is a grid-cell size. Using Eq. (56), we evaluate the deviation of the numerical results u(η i ; N cell ) from the exact solutions. At the same time we can know the convergence speed to the exact solution of numerical algorithm. In Fig. 7, the L1 norms of energy density and w η are shown. In both energy density and velocity, the values of the L1 norm of the code without the source terms are smaller than those with the source terms, which means that the code without the source term has smaller numerical viscosity than that with the source terms. As expected, the existence of the source terms produces more artificial viscosity.
If the initial discontinuity is set at the larger rapidity, this makes the velocity slope at rarefaction wave larger and gives more severe problems to the codes in Milne coordinates. We find that the code without the source terms is more stable than that with the source terms. For example, if we set the initial discontinuity at η = 1.7, we find that numerical instability occurs in the code with the source terms.

Landau-Khalatnikov solution
We employ the Landau-Khalatnikov solution [44,45] as a one-dimensional numerical test problem. The initial condition of it is expressed by a thin slab of hot and dense matter created after the collisions, which is the same as the problem discussed in Sect. 4.1.2. In the expansion of the slab of matter, two rarefaction waves travel into the slab from both sides and start to overlap at the center of the slab. The region where rarefaction waves overlap is described by the Landau-Khalatnikov solution. The asymptotic form of the Landau-Khalatnikov solution for sufficiently later time τ ≫ ∆ is written by and w η = 0, where ∆ is the thickness of the slab. The asymptotic solution Eq. (57) is used for an investigation of the rapidity distributions of the produced particles at RHIC [49][50][51][52][53][54].
In the numerical calculation, we start the simulation at τ 0 = 500 fm with the initial condition given by Eq. (57), where e 0 and the thickness size are set to 10 GeV/fm 3 and ∆ = 0.5 fm, respectively. The numerical calculation is performed with the grid size ∆ η = 0.1, the time-step size ∆ τ = 0.1 fm and ∆ τ = 0.1τ 0 ∆ η = 5 fm, which is determined by the CFL condition in Sect. 3.3. Figure 8 shows the energy density distributions at τ = 510, 600, 700, and 1000 fm together with the analytical solution. Calculations with the time-step sizes ∆ τ = 5 fm and 0.1 fm can explain the analytical solution, which suggests that the computational time can be saved if the time-step size is determined by the CFL condition in Sect.3.3. There is a small deviation between numerical calculations and the analytical solution at large |η|, which implies that the asymptotic form of the Landau-Khalatnikov solution Eq. (57) cannot be applicable at large rapidity [45].

Propagation of longitudinal fluctuations around Bjorken flow
The longitudinal fluctuation in particle distributions and collective flows is one of the interesting topics in high-energy heavy-ion collisions at LHC [55]. For instance, the propagation of fluctuations around Bjorken flow in heavy-ion collisions are investigated from linear analyses [7,56]. The small longitudinal fluctuations around Bjorken flow propagate according to the following linearized equations: Here e B = e 0 τ 0 τ 1+λ is the energy density from Bjorken's scaling solution [57] for the equation of state p = λ e 3 . Since the background is rapidity independent, we can obtain solutions with a definite wave number δ e, δ w η ∝ e ikη . The solutions consist of two modes, as Eq. (60)  Here we compare the analytical solutions and the numerical calculation with our hydrodynamic code for two cases, D > 0 and D < 0. In both cases, we choose initial conditions so that we can single out a particular mode of attenuation (D > 0) and propagation (D < 0). To be specific, for D > 0 we choose where θ is defined by θ ≡ 1 2 √ −Dlog(τ/τ 0 ). The phase velocity of the fluctuation is √ −D/(2kτ). We set the initial time to τ 0 = 1 fm and use the ideal gas equation of state, e = 3p and λ = 1/3. The initial energy density e 0 in Bjorken flow is given by e 0 = 1000 fm −4 . In the case of D > 0, we choose k = 0.5, D = 0.111 and A = 0.1 fm −4 . In the case of D < 0, we choose k = 2π, D = −52.193 and A = 0.1 fm −4 . Figures 9 and 10 show the analytical and numerical results of the fluctuations of the energy density and the velocity around Bjorken flow for D > 0 and D < 0. In the numerical calculation, we use the periodic boundary condition. When the value of D is positive, the fluctuation does not propagate and its amplitude decreases with time ( Fig. 9). On the other hand, when the value of D is negative, the fluctuation propagates and its amplitude decreases with time (Fig. 10). However, if we chose different initial conditions such that more than one mode were involved, the amplitude of the fluctuation would at first grow and then reduce for both D > 0 and D < 0 due to the interference of two modes. The similar amplification of the fluctuations around Bjorken flow is reported in Ref. [56]. Our numerical results show good agreement with the analytical solutions.

Gubser flow
An analytic solution to the relativistic, conformally invariant Navier-Stokes equation is constructed based on symmetry considerations in (τ, η, x ⊥ , φ ) coordinate system [46,47]. The Gubser flow and its related solutions are utilized for checking or improvement of hydrodynamic codes [12,[58][59][60][61][62][63]. The solution is a generalization of   Bjorken flow where the medium expands both longitudinally and radially, which gives us a realistic description of the space-time evolution of high-energy heavy-ion collisions.
According to the solution for inviscid fluid, the transverse or radial velocity w ⊥ and energy density e are given by whereê 0 is a dimensionless integration constant, q is an arbitrary dimensional constant with unit of inverse length of the system size [46,47]. We compare our numerical calculations and the analytical solution in Figs.11 and 12. In our numerical calculation, the parameters are set to q = 1 fm −1 andê 0 = 400. The hydrodynamic expansion starts at τ 0 = 1 fm. The numerical simulation is performed with the grid size ∆ x = ∆ y = 0.05 fm, ∆ η = 0.1 and the time-step size ∆ τ = 0.1∆ x. First the consistency between our numerical calculations and the analytical solution suggests the Strange splitting method in Sect. 3.2.2 works correctly. In the Gubser flow, the existence of the initial transverse flow is an origin of the strong transverse flow at later time. To reproduce the strong transverse flow in the Gubser flow, we find that careful choice of the interpolation method is required in the numerical calculation. For example, if the second-order interpolation method is employed, the energy density around x = 0 fm of the numerical calculations is larger than that of the analytical solutions. However, if we use the PPM interpolation method, we can reproduce the analytical solutions numerically. It is discussed in Appendix E.

Conservation property
We check the energy and momentum conservation of our code. The conserved quantities in our algorithm in Milne coordinates are given by τT τν (ν = t, x, y, z). Their time evolution in our algorithm is schematically written by where F ν i represents the flux of the conserved quantities which flow into and out of i-th grid cell during ∆ τ. Equation (67) contains two steps: Lagrange and remap steps. Integrating Eq. (67) on all spatial grids, we obtain which suggests that the total variation of the conserved quantities depends on the amount of inflow and outflow from the boundary. If the equations with the source terms, Eqs. (16)- (19) are used in the code, the right-hand side of  Table 1 The violation of the total energy and momentum conservation.
Eq. (68) has an additional term from the source terms, which can spoil the conservation property and affects the numerical accuracy in application to physical problems.
Here we focus on the effects of existence of the source terms in Milne coordinates on conservation property. We perform our numerical calculation with the initial energy density and flow distributions which are usually used in study of the relativistic heavy-ion collisions, where Y b = 5.3 is the beam rapidity, σ η = 2.1 and η flat = 2.6 show the size of the flat structure of the initial energy density distribution in the rapidity, and e 0 = 30 GeV/fm 3 is the maximum value of the energy density. We choose a typical parameter set which is tuned for the RHIC collision energy [60,64]. To discuss the source-term effect in Milne coordinates clearly, we carry out the numerical calculation only in the rapidity direction. In Ref. [60] one calculates the total energy and entropy in the beginning and at the end of 3D hydrodynamic evolution, using the Glauber model with a limited rapidity profile for an initial condition. From comparison between them they find that the energy is conserved on a level of better than 3 % in their code.
To check the conservation property of our code, we evaluate the total energy and momentum in the hydrodynamic expansion, We carry out a numerical calculation from τ 0 = 1 fm to τ = 10 fm on ∆ η = 0.02, 0.1, 0.2, and 0.5 grid sizes with ∆ τ = 0.1τ 0 ∆ η time-step size. The total momentum in the beginning is vanishing in the numerical precision. We observe that in the case where the source terms are explicitly included, the total energy and momentum increase with the proper time τ monotonically on ∆ η = 0.5. On other grid sizes, first the total energy and momentum increase with τ and after some time steps they start to decrease. If there are no source terms, sometime the total energy and momentum increase and at other times they decrease. This behavior suggests that the simple comparison between total conserved quantities in the beginning and those at the end is not suitable for an investigation of the conservation property. Instead of the simple comparison, we evaluate the violation of the energy and momentum conservation at each time step and sum it from the beginning to the end, We show the calculated results of them in Table 1. We find that the numerical calculation based on the equations without the source terms keeps the energy and momentum conservation with high accuracy compared with that with source terms on every grid size. In the case where the source terms are explicitly included, enough numerical accuracy is still kept, but the amount of the violation of the conservation property increases with grid size. On the other hand, in the numerical algorithm with the conservative form, it does not depend on the grid size. Even on the the course grid, the conservation property is kept with very high accuracy. Next we investigate the total energy and momentum conservation in hydrodynamic evolution which starts from a fluctuating initial condition. We add the fluctuation to the energy density and velocity distributions in Eqs. (69) and (70), where e flat is given by Eq. (69) and values of η e n and η v n are chosen between η = −Y b and Y b at random. We set δ e n = 0.05 and δ w η = 0.05 fm −1 for all n. We carry out numerical calculations with the grid size ∆ η = 0.2, which is often chosen in calculations of high-energy heavy-ion collisions [60]. We set the time-step size equal to ∆ τ = 0.1τ 0 ∆ η. In Fig. 13 the energy density and velocity distributions at τ = 10 fm with and without the source terms are shown. At the mid-rapidity the energy density and the flow distributions of numerical calculations with and without the source terms are consistent with each other. In the region of |η| > 4, however, the differences between them are observed in the small structure of both distributions. The growth of the velocity to the vacuum |η| ∼ 8 gives a difficulty of numerical calculation and becomes the reasons for the differences. The deviation from the energy and momentum conservation are listed in Table 2. For both cases, we find that they are around ten times as large as those with the smoothed initial condition. Nevertheless, the code based on the conservative form keeps conservation property with high accuracy. On the other hand, in the code with the source terms a few % deviation from the energy and momentum conservation appears, which is still acceptable. In the code with the source terms, numerical calculation with fine grid size is indispensable for the energy and momentum conservation. There exist other ingredients which can cause additional error, and violation of the conservation property originates  Table 2 The violation of the total energy and momentum conservation with fluctuating initial conditions. The initial total energy and momentum are E 0 = 2224 GeV and M 0 = −94 GeV.
from the geometric source term, for instance shock waves and jets in medium [65][66][67]. In addition, the existence of the viscosity can be the origin of the breakdown of the conservation property [60]. To avoid such problems, we need to construct the codes based on constitutive equations with the conservative form or perform numerical calculations on sufficiently fine grids.

Summary
We constructed a new Godunov type relativistic hydrodynamic code in Milne coordinates based on the algorithm in Cartesian coordinates [23]. We evaluated the flux terms, using the numerical solution of the Riemann problem with the initial condition at the constant proper time τ. We checked the correctness of our algorithm from the comparison between numerical calculations and analytical solutions of shock tube, expansion of matter into the vacuum, the Landau-Khalatnikov solution, propagation of fluctuation around Bjorken flow and the Gubser flow. We investigated the energy and momentum conservation of our code from a calculation of the longitudinal hydrodynamic expansion with an initial condition for high-energy heavy-ion collisions.
In particular, we focused on the effects of the source terms in relativistic numerical hydrodynamics in Milne coordinates on stability and numerical viscosity. We analyzed those effects in the test problems of expansion into the vacuum and the conservation property. In expansion of matter into the vacuum, we showed that numerical results from the code without the source terms are closer to the analytical solution compared with that with source terms. Besides, the code without the source terms is more stable and has less numerical viscosity than the code with the source terms. In addition, we observed that the code written in the conservative form keeps the conservation property with high accuracy in the expansion from the fluctuating initial longitudinal profile for high-energy heavy-ion collisions, even on a coarse grid.
Our algorithm is easily extended to the code with the QCD equation of state and finite viscosities [23,24]. After that, we shall employ our hydrodynamic code to investigate experimental results at RHIC and LHC and understand the detailed QGP bulk property using a reliable 3D relativistic viscous hydrodynamic expansion with small numerical viscosity.

Acknowledgments
The work of CN is supported by the JSPS Grant-in-Aid for Scientific Research (S) No. 26220707 and US Department of Energy Grant DE-FG02-05ER41367. Y. A. is partially supported by JSPS Postdoctoral Fellowships for Research Abroad.

Appendix A: Riemann problem in Milne coordinates
We show that the two initial problems Eqs. (26) and (27) have the same analytic solution. First, we prove that the hydrodynamic state which satisfies the condition ∂ η V V V = 0 at some proper time keeps ∂ τ V V V = 0. In other words, if the hydrodynamic state V V V is uniform at some initial proper time, it remains the initial uniform state at all times. This state corresponds to the initial condition of the Riemann problem in Milne coordinates Eq. (27). For the (1 + 1)-dimensional case, we rewrite the energy conservations, Eqs. (24) and (25), as Inserting the conditions, ∂ η V V V = 0, namely ∂ η T tt = 0, ∂ η T tz = 0 and ∂ η T zz = 0, into Eqs.  27), the hydrodynamic state outside the light cone of the discontinuity (τ 0 , η i ) remains V V V L or V V V R , based on the fact that the signal velocity in the ideal hydrodynamic equation is smaller than the speed of light. In Fig. 14 (b), the initial condition of the Riemann problem is set on the hyperbolic curve τ = τ 0 with the discontinuity at η i , the hydrodynamic state outside the light cone of the discontinuity is given by V V V L or V V V R . Figure  14 (a) shows the initial condition of the Riemann problem in Cartesian coordinates with the same discontinuity point at (t 0 , z i ), which corresponds to (τ 0 , η i ) in Fig. 14 (b). From a comparison between Figs. 14 (a) and (b), the analytical solutions for initial-condition problems which are given by Eqs. (26) and (27) satisfy the same boundary conditions on the light cone at the discontinuity (τ 0 , η i ). We conclude

Appendix B: Discretized hydrodynamic equation with Lagrange step
We show the detailed calculation of numerical flux terms Eqs. (35)-(37). We represent the velocity and pressure of the grid-cell boundary at η i as V z i and P i , respectively. First we shift the grid-cell boundary at η i to η = 0 by the Lorentz boost transformation to derive the moving distance δ η i (τ ′ ) of the grid-cell boundary during τ ′ (0 ≤ τ ′ ≤ ∆ τ). Under the Lorentz boost transformation, δ η i (τ ′ ) is invariant. The velocity of the grid-cell boundary V z 0,i after the Lorentz boost transformation is related to V z i by The moving distance δ η i (τ ′ ) of the grid-cell boundary from η = 0 during τ ′ is given by Eq. (33) with Eq. (34). Equation (34) is derived from the simultaneous equations z = V z 0,i (t − τ n ) and (τ n + τ ′ ) 2 = t 2 − z 2 with t = τ n + δt(τ ′ ).
the calculation, we use the same parameters as those in Sect. 4.4. The numerical calculations with the third-order interpolation procedure reproduce the analytical solutions. On the other hand, the energy density with the second-order interpolation procedure is slightly larger than that of the analytical solutions. In particular, we observe the deviation from the analytical solution of energy density in |x| < 2 fm (|x| < 5 fm) at τ = 2 (τ = 5 fm). We find the same behavior in other second-order interpolation procedures, minmod and superbee limiters. However, in the case of the one-dimensional expansion, even if the strong expansion exists, the numerical calculation with the second-order interpolation procedures shows good agreement with the analytical solution. We observe the deviation between the analytical solution and the numeral results in multidimensional calculation, which suggests that the operator splitting method is also a possible key issue for the problem. In Fig. 16, we show the transverse velocities at τ = 2 and τ = 5 fm. The gradient of the transverse velocity increases rapidly up to x ∼ 2 (x ∼ 5) fm at τ = 2 (τ = 5) fm, where the value of the transverse velocity is slightly smaller than that of the analytical solutions, which implies that the secondorder interpolation schemes do not satisfy the description of such a rapid expansion. The inadequate velocity growth causes the delay of the decrease of the energy density which is observed in Fig. 15.
In Ref. [61] one points out the importance of adjusting the flux limiter in the algorithm (KT algorithm), using the relativistic viscous hydrodynamics. One shows that a free parameter ξ in the van Leer minmod filter is fixed from a comparison between the solutions of the shear-stress tensor from Gubser flow and the numerical calculations.