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

TL;DR

This paper introduces a new relativistic viscous hydrodynamics code optimized for Milne coordinates, combining advanced numerical methods, and applies it to study the Kelvin-Helmholtz instability in high-energy heavy-ion collisions.

Contribution

The paper presents a novel relativistic viscous hydrodynamics code with a unique splitting method and validation, applied to Kelvin-Helmholtz instability analysis in heavy-ion collisions.

Findings

Validated the code against analytical solutions and known theories.

Demonstrated the code's capability to simulate Kelvin-Helmholtz instability.

Provided insights into instability development in high-energy collisions.

Abstract

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.

Equations

$$N_{;\mu }^{\mu }$$

$$T_{;\mu }^{\mu \nu }$$

$$N^{\mu }$$

$$T^{\mu \nu }$$

$${\Delta }^{\mu \nu }=g^{\mu \nu }-u^{\mu }{u}^{\nu },$$

$$N^{\mu }$$

$$T^{\mu \nu }$$

$$D{A}^{{\mu }_1\cdots {\mu }_n}\equiv u^{\beta }{A}_{;\beta }^{{\mu }_1\cdots {\mu }_n},$$

$${\nabla }_{\alpha }{A}^{{\mu }_1\cdots {\mu }_n}\equiv {\Delta }_{\alpha }^{\beta }{A}_{;\beta }^{{\mu }_1\cdots {\mu }_n},$$

$${\Delta }_{\alpha }^{\mu }D{n}^{\alpha }$$

$${\Delta }_{\alpha }^{\mu }{\Delta }_{\beta }^{\nu }D{\pi }^{\alpha \beta }$$

$$D\Pi$$

$$n_{\mathrm{NS}}^{\mu }$$

$${\pi }_{\mathrm{NS}}^{\mu \nu }$$

$${\Pi }_{NS}$$

$${\partial }_{\tau }{N}^{\tau }+{\partial }_i{N}^i=-\frac{1}{\tau }{N}^{\tau },$$

$${\partial }_{\tau }{T}^{\tau \nu }+{\partial }_i{T}^{i\nu }=S^{\nu },$$

$$S^{\nu }=\left(-\frac{1}{\tau }{T}^{\tau \tau }-\tau T^{\eta \eta },-\frac{1}{\tau }{T}^{\tau x},-\frac{1}{\tau }{T}^{\tau y},-\frac{3}{\tau }{T}^{\tau \eta }\right).$$

$$({\partial }_{\tau }+v^i{\partial }_i)n^{\mu }$$

$$({\partial }_{\tau }+v^i{\partial }_i){\pi }^{\mu \nu }$$

$$({\partial }_{\tau }+v^i{\partial }_i)\Pi$$

$$J_n^{\tau }=\tau v^{\eta }{n}^{\eta },J_n^j=0,$$

$$J_n^{\eta }=\frac{1}{\tau }{v}^{\eta }{n}^{\tau }+\frac{1}{\tau }{n}^{\eta },$$

$$J_{\pi }^{\tau \tau }=2\tau v^{\eta }{\pi }^{\tau \eta },J_{\pi }^{\tau j}=\tau v^{\eta }{\pi }^{j\eta },$$

$$J_{\pi }^{xx}=J_{\pi }^{yy}=J_{\pi }^{xy}=0,$$

$$J_{\pi }^{i\eta }=\frac{1}{\tau }{\pi }^{i\eta }+\frac{1}{\tau }{v}^{\eta }{\pi }^{i\tau },$$

$$J_{\pi }^{\eta \eta }=\frac{2}{\tau }{v}^{\eta }{\pi }^{\tau \eta }+\frac{2}{\tau }{\pi }^{\eta \eta },$$

$$K_n^{\mu }=n^{\lambda }{v}^{\mu }D{u}_{\lambda },$$

$$K_{\pi }^{\mu \nu }=({\pi }^{\lambda \mu }{v}^{\nu }+{\pi }^{\lambda \nu }{v}^{\mu })D{u}_{\lambda },$$

