Stochastic hydrodynamics and long time tails of an expanding conformal charged fluid

We investigate the impact of hydrodynamic fluctuations on correlation functions in a scale invariant fluid with a conserved $U(1)$ charge. The kinetic equations for the two-point functions of pressure, momentum and heat energy densities are derived within the framework of stochastic hydrodynamics. The leading non-analytic contributions to the energy-momentum tensor as well as the $U(1)$ current are determined from the solutions to these kinetic equations. In the case of a static homogeneous background we show that the long time tails obtained from hydro-kinetic equations reproduce the one-loop results derived from statistical field theory. We use these results to establish bounds on transport coefficients. We generalize the stochastic equation to a background flow undergoing Bjorken expansion. We compute the leading fractional power $\mathcal{O}((\tau T)^{-3/2})$ correction to the $U(1)$ current and compare with the first order gradient term.


I. INTRODUCTION
The hydrodynamic description of relativistic heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) has been extremely successful [1][2][3]. State-of-the-art models contain not only higher order gradient terms in relativistic fluid dynamics, but also final state kinetic theory afterburners, and initial state models that account for event-by-event fluctuations. There is strong evidence for the role of initial state fluctuations from the observation of odd azimuthal moments of hydrodynamic flow [4].
The role of thermodynamic fluctuations during the evolution of the system is much less explored. These fluctuations arise from the fact that fluid dynamics is a coarse-grained description, so that at any coarse graining scale unresolved microscopic degrees of freedom lead to fluctuations of the macroscopic variables. These effects are interesting, because hydrodynamics is a nonlinear theory, so that couplings between hydrodynamic modes lead to novel phenomena beyond just Gaussian fluctuations in the hydrodynamic variables. For example, it is known that fluctuations lead to hydrodynamic "tails", which are non-analytic contributions to the time (or frequency) dependence of correlation functions [5][6][7]. In the gradient expansion these terms are formally more important than corrections beyond the Navier-Stokes approximation. Hydrodynamic tails also imply non-trivial bounds on the transport coefficients that are completely independent of ideas about holographic dualities.
More importantly, hydrodynamic fluctuations dominate the dynamics near a critical point, and any fully dynamical description of fluctuations of conserved charges in a heavy ion collision must include fluctuations if the hydrodynamic evolution equations. This implies that the study of hydrodynamic fluctuations is crucial for interpreting the results from the RHIC beam energy scan program.
There are a variety of techniques for taking int ac-count hydrodynamic fluctuations. In homogeneous systems fluctuations have been studied using diagrammatic or effective action methods [5,8]. Combined with the epsilon expansion, these methods have led to the determination of dynamical critical exponents [9]. In a complicated environment like a heavy ion collision these methods are less convenient. An alternative is to numerically study stochastic hydrodynamic equations. There is some progress in this effort [10], but the calculations are difficult, and the results are not always simple to interpret. Several authors have derived deterministic equations for the time evolution of hydrodynamic n-point functions [11][12][13]. These methods rely on linearizing the fluid dynamic equations, and on truncations in the number of moments, but the resulting equations are easier to evolve in time, and the comparison to analytical results for homogeneous systems is more straightforward.
In this work we adopt the latter approach, extending a formalism developed by Akamatsu et al. [11] to two-point functions involving a conserved number density. This is an important problem, because the conserved baryon density (or the entropy per baryon number) is believed to be the order parameter for a possible critical point in the QCD phase diagram. In the present work, we will restrict ourselves to the two-point function of a scale invariant fluid, leaving the role of a critical equation of state to a future study. Estimates of the role of critical fluctuations have recently appeared in [14]. As indicated above we will linearize the equations of motion, an expand in powers of fluctuating fields and external sources. We provide bounds on the diffusion constant, and derive the tails of diffusive correlators in a Bjorken backgound.
We use the following conventions: We use a bar to denote background values of the hydrodynamic variables. Four dimensional indices are denoted by Greek letters µ, ν, ... = 0, 1, 2, 3; three dimensional indices are labeled by Latin letters i, j, ... = 1, 2, 3. Three dimensional spatial vectors are denoted by bold letters (x,k). The signature of the metric is taken to be 'mostly-plus'. A four dimensional vector is denoted by x µ = (x 0 , x). The fluid velocity is a time-like vector with u µ u µ = −1. We denote the spatial measure of the integral d 3 x √ −g ≡ x while the three dimensional measure in momentum space is d 3 k/[(2π) 3 √ −g] ≡ k .

