Causality of fluid dynamics for high-energy nuclear collisions

Dissipative relativistic fluid dynamics is not always causal and can favor superluminal signal propagation under certain circumstances. On the other hand, high-energy nuclear collisions have a microscopic description in terms of QCD and are expected to follow the causality principle of special relativity. We discuss under which conditions the fluid evolutions for a radial expansion are hyperbolic and how the properties of the solutions are encoded in the associated characteristic curves. The expansion dynamics is causal in relativistic sense if the characteristic velocities are smaller than the speed of light. We obtain a concrete inequality from this constraint and discuss how it can be violated for certain initial conditions. We argue that causality poses a bound to the applicability of relativistic fluid dynamics. }


Introduction
Relativistic fluid dynamics has been used with success as an effective theory of QCD in the regime of high energy density probed by relativistic heavy ion collisions [1][2][3][4][5]. It is seen as an expansion around local thermal equilibrium states in terms of derivatives. The lowest order in this derivative expansion is ideal fluid dynamics and describes in particular equilibrium states. Higher orders in a formal derivative expansion describe viscous and other transport corrections. The microscopic physics of QCD enters the formalism in terms of the thermodynamic equation of state and the transport properties such as viscosities, conductivities or relaxation times.
One of the most puzzling properties of relativistic fluid dynamics is that it is not always causal. More specific, it has been recognized early that relativistic versions of the Navier-Stokes equation, formulated either in the Landau frame or in the Eckart frame, allow for signal propagation with arbiltrarilly large velocity. This can be traced back to the fact that the non-relativistic Navier-Stokes equation is a parabolic equation (and not hyperbolic as for example the Klein-Gordon equation). Hyperbolic extensions of this setup have been proposed, most prominently by Israel and Stewart [7] (see also refs. [8,[23][24][25] for related discussions). The setup of Israel and Stewart goes beyond the Chapman-Enskog expansion and shear stress, bulk viscous pressure and diffusion current are dynamic fields with their own evolution equations. It was subsequently shown by Hiscock and Lindblom [26] that for appropriate values of the velocity of sound and transport properties, the theory of Israel and Stewart has equilibrium states around which small perturbations evolve in a relativistically causal way, indeed. They also investigated the stability of these linear perturbations and found that the conditions for causality also lead to linear stability of the homogeneous equilibrium states.
A generalization of Isreal-Stewart theory that is more complete in the sense that it takes all possible terms at second order in derivatives into account, has been put forward by Denicol, Niemi, Molnar and Rischke [9]; we will review their results in section 3. The relation between causality and linear stability for perturbations around equilibrium states was investigated in more detail in refs. [27,28]. It was confirmed that causality, in the sense of an asymptotic group velocity that is smaller than the speed of light, implies linear stability.
So far, investigations of causality were restricted to small perturbations of thermal equilibrium states. It is certainly necessary that such perturbations behave in a causal way, however, it is not sufficient to guarantee causality in more general situations. In fact, in the absence of a more general theoretical argument, causality needs to be established for every solution of the non-linear relativistic fluid equations of motion, case by case. This might be different only if the theory is organized in a scheme different from the conventional derivative expansion, for example as in ref. [29].
In the present work, we will consider a class of such solutions of relativistic fluid dynamics describing a fireball produced by high-energy nuclear collisions with longitudinal and transverse expansion. We discuss the corresponding solution of the equations proposed in ref. [9] and we discuss the issue of causality. Our main result will be concrete relations in the form of inequalities that tell whether perturbations propagating in the radial direction around the full, non-linear solutions of the fluid dynamic equations, behave in a causal way. This provides a non-linear form of the causality constraint and goes beyond the causality constraint for perturbations around equilibrium states. In particular, it turns out that the causality bounds are not satisfied for arbitrary initial conditions. We discuss what kind of initial conditions are save in this respect.
Throughout the manuscript we will work with natural units where c = = k B = 1 and we will take the signature of the Minkowski metric to be (−, +, +, +).
2 Hyperbolic partial differential equations and causality 2.1 Sets of hyperbolic partial differential equations In the present section we recall some mathematical knowledge about partial differential equations of the type relevant for dissipative relativistic fluid dynamics. This will lay the ground for a proper discussion of causality. We follow here mainly ref. [33].
Let us consider a set of first order partial differential equations of the type We will see that the fluid equations of motion for a longitudinally and radially expanding fluid can be expressed in this form under rather general conditions. The independent variables or coordinates are a time variable x 0 and a spatial (radial) variable x 1 . The dependent variables or fields are collected into the vector structure Φ j with index j = 1, . . . , n. The coefficient matrices A ij , B ij as well as the "source terms" C i are functions of the fields Φ and of the coordinates (x 0 , x 1 ) but do not depend on derivatives of Φ. Let us recall that the system of partial differential equations (2.1) would be called linear if A ij and B ij were independent of Φ and C i would depend at most linearly on Φ. If A ij and B ij were independent of Φ and the source terms C i some (non-linear) function of Φ, the system would be called semi-linear. In the more general case (relevant to us) where A ij and B ij are non-linear functions of Φ one calls (2.1) a quasilinear set of partial differential equations of first order. The equations (2.1) would be called homogeneous for C i (Φ) = 0.
We also assume that the matrix A ij in front of the time derivative is non-singular, det(A) = 0, and can be inverted. In the following we consider Cauchy's initial value problem, formulated as follows. We start with some initial configuration on some curve C. This curve represents a Cauchy surface in the case of 1+1 dimensions and will in practice be for example a line of constant time x 0 = const. More general, we assume that the curve C is specified by an equation ϕ(x 0 , x 1 ) = 0 with ∂ 0 ϕ(x 0 , x 1 ) = 0. We now ask whether the information given in terms of the field values Φ j (x) on the curve C is sufficient to determine the first derivatives of Φ j via the system of equations (2.1).
We first note that with Φ j (x) on C, we have also all the information to determine the internal derivatives 1 or derivatives along i. e. we can take D j to be also known along C. We can solve this equation for ∂ 1 Φ j and, using λ This is now a linear set of equations for the time derivatives ∂ 0 Φ j . Accordingly, a necessary and sufficient condition for the first derivatives to be uniquely determined by (2.1) along the curve C is given by Here, Q is known as the characteristic determinant. If Q = 0 along the curve C, the latter is called free. The fields Φ j along such curves can be continued into a "strip" where they solve (2.1). By going in this way from one Cauchy surface to the next, one can construct solutions of (2.1). If λ (m) is a real solution of the algebraic equation Q = 0, the solution of the differential equation This is obviously the case for α 0 = −∂1ϕ and α 1 = ∂0ϕ.
is called a characteristic curve. The real eigenvalue λ (m) is known as a characteristic velocity. Note that characteristic curves are not free in the sense defined above. For a given real solution λ (m) of Q = 0 one can find a corresponding (left) eigenvector w The system of n equations (2.1) is called hyperbolic if n linearly independent eigenvectors w (m) i with corresponding real eigenvectors λ (m) can be found. Note that the characteristic velocities λ (m) might be degenerate. If they were all different from each other, the system would be called totally hyperbolic.
It is instructive to use the left eigenvectors for an alternative formulation of the set of equations (2.1). In fact, they can be used to find so-called Riemann invariant or characteristic variables [34]. To this end we start from (2.1) in the form Contracting with the left eigenvectors leads to w (m) If one now introduces new variables J (m) such that Interestingly, we have now obtained a set of equations which is formulated in terms of derivatives along the characteristic curves labeled by the index m with corresponding velocities λ (m) . Without the inhomogeneous term ∼ C, the variables J (m) would actually be conserved along the characteristic curves dx 1 /dx 0 = λ (m) . This explains the name Riemann invariants or characteristic variables. The inhomogeneous terms lead to a modified behavior which typically results in an additional damping. Note that the causality structure of the system of partial differential equations is particularly transparent in the characteristic form (2.10). In each infinitesimal time step, information encoded in the spatial dependence of the variables J (m) is transported with the characteristic velocity λ (m) .
The initial conditions must be given on a Cauchy curve (or a Cauchy surface for more than one space dimension) which is free and therefore non-parallel to any of the characteristic curves. The Cauchy problem is then well posed and has (locally) a unique solution. Colloquially speaking, the initial conditions are specified on a Cauchy curve, while they are propagated along the characteristic curves.

Causality
The notion of causality can be made more precise for small (linear) perturbations around a given solution of the set of equations (2.1). For a given space-time point P = (x 0 , x 1 ) one can show that small changes in the initial conditions outside of a certain region in the past of P cannot change the solution at the point P . The region in the past of P where small changes in the initial conditions can influence the solution at P is called the domain of dependence Γ d . In a similar way, small changes at P can only affect a certain region in the future of P which is called domain of influence Γ i . The domain of dependence and domain of influence generalize the concepts of a past and future light cone familiar from electromagnetism to the present situation of non-linear (but quasi-linear) partial differential equations.
The domain of dependence Γ d and domain of influence Γ i are bounded by the characteristic curves with minimal and maximal characteristic velocity. This illustrates again the important role of the characteristic curves for the causal structure. We discuss this further in section 4.3 and in figure 4 we illustrate the regions Γ d and Γ i for a given point in the space-time history of a heavy ion collision. In appendix C we recall briefly the derivation of certain inequalities (so-called energy inequalities) which allow to prove the above statements about the domain of dependence and domain of influence for linear perturbations.
If the characteristic velocities are bounded in magnitude |λ (m) | ≤ v max , the domain of dependence and domain of influence of a certain point P are bound to lie in the interior of cones with velocity v max , similar to light cones. As long as there is some finite maximal velocity v max , one can in principle see this as a causal structure, similar to the causal structure of special (and general) relativity with v max replacing the velocity of light c.
In principle, if one would study a relativistic fluid in the absence of any other physical phenomena such as electromagnetism or gravity, it might be acceptable to have a causal structure with a maximal characteristic velocity v max larger than the speed of light. Note, however that also the fluid velocity could become as large as v max . However, in practice we are not interested in such an isolated situation but want to study a QCD fluid that interacts both with electromagnetic fields and (at least in principle) gravitational fields. Moreover, on a more microscopic level this fluid is governed by the laws of QCD as a quantum field theory, and the latter certainly imply an upper bound on the maximal velocity of signal propagation being given by the velocity of light, in agreement with the theory of relativity. It is for these physics reasons, that we demand for an effective (macroscopic) theory of a relativistic QCD fluid to follow the causality principle of special (and general) relativity according to which the maximal velocity of signal propagation, and therefore the maximal characteristic velocity, must be bounded from above by the velocity of light, |λ (m) | ≤ c = 1. (2.11) The last equality holds in our system of natural units.

Hyperbolic relativistic fluid equations to second order
The equations of motion for a relativistic fluid are the ones of energy-momentum conservation plus possibly additional conservation laws (such as for baryon number) and constitutive relations. We consider here a fluid where baryon number and other conserved quantities can be dropped and the only relevant conservation law is for energy and momentum, The energy-momentum tensor can be decomposed as where the fluid velocity u µ is defined in the Landau frame as the (unique) time-like eigenvector of T µ ν and the energy density is the corresponding eigenvalue.
The pressure p is related to by the same relation as in thermal equilibrium, the thermodynamic equation of state p = p( ) while π bulk is the bulk viscous pressure. The projector orthogonal to the fluid velocity is given by ∆ µ ν = u µ u ν + δ µ ν . Finally, the shear stress is symmetric π µν = π νµ , traceless π µ µ = 0, and transverse to the fluid velocity, u µ π µ ν = 0. The fluid velocity itself is normalized to u µ u µ = −1. The conservation law leads to the following evolution equations for the energy density and fluid velocity respectively, These equations must be supplemented by constitutive relations for the shear stress π µ ν and bulk viscous pressure π bulk , either in the form of constraint equations or as dynamical laws.
In what follows we shall work in the framework of a gradient expansion and include terms that are formally up to second order in gradients of temperature and fluid velocity. This setup corresponds to an expansion around local equilibrium configurations where the zeroth order corresponds to a locally equilibrated state characterized by temperature T (x) and fluid velocity u µ (x). One can also do the organization in terms of appropriately defined Knudsen and Reynolds numbers. The Knudsen number Kn is generically defined as a ratio between a microscopic scale, such as the mean free path, and a macroscopic scale of the fluid, such as the length over which the fluid velocity changes. The Reynolds number corresponds typically to the ratio of this macroscopic length scale to the dissipation scale, i. e. the length scale of perturbations that are efficiently damped by viscosity effects. In the present situation, it is convenient to define the inverse Reynolds number Re −1 as the ratio of dissipative fields such as π µν and corresponding equilibrium fields such as pressure or enthalpy density.
The most general equation for the shear stress π µ ν and bulk viscous pressure π bulk at second order in Knudsen number Kn and in inverse Reynolds number Re −1 have been obtained in ref. [9]. We are here interested in situations without a conserved particle number current. In this case, the evolution equation for the shear stress becomes We have used here the projector to the symmetric, transverse and trace-less part of a tensor, We also use the abbreviations Similarly for π bulk one finds the evolution equation The most important term in (3.4) is the one proportional to shear viscosity η and to first order in gradients one would obtain the Navier-Stokes result π µ ν = −2ησ µ ν . Similarly, (3.7) would give to first order in gradients, the bulk viscous pressure π bulk = −ζ∇ ρ u ρ . At second order in gradients, in particular the relaxation times τ shear and τ bulk come in and eqs. (3.4) and (3.7) mainly describe the dynamical relaxation of π µ ν and π bulk towards their Navier-Stokes values.
There are also additional, non-linear terms of second order with additional transport coefficients. More specific, the coefficients τ ππ , δ ππ , λ πΠ , δ ΠΠ and λ Ππ are formally of order O(Kn Re −1 ); while ϕ 7 , ϕ 6 , ϕ 1 and ϕ 3 are formally of order O(Re −2 ). All these terms can be understood as non-linear modifications of how π µ ν and π bulk relax to their asymptotic Navier-Stokes or equilibrium values. The terms of order O(Kn Re −1 ) contain one space-or time derivative acting on the dynamical fields, which may be taken to be energy density or temperature T , three independent components of fluid velocity u µ , five independent components of the shear stress π µ ν and the bulk viscous pressure π bulk . The terms of order O(Re −2 ) contain actually no derivatives.
Note that we have dropped in (3.4) contributions of order O(Kn 2 ), that correspond to non-linear terms in the transverse gradients, like for example (∇ ρ u ρ ) 2 or ω µ λ ω λν . Such terms are non-linear in derivatives of the fluid velocity and temperature and in principle appear naturally in a gradient expansion scheme. Corresponding transport coefficients have been computed in ref. [30,31]. However the problem with these terms is that they would transform the set of hyperbolic equations into a parabolic or mixed set of partial differential equations which would in general be neither causal (in the relativistic sense) nor stable with respect to small perturbations. As they stand, eqs. (3.4) and (3.7) constitute the most general set of equations for the evolution of shear stress and bulk viscous pressure at second order in the formal gradient expansion around local equilibrium but with only first order derivatives acting on the extended fluid fields , u µ , π µ ν and π bulk . These equations are hyperbolic and have therefore a chance to be causal and stable with respect to linear perturbations.

Coordinate system
We will now consider a more concrete situation of a relativistic nuclear collisions at high energies with a particular set of symmetries. We discuss here first coordinate systems that are particularly useful for our purpose.
Cartesian coordinates may be chosen such that the z-axis agrees with the beam axis, that the collision of the strongly Lorentz contracted nuclei takes place at t = z = 0, and that the center of the overlap region is at the origin of the transverse plane x = y = 0. We follow Bjorken in assuming that the dynamics is invariant under longitudinal boosts. This should be a good approximation at large collision energy and in the central longitudinal region.
It is therefore convenient to introduce the proper time τ = √ t 2 − z 2 and longitudinal rapidity η such that t = τ cosh(η) and z = τ sinh(η). Moreover, we introduce polar coordinates in the transverse plane such that x = r cos(φ) and y = r sin(φ). In these coordinates, the Minkowski space metric becomes (3.8) Below we will mainly concentrate on central collisions where the initial conditions are invariant under translations in rapidity, η → η + ∆η, and the azimuthal angle, φ → φ + ∆φ, corresponding to Bjorken boosts and azimuthal rotations, respectively. By symmetry reasons, this will then also be the case at later times. The dynamics is therefore effectively reduced from 3+1 space and time dimensions to the 1+1 dimensional subspace of proper time τ and transverse radius r.
To discuss the causal structure of the transverse expansion dynamics it will be convenient to introduce another parametrization that is related to τ and r by a conformal mapping similar to a Penrose diagram [32]. We write where the new coordinates σ and ρ replace τ and r. The monotonic function h(x) may be chosen to be a map from the interval (−1, 1) to (−∞, ∞) with the following properties As an example take with some length R > 0 and exponent α > 0. In this case, the infinite sector 0 ≤ τ, r < ∞ is parametrized by a finite region in terms of σ and ρ namely 0 ≤ σ < 1 and 0 ≤ ρ < 1 − σ. This coordinate system may also be useful for a numerical treatment because only a finite coordinate region must be considered.
In figure 1(a) we display the function h(x) in eq. (3.11) and in fig. 1(b) we show the coordinate lines of constant r and τ as a function of ρ and σ. Note that the point ρ = 1, σ = 0 corresponds to spatial infinity r → ∞ with fixed time τ and is labeled by i 0 in fig. 1(b). Similarly, the point σ = 1, ρ = 0 corresponds to τ → ∞ with fixed r, i.e. timelike infinity and is labeled by i + in figure 1(b). Finally, the line σ = 1 − ρ corresponds to τ → ∞, r → ∞ with fixed ratio τ /r. This corresponds to lightlike infinity and is labeled by J + in figure 1(b). The Minkowski space metric becomes Note that for dφ = dη = 0 the metric (3.12) is indeed related to the metric (3.8) by a conformal transformation, in particular light rays dτ = ±dr preserve their form and become dσ = ±dρ. This feature is particularly useful to investigate the issue of causality. We note as a side remark that (3.9) is not a conformal transformation in the full four-dimensional sense. Note also that the metric (3.8) and (3.12) are special cases of the ansatz g µν = diag(−g 11 , g 11 , g 22 , g 33 ), (3.13) with g 11 , g 22 and g 33 being functions of the time coordinate x 0 and radial coordinate x 1 . In the following it will often be convenient to work with this general case and restrict to the special cases of Bjorken coordinates τ , r, φ and η or the conformally related coordinates σ, ρ, φ and η at other places. We collect useful formulas for these coordinate choices, such as Christoffel symbols, projection operators etc. in appendix A.

Evolution equations for radial expansion
Assuming longitudinal boost invariance as well as azimuthal rotation symmetry, the fluid dynamic equations of motion become partial differential equations with one time coordinate x 0 and one radial coordinate x 1 . For the set of equations (3.3), (3.4) and (3.7), the equations are generically of hyperbolic type, and can be formulated as a first order system of equations for conveniently chosen independent fields. For our numerical treatment we will take as independent fields The connection between τ, r and σ, ρ is a conformal map in the two-dimensional space, similar to a Penrose diagram. The infinite half-plane 0 ≤ τ, r ≤ ∞ is covered by a finite region in terms of σ and ρ. The point i 0 at σ = 0 and ρ = 1 corresponds to spatial infinity r → ∞ at fixed τ . Similarly i + at σ = 1 and ρ = 0 corresponds to timelike infinity τ → ∞ with fixed finite r. The connecting line J + corresponds to lightlike infinity τ → ∞ and r → ∞ at fixed ratio τ /r.
• the radial fluid rapidity χ in terms of which the properly normalized fluid velocity is given by • the temperature T or actually its logarithm ln(T ), • two independent components of shear stress normalized by the enthalpy densityπ φ φ = π 3 3 /( + p) andπ η η = π 4 4 /( + p) to which the other components are related by eq. (A.9), • and the bulk viscous pressure normalized by enthalpy densityπ bulk = π bulk /( + p).
All these fluid variables are dimensionless and unconstrained fields, for which evolution equations follow from the set of equations (3.3), (3.4) and (3.7). These equations can be formulated as a quasi-linear differential equation with the combined field Φ = (χ, ln T,π φ φ ,π η η ,π bulk ) and coefficients A ij , B ij and C i that are functions of Φ but not of its derivatives. More concrete, we present the full set of equations of motion in appendix B and we obtain in particular the matrix associated to the time derivative 15) and the one associated to the radial derivative We have used here the abbreviations (3.17) The "source terms" C i are somewhat lengthy and we do not display them explicitly. However, they can easily be read off from the equations in appendix B. Together with initial values provided on an appropriate Cauchy surface (for example τ = τ 0 ), the differential equation (3.14) specifies the solution of the fluid equations completely.

