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

The European Physical Journal C, Jun 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-017-4944-0.pdf

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

Eur. Phys. J. C A new relativistic viscous hydrodynamics code and its application to the Kelvin-Helmholtz instability in high-energy heavy-ion collisions Kazuhisa Okamoto 2 Chiho Nonaka 0 1 2 0 Department of Physics, Duke University , Durham, NC 27708 , USA 1 Kobayashi-Maskawa Institute for the Origin of Particles and the Universe (KMI), Nagoya University , Nagoya 464-8602 , Japan 2 Department of Physics, Nagoya University , Nagoya 464-8602 , Japan 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. - 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 mechanism and the QGP bulk property is elucidated from comparison between hydrodynamic calculation and experimental data. a e-mail: b e-mail: 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–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 spacetime 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 heavy-ion collisions. Hydrodynamic instability and turbulent flow are discussed in Refs. [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. 2 Relativistic viscous hydrodynamic equations The relativistic hydrodynamics is based on the conservation equations, (1) T μν = euμuν − p μν , 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 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 with the metric tensor gμν . Here the uμ is determined uniquely. where τn, τπ , and τ are relaxation times, Inμ, Iπμν , and I represent second-order terms. nμNS, πNμSν , and NS are the 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 energymomentum tensor of the viscous fluid are decomposed as T μν = euμuν − ( p + 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–34]. The problem can be resolved by introducing the second-order 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–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–45]. In addition to the framework of the Israel–Stewart theory, other approaches such as the AdS/CFT correspondence [4,46–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 D Aμ1···μn ≡ uβ Aμβ1···μn , ; ∇α Aμ1···μn ≡ βα A;μβ1···μn , respectively. For example, in the second-order Israel–Stewart formalism the constitutive equations of the viscous tensors are given by μα νβ Dπ αβ = − τπ ( − NS) − I , 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 uμ = γ (1, vx , v y , vη), where vi (i = x , y, η) and γ = (1 − vx 2 − v y 2 − τ 2vη2)−1/2 are the three-velocity and the Lorentz factor, respectively. The conservation equations Eqs. (1) and (2) are explicitly written as 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 1 (∂τ + vi ∂i )nμ = − γ τn (nμ − nμNS) − Inμ − Jnμ − Knμ, (π μν −πNμνS)− Iπμν − Jπμν − Kπμν , NS) − I , where τn, τη and τ are the relaxation times, and the secondorder terms are defined by Jnτ = τ vηnη, Jnj = 0, Jnη = τ1 vηnτ + τ1 nη, J τ τ π = τ vηπ jη, π = 2τ vηπ τ η, J τ j j = x , y and λ = τ, x , y, η. Here, Jnτ and J μν are the π geometric source terms which come from the convective time derivative of nμ and π μν respectively. Knμ and Kπμν ensures μ the constraints nμuμ = 0, π μν uμ = 0 and πμ = 0. 3 Numerical algorithm In this section, we present our numerical algorithm for solving the relativistic viscous hydrodynamic equations in the Milne coordinates. 3.1 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: N μ = Niμd + Nvμis and T μν = Tiμdν + Tvμisν , where Niμd ≡ nuμ, Nvμis ≡ nμ, Tiμdν ≡ euμuν − p μν and Tvμisν ≡ π μν − μν . The subscripts “id” and “vis” mean the ideal part and the viscous part, respectively. The equations of the ideal part are expressed by where Siνd = −Tiτdτ /τ − τ Tiηdη, −Tiτdx /τ, −Tiτdy /τ, −3Tiτdη/τ . They are nothing but the usual ideal hydrodynamic equations in the Milne coordinates. On the other hand, the equations of the viscous part are given by i ∂τ (Niτd + Nvτis) + ∂i Nvis = − ∂τ (Tiτdν + Tvτisν ) + ∂i Tviiνs = Svνis, where Svνis = −Tvτisτ /τ − τ Tvηisη, −Tvτisx /τ, −Tvτisy /τ, −3Tvτisη/τ . 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, = 0, the relaxation equations, 1 = − γ τn 1 = − γ τη and the equations with source terms, ∂τ nμ = −Inμ − Jnμ − Knμ, ∂τ π μν = −Iπμν − Jπμν − Kπμν , = −I . 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 the evolutions of the viscous tensors, which means that a small relaxation time gives us a 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), nμ(τ ) = (n0μ − nμNS)exp − π μν (τ ) = (π0μν − πNμSν )exp − τγ−τπτ0 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 nx , n y , nη, π x x , π yy , π ηη, π x y , π 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. 3.2 Numerical implementation The decomposed hydrodynamic equations Eqs. (30)–(42) are solved by the following procedure. Here, we represent a conserved variable as U = U id + U vis, where U id ≡ (Niτd, Tiτdν ) and U vis ≡ (Nvτis, Tvτisν )(ν = τ, x , y, η). Fluid and dissipative variables are described by V id ≡ (n, p, vi ) and V vis ≡ (ni , π i j , )(i, j = x , y, η), respectively. 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 id(τ ) is evolved into U i∗d(τ + τ ), where the asterisk indicates a variable evolved only in the ideal part. 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 i∗d(τ + τ ) from U i∗d(τ + τ ) with the algorithm for recovery of the primitive variables V id from the conserved variables U id [25]. Second, we solve the constitutive equations of the viscous tensors Eqs. (34)–(42) to obtain V vis(τ + τ ). The convection equations Eqs. (34)–(36), the relaxation equations Eqs. (37)–(39) and Eqs. (40)–(42) are solved by the upwind scheme, the PES method and the predictor corrector method, respectively. The Navier–Stokes terms nμNS, πNμSν , NS and the second-order terms I , K in the right-hand sides of Eqs. (37)–(42) contain not only the spatial derivatives of fluid variables but also the time derivatives of them. The time derivatives in the right-hand sides of Eqs. (37)–(42) are obtained by ∂τ V id = (V i∗d(τ + τ ) − V id(τ ))/ τ. Here we keep the middle time-step value of the viscous tensor V vis(τ + τ /2) for the next step. Next, the conserved variables U i∗d(τ + τ ) and U vis(τ ) are evolved into U id(τ + τ ) and U vis(τ + τ ) by the viscous part of conservation equations Eqs. (32) and (33). Then we recover the fluid variables V id(τ + τ ) from conserved variables U vis(τ + τ ) [53]. We keep the middle time-step value V id(τ + /2). To achieve the second-order accurate in time, we repeat the above whole steps using the middle time-step values V id(τ + τ /2) and 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 the 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–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. 4 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 the 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. 4.1 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 = T0(τ0/τ )1/3, where τ0 and T0 are the initial proper time and temperature, respectively [51]. In the viscous fluid, the non-vanishing components of viscous tensor π μν are π x x , π yy , and π ηη. From the symmetries of the system, the relation 2π x x = 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 If η/s is constant, Eqs. (46) and (47) give the time evolution of the temperature, T = 1 − The numerical calculation is carried out on the space-grid size η = 0.1 with the time-step size τ = 0.1τ0 η. The initial temperature T0 and the proper time τ0 are set to T0 = 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, Fig. 1 The numerical and analytical results of the time evolution of the temperature in the Bjorken flow with and without shear viscosity. The viscosity to entropy density ratio is η/s = 0 and 0.2 Fig. 2 The numerical and analytical results of the time evolution of the Bulk pressure in the Bjorken flow. The bulk viscosity is ζ = 1000 MeV/fm2 ( − 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/fm2, τ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. 4.2 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, u⊥ 2q2τ x⊥ v⊥ = uτ = 1 + q2τ 2 + q2 x 2 ⊥ 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 τη = 5η/(T s). The numerical simulation starts at τ0 = 1 fm. The time-step size and the space-grid size in numerical simulation are set to τ = 0.1 x and ( x , y, η) = (0.05 fm, 0.05 fm, 0.1), respectively. Figure 3 shows the numerical results and the semi-analytic solutions of temperature and x component of fluid velocity as a function of x at τ = 1.2, 2 and 3 fm. The numerical results are consistent with the semi-analytic solutions. In our previous test calculation of the ideal Gubser flow the temperature and the fluid velocity follow the analytic solution until τ = 7 fm [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 π x x , π yy , π ηη and π x y at τ = 1.2, 2 and 3 fm are presented together with the semi-analytic solutions. Here the profile of π x y is shown along a line x = y, since the value of π x y vanishes on the x and y axes. The shear tensors π x x , π yy and π ηη in our numerical calculations show good agreement with the semi-analytic solutions. However, in π x y 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 Fig. 4 Comparison between the solutions for shear tensors π x x (top left), π yy (top right), τ 2π ηη (bottom left) and π xy (bottom right) from the Gubser flow and our numerical calculation as a function of x . The solid lines stand for the semi-analytic solutions and the pluses stand for numerical results 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 Appendices A and B. Figure 5 shows the semi-analytic solutions and numerical results of the shear tensors π x x , π yy , τ 2π ηη and π x y at τ = 3 fm. In π x x , π yy and τ 2π ηη, the results of all numerical schemes are reasonably consistent with the semi-analytic solutions. In addition, differences among them are small. However, in π x y 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. 5 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 the 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 heavyion 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’s scaling solution eB = e0(τ0/τ )4/3 and vηB = 0. Shifting the Bjorken scaling solution to ± z(= 0.3 fm) in the z direction, we obtain the energy density eU (eD) and the η component of velocity of the flow flux vUη (vDη) in x > 0 (x < 0), = e0 = e0 z coshη z) = τ 2 1 − τz sinhη , z coshη z) = − τ 2 1 + τz sinhη , where τ0 and e0 are the initial time and the energy density, respectively. Figure 6 shows the energy densities eU and eD, and the η component of velocity of the flow flux vUη and vDη in x > 0 and x < 0. The energy density and the η components of velocity of the flow flux are dependent on η in x > 0 and x < 0, because of the translational transformation of Bjorken’s scaling solution in the z direction. Importantly, one can see that the shear flow is created between the two initial flow fluxes. Furthermore, we put the fluctuation xb = 0.01sin(2π η/λ) with a wavelength λ along the boundary between the flow fluxes. Finally our initial energy density and flow velocity are written by Fig. 6 The energy density (top panel) and the rapidity component of velocity (bottom panel) of the translated Bjorken flow. The results at τ = τ0 = 1fm are shown. The red solid lines indicate the profile of Eqs. (54) and (56). The blue dotted lines indicate the profile of Eqs. (55) and (57) where the energy density and flow velocity around the boundary are connected from eU and vUη to eD 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 e0 = 741 GeV/fm3(T0 = 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 , ∂ ux The arrows stand for the velocity field in (τ vη, vx ). 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 = 3 p. ’ideal/xz0000.dat’ u 3:1:(-$8) -wy(fm-1) Fig. 7 Initial condition for the shear flow with the Bjorken expansion. The color profile show the distribution of vorticity −wy . The arrows indicate the three-fluid vector (τ vη, vx ) 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 Eqs. (58) and (59) and Fig. 6 suggest that eU is larger than eD and |vz | in x > 0 is larger than that in x > 0. The eU decreases more slowly than eD does due to the time dilation from larger |vz |. 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 the 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 Sect. 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 ’eta001/xz3000.dat’ u 3:1:(-$8) ’eta001/xz6000.dat’ u 3:1:(-$8) ’ideal/xz3000.dat’ u 3:1:(-$8) ’ideal/xz6000.dat’ u 3:1:(-$8) 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 long wavelength 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. 6 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 We define space averages of an interpolation function, Fi,R (σi ) and Fi,L (σi ), Fi,R = σi x xi+1/2−σi x where a I (x ) is an interpolation function of a(x ) and σi = |ui | t / x . Here we use the sound velocity (the fluid velocity) for ui in the conservation equation (the convection equation). We utilize Fi,R (σi ) and Fi+1,L (σi+1) for the initial condition of the Riemann problem at the cell interface xi+1/2 in the conservation equation. In the convection equation, Fi,R (σi ) or Fi+1,L (σi+1) corresponds to the numerical flux passing through the cell boundary xi+1/2 (Appendix B). In the linear interpolation, Fi,R (σi ) and Fi,L (σi ) are expressed by Appendix A: Interpolation procedures Appendix A.2: Piecewise parabolic method (PPM) [55–57] 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. Acknowledgements We would like to thank Berndt Mueller and Rainer Fries for valuable discussions. C.N. is grateful for the hospitality of the members of Cyclotron Institute and Department of Physics and Astronomy in Texas A&M. The work of C.N. is supported by the JSPS Grant-in-Aid for Scientific Research (S) No. 26220707 and US Department of Energy grant DE-FG02-05ER41367. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP3. 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 i th cell by xi and the boundary between the i th and the i + 1th cell by xi+1/2. We assume that we have the average value ai of the quantity a(x ) in the cell (xi−1/2, xi+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, aR,i = limx→xi+1/2 a(x ) and aL,i = limx→xx−1/2 a(x ) from the average value ai . 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, aR,i = ai + ai /2, aL,i = ai − In the MC limiter [54], ai is given by = 0 otherwise. × sign(ai+1 − ai−1) if (ai+1 − ai )(ai − ai−1) > 0, (A.2) First, we calculate the interpolated values of a(x ) at cell interfaces using fourth-order interpolation: 7 1 ai+1/2 = 12 (ai + ai+1) − 12 (ai−1 + ai+2). If the condition min(ai , ai+1) ≤ ai+1/2 ≤ max(ai , ai+1) is not satisfied, ai+1/2 is limited as follows: (D2a)i+1/2 = (D2a)i+1/2,L = (D2a)i+1/2,R = If the signs of (D2a)i+1/2, (D2a)i+1/2,R and (D2a)i+1/2,L are all the same, |(D2a)i+1/2| sign((D2a)i+1/2), otherwise, (D2a)i+1/2,lim = 0. Then the modified values of ai+1/2 read 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 aL,i+1 = aR,i = ai+1/2. We perform the flattening algorithm near strong shocks to prevent numerical oscillations, aR,i → ai fi + aR,i (1 − fi ), aL,i → ai fi + aL,i (1 − fi ). The flattening parameter fi is fixed by fi = max( f˜i , f˜i+si ), where s j = +1 for pi+1 − pi−1 > 0 and s j = −1 for pi+1 − pi−1 < 0, f˜i = min 1, wi max 0, The constant wi is chosen by pi+2 − pi−2 pi+1 − pi−1 − w(1) > , vi−1 > vi+1, 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 ai,R and ai,L to ensure the interpolated function remains monotonic. If (ai,R − ai )(ai − ai,L ) ≤ 0 or (ai−1 − ai )(ai − ai+1) ≤ 0, the i th cell contains a local extremum. The values of ai,R and ai,L are modified as follows: (D2a)i = − 2a6,i , (D2a)i,C = (D2a)i,L = (D2a)i,R = where a6,i = 6ai − 3(ai,L + ai,R ). If (D2a)i and (D2a)i,{L,C,R} have the same sign, |(D2a)i |)sign((D2a)i ), otherwise, (D2a)i,lim = 0. Then we obtain (D2a)i,lim , ai,R → ai + (ai,R − ai ) (D2a)i (D2a)i,lim . ai,L → ai + (ai,L − ai ) (D2a)i If (D2a)i = 0, we set the second term of Eqs. (A.22) and (A.23) to be zero. In the last limiter, the values of ai,R and ai,L are modified as The space averages of a parabolic interpolant are written ai,R − ai,L − ai,R − ai,L + Again, Fi,R (σi ) and Fi+1,L (σi+1) are used for the initial condition of the Riemann problem at the cell interface xi+1/2 in the conservation equation. In the convection equation, Fi,R (σi ) or Fi+1,L (σi+1) corresponds to the numerical flux passing through the cell boundary xi+1/2 (Appendix B). Appendix B: Numerical schemes for convection equations Appendix B.1: High-resolution upwind method We consider the one-dimensional convection equation, = 0. In the high-resolution upwind method, we obtain the solution of the convection equation Eq. (B.28), ain+1 = ain − where ain is the value of a(t, x ) at (t, x ) = (t n, xi ), ain+1 is the value of a at next time step t = t n+1 = t n + t . The numerical flux ain++11//22 reads We evaluate the Fi,R (σi ) and Fi,L (σi ), using the MC limiter (Eqs. (A.5) and (A.6)) or the PPM (Eqs. (A.26) and (A.27)). In the case of multidimensional problems, we employ the Strang splitting method [52]. Using the operator Lik , which represents one-dimensional evolution in the i direction during the time k t , we express the two-dimensional expansion in the (x , y) coordinates as an+1 = L 1x/2 L1y L 1x/2an. Similarly the three-dimensional expansion in (x , y, z) coordinates is written by × L 1z/3 L 1x/6 L1y/3 L 1z/6 L 1x/6an. Appendix B.2: Corner transport upwind (CTU) scheme [63] We consider two-dimensional convection equation, ∂a(t, x , y) ∂t ∂a(t, x , y) ∂ x ∂a(t, x , y) ∂ y In the CTU, the solution of the convection equation Eq. (B.33) reads − max(ui, j , 0) (ain++11//22, j − ain−+11//22, j ) (ain, +j+1/12/2 − ain, +j−1/12/2), ain, +j1 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 an+1/2 n i+1/2, j = ai, j + − max(vi, j , 0) − min(vi, j , 0) − max(vi+1, j , 0) − min(vi+1, j , 0) (ain, j+1 − ain, j ) if ui, j ≥ 0, = 0. − max(ui, j+1, 0) − min(ui, j+1, 0) − min(ui, j , 0) (ain+1, j − ain, j ) if vi, j ≥ 0, 2 y (ain+1, j+1 − ain, j+1) Here we evaluate the variation of a(t, x , y) in the x ( x ai, j ) and y direction ( y ai, j ) using the MC limiter (A.2). 1. T.D. Lee , Quark gluon plasma, new discoveries at RHIC: case for the strongly interacting quark-gluon plasma . Nucl. Phys. A 750 , 1 - 172 ( 2005 ) 2. J.M. Stewart , Proc. R. Soc . Lond . A 357 , 59 ( 1977 ). doi:10.1016/ 0003 -49 1 3. W. Israel , J. M. Stewart , Ann. Phys. (N.Y.) 118 , 341 ( 1979 ). doi:10. 1016/ 0003 -49 1 4. R. Baier , P. Romatschke , D.T. Son , A.O. Starinets , M.A. Stephanov , JHEP 04 , 100 ( 2008 ). doi:10.1088/ 1126 - 6708 /2008/ 04/100. arXiv: 0712 .2451 [hep-th] 5. H. Song , U.W. Heinz , Phys. Rev. C 77 , 064901 ( 2008 ). doi:10. 1103/PhysRevC.77.064901. arXiv: 0712 .3715 [nucl-th] 6. K. Dusling , D. Teaney , Phys. Rev. C 77 , 034905 ( 2008 ). doi:10. 1103/PhysRevC.77.034905. arXiv: 0710 .5932 [nucl-th] 7. M. Luzum , P. Romatschke , Phys. Rev. C 78 , 034915 ( 2008 ). doi:10. 1103/PhysRevC.78.034915. arXiv: 0804 .4015 [nucl-th] 8. G.S. Denicol , T. Kodama , T. Koide , P. Mota , Phys. Rev. C 80 , 064901 ( 2009 ). doi:10.1103/PhysRevC.80.064901. arXiv: 0903 .3595 [hep-ph] 9. B. Schenke , S. Jeon , C. Gale , Phys. Rev. Lett . 106 , 042301 ( 2011 ). doi:10.1103/PhysRevLett.106.042301. arXiv: 1009 .3244 [hep-ph] 10. P. Bozek , Phys. Rev. C 85 , 034901 ( 2012 ). doi:10.1103/PhysRevC. 85.034901. arXiv: 1110 .6742 [nucl-th] 11. V. Roy , A.K. Chaudhuri , Phys. Rev. C 85 , 024909 ( 2012 ). doi:10. 1103/PhysRevC.85.024909. arXiv: 1109 .1630 [nucl-th] 12. I. Karpenko , P. Huovinen , M. Bleicher , Comput. Phys. Commun . 185 , 3016 ( 2014 ). doi:10.1016/j.cpc. 2014 .07.010. arXiv: 1312 .4160 [nucl-th] 13. J. Noronha-Hostler , J. Noronha , F. Grassi , Phys. Rev. C 90 , 034907 ( 2014 ). doi:10.1103/PhysRevC.90.034907. arXiv: 1406 .3333 [nucl-th] 14. L.G. Pang , Y. Hatta , X.N. Wang , B.W. Xiao , Phys. Rev. D 91 , 074027 ( 2015 ). doi:10.1103/PhysRevD.91.074027. arXiv: 1411 .7767 [hep-ph] 15. C. Gale , S. Jeon , B. Schenke , P. Tribedy , R. Venugopalan , Phys. Rev. Lett . 110 , 012302 ( 2013 ). doi:10.1103/PhysRevLett.110. 012302. arXiv: 1209 .6330 [nucl-th] 16. G. Aad et al., ATLAS Collaboration, Phys. Rev. C 92 , 034903 ( 2015 ). doi:10.1103/PhysRevC.92.034903. arXiv: 1504 .01289 [hep-ex] 17. J. Adam et al., ALICE Collaboration, Phys. Rev. Lett 117 , 182301 ( 2016 ). doi:10.1103/PhysRevLett.117.182301. arXiv: 1604 .07663 [nucl-ex] 18. G. Aad et al., ATLAS Collaboration, Phys. Rev. C 90 , 024905 ( 2014 ). doi:10.1103/PhysRevC.90.024905. arXiv: 1403 .0489 [hepex] 19. V. Khachatryan et al., CMS Collaboration, Phys. Rev. C 92 , 034911 ( 2015 ). doi:10.1103/PhysRevC.92.034911. arXiv: 1503 .01692 [nucl-ex] 20. L. Yan , J.-Y. Ollitrault , Phys. Lett . B 744 , 744 , 82 ( 2015 ). doi:10. 1016/j.physletb. 2015 .03.040. arXiv: 1502 .02502 [nucl-th] 21. L. Adamczyk et al. , STAR Collaboration , arXiv:1701.06496 [nuclex] 22. L. Adamczyk et al. , STAR Collaboration , arXiv:1701.06497 [nuclex] 23. A. Mignone , T. Plewa , G. Bodo , ApJS 160 , 199 ( 2005 ). doi:10. 1086/430905. arXiv:astro-ph/ 0505200 24. K, Okamoto , Y. Akamatsu , C. Nonaka , Eur. Phys. J. C (2016) 76 : 579. doi:10.1140/epjc/s10052- 016 - 4433 -x. arXiv: 1607 .03630 [nucl-th] 25. Y. Akamatsu , S. Inutsuka , C. Nonaka , M. Takamoto , J. Comput . Phys. 256 , 34 ( 2014 ). doi:10.1016/j.jcp. 2013 .08.047. arXiv: 1302 .1665 [nucl-th] 26. P. Romatschke , Theor. Phys. Suppl . 174 , 137 ( 2008 ). doi:10.1143/ PTPS.174.137. arXiv: 0710 .0016 [nucl-th] 27. S. Floerchinger , U.A. Wiedemann , J. High Energy Phys . 11 , 100 ( 2011 ). doi:10.1007/ JHEP11(2011)100 . arXiv: 1108 .5535 [nuclth] 28. L. P. Csernai , D. D. Strottman , Cs. Anderlik, Phys. Rev. C 85 , 054901 ( 2012 ). doi:10.1103/PhysRevC.85.054901. arXiv: 1112 .4287 [nucl-th] 29. C. Eckart , Phys. Rev . 58 , 919 ( 1940 ). doi:10.1103/PhysRev.58.919 30. L.D. Landau , E.M. Lifshitz , Fluid Mechanics , 2nd edn, Course of Theoretical Physics , vol. 6 ( Butterworth-Heinemann , Oxford, 1987 ) 31. A. Andronic , P. Braun-Munzinger , K. Redlich , J. Stachel , Nucl. Phys. A 904 , 535c ( 2013 ). doi:10.1016/j.nuclphysa. 2013 .02.070. arXiv: 1210 .7724 [nucl-th] 32. W.A. Hiscock , L. Lindblom, Ann. Phys. (NY) 151 , 466 ( 1983 ). doi:10.1016/ 0003 -4 9 33. W.A. Hiscock , L. Lindblom , Phys. Rev. D 31 , 725 ( 1985 ). doi:10. 1103/PhysRevD.31.725 34. S. Pu , T. Koide , D.H. Rischke , Phys. Rev. D 81 , 114039 ( 2010 ). doi:10.1103/PhysRevD.81.114039. arXiv: 0907 .3906 [hep-ph] 35. P. Huovinen , D. Molnar , Phys. Rev. C 79 , 014906 ( 2009 ). doi:10. 1103/PhysRevC.79.014906. arXiv: 0808 .0953 [nucl-th] 36. I. Bouras , E. Molnar , H. Niemi , Z. Xu , A. El , O. Fochler , C. Greiner , D.H. Rischke , Phys. Rev. C 82 , 024910 ( 2010 ). doi:10. 1103/PhysRevC.82.024910 37. D. Molnar , P. Huovinen , Nucl. Phys. A 830 , 475c ( 2009 ). doi:10. 1016/j.nuclphysa. 2009 .10.104. arXiv: 0907 .5014 [nucl-th] 38. M. Takamoto , S. Inutsuka , Physica (Amsterdam) 389, 4580 ( 2010 ). doi:10.1016/j.jcp. 2011 .05.030. arXiv: 1106 .1732 [astro-ph] 39. A. El , Z. Xu , C. Greiner , Phys. Rev. C 81 , 041901 ( 2010 ). doi:10. 1103/PhysRevC.81.041901. arXiv: 0907 .4500 [hep-ph] 40. A. Muronga , Phys. Rev. C 76 , 014910 ( 2007 ). doi:10.1103/ PhysRevC.76.014910. arXiv:nucl-th/0611091 41. B. Betz , D. Henkel , D.H. Rischke , J. Phys . G 36 , 064029 ( 2009 ). doi:10.1051/epjconf/20111307005. arXiv: 1012 .5772 [nucl-th] 42. G.S. Denicol , T. Koide , D.H. Rischke , Phys. Rev. Lett . 105 , 162501 ( 2010 ). doi:10.1103/PhysRevLett.105.162501. arXiv: 1004 .5013 [nucl-th] 43. E. Calzetta , J. Peralta-Ramos , Phys. Rev. D 82 , 106003 ( 2010 ). doi:10.1103/PhysRevD.82.106003. arXiv: 1009 .2400 [hep-ph] 44. G.S. Denicol , H. Niemi , E. Molnar , D.H. Rischke , Phys. Rev. D 85 , 114047 ( 2012 ). doi:10.1103/PhysRevD.85.114047. arXiv: 1202 .4551 [nucl-th] 45. A. Jaiswal , Phys. Rev. C 87 , 051901 ( 2013 ). doi:10.1103/ PhysRevC.87.051901. arXiv: 1302 .6311 [nucl-th] 46. S. Bhattacharyya , V.E. Hubeny , S. Minwalla , M. Rangamani , J. High Energy Phys . 02 , 045 ( 2008 ). doi:10.1088/ 1126 - 6708 /2008/ 02/045. arXiv: 0712 .2456 [hep-th] 47. M. Natsuume , T. Okamura , Phys. Rev. D 77 , 066014 ( 2008 ). doi:10. 1103/PhysRevD.77.066014. arXiv: 0712 .2916 [hep-th] 48. P. Romatschke , Class. Quantum Gravity 27 , 025006 ( 2010 ). doi:10. 1088/ 0264 - 9381 /27/2/025006. arXiv: 0906 .4787 [hep-th] 49. K. Tsumura , T. Kunihiro , Eur. Phys. J. A 48 , 162 ( 2012 ). doi:10. 1140/epja/i2012- 12162 -x. arXiv: 1206 . 1929 [nucl-th] 50. K. Tsumura , Y. Kikuchi , T. Kunihiro , Phys. Rev. D 92 , 085048 ( 2015 ). doi:10.1103/PhysRevD.92.085048. arXiv: 1506 .00846 [hep-ph] 51. J.D. Bjorken, Phys. Rev. D 27 , 140 ( 1983 ). doi:10.1103/PhysRevD. 27.140 52. G. Strang , SIAM J. Numer . Anal. 5 , 506 ( 1968 ). doi:10.1137/ 0705041 53. M. Takamoto , S. Inutsuka , J. Comput . Phys. 11 , 38 ( 2011 ). doi:10. 1016/j.jcp. 2011 .05.030. arXiv: 1106 .1732 [astro-ph] 54. B. Van Leer , J. Comput. Phys . 32 , 101 ( 1979 ). doi:10.1016/ 0021 -999 1 55. P. Colella , P.R. Woodward , J. Comput . Phys. 54 , 174 ( 1984 ). doi:10. 1016/ 0021 -9991( 8 56. J.M. Marti , E. Müller , J. Comput . Phys. 123 , 1 ( 1996 ). doi:10.1006/ jcph.1996.0001 57. P. Colella , M.D. Sekora , J. Comput . Phys. 227 , 7069 ( 2008 ). doi:10. 1016/j.jcp. 2008 .03.034 58. H. Marrochio , J. Noronha , G.S. Denicol , M. Luzum , S. Jeon , C. Gale , Phys. Rev. C 91 , 014903 ( 2015 ). doi:10.1103/PhysRevC.91. 014903. arXiv: 1307 .6130 [nucl-th] 59. S.S. Gubser , Phys. Rev. D 82 , 085027 ( 2010 ). doi:10.1103/ PhysRevD.82.085027. arXiv: 1006 .0006 [hep-th] 60. S.S. Gubser , A. Yarom , Nucl. Phys . B 846 , 469 ( 2011 ). doi:10. 1016/j.nuclphysb. 2011 .01.012. arXiv: 1012 .1314 [hep-th] 61. C. Shen , Z. Qiu , H. Song , J. Bernhard , S. Bass , U. Heinz , Comput. Phys. Commun . 199 , 61 ( 2016 ). doi:10.1016/j.cpc. 2015 .08.039. arXiv: 1409 .8164 [nucl-th] 62. F. Becattini , G. Inghirami , V. Rolando , A. Beraudo , L. Del Zanna , A. De Pace , M. Nardi , G. Pagliara , V. Chandra , Eur. Phys. J. C 75 , 406 ( 2015 ). doi:10.1140/epjc/s10052- 015 - 3624 -1. arXiv: 1501 .04468 [nucl-th] 63. P. Collela , J. Comput . Phys. 87 , 171 ( 1990 ). doi:10.1016/ 0021 -9991(90) 90233 - Q 64. P.G. Drazin , W.H. Reid , Hydrodynamic Stability (Cambridge University Press, London, 1981 ) 65. G. Chen , R.J. Fries , J.I. Kapusta , Y. Li , Phys. Rev. C 92 , 064912 ( 2015 ). doi:10.1103/PhysRevC.92.064912. arXiv: 1507 .03524 [nucl-th] 66. R. J. Fries , G. Chen , S. Somanathan , private communication


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-017-4944-0.pdf

Kazuhisa Okamoto, Chiho Nonaka. A new relativistic viscous hydrodynamics code and its application to the Kelvin–Helmholtz instability in high-energy heavy-ion collisions, The European Physical Journal C, 2017, DOI: 10.1140/epjc/s10052-017-4944-0