II. RELATIVISTIC HYDRODYNAMIC FLUCTUATIONS OF A CONFORMAL CHARGED FLUID
In this section we describe the theoretical framework for studying the evolution of hydrodynamic fluctuations on top of an evolving background for a charged fluid. We consider a relativistic fluid with an unbroken global U (1) symmetry (such as baryon number) where parity is conserved. The local conservation laws for this system are [15,16] where T µν is the energy-momentum tensor of the fluid, J µ is the U (1) current and F µν = ∂ [µ A ν] is the field strength tensor of a gauge field A µ coupled to the conserved current. The energy momentum tensor couples to the metric g µν , and the covariant derivative associated with g µν is denoted by D µ . In the Landau frame where u µ T µν = ǫu ν (being u µ the time-like vector denoting the fluid velocity and ǫ the energy density) the energy-momentum tensor T µν and the current J µ can be decomposed as In the previous equations the macroscopic quantities of the fluid are the isotropic pressure p, the viscous stress tensor π µν , the particle density n and the diffusive current v µ . Fluctuating contributions to the energymomentum tensor and current are denoted by S µν and I µ . We introduce the tensor ∆ µν = g µν + u µ u ν which is orthogonal to the fluid velocity. In order to solve the conservation laws of stochastic hydrodynamics (1) one requires an equation of state and the constitutive relations that express the conserved currents in terms of hydrodynamic variables. Within the Navier-Stokes approximation and for conformal systems the constitutive relations for the stress tensor π µν and the particle diffusion current v µ read as respectively where σ µν is the Navier-Stokes shear tensor and E µ is an external electric field. σ µν and E µ are defined as follows where W µν = ∆ µν αβ A αβ with Note that ∆ µν αβ is a 2-rank tensor which is traceless, symmetric, and orthogonal to the fluid velocity u µ . The bare transport coefficients associated with charge conductivity and shear viscosity are denoted by σ 0 and η 0 respectively. In Eqs. (3) we introduce the projector orthogonal to the fluid velocity ∆ µν = g µν + u µ u ν .
The set of equations (1) are not complete unless we specify the equation of state as well as the functional form of the random sources S µν and I µ . At leading order in the low momentum (gradient) expansion we can model the random sources as delta-correlated white noise [17][18][19], i.e., These approximations will allow us to recast the problem of solving the equations of stochastic hydrodynamics onto a set of coupled Langevin equations. In this work we consider a conformally invariant fluid whose equation of state (EOS) is where g is an arbitrary function. For this particular equation of state, the conformal symmetry implies the ratio µ/T is invariant under Weyl rescaling [20]. Different thermodynamic relations associated to the EOS (7) are listed in App. A.

III. HYDRODYNAMIC FLUCTUATIONS AROUND A STATIC BACKGROUND
Hydrodynamic tails appear in the long time behavior of response functions. We compute these response functions by coupling the hydrodynamic evolution in the presence of noise terms to external fields, and then determine the variational derivative of currents with respect to the sources. Specifically, the metric can be used to study response functions of the stress tensor, and the external gauge field determines the current response.
In order to study diffusive terms in the shear response in a static background it is sufficient to consider a metric of the type g µν = η µν +h µν where the metric perturbation is parametrized as where h s is a time-dependent scalar function. The current response is studied by considering an external gauge field of the form In a static background the hydrodynamic equations for the background hydrodynamic are The solutions of these two equations are given in terms of static background configurations, i.e.,ǭ = ǫ 0 andn = n 0 respectively. In the following we will also make use of the heat energy density [21,22] q =ǭ −wn ,w = (ǭ +p)/n , wherew is the enthalpy per particle. From the first law of thermodynamics the variation of the heat energy density represents the change of the entropy per particle 1 δq ≡ Tn δ s n .
Given the conservation of the entropy (or the heat energy density) one can define the relation between the pressure and the energy density as where we recognize v 2 a as the (adiabatic) speed of sound. For the EOS (7) β 1 = 1/3, β 2 = 0 so v 2 a = 1/3 as expected (see App. A).