Characteristics of the differential equations
As we have reviewed in section 2, quasi-linear systems of differential equations such as (3.14) can be characterized in terms of their characteristics. The characteristic curves can be seen as the lines along which information is transported, for example small perturbations, discontinuities, defects or shocks. The system is causal in the relativistic sense precisely if the characteristic velocities are smaller than (or, as a limit, equal to) the velocity of light.
More specific, the solution of the differential equation (3.14) at a particular space-time point x has a domain of dependence in the past of that point bounded by the characteristics with largest and smallest characteristic velocities [33]. This will be illustrated in more detail below.
In order to find the characteristics of eq. (3.14) we need to solve the eigenvalue problem where the eigenvalues λ (n) correspond to the characteristic velocities and w A direct calculation shows that they are given by where v = tanh(χ) is the fluid velocity andc is a modified sound velocity which we discuss below. Note that the characteristic velocities λ (1) and λ (2) correspond to "relativistic sums" of the fluid velocity v and the modified speed of soundc. Causality |λ (n) | ≤ 1 is guaranteed when |v| ≤ 1 and |c| ≤ 1.
Physically, we expect that the two characteristics with velocities λ (1) and λ (2) describe generalized sound propagation including viscous as well as non-linear effects. In contrast, the characteristics with velocities λ (3) , λ (4) and λ (5) should be understood as non-linear generalizations of diffusive modes for which information is propagated along the fluid flow lines.
Let us now discuss the modified sound velocityc. It can be written as where the ideal fluid velocity of sound is determined by the thermodynamic relation and the viscous and non-linear modification is parametrized by the combination (3.23) The second equation uses the dimensionless variables introduced in appendix A and the third equation uses the abbreviations (3.17). We note that causality requires large enough relaxation times τ shear and τ bulk for given shear viscosity η and bulk viscosity ζ. In fact, this was already known from the analysis of small perturbations around thermal equilibrium states by Hiscock and Lindblom [26]. In contrast to ref. [26], our analysis applies for the general, non-linear evolution in the specific situation of azimuthally symmetric and Bjorken boost invariant heavy ion collisions.
Note that the causality constraint formulated here is not only a constraint on thermodynamic and transport properties. At the non-linear level (in deviations from equilibrium) it involves the shear stress as parametrized in terms of π φ φ and π η η as well as the bulk viscous pressure π bulk . One needs to check for a specific solution of the fluid equations whether the causality constraint |c| ≤ 1 is satisfied, and we will do so for a typical heavy ion collision below. A set of left eigenvectors corresponding to (3.18)  4 Causality of the radial expansion