$${\partial }_{\tau }{N}_{\mathrm{id}}^{\tau }+{\partial }_i{N}_{\mathrm{id}}^i=-\frac{1}{\tau }{N}_{\mathrm{id}}^{\tau },$$

$${\partial }_{\tau }{T}_{\mathrm{id}}^{\tau \nu }+{\partial }_i{T}_{\mathrm{id}}^{i\nu }=S_{\mathrm{id}}^{\nu },$$

$${\partial }_{\tau }(N_{\mathrm{id}}^{\tau }+N_{\mathrm{vis}}^{\tau })+{\partial }_i{N}_{\mathrm{vis}}^i=-\frac{1}{\tau }{N}_{\mathrm{vis}}^{\tau },$$

$${\partial }_{\tau }(T_{\mathrm{id}}^{\tau \nu }+T_{\mathrm{vis}}^{\tau \nu })+{\partial }_i{T}_{\mathrm{vis}}^{i\nu }=S_{\mathrm{vis}}^{\nu },$$

$$({\partial }_{\tau }+v^i{\partial }_i)n^{\mu }$$

$$({\partial }_{\tau }+v^i{\partial }_i){\pi }^{\mu \nu }$$

$$({\partial }_{\tau }+v^i{\partial }_i)\Pi$$

$${\partial }_{\tau }{n}^{\mu }$$

$${\partial }_{\tau }{\pi }^{\mu \nu }$$

$${\partial }_{\tau }\Pi$$

$${\partial }_{\tau }{n}^{\mu }$$

$${\partial }_{\tau }{\pi }^{\mu \nu }$$

$${\partial }_{\tau }\Pi$$

$$n^{\mu }(\tau )$$

Full text

††journal: Eur. Phys. J. C

1 Introduction

Since the success of production of the strongly interacting quark-gluon plasma (QGP) at Relativistic Heavy Ion Collider (RHIC) QGP-RHIC , 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 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 Stewart1977 ; Israel1979 and conformal hydrodynamics Baier2007 are often used. Solving them numerically, study of experimental data of high-energy heavy-ion collisions is performed Song2007 ; Dusling2007 ; Luzum2008 ; Denicol2009 ; Schenke2011 ; Bozek2011 ; Roy2012 ; Karpenko2014 ; Noronha2014 ; Pang2014 .

Now the relativistic viscous hydrodynamic model can explain not only the elliptic flow but also higher harmonics Gale2013 . 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 ATLAS2015 ; ALICE2016 , event plane correlation ATLAS2014 ; CMS2015 , non-linearity of higher flow harmonics Yan2015 and three particle correlation STAR2017a ; STAR2017b . At the same time, we can investigate the QGP property further using information of (3+1)-dimensional space-time evolutions CMS2015 ; STAR2017a ; STAR2017b . 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 Mignone2005 is employed Okamoto2016 . It is stable even with small numerical viscosity Akamatsu2014 . 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 Ref. Romatschke2008 ; Floerchinger2011 and the possibility of KH instability is argued in Ref. Csernai2012 . 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,

[TABLE]

where $N^{\mu}$ is the net charge current and $T^{\mu\nu}$ is the energy-momentum tensor. In the case of ideal fluid, the net charge current and energy-momentum tensor are given by

[TABLE]

where $n$ is the net charge density, $e$ is the energy density, $p$ is the pressure and $u^{\mu}$ is the fluid four-velocity which satisfies the normalization $u^{\mu}u_{\mu}=1$. $\Delta^{\mu\nu}$ is the orthogonal projection tensor to $u^{\mu}$, which is defined by

[TABLE]

with the metric tensor $g^{\mu\nu}$. Here the $u^{\mu}$ is determined uniquely.

On the other hand, in dissipative flow, there are several possible choices to determine $u^{\mu}$. For example, one can assign the $u^{\mu}$ as net charge flow (Eckart frame Eckart1940 ) or as energy flow (Landau frame Landau1987 ). The decomposition of $N^{\mu}$ and $T^{\mu\nu}$ in viscous fluid depends on the choice of $u^{\mu}$. 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 Andronic2013 .

