A new relativistic viscous hydrodynamics code and its application to the Kelvin-Helmholtz instability in high-energy heavy-ion collisions

We construct a new relativistic viscous hydrodynamics code optimized in the Milne coordinates. We split the conservation equations into an ideal part and a viscous part, using the Strang spitting method. In the code a Riemann solver based on the two-shock approximation is utilized for the ideal part and the Piecewise Exact Solution (PES) method is applied for the viscous part. We check the validity of our numerical calculations by comparing analytical solutions, the viscous Bjorken's flow and the Israel-Stewart theory in Gubser flow regime. Using the code, we discuss possible development of the Kelvin-Helmholtz instability in high-energy heavy-ion collisions.


Introduction
Since the success of production of the strongly interacting quark-gluon plasma (QGP) at Relativistic Heavy Ion Collider (RHIC) [1], relativistic viscous hydrodynamic model has been one of promising phenomenological models. Now at RHIC as well as at the Large Hadron Collider (LHC) high-energy heavy-ion collisions are carried out. The strong collective dynamics observed in experimental data at RHIC and the LHC provides us with a clue of understanding the QCD matter. A relativistic hydrodynamic model is suitable for description of space-time evolution of strongly interacting QCD matter produced after collisions. Besides, it has a close relation to an equation of state and transport coefficients of the QGP. The QCD phase transition a e-mail: okamoto@hken.phys.nagoya-u.ac.jp b e-mail: nonaka@hken.phys.nagoya-u.ac.jp mechanism and the QGP bulk property is elucidated from comparison between hydrodynamic calculation and experimental data.
A relativistic viscous hydrodynamic model plays an important role in the quantitative understanding of the QGP bulk property. However, introducing viscosity effect into the framework of relativistic hydrodynamics is not an easy task, because of the acausality problem. There is not the unique way to extract the second-order relativistic viscous hydrodynamic equation. In high-energy heavy-ion collisions, currently the Israel-Stewart theory [2,3] and conformal hydrodynamics [4] are often used. Solving them numerically, study of experimental data of high-energy heavy-ion collisions is performed [5][6][7][8][9][10][11][12][13][14].
Now the relativistic viscous hydrodynamic model can explain not only the elliptic flow but also higher harmonics [15]. In particular, analyses of the higher harmonics bring us progress of understanding of the QGP, because it is more sensitive to the QGP bulk property. Furthermore, a lot of experimental data are reported; correlation between flow harmonics [16,17], event plane correlation [18,19], non-linearity of higher flow harmonics [20] and three particle correlation [21,22]. At the same time, we can investigate the QGP property further using information of (3+1)-dimensional space-time evolutions [19,21,22]. The rich experimental data realizes investigation of both shear and bulk viscosities and even their temperature dependence.
We need to perform numerical calculations for relativistic viscous hydrodynamics with high accuracy, to achieve the quantitative analyses of the transport coefficients of the QGP from comparison with high statistics and high precision experimental data. For example, the following features in numerical calculations are demanded: A fluctuating initial condition is correctly captured and numerical viscosity which is needed for stability of calculation is much smaller than physical viscosity. Furthermore, time evolution of the viscous stress tensor is sensitive to numerical scheme, because it consists of time and space derivatives of hydrodynamic variables.
Here we present a new relativistic viscous hydrodynamics code optimized in the Milne coordinates. The code is developed based on our algorithm of the ideal fluid in which a Riemann solver with the two shock approximation [23] is employed [24]. It is stable even with small numerical viscosity [25]. We shall show comparison between numerical calculations and analytic solutions of viscous Bjorken's flow and the Israel-Stewart theory in Gubser flow regime.
Using the code, we shall discuss possible development of Kelvin-Helmholtz (KH) instability in high-energy heavyion collisions. Hydrodynamic instability and turbulent flow are discussed in Ref. [26,27] and the possibility of KH instability is argued in Ref. [28]. The hydrodynamic instability is affected by a viscosity effect, which suggests that the numerical code with less numerical viscosity is indispensable for study of it.
This paper is organized as follows. We begin in Sect. 2 by showing the relativistic viscous hydrodynamic equations briefly. In Sect. 3 we explain the numerical algorithm; Strang splitting method and numerical implementation. We check the validity of our code comparing analytic solutions of viscous Bjorken flow and the Israel-Stewart theory in the Gubser flow regime in Sect. 4. In Sect. 5, we discuss the possible development of KH instability in high-energy heavy-ion collisions. We end in Sect. 6 with our conclusions.