Causal initial condition and evolution
Let us now investigate the condition for a causal evolution in more detail. The characteristic velocities in (3.20) are below the speed of light and the system of first order hyperbolic equations (3.14) is accordingly causal, precisely if the modified sound velocityc in (3.21) is bounded by |c| ≤ 1. This in turn is a condition that involves thermodynamic quantities such as the ideal fluid velocity of sound (3.22), (ratios of) transport coefficients such as η/τ shear but also the components of shear stress π η η and π φ φ as well as the viscous pressure π bulk . The initialization of the fluid evolution includes also the specification of initial values for the shear stress and bulk viscous pressure. This means that the causality bounds |c| ≤ 1 is also a condition for a viable initial conditions.
To investigate how important this causality bound is in practice, we will now determine the modified sound velocityc for typical initial conditions relevant to high energy nuclear collisions.
For the numerical investigation we use the QCD equation of state provided by [35] and we adopt a temperature dependent shear-viscosity [36] calculated for Yang-Mills theory. We neglect the bulk viscous pressure. Among the second order transport coefficient we keep only τ shear and δ ππ , while other second order coefficients do not enter or are neglected. We have chosen the value of this transport coefficient such that as obtained from AdS/CFT calcualtions [11]. We initialize fluid dynamics on the hypersurface τ = τ 0 with vanishing radial fluid velocity χ(τ 0 ) = 0 and we choose the initial temperature profile according to the model discussed in ref [37], with a maximal temperature in the center of the fireball of T max = 0.4 GeV.
For the initial condition of the shear stress we compare four possibilities: we have setπ φ φ =π η η = 0 at the initialization time τ 0 = 0.6 fm/c. The causality constraints |c| ≤ 1 is satisfied there as well as at subsequent times. In (b) and (c) we have chosen Navier-Stokes initial conditions according to equation (4.2) at the initialization time τ 0 = 0.6 fm/c and τ 0 = 0.1 fm/c, respectively. One observes that in the former case causality is violated at early times and for large radii and in the latter case for early times at all radii. For later times one finds decreasingc such that |c| ≤ 1 becomes valid. In (d) we have used a modification of Navier-Stokes initial condition according to equation (4.3). The causality constraint is satisfied then at all times τ and radii r.