In the Landau frame, the net charge current and the energy-momentum tensor of the viscous fluid are decomposed as

[TABLE]

where $n^{\mu}$ is the charge diffusion current, $\Pi$ is the bulk pressure, and $\pi^{\mu\nu}$ is the shear tensor Landau1987 . The relativistic extension of Navier-Stokes theory in non-relativistic fluid usually has a problem of acausality and instability Hiscock1983 ; Hiscock1985 ; Pu2010 . 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 Stewart1977 ; Israel1979 . However, the original Israel-Stewart theory does not reproduce the results of the kinetic equation quantitatively Huovinen2009 ; Bouras2010 ; Molnar2009 ; Takamoto2010 ; Xu2010 . The construction of second-order relativistic viscous hydrodynamic equations is still under investigation. The extension of the Israel-Stewart theory is also proposed Muronga2007 ; Betz2009 ; Denicol2010 ; Calzetta2010 ; Denicol2012 ; Jaiswal2013 . In addition to the framework of the Israel-Stewart theory, other approaches such as the AdS/CFT correspondence Baier2007 ; Natsuume2008 ; Bhatta2008 ; Romatschke2010 and renormalization group method are applied to the construction of causal relativistic hydrodynamics Tsumura2012 ; Tsumura2015 .

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 $\nabla^{\mu}$, which are defined by

[TABLE]

respectively. For example, in the second-order Israel-Stewart formalism the constitutive equations of the viscous tensors are given by

[TABLE]

where $\tau_{n},\tau_{\pi}$, and $\tau_{\Pi}$ are relaxation times, $I^{\mu}_{n},I^{\mu\nu}_{\pi},$ and $I_{\Pi}$ represent second-order terms. $n^{\mu}_{\rm NS},\pi^{\mu\nu}_{\rm NS}$, and $\Pi_{\rm NS}$ are the Navier-Stokes value of viscous tensors written as

[TABLE]

where $T$ is the temperature, $\mu$ is the chemical potential, $\theta\equiv u^{\mu}_{;\mu}$ is the expansion scalar, $\sigma$ is the charge conductivity, $\eta$ is the shear viscosity, and $\Pi$ is the bulk viscosity.

We construct a relativistic viscous hydrodynamics code in the Milne coordinates $(\tau,x,y,\eta)$ which is optimized for description of the strong longitudinal expansion Bjorken1983 at RHIC and the LHC. In the Milne coordinates, the metric tensor is given by $g^{\mu\nu}={\rm diag}(1,-1,-1,-1/\tau^{2})$ and the fluid four-velocity has the form $u^{\mu}=\gamma(1,v^{x},v^{y},v^{\eta})$, where $v^{i}(i=x,y,\eta)$ and $\gamma=(1-{v^{x}}^{2}-{v^{y}}^{2}-\tau^{2}{v^{\eta}}^{2})^{-1/2}$ are the three-velocity and the Lorentz factor, respectively. The conservation equations Eqs. (1) and (2) are explicitly written as

[TABLE]

where $i=x,y,\eta$ and the right-hand sides of them represent geometric source terms. $S^{\nu}$ is given by

[TABLE]

The constitutive equations Eqs. (10), (11), and (12) in Milne coordinates read

[TABLE]

where $\tau_{n}$, $\tau_{\eta}$ and $\tau_{\Pi}$ are the relaxation times, and the second-order terms are defined by

[TABLE]

