The Existence and Stability of Relativistic Shock Waves: General Criteria and Numerical Simulations for a Non-Convex Equation of State

A small viscosity approach to discontinuous flows is discussed in relativis-tic hydrodynamics with a general (possibly, non-convex) equation of state that typically occurs in the domains of phase transitions. Different forms of criteria for the existence and stability of relativistic shock waves, such as evolutionarity conditions, entropy criterion and corrugation stability conditions are compared with the requirement of the existence of shock viscous profile. The latter is shown to be most restrictive in case of a single-valued shock adiabat expressed as a function of pressure. One-dimensional numerical simulations with artificial viscosity for a simple piecewise-linear equation of state are carried out to illustrate the criteria in the case of pla-nar and spherical shock waves. The effect of a phase transition domain on the shock amplitude in the process of a hydrodynamical spherical collapse is demonstrated.


Introduction
A number of applications of relativistic hydrodynamics deal with shock wave propagation in a matter at extremely high densities.This problems concern, e.g., the relativistic stellar collapse, the theory of early Universe, hydrodynamical models of elementary particle production in nucleus-nucleus collisions.Such a consideration usually requires taking into account various phase transitions that may occur in a super-dense matter under appropriate conditions.The domain of phase transitions typically has essential peculiarities from the viewpoint of hydrodynamics because of the possible non-convexity of the equation of state.Namely, this requires special attention to the existence and stability properties of discontinuous solutions 1,2].The main problem is that standard equations for an ideal uid ow with the general equation of state do not restrict properly the discontinuous solutions and the question is how to distinguish the physical solutions obtained either by analytical or numerical algorithms from those that cannot be realized in nature.This problem is well known in classical hydrodynamics 1,2] and needs to be investigated in a relativistic case.
The equations for a uid ow of relativistic hydrodynamics are de ned by conservation laws for energy-momentum and baryonic charge (or another conserved charge) combined with the equation of state.However, in the presence of discontinuities (shock waves) these equations must be complemented by additional restrictions in order to de ne correctly a generalized solution and to provide its uniqueness.Otherwise, there may be di erent solutions satisfying the same discontinuous initial data.This situation is typical of the quasilinear equations theory 1,2].The criteria for the existence of a viscous pro le (EVP), entropy criterion (EnC) and evolutionarity criterion (EvC) may be considered as the most famous examples of the above restrictions 1{3].EVP selects the discontinuous solutions that can be obtained from continuous solutions of hydrodynamical equations with a non-zero viscosity in the small-viscosity limit.The physical background of this approach involves a discussion of validity of a viscous hydrodynamics approximation to describe the shock wave structure.This approximation is not always physically admissible and this diminishes the role of an EVP criterion.Nevertheless, a small viscosity model is physically well understandable and may be useful in phenomenological schemes of relativistic hydrodynamics (see, e.g., 4,5]) dealing with a superdense matter, where investigation of the shock structure at a microlevel is rather complicated.EnC prohibits unphysical shocks with the decreasing entropy.In many cases it is the only criterion that may be taken into account, but in the case of a non-convex equation of state it is not su cient 1,2].EvC may be considered as a general requirement that provides a correct construction of the solution by the determination of hydrodynamical parameters on the shock 3].The solutions that do not satisfy these requirements are considered to be unphysical and must be rejected on a formal level 1{3].An important practical question concerns the validity of numerical algorithms and their correspondence to the above criteria.
Another set of criteria concerns the stability of shock waves.According to the general notion of stability these criteria do not prohibit (mathematically) the existence of certain solutions.However, they show that some solutions are destroyed with time by exponentially growing small perturbations and, therefore, they cannot be realized in nature.The corrugation stability criterion 6,7] examines such perturbations of hydrodynamical quantities in a continuous ow and three-dimensional deformations of the discontinuity front itself.
In classical hydrodynamics much attention has been paid to compare di erent criteria for the existence and stability of shock waves.The aim of this paper is to outline some results concerning such a comparison in relativistic hydrodynamics.
The plan of the paper is as follows.In section 2 we review the conditions for EVP of relativisic shock waves with some generalization that allows one to compare di erent types of viscosity.In section 3 the EVP is compared with the corrugation stability conditions.Section 4 deals with numerical simulations of the shocks in the case of plane and spherical symmetry for a simple piecewise-linear equation of state that models a phase transition.We compare the numerical and analytical results in a plane case and study some qualitative e ects due to the domain of phase transition in a spherical collapse.Section 5 summarises the results of the paper.