(d) A variant of the Navier-Stokes initial condition specified by
and similarly forπ φ φ .
For the four choices of initial conditions above we show the velocity parameterc as a function of the radius r in figure 2. We also plotc at later times τ as obtained from solving equations (3.14) numerically.
One observe that the initial condition (a) leads to |c| ≤ 1 for all times τ and radii r. Relativistic causality is indeed satisfied here. In contrast, the Navier-Stokes initial conditions (b) and (c) are not viable from a causality point of view. With the initialization time τ 0 = 0.6 fm/c corresponding to (b), the velocityc exceed the velocity of light for larger radii r at early time. At later time τ the shear stress relaxes andc decreases below the velocity of light. If the initialization time is chosen as τ 0 = 0.1 fm/c corresponding to (c), the causality bound is actually violated at all radii r for early time τ . The Navier-Stokes initial condition is therefore not viable at early times τ in the second order approximation for fluid dynamics (3.4), althought the solution relaxes towards smaller value ofc at late time. Finally, case (d) corresponds to a generalization of Navier-Stokes initial condition that involves also the additional transport coefficient λ 2 . As one can read off from 2(d), when used at the initialization time τ 0 = 0.6 fm/c, this leads to |c| ≤ 1 at all relevant times and radii r.
The violation of the causality bound for Navier-Stokes initial condition at early time can be also be highlighted directly from equation (3.23).The initial conditions can be written in a simple form and the causality constraint on the modified sound velocity |c| < 1 leads to the condition where we have used the abbreviation Γ = 2η 3( + p)τ shear .
In general, the right-hand side is a function of the initial temperature, but in order to have an estimation of magnitude of the bound, we can approximate the speed of sound by c 2 s = 1/3, while (4.1) leads to Γ ≈ 2/9 and for δ ππ we can take the conformal value δ ππ = 4/3τ shear . Inserting this values into the inequality (4.4) we obtain τ 0 τ shear ≥ 2 9 , (4.5) and restoring the units we have This simple estimation explains the differences between the initial condition (b) and (c) in the center of the fireball, where the first is initialized at τ 0 = 0.6 fm/c, consequently the causality condition is satisfied while (c) is initialized at τ 0 = 0.1 fm/c and the causality bound is violated.