$j=x,y$ and $\lambda=\tau,x,y,\eta$. Here, $J^{\tau}_{n}$ and $J^{\mu\nu}_{\pi}$ are the geometric source terms which come from the convective time derivative of $n^{\mu}$ and $\pi^{\mu\nu}$ respectively. $K^{\mu}_{n}$ and $K^{\mu\nu}_{\pi}$ ensures the constraints $n^{\mu}u_{\mu}=0,\pi^{\mu\nu}u_{\mu}=0$ and $\pi^{\mu}_{\mu}=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 Strang1968 . Specifically, the net charge current and the energy-momentum tensor are divided as follows: $N^{\mu}=N^{\mu}_{\rm id}+N^{\mu}_{\rm vis}$ and $T^{\mu\nu}=T^{\mu\nu}_{\rm id}+T^{\mu\nu}_{\rm vis}$, where $N^{\mu}_{\rm id}\equiv nu^{\mu},N^{\mu}_{\rm vis}\equiv n^{\mu}$, $T^{\mu\nu}_{\rm id}\equiv eu^{\mu}u^{\nu}-p\Delta^{\mu\nu}$ and $T^{\mu\nu}_{\rm vis}\equiv\pi^{\mu\nu}-\Pi\Delta^{\mu\nu}$. The subscripts “id” and “vis” mean the ideal part and the viscous part, respectively. The equations of the ideal part are expressed by

[TABLE]

where $S^{\nu}_{\rm id}=\left(-T^{\tau\tau}_{\rm id}/\tau-\tau T^{\eta\eta}_{\rm id},\>-T^{\tau x}_{\rm id}/\tau,\>-T^{\tau y}_{\rm id}/\tau,\>-3T^{\tau\eta}_{\rm id}/\tau\right)$. 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

[TABLE]

where $S^{\nu}_{\rm vis}=\left(-T^{\tau\tau}_{\rm vis}/\tau-\tau T^{\eta\eta}_{\rm vis},\>-T^{\tau x}_{\rm vis}/\tau,\>-T^{\tau y}_{\rm vis}/\tau,\>-3T^{\tau\eta}_{\rm vis}/\tau\right)$. 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,

[TABLE]

the relaxation equations,

[TABLE]

and the equations with source terms,

[TABLE]

In numerical simulation of relativistic hydrodynamic equation, a time-step size $\Delta\tau$ is usually determined by the Courant-Friedrichs-Lewy (CFL) condition. However, in the relativistic dissipative hydrodynamics, one needs to determine the value of $\Delta\tau$ carefully. The relaxation times $\tau_{n},\tau_{\eta}$, and $\tau_{\Pi}$ 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 $\Delta\tau$ than the CFL condition does. If the relaxation times $\tau_{n},\tau_{\eta}$, and $\tau_{\Pi}$ are much shorter than the fluid timescale $\tau_{\rm fluid}$, the time-step size $\Delta\tau$ 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 Inutsuka2011 , instead of using a simple explicit scheme. In the PES method, formal solutions of Eqs. (37)-(39),

[TABLE]

can be used. On the other hands, if the relaxation times are larger than $\Delta\tau$ determined by the CFL condition, the PES method is not applied Karpenko2014 .

In our algorithm, we solve the time evolution of $n^{x},n^{y},n^{\eta},\pi^{xx},\pi^{yy},\pi^{\eta\eta},\pi^{xy},\pi^{y\eta},\pi^{\eta x}$ and $\Pi$ directly. Other components of viscous tensors $n^{\tau},\pi^{\tau\tau},\pi^{\tau x},\pi^{\tau y}$ and $\pi^{\tau\eta}$ are derived from the orthogonality conditions $n^{\mu}u_{\mu}=0$ and $\pi^{\mu\nu}u_{\nu}=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 $\bm{U}=\bm{U}_{\rm id}+\bm{U}_{\rm vis}$, where $\bm{U}_{\rm id}\equiv(N^{\tau}_{\rm id},T^{\tau\nu}_{\rm id})$ and $\bm{U}_{\rm vis}\equiv(N^{\tau}_{\rm vis},T^{\tau\nu}_{\rm vis})(\nu=\tau,x,y,\eta)$. Fluid and dissipative variables are described by $\bm{V}_{\rm id}\equiv(n,p,v^{i})$ and $\bm{V}_{\rm vis}\equiv(n^{i},\pi^{ij},\Pi)(i,j=x,y,\eta)$, respectively.