Existence of a viscous profile of a shock wave
A detailed investigation of the EVP criterion for relativistic shocks in an ideal uid with a general equation of state has been carried out in 4,5].It is important to note that the results of 4,5] do not use supposition about the convexity of Poisson adiabats, that is one of the Bethe{Weyl conditions (see, e.g., 1]) for a normal medium.In this section we reconsider some of these results in order to compare them with the stability criteria and with numerical algorithms involving the Neumann{Richtmyer arti cial viscosity.
The equations of motion of an ideal relativistic uid can be written as conservation equations of the energy-momentum tensor 3] T = (" + p) u u ?pg ; (1) where u is 4-velocity of the uid, fg g = diag(1; ?1; ?1; ?1), p is pressure and " is the proper energy density.The pressure is related to " and to the baryon number density n (or density of some other conserved charge) by a su ciently smooth equation of state p = p ("; V ), V = 1=n being the speci c volume.
The conservation equations for T must be complemented by conservation equations for the baryonic charge @ (nu ) = 0: ( In a discontinuous ow the hydrodynamic quantities on both sides of the shock are related through the integral form of conservation laws.However, relations on the discontinuities must be complemented by certain restrictions that are necessary to provide the uniqueness of discontinuous solutions 1,2].These restrictions may be obtained by analyzing the structure of a shock transition in the presence of small dissipative e ects, such as viscosity.In the presence of viscosity the conservation equation for T must be replaced by where = R visc , visc = (u ; + u ; ?u u u ; ?u u u ; ) + ? 2 3 @u @x (g ?u u ) , (4) is the Landau{Lifshits relativistic viscosity tensor 3]; and the multiplyer R, that makes equation (3) somewhat di erent from that used in 4,5], may be any positive function of ow parameters to be speci ed later.Under the given equation of state and appropriate initial conditions, equations ( 2) and (3) de ne solutions for the viscous uid dynamics.To obtain the discontinuous ideal uid motion we suppose that some limit of such a dissipative solution exists for !0, !0. This limit is considered as a generalized solution of the conservation equations for an ideal uid.This is one of a number of possible methods to de ne discontinuous solutions; there may be, e.g., other dissipative e ects besides viscosity; one may also use di erent de nitions of a generalized solution of hydrodynamic equations of motion 3].However, the approach based on the viscosity tensor 3] appears to be the simplest and physically well understandable.It is su cient to obtain a continuous pro le of the shock, the resulting restrictions being more stringent, e.g., than the evolutionarity conditions.
Thus, we expect that in the limit of small viscosity the solutions of hydrodynamical equations yield discontinuous ows describing shock waves.In the case of a nonzero viscosity the stationary shock wave propagating in a space-like direction l is locally represented in a proper frame by a stationary viscous ow depending upon the only variable x = x l .We choose the coordinates in such a way that fl g = f0; 1; 0; 0g, x = x 1 .Further, without loss of generality, one may also put u 2 u 3 0 due to the choice of the reference frame.
The equations ( 1), ( 2) yield T 1 + 1 = T 1 0 = const, (5) nu 1 = j = const; (6) where j = n 0 u 1 0 and T 1 0 are asymptotic quantities either for x ! 1 or for x !?1.Here we assume that for x !?1, all the quantities tend to some constant values marked by \0" (the state ahead of the shock) and for x ! 1, the quantities tend to constant values marked by \1" (behind the shock).Because !0 for x ! 1, the solution of equations ( 5), (6) (if it exists) satis es the conservation relations that connect hydrodynamical quantities on the both sides of the shock.
However, there arises the question whether there is a continuous one-dimensional ow represented by the solution of ( 5), ( 6) that satis es the both above asymptotics for x ! 1 and x !?1.The existence conditions for a continuous solution of this boundary value problem are then interpreted as admissibility conditions for shock transition u (0) ; V 0 ; p 0 !u (1) ; V 1 ; p 1 with the corresponding states satisfying the relations on the shock wave.
Equations ( 5), ( 6) can be reduced, analogously to 4,5], to the rst order ordinary di erential equation R + 4  3 e " (V ) = T 1 (0) u =u 1 : (9) The right-hand side of equation ( 7) equals zero for V = V 0 and V = V 1 .Now one must analyse di erent choices of R. a) R 1.This case has been studied in 4,5].Let, e.g., V 1 > V 0 (the reverse case yields the same results), and the r.h.s. of ( 7) is smooth and positive for V 2 (V 0 ; V 1 ).Then, there is a smooth solution V (x) of ( 7) that tends to V 0 for x !?1 and to V 1 for x ! 1.If the right-hand side of (7) changes its sign at some point between V 0 and V 1 , then there is no continuous solution of ( 7) in (?1; 1) connecting V 0 and V 1 .As a result, we obtain the following criterion of admissibility of a stationary shock transition: ) for all V between V 0 and V 1 .The cases when the left-hand side of (10) has zeros within (V 0 ; V 1 ) but does not have negative values (e.g. in the case of Joquet points) may be also included, but this requires some additional consideration.
b) The hydrodynamical numerical calculations often deal with an arti cial viscosity of the Neumann{Richtmyer type 1,11].We shall take this into account by putting R = C dV dx in (4), C being a positive constant.Let again, for de niteness, V 1 > V 0 , and one has a strict inequality in (10).In this case, expanding equation (7) in the neighborhood of V 1 and V 0 , one obtains dV=dx / (V 1 ?V ) 1=2 and an analogous relation for V 0 .Using this it is easy to see that we have a nite interval x 0 ; x 1 ] for the smooth solution V (x), such that V (x 0 ) = V 0 , V (x 1 ) = V 1 ( nite width of the shock front, cf.1]).Outside this interval the solution may be continued by putting V (x) V 0 for x < x 0 , V (x) V 1 for x > x > x 1 .
On the other hand, if there is point V 0 2 (V 0 ; V 1 ) where the right-hand side of (7) changes its sign, then the solutions of (7) must be increasing for V < V 0 and decreasing for V > V 0 .Therefore, there is no continuous solution connecting states \0" and \1".Therefore, we again come to criterion (10).
It is shown in 4,5] that inequality (10) is equivalent to another one (" 1 ?" 0 ) (e p (") ?p H (")) > 0 (11) expressed in terms of the Hugoniot{Taub shock adiabat p H (") in plane of thermodynamical variables p ? "; the pattern function e p (") = A ? C= (A + ") represents curve (8) passing through the initial and the nal state in the p ? " plane of thermodynamical parameters (this condition de nes constants A and C).Inequality (11) must hold for all the values of energy density " between states \0" and \1", ahead of and behind the shock, respectively.This form of the EVP criterion works in the case of the single-valued Hugoniot{Taub adiabat p H ("). If this is not the case, one may use another form 4,5] of the criterion that uses an equation of state only.An equivalent assertion in terms of shock adiabats may be formulated as follows.
The compression (rarefaction) shock transition has a viscous pro le on the condition that pattern curve e p (") lies above (below) the Hugoniot{Taub adiabat in the vicinity of the initial point and does not intersect it between the initial and the nal states.
It has been shown 4,5] that EnC follows from EVP, but the reverse statement is not valid.The evolutionarity conditions for relativistic shocks yield the inequalities for uid velocities and sound speeds v 0 ; v 1 ; c 0 ; c 1 ahead of and behind the shock v 0 > c 0 ; v 1 6 c 1 , (12) that have the same form as for nonrelativistic ones 3].Inequalities (12) can be also be obtained as a consequence of EVP 4,14].Evidently, EVP imposes more powerful restrictions upon the parameters of the shock because (10) or (11) concerns all the states on the Hugoniot{Taub adiabat, not only in the neighbourhood of \0" and \1".However, in the case of a normal equation of state (that allows only compression shocks) these criteria are equivalent.

Shock wave instabilities and viscous profiles
Consider now the corrugation stability conditions 6,7] dealing with perturbations of the solutions of hydrodynamical equations.The shock transition \0"!\1" is corrugationally unstable 6,7] if at the nal state p = p 1 , = 1 either m 2 @ @p H < ?1 where M = v 1 =c 1 is the Mach number behind the shock; v is the uid velocity with respect to the shock front, m= v= p 1 ?v 2 is a conserved ux, index H means that the derivatives are taken along the Hugoniot{Taub adiabat.
On account of the relation m 2 = (p ?p 0 ) = ( 0 ? ) we see that by virtue of ( 14) the relative disposition of the shock adiabat and the Rayleigh line in the neighbourhood of the nal point p = p 1 contradicts (13).Then condition ( 14) is not satis ed for shocks satisfying EVP.
To analyse the second inequality (15) we pass on to variables u; p, where u is the uid velocity ahead of the shock with respect to the uid behind the shock, u 2 = (p ?p 0 ) (" ?" 0 ) (" 0 + p) (" + p 0 ) . (16) Let p H (u) represent a shock adiabat in the p?u plane.We shall now analyze inequality (15) in view of the requirement that p H (u) is a single-valued function.This is not a physical requirement; and, indeed, there are examples of non-single-valued adiabats 2].However, the case of single valued adiabats p H (u) covers a signicant number of physical applications; and in the opposite case the investigation of discontinuous solutions requires additional information of a non-hydrodynamical nature 2], 9].
Due to this equation we have the derivative (@"=@p) H that is used to obtain (@u=@p) H from (16); then, substitution into (15) yields @u @p H (1 ?v 0 v 1 ) 2 < 0 for p = p 1 . (17 On the other hand, direct calculation of this derivative from (17) for the initial point yields (@u=@p) H = 1= c 0 (" 0 + p 0 )] > 0. So, this derivative changes its sign between the initial and the nal states, and, therefore, dependence u H (p) is not monotonous.This conclusion is the same as in classical hydrodynamics.In other words, shock adiabat p H (u) is not single-valued.
Therefore, either the corrugationally unstable shocks do not satisfy EVP, or they are non-single-valued in the above sense.
So, we infer that shocks for the single-valued Hugoniot-Taub adiabat segments satisfying criterion (13) are evolutionary, corrugationally stable and satisfy the entropy criterion.Discontinuities that correspond to the segments of the shock adiabat where ( 13) is not valid break up, although their existence is not prohibited by conservation laws.
If the Hugoniot{Taub adiabat p H (u) is not single-valued, this statement generally does not hold.Such shocks can break up in a non-unique way 2].Even stable shocks may break up for non-single-valued Hugoniot{Taub adiabats.Consideration of shocks in this case needs additional information of a non-hydrodynamical nature 2,9].