Characteristic curves
It is instructive to study in more detail the characteristic curves defined by the characteristic velocities λ (m) in (3.20). More specific, the characteristic curves are defined as the solution of the differential equation  Figure 3. Characteristic curves in the coordinates σ, ρ related to Bjorken time τ and radius r by the conformal map (3.9) (see also fig. 3.2). In (a) we plot characteristic curves with characteristic velocity λ (1) in (3.20), in (b) with characteristic velocity λ (2) and in (c) with characteristic velocity λ (3) , which corresponds to the fluid velocity. For better orientation we also show curves corresponding to constant temperature as dashed lines. For a discussion of the significance of i 0 , i + and J + see figure 3.2. To calculate these curves we have used initial conditions for the shear stress according to eq. (4.3) .
For our purpose, it is particularly convenient to explore these curves in the σ, ρ coordinate system introduced in section (3.2). We show the result in fig. 3. The diagram in figure 3(a) shows the characteristic curves with velocity λ (1) , while figure 3(b) shows those with characteristic velocity λ (2) . Causality demands that these velocity be smaller than unity which is indeed the case. Figure 3(c) shows the characteristic curves corresponding to λ (3) = v. These curves can therefore also be understood as the fluid flow lines. For better  Figure 4. Domain of dependence Γ d and domain of influence Γ i of a space-time point x. The domain of dependence Γ d is the region in the past of x between the characteristic curves with velocities λ (1) and λ (2) , respectively. Hypothetical changes of field values in that region could modify the solution at the point x. In turn, the domain of influence Γ i is the region in the future of x bounded by the characteristic curves with velocities λ (1) and λ (2) where field values could be influenced by small changes at the point x. For better orientation we have also show the curve of constant temperature T = 0.145 GeV (dashed line) and the future and the past light cones of x (dotted lines). orientation we also show curves of constant temperature as the dashed lines in figure 3.

