Shock propagation and stability in causal dissipative hydrodynamics
aa r X i v : . [ h e p - ph ] M a y Shock propagation and stability in causal dissipative hydrodynamics
G.S.Denicol, T. Kodama, T. Koide, and Ph. Mota
Instituto de F´ısica, Universidade Federal do Rio de Janeiro, C. P. 68528, 21945-970, Rio de Janeiro, Brasil
We studied the shock propagation and its stability with the causal dissipative hydrodynamics in1+1 dimensional systems. We show that the presence of the usual viscosity is not enough to stabilizethe solution. This problem is solved by introducing an additional viscosity which is related to thecoarse-graining scale of the theory.
PACS numbers: 47.10.-g,25.75.-q
I. INTRODUCTION
It is by now widely accepted that the basic features ofcollective motions in relativistic heavy-ion collisions canbe well described as those of a hydrodynamical motion ofan (almost) ideal fluid [1]. Several studies on the effectsof viscosity are available in literatures and support sucha vision [2, 3, 4, 5]. In addition, the possibility of theexistence of the lower bound of the shear viscosity co-efficient has been discussed by assuming the AdS/CFTcorrespondence [6].However, strictly speaking there are still several openquestions in hydrodynamic approaches of the heavy-ioncollisions. Hydrodynamic observables are not so restric-tive to determine uniquely many unknown factors, suchas initial condition, equation of state, dissipation mech-anisms, and so on. In particular, even the theory of arelativistic dissipative fluid is not yet well understood.Therefore, in order to conclude that the matter createdin relativistic heavy-ion collisions really behaves just asan ideal fluid, we need to investigate the effect of dissi-pation more carefully.The difficulty of the construction of the relativisticdissipative hydrodynamics is because of the problem ofacausality and instability: a naive relativistic general-ization of the Navier-Stokes equation present an infinitepropagation speed of pulse signals and the solution is un-stable for small perturbations. To solve this problem, oneshould, for example, take into account memory effects byintroducing a relaxation time [5].So far, there are several attempts to study the effectof dissipation to relativistic heavy-ion collisions by im-plementing the 1+1 and 2+1 dimensional calculations[4]. In these studies, they mainly deal with cases wherethe deviation of the ideal fluid is rather small. On theother hand, there is an evidence that the bulk viscos-ity becomes large in the vicinity of the critical point[7].In addition, to know the limitation of the theory itself,we should investigate the behavior of solutions for largeviscosity as well. In particular, for higher energies suchas in the LHC regime, we expect that the relative im-portance of the viscosity becomes significant. This pointis essential since the problem of acausality or instabili-ties is directly related with the size of viscosity and therelaxation time mentioned above.Other interesting aspects of viscous fluid dynamics ap- pear when discontinuities emerge during the time evolu-tion or already exist in the initial condition. In the usualapplication of hydrodynamics, only very smooth initialdistributions, both in energy and velocity, have been ap-plied. In such cases, the fluid profile usually remainssmooth in time and no special attention is required forthe treatment of discontinuities.However, it is sometimes necessary to discuss the ex-treme cases which involve discontinuities. For example,the Landau type initial condition is often discussed as aninteresting possibility of meson production mechanismwith null initial velocity in the p-p collisions and A-Acollisions [8]. Furthermore, the shock phenomena, whichare typical discontinuous propagation of hydrodynamicalvariables, also may occur in the heavy-ion collisions byhigh energy jet propagations in the QGP [3, 9]. We alsoexpect a formation of shock wave in the region near thecoexisting phase, if the QCD-hadron phase transition isof first order, where the velocity of sound vanishes (orbecomes very small). Shock phenomena in relativisticheavy ion collisions, if any, are particularly interestingsince they would furnish genuine hydrodynamical signals.It is well-known that to deal with dynamical discon-tinuities, such as shocks, is not a simple problem. Wehave to introduce some specific techniques such as theGodunov method or artificial viscosity (pseudo-viscosity)[10] to achieve physically meaningful results interpolatingsmoothly the discontinuity. For a causal dissipative hy-drodynamics, there are a few works where the dynamicsof shock discontinuities are discussed [11, 12]. However,to the author’s knowledge, detailed numerical study of acausal dissipative fluid dynamics involving discontinuitieshas not yet been done.In this paper, we study in detail the dynamics of a vis-cous fluid. For the sake of simplicity, we concentrate our-selves to 1+1 dimensional systems. The basic objectiveis to analyze the problem of instabilities associated withcausality through several numerical examples, in partic-ular, which contain discontinuities or shock phenomena.We also pay attention to the origin of the so-called “ar-tificial viscosity” for the numerical calculations of shockphenomena in the framework of causal viscous hydrody-namics. We investigate its role and interpret its origin interms of the scale of the coarse graining introduced in thehydrodynamical theory. To clarify this point, we use thesmoothed particle hydrodynamics (SPH) as the numer-ical method instead of the commonly used space-fixedgrid methods. Furthermore, having in mind the collec-tive QGP motion in the LHC energy regime, we restrictourselves to systems described by an equation of state ofa baryon-number free, massless relativistic gas.In the next section, we briefly review the causal dis-sipative hydrodynamics and describe how the memoryfunction is introduced. To check causality and stabilityof our theory, we discuss the dispersion relation for thepropagation of a perturbative plane-wave. In Sec. III, wedescribe the SPH method applied to our problem whichintroduces the coarse-graining scale h . In Sec. IV, weshow several examples to reveal the effects of viscosity,in particular, the case of steady shock wave propagationinduced by a large pressure gradient in the initial condi-tion. In the first example, we show the universal relationbetween pressure and viscosity in the fluid expansion tovacuum. This relation is satisfied independently of ini-tial conditions and equation of states. In these examples,some quick oscillating modes appear in the dynamics in-volving shocks, leading to instabilities of the numericalsolutions. That is, the normal viscosity does not neces-sarily damp all the high frequency modes. The presenceof such modes indicates the necessity of a new ingredientfor the theory to be physically consistent. In Sec. V, weintroduce an additional viscosity, having a different scaleof the memory function associated with the coarse grain-ing size, which solves this problem. We show severalexamples where the dynamics of shock phenomena aredescribed satisfactorily in this scheme. Finally in Sec.VI,we summarize our work on the analysis of a causal dissi-pative hydrodynamics in 1+1 and discuss the problemsstill open in such theories. II. RELATIVISTIC DISSIPATIVEHYDRODYNAMICS
Various theories have been proposed to incorporate dis-sipation consistent with causality and stability; the di-vergence type theory [13], the Israel-Stewart theory [14]and its extension based on the extended irreversible ther-modynamics [15], Carter’s theory [16], ¨Ottinger-Grmelaformulation [17] and the memory function method [5].Here we briefly review the memory function method[5] to obtain a causal dissipative hydrodynamics for 1+1dimensional case. That is, we ignore the motion in thetransverse direction and concentrate on the longitudinaldynamics.As was mentioned in the introduction, we just considerthe case of vanishing baryon chemical potential. In thiscase, the hydrodynamical equation of motion can be writ-ten only as the conservation of the energy-momentumtensor, ∂ µ T µν = 0 . (1)together with the thermodynamical relations. We use theLandau definition for the local rest frame and assume, as usual, that the thermodynamic relations are valid in thisframe. Then the energy-momentum tensor is expressedas T µν = ( ε + p + Π) u µ u ν − ( p + Π) g µν , (2)where, ε , p , u µ and Π are, respectively, the energy den-sity, pressure, four velocity and bulk viscosity.In the presence of these irreversible currents, the en-tropy is not conserved anymore. Instead, from Eq.(1),we have [18] ∂ µ σ µ = − T Π ∂ µ u µ , (3)where the entropy four-flux is identified by Landau andLifshitz as σ µ = su µ . (4)In irreversible thermodynamics, it is interpreted thatentropy production is the sum of the products of thermo-dynamic forces and irreversible currents. From Eq. (3),we define the thermodynamic force as F = ∂ α u α . (5)To satisfy the second law of thermodynamics locally, andhence the positiveness of the entropy production, Landauproposed that the irreversible current should be propor-tional to the thermodynamic force [18],Π = − ζF = − ζ∂ α u α , (6)where ζ is the viscosity coefficient. However, it is knownthat the derived equations have the problem of acausalityand instability [19, 20, 21]. To solve these difficulties, weintroduce a memory effect to the irreversible current byusing a memory function[5]. One of the simplest formsof the memory function is G ( τ, τ ′ ) → τ R ( τ ′ ) e − R ττ ′ τR ( τ ′′ ) dτ ′′ . (7)The relaxation time τ R ( τ ) is, in general, a function of thelocal proper time τ = τ ( ~r, t ) through the thermodynam-ical quantities. Then, the irreversible current is modifiedas follows;Π ( τ ) = − Z ττ dτ ′ G ( τ, τ ′ ) ζ∂ α u α ( τ ′ ) + e − ( τ − τ ) /τ R Π , (8)where Π is the initial value given at τ .Because of the modification of the relation between theirreversible currents and the thermodynamic forces, thealgebraic positivity of the second law of thermodynamicsis not satisfied. However, we checked that the second lawof thermodynamics is not violated in all examples dis-cussed in this paper numerically. It is worth mentioningthat the algebraic positivity is not automatically satis-fied even in the Israel-Stewart theory (See the discussionbelow Eq. (2.31) in [14]). To solve this problem, the con-cepts of thermodynamics have to be extended. See [15]for details.The integral expression (8) are equivalent to the fol-lowing differential equation,Π = − ζ∂ α u α − τ R d Π dτ , (9)where d/dτ = u µ ∂ µ is the total derivative with respect tothe proper time. In the 1+1 dimensional case discussedhere, the equations derived above are equivalent to thoseof the Israel-Stewart theory. Thus, the conclusions in thispaper is applicable also to the Israel-Stewart theory. A. Propagation speed of signals
We discuss the propagation speed of the 1+1 dimen-sional system. For this purpose, we consider a small per-turbation of the three independent variables in the formof a plane wave, δǫδU δ Π ∝ e iωt − ikx , propagating in a hydrostatic equilibrated background.Then, the linearized hydrodynamic equation for theseperturbations should satisfy iω − ik ( ǫ + p ) 0 α ( − ik ) iω ( ǫ + p ) − ik − ikζ τ R γiω δǫδU δ Π = 0 , (10)where α ≡ dp/dε is square of the velocity of sound in theabsence of the viscosity.To have the non-trivial solution for the perturbation,the determinant of the 3 × ω should satisfy the fol-lowing dispersion relation, ω − αk = i ζ ( ǫ + p ) ωk iτ R ω . (11)whose solution for ω can be written as ω = x + A x + i τ R , where x = ( − i ) / vuut B s(cid:18) B (cid:19) + (cid:18) A (cid:19) , (12) A = (cid:18) b + α (cid:19) k − τ R , (13) B = 13 τ R (cid:18) α − b (cid:19) k + 2(3 τ R ) , (14) FIG. 1: The real part of the frequency ω as a function of k .There are two propagating modes (dashed lines A and B) andone non-propagating mode (solid line C).FIG. 2: The imaginary part of the frequency ω as a functionof k . The two propagating modes A and B shown in Fig.1 are degenerated (dashed lines). The solid line is the non-propagating mode C. with 1 b ≡ ζτ R ε + p . (15)The asymptotic forms of the dispersion relation forlarge k are ω = ( ± k q b + α + i τ R (1+ αb ) + O (cid:0) k − (cid:1) i αbτ R (1+ αb ) + O (cid:0) k − (cid:1) , (16)whereas for small k , ω = ( ± k √ α − ( α +1 / b ) τ R k /b ( − τ R k /b ) + ik τ R ( b − τ R k ) i/τ R . (17)We have three solutions; two of them are propagatingmodes and the remaining one is a non-propagating mode( ω pure imaginary). From Eqs.(16) and (17) we can seethat the sound velocity of propagating modes reduces tothat of the ideal one for small k , on the other hand, forlarge k , the sound velocity is given by p /b + α . In Figs.1 and 2, we show, respectively, the real and imaginaryparts of ω from these dispersion relations as functions ofmomentum k for a = 0 . , b = 6 , and the temperature T =200 M eV . In Fig. 1, the lower two curves (solid and dot-dashed lines) correspond to the upper line of Eq.(16) andthe horizontal line (dotted) corresponds to the second lineof Eq.(16). As we see from Fig. 2, the imaginary part of ω is always positive for all k and converges to constantvalues. Positivity of the imaginary part guarantees thestability of the plane-wave perturbation.We assume that the propagation of physical quantitiesfor the propagating modes are characterized by the groupvelocity. Then, from this relation, the maximum velocityin such a theory is determined as v M = p /b + α. (18)It should be noted that Eq. (18) gives the velocity ofsound in the causal hydrodynamics. As a matter of fact,in the vanishing ζ , this definition coincides with the ve-locity of sound of ideal fluid. It is clear that the groupvelocity becomes infinite at the limit of τ R →
0. In sucha situation, always some portion of the matter tries topropagate with velocity larger than the velocity of thelight for any initial condition. However, since the hy-drodynamic equation is covariant, the presence of thelight-cone singularity forbids such a propagation. Thenthe matter tends to accumulate at the light-cone and insuch situation the linearized wave analysis in the homo-geneous back ground at rest will breakdown. This kindof conflict between the causality and relativistic covari-ance leads eventually to the instabilities of the solutionnear the light-cone [19, 20, 21].As for the non-propagating mode, we find that theimaginary part becomes a constant for large k. That is,the large k component just damps exponentially. Thisis different from the case of a diffusion process wherethe non-propagating mode behaves as k in the large k limit, which leads to an infinite propagation speed. Inthis sense, our theory is causal. B. Parameters
There are so far two approaches to estimate the trans-port coefficients; kinetic approach and microscopic ap-proach. The calculation of the bulk viscosity coefficients are much involved in the kinetic approach. We cannot usethe Boltzmann equation since it contains only the infor-mation of two-body collisions. We have to use the (mod-ified) Enskog model or the Bogoliubov-Cho-Uhlenbeckequation where the multiple collision effects are included.In the microscopic approach, it is known that the trans-port coefficients are calculated by using the Green-Kubo-Nakano (GKN) formula. However, the GKN formula isthe formula for the Newtonian fluid like the Navier-Stokesfluids, and hence we cannot use for the transport coef-ficients of the causal dissipative hydrodynamics becauseit is a non-Newtonian fluid. To calculate the transportcoefficients of the causal dissipative hydrodynamics, wehave to derive a new formula. One possibility of sucha new formula is proposed in [22, 23, 24]. In any case,no reliable estimate for the bulk viscosity coefficient isavailable.In the present analysis, we will not deal with a precisequantitative description of the behavior of the mattercreated in heavy-ion collisions but rather interested inqualitative role of viscosity. Thus, we assume a simpleexpression for the bulk viscosity coefficient ζ as usuallyadopted for the shear viscosity, ζ = as, where s is the entropy density in the local rest frame.Another important parameter of the theory is the re-laxation time τ R . As was mentioned in the previous sec-tion, causality constraints the relation between the twoparameters, ζ and τ R through v M ≤ b as constant. This determines the relaxation timeas τ R = ζǫ + p b. The later examples are presented in terms of these pa-rameters a and b . However, there is no theoretical reasonthat a and b are constants, but they may depend on ther-modynamical quantities (see the later discussion). In thefollowing calculations, we use two values for the param-eter a , 0 . , with fixed value of b = 6. III. SMOOTHED PARTICLESHYDRODYNAMICS
To solve numerically the relativistic hydrodynamicequations we use the Smoothed Particle Hydrodynamic(SPH) method. This method was initially introduced[25] for application in astrophysics. Using variational ap-proach, this method was extended to the application forheavy ion collisions [26].The original idea of the SPH method is to obtain anapproximate solution of hydrodynamics by parameter-izing the fluid into a set of “effective particles”. How-ever, in the application to the heavy ion dynamics, theSPH method is not a mere mathematical discretizationscheme, but can be interpreted as a physical model ofthe collective motion in terms of a finite set of dynamicalvariables.To see this, let us consider a distribution a ( r , t ) of anyextensive physical quantity, A . In a system like the hotand dense matter created in heavy ion collisions, the be-havior of a ( r , t ) contains the effects of whole microscopicdegrees of freedom. We are not interested in the ex-tremely short wavelength behavior of a ( r , t ) but rather inglobal behaviors which are related directly to the exper-imental observables. Therefore, we would like to intro-duce a coarse-graining procedure for a . To do this, we in-troduce the kernel function W ( r − ˜r , h ) which maps theoriginal distribution a to a coarse-grained version a CG as, a CG ( r , t ) = Z a ( ˜r , t ) W ( r − ˜r , h ) d ˜r (19)where W is normalized, Z W ( ˜r , h ) d ˜r = , (20)and has a bounded support of the scale of h , W ( r , h ) → , | r | & h, satisfying lim h → W ( ˜r , h ) = δ ( ˜r ) . Here, h is a typical length scale for the coarse-grainingin the sense that the kernel function W introduces a cut-off in short wavelength of the order of h . Thus we willtake this value as the scale of coarse graining in the QCDdynamics (i.e., the mean-free path of partons) to obtainthe hydrodynamics of QGP ( h ≃ . f m ).The second step is to approximate this coarse graineddistribution a CG ( r , t ) by replacing the integral in Eq.(19)by a summation over a finite and discrete set of points, { r α ( t ) , α = 1 , .., N SP H } ,a SP H ( r , t ) = N SPH X α =1 A α ( t ) W ( | r − r α ( t ) | ) . (21)If the choice of { A α ( t ) , α = 1 , .., N SP H } and { r α ( t ) , α = 1 , .., N SP H } are appropriate, the aboveexpression should converge to the coarse-graineddistribution a CG for large N SP H . Parameters { A α ( t ) , α = 1 , .., N SP H } and { r α ( t ) , α = 1 , .., N SP H } should be determined from the dynamics of the system.In practice, we first choose the reference density σ ∗ which is conserved, ∂σ ∗ ∂t + ∇ · j = 0 , (22) where ~j is the current associated with the density σ ∗ . Then, we note that the following ansatzs, σ ∗ SP H ( r , t ) = N SPH X α =1 ν α W ( | r − r α ( t ) | ) , j SP H ( r , t ) = N SPH X α =1 ν α d r α ( t ) dt W ( | r − r α ( t ) | ) , satisfies the equation, ∂σ ∗ SP H ∂t + ∇ · j SP H = 0 , where ν α ´ s are constant. By using the normalization of W, Eq.(20), we have Z SP H σ ∗ ( r , t ) d r = N SPH X α =1 ν α . Then we can interpret the quantity ν α as the conservedquantity attached at the point r = r α ( t ). Therefore,the distribution σ ∗ SP H ( r , t ) is a sum of small piece-wisedistribution, carrying the density, ν α W ( | r − r α ( t ) | ) . These pieces are referred to as ”SPH-particles”.Using the above reference density and the extensivenature of A, we can write A α in Eq.(21) as A α ( t ) = ν α a ( r α , t ) σ ∗ ( r α , t )which represents the quantity A carried by the SPH par-ticle at the position r = r α ( t ). In fact, the total amountof A of the system at the instant t is given by A ( t ) = N SPH X α =1 A α ( t ) . In the previous works the entropy density is chosen asthe reference density and the dynamics of the parameters { r α ( t ) , α = 1 , .., N SP H } are determined from the varia-tional principle from the action of ideal hydrodynamics[26]. Thus, the SPH particle coordinates and their timederivatives are considered as the variational parameterswhich optimize dynamically the action of the system.The entropy density is, however, conserved only forthe motion of an ideal fluid. When we discuss the behav-ior of viscous fluids, where the entropy is not conserved,we cannot use the entropy density as the reference SPHdensity. Thus we introduce a new conserved quantity, thespecific proper density σ , which is defined by the flow ofthe fluid, ∂ µ ( σu µ ) = 0 , (23)and we will use it as the reference density for viscousfluids. Here, the four-velocity u µ is defined in terms ofthe local rest frame of the energy flow (Landau frame).The specific density is expressed in the SPH form as σ ∗ ( r , t ) = N SPH X α =1 ν α W ( | r − r α ( t ) | ) , where σ ∗ = σu is the specific density in the laboratoryframe and ν α is the inverse of the specific volume of theSPH particle α, and is chosen as an arbitrary constant.Final results do not depend on this choice and we set ν α = 1 for simplicity. As for the kernel W ( r ), we use thespline function [25].Strictly speaking, this procedure is only possible pro-vided that the lines of flow in space defined by the veloc-ity field u µ do not cross each other during the evolutionin time. That is, if there appear turbulences or singular-ities in the flow lines, the above definition of Lagrangecoordinates fails. However, if the size of h is consistentlychosen as the size of coarse graining of the underlyingmicroscopic theory, the flux line calculated using this h should not cross.Now we apply this method to the causal dissipativehydrodynamics in 1 + 1 dimension. We have to solve theevolution equation of the viscosity in the SPH scheme.To do so, we express the viscosity asΠ = N SPH X α =1 ν α Π α σ ∗ α W ( | r − r α ( t ) | ) , (24)Time evolution of the term Π α can be calculated as γ α d Π α dt = − ζτ R ( ∂ µ u µ ) α − τ R Π α (25)where γ α is the Lorentz factor of the α − th particle. Atthe same time, using the SPH expression for the entropydensity s ∗ in the observable frame, s ∗ = N SPH X α =1 ν α (cid:16) sσ (cid:17) α W ( | r − r α ( t ) | ) , (26)and using Eq.(3) we find, ddt (cid:16) sσ (cid:17) α = − T Π α σ ∗ α ( ∂ µ u µ ) α . where s = s ∗ /u is the proper entropy density. In the fol-lowing, we denote the quantity in the observable framewith the asterisk. In the above expressions, relaxationtime τ R , viscosity coefficient ζ and temperature T arefunctions of space and time, so that they should be eval-uated at the position of each particle α. Finally, we need to express the momentum conserva-tion equation by the SPH variables. We write the spacecomponent of Eq.(1) in terms of the reference density, σ ddτ (cid:18) ǫ + p + Π σ u i (cid:19) + ∂ i ( p + Π) = 0 . (27) It should be noted that there exist ambiguities withinthe resolution of the coarse-graining size h to express theequation of motion in the SPH form. However, in theideal fluid, the SPH equation of motion can be derivedby the variational method uniquely. Thus, we obtain theequation of motion by using the same SPH parametriza-tion to Eq.(27) , σ α ddτ α (cid:18) ǫ α + p α + Π α σ α u iα (cid:19) = N SPH X β =1 ν β σ ∗ α p β + Π β (cid:16) σ ∗ β (cid:17) + p α + Π α ( σ ∗ α ) ∂ i W ( | r α − r β ( t ) | ) , (28)where the right hand side of Eq.(28) corresponds to theterm ∂ i ( p + Π) written in terms of the SPH parametriza-tion. We remark that in the case of vanishing viscosityour result is reduced to the expression derived with vari-ational principle for ideal fluids.By separating the acceleration and force terms inEq.(28), we obtain our final expression of the equationof motion for the SPH particles, M α d~u α dt = ~F α , where M ij = γ ( ǫ + p + Π) δ ij + Au i u j , (29) F j = − ∂ j ( p + Π) + Bu j , (30)with A = 1 γ (cid:20) ε + p + Π − ∂∂s ( ε + p ) (cid:18) s + Π T (cid:19) − ζτ R (cid:21) , (31) B = A γ σ ∗ dσ ∗ dt + Π τ R . (32)These set of equations define the coarse-grained dynamicsto represent the continuity equation for the energy andmomentum tensor of a relativistic fluid, together with anirreversible mechanism which converts a part of collectivekinetic energy (the motion of SPH particles) into internalheat of the fluid.For the sake of book-keeping, we summarize the prac-tical algorithm of calculating the SPH method to solvenumerically the causal dissipative hydrodynamics. Thedynamics is described by the following variables corre-sponding to the quantities attached to the each SPH par-ticle: n r α , u α , ν α , (cid:16) sσ (cid:17) α , Π α ; α = 1 , .., N SP H o . At the initial time, their values are determined accord-ing to the initial condition. The entropy density profileand the bulk viscosity are then obtained with the in-terpolations in Eq.(26) and Eq.(24), respectively. Theenergy density, pressure and temperature are calculatedwith the equation of state. The time evolution of thesequantities are calculated by solving the equations derivedpreviously, d~u α dt = M − α ~F α ,γ α d Π α dt = − ζτ R ( ∂ µ u µ ) α − τ R Π α ,ddt (cid:16) sσ (cid:17) α = − T Π α σ ∗ α ( ∂ µ u µ ) α Here the inverse matrix of M is calculated explicitly as M − = 1 γ ( ǫ + p + Π) b − Aγ ( ǫ + p + Π) ( γ ( ǫ + p + Π) + ( γ − A ) ~u ~u T . and the four-divergence of the velocity is calculated as( ∂ µ u µ ) a = − γ α σ ∗ α (cid:18) dσ ∗ dt (cid:19) α + 1 γ α ~u α · d~u α dt with dσ ∗ dt = 1 σ ∗ ~j · ∇ σ ∗ − ∇ · ~j. IV. EXAMPLESA. Expansion to the vacuum and stationaryboundary
Let us consider first the Landau type initial condition.For the ideal case, this example has already been dis-cussed in [27] and the SPH scheme works very well. Whenwe introduce the viscosity, we found that the sharp dis-continuity in the boundary leads to undesirable instabil-ities. As we will discuss later, the origin of these insta-bilities is due to the presence of a space discontinuity.When we relax such a steep boundary by replacing theinitial distribution consistent with the SPH scale used,the above instabilities disappears. In this subsection, weuse basically the Landau type initial distribution withthe temperature 200
M eV , distributed uniformly withinthe range x ∈ [ − , f m. To relax the sharp boundaries,we add the surface thickness of 10 h, where h = 0 . f m .For simplicity, we take both vanishing initial velocity andviscosity, u ( x,
0) = 0 , Π( x,
0) = 0 . As we mentioned in the introduction, we use the equationof state that of massless ideal gas, p = ǫ , (33) FIG. 3: The time evolutions of the proper entropy density forideal fluid (dotted line) and the viscous fluid of a = 0 . which in turn, ǫ = Cs / , T = 43 Cs / , where C is a constant related to the Stephan-Boltzmannconstant of massless 3 flavor quark-gluon gas.The temporal evolution of the density profile for theentropy at t = 5 ,
10 and 15 fm is shown in Fig.3. Thesolid line represents the results of the causal hydrody-namics with small viscosity, a = 0 . b = 6. For thesake of comparison, we show the time evolution of thesame initial condition for the ideal fluid with the dashedlines.For the massless ideal fluid, we know that the prop-agation speed into the vacuum should be the speed oflight. On the other hand, one can see that the propaga-tion speed into the vacuum for the viscous fluid is slowerthan that of the ideal fluid. This is because the viscosityacts as an attractive interaction during the expansion ofthe fluid. Thus it takes more time to achieve the speedof light in the causal dissipative hydrodynamics. On theother hand, we can also observe that the propagation ofrarefaction wave into the matter of the viscous fluid isfaster than that of the ideal case. This is what we expectfrom Eq.(18).Interestingly enough, we observe that the behavior ofthe boundary of the viscous fluid seems to be a kind ofa stationary wave. To understand this phenomenon, letus introduce the relative coordinate y = x − v s t witha velocity of the stationary wave v s and suppose thatthe energy-momentum tensor depends only on y near theboundary. Then, from the equation of continuity of the FIG. 4: The time evolution of the pressure P (dotted) andviscosity Π (solid) of the fluid. At the boundary, we can seethe relation P = − Π. The in-set shows the details of theregion where the two curves starts to coincide. energy-momentum tensor, we have ddy [ − v s T + T ] = 0 , (34) ddy [ − v s T + T ] = 0 . (35)where, T µν = ( ε + P + Π) u µ u ν − g µν ( p + Π) . and u µ = ( γ, γv ) . From the boundary condition at thevacuum where T µν vanishes, we get the following twoconditions to define the stationary wave, v s ( T + T ) = T (1 + v s ) , (36) T = v s T . (37)From Eq. (36), we find v s = v or 1 /v . The first solutionis consistent of our stationary wave assumption, showingthat the fluid propagates to the vacuum with the veloc-ity v . The second solution leads to a trivial situation ǫ = P = 0. When we substitute v s = v into Eq.(37) weget P = − Π. Thus, if the boundary of the fluid prop-agates as a stationary wave, then the pressure and theviscosity should satisfy this relation. In fact, this relationcorresponds to the case where the acceleration vanishesfor the fluid motion. In Fig.4, we plotted the profiles of P and Π and we can see clearly the relation P = − Π issatisfied near the boundary (see the in-set). In this fig-ure, we notice that there appear quick oscillation modeswith the wavelength of the order of few h. We will discussthis point in the following sections.
FIG. 5: The shock formation in the proper entropy profile(solid) of the viscous fluid of a = 0 .
1, starting from the dis-continuous initial condition (dotted).
We have checked for various combinations of parame-ters and we found that the above feature of the appear-ance of a stationary wave near the boundary is universalfor the free expansion of viscous matter into the vacuum.This fact is very important because the viscosity Π isalways the same order as the pressure P near the bound-ary to the vacuum, showing that the viscous effect cannot be treated as a small correction to the equilibriumthermodynamical quantities. B. Density discontinuity and shock wave
One interesting question is whether the dissipative hy-drodynamics can describe the formation of a shock wave.We know that for an ideal fluid, when a shock is formed,we need the so-called artificial viscosity for smoothingthe shock region. One might argue that when the realviscosity is present, we do not need to introduce such anartificial viscosity. In Fig. 5, we show the time evolutionof the entropy density of a viscous fluid of a = 0 . h , where we see a quick oscilla-tion of large amplitude. It is interesting to observe thatthe wavelength of these oscillations is exactly the orderof h of the SPH scheme.The physical reason for these oscillations of the SPH FIG. 6: The behavior of SPH solution for the ideal fluid inthe presence of shock wave corresponding to Fig. 5. particles is clear. From the well-known Hugoniot-Ranking relation, we know that at the shock front, thereshould exist a production of entropy [18]. However, for anideal fluid, the total entropy is conserved so that the SPHparticles carry the extra kinetic energy corresponding tothe entropy production at the shock. These extra oscilla-tions propagate with the smallest wavelength (order of h ). The presence of viscosity can damp these oscillationsand turns the kinetic energy of the SPH particles to in-ternal heat of the fluid, recovering the Hugoniot-Rankingrelation. This is the case of the example in Fig. 5.However, this is not always the case. For example,when a becomes large, the relaxation time τ R increases,because of causality. Then the time scale for the dampingbecomes comparable to the evolution time scale generat-ing an oscillatory behavior. Such example is shown inFig. 7 with a = 1. In addition, when we look preciselythe density profile shown in the example of Fig. 5, wenote the existence of quick oscillation modes with thewavelength of the order of h, although in this example,their amplitude does not increase in time.In the SPH calculation, by assumption, h is the scale ofcoarse graining and the resulting dynamics should alwayshave larger wavelength than this. Thus, it is clear thatthe appearance of such rapidly oscillating modes of shortwavelength of the order of h is a signature of an inconsis-tency of the theory with this finite coarse-graining scale.Such a situation also happens in a simple expansion ofthe fluid into the vacuum in the presence of viscosity. Asmentioned in the beginning of this section, if we take areal sharp discontinuous Landau initial condition, sim-ilar rapidly oscillating modes appear at the expandingboundary. This is because a real discontinuity is not com-patible with a finite coarse-graining scale h. See also thesmall oscillations in Fig. 4. When we discuss the shock
FIG. 7: The shock formation in the proper entropy profile(solid) of the viscous fluid of a = 1, starting from the discon-tinuous initial condition (dotted). phenomena using an ideal fluid, then the shock front hasnull thickness and it should be treated as a singularityof the theory. This is because the usual hydrodynamicsassume a null coarse-graining scale.Note that the usual viscosity is also obtained assuminga vanishing coarse-graining scale. To be consistent withthe coarse-graining scheme of the theory, there shouldbe an additional dissipation mechanism which is relatedto the coarse-graining scale. Such a mechanism shouldvanish in the limit of h/L ≪ , where L is the typicalobservable scale of the dynamics. V. VISCOSITY ASSOCIATED WITH COARSEGRAINING SCALE
As discussed above, we need an additional viscosity,whose scale is determined by h. Such a viscosity alsoshould obey the requirement of causality. Thus we pro-pose Π
T ot = Π + Π h , where as before τ ( h ) R d Π h dτ = − Π h − ζ ( h ) ∂ µ u µ with ζ ( h ) = a ( h ) ( ε + P ) h,τ ( h ) R = ζ ( h ) ǫ + p b ( h ) . The coefficients a ( h ) and b ( h ) are the numbers of the or-der of unity. Strictly speaking, these numbers should be0determined from a microscopic theory by incorporatingthe effect of the coarse-graining scale appropriately. Herewe take them as phenomenological parameters and foundthat the values, a ( h ) = 1 / , b ( h ) = 2 , (38)can eliminate undesirable oscillations in the SPH motion.When this additional viscosity is included, the viscosityterm Π in Eqs.(29) and (30) should be replaced by Π T ot and the coefficients A and B in Eqs.(31) and (32) aremodified as A → γ " ε + p + Π T ot − ∂ ( ε + p ) ∂s (cid:18) s + Π T ot T (cid:19) − ζτ R − ζ ( h ) τ ( h ) R , (39) B = A γ σ ∗ dσ ∗ dt + Π τ R + Π ( h ) τ ( h ) R . (40)Consequently, the dispersion relation for the linearizedsound wave, Eq.(11) is modified as ω − αk = i ζ ( ǫ + p ) ωk iτ R ω + i ζ ( h ) ( ǫ + p ) ωk iτ ( h ) R ω . (41)In Figs. 8 and 9, we show the real and imaginaryparts, respectively, of the frequency ω as function of k ,calculated from Eq.(41) for a = 0 . , b = 6 and T =200 M eV . This time, there are 4 solutions, but only 2of them are propagating modes and the other 2 are non-propagating modes. All them have positive imaginarypart so that our equation is stable to a linear perturbationaround a hydrostatic state.Under the normal condition, τ ( h ) R /τ R ≪ , ζ ( h ) ≪ ζ, the above equation gives the same dispersion relation asthe previous case, as far as τ ( h ) R ω < , since the secondterm of the right-hand side of Eq.(41) can be neglectedcompared to the first term. On the other hand, for thevery large k limit, the asymptotic form of the dispersionrelation becomes ω → ± r α + 1 b + 1 b ( h ) k + i
12 1 / ( bτ R ) + 1 / ( b ( h ) τ ( h ) R ) α + 1 /b + 1 /b ( h ) . The most effective choice of the additional viscosity isobtained for the smallest possible value of τ ( h ) R , preservingcausality. This determines the value of b ( h ) as b ( h ) = 11 − α − /b . The choice Eq.(38) corresponds to the case b = 6 and α = 1 / . The difference of the two regimes determinedby each viscosity can be seen in these figures.The additional viscosity introduced here can be consis-tent with the stability and causality of the theory. Withthis procedure, the rapid oscillation modes associatedwith the degrees of freedom beyond the applicability ofthe coarse-grained theory can be naturally eliminated.
FIG. 8: The real part of the frequency ω as a function of k .There are two propagating modes A, B (dashed lines) andtwo non-propagating modes C and D which are degenerated(solid).FIG. 9: The imaginary part of the frequency ω as a functionof k . The two propagating modes A and B are degenrated(dashed line). The non-propagating modes C and D havedifferent k dependence. A. “Double Shock” Phenomena
We have already shown that, without the additionalviscosity, the SPH time evolution develops rapidly oscil-lating modes with large amplitude whose wavelength areof the order of h. In Fig. 10, we show the result of thesame time evolution with the additional viscosity. We seethat this new viscosity successfully smears out the quickoscillating modes.1
FIG. 10: The same example as Fig. 7, calculated with theadditional viscosity. Quick osculations at the shock region areremoved.FIG. 11: The velocity profile corresponding to Fig. 10.
A very interesting feature of this example, in Fig. 10, isthe appearance of non-trivial structure of discontinuitieswhen compared to the usual shock discontinuity seen inFig. 4. By changing the viscosity parameter a from thesmall value to the large ones, we found that the discon-tinuity associated with the usual shock corresponds tothe discontinuity A shown in this figure. The new dis-continuity B starts to appear at some critical value ofparameter a. This can be interpreted as the transporta-tion of the discontinuity in the initial distribution by acausal diffusion mechanism, as is known in the case ofmatter transportation by a telegraph equation (See Ap-pendix ). In fact, the velocity profile for this apparent“double shock” phenomena does not have any disconti-nuity at the location of B as is shown in Fig. 11.
FIG. 12: The shock formation of the proper entropy densityof the ideal fluid, starting from the homogeneous initial con-dition (dotted).
B. Shock wave Formation by quickly moving fluidcomponent
Another typical example where shock discontinuitiesappear is that when a component of fluid is acceleratedby an external force and achieves a velocity greater thanthat of the sound in the medium. Such phenomena areparticularly interesting in heavy ion physics. When highenergy partonic jets punch out the thermally equilibratedQGP, they may transfer a large amount of momentumand energy, dragging the piece of the fluid with high ve-locity. Such a scenario is often discussed with the ob-served angular correlation of produced hadrons, connect-ing to the formation of Mach cones[3, 9]. Of course, theMach cones do not exist in 1+1 dimensional systems,but it is interesting to study the dynamical formation ofa shock discontinuity in our theory. If we calculate suchsituation without the additional viscosity, the formationof shock leads to quick unphysical oscillations as shownin Fig. 12. Here, a small part of the fluid has an finite ini-tial velocity (Lorentz factor γ ≈
2) with constant entropydensity in the local rest frame. However as is discussedin the previous section, the inclusion of the additionalviscosity can eliminate these unphysical oscillations asshown in Fig. 13, confirming the efficiency of our addi-tional viscosity. Differently from the case of discontinuityinitial condition shown in Fig. 4, we found that the shockfront created by a rapidly moving fluid element cannotbe smoothed out even in the presence of normal viscos-ity. The effect of piling up the matter at the shock frontgenerates a steep density variation whose wavelength be-comes eventually smaller than the coarse-graining scale h. Therefore, the introduction of the additional viscosityis essential to calculate the dynamical shock formationprocesses.In Figs. 14 and 15, we show the entropy density pro-2
FIG. 13: The same as Fig. 12, calculated with the additionalviscosity.FIG. 14: The same as Fig. 13 with a = 0 . a = 1. FIG. 16: The corresponding velocity profile of the viscousfluid of a = 1. file of the same situation for different values of viscouscoefficient a with additional viscosity. It is interesting tonote that a rarefaction wave is formed and propagatesbackwards with a smaller velocity, as is shown in thesefigures, although the fluid velocity is positive everywhere,as is shown in Fig. 16. Note that, another density dis-continuity is formed in the rarefaction wave part whichbehaves exactly in the same manner as the case of shockfront created by the initial density discontinuity. VI. DISCUSSION AND CONCLUDINGREMARKS
In this paper, we studied the causal dissipative hydro-dynamics in 1+1 dimensional systems. To clarify theeffects of viscosity and its relation to the coarse-grainingaspect of the theory, we apply the SPH formulation torepresent the macroscopic variables. Here, we considerthe SPH scale parameter h as the scale of coarse-graining.We have shown that once the viscosity is introducedappropriately, then the expansion of the fluid into thevacuum should form a steady wave and there the univer-sal relation p + Π = 0 should be satisfied. We also stud-ied the various situations where the shock discontinuitiesemerge. We argue that for the consistency of the causaldissipative hydrodynamics, we have to introduce an addi-tional viscosity associated with the coarse-graining scaleof the theory. This is because, the normal viscosity ob-tained from microscopic theories, such as kinetic equa-tions or the Green-Kubo-Nakano formula, is associatedwith the bulk properties of the matter and is calculated inthe vanishing coarse-graining scale h → . To be consis-tent with the coarse-graining scheme of the theory, thereshould be additional dissipation mechanisms that disap-pears in the limit of vanishing h. We proposed a schemewhich contemplates such an additional viscosity, keeping3causality and stability of the theory.Most of the results in this study will be more relevantwhen the role of viscosity becomes effective. Such situ-ations will be realized in the coming LHC experimentswhere we may expect that a large inhomogeneity in thevelocity profile can be created. The application of thepresent theory for the LHC energy is now in progress.One might think that the problems discussed heremight be out of the range of applicability of hydrody-namics in the usual argument based on the smallness ofthe Knudsen number. However, as we already pointedout, it is well-known that hydrodynamics is still appliedeven for cases where the deviation from equilibrium isapparently large. In fact, the microscopic derivationof hydrodynamics is still an open question. Interestingworks are still under progress to justify the applicabil-ity of hydrodynamics based on the asymptotic theory,the fluctuation theory and so on. In our opinion, whilewe do not establish a precise theoretical criteria for itsuse, hydrodynamics should be explored independent ofsuch formal limitations, since it contains very importantingredients to describe phenomenologically the collectiveaspects of matter. As a matter of fact, if we had really torestrict ourselves to the region where the hydrodynamicsapproach is clearly available, we would never discuss thephase transition dynamics in this scheme [28], so that allthe ideal fluid model for heavy-ion collisions would be-come meaningless. In the same sense, we consider thatit is significant to study the causal dissipative hydrody-namics comprehensively.In this paper, we discussed the problem of causalityand stability of our theory, and concluded that the theoryis stable for the small linear perturbations. To be precise,we did not prove that the theory is stable for nonlinearperturbations, although the numerical solutions exam-ined here did not show any instabilities. However, if thetheory is shown to be unstable in some regime, then wehave to construct another stable theory. Analysis on thisline is under investigation and will be reported in anotherpaper. The effect of the equation of state containing phase transitions also will be discussed in future.This work has been supported by CNPq, FAPERJ,CAPES and PRONEX.
APPENDIX A: TELEGRAPH EQUATION
We discuss the behavior of the causal diffusion equa-tion, which is described by τ ∂ ∂t n + ∂∂t n − D ∂ ∂x n = 0 . As is well-known, the maximum propagation speed ofsignals described by the equation is v = p D/τ . As isdiscussed in [29], the solution of the causal diffusion equa-tion is given by n ( x, t ) = 12 e − t τ [ n ( x + vt ) + n ( x − vt )]+ e − t τ Z x + vtx − vt dx (r Dτ I " r Dτ p v t − ( x − x ) + 12 r τD ∂∂t I " r Dτ p v t − ( x − x ) n ( x )+ 12 r τD e − t τ Z x + vtx − vt dx I " r Dτ p v t − ( x − x ) × ∂∂t n ( x ) , (A1)where n ( x ) is the initial distribution and I ( x ) is themodified Bessel function. The first term of the solutionindicates that the initial distribution n ( x ) is separatedinto two fragments, which travel to the opposite direc-tions with the velocity v . Thus, in the shorter time scalethan the relaxation time τ , the memory of the initialdistribution profile survives near the boundary. [1] See for example, P. Huovinen and P.V. Ruuskanen, Ann.Rev. Nucl. Part. Sci. , 163 (2006); Jean-Yves Olli-trault, Euro. J. Phys. , 275 (2008) and referencestherein.[2] D. Teany, Phys. Rev. C68 , 034913 (2003); P. Van, T. S.Biro, Eur. Phys. J. ST , 201 (2008).[3] J. Noronha, G. Torrieri and M. Gyulassy,arXiv:0712.1053; B. M¨uller and J. Ruppert,arXiv:0802.2254.[4] A. Muronga, Phys. Rev. Lett. 88, 062302 (2002) [Erra-tum ibid. 89, 159901 (2002)]; Phys. Rev. C 76, 014909(2007) ; P. Romatschke and U. Romatschke, Phys. Rev.Lett. 99, 172301 (2007) ; H. Song and U. W. Heinz,arXiv:0712.3715; A. K. Chaudhuri, arXiv:0801.3180; K.Dusling and D. Teaney, arXiv:0710.5932; R. S. Bhaleraoand S. Gupta, Phys. Rev.
C77 , 014902 ; A. Dumitru E. Mol´ar and Y. Nara, Phys. Rev.
C76 , 024910 (2007); S.Pratt, Phys. Rev.
C77 , 024910 (2008); and referencestherein.[5] T. Koide, G.S. Denicol, Philipe Mota and T. Kodama,Phys. Rev.
C75 , 034909 (2007).[6] P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev.Lett. 94, 111601 (2005).[7] F. Karsch, D. Kharzeev and K. Tuchin, arXiv:0711.0914.[8] See for example P. Steinberg, nucl-ex/0702019 and refer-ences therein.[9] H. St¨ocker, Nucl. Phys. A 750 (2005) 121; J. Casalderrey-Solana, E. V. Shuryak and D. Teaney, J. Phys. Conf. Ser.27, 22 (2005).[10] See for example D. H. Rischke, S. Bernard and J. A.Maruhn, Nucl. Phys. A , 346 (1995); D. H. Rischke,Y. P¨urs¨un and J. A. Maruhn, Nucl. Phys. A , 383 (1995); For a classical paper, J. VonNeumann and R. D.Richtmyer, J. Ap. Phys. , 232 (1950).[11] D. Jou and D. Pav´on, Phys. Rev. A44 , 6496 (1991).[12] T. Ruggeri, Phys. Rev.
E47 , 4135 (1993).[13] I. M¨uller, Z. Phys. 198, 329 (1967); as a review paper,see I. M¨uller, Living Rev. Relativity , 1 (1999).[14] W. Israel and J. M. Stewart, Ann. Phys. (N.Y.) , 341(1979).[15] As a review paper, see D. Jou, J. Casas-V´azquez and G.Lebon, Rep. Prog. Phys. , 1105 (1988); ibid , 1035(1999).[16] B. Carter, Proc. R. Soc. London, Ser A, , 45 (1991);as a review paper, see N. Andersson and G. L. Comer,Living Rev. Relativity , 1 (2007).[17] M. Grmela and H. C. ¨Ottinger, Phys. Rev. E56 , 6620(1997).[18] L. D. Landau and E. M. Lifshitz,
Fluid Mechanics , (Perg-amon; Addison-Wesley, London, U.K.; Reading, U.S.A.,1959).[19] W. A. Hiscjock and L. Lindblom,Ann. Phys. (N.Y.), , 466 (1983); Phys. Rev.
D31 , 725 (1985); Phys. Rev.
D35 , 3723 (1987); Phys. Lett.
A131 , 509 (1988).[20] H. Kouno, M. Maruyama, F. Takagi and K. Saito, Phys.Rev.
D41 , 2903 (1990).[21] G. Torrieri and I Mishustin, arXiv:0805.0442.[22] T. Koide, Phys.Rev.
E72 , 026135 (2005).[23] T. Koide, Phys. Rev.
E75 , 060103(R) (2007).[24] T. Koide and T. Kodama in preparation.[25] L.B. Lucy ApJ. , 1013 (1977), J.J. Monaghan, Annu.Rev. Astron. Astrophys. , 543 (1992).[26] C. E. Aguiar, T. Kodama, T. Osada and Y. Hama, J.Phys. G27 , 75 (2001).[27] Y. Hama, T. Kodama and O. Socolowski, Braz. J. Phys. , 24 (2005).[28] See, for example, the discussion in the introduction in T.Koide and M. Maruyama, Nucl. Phys. A742 , 95 (2004).[29] P. M. Morse and H. Feshbach,