A. The Navier Stokes-Langevin equations
In this section we will derive the evolution equation for hydrodynamic fluctuations in the presence of noise, the linearized Navier-Stokes Langevin equation. We begin by linearizing the equations of motion of relativistic fluid dynamics around a static background. We write the hydrodynamic variables as ǫ =ǭ + δǫ , p =p + δp , n =n + δn , where δǫ, δp and δn are perturbations in the energy density, pressure and particle density, respectively. We will write the equations of motion in a mixed representation where we perform a Fourier transform with respect to 1 Eq. (12) considers that the total number of particles N is fix and thus, all the thermodynamic derivatives must be taken at fixed N [21,22]. spatial variables 2 . The linearized stochastic hydrodynamic equations are Here, we have defined the electric field by E i (x 0 , k) = F i0ū 0 = F 0i , and we denote k 2 = |k| 2 . We have also introduced the momentum diffusion constant whereH is the enthalpy density. In Eq. (16c) we did not keep higher order terms in frequency and wave number. We have also defined the thermodynamic quantities The natural degrees of freedom of Eqs. (16) are the variations of the energy density δǫ, particle density δn and fluid velocity δu i . However we are free to use any other set of variables when analyzing Eqs. (16). It is more convenient to use the heat energy density δq, pressure δp/v a and momentum flux density g i ≡ T 0i . These variables have the property that in thermal equilibrium their fluctuations are statistically independent [17,21,22]. For the conformal EOS (7) and close to the thermal equilibrium we have ϕ a = M ab ψ b with ϕ = (δǫ, δu i , δn) and ψ = (δp/v a , g i , δq) and the matrix M ab given by The momentum density is a three-dimensional vector. It is convenient to express this vector in a basis spanned by the vectorsk = k/k andê Ti (i=1,2) defined by [11] k = (sin θ cos φ, sin θ sin φ, cos θ) , (20a) e T1 = (− sin φ, cos φ, 0) , (20b) e T2 = (cos θ cos φ, cos θ sin φ, − sin θ) .
(20c) 2 The spatial Fourier transform is defined as In this basis the momentum flux density g i is decomposed as Using (19) and (16) we obtain the Navier-Stokes-Langevin (NSL) equations in the form where is a five dimensional vector representing the fluctuating hydrodynamic fields. The acoustic and diffusive matrices, denoted by A and D respectively, are where we have defined the diffusion coefficient D 0 = σ 0 α 2 . The NSL equation describes five hydrodynamic modes, and the matrices A and D determine the dispersion relation of these modes at O(k) and O(k 2 ) in a expansion the wave number k. At leading order (LO) the acoustic matrix A describes two sound modes with eigenvalues λ ± = ± i v a , one diffusive (d) mode and two shear modes (T 1 , T 2 ) [5,[21][22][23]. At this order in the expansion the eigenvalues of the diffusive and shear modes are λ A = 0 (A = d, T 1 , T 2 ). The absence of dissipation at LO implies that δq is a conserved quantity [22].
At next-to-leading order (NLO) non-diagonal terms in the matrix D lead to mixing between the pressure, the heat energy density, and the longitudinal momentum density. Both the momentum density and the heat energy receive real (diffusive) corrections O(k 2 ). As a result the dynamics of the modes (T 1 , T 2 ) remains purely diffusive, while (δp/v a , g L , δq) exhibit a combination of diffusive and reactive behaviour.
The source matrix P ab in Eq. (22) takes into account the non-linear couplings between the fluctuating hydrodynamic fields and the external sources. Its explicit components are Finally the vector ξ a carries the information of the random noise sources. Its components are The solutions of the NSL equations (22) determine the linear evolution of variations of the pressure, heat energy, and momentum flux densities in response to noise and external sources.

IV. EVOLUTION EQUATIONS FOR TWO-POINT CORRELATION FUNCTIONS
The general solution of the NSL equations (22) can be written as where ψ b (0) is the initial condition for the hydrodynamic fluctuating fields at the initial time x 0 = 0 and the matrix operator This operator is a solution of the differential equation The discussion of the mathematical properties of the evolution operators (28) and their application to stochastic processes in statistical mechanics can be found in Refs. [24,25]. If one performs a stochastic average over the solution (27) in the absence of external sources then one finds that ψ a (x 0 , k) = 0. If instead one considers an ensemble average over different initial conditions then ψ a (x 0 , k) = U ab (x 0 ) ψ b (0) describes the hydrodynamic response to initial state fluctuations. This is an interesting problem in its own right [26][27][28][29][30][31][32][33][34][35], but not the main focus of the present study.
Information on hydrodynamic fluctuations is contained in two-point and higher n-point correlation functions. We define the symmetric two point correlation function at equal times as 3 The reality constraint of the hydrodynamic fluctuating fields implies φ * (x 0 , k) = φ(x 0 , −k) so that C ab is a real symmetric matrix. Using the equation of motion for the hydrodynamic variables ψ a , Eq. (22), we can derive an evolution equation for the two-point function. We find 3 In an effective hydrodynamic theory where the fields are purely classical there is no distinction between the symmetrized correlator 1 2 {B(t), A(0)} and the Wightman correlator B(t)A(0) [36].
where N denotes the two-point correlation function of the noise The noise correlator was originally defined in Eq. (31).
Here, we have expressed N in terms of the diffusive matrix D ab and the susceptibility matrix χ ab . We define χ ab in Appendix A.
The equation of motion for the correlation functions C ab mixes different channels. Following the approach of Akamatsu et. al. [11] we write the ODE's (30) in a basis spanned by the eigenvectors of the acoustic matrix A. The acoustic matrix A (24a) is a 5 × 5 antihermitian matrix (i.e. A † = −A) and the solutions to its eigenproblem The eigenvalues of A are purely imaginary or zero. In our case these are given by In general, the eigenvectors of a non-hermitian matrix need not be mutually orthogonal, and it is necessary to construct a set of left and right eigenvectors. Nonetheless, the antihermitian property of A ensures that the eigenvectors are orthogonal. The eigenvectors φ i of A (24a) are given by wherek andê T l are given by Eqs. (20). In the eigenmode basis Eqs. (30) take the form where are the eigenvalues of the acoustic matrix A. In components the equations of mo-tion can be written as with C 0 =TH. In a recent work Akamatsu et. al. [14] studied the dynamics of Eqs. (35a)-(35d) for a generic EOS. For the purpose of determining hydrodynamic tails we are interested in the long-time, low-frequency behavior of the correlation functions. In the long-time limit diagonal correlation functions relax to equilibrium states fixed by fluctuation-dissipation relations. External sources lead to small deviations that can be studied perturbatively. The solutions of Eqs. (35a)-(35b) are given by where The correlation function C dd does not directly couple to a source and is given by Off-diagonal correlation functions relax to zero, but external sources, coupled to the equilibrium behavior of C T l T l and C dd , can drive them off equilibrium. The asymptotic behavior of C dT l is where we have defined χ dT l = χ T l T l + χ dd .

V. MODIFICATION OF THE RESPONSE FUNCTIONS DUE TO HYDRODYNAMIC FLUCTUATIONS
The response functions of the energy-momentum tensor and the charge current are defined by For conformal systems the response functions G µν,λσ R and G µν R can be decomposed as [37][38][39] G µν,λσ where G R S and G R J are scalar functions while P µν and H µν,λσ (for d=4) are given by The operator H µν,λσ is a projector onto conserved traceless symmetric tensors which satisfies η µν H µν,λσ = 0 [37][38][39]. In the limit k → 0 the scalar functions G B , G S and G J can be written in terms of δT ij and δJ i as follows In hydrodynamics we can read off the variation of the spatial components of the energy-momentum tensor (2a) and the charge current (2b) from the constitutive relations. Taking into account both the effects of external fields and the variations of the hydrodynamic variable we obtain where we include corrections up to second order in fluctuations. These terms describe the non-linear modecoupling between the currents and the hydrodynamic variables. Note that the second order terms are determined by the symmetric two-point functions (29) studied in the previous section [5,11,25,36,40] Here, we have expressed the symmetric correlation matrix C ab (29) in terms of the eigenmodes of the acoustic matrix A. This provides a transparent interpretation of different contribution to the currents in terms of hydrodynamic modes [6,7,11,17,21,22,36,41].

VI. LONG TIME TAILS: STATIC CASE
As a first application of this formalism we compute the low frequency behavior of the response in the static case. For harmonic perturbations h s , δA ∼ e iωx 0 Eqs. (44) give In Eq. (46) the non-linear contributions of the hydrodynamic fluctuations are encoded in the correlators g i g j , (δp/v a ) g i and δq g i . In terms of the eigenmodes of the acoustic matrix (33) these correlators are given by Here, we have neglected terms that correspond to rapidly oscillating modes [11,42]. In the language of the diagrammatic approach these are diagrams that have poles far away from the real axis in frequency space [6,7,36,41]. We observe that the leading term in the stress tensor response, Eq. (46a) and (47a), only contains diagonal correlation functions [11,42]. The current response, on the other hand, is sensitive to off-diagonal correlation functions, see Eq. (46b) and (47c). This is interesting, because in the diagrammatic approach the current response function is factorized into diagonal propagators, see Kovtun and Yaffe [36] as well as App. B. We will show that the result (for the modes taken into account in both works) nevertheless agrees.
We have all the ingredients to calculate the response functions (43). The stress tensor response G R S (43a) was obtained by Akamatsu et. al. [11] and we simply state their result where we have used v 2 a = 1/3 as well as χ AA = 2 for A = ± and χ AA = 1 for A = T l . In Eq. (48) we have regularized UV divergent integrals using a sharp cutoff Λ in momentum space. There are two divergent terms, a O(Λ 3 ) term that can be viewed as a correction to the pressure, and a O(Λ) contribution that modifies the shear viscosity. Stochastic hydrodynamics is renormalizable in the sense that divergent terms can be absorbed in the parameters of the theory, the equation of state and the transport coefficients. The main result is the leading nonanalytic term O(ω 3/2 ), which is independent of the cutoff. This term corresponds to a t −3/2 decay of the real time correlation function and is known as a hydrodynamic tail. The tail of the relativistic stress tensor correlation function was computed in [36,41], and checked using the formalism employed in this work in [11].
Using the same methods we can determine the current response function G R J (43b). We find where we have approximated χ dT l = T c p /H, where c p is the specific heat at constant pressure, to facilitate the comparison with diagrammatic calculations, see App. B. Also see App. A for the definition of thermodynamic quantities. We note that the current response does not contain a cubic divergence. There is a linear divergence which can be absorbed into the conductivity. We find a non-analytic contribution, which is again of the form ω 3/2 . It receives contributions from sound modes, proportional to γ −3/2 η0 , and from the mixing of diffusive heat and shear waves, governed by (D + γ η0 ) −3/2 . Note that in a conformal fluid sound attenuation is purely governed by γ η0 . The purely diffusive terms were previously com-puted by Kovtun and Yaffe [36] 4 . We provide a more detailed comparison with their work in App. B.

Spectral functions
The imaginary part of the response function determines the spectral function, and it can be used to define frequency dependent diffusion constants. We find where the renormalized momentum diffusion constant γ ren. η and renormalized charge diffusion constant D ren. are given by .
We note that the spectral function has a ω 1/2 nonanalyticity near ω = 0, and that the spectral functions are enhanced in the small frequency limit. We also note that the negative sign of the ω 1/2 terms implies that the zero-frequency diffusion constants cannot be arbitrarily small.

Lower bounds of the transport coefficients
The renormalized transport coefficients, Eqs. (51), are given as the sum of a bare contribution and a loop correction which scales inversely with the bare coupling. This implies that the renormalized transport coefficients cannot be arbitrarily small. The physical mechanism underlying this bound is the fact that even if the microscopic diffusion constant is small there is always a macroscopic mechanism, associated with the propagation of hydrodynamic modes, that leads to diffusion.
More formally, we can view this bound as arising from a matching procedure. Consider a UV scale Λ at which the long wavelenght hydrodynamic theory is matched to a microscopic theory with transport coefficients γ η0 and D 0 . Consistency requires that both γ η0 and D 0 are positive. Contributions to the transport coefficients from scales below Λ can be computed in stochastic hydrodynamics, and lead to terms that scale with a positive power of the cutoff, and inversely with the bare transport coefficient. This means that for any value of the UV scale there is a bound on the transport coefficient, and that the bound is stronger the larger the cutoff can be chosen.
The largest possible cutoff scale is determined by the regime of validity of fluid dynamics. Consider a diffusive shear mode. The decay rate is proportional to ω ∼ γ η k 2 . For hydrodynamics to be valid this rate has to be smaller than the next order gradient correction, the shear relaxation rate τ −1 π . This means γ η k 2 < τ −1 π and Λ = (γ η τ π ) −1/2 . In microscopic theories, both in the weak and strong coupling limit, the diffusion constant and the relaxation time are driven by similar physical processes, and τ π = γ η λ(p(T, µ)) [20,[43][44][45][46][47][48], where λ is a function of the thermodynamic pressure p. We conclude that Similar estimates can be obtained by considering the sound and heat diffusion channels [5,7], and we will assume that Eq. (52) can be applied in all channels. For Λ = Λ max. the physical values of the sound attenuation length and the diffusion coefficient, Eqs. (51) read as By extremizing the first of these expressions respect to γ η we get the following lower bound for γ η γ η ≥ γ min.
This result was found previously by Kovtun et. al. [41] for vanishing chemical potential so λ ≡ λ(p(T )) in Eq. (4.3) of Ref. [41]. Now the lower bound of the diffusion coefficient is obtained by considering the bare diffusion coefficient to vanish exactly D 0 = 0 in Eq. (53b) such that Therefore, the lower bound of the diffusion coefficient depend exclusively on γ η0 . The lower bounds of the attenuation transport coefficients also depend on the choice of the equation of state.

VII. LONG TIME TAILS: BJORKEN EXPANSION CASE
We now come to the central problem studied in this work, extending the calculation of diffusive tails to an expanding background. We will consider a conformal U(1) charged fluid undergoing longitudinal boost invariant expansion, known as Bjorken flow. The symmetries of the Bjorken flow are manifest in Milne coordinates where τ = (x 0 ) 2 − (x 3 ) 2 is the proper time, ς = arctanh(x 3 /x 0 ) is spatial rapidity, and the metric is g µν = diag. (−1, 1, 1, τ 2 ).
In the absence of external sources and noise the conservation laws provide the evolution equations for the background hydrodynamic fields. Bjorken flow corresponds to the well known equations where we have neglected dissipative terms. For the conformal EOS (7) the hydrodynamic background fields evolve according tō We also find the background pressurep(τ ) = v 2 aǭ (τ ) as well as the heat energy densityq =ǭ(τ ) −w(τ )n(τ ).
The evolution equations for the hydrodynamic variables in the presence of external fields and noise can be obtained using the procedure discussed in Sect. III A. We again adopt a mixed representation and write define spatial momenta k ≡ (k, k ς ). The spatial Fourier transform is defined as with κ = k ς . The linearized stochastic hydrodynamic equations are where the electric field is E i = F iτ u τ ≡ F τ i . We adopt the basis and write the Navier-Stokes-Langevin equations in the matrix form 5 (62) 5 Here we assume that the background is slowly evolving, so that ∂τ δq = ∂τ δǫ −w∂τ δn.
Here, the acoustic and diffusive matrices are and the source matrix and noise vector are given by where E = (E ⊥ , τ E ς ), 1 2 is the identity in two dimensions, κ = τ 2 k ς and we define Both − → K andK are time depending vectors [11]. Moreover, the ISO(2) symmetry of the background hydrodynamic fields forbids the azimuthal angle ϕ K to be time dependent.
The acoustic matrix has the same structure as its static counterpart (24a), and the set of orthonormal eigenvectors is given by [11] where we have defined The evolution equations for the symmetric correlation functions are determined by following the procedure discussed in Sect. IV. In the following we will make use of the following components of C AB withC 0 = T (τ )H(τ ). We have also used the noise correlator The factor 1/τ stems from the Jacobian transformation of the space-time coordinates in Eq. (6). In order to find the asymptotic solutions we have to specify a functional form for the gauge field A µ . Since we are interested in the propagation of fluctuations that explicitly break the Bjorken symmetry, it is enough to consider an electric field that is a constant vector, i.e. the temporal component of the gauge field is a linear function that grows linearly with the distance A µ = (E · x, 0, 0, 0). The equations of motion for the Bjorken flow have a similar mathematical structure as their static counterparts (35) so one can also obtain time-dependent asymptotic solutions. The diagonal components of the symmetric matrix C AB admit the following asymptotic solution where We proceed as in Sect. IV and obtain the asymptotic solution of the non-diagonal component C dT l In the remainder of this section we outline the calculation of long time tails of the response functions in a Bjorken background. Some of the details are relegated to App. C.
In the static case we observed that the energy-momentum tails (48) are not affected by the chemical potential in the case of a conformal equation of state (7). This is also the case for a Bjorken expansion, and we shall not repeat the calculation of Akamatsu et al. [11]. Instead, we will focus on stochastic corrections to the particle current.
For the Bjorken flow the particle current in stochastic fluid dynamics is where we have kept terms up to second order in fluctuating variables, and we have used that δu τ = (δu i ) 2 /2. The correlation function φ a φ b can be written in terms of the eigenstates of the acoustic matrix as (74) Using the asymptotic solutions for the correlation func-tion (70)-(72) we obtain the currents In Eq. (75a) we recognized that the quadratic perturbations of the momentum flux are uniquely related with δT τ τ (see Eq. (67a) in Ref. [11]). The terms ∆J µ are the finite residual contributions to the particle current which are cutoff independent. The ultraviolet divergences are absorbed into the renormalized macroscopic hydrodynamic variables and parameters, i.e., where the previous expressions are determined by comparing orders with the same power law τ dependence. The renormalized diffusion coefficient matches exactly the static homogeneous one (51b). This also holds for the renormalized sound attenuation lenght γ η0 [11]. The second terms in the RHS of Eqs. (76) correspond to the bare and cutoff dependent terms. The cutoff dependent pieces, O(Λ) and O(Λ 3 ) respectively, are the contributions of the hydrodynamic fluctuating fields below the cutoff Λ. The cutoff dependent terms can be obtained directly from the asymptotic solutions of the symmetric matrix C AB , Eqs. (70)-(72). The bare quantities captures the macroscopic thermodynamic equilibrium properties of the fastest modes which have been already integrated out above the cutoff scale Λ.

A. Finite residual contributions of the particle current
In this section we determine finite cutoff-independent contributions to the particle current. We follow closely the procedure outlined in [11]. We compute the general solution for the correlation function C AB and focus on those terms that are independent of the initial conditions. We remove the leading divergence identified in the previous section, and compute the remaining finite terms. The technical details of this calculation are presented in appendix C.
From the formalism discussed in Sect. IV the general solution of the evolution equation for the correlation matrix C AB is with U(τ, τ ′ , K) = e − τ τ ′ ds f (s,K) where the term f (τ, K) corresponds to the linear terms of the correlation function in Eqs. (68). f (τ, K) takes into account the dissipative and external forcing sources due to either the expansion rate or the constant electric field. The term g(τ, K) is either the equilibrium value of the correlation function (for the diagonal correlators C AA ) or the couplings between the external electric field and correlators C dd and C T l T l (for the non-diagonal correlators C dT l ). As we pointed out previously, the first line of the solution (77) describes the decay of the initial condition for the correlator while the second line is the genuine contribution coming from hydrodynamic fluctuations [49]. We shall focus in the rest of this section only on the hydrodynamic fluctuating contributions.
Using the exact solution (77) and subtracting the divergent piece, we obtain the finite contribution to the temporal and spatial components of the particle current The coefficients appearing in the RHS of the previous expressions were calculated numerically as indicated in App. C.

Diagrammatic
Hydrokinetics Hydrokinetics In this work we have studied the role of hydrodynamic fluctuations in a conformal fluid with a conserved U (1) charge. The hydrokinetic formalism [11,14,42] has been extended to the regime of finite chemical potential and/or baryon density. The out-of-equilibrium hydrodynamic fluctuating contributions to the particle current renormalize the particle density and the diffusion coeffi-cient. We have verified that the one-loop correction to the Green functions derived from the hydrokinetic and diagrammatic formalisms match exactly in the case of a static homogeneous background. The long time tails of the particle current are generated by contributions associated with sound modes and a mixture of shear and diffusive modes. The sound mode contribution is pro-portional to γ −3/2 η while the shear-diffusive one is proportional to (D + γ η ) −3/2 .
Hydrodynamic fluctuations imply the existence of lower bounds on the transport coefficients. These bounds are purely classical and they are not universal. In fact, the lower bounds depend on the thermodynamic properties and the typical momentum scale at which hydrodynamics breaks down. The lower bound on the diffusion coefficient is enhanced when the value of the ratio of shear viscosity to entropy density is small and thus, the experimental determination of this type of transport coefficient has the potential to determine the impact of hydrodynamic fluctuations in relativistic heavy ion collisions.
In the case of a background undergoing boost invariant Bjorken flow the constitutive relations for the particle current in the presence of hydrodynamic fluctuations are (79c) where the particle density as well as the shear and diffusive transport coefficients are renormalized quantities.
In Table I we summarize the properties of hydrodynamic fluctuations obtained from the different approaches in the static and expanding cases. The renormalized transport coefficients obtained in the diagrammatic approach coincide with the hydrokinetic approach in both homogeneous and expanding fluids. There is a mismatch between the hydrokinetic and diagrammatric results for the cubic divergences that enter the pressure and particle density [6,41]. This mismatch arises from the approximations made in the propagators of the sound modes, which affect the UV behavior of the loop integrals.
We are now in a position to estimate the relative size of fluctuation corrections to the U (1) current. For simplicity we will focus on the spatial components J x /E x = J y /E y . As an estimate of the renormalized conductivity σ = D/α 2 we will use the lower bound (55) derived in Sect. VI 2. This has the virtue that the gradient term and the fluctuation correction scale in the same way with the enthalpy per particlew. We use λ = 1/4 as suggested by simple kinetic estimates. Both the bound on the gradient term and the fluctuation corrections are enhanced by a term proportional to c p . For simplicity we ignore these terms. We express η/s in units of 1/(4π), and the dimensionless proper time τ T in units of a typical freezeout time τ T ≃ 4.5. We find indicating that fluctuations are important, but do not overwhelm the leading term in the gradient expansion at freezeout.
There are a number of issues that can be studied in the future. On the technical side it would be interesting to study convergence properties and higher order asymp-totics of the long time expansion in the presence of fluctuations. This is related to resurgent asymptotics, which has received a considerable amount of attention in the context of both kinetic theory and holographic dualities. From a phenomenological point of view it is crucial to extend this formalism to include critical behavior [14,[50][51][52][53]. Finally it would be interesting to study fluctuations in time-dependent backgrounds in non-relativistic fluids with possible applications to ultracold atomic gases or graphene [6,7,[54][55][56].
are determined from the pressure through the following thermodynamical expressions The energy density is obtained from the Gibbs-Duhem relation Let us then calculate the specific heat at constant volume The isothermal and adiabatic compressibilities are The thermal expansion coefficient is Notice that the thermodynamical relations are satisfied. One can write the adiabatic speed of sound v 2 a in terms of the coefficients calculated above. As a result one gets The isothermal speed of sound is The coefficient α 1 (18a) can be rewritten as Now, the coefficient α 2 (18b) is During the text we find the combination α 1 + α 2 /w which can be written as We can also obtain the coefficient β 2 The susceptibility matrix χ ab relates the variations between the hydrodynamical fields ϕ a = (δǫ, g i , δn) (i=1,2,3) and the variations of the temperature δT , fluid velocity δu i and chemical potential δµ throught the relation ϕ a = χ ab ρ b with ρ a = (δT /T, δu i , δµ− (µ/T ) δT ) [5] and χ ab given by It follows that det χ = H 3 (χ 11 χ 33 − χ 13 χ 31 ). The inverse susceptibility matrix χ −1 ab which is Onsager relations [57,58] require the susceptibility matrix χ ab (A16) to be symmetric. This condition is satisfied in Eq. (A16) because in the gran canonical potential T (∂n/∂T ) µ/T = (∂ǫ/∂µ) T so χ 15 = χ 51 . One can also show that χ 11 , χ 55 ≥ 0 as well as det χ ≥ 0 (see discussion in Sect. 2.5 of Ref. [5]). Furthermore, for time-reversal invariant theories the Maxwell relations together with the Gibbs-Duhem relation lead to the following identities [5] β 1 χ 11 + β 2 χ 51 = H , (A18a) β 1 χ 15 + β 2 χ 55 = n , (A18b) α 1 χ 11 + α 2 χ 51 = 0 , (A18c) α 1 χ 15 + α 2 χ 55 = 1 . (A18d) Thus one can express the thermodynamic coefficients α 1,2 , Eqs. (A11)-(18b), and β 1,2 , Eqs. (A14)-(A15), in terms of the coefficients of the susceptibility matrix χ ab and viceversa. Now for the conformal EOS given by the energy, particle density and entropy densities are given by respectively where g ′ (µ/T ) is the derivative respect to the argument. It follows that ǫ = 3p so v 2 a = 1/3 and thus, the enthalpy density H = ǫ + p = 4p. Using the results of the previous sections we calculate the thermodynamical coefficients β 1 , β 2 , α 1 and α 2 for the EOS (A19) which gives us respectively It is straighforward to show that for the EOS (A19) the combination α 1 + α 2 /w = 0.
which measures the deviation of the full response function from its equilibrium value. The equations for the different response ratio functions read as with A = ±, T l , the terms P AA and D AA are given by Eqs. (71) and we use that ∂ τ [C 0 (τ )/τ ] ≈ −2 1 + v 2 a C 0 (τ )/τ 2 . The general solutions for R AA are where K 0 ≡ K(τ 0 ) and G AA (τ, τ 0 , K) is the Green function. An explicit calculation of the Green functions yields to the following expressions with t = τ ′ /τ , r s = (γ η (τ )τ ) 1/2 K and r d = (D(τ )τ ) 1/2 K. Moreover we define the following functions In the derivation of Eqs. (C4) we used the Bjorken scaling of the temperature and finite chemical potential scale like τ −1/3 in the ideal fluid limit [26]. We approximate K · E R ±± ≈ K · E in the LHS of Eq. (C2) when A = ±. This is justified since we are interested only in the linear response of the . Thus, this term enters as an independent source in the LHS of the corresponding equation for R ±± . By following a similar procedure we find that the differential equation for the response ratio functions R dT l = C dT l / (C 0 /τ ) are with P dT l given by The solution of Eq. (C6) is dτ ′ (ê T l · E) w(τ ′ ) G dT l (τ ′ , τ 0 , K) (R T l T l + R dd ) .

(C8)
The explicit form of the Green functions G dT l are .
with K dT = [(D +γ η )τ ] −1/2 . In order to derive the previous expressions we used explicitly that the Wiedemann-Franz law holds, i.e., D ∼ γ η [59,60]. Now we define the residual function R (r) AB as the difference between the exact solution (C3) [11] and its leading asymptotic series term (which reads off directly from Eqs. (70) or Eq. (72)), i.e., (C10) In terms of these functions the finite residual contributions to the particle current components are given by the following expressions In the following we present the generic method used in this work to obtain the residual contributions to the par-ticle current. We calculate explicitly the mixed diffusivesound contribution appearing in Eqs. (C11a), i.e., with F (τ, τ 0 ) given by (C13) In Eq. (C12) we replace the asymptotic solutions for R T1T1 and R dd which for our purposes is enough since one is interested in the long time asymptotic behaviour. The integrand in Eq. (C13) is divergent since B(t, θ K ) ≈ 1 − t when t → 1 [11]. The finite and linear divergent pieces of the integral (C13) are obtained by rewriting it as follows Eq. (C13) = The integral in the first line is fully convergent while the second one is not and thus, an arbitrary cutoff parameter δ was introduced in order to regulate the divergence. The former regularized integral gives an analytical expression which can be Taylor series expanded and its the leading order O(δ/K dT ) contribution read as By using Eqs. (C15) and (C16) into Eq. (C12) we finally get Therefore we conclude that the linear divergence cancels exactly with the associate late time divergence in Eq. (C12) iff This simple relation resembles the usual M S renormal-ization scheme [61].
By using the same regularization procedure for each integral in Eq. (C11) we calculate the residual hydrodynamical fluctuating contribution to the current J µ . Below we list the numerical values of the relevant integrals needed in this calculation The sum of the different contributions give rise to the numerical values of the particle current quoted in Eq. (78).
For the residual contribution of ∆J τ (78a) we used the results listed in Table 1 of Ref. [11].