Domain of dependence and domain of influence
The characteristics curves with the smallest and the largest velocity form a boundary of the region of dependence in the past of a given space-time point x. We illustrate this in figure Fig 4 for a point on the line of constant temperature T = 0.145 GeV. The characteristic curves with velocity λ (1) and λ (2) going through this point are the boundaries of the region Γ d , the region of dependence. This is the region inside which a (hypothetical) change of the field Φ, in sense of a change of initial condition, can influence the value of Φ at the point x. On the other side, changing the initial condition outside of this region has no influence on Φ(x), as a consequence of causality, see also appendix C. For clarity we also show the past and future light cone of the point x and one can observe that indeed the domain of dependence is a subregion of the past light cone.
In a similar way, one can define the domain of influence of the point x as the region in its future, bounded by the characteristic curves with smallest and largest velocity, respectively. As the name suggests, this is the region where field values can be influenced by a small change in the value of Φ(x) at the point x. Again this region is a subregion of the corresponding light cone as required by causality.

Causality as a bound to applicability of fluids dynamics
As discussed in section 1 and 3, fluid dynamics is an effective theory for the dynamics of a quantum field theory that becomes valid for slow enough dynamics and small enough spatial gradients.
Although the framework of Israel-Stewart theory [7] or the more general DNMR [9] setup are going beyond a strict (Chapman-Enskog) derivative expansion in the thermodynamic fields and fluid velocity,  Figure 5. Ratio of different pressure P η /P r according to (4.8) in the central region of the fireball r = 0 as a function of Bjorken time τ . We compare different initial condition forπ η η andπ φ φ within the range allowed by causality. For panel 5(a) we have set τ shear = 3η/( + p), for the panel 5(b) τ shear = 30η/( + p).
it is clear that the expansion assumes small Knudsen number and inverse Reynolds number according to the definition in section 3. In particular, this assumes that the components of the shear stress and bulk viscous pressure are small compared to the thermodynamic enthalpy density.
It is often stated in literature that the application of the fluid approximation needs an approximately isotropic energy-momentum tensor T µν . Indeed, deviations from isotropy in the decomposition (3.3) are parametrized by the shear stress. At least superficially, a small inverse Reynolds number Re −1 implies also a close-to-isotropic energy-momentum tensor. However, is not clear where the bound of applicability is precisely. Oftentimes, this is discussed in terms of ratios of different "pressure components", for example (assuming π bulk = 0 for simplicity) , (4.8) and it is believed that such ratios should not deviate too much from unity in order to not leave the range of applicability of fluid dynamics. However, it is not so clear where precisely the boundary is. In particular for p the above ratio (4.8) could actually deviate for unity substantially while the ratiosπ η η = π η η /( + p) andπ φ φ = π φ φ /( + p) might be below one. A more severe bound forπ φ φ andπ η η comes from the causality constraint |c| ≤ 1 in terms of the relations (3.21), (3.22) and (3.23). It is interesting to investigate the range for the ratios such as the in eq. (4.8) for solution of fluids dynamics that satisfy the causality bound. We have done this in fig. 5(a) where we plot the ratio P η /P r in the center of the fireball for different initial conditions ofπ φ φ andπ η η within the causality bound. One observes that a rather large range of values is allowed by causality and that after a rather short time of order 1 fm/c they approach the same attractor solution which then evolves towards late time. The speed of approach to this attractor is actually governed by the shear stress relaxation time τ shear . This can be see by comparison to fig. 5(b) where we have chosen τ shear to be larger by a factor 10. The approach to the attractor is accordingly slower. In Fig. 5 we have chosen the initialization time τ 0 = 0.6 fm/c; for earlier initialization times τ 0 the range of P η /P r allowed by causality is even larger.