Relativistic viscous hydrodynamic equations
The relativistic hydrodynamics is based on the conservation equations, where N µ is the net charge current and T µν is the energymomentum tensor. In the case of ideal fluid, the net charge current and energy-momentum tensor are given by where n is the net charge density, e is the energy density, p is the pressure and u µ is the fluid four-velocity which satisfies the normalization u µ u µ = 1. ∆ µν is the orthogonal projection tensor to u µ , which is defined by with the metric tensor g µν . Here the u µ is determined uniquely.
On the other hand, in dissipative flow, there are several possible choices to determine u µ . For example, one can assign the u µ as net charge flow (Eckart frame [29]) or as energy flow (Landau frame [30]). The decomposition of N µ and T µν in viscous fluid depends on the choice of u µ . Here we choose the Landau frame for relativistic viscous hydrodynamic equations, because we focus on the high-energy heavy-ion collisions as RHIC and the LHC where the net baryon number is very small [31].
In the Landau frame, the net charge current and the energy-momentum tensor of the viscous fluid are decomposed as where n µ is the charge diffusion current, Π is the bulk pressure, and π µν is the shear tensor [30]. The relativistic extension of Navier-Stokes theory in non-relativistic fluid usually has a problem of acausality and instability [32][33][34].
The problem can be resolved by introducing the secondorder terms of the viscous tensor and the derivative of fluid variables into the hydrodynamic equations [2,3]. However, the original Israel-Stewart theory does not reproduce the results of the kinetic equation quantitatively [35][36][37][38][39]. The construction of second-order relativistic viscous hydrodynamic equations is still under investigation. The extension of the Israel-Stewart theory is also proposed [40][41][42][43][44][45]. In addition to the framework of the Israel-Stewart theory, other approaches such as the AdS/CFT correspondence [4,[46][47][48] and renormalization group method are applied to the construction of causal relativistic hydrodynamics [49,50].
In the second-order viscous hydrodynamics, additional equations for evolution of the viscous tensors are needed. Here, we introduce the convective time derivative D and the spatial gradient operator ∇ µ , which are defined by respectively. For example, in the second-order Israel-Stewart formalism the constitutive equations of the viscous tensors are given by where τ n , τ π , and τ Π are relaxation times, I µ n , I µν π , and I Π represent second-order terms. n µ NS , π µν NS , and Π NS are the Navier-Stokes value of viscous tensors written as where T is the temperature, µ is the chemical potential, θ ≡ u µ ;µ is the expansion scalar, σ is the charge conductivity, η is the shear viscosity, and Π is the bulk viscosity.
We construct a relativistic viscous hydrodynamics code in the Milne coordinates (τ, x, y, η) which is optimized for description of the strong longitudinal expansion [51] at RHIC and the LHC. In the Milne coordinates, the metric tensor is given by g µν = diag(1, −1, −1, −1/τ 2 ) and the fluid four-velocity has the form where i = x, y, η and the right-hand sides of them represent geometric source terms. S ν is given by The constitutive equations Eqs. (10), (11), and (12) in Milne coordinates read where τ n , τ η and τ Π are the relaxation times, and the secondorder terms are defined by J ττ π = 2τv η π τη , J τ j π = τv η π jη , j = x, y and λ = τ, x, y, η. Here, J τ n and J µν π are the geometric source terms which come from the convective time derivative of n µ and π µν respectively. K µ n and K µν π ensures the constraints n µ u µ = 0, π µν u µ = 0 and π µ µ = 0.

Numerical algorithm
In this section, we present our numerical algorithm for solving the relativistic viscous hydrodynamic equations in the Milne coordinates.

Strang splitting method
In our algorithm, we split the conservation equations Eqs. (16) and (17) into two parts, an ideal part and a viscous part using the Strang splitting method [52]. Specifically, the net charge current and the energy-momentum tensor are divided as follows: vis ≡ π µν − Π ∆ µν . The subscripts "id" and "vis" mean the ideal part and the viscous part, respectively. The equations of the ideal part are expressed by where S ν id = −T ττ id /τ − τT ηη id , −T τx id /τ, −T τy id /τ, −3T τη id /τ . They are nothing but usual ideal hydrodynamic equations in the Milne coordinates. On the other hand, the equations of the viscous part are given by where S ν vis = −T ττ vis /τ − τT ηη vis , −T τx vis /τ, −T τy vis /τ, −3T τη vis /τ . They give viscous corrections to the evolution of the ideal fluid.
The Strang splitting technique is also applied to evaluate the constitutive equations of the viscous tensors Eqs. (19)- (21). We decompose the constitutive equations into the following three parts; the convection equations, the relaxation equations, and the equations with source terms, In numerical simulation of relativistic hydrodynamic equation, a time-step size ∆ τ is usually determined by the Courant-Friedrichs-Lewy (CFL) condition. However, in the relativistic dissipative hydrodynamics, one needs to determine the value of ∆ τ carefully. The relaxation times τ n , τ η , and τ Π in the constitutive equations show the characteristic timescale of evolutions of the viscous tensors, which means that a small relaxation time gives us more restrictive condition to ∆ τ than the CFL condition does. If the relaxation times τ n , τ η , and τ Π are much shorter than the fluid timescale τ fluid , the time-step size ∆ τ should be smaller than the relaxation timescale, which makes the computational cost increase. To avoid this problem, we use the Piecewise Exact Solution (PES) method [53], instead of using a simple explicit scheme. In the PES method, formal solutions of Eqs. (37)-(39), can be used. On the other hands, if the relaxation times are larger than ∆ τ determined by the CFL condition, the PES method is not applied [12]. In our algorithm, we solve the time evolution of n x , n y , n η , π xx , π yy , π ηη , π xy , π yη , π ηx and Π directly. Other components of viscous tensors n τ , π ττ , π τx , π τy and π τη are derived from the orthogonality conditions n µ u µ = 0 and π µν u ν = 0.

Numerical implementation
The decomposed hydrodynamic equations Eqs. (30)- (42) are solved by the following procedure. Here, we represent a conserved variable as First, we solve the ideal part of the conservation equations Eqs. (30) and (31) using the Riemann solver [24]. In this step, the conserved variable U U U id (τ) is evolved into U U U * id (τ + ∆ τ), where the asterisk indicates a variable evolved only in the ideal part. V V V id (τ) is used to evaluate the numerical flux and the geometric source terms in Eqs. (30) and (31). We calculate the fluid variable V V V * id (τ + ∆ τ) from U U U * id (τ + ∆ τ) with the algorithm for recovery of the primitive variables V V V id from the conserved variables U U U id [25].
Here we keep the middle time-step value of the viscous tensor V V V vis (τ + ∆ τ/2) for the next step.
Next, the conserved variables U U U * id (τ + ∆ τ) and U U U vis (τ) are evolved into U U U id (τ + ∆ τ) and U U U vis (τ + ∆ τ) by the viscous part of conservation equation Eqs. (32) and (33). Then we recover the fluid variables V V V id (τ + ∆ τ) from conserved variables U U U vis (τ + ∆ τ) [53]. We keep the middle To achieve the second-order accurate in time, we repeat the above whole steps using the middle time-step values V V V id (τ + ∆ τ/2) and V V V vis (τ + ∆ τ/2). However, we find that numerical errors arise mainly from the constitutive equations Eqs. (34)- (42). Therefore we carry out numerical calculation in the second-order accurate in time only in constitutive equations Eqs. (34)- (42) and the viscous part of conservation equations Eqs. (32) and (33).
Throughout all above steps, we evaluate space derivative terms using the MC limiter [54] for the second-order accurate in space or the piecewise parabolic method (PPM) [55][56][57] for the third-order accurate in space. We shall give the explicit expressions of the interpolation procedures, the MC limiter and the PPM in Appendix A.

Numerical tests
We check the correctness of our code in the following test problems; the viscous Bjorken flow for one-dimensional expansion and the Israel-Stewart theory in Gubser flow regime [58] for the three-dimensional calculation. We use the ideal massless gas equation of state, p = e/3 and set the net charge to be vanishing.

Viscous Bjorken flow
The Bjorken flow is one of the simplest one-dimensional test problems for the code which is optimized in the Milne coordinates. In the ideal fluid, the time evolution of temperature follows T = T 0 (τ 0 /τ) 1/3 , where τ 0 and T 0 are the initial proper time and temperature, respectively [51]. In the viscous fluid, the non-vanishing components of viscous tensor π µν are π xx , π yy , and π ηη . From the symmetries of the system, the relation 2π xx = 2π yy = −τ 2 π ηη holds. First, we focus on the shear viscosity effects at the Navier-Stokes limit. In the Navier-Stokes limit, the relativistic viscous hydrodynamic equation with the boost invariance is written as where π ηη NS is the Navier-Stokes value of shear tensor, If η/s is constant, Eqs. (46) and (47) give the time evolution of the temperature, The numerical calculation is carried out on the spacegrid size ∆ η = 0.1 with the time-step size ∆ τ = 0.1τ 0 ∆ η. The initial temperature T 0 and the proper time τ 0 are set to T 0 = 300 MeV and τ 0 = 1 fm, respectively. We set the relaxation time to be τ η = 0.0001 fm as the Navier-Stokes limit. Since τ η is smaller than ∆ τ, the PES method is applied to solve the relaxation equations Eqs. (37)- (39). Figure 1 shows the analytical and numerical results of the Bjorken flow with and without shear viscosity. In the case of finite shear viscosity, the temperature decreases with proper time more slowly, compared to that of the ideal fluid. In both cases, our numerical results show good agreement with the analytical solutions.
Next, we check the time evolution of the bulk pressure in the viscous Bjorken's flow. Ignoring the second-order terms I Π in Eq. (21), we write the relaxation equation of the bulk pressure, with the Navier-Stokes value of the bulk pressure Π NS = ζ /τ. If we assume ζ and τ Π are constant, we obtain the analytical solution of Eq. (49), where Π 0 is the initial value of the bulk pressure and Ei(x) is the exponential integral function.
In the numerical calculation, we set Π 0 = 0, ζ = 1000 MeV/fm 2 , τ 0 = 1 fm, and τ Π = 1 fm. The space-grid size ∆ η = 0.1 and the time-step size ∆ τ = 0.1τ 0 ∆ η are utilized. Figure 2 shows the analytical and numerical results of the time evolution of the bulk pressure in the Bjorken flow. Our numerical calculation is consistent with the analytical solution.

Israel-Stewart theory in the Gubser flow regime
Based on the symmetry arguments developed by Gubser [59,60], a semi-analytic solution of the Israel-Stewart theory in the Gubser flow regime is obtained [58]. The semi-analytic solution is a useful test problem for the code of relativistic viscous hydrodynamics which is developed for application to the high-energy heavy-ion collisions [13,14,58,61,62]. The velocity profile of the semi-analytic solution is the same as that of the ideal Gubser flow, where q is an arbitrary dimensional constant with unit of inverse length of the system size and set to q = 1 in comparison with numerical computation. The solutions of the temperature and the shear tensors are derived by solving a set of two ordinary differential equations numerically [58]. The second-order terms and the relaxation time in Eq. (20) are given by where c is a constant [58]. We carry out the numerical calculation with the finite shear viscosity η/s = 0.2. We set the relaxation time to   [24]. On the other hand, in the finite viscosity calculation the difference between the numerical calculation and the semi-analytic solution appears after τ = 4 fm.
In Fig. 4 the numerical results of the shear tensors π xx , π yy , π ηη and π xy at τ = 1.2, 2 and 3 fm are presented together with the semi-analytic solutions. Here the profile of π xy is shown along a line x = y, since the value of π xy vanishes on the x and y axes. The shear tensors π xx , π yy and π ηη in our numerical calculations show good agreement with the semi-analytic solutions. However, in π xy the deviation from the semi-analytic solution starts to appear at τ = 2 fm and grows at later time.
Since in the Israel-Stewart theory the second-order terms in π µν become small compared with the first-order terms, choice of numerical scheme for evaluation of the convection term in Eq. (35) is important. For example, in Ref. [58], they show that adjustment of the flux limiter which controls possible artificial oscillation in a higher order discretization scheme is crucial for good agreement with the semi-analytic solution. Here we employ the PPM for solving the convection part numerically, instead of the MC limiter. In the case of three-dimensional calculation we use the dimensional splitting method [24]. We find that the Corner Transport Upwind (CTU) scheme [63] which is a three-dimensional unsplit method, realizes good agreement of the semi-analytic solution even with the MC limiter.
We discuss the numerical scheme dependence on the shear tensors in solving the convection term in Eq. (35). We compare the three numerical schemes; a dimensional splitting method with the MC limiter, a dimensional splitting method with the PPM and the CTU method [63] with the MC limiter for three-dimensional unsplit method. We shall explain the details of each scheme in Appendix A and Appendix B. Figure 5 shows the semi-analytic solutions and numerical results of the shear tensors π xx , π yy , τ 2 π ηη and π xy at τ = 3 fm. In π xx , π yy and τ 2 π ηη , results of all numerical schemes are reasonably consistent with the semi-analytic solutions. In addition, differences among them are small. However, in π xy we can clearly see the scheme difference. In the solution obtained with the MC limiter, the large deviation from the semi-analytic solution at the peak around x = 2 fm appears, whereas the PPM and the CTU methods keep the good agreement with the semi-analytic solution. The CTU method can achieve the high numerical accuracy with the second-order accurate in space, but it needs the more computer memory than the dimensional splitting method with the PPM does. Therefore we employ the dimensional splitting method with the PPM for solving the convection term.

Kelvin-Helmholtz instability in Bjorken expansion
We discuss the possible development of the KH instability in relativistic heavy-ion collisions. The KH instability is one of the hydrodynamic instabilities. It occurs on the interface between two horizontal streams which have different velocities [64]. If it takes place, perturbations to the interface between fluids grow and result in vortex formation. In heavy-ion collisions, the color-flux tube structure in initial condition can be an origin of the KH instability; fluctuations in the longitudinal direction are amplified with the KH instability, however, vortex formation is not observed [28]. Recently initial fluctuations and QGP expansion not only in the transverse direction but also in the longitudinal direction have attracted interest [19,21,22]. Using the new relativistic viscous hydrodynamics code which has small numerical viscosity, we investigate the KH instability and vortex formation in heavy-ion collisions.
For simplicity, we focus on hydrodynamic expansion in the (x, η) plane. The heavy ion accelerated with high-energy still has about 1 fm width in the longitudinal direction (z direction) due to the uncertainty principle. In other words, a thin disk composed of large-x partons is covered by a cloud of small-x partons. As a result, in the high-energy heavy-ion collisions parton-parton interactions may take place in the area within around 1 fm from z = 0 fm. Then if we consider the color-flux tube structure in the initial condition, each color-flux tube may evolve from a different interaction point in |z| < 1 fm.
Suppose that two initial flow fluxes are located in x > 0 and x < 0 which represent two color-flux tubes starting to expand at z = ∆ z and z = −∆ z, respectively. Energy density and η component of velocity of the flow flux are assumed to be described by Bjorken where τ 0 and e 0 are the initial time and the energy density, respectively. Figure where the energy density and flow velocity around the boundary are connected from e U and v η U to e D and v η D smoothly with the parameter ∆ . Here, we set the wavelength of a fluctuation and the width of boundary between two fluid fluxes to λ = 0.4 and ∆ = 0.02 fm, respectively. Focusing the hot spots in a fluctuating initial condition at the LHC, we fix the initial energy density (temperature) to e 0 = 741 GeV/fm 3 (T 0 = 800 MeV). Figure 7 shows the velocity field and profile of the vorticity w y of the initial condition Eqs. (58) and (59). Here the definition of the vorticity w y , The arrows stand for the velocity field in (τv η , v x ). We start the numerical calculation on the grid (∆ x, ∆ η) = (0.005 fm, 0.00625) at τ 0 = 1 fm with time-step size ∆ τ = 0.2∆ x. We use the ideal gas equation of state e = 3p.
First we argue on the KH instability in the ideal fluid. We find a starting vortex formed around the boundary at τ ∼ 3 fm. In Fig. 8 the velocity field and the profile of the vorticity w y at τ =4 and 7 fm are shown. We observe that the boundary with the two vortexes tilt toward negative x. The initial conditions Eq. (58) and (59) and Fig. 6 suggest that e U is larger than e D and |v z | in x > 0 is larger than that in x > 0. The e U decreases more slowly than e D does due to the time dilation from larger |v z |. The energy density and the flow differences between x > 0 and x < 0 cause the flow in the negative x direction. The two vortices expand with time and their sizes grow because of existence of the Bjorken flow. As a result, the intensity of the vortices becomes small. The larger the difference of velocity in the shear flow is, the faster the growth of instability is. That is why the development of vortex at η ∼ 0.6 is faster than that at η ∼ 0.2. The fluctuation with a longer wavelength grows slower in the KH instability than that with a shorter wavelength does. If we set the wavelength λ to λ > 0.5 in the region |η| < 0.8, the growth of fluctuation is too slow to form the vortex and the fluctuation is smeared with the Bjorken flow. However, at the forward rapidity, a fluctuation with a long wavelength can survive to form a vortex.
Next we discuss the KH instability with finite viscosity. We employ the same values of the second-order term and the relaxation time in Eq. (20) as those in Sec. 4.2. The shear viscosity is set to η/s = 0.01. Figure 9 shows the numerical results of KH instability at τ = 4 and 7 fm. In contrast to Fig. 8, we cannot find the clear vortex but small and vague enhancement of vorticity around η ∼ 0.2 and 0.6. Again we can see that the flow in the negative x direction is produced. The width between two fluid fluxes expands and the fluctuation is washed away before it forms a vortex because of the viscosity effect. The KH instability is not developed. In viscous fluid, a small size vortex compared with the Kolmogorov length scale is smeared by the viscosity and cannot exist. The fluctuation with the wavelength λ = 0.4 at τ 0 = 1 fm may be smaller than the Kolmogorov length scale. In the mid rapidity, |η| < 0.8, a fluctuation with longwave length disappears due to the Bjorken flow and a fluctuation with short wavelength is smeared by the viscosity. However, because at forward rapidity a fluctuation with long wavelength grows faster, there may be a chance that the KH instability occurs. Or if the longitudinal flow is smaller than Bjorken's flow, a fluctuation with long wavelength survives and can form a vortex. The existence of the KH instability depends on the viscosity and the flow distribution.

Summary
In this paper we have developed the new relativistic viscous hydrodynamics code. In the code, we employed the Milne coordinates which are suitable for the initial strong longitudinal expansion at high-energy heavy-ion collisions. After the brief explanation of the relativistic viscous hydrodynamic equations, we showed the numerical algorithm of the code which has the ideal part and the viscous part. For the ideal part we employed the Riemann solver with the two shock approximation which achieves stable calculation even with the small numerical viscosity [25] and for the viscous part we utilized the PES method [53]. Because we found that the order of accurate in space in the convection part of the viscous part is important, we applied the PPM instead of the MC limiter to the convection part. Next we examined the validity of our code using two test calculations; the viscous Bjorken flow for the one-dimensional test and the Israel-Stewart theory in the Gubser flow regime for the three-dimensional test. In both tests, our numerical calculations showed good agreement with analytical solutions. Besides, we pointed out that in the Gubser flow the shear tensors are sensitive to numerical scheme. Finally, we discussed the possible vortex formation through the KH instability in high-energy heavy-ion collisions. We focused on the mid rapidity and started the numerical calculations with the simple initial conditions inspired by the color-flux tube structure of hot spots in fluctuating initial conditions. In the case of the ideal fluid we found the vortex formation after τ ∼ 3 fm, however, we did not observe the vortex formation in the viscous fluid even with very small viscosity. To obtain a more conclusive result for the vortex formation in high-energy heavy-ion collisions, we need to use the more realistic initial conditions. For example, the existence of shear flow is found in the initial condition based on the Color Glass Condensate [65,66]. In addition, the effect of deviation from the Bjorken flow in a realistic initial condition is also important. Furthermore we shall apply our new code to analyses of experimental data at RHIC and the LHC; correlation between flow harmonics [16,17], event plane correlation [18,19], non-linearity of higher flow harmonics [20] and three particle correlation [21,22]. A comprehensive investigation of experimental data with the accurate numerical method of the relativistic viscous hydrodynamics gives us deep insight of QCD matter.
We give the explicit expressions for the interpolation procedure used in our relativistic viscous hydrodynamics code. We use the MC limiter for the second-order accurate in space and the PPM for the third-order accurate in space. We denote the center of the ith cell by x i and the boundary between the ith and the i + 1th cell by x i+1/2 . We assume that we have the average value a i of the quantity a(x) in the cell (x i−1/2 , x i+1/2 ) where a(x) stands for fluid variables and viscous tensors. In the interpolation procedure, we evaluate the values of a(x) at the right and left interfaces, a R,i = lim x→x i+1/2 a(x) and a L,i = lim x→x x−1/2 a(x) from the average value a i .

Appendix A.1: MC limiter
The second-order accuracy in space is achieved by the linear interpolation. In the second-order interpolation, we evaluate the interpolated values of a(x) at right and left interfaces, In the MC limiter [54], ∆ a i is given by We define space averages of an interpolation function, F i,R (σ i ) and F i,L (σ i ), where a I (x) is an interpolation function of a(x) and σ i = |u i |∆t/∆ x. Here we use the sound velocity (the fluid velocity) for u i in the conservation equation (the convection equation). We utilize F i,R (σ i ) and F i+1,L (σ i+1 ) for the initial condition of the Riemann problem at the cell interface x i+1/2 in the conservation equation. In the convection equation, F i,R (σ i ) or F i+1,L (σ i+1 ) corresponds to the numerical flux passing through the cell boundary x i+1/2 (Appendix B). In the linear interpolation, F i,R (σ i ) and F i,L (σ i ) are expressed by Appendix A.2: Piecewise Parabolic Method (PPM) [55][56][57] First, we calculate interpolated values of a(x) at cell interfaces using forth-order interpolation.
If the condition min(a i , a i+1 ) ≤ a i+1/2 ≤ max(a i , a i+1 ) is not satisfied, a i+1/2 is limited as follows: where C > 1 is a constant. We set C to C = 1.25 [57]. Then the interpolated values of a(x) at right and left interfaces are initiated as a L,i+1 = a R,i = a i+1/2 . We perform the flattening algorithm near strong shocks to prevent numerical oscillations, 14) The flattening parameter f i is fixed by where s j = +1 for p i+1 − p i−1 > 0 and s j = −1 for p i+1 − p i−1 < 0, (2) .
The constant w i is chosen by The parameters are set to ε = 1, w (1) = 0.52, and w (2) = 10 [56]. The flattening algorithm is applied for conservation equations. Furthermore, we modify the values of a i,R and a i,L to ensure the interpolated function remains monotonic. If (a i,R − a i )(a i − a i,L ) ≤ 0 or (a i−1 − a i )(a i − a i+1 ) ≤ 0, the ith cell contains a local extremum. The values of a i,R and a i,L are modified as follows:    where a n i is the value of a(t, x) at (t, x) = (t n , x i ), a n+1 i is the value of a at next time step t = t n+1 = t n + ∆t. The numerical flux a n+1/2 i+1/2 reads We evaluate the F i,R (σ i ) and In the case of multidimensional problems, we employ the Strang splitting method [52]. Using the operator L k i , which represents one-dimensional evolution in the i direction during the time k∆t, we express two-dimensional expansion in the (x, y) coordinates as a n+1 = L We consider two-dimensional convection equation, ∂ a(t, x, y) ∂t + u(x, y) ∂ a(t, x, y) ∂ x + v(x, y) ∂ a(t, x, y) ∂ y = 0. where a n i, j is the value of a(t, x, y) at t = t n , x = x i , y = y j , a n+1 i, j is the value of a(t, x, y) at next time step t = t n+1 = t n + ∆t, the second and third terms stand for the numerical flux passing through the cell boundary. The numerical flux is given by y (a n i, j − a n i, j−1 ) − min(v i, j , 0) ∆t 2∆ y (a n i, j+1 − a n i, j ) if u i, j ≥ 0, y (a n i+1, j − a n i+1, j−1 ) − min(v i+1, j , 0) ∆t 2∆ y (a n i+1, j+1 − a n i+1, j ) if u i, j < 0, (B.35) x (a n i, j − a n i−1, j ) − min(u i, j , 0) ∆t 2∆ x (a n i+1, j − a n i, j ) if v i, j ≥ 0, ∆ y a i, j+1 ∆ y − max(u i, j+1 , 0) ∆t 2∆ y (a n i, j+1 − a n i−1, j+1 ) − min(u i, j+1 , 0) ∆t 2∆ y (a n i+1, j+1 − a n i, j+1 ) if v i, j < 0. (B.36) Here we evaluate the variation of a(t, x, y) in the x (∆ x a i, j ) and y direction (∆ y a i, j ) using the MC limiter (A.2).