Numerical simulations of shock waves
We have shown in section 2 that introduction of either the Landau-Liftshits viscosity or an analytical counterpart of the Neumann-Richtmyer viscosity leads to similar criteria for the shock existence in a one-dimensional ow.However, strictly speaking, when we pass on to a nite-di erence analogues of hydrodynamical equations, such a comparison needs special mathematical investigation of a concrete numerical algorithm.In this connection we outline the results of numerical simulations of discontinuous ows in order to illustrate the correspondence between the above EVP criterion and an algorithm involving arti cial viscosity.In this section we also apply a numerical algorithm to study the formation of a shock wave in a hydrodynamical spherical collapse.
To perform the numerical calculations we have worked out a computer program that realizes an absolutely conservative numerical algorithm for relativistic hydrodynamical equations using small viscosity of the Neumann{Richtmyer type.The principle of absolute conservation means that not only conservation laws are approximated numerically but so are all their di erential consequences.Such a property permits us to work with rather rough grids (which do not require powerful hardware) in order to obtain qualitatively good results by means of personal computers.For a detailed description of the method for nonrelativistic equations see 11]; we have modi ed this approach to the relativistic one-dimensional ow in the case of planar and spherical symmetry.
A. In a planar case the problem we deal with is essentially a break-up problem of arbitrary discontinuity which is analogous to the well-known classical problem 1].An algorithm of an analytical solution for a relativistic version of this problem can be found in 13,14]; a detailed analytical investigation of discontinuous initial data break-ups for the above equation of state is presented in 14].The hydrodynamical ows of the breakup problem depend upon the only similarity variable x=t (t is time, x { a spatial variable).The analytical solutions for discontinuous initial data have been compared with the results of computer simulations.The results shown in gures 1 and 2 demonstrate that the initial con gurations not satisfying (13) on discontinuities are disrupted into smaller shocks and simple waves.The results also show the existence of admissible (satisfying ( 13)) rarefaction shocks compatible with analytical calculations.B. In the case of spherical symmetry we cannot construct an analytic solution for the equation of state (18) with the same degree of completeness as in the planar case.To check up the program in this case we use a speci c self-similar solution describing evaporation of a liquid drop through a converging rarefaction shock wave.The solution depends upon the only variable = r=t (r being a radial variable).The trajectory of the converging shock is = ?0 .
Outside the drop (j j > 0 ) the solution is described by an ordinary di erential equation (see, e.g., 15]) that is solved numerically.The initial values for v and " in the hadron phase are obtained from " 0 and 0 through conservation relations that relate the hydrodynamical quantities on the both sides of the shock.Then this result is compared with the results of numerical simulations using arti cial viscosity.These are also in good agreement (see gure 3).C. The following question concerns the application of the above algorithm to the study of a role of the shelf part (modelling the phase transition domain) in the equation of state (18) in a hydrodynamical spherical collapse.As initial conditions for the converging rarefaction shock wave we take a spherically symmetric con guration with energy density " (r; 0) = 1= (2 + r 2 ) and with homological conditions on velocity v (r; 0) = ?r.We consider a class of equations of state (18) with di erent sizes of the \shelf" " = " 2 ?" 1 and our aim is to study the qualitative dependence of the ow upon it.The results of the calculations are depicted in gure 4.These gures show the sequential time slices of energy density " (r; t) as a function of the Eulerian coordinate.and c) correspond to the values of the \shelf" " equal to 0,5,17, respectively.The wave-like structures and oscillations in the numerical solution are the consequence of the grid roughness.The \tails" of the solutions at the right-hand side of the curves correspond to a rarefaction wave of the outer layers in view of a nonequillibrium initial con guration that is cut o at " = 0:05: As it can be seen from the gures, the shelf part of the EOS causes an emergence of a shock wave moving from the centre of the con guration to the outer layers.This takes place at the moment when the energy density reaches " 2 .This is not surprising because the sti ening of the equation of state at this point acts as an elastic piston.Such a shock emerging in a star core in the collapse process will press on the outer layers of the star as a piston and may contribute to the envelope overthrowing in supernovae explosions.Dependence of the intensity of the shock wave upon the size of a phase transition, i.e. on the size of the \shelf", is depicted in gure 5.As we see, the larger ", the stronger is the shock, i.e. the stronger is the e ect of the piston on the outer envelope.