First, we solve the ideal part of the conservation equations Eqs. (30) and (31) using the Riemann solver Okamoto2016 . In this step, the conserved variable $\bm{U}_{\rm id}(\tau)$ is evolved into $\bm{U}_{\rm id}^{*}(\tau+\Delta\tau)$, where the asterisk indicates a variable evolved only in the ideal part. $\bm{V}_{\rm id}(\tau)$ is used to evaluate the numerical flux and the geometric source terms in Eqs. (30) and (31). We calculate the fluid variable $\bm{V}_{\rm id}^{*}(\tau+\Delta\tau)$ from $\bm{U}^{*}_{\rm id}(\tau+\Delta\tau)$ with the algorithm for recovery of the primitive variables $\bm{V}_{\rm id}$ from the conserved variables $\bm{U}_{\rm id}$ Akamatsu2014 .

Second, we solve the constitutive equations of the viscous tensors Eqs. (34)-(42) to obtain $\bm{V}_{\rm vis}(\tau+\Delta\tau)$. 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^{\mu}_{\rm NS},\pi^{\mu\nu}_{\rm NS},\Pi_{\rm 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 $\partial_{\tau}\bm{V}_{\rm id}=(\bm{V}^{*}_{\rm id}(\tau+\Delta\tau)-\bm{V}_{\rm id}(\tau))/\Delta\tau.$ Here we keep the middle time-step value of the viscous tensor $\bm{V}_{\rm vis}(\tau+\Delta\tau/2)$ for the next step.

Next, the conserved variables $\bm{U}^{*}_{\rm id}(\tau+\Delta\tau)$ and $\bm{U}_{\rm vis}(\tau)$ are evolved into $\bm{U}_{\rm id}(\tau+\Delta\tau)$ and $\bm{U}_{\rm vis}(\tau+\Delta\tau)$ by the viscous part of conservation equation Eqs. (32) and (33). Then we recover the fluid variables $\bm{V}_{\rm id}(\tau+\Delta\tau)$ from conserved variables $\bm{U}_{\rm vis}(\tau+\Delta\tau)$ Inutsuka2011 . We keep the middle time-step value $\bm{V}_{\rm id}(\tau+\Delta/2)$.

To achieve the second-order accurate in time, we repeat the above whole steps using the middle time-step values $\bm{V}_{\rm id}(\tau+\Delta\tau/2)$ and $\bm{V}_{\rm vis}(\tau+\Delta\tau/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 Van1979 for the second-order accurate in space or the piecewise parabolic method (PPM) Colella1984 ; Marti1996 ; Colella2008 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 Gubser flow regime Marrochio2015 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=T_{0}(\tau_{0}/\tau)^{1/3}$, where $\tau_{0}$ and $T_{0}$ are the initial proper time and temperature, respectively Bjorken1983 . In the viscous fluid, the non-vanishing components of viscous tensor $\pi^{\mu\nu}$ are $\pi^{xx}$, $\pi^{yy}$, and $\pi^{\eta\eta}$. From the symmetries of the system, the relation $2\pi^{xx}=2\pi^{yy}=-\tau^{2}\pi^{\eta\eta}$ 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

[TABLE]

where $\pi^{\eta\eta}_{\rm NS}$ is the Navier-Stokes value of shear tensor,

[TABLE]

If $\eta/s$ is constant, Eqs. (46) and (47) give the time evolution of the temperature,

[TABLE]

The numerical calculation is carried out on the space-grid size $\Delta\eta=0.1$ with the time-step size $\Delta\tau=0.1\tau_{0}\Delta\eta$. The initial temperature $T_{0}$ and the proper time $\tau_{0}$ are set to $T_{0}=300$ MeV and $\tau_{0}=1$ fm, respectively. We set the relaxation time to be $\tau_{\eta}=0.0001$ fm as the Navier-Stokes limit. Since $\tau_{\eta}$ is smaller than $\Delta\tau$, 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_{\Pi}$ in Eq.  (21), we write the relaxation equation of the bulk pressure,

[TABLE]

with the Navier-Stokes value of the bulk pressure $\Pi_{\rm NS}=\zeta/\tau$. If we assume $\zeta$ and $\tau_{\Pi}$ are constant, we obtain the analytical solution of Eq. (49),

[TABLE]

where $\Pi_{0}$ is the initial value of the bulk pressure and ${\rm Ei}(x)$ is the exponential integral function.

In the numerical calculation, we set $\Pi_{0}=0$, $\zeta=1000$ MeV/fm2, $\tau_{0}=1$ fm, and $\tau_{\Pi}=1$ fm. The space-grid size $\Delta\eta=0.1$ and the time-step size $\Delta\tau=0.1\tau_{0}\Delta\eta$ 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 Gubser2010 ; Gubser2011 , a semi-analytic solution of the Israel-Stewart theory in the Gubser flow regime is obtained Marrochio2015 . 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 Marrochio2015 ; Noronha2014 ; Shen2014 ; Pang2014 ; ECHO2015 . The velocity profile of the semi-analytic solution is the same as that of the ideal Gubser flow,

[TABLE]

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 Marrochio2015 . The second-order terms and the relaxation time in Eq. (20) are given by

[TABLE]

where $c$ is a constant Marrochio2015 .

We carry out the numerical calculation with the finite shear viscosity $\eta/s=0.2$. We set the relaxation time to $\tau_{\eta}=5\eta/(Ts)$. The numerical simulation starts at $\tau_{0}=1$ fm. The time-step size and the space-grid size in numerical simulation are set to $\Delta\tau=0.1\Delta x$ and $(\Delta x,\Delta y,\Delta\eta)=(0.05{\rm\hskip 2.84526ptfm},0.05{\rm\hskip 2.84526ptfm},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 $\tau=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 $\tau=7$ fm Okamoto2016 . On the other hand, in the finite viscosity calculation the difference between the numerical calculation and the semi-analytic solution appears after $\tau=4$ fm.

In Fig. 4 the numerical results of the shear tensors $\pi^{xx}$, $\pi^{yy}$, $\pi^{\eta\eta}$ and $\pi^{xy}$ at $\tau=1.2$, 2 and 3 fm are presented together with the semi-analytic solutions. Here the profile of $\pi^{xy}$ is shown along a line $x=y$, since the value of $\pi^{xy}$ vanishes on the $x$ and $y$ axes. The shear tensors $\pi^{xx}$, $\pi^{yy}$ and $\pi^{\eta\eta}$ in our numerical calculations show good agreement with the semi-analytic solutions. However, in $\pi^{xy}$ the deviation from the semi-analytic solution starts to appear at $\tau=2$ fm and grows at later time.

Since in the Israel-Stewart theory the second-order terms in $\pi^{\mu\nu}$ 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.Marrochio2015 , 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 Okamoto2016 . We find that the Corner Transport Upwind (CTU) scheme Colella1990 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 Colella1990 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 $\pi^{xx}$, $\pi^{yy}$, $\tau^{2}\pi^{\eta\eta}$ and $\pi^{xy}$ at $\tau=3$ fm. In $\pi^{xx}$, $\pi^{yy}$ and $\tau^{2}\pi^{\eta\eta}$, results of all numerical schemes are reasonably consistent with the semi-analytic solutions. In addition, differences among them are small. However, in $\pi^{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.

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 Drazin1981 . 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 Csernai2012 . Recently initial fluctuations and QGP expansion not only in the transverse direction but also in the longitudinal direction have attracted interest CMS2015 ; STAR2017a ; STAR2017b . 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,\eta)$ 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=\Delta z$ and $z=-\Delta z$, respectively. Energy density and $\eta$ component of velocity of the flow flux are assumed to be described by Bjorken’s scaling solution $e_{B}=e_{0}(\tau_{0}/\tau)^{4/3}$ and $v^{\eta}_{B}=0$. Shifting the Bjorken scaling solution to $\pm\Delta z(=0.3$ fm) in the $z$ direction, we obtain energy density $e_{\rm U}$ ($e_{\rm D}$) and the $\eta$ component of velocity of the flow flux $v^{\eta}_{\rm U}$ ($v^{\eta}_{\rm D}$) in $x>0$ ($x<0$),

[TABLE]

[TABLE]

where $\tau_{0}$ and $e_{0}$ are the initial time and the energy density, respectively. Figure 6 shows the energy densities $e_{\rm U}$ and $e_{\rm D}$, and the $\eta$ component of velocity of the flow flux $v^{\eta}_{\rm U}$ and $v^{\eta}_{\rm D}$ in $x>0$ and $x<0$. The energy density and the $\eta$ components of velocity of the flow flux are dependent on $\eta$ 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 $x_{b}=0.01{\rm sin}(2\pi\eta/\lambda)$ with a wavelength $\lambda$ along the boundary between the flow fluxes. Finally our initial energy density and flow velocity are written by

[TABLE]

where the energy density and flow velocity around the boundary are connected from $e_{U}$ and $v^{\eta}_{U}$ to $e_{D}$ and $v^{\eta}_{D}$ smoothly with the parameter $\Delta$. Here, we set the wavelength of a fluctuation and the width of boundary between two fluid fluxes to $\lambda=0.4$ and $\Delta=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/fm3($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}$,

[TABLE]

The arrows stand for the velocity field in $(\tau v^{\eta},v^{x})$.

We start the numerical calculation on the grid $(\Delta x,\Delta\eta)=(0.005{\rm\hskip 2.84526ptfm},0.00625)$ at $\tau_{0}=1$ fm with time-step size $\Delta\tau=0.2\Delta 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 $\tau\sim 3$ fm. In Fig. 8 the velocity field and the profile of the vorticity $w^{y}$ at $\tau=$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 $\eta\sim 0.6$ is faster than that at $\eta\sim 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 $\lambda$ to $\lambda>0.5$ in the region $|\eta|<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 $\eta/s=0.01$. Figure 9 shows the numerical results of KH instability at $\tau=4$ and 7 fm. In contrast to Fig. 8, we cannot find the clear vortex but small and vague enhancement of vorticity around $\eta\sim 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 $\lambda=0.4$ at $\tau_{0}=1$ fm may be smaller than the Kolmogorov length scale. In the mid rapidity, $|\eta|<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.

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 Akamatsu2014 and for the viscous part we utilized the PES method Inutsuka2011 . 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 $\tau\sim 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 Fries1 ; Fries2 . 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 ATLAS2015 ; ALICE2016 , event plane correlation ATLAS2014 ; CMS2015 , non-linearity of higher flow harmonics Yan2015 and three particle correlation STAR2017a ; STAR2017b . A comprehensive investigation of experimental data with the accurate numerical method of the relativistic viscous hydrodynamics gives us deep insight of QCD matter.

Acknowledgments

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.

Appendix Appendix A Interpolation procedures

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 $x_{i}$ and the boundary between the $i$th and the $i+1$th 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}={\rm lim}_{x\rightarrow x_{i+1/2}}a(x)$ and $a_{L,i}={\rm lim}_{x\rightarrow 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,

[TABLE]

In the MC limiter Van1979 , $\Delta a_{i}$ is given by

[TABLE]

We define space averages of an interpolation function, $F_{i,R}(\sigma_{i})$ and $F_{i,L}(\sigma_{i})$,

[TABLE]

where $a^{I}(x)$ is an interpolation function of $a(x)$ and $\sigma_{i}=|u_{i}|\Delta t/\Delta 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}(\sigma_{i})$ and $F_{i+1,L}(\sigma_{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}(\sigma_{i})$ or $F_{i+1,L}(\sigma_{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}(\sigma_{i})$ and $F_{i,L}(\sigma_{i})$ are expressed by

[TABLE]

Appendix A.2 Piecewise Parabolic Method (PPM)Colella1984 ; Marti1996 ; Colella2008

First, we calculate interpolated values of $a(x)$ at cell interfaces using forth-order interpolation.

[TABLE]

If the condition ${\rm min}(a_{i},a_{i+1})\leq a_{i+1/2}\leq{\rm max}(a_{i},a_{i+1})$ is not satisfied, $a_{i+1/2}$ is limited as follows:

[TABLE]

If the signs of $(D^{2}a)_{i+1/2},(D^{2}a)_{i+1/2,R}$ and $(D^{2}a)_{i+1/2,L}$ are all the same,

[TABLE]

otherwise, $(D^{2}a)_{i+1/2,{\rm lim}}=0$. Then the modified values of $a_{i+1/2}$ read

[TABLE]

where $C>1$ is a constant. We set $C$ to $C=1.25$ Colella2008 . 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,

[TABLE]

The flattening parameter $f_{i}$ is fixed by $f_{i}={\rm max}(\tilde{f}_{i},\tilde{f}_{i+s_{i}})$, where $s_{j}=+1$ for $p_{i+1}-p_{i-1}>0$ and $s_{j}=-1$ for $p_{i+1}-p_{i-1}<0$,

[TABLE]

The constant $w_{i}$ is chosen by

[TABLE]

The parameters are set to $\epsilon=1$, $w^{(1)}=0.52$, and $w^{(2)}=10$ Marti1996 . 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})\leq 0$ or $(a_{i-1}-a_{i})(a_{i}-a_{i+1})\leq 0$, the $i$th cell contains a local extremum. The values of $a_{i,R}$ and $a_{i,L}$ are modified as follows:

[TABLE]

where $a_{6,i}=6a_{i}-3(a_{i,L}+a_{i,R})$. If $(D^{2}a)_{i}$ and $(D^{2}a)_{i,\{L,C,R\}}$ have the same sign,

[TABLE]

otherwise, $(D^{2}a)_{i,{\rm lim}}=0$. Then we obtain

[TABLE]

If $(D^{2}a)_{i}=0$, we set the second term of Eqs. (82) and (83) to be zero. In the last limiter, the values of $a_{i,R}$ and $a_{i,L}$ are modified as

[TABLE]

The space averages of a parabolic interpolant are written

[TABLE]

Again, $F_{i,R}(\sigma_{i})$ and $F_{i+1,L}(\sigma_{i+1})$ are used 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}(\sigma_{i})$ or $F_{i+1,L}(\sigma_{i+1})$ corresponds to the numerical flux passing through the cell boundary $x_{i+1/2}$ (Appendix B).

Appendix Appendix B Numerical schemes for convection equations

Appendix B.1 High-resolution upwind method

We consider the one-dimensional convection equation,

[TABLE]

In the high-resolution upwind method, we obtain the solution of the convection equation Eq.(88),

[TABLE]

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}+\Delta t$. The numerical flux $a^{n+1/2}_{i+1/2}$ reads

[TABLE]

We evaluate the $F_{i,R}(\sigma_{i})$ and $F_{i,L}(\sigma_{i})$, using the MC limiter (Eqs. (65) and (66)) or the PPM (Eqs.(86) and (87)).

In the case of multidimensional problems, we employ the Strang splitting method Strang1968 . Using the operator $L^{k}_{i}$, which represents one-dimensional evolution in the $i$ direction during the time $k\Delta t$, we express two-dimensional expansion in the $(x,y)$ coordinates as

[TABLE]

Similarly the three-dimensional expansion in $(x,y,z)$ coordinates is written by

[TABLE]

Appendix B.2 Corner transport upwind (CTU) scheme Colella1990

We consider two-dimensional convection equation,

[TABLE]

In the CTU, the solution of the convection equation Eq.(93) reads

[TABLE]

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}+\Delta t$, the second and third terms stand for the numerical flux passing through the cell boundary. The numerical flux is given by

[TABLE]

Here we evaluate the variation of $a(t,x,y)$ in the $x$ ($\Delta^{x}a_{i,j}$) and $y$ direction ($\Delta^{y}a_{i,j}$) using the MC limiter (62).