Fluid dynamic propagation of initial baryon number perturbations on a Bjorken flow background

Baryon number density perturbations offer a possible route to experimentally measure baryon number susceptibilities and heat conductivity of the quark gluon plasma. We study the fluid dynamical evolution of local and event-by-event fluctuations of baryon number density, flow velocity and energy density on top of a (generalized) Bjorken expansion. To that end we use a background-fluctuation splitting and a Bessel-Fourier decomposition for the fluctuating part of the fluid dynamical fields with respect to the azimuthal angle, the radius in the transverse plane and rapidity. We examine how the time evolution of linear perturbations depends on the equation of state as well as on shear viscosity, bulk viscosity and heat conductivity for modes with different azimuthal, radial and rapidity wave numbers. Finally we discuss how this information is accessible to experiments in terms of the transverse and rapidity dependence of correlation functions for baryonic particles in high energy nuclear collisions.


I. INTRODUCTION
One of the most important goals of the experimental program of high energy nuclear collisions is to determine the transport and thermodynamical properties of QCD as a function of temperature T and baryon chemical potential µ. During the past few decades, the experimental data measured at the Relativistic Heavy Ion Collider (RHIC) at the Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN in Geneva, has shown collective behavior of the QCD matter created after the collision of heavy nuclei at high energies [1][2][3][4][5][6][7]. The low momentum region of the transverse hadron spectra and the two particle correlation functions are well described by relativistic viscous fluid dynamics with a very small value of the shear viscosity over entropy ratio 1 . These results have been taken as evidence for the production of an almost perfect liquid, a strongly coupled quark gluon plasma.
The hydrodynamic modeling of heavy ion collisions solves on an event-by-event basis the relativistic fluid equations corresponding to energy-momentum conservation laws together with the so called constitutive relations for the shear viscous tensor and bulk pressure. Within this approach, little attention has been paid to the possible role of the baryon density n and/or baryon chemical potential µ. At high energies, this is justified because n and µ are very small, at least in the midrapidity region. However, interesting physics could be probed by investigating event-by-event fluctuations in the local baryon number density.
Baryon number fluctuations have been mainly discussed in the context of heavy ion collisions at lower energy where larger values of µ can be realized. Interesting features of the QCD phase diagram can emerge there [10]. Different effective models have predicted the existence of a first-order phase boundary that separates hadronic matter from the quark gluon plasma at larger values of the baryon chemical potential. This boundary comes to an end at some critical values of the temperature T c and baryon chemical potential µ c . Right now there is no conclusive evidence for the location of a critical point in the T − µ plane from lattice QCD calculations at finite baryon density [11].
On the other hand, in heavy ion collisions it has been proposed to study second and higher order cumulants of particle multiplicity distributions as a function of the center of mass energy √ s [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. From thermodynamic considerations, it is expected that these moments scale with the correlation length which is expected to become large near the QCD critical point [16-18, 26, 28, 32]. Possible signs of the critical point have been measured at RHIC but at present these do not provide a conclusive evidence [33][34][35][36]. If the expanding fireball of nuclear matter passes through a critical region (close to a critical point), one can extract information about the equation of state and the critical behavior of transport coefficients from the particle spectra formed at the freeze-out surface. It is important to determine whether the possible signatures of the critical point can survive the entire evolution of the expanding fireball.
In the fluid dynamic framework, different aspects of the evolution of the fireball can change the pattern expected from purely thermodynamic considerations. Thermodynamic fluctuations are in principle part of a fluid dynamic description, at least in an extended sense where one accounts also for noise. Fluctuations evolve in time and space during the expansion of the fireball and thus, these are indeed effected by the equation of state and specially the transport coefficients such as the viscosities and conductivities. Close to equilibrium, there is also a deep theoretical connection between thermodynamic fluctuations in fluid dynamic fields and dissipative transport properties as stated by the fluctuation-dissipation theorem; see, e. g., Refs. [37][38][39]. In the vicinity of the critical point, heat conductivity κ as well as the shear and bulk viscosities η and ζ show critical behavior [40][41][42][43].
Besides genuine thermodynamic fluctuations (or noise), there is another possible source of fluctuations in the fluid dynamic approach to heavy ion collisions. These are the fluctuations already present in the initial state when the fluid dynamic treatment becomes valid. Their origin can be either the substructure of the colliding nuclei or the farfrom-equilibrium dynamics preceding a fluid dynamic regime. This kind of initial state perturbation is particularly important for energy and/or entropy density. Fluctuations in the geometric distribution of nucleons within a nucleus lead to initial density perturbations which -after a fluid dynamical evolution -determine the spectrum of harmonic flow coefficients and, for example the form of the two-particle correlation function ("the ridge") in heavy ion collisions.
In a very similar way to fluctuations in the initial energy density, one can also expect, for example from a Glaubertype description of the initial state, initial fluctuations in the baryon number density. Indeed, baryon number density carried by protons and neutrons is presumably not distributed homogeneously within a nucleus and fluctuates locally and from event to event. In addition, the baryons and anti-baryons produced by pair production directly after the collision are subject to some local and event-by-event fluctuations [44].
In order to discriminate the effects associated to the thermodynamic fluctuations from the initial state fluctuations, it is necessary to understand their space-time evolution. In the present work we will concentrate mainly on the dynamics of initial state fluctuations although parts of our formalism are relevant also for the evolution of thermodynamic fluctuations. Initial state fluctuations are interesting on their own. For instance, the evolution of the initial perturbations of energy density depends on the viscosities, in particular shear viscosity. In a similar way, the evolution of baryon number density depends on heat conductivity (in the Landau frame one may see heat conductivity equivalently as baryon number diffusion). If one has a theoretical understanding of initial state perturbations in baryon number density and their fluid dynamic evolution, it is possible to study their consequences for particle spectra at freeze-out. Provided possible signals are large enough to be seen within the constraints set by finite statistics, there could be a possibility to constrain the heat conductivity of the quark gluon plasma from experimental data. This would be very interesting for not only low energy collision experiments which aim at exploring the QCD phase diagram, but also at RHIC and LHC energies where baryon number diffusion could be another characteristic of the quark-gluon plasma.
As a first step in this direction we study here the fluid dynamic propagation of local and event-by-event fluctuations of the baryon number density, flow velocity and energy density. These fluctuations propagate on top of a hydrodynamical background which for simplicity, we consider to be described by Bjorken's model [45] (which includes finite baryon number density). In order to study the fluid dynamic propagation of perturbations we use a background-fluctuation splitting and a Bessel-Fourier decomposition for the fluctuating part of the fluid fields [46][47][48][49][50][51][52][53]. We derive the evolution equations of the linear fluctuations and solve them for different initial conditions, values of the transport coefficients and equation of state.
This work is organized as follows. In Sec. II we review briefly the theory of relativistic fluid dynamics at finite chemical potential putting emphasis on the role of the equation of state and current estimates of the transport coefficients in the strong and weakly coupling regimes. The main features of the temporal evolution of the background fields are discussed in Sec. III. In Sec. IV we formulate the theory of linear perturbations on top of this evolving background and discuss numerical solutions. In Sec. V we draw some conclusions for a potential experimental observable, the correlation function of net baryon number as a function of azimuthal angles and rapidity. General conclusions are presented in Sec. VI. Some technical details of our calculations are presented in Appendixes A and B, respectively.

II. RELATIVISTIC FLUID DYNAMICS WITH A GLOBALLY CONSERVED CHARGE
We consider a relativistic fluid with one globally conserved quantum number current (baryonic number for our purposes). The energy-momentum tensor and number current are Here, is the energy density, u µ is the fluid velocity, π µν is the shear stress tensor, π bulk is the bulk viscous pressure, n is the particle density and ν µ is the particle diffusion current. We choose the signature of the metric g µν to be (−, +, +, +) and the projector orthogonal to the fluid velocity is The fluid velocity is normalized to u µ u µ = −1. We work in the Landau frame where the fluid velocity is chosen such that u µ T µν = − u ν . The shear stress tensor is transverse to the fluid velocity, The shear stress tensor is also symmetric and traceless. The particle number density is defined by n = −u µ N µ such that the diffusion current is orthogonal to the fluid velocity, u µ ν µ = 0. It is clear that an arbitrary (symmetric) energy-momentum tensor T µν (with a time-like eigenvector) and current N µ can be written in the above form. The decomposition becomes unique by requiring that the pressure p is related to the energy density and the baryon density n by the same relation as in thermodynamic equilibrium, i.e. by an equation of state p = p( , n) .

(5)
Here we have introduced the comoving derivative defined as D = u µ ∇ µ .
To close the evolution equations (5) one needs expressions for π bulk , π µν and ν µ . Within the formalism of fluid dynamics one writes these objects as a derivative expansion in terms of the fluid velocity u µ and thermodynamic variables , n. In the present work we concentrate for simplicity on the first order of this expansion. One should keep in mind that terms of second order are expected to improve the results quantitatively and are in general needed for an acceptable causal structure and linear stability [54,55].
The constitutive relation for the shear stress is where η is the shear viscosity transport coefficient. The bulk viscous pressure is obtained from the following expression where ζ is the bulk viscosity and θ is the expansion scalar. Finally, the particle diffusion current is where κ is the heat conductivity. In the last equation we have introduced the chemical potential µ, which is conjugate to the baryon density n, and the temperature T . In summary, the hydrodynamic equations at this stage involve the fluid velocity u µ (with three independent components), the energy density , pressure p, baryon density n, baryon chemical potential µ, temperature T as well as the shear viscosity η, bulk viscosity ζ and the thermal conductivity κ. Only two thermodynamic variables are independent and they also determine the transport properties η, ζ and κ. In a non-equilibrium situation only energy density = u µ u ν T µν and baryon number density n = −u µ N µ are directly related to the physical energy-momentum tensor T µν and number current N µ . All other thermodynamic variables are defined indirectly via their relation to and n in thermal equilibrium.
For the practical calculations one is in principle free to use any set of independent thermodynamic variables. The form of Eqs. (5) suggests the use of the energy density and baryon density n. However, because most microscopic calculations are done in the grand canonical ensemble, the thermodynamic equation of state and the transport coefficients are usually obtained as a function of the temperature T and chemical potential µ, for example p = p(T, µ). Thus, it can be advantageous to use T and µ as independent variables in fluid dynamics, as well. This avoids the inversion of functions which can be numerically difficult. One should keep in mind that T and µ in a non-equilibrium situation are defined via their relation to and n. Eq. (5) can be transformed using thermodynamic relations compiled in Appendix A. The evolution equation for energy density becomes where we have now used the constitutive relations (6) and (7). The evolution equation for the fluid velocity is now of the form and finally, the particle number conservation law becomes Note that Eqs. (9) and (11) form a linear system of equations that can be solved for DT = u α ∂ α T and Dµ = u α ∂ α µ as long as To solve the fluid dynamic equations we will use a background-fluctuation splitting. To this end we write the fluid dynamic fields as and similar for the other fields. We are interested in perturbations δu µ , δ etc. that are small enough so that only linear terms in the evolution equations need to be kept. The background fieldsū µ ,¯ etc. satisfy the fluid dynamic equations (5) while the perturbations satisfy linear equations that depend on the background solution. We derive these linearized equations for arbitrary background fields in Appendix B. The structure of the linearized equations permits us to simply use δ , δn and three independent components of the fluid velocity as variables (the fourth component of the fluid velocity follows from the constraintū µ δu µ = 0). However, all the background-dependent thermodynamic quantities can be expressed in terms ofT andμ. Useful thermodynamic relations for this purpose are compiled in Appendix A.
In the rest of this section we briefly discuss some simple parametrization of the thermodynamic equations of state p(T, µ) and transport properties η(T, µ), ζ(T, µ), κ(T, µ). We emphasize that our formalism can be used for an arbitrary form of these functions once these have been determined from a particular microscopic description.

A. Equation of state
The fluid hydrodynamical equations require an equation of state (EOS) p(T, µ) as an input. In principle, the equation of state can be calculated from the inherent quantum field theory associated to a particular system but this is a formidable task. In recent years there have been important advances to determine analytically and numerically the thermodynamical properties of QCD at high temperatures and chemical potential by considering effective thermal field theories [56][57][58][59][60][61][62] while in the low temperature and chemical potential regimes one expects that a non-interacting hadron resonance gas provides a reasonably good approximation [63].
At intermediate temperatures, non-perturbative methods are needed to describe the transition which separates the hadronic, confined phase and the quark-gluon plasma (QGP) phase. While several studies of lattice QCD simulations are available at the moment at vanishing chemical potential µ = 0 2 , at µ > 0 lattice simulations are not possible due to the sign problem. However, different alternatives have been studied in order to circumvent this problem such as reweighting [64], Taylor-expansion in µ [65][66][67][68][69][70][71][72], analytic continuation from imaginary µ [73], the density of states method, or using the canonical ensemble. Of course, each of these methods have their advantages and disadvantages.
One of the main goals in the analysis of fluid dynamic fluctuations and their propagation is to provide a phenomenological determination of the equation of state (or at least of some of its properties). In the derivation of the evolution equations for the background and fluctuating fields we shall keep the equation of state p(T, µ) unspecified as far as possible in analytic expressions. For some numerical calculations and illustrations we use the simplest possible case, a non-interacting gas of N F massless quarks that come in N C colors and N 2 C − 1 gluons, where we use the abbreviations The baryon chemical potential µ measures the net baryon density of the system. In our convention, quarks carry baryon number charge 1/3 and anti-quarks −1/3. Corrections to the ideal EOS arise as a consequence of interactions and the breaking of conformal invariance by dimensional transmutation and non-zero quark masses. They are most important at low temperatures. We follow here the Wuppertal collaboration which has parametrized the QCD equation of state for finite chemical potential in Transport coefficient Weakly-coupled QCD Strongly-coupled theories  [75][76][77][78][79] and strongly coupled theories with holographic duals [80][81][82]. See text for discussion.
terms of a Taylor expansion [71]. The leading order expression for the trace anomaly or QCD interaction measure where I(T, 0) is the interaction measured at µ = 0 and χ 2 (T ) is the leading-order Taylor coefficient. Both terms, I(T, 0) and χ 2 (T ), can be parametrized analytically as [71,74] I(T, 0) where t = T /(0.  [71,74]. The pressure at finite µ is given by At µ = 0 the relation between the pressure and the trace anomaly (17a) is All other thermodynamic quantities can be derived from p(T, µ) using the standard relations (compiled in appendix A). The equation of state (18) with the above parametrization is valid for small chemical potentials µ/T < 3 in the temperature window 0 < T < 400 MeV. We will use Eq. (18) to study the influence of the EOS for the dynamics of the background fluid dynamic fields.

B. Transport coefficients
In addition to the thermodynamic equation of state, the fluid dynamical description needs as an input transport coefficients. These can either be determined experimentally, or, if a microscopic underlying theory is known, they can at least in principle be calculated as a function of the thermodynamic variables via Kubo relations. In this section we briefly summarize the current theoretical knowledge for the shear and bulk viscosities and thermal conductivity of QCD and related theories, both in weakly and strongly coupled regimes 3 .

C. Weak coupling regime
When the interaction strength is small, effective thermal field theory methods allows us to calculate the transport coefficients. For weakly coupled QCD in the high temperature and vanishing chemical potential regime, the leading logarithmic result for the shear viscosity is [75][76][77] where g is the strong coupling constant. In the previous expression k is a constant that depends on the number of fermions species [75][76][77]. Arnold et al. showed that at leading log accuracy and for high temperatures with vanishing chemical potential there is an approximate scaling between the shear (η) and bulk (ζ) viscosities for weakly coupled QCD [77] where c 2 s = dp/d is the speed of sound. A similar expression was first derived by Weinberg for a gas of photons [84]. To date there is no complete leading logarithmic calculation of the heat conductivity κ(T, µ) and so far only two estimates of κ(T, µ) have been provided in the literature for different kinematic regions of the T − µ plane [78,79] where F (T, m D ) is a function that depends on the temperature and the Debye screening mass m D (see Ref. [78] for details). In the case of small chemical potential, the proportionality constant C depends on the number of flavors and the gauge group [79]. In the limit where µ → 0 the heat conductivity κ ∼ µ −2 is divergent. However, the particle diffusion current (8) remains finite [79]. In the context of relativistic kinetic theory, some general expressions for the transport coefficients with constant cross section or within the relaxation time approximation have been derived recently [85][86][87][88][89][90]. However, these calculations do not take into account the quantum screening effects of the QCD plasma.
Despite relatively large uncertainties, experimental results indicate that the value of the shear viscosity over the entropy ratio η/s is smaller than the one calculated from weakly coupled QCD (20) [8,9]. For the case of the bulk viscosity the situation is less clear: the uncertainties in its experimental determination are even larger (see Ref. [91] and references therein). In addition, there are no experimental constraints for the value of heat conductivity in high energy-nuclear collisions so far.

D. Strong coupling regime
From the previous discussion it is clear that at this moment perturbative QCD calculations of the transport coefficients are not completely under control for all the possible physical values of the temperature and chemical potential. On the other side, there are certain classes of strongly interacting theories where transport coefficients can be determined for almost all values of T and µ. These are field theories with known gravitational duals where the computations can be done via the anti-de Sitter/conformal field theory (AdS/CFT) correspondence. Despite the fact that those theories are not equivalent to QCD, they share some qualitative aspects with it and thus, these theories might provide some guidance in the regimes where pQCD calculations are not reliable. We take here a pragmatical approach and consider the estimates of the transport coefficients based on holographic calculations as toy models which allow us to study the propagation of perturbations in fluid dynamic fields. For large t'Hooft coupling and for N = 4 SYM theory, holographic methods give the well known result [80], This result holds also for any holographic theory at sufficiently large coupling and number of colors as long as the theory is spatially isotropic. This relation for η/s holds even in the presence of non-zero chemical potential [81]. Initially this result was conjectured to be an universal lower bound but today there is evidence showing that this relation does not hold in general [92][93][94][95][96][97][98][99]. Incidentally, the value of the shear viscosity extracted from experiments in high energy nuclear collisions is closer to the one predicted for strongly coupled theories (23) than the one calculated in weakly coupled QCD (20) (see Ref. [8] for a recent review).
The shear viscosity has also been calculated for pure Yang-Mills theory using lattice gauge theory for specific values of temperature [100,101]. The estimated values for η/s are somewhat above the AdS/CFT values. Similarly, η/s as a function of temperature for vanishing baryon chemical potential has also been estimated for Yang-Mills theory as well as QCD by using diagrammatic functional relations and gluon spectral functions obtained by numerical analytic continuation from Euclidean quantum field theory [102,103]. The minimal value for QCD was found to be η/s ≈ 0.17 For holographic theories that deviate from conformal behavior the bulk viscosity has also been calculated [82] As in the case of the shear viscosity value (23) this relation holds for certain theories with finite chemical potential [105] but it is not an universal bound [106]. By comparing the scalings between ζ and η, Eqs. (24) and (21), one observes that they differ in the strong and weak coupling regime. This mismatch between both parametrizations is currently not understood. In the case of the thermal conductivity κ, the calculations for strongly coupled plasmas with finite chemical potential give the following result [81] κ which is an analog of the Wiedemann-Franz law [107] 5 . As in the weakly coupled case (22), the heat conductivity is divergent ∼ µ −2 while the particle diffusion current (8) is finite. Recently the temperature-dependence of the first and second order transport coefficients have been studied in a particular holographic model [109]. We summarize the discussion presented in this section in Table I, where we show the estimates of the transport coefficients in both strong and weak coupling regimes. Mainly for reasons of simplicity, we shall concentrate here on the parametrizations of the transport coefficients in the strong coupling regime Eqs. (23), (24) and (25) for our numerical calculatons. Another advantage of using the parametrization of strongly coupled theories is that both transport coefficients, the bulk viscosity ζ and the heat conductivity κ, are proportional to the shear viscosity η and thus, one can not only study the effect of the dissipative corrections but also one can investigate the 'weak' and 'strong' regimes by varying the values of η/s. We keep the functions η(T, µ), ζ(T, µ), κ(T, µ) unspecified as far as possible in our analytic calculations.

III. BJORKEN BOOST INVARIANT SOLUTION
In this section we study the solutions of the fluid dynamical equations for a quark-gluon plasma undergoing boost invariant longitudinal expansion. We assume translational and rotational symmetry in the transverse plane and arrive at a simple model for the early stages of a heavy ion collision first studied by Bjorken [45]. Our analysis is extended to the case where there is a non-vanishing baryon number density. The relatively simple homogeneous solutions will also serve as a background for a more elaborate discussion of perturbations around it in Sec. IV.
It is convenient to change from Cartesian coordinates x µ = (t, 3 is the longitudinal proper time, η = arctanh(x 3 /t) is the longitudinal (space) rapidity and r and φ are the usual polar coordinates in the transverse plane. The metric in the Milne coordinates is g µν = diag(−1, 1, r 2 , τ 2 ). The main advantage of using these coordinate systems is that the symmetries of the Bjorken solution are explicitly manifest. Specifically, the symmetry group ISO(2)⊗SO(1, 1)⊗Z 2 consists of translations and rotations in the transverse plane, longitudinal boosts η → η + ∆η and reflections η → −η [110]. The Bjorken flow velocity u µ = (1, 0, 0, 0) is the only invariant unit vector and the symmetry also implies that all fluid dynamic fields depend only on the longitudinal proper time τ [45].
From Eqs. (5) one finds that the evolution equations for energy density and particle number density are We have used here the Christoffel symbols of the Milne coordinate system. The non-vanishing ones are Γ η τ η = with σ µν σ µν = 2 3τ 2 . The expansion scalar is θ = 1 τ and the projector orthogonal to the fluid velocity is ∆ µ ν = diag(0, 1, 1, 1). The particle diffusion current ν µ (8) is a vector orthogonal to u µ and therefore vanishes exactly for the Bjorken flow.
While the particle number density is simply diluted by the one-dimensional expansion, the evolution of energy density in (26) contains an additional loss term from the thermodynamic work done by the expansion and a gain term from shear and bulk viscous effects. After the variable change to T and µ eq. (26) becomes We observe that the size of viscous corrections to an isentropic expansion is determined by the parameter Formally, the gradient expansion underlying viscous fluid dynamics can be used for γ 1. Note that for a given thermodynamic equation of state p(T, µ) and viscosities η(T, µ), ζ(T, µ) one can solve the two coupled ordinary differential equations (27).
In the remainder of this section we discuss as a simple illustrative example the equation of state of an ideal gas of massless quarks and gluons in Eq. (14). The evolution equations (27) for the temperature T and chemical potential µ become where the coefficients a 1 , a 2 and a 3 are given in Eq. (15). Note that we use conventions where µ is the chemical potential for baryons, and the chemical potential for quarks is µ q = µ/3. Let us first discuss some interesting limiting cases of Eqs. (29): 1. Ideal fluid dynamic expansion. When shear and bulk viscosities vanish, η = ζ = 0, the temperature and the chemical potential decouple from each other. This allows us to solve Eqs. (29) exactly, which gives The scaling solution of the temperature is not modified by the presence of the chemical potential and it coincides with the well known result found by Bjorken [45].

2.
Vanishing chemical potential. The point with µ = 0 corresponds to a (partial) fixed point of the evolution equations (29) with extended symmetry (baryon number parity). The evolution equation for temperature becomes where γ is given by Eq. (28). For vanishing bulk viscosity, ζ = 0, and constant ratio η/s, the exact solution to the previous equation is [111][112][113][114] T Viscous corrections are relevant only at early times where velocity gradients are large while at late times these are suppressed and thus, T (τ ) ∼ τ −1/3 .

3.
Small chemical potential. For µ/T 1 the dynamics of T is approximmately determined by Eq. (31) while the evolution equation for µ is The viscous effects (encoded in the parameter γ) have the tendency to accelerate the decrease of µ due to the expansion. This is in contrast to the temperature where viscosity has the opposite effect. To lowest order in η/s, the solution of (33) is 4. Small temperature. For T /µ 1 the evolution equation for the chemical potential is the one of eq. (30b) with a simple scaling solution. For the temperature we obtain to lowest order in T /µ which has a solution similar to Eq. (32) when η/s and ζ/s have constant values. If one chooses T (τ 0 ) = 0 as initial condition the solution to Eq. (35) becomes Even if the temperature vanishes initially, the system is heated up due to shear and bulk dissipative effects. In contrast to µ = 0, vanishing temperature T = 0 does not correspond to a (partial) fixed point of the evolution.
Let us now consider the evolution equations (29) in the general case where we find their solution numerically. In First we discuss the properties of the numerical solutions of Eqs. (27) for the ideal (and massless) EOS (14). For both variables, T and µ, the effect of viscous corrections are more relevant during the early stages of the expansion while at late times their effects are negligible as expected. In the left panel of Fig. 1 we see that the viscosity reduces the effect of the longitudinal expansion on the temperature. This is simply the expected heating by dissipative effects. At the final time τ f = 10 fm/c the temperature is larger by values of the order of 10% for η/s = 2/(4π) compared to the ideal fluid expansion. For the chemical potential we find that the inclusion of dissipative corrections has the opposite effect, i.e. the chemical potentials decrease faster in the viscous case. This is clearly seen in the right panel of Fig. 1 when comparing the final values of the chemical potential µ(τ f ). The changes with respect to the ideal fluid expansion are also somewhat larger, of the order of 15% for η/s = 2/(4π).
When using the lattice-based EOS (18) we find that the numerical solutions of Eqs. (27) for T and µ are qualitatively similar to the ones obtained from the ideal EOS during the early stages of the evolution. As a function of time, the temperature is always decreasing and the dissipative corrections are larger at early times than at late times. The chemical potential decreases faster for larger values of the shear viscosity. For the lattice EOS, the changes induced by the dissipative corrections are on the order of 8 − 15%.
Interestingly, the evolution of µ with time differs substantially between the two choices for the equation of state. In the right panel of Fig. 1 one observes that the decrease with time is much weaker for the lattice EOS than for the ideal EOS. At the freeze-out time τ f = 10 fm/c and for vanishing η/s, one has µ(τ f ) ≈ 0.29 GeV for the lattice EOS while µ(τ f ) ≈ 0.12 GeV for the ideal EOS. The difference between those values increases slightly for finite values of η/s. Moreover, at late times µ increases slowly (and somewhat more for larger values of η/s). Our numerical results show also that when using the lattice EOS the values of the temperature are somewhat larger than for the ideal EOS specially at late times.
In Fig. 2 we show the Bjorken flow trajectories in the plane of chemical potential µ and temperature T for the ideal EOS (14)  In both panels we compare the viscous effects by choosing η/s = 2/(4π) (dashed lines) to the the case of vanishing viscosity, η/s = 0 (solid lines). All lines end at fixed final time τ f = 10 fm/c. Note that we use conventions where µ is the chemical potential for baryons, the chemical potential for quarks is µq = µ/3. vary the shear viscosity to entropy η/s = 2/(4π) (dashed lines) and η/s = 0 (solid lines). All trajectories end at fixed final time τ = 10 fm/c.
For the ideal EOS (left panel of Fig. 2) we observe that the viscosity weakens the effect of the expansion on the temperature T while it does the opposite for the chemical potential µ and thus the trajectories end at larger values of T and smaller vales of µ for non-zero η/s. This is in agreement with the previous discussion of the temporal evolution of T and µ. For the lattice EOS (right panel of Fig. 2) we observe similar trajectories for small initial values of µ(τ 0 ). For larger values of µ(τ 0 ), the trajectories start to bend towards larger values of µ while they continue to decrease towards lower values of T . This behavior is understood from the previous discussion, as well.
In summary, the time-evolution of temperature and chemical potential for a Bjorken expansion is given by Eqs. (27) for an arbitrary EOS. The evolution of µ as a function of time is quite sensitive to the choice of the EOS. The effect of the viscosity is relatively small. This is actually expected for the homogeneous background while we expect more prominent dissipative effects for non-homogeneous perturbations around it. 6 We turn to those in the next section.

IV. FLUCTUATIONS AROUND BJORKEN FLOW
After having studied the solution of the hydrodynamic evolution equations with Bjorken boost invariance and transverse translational symmetries we study now the evolution of fluctuations or deviations from that solution. We will concentrate here on deviations that are small enough in magnitude to describe their evolution by linearized evolution equations. In other words, we write the fluid dynamic fields as whereū µ ,¯ ,n is the Bjorken-type solution discussed in the previous section. The linearized evolution equations for the perturbations δu µ , δ , δn are discussed for a generic background solution and arbitrary coordinate system in Appendix B. If one specializes to the Bjorken background and the coordinate system (τ, r, φ, η), the independent fluid dynamic fields are in the first order formalism δ , δn, δu r , δu φ and δu η . (We take the background fluid velocityū µ and the full fluid velocity u µ =ū µ + δu µ to be normalized, u µ u µ =ū µū µ = −1, such that one has δu τ = 0 at linear order in perturbations). Equation (B3) yields the following equation for the perturbation in energy density (each hydrodynamical fluctuating field depends on (τ, r, φ, η) which we suppress for better readability) The thermodynamic derivatives like (∂p/∂ ) n , etc., are to be evaluated here on the background solution and similarly the transport coefficients and their derivatives. The evolution equation for the perturbation in baryon number density is The derivative operator of second order that appears in the last two lines in front of δ and δn, respectively, is the Laplace operator in the spatial coordinates r, φ and η.
The fluid velocity in the radial direction is determined by the following evolution equation the one in the azimuthal direction by and finally the fluid velocity component in the rapidity direction is governed by Equations (38) -(42) are hyperbolic coupled linear differential equations for the variables δ , δn, δu r , δu φ and δu η : They contain only first order derivatives with respect to the time coordinate τ but up to second order derivatives with respect to the spatial coordinates r, φ and η. In the second order gradient expansion the equations would be elliptical but also contain more degrees of freedom and transport coefficients. In order to analyze the differential equations (38) - (42) it is convenient to use a Bessel-Fourier transformation. For the perturbation in energy density this reads Since δ (τ, r, φ, η) ∈ R and J −m (kr) = (−1) m J m (kr) one has δ * (τ, k, m, q) = (−1) m δ (τ, k, −m, −q).
For the baryon number density fluctuation δn and the rapidity component of the fluid velocity δu η one can use the same expansion. For the fluid velocity components δu r and δu φ we write instead with δu + * (τ, r, φ, η) = δu − (τ, r, φ, η). We expand δu − (τ, r, φ, η) and δu + (τ, r, φ, η) similar to Eq. (43) but replace J m (kr) by J m−1 (kr) and J m+1 (kr), respectively. The reality constraint becomes In terms of the Bessel-Fourier transformed variables one can easily perform the spatial derivatives in Eqs.  The Bessel expansion we use in Eqs. (43) contains an integral over all (positive) values of k. This expansion, also known as the Hankel transformation, is appropriate for functions on the open interval r ∈ (0, ∞). More realistically, the energy distribution in a heavy ion collision is non-zero only on a compact interval (0, R) with some radius R that depends on time during the expansion of the fireball and it is of the order of R ∼ 10 fm. On such a compact interval the Bessel expansion becomes discrete, in the sense that the integral over k is replaced by a sum over a discrete subset. For example, the boundary condition δ = 0 at r = R leads to the values k where ρ(r) is a monotonous function into the interval (0, 1) and a particularly useful choice for ρ(r) is discussed in Appendix A of Ref. [51].
The evolution equation for the perturbation in energy density, Eq. (38) becomes in Bessel-Fourier space (all perturbation functions have now the argument (τ, k, m, q) that we suppress for better readability) Similarly, the evolution equation for the perturbation in baryon number density becomes Let us now turn to the perturbations in the fluid velocity. Equations (40) and (41) lead to the following equations for δu + and δu − in Bessel-Fourier space and for the rapidity component we find from Eq. (42) Note that eqs. (49) -(52) are now coupled ordinary differential equations. All spatial derivatives have become algebraic and one can directly integrate for the time dependent perturbations δ (τ, k, m, q) etc. To construct such a solution one needs as an input the background or Bjorken solution forT (τ ) andμ(τ ) as well as the relations that express all other thermodynamic densities (¯ ,p,n etc.), transport coefficients (ζ,η,κ) and derivatives ((∂p/∂ ) n , (∂p/∂n) , (∂ζ/∂ ) n etc.) in terms of the independent thermodynamic variablesT andμ. Let us first discuss some limiting cases of Eqs. (49) -(52) with extended symmetries.

A. Statistical baryon number conjugation symmetry
If the baryon number density vanishes in the background solution, i. e.n =μ = 0, one has an extended symmetry namely baryon-anti-baryon or baryon number conjugation symmetry corresponding to , n → −n. Odd derivatives such as (∂p/∂n) or (∂η/∂n) have to vanish and one finds that δn decouples from the equations for δ in Eq. (49) and the perturbations of fluid velocity in Eqs. (51) and (52). However, this does not imply that δn has to vanish as well. Locally and event-by-event one may have a non-zero baryon number density. The evolution equation for this perturbation is obtained from Eq. (50) as The second term on the left hand side accounts simply for the dilution due to the longitudinal expansion while the third term is a diffusion term due to heat conductivity. Note thatκ is expected to be singular in the limitn → 0 in such a way that the combination of terms that multiplies (k 2 + q 2 τ 2 ) δn remains finite [79]. Therefore, the diffusion term indeed plays a role for the evolution of perturbations δn.
Equation (53) can be directly integrated and its solution reads as where the integrals depend on the heat conductivity and thermodynamic quantities on the background Bjorken solution. While the integral I 1 is typically dominated by late times τ (for example for the ideal thermodynamic equation of state (14), heat conductivity of the form (22) and Bjorken expansion as in Eq. (31)), the integral I 2 is dominated by early times τ ≈ τ 0 . Moreover, for fast thermalization τ 0 → 0 one has formally I 2 → ∞ such that in reality it might be rather large. Modes with q = 0 are therefore strongly damped by dissipative effects of heat conductivity. The evolution equations for the perturbations in energy density δ and fluid velocity are independent of δn. Their solution has already been discussed in a similar setup in Ref. [46].

B. Exact Bjorken boost symmetry
The evolution equations for perturbations (49) -(52) simplify also in a situation where Bjorken boost invariance is realized as an exact symmetry instead of only on a statistical level. In that case one has δu η = 0 and the perturbations δ , δn etc. vanish except for q = 0. Equation (49) becomes and similarly Eq. (50) becomes One observes that (56) and (57) depend on δu + and δu − only via the combination (δu + − δu − ) / √ 2, for which one obtains from Eq. (51), Equations (56) - (58) together with the information about background quantities form a closed system that describes the analog of sound propagation and baryon number diffusion in the transverse plane of a longitudinally expanding fireball. The orthogonal combination of fluid velocity perturbations δu + + δu − is a shear mode with purely dissipative behavior (equation not shown). It is interesting to compare these equations to the ones that govern perturbations in a static medium. In that case all terms that involve explicit factors 1/τ or derivatives of background quantities with respect to τ vanish. For example, the analog of Eq. (56) is while the analog of Eq. (57) is = 1.30 fm −1 (blue curves). We compare two different values of the ratio of shear viscosity to entropy density η/s = 1/(4π) (left column) and (b) η/s = 10/(4π) (right column). Heat conductivity is related to this by eq. (25). We use T0 = 0.5 GeV, µ0 =0.05 GeV, τ0= 1 fm/c, τ f =10 fm/c and for the initial values of the hydrodynamic fluctuations we choose δ (τ0) = 0, δn(τ0) = δu + (τ0) = δu − (τ0)=0. We denote . See text for further details. and the analog of Eq. (58) is  and have to be integrated numerically for a given background solution and wavenumber k. This has already been discussed in Ref. [46]. In Figs. 4, 5 and 6 we show numerical solutions of the evolution equations (56)-(58) for the ideal EOS (14). For the background fields we employ the scaling solution (30). We compare the numerical results for different initial conditions and two different values of the ratio of shear viscosity to entropy density and assume ζ = 0 for simplicity. More precisely, the left columns of Figs. 4, 5 and 6 correspond to η/s = 1/(4π), the right columns to η/s = 10/(4π). In all cases, the heat conductivity is taken to be related to the shear viscosity by (25). We also compare different values of the radial wavenumber k = k . In all cases, the modes with larger k are damped more quickly as expected. In order to simplify the notation we use the abbreviation ∆ − ≡ δu + − δu − in Figs. 4, 5 and 6.
In Fig. 4 we have chosen initial conditions with non-vanishing perturbations in energy density δ (τ ) = δ 0 while the perturbations in baryon number density δn and fluid velocity δu + , δu − vanish initially. The pressure gradients associated with δ induce sound waves with the typical oscillating behavior between δ and δu + − δu − , modified by the longitudinal expansion and viscous damping. As expected, the oscillation frequency is larger for larger radial wave numbers k. The perturbation in energy density δ induces also a small perturbation in baryon number density δn at times τ > τ 0 . This is due to the linear mixing between the different fluctuating fields (δ , δu + and δu − ) for non-vanishing background baryon chemical potential (we choose µ 0 = 0.05 GeV). Forμ =n = 0, the evolution equation for δn would decouple from the other fluctuating fields as we discussed in the previous section. Because we solve linearized equations for the perturbations, the solution scales linearly with the initial value δ 0 .
In Fig. 5 we initialize with non-vanishing perturbation in the baryon number density δn(τ 0 ) = δn 0 but set δ (τ 0 ) = δu + (τ 0 ) = δu − (τ 0 ) = 0. The mode excited in this way has essentially diffusive behavior. This is most clearly seen in the intermediate panel which shows the temporal evolution of δn/δn 0 . There are no oscillations seen but simply a decay in amplitude which is faster for large values of k. This decay is mainly a consequence of heat conductivity (or equivalently, baryon number diffusion). In addition to the baryon number density perturbation, also a (small) perturbation in δ and δu + − δu − is excited for τ > τ 0 . This is again a consequence of the non-vanishing baryon number density in the background. The behavior of these perturbations is oscillatory, i.e., of sound type.
In Fig. 6 we choose initial conditions with ∆ − 0 = δu + 0 − δu − 0 = 0 while the perturbations δ and δn vanish initially. This results again in sound propagation of the typical oscillating type. In Fig. 6 we also show the behavior of perturbations in the orthogonal combination ∆ + = δu + + δu − which is a shear mode whose decay rate is determined by shear viscosity η. The shear viscosity dependence of the decay rate for this particular shear mode can be obtained directly from the corresponding evolution equation (51).
In Fig. 7 we show the final amplitude of the perturbations in energy density (left column) and particle density (right panel) at τ f = 10 fm/c as a function of the k-wave number in units of the initial weight at time τ 0 =1 fm/c for η/s = 1/(4π) and different initial conditions of the perturbations of the fluctuating fields. This plot shows that some modes of the initial perturbations characterized by the k-wave number indeed survive the entire evolution of the system and at the same time, it also indicates the distribution of of the surviving modes at the scales of time relevant for the freeze-out surface. 7 . The uppermost panel corresponds to the non-zero value for the initial perturbation in energy density δ (τ 0 ) = δ 0 while the remaining fluctuating fields, δn 0 , δu + 0 and δu − 0 , vanish exactly. The middle panel corresponds to the case where δn(τ 0 ) = δn 0 and δ 0 = δu + 0 = δu − 0 = 0. The bottom panel corresponds to ∆ − 0 = δu + 0 − δu − 0 = 0 and δ 0 = δn 0 = 0. For the sound wave type initial conditions (δ 0 = 0 or ∆ − 0 = 0) the size of the amplitudes at τ f present a damped oscillatory behavior while for the case when δn 0 = 0 the fluctuation of the δ and ∆ − present an oscillatory behavior while δn shows a exponential type decay which is typical to diffusive processes. In all the cases we observe that essentially none of the modes survive for values of k ≥ 2 fm −1 .
We conclude this subsection by emphasizing again the observation that perturbations in baryon number density have a diffusive time evolution with a dissipation rate determined by the heat conductivity. For typical values corresponding to strong coupling behavior, the damping is rather strong but at least the modes with the smallest radial and azimuthal wave-numbers (small values of m and l) are not dissipated completely and could have experimentally observables consequences.

C. Exact transverse translation and rotation symmetry
One can also consider a situation with exact symmetry under translations and rotations in the transverse plane. In that case only perturbations with k = 0 are possible and the fluid velocities in transverse directions have to vanish, δu + = δu − = 0. Again Eqs. (49)-(52) simplify substantially, albeit not to a point where they can be integrated directly. Specifically, Eq. (49) becomes and the evolution equation for the perturbation in baryon number density (50) becomes ∂η ∂n iq τ 2 δn + ζ + 4 3η These equations simplify further in a situation with vanishing baryon number density, the numerical solution for this situation has already been discussed in Ref. [46]. The solution of the fluctuating fields is in general complex but subject to the reality constraints δ * (τ, q) = δ (τ, −q) and similar for δn and δu η .
In order to gain some qualitative insights let us consider a simple equation of state = 3p while settingn =μ = 0 and neglecting the effects of viscosities where they are sub-leading compared to other background terms. One can then derive for the variable δ = δ /¯ the equation More general, the set of equations (62), (63) and (64) describe also baryon number density waves and diffusion in the longitudinal direction. We show numerical solutions to the evolution equations (62) -(64) in Figs. 8, 9 and 10 for different initial conditions. As we proceed in Sect. IV B we compare two different values of the ration of shear viscosity to entropy density η/s = 1/(4π) (left panel) and η/s = 10/(4π) (right panel). For the background fields we use again the scaling solution (30).
For Fig. 8 we choose only δ to be non-zero initially. Compared with the behavior of the transverse sound waves or sound waves in a static medium discussed in the previous section, the evolution of the resulting longitudinal sound waves is completely different. In particular, no proper oscillations are visible during the entire temporal evolution. Rather, one observes a decay in amplitude, in particular at early times. This effect of the longitudinal expansion is particularly strong for large values of q. At later times the damping actually weakens to the extent that amplitudes remain non-zero at the final time. Interestingly, the influence of viscosity on the time evolution of longitudinal perturbations is relatively weak. Some quantitative differences are of course visible between the left and right panels of Fig. 8 but qualitatively, the evolution is surprisingly similar.  In Figure 10 we choose δu η 0 = 0. As in the previous two situations, there is no proper oscillation visible for the time interval shown.
Finally, Fig. 11 shows the final amplitude of perturbations at τ f = 10 fm/c in units of the initial weight at time τ 0 as a function of the longitudinal q-wavenumber. Figure 11 is obtained by choosing a non-vanishing value initially for one particular fluctuating field while the remaining fluctuating fields are initially set to zero. The top, middle and bottom panels of Fig. 11 corresponds to δ 0 = 0, δn 0 = 0 and δu − 0 = 0, respectively. We observe that the amplitude of the fluctuating modes goes asymptotically to zero for q ≥ 25 which corresponds to a small window in the rapidity variable (i.e., ∆η ∼ (∆q) −1 ). We expect that modes with intermediate and large q would be damped stronger for earlier initialization time τ 0 .
Finally, in a situation where Bjorken boost invariance as well as translations and rotations in the transverse plane are realized exactly, i. e., δu η = δu + = δu − = k = q = 0, Eqs. (49) and (50) reduce simply to a linearized version of the Bjorken expansion in Eq. (26) as it has to be.

V. THE TWO POINT CORRELATION FUNCTION OF BARYONIC PARTICLES
In this section we discuss the possibility to access the information about perturbations in the baryon number density experimentally by measuring a correlation function of the net number of baryons (baryons minus anti-baryons) as a function of the rapidity and azimuthal angle. We concentrate for simplicity on the case of vanishing background baryon number density as discussed in Sect. IV A.
Perturbations in baryon number density in position space as described by Eq. (54) are not directly accessible to experiments. However, a fluctuating baryon number density and chemical potential on the kinetic freeze-out surface has an influence on the distribution of particles with non-zero baryon number in momentum space. This concerns in particular protons but also resonances with non-vanishing baryon number. Similar as for flow observables, there is a direct link between different harmonics in azimuthal angle and rapidity in the fluid dynamic description and the corresponding harmonics in the momentum space particle distribution. Thus, we can partly access the physical information contained in Eq. (54). As an example, we consider a connected two-point correlation function of the type 8 which measures correlations of baryonic particles (i.e. the number of baryons minus anti-baryons) as a function of the difference between (particle momentum) azimuthal angles φ 1 − φ 2 and (particle momentum) rapidities η 1 − η 2 . In Eq. (66), n Baryons (φ, η) is the number of baryons minus anti-baryons as found in the detector in a particular bin in azimuthal angle φ and rapidity η 9 . We also introduce the Fourier representation The correlation function in Eqs. (66) and (67) is determined by a combination of initial conditions (set at the time where a fluid dynamic description becomes valid) and response functions that describe how baryon number density perturbations propagate in the fluid dynamic regime and how they influence the particle distributions at freeze-out. In the following we discuss both parts in a bit more detail. First, the initial state after a heavy ion collision (and after the early non-equilibrium dynamics) at the time τ 0 when a fluid dynamic description becomes valid, is characterized by a fluctuating baryon number density δn(τ 0 , r, φ, η) around some average or expectation valuen(τ 0 , r). (The latter might be rather small at LHC and upper RHIC energies and we neglect it in the following.) For the fluctuating part we use a Bessel-Fourier decomposition An event-by-event ensemble of initial conditions conditions for the baryon number density can be characterized in terms of the weights δn (m) l (q). For example, the two-mode correlation function is δnδn;l1,l2 (q).
We have assumed here that the ensemble of initial conditions is statistically symmetric under azimuthal rotations and longitudinal boosts leading to the factors δ m1+m2,0 and 2πδ(q 1 + q 2 ) on the right hand side of Eq. (69). For a single event with baryon number perturbation as in (68), the baryon number distribution in momentum space after kinetic freeze-out will be proportional to the weights δn   8 The brackets · · · in Eq. (66) denote an event average O(x, y) = lim y) . 9 There is a complication due to the fact that neutrons cannot be measured experimentally. Further studies are needed in order to quantify whether this presents a problem for observables as in (66) and if so, how these can be overcome. Also, one should estimate possible contributions to Eq. (66) from sources other than fluid dynamics, such as resonance decays.
with linear baryon number response function S Baryons;(m)l (q). The object on the left hand side of Eq. (70) is the Bessel-Fourier weight of the (momentum space) distribution of the number of baryons minus anti-baryons. The correlation function on the right hand side of (67) can be written as δnδn;l1,l2 (q).
For a more detailed discussion of the response function formalism briefly introduced here we refer to Ref. [51]. The linear response functions S Baryon;(m)l (q) are in particular also affected by heat conductivity. More specific, the analog of the factor exp(−k 2 I 1 − q 2 I 2 ) in a situation with realistic transverse dependence and radial flow leads to a suppression of modes with q 2 > 0 and large values of m and/or the radial wave number l. Qualitatively, one expects that the scale for the suppression in the transverse direction is set by the (time dependent) radius R of the fireball. Moreover, the l'th zero crossings z where on the right hand sideCκ =0 Baryon (m, q) would be the corresponding correlation function in the (somewhat hypothetical) situation of vanishing heat conductivity and the dissipative attenuation terms can be roughly estimated as Going now back to the two-particle correlation function (66), the exponential suppression factor in Eq. (72) implies for large I 2 long range correlations with respect to the rapidity difference η 1 − η 2 , with a decay that is determined by the value of I 2 (except ifCκ =0 Baryon (m, q) has a very strong decay with q already) and a similar, although weaker, effect with respect to the azimuthal wavenumber m. In order to make our qualitative statements more precise, it is necessary to generalize the calculations described here to a more realistic background. A more realistic background would have a realistic transverse profile and expansion in addition to the longitudinal (boost-invariant) expansion. Moreover, one also has to perform more detailed studies of the initial conditions and kinetic freeze-out, that both affectCκ =0 Baryon (m, q).

VI. CONCLUSIONS
We have studied solutions of the fluid equations describing relativistic heavy ion collisions in the presence of a globally conserved quantum number (baryon number) using a background-fluctuating splitting. For the background we have assumed Bjorken boost and transverse translation and rotation invariance. This generalizes Bjorken's original solution to non-vanishing baryon number density as well as shear and bulk viscosities. Heat conductivity does not play a role on the background equations since the diffusion current vanishes exactly due to the symmetries of the Bjorken flow.
We derived evolution equations for the perturbations around this background solution. While the amplitude of these perturbations was assumed to be small, such that linearized equations could be used, the formalism allows us to treat perturbations with arbitrary dependence on the transverse coordinates and rapidity. Technically, this is done by employing a Bessel-Fourier expansion. The partial differential equations of relativistic fluid dynamics become ordinary differential equations for the different modes that are characterized by radial, azimuthal and rapidity wave numbers. The evolution of these perturbations is governed by the thermodynamic properties encoded in the equation of state p(T, µ) as well as the transport properties (i.e., shear viscosity η(T, µ), bulk viscosity ζ(T, µ) and heat conductivity κ(T, µ) in the first order formalism we use).
Generically, one finds that perturbations with large wave numbers are damped more quickly by the dissipative processes, as expected. The dissipation of different modes depends on time in a different way and, in particular, deviations from Bjorken boost symmetry show a fast damping at early times. In principle, it might be possible to use these dependencies to probe transport and thermodynamic properties at different times in the evolution history and therefore for different temperatures of the quark-gluon plasma produced in a heavy ion collision.
In order to make more quantitative statements, one must take a realistic transverse density profile and expansion into account, of course. This has been done for perturbations with exact Bjorken boost symmetry and vanishing baryon number in Refs. [47,48,50,51]. In the present paper we have concentrated mainly on the evolution of perturbations in baryon number density. They have diffusion-type evolution governed by the longitudinal expansion and heat conductivity. (In the Landau frame, heat conductivity can in fact be understood as baryon number diffusion.) There are characteristic differences in the dependencies on longitudinal and transverse wave numbers. More specific, baryon number perturbations are quickly "flattened out" in the longitudinal direction at early times.
In principle, the information on baryon number perturbations is accessible experimentally via two-point (and higher order) correlation functions of particles with non-zero baryon number, as a function of the difference in azimuthal angles and rapidities. Based on the evolution equations for perturbations, we expect long-range correlations in rapidity (a "baryon number ridge"). For a more detailed theoretical picture one needs a better description of the local event-byevent fluctuations in baryon number density at the initial time when fluid dynamics becomes valid. Also, one should take a realistic transverse expansion into account and study the implications of baryon number perturbations at the kinetic freeze-out. (Formulas needed for this have already been derived in Ref. [50].) It would be very interesting to study net-baryon number correlations experimentally, as well as theoretically in more detail, and thereby constrain heat conductivity as another property of the quark-gluon plasma. In this appendix we compile some thermodynamic relations in the grand canonical ensemble that we found useful in the context of relativistic fluid dynamics with a conserved charge. We start from the pressure p(T, µ), which is related to the thermodynamic potential of the grand canonical ensemble (the Landau potential) by p = −Ω/V . The differential of pressure is dp = sdT + ndµ . (A1) All thermodynamic quantities can be obtained from this and the Gibbs-Duhem relation + p = T s + µn, for example In the following we will sometimes drop the subscripts with the convention that pressure is evaluated as a function of T and µ unless indicated otherwise. Also we find it useful to express all susceptibilities in terms of the pressure and its derivatives. This avoids ambiguities and realizes Maxwell's relations automatically. For example, the energy density is obtained then as = −p + T ∂p ∂T + µ ∂p ∂µ . (A3) Its differential, as well as the one for density, are d = T ∂ 2 p ∂T 2 + µ ∂ 2 p ∂T ∂µ dT + T ∂ 2 p ∂T ∂µ + µ ∂ 2 p ∂µ 2 dµ , and a modified sound velocity at fixed particle densitỹ Both sound velocities agree for vanishing baryon number density, n = 0. Note that the usual relations where, similarly to pressure p(T, µ), the bulk viscosity ζ(T, µ) and shear viscosity η(T, µ) are functions of T and µ on the right hand side.