Conclusions
We have compared di erent known criteria for the existence and stability of relativistic shock waves: the existence of a viscous pro le, entropy criterion, evolutionarity criterion and corrugation stability.These criteria have di erent senses and levels of physical justi cation and di erently restrict discontinuous ows.The EnC and EvC follow from very general requirements, however, in the case of nonconvex equations of state they do not properly restrict discontinuous solutions of hydrodynamical equations.The solutions may be speci ed due to EVP requirements; in particular, this criterion does work in the relativistic problem of decay of arbitrary discontinuity 14].However, justi cation of this criterion is less general from the physical point of view.This criterion was obtained in 4,5] by using the phenomenological relativistic viscosity tensor 3] (in a small viscosity limit) for a one-dimensional picture of the shock wave, so one may enquire about other models of the relativistic shock structure yielding other ways of studying discontinuous solutions.Here we do not discuss possible physical situations related to this model, but concentrate on formal results.
Though the relativistic EVP criterion of 4,5] is essentially one-dimensional and the corrugation stability concerns three-dimensional perturbations of the uid ow, in the case of single-valued shock adiabats (in the p ? u plane) the requirement of EVP is shown to be most restrictive: the shocks satisfying this criterion also satisfy EvC, EnC and are corrugationally stable.This is analogous to the non-relativistic case 2].Though the requirement of a single-valued shock adiabat is not a physical one 2], this case covers a signi cant number of physical applications.If the shock adiabat is not single-valued, then there may exist corrugationally unstable shocks that satisfy EVP, EvC and EnC 2,9] and their consideration requires additional information of a non-hydrodynamical nature.Note that the results are valid for a wide class of equations of state including those having anomalous domains that are typical in the neighbourhoods of phase transitions.
It is also important to note that the results of the computer simulation based on widely used arti cial viscosity schemes agree with the EVP criterion.This has been illustrated by numerical simulations using an absolutely conservative scheme with arti cial viscosity for special one-dimensional relativistic ows with a simple piecewise-linear equation of state.The results of these simulations show an essential in uence of the model phase transition domain on the amplitude of the shock arising in the spherical hydrodynamical collapse.This may be of interest for more detailed considerations of a relativistic stellar collapse in the supernova ares theory.

Figure 1 .
Figure 1.Numerical simulations using the Neumann{Richtmyer viscosity (crosses) versus analytical solutions (solid line).The picture describes an energy density prole that evolves from an \unstable" plane compression shock.

Figure 2 .
Figure 2. Energy density pro le evolving from the desruption of an unstable plane rarefaction shock.

Figure 3 .
Figure 3. Spherically symmetrical discontinuous ow in the case of a converging rarefaction shock.The solid line shows the solution obtained by numerical integration of ordinary di erential equations corresponding to the self-similarity case; the crosses show the results of an arti cial viscosity algorithm.

Figure 4 .
Figure 4. Energy density pro les in the case of a spherical collapse for di erent time moments.The gures a), b) and c) correspond to the values of the \shelf"" equal to 0,5,17, respectively.The wave-like structures and oscillations in the numerical solution are the consequence of the grid roughness.The \tails" of the solutions at the right-hand side of the curves correspond to a rarefaction wave of the outer layers in view of a nonequillibrium initial con guration that is cut o at " = 0:05:

Figure 5 .
Figure 5. Numerical simulations using the Neumann{Richtmyer viscosity (crosses) versus analytical solutions (solid line).The picture describes an energy density prole that evolves from an \unstable" plane compression shock.