A new relativistic viscous hydrodynamics code and its application to the Kelvin–Helmholtz instability in highenergy heavyion collisions
Eur. Phys. J. C
A new relativistic viscous hydrodynamics code and its application to the KelvinHelmholtz instability in highenergy heavyion collisions
Kazuhisa Okamoto 2
Chiho Nonaka 0 1 2
0 Department of Physics, Duke University , Durham, NC 27708 , USA
1 KobayashiMaskawa Institute for the Origin of Particles and the Universe (KMI), Nagoya University , Nagoya 4648602 , Japan
2 Department of Physics, Nagoya University , Nagoya 4648602 , 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 twoshock 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 IsraelStewart theory in Gubser flow regime. Using the code, we discuss possible development of the KelvinHelmholtz instability in highenergy heavyion 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)
highenergy heavyion 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 spacetime 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 email:
b email:
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 secondorder relativistic viscous
hydrodynamic equation. In highenergy heavyion collisions,
currently the Israel–Stewart theory [2,3] and conformal
hydrodynamics [4] are often used. Solving them numerically, study
of experimental data of highenergy heavyion 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], nonlinearity
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 twoshock 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 highenergy heavyion
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 highenergy heavyion
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 fourvelocity 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 energymomentum 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 secondorder 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 highenergy
heavyion 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 nonrelativistic fluid usually
has a problem of acausality and instability [32–34]. The
problem can be resolved by introducing the secondorder terms
of the viscous tensor and the derivative of fluid variables
into the hydrodynamic equations [2,3]. However, the
original Israel–Stewart theory does not reproduce the results of the
kinetic equation quantitatively [35–39]. The construction of
secondorder 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 secondorder 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 secondorder 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
fourvelocity 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 threevelocity and the Lorentz factor, respectively. The
conservation equations Eqs. (1) and (2) are explicitly written
as
where i = x , y, η and the righthand 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 energymomentum 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 timestep 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 timestep 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 secondorder terms I , K in the righthand 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 righthand sides of Eqs. (37)–(42) are obtained
by ∂τ V id = (V i∗d(τ + τ ) − V id(τ ))/ τ. Here we keep the
middle timestep 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 timestep
value V id(τ + /2).
To achieve the secondorder accurate in time, we repeat
the above whole steps using the middle timestep 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 secondorder 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 secondorder accurate
in space or the piecewise parabolic method (PPM) [55–57]
for the thirdorder 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 onedimensional
expansion and the Israel–Stewart theory in the Gubser flow
regime [58] for the threedimensional 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 onedimensional 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 nonvanishing 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 spacegrid size
η = 0.1 with the timestep 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 secondorder 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 spacegrid size
η = 0.1 and the timestep 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 semianalytic solution of the Israel–Stewart theory in
the Gubser flow regime is obtained [58]. The semianalytic
solution is a useful test problem for the code of relativistic
viscous hydrodynamics which is developed for application to
the highenergy heavyion collisions [13, 14, 58, 61, 62]. The
velocity profile of the semianalytic 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 secondorder 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 timestep size and the spacegrid 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 semianalytic
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 semianalytic 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 semianalytic 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 semianalytic 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 semianalytic solutions. However, in π x y the deviation
from the semianalytic solution starts to appear at τ = 2 fm
and grows at later time.
Since in the Israel–Stewart theory the secondorder terms
in π μν become small compared with the firstorder 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 semianalytic
solution. Here we employ the PPM for solving the
convection part numerically, instead of the MC limiter. In the case of
threedimensional calculation we use the dimensional
splitting method [24]. We find that the Corner Transport Upwind
(CTU) scheme [63], which is a threedimensional unsplit
method, realizes good agreement of the semianalytic
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 semianalytic solutions and the pluses stand for numerical
results
with the PPM and the CTU method [63] with the MC
limiter for threedimensional unsplit method. We shall explain
the details of each scheme in Appendices A and B. Figure 5
shows the semianalytic 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 semianalytic 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
semianalytic solution at the peak around x = 2 fm appears,
whereas the PPM and the CTU methods keep the good
agreement with the semianalytic solution. The CTU method can
achieve the high numerical accuracy with the secondorder
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 heavyion 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 heavyion collisions, the colorflux 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 heavyion
collisions.
For simplicity, we focus on hydrodynamic expansion in
the (x , η) plane. The heavy ion accelerated with highenergy
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 largex partons is covered by a cloud
of smallx partons. As a result, in the highenergy
heavyion collisions parton–parton interactions may take place in
the area within around 1 fm from z = 0 fm. Then if we
consider the colorflux tube structure in the initial condition,
each colorflux 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 colorflux 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 timestep size
0.2 x . We use the ideal gas equation of state e = 3 p.
’ideal/xz0000.dat’ u 3:1:($8)
wy(fm1)
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 threefluid 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 secondorder 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 highenergy heavyion 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 twoshock
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 onedimensional test and the Israel–Stewart
theory in the Gubser flow regime for the threedimensional 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 highenergy heavyion
collisions. We focused on the mid rapidity and started the
numerical calculations with the simple initial conditions inspired
by the colorflux 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 highenergy heavyion 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], nonlinearity 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 GrantinAid for Scientific Research (S) No. 26220707 and US
Department of Energy grant DEFG0205ER41367.
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 secondorder accurate in
space and the PPM for the thirdorder 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 secondorder accuracy in space is achieved by the linear
interpolation. In the secondorder 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 fourthorder 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: Highresolution upwind method
We consider the onedimensional convection equation,
= 0.
In the highresolution 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 onedimensional evolution in the i direction
during the time k t , we express the twodimensional expansion
in the (x , y) coordinates as
an+1 = L 1x/2 L1y L 1x/2an.
Similarly the threedimensional 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 twodimensional 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 quarkgluon 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 [hepth]
5. H. Song , U.W. Heinz , Phys. Rev. C 77 , 064901 ( 2008 ). doi:10. 1103/PhysRevC.77.064901. arXiv: 0712 .3715 [nuclth]
6. K. Dusling , D. Teaney , Phys. Rev. C 77 , 034905 ( 2008 ). doi:10. 1103/PhysRevC.77.034905. arXiv: 0710 .5932 [nuclth]
7. M. Luzum , P. Romatschke , Phys. Rev. C 78 , 034915 ( 2008 ). doi:10. 1103/PhysRevC.78.034915. arXiv: 0804 .4015 [nuclth]
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 [hepph]
9. B. Schenke , S. Jeon , C. Gale , Phys. Rev. Lett . 106 , 042301 ( 2011 ). doi:10.1103/PhysRevLett.106.042301. arXiv: 1009 .3244 [hepph]
10. P. Bozek , Phys. Rev. C 85 , 034901 ( 2012 ). doi:10.1103/PhysRevC. 85.034901. arXiv: 1110 .6742 [nuclth]
11. V. Roy , A.K. Chaudhuri , Phys. Rev. C 85 , 024909 ( 2012 ). doi:10. 1103/PhysRevC.85.024909. arXiv: 1109 .1630 [nuclth]
12. I. Karpenko , P. Huovinen , M. Bleicher , Comput. Phys. Commun . 185 , 3016 ( 2014 ). doi:10.1016/j.cpc. 2014 .07.010. arXiv: 1312 .4160 [nuclth]
13. J. NoronhaHostler , J. Noronha , F. Grassi , Phys. Rev. C 90 , 034907 ( 2014 ). doi:10.1103/PhysRevC.90.034907. arXiv: 1406 .3333 [nuclth]
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 [hepph]
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 [nuclth]
16. G. Aad et al., ATLAS Collaboration, Phys. Rev. C 92 , 034903 ( 2015 ). doi:10.1103/PhysRevC.92.034903. arXiv: 1504 .01289 [hepex]
17. J. Adam et al., ALICE Collaboration, Phys. Rev. Lett 117 , 182301 ( 2016 ). doi:10.1103/PhysRevLett.117.182301. arXiv: 1604 .07663 [nuclex]
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 [nuclex]
20. L. Yan , J.Y. Ollitrault , Phys. Lett . B 744 , 744 , 82 ( 2015 ). doi:10. 1016/j.physletb. 2015 .03.040. arXiv: 1502 .02502 [nuclth]
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:astroph/ 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 [nuclth]
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 [nuclth]
26. P. Romatschke , Theor. Phys. Suppl . 174 , 137 ( 2008 ). doi:10.1143/ PTPS.174.137. arXiv: 0710 .0016 [nuclth]
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 [nuclth]
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 ( ButterworthHeinemann , Oxford, 1987 )
31. A. Andronic , P. BraunMunzinger , K. Redlich , J. Stachel , Nucl. Phys. A 904 , 535c ( 2013 ). doi:10.1016/j.nuclphysa. 2013 .02.070. arXiv: 1210 .7724 [nuclth]
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 [hepph]
35. P. Huovinen , D. Molnar , Phys. Rev. C 79 , 014906 ( 2009 ). doi:10. 1103/PhysRevC.79.014906. arXiv: 0808 .0953 [nuclth]
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 [nuclth]
38. M. Takamoto , S. Inutsuka , Physica (Amsterdam) 389, 4580 ( 2010 ). doi:10.1016/j.jcp. 2011 .05.030. arXiv: 1106 .1732 [astroph]
39. A. El , Z. Xu , C. Greiner , Phys. Rev. C 81 , 041901 ( 2010 ). doi:10. 1103/PhysRevC.81.041901. arXiv: 0907 .4500 [hepph]
40. A. Muronga , Phys. Rev. C 76 , 014910 ( 2007 ). doi:10.1103/ PhysRevC.76.014910. arXiv:nuclth/0611091
41. B. Betz , D. Henkel , D.H. Rischke , J. Phys . G 36 , 064029 ( 2009 ). doi:10.1051/epjconf/20111307005. arXiv: 1012 .5772 [nuclth]
42. G.S. Denicol , T. Koide , D.H. Rischke , Phys. Rev. Lett . 105 , 162501 ( 2010 ). doi:10.1103/PhysRevLett.105.162501. arXiv: 1004 .5013 [nuclth]
43. E. Calzetta , J. PeraltaRamos , Phys. Rev. D 82 , 106003 ( 2010 ). doi:10.1103/PhysRevD.82.106003. arXiv: 1009 .2400 [hepph]
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 [nuclth]
45. A. Jaiswal , Phys. Rev. C 87 , 051901 ( 2013 ). doi:10.1103/ PhysRevC.87.051901. arXiv: 1302 .6311 [nuclth]
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 [hepth]
47. M. Natsuume , T. Okamura , Phys. Rev. D 77 , 066014 ( 2008 ). doi:10. 1103/PhysRevD.77.066014. arXiv: 0712 .2916 [hepth]
48. P. Romatschke , Class. Quantum Gravity 27 , 025006 ( 2010 ). doi:10. 1088/ 0264  9381 /27/2/025006. arXiv: 0906 .4787 [hepth]
49. K. Tsumura , T. Kunihiro , Eur. Phys. J. A 48 , 162 ( 2012 ). doi:10. 1140/epja/i2012 12162 x. arXiv: 1206 . 1929 [nuclth]
50. K. Tsumura , Y. Kikuchi , T. Kunihiro , Phys. Rev. D 92 , 085048 ( 2015 ). doi:10.1103/PhysRevD.92.085048. arXiv: 1506 .00846 [hepph]
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 [astroph]
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 [nuclth]
59. S.S. Gubser , Phys. Rev. D 82 , 085027 ( 2010 ). doi:10.1103/ PhysRevD.82.085027. arXiv: 1006 .0006 [hepth]
60. S.S. Gubser , A. Yarom , Nucl. Phys . B 846 , 469 ( 2011 ). doi:10. 1016/j.nuclphysb. 2011 .01.012. arXiv: 1012 .1314 [hepth]
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 [nuclth]
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 [nuclth]
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 [nuclth]
66. R. J. Fries , G. Chen , S. Somanathan , private communication