Conclusions
We have discussed here the causality of relativistic fluid dynamics for high-energy nuclear collisions. The fluid equations we use contain all possible terms up to second order in a formal gradient expansion compatible with symmetry and with hyperbolicity [9]. We concentrate on an expanding fireball with longitudinal Bjorken boost and azimuthal rotation symmetry. The causality of the dynamics in the reduced configuration space of Bjorken time τ and radius r can be conveniently discussed in terms of characteristics. We found five characteristic velocities for the five independent fluid fields. Three of them are degenerate and corresponds to the fluid velocity. The remaining two differ from this by a modified sound velocity (see eq. (3.21), (3.22) and (3.23)) which contains the standard ideal fluid velocity of sound but also nonlinear modifications due to the dissipative terms. In particular, the modified sound velocity depends not only on thermodynamic variables and transport properties, but also on the components of the shear stress tensor and bulk viscous pressure. Causality is in this sense state-dependent.
In a subsequent step we investigate whether the modified sound velocity remains below the speed of light for characteristic situations that occur in the context of high-energy nuclear collisions. Depending on the initial condition, we found that the causality constraint can be violated at early times, which indicates a violation of the applicability conditions for the fluid approximation. In particular, such violations occur for so-called Navier-Stokes initial conditions that have often been used in phenomenological investigations. We also show that modifications of these initial conditions lead to a dynamical evolution that does not violate the causality constraint. Finally, we investigate the ratio of longitudinal and radial pressure in the center of the fireball for a range of initial conditions that does not lead to any violations of causality. We find that this ratio can vary substantially at early times but approaches an attractor solution quickly at late times where the "speed of approach" is determined by the shear stress relaxation time as expected.
The most important conclusion from our findings is that causality poses a bound on the applicability of relativistic fluid dynamics. This helps for concrete phenomenological applications of the formalism but also improves the more conceptual understanding of its foundations.
The coordinates σ and ρ are related to the more familiar longitudinal proper time τ = √ t 2 − z 2 and transverse radius r = x 2 + y 2 by The non-vanishing Christoffel symbols are and the square root of the metric determinant is Note that in terms of σ and ρ, the half plane 0 ≤ τ, r < ∞ corresponds to 0 ≤ σ, ρ < 1 with ρ + σ < 1. However, it is sometimes useful to analytically continue to negative ρ with |ρ| < 1 − σ. One can obviously recover the original Bjorken coordinate system with the trivial choice h(x) = x such that σ = τ and ρ = r. In a 1+1 dimensional situation it is convenient to parametrize the fluid velocity as with radial fluid rapidity field χ. The projector orthogonal to the fluid velocity is The fluid velocity divergence is (A.7) We also need the tensor ∇ µ u ν . It has the following non-vanishing components √ g 11 ∇ 0 u 0 =∂ 0 cosh(χ) + sinh(χ) ∂ 1 ln ( √ g 11 ) , √ g 11 ∇ 0 u 1 =∂ 0 sinh(χ) + cosh(χ) ∂ 1 ln ( √ g 11 ) , (A.8) The shear stress π µν has only two independent components and for convenience they can be chosen as the dimensionless ratiosπ η η = π 3 3 /( + p) andπ φ φ = π 2 2 /( + p). The other non-vanishing components are related to these by We also parametrize the bulk viscous pressure by a dimensionless variable as π bulk =π bulk ( + p).

B Equations of motion
In this appendix we compile the relativistic fluid dynamic equations of motion for a central heavy ion collision with longitudinal and transverse expansion. The equations are partial differential equations involving a time coordinate x 0 and a radial coordinate x 1 . The evolution equations for energy density and fluid velocity follow directly from energy-momentum conservation. The energy density, pressure and other thermodynamic and transport properties close to local equilibrium can be expressed in terms of a single independent thermodynamic variable, which we take to be temperature.

C Energy inequalities and uniqueness of the solution
In this appendix we explain in more details the notion of causality of hyperbolic equations and its relation with characteristic curves. As we discuss in the main text, the evolution is causal if the domain of dependence Γ d of a given point P = (x 0 , x 1 ) is contained in the past light cone of that point. The domain of dependence Γ d is defined as the region in the past of the point P bounded by the two extreme characteristic curves, i.e. the two integral curves with maximal and minimal velocities. Any perturbation of the initial data that vanishes in this region does not change the value of the solution at the point P . This property is due to the existence of a so-called "energy inequality" that bounds the spatial average of the solution to the spatial average of the initial data. To be more specific consider a perturbation δΦ, solution of the linearized version of the partial differential equation (3.14) around a generic background Φ 0 , and consider the domain of dependence of the point P shown in figure 6: C l is the characteristic curves with largest velocity and C r with the smallest velocity, the line x 0 = 0 is the where the initial data are provided and where A l , A r , P l and P r are the intersection points of the characteristic curves with constant time lines shown int the figure 6, and µ is some positive constant. The first consequence of this inequality is that the solution of the linearized equation δΦ can grow at most exponential with respect to the initial data. The second consequence, and most interesting for our purpose, is that if the initial data vanish, i.e. for all times x 0 = t, and in particular at the point P , therefore δΦ(P ) = 0. As a consequence, any perturbation δφ outside of the domain of dependence can not modify the value of the solution at the point P , to linear order.
In the rest of this appendix we recall briefly the proof of the energy inequalities that can be found in [33] applied to equation (3.14). Consider a small perturbation δΦ around a generic solution Φ 0 (x) of (3.14), ∂ 0 δΦ j + M jk ∂ 1 δΦ k + S jk δΦ k = 0. .

(C.5)
This set of equations is linear and hyperbolic and the matrix M has a set of left eigenvector given by (3.24). Eq. (C.4) is not symmetric-hyperbolic because M in general is not a symmetric matrix. However, The source term can be expressed in terms of δJ (m) as well using the right eigenvectors r where Λ = diag(λ (1) , λ (2) , · · · , λ (n) ) and the scalar product (·, ·) is the Euclidean one, defined on the vector space of the perturbations. The unknown vector δJ can be rescaled, δJ = e µx 0 δJ which leads to The quadratic form on the right hand side can be taken as a negative definite, if we select a sufficiently large value for the constant µ. Consequently we can write an inequities for the left hand side, 1 2 ∂ 0 (δJ, δJ) + 1 2 ∂ 1 (δJ, ΛδJ) ≤ 0. (C.10) For a given point P = (x 0 , x 1 ) is possible to compute the characteristic curves x 1 m (x 0 ) solving the differential equation dx 1 /dx 0 = λ (m) , this is possible because the eigenvalue λ (m) depends only on the background solution. Consider now the trapezoidal domain Γ bound by the lines x 0 = 0, x 0 = t and the two outer characteristic curves that arrive at the point P (see figure 6). Integrating the inequality (C.10) over this domain, we obtain the inequality Parametrizing the curves C l and C r with x 0 and considering the different orientations of the two curves C l and C r respect to this parameter the righthand side of the inequality (C.12) can be written as The two characteristic curves were chosen such that C l has the maximal velocity and C r the minimal velocity, which corresponds to the extreme eigenvalues of Λ,