Pseudorapidity correlations in heavy ion collisions from viscous fluid dynamics

We demonstrate by explicit calculations in 3+1 dimensional viscous relativistic fluid dynamics how two-particle pseudorapidity correlation functions in heavy ion collisions at the LHC and RHIC depend on the number of particle producing sources and the transport properties of the produced medium. In particular, we present results for the Legendre coefficients of the two-particle pseudorapidity correlation function in Pb+Pb collisions at 2760 GeV and Au+Au collisions at 200 GeV from viscous hydrodynamics with three dimensionally fluctuating initial conditions. Our results suggest that these coefficients provide important constraints on initial state fluctuations and the transport properties of the quark gluon plasma.


INTRODUCTION
The study of the Fourier coefficients v n of the azimuthal particle distribution in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) has led to great insight into the initial state fluctuations and hydrodynamic evolution of the produced system [1]. In particular, it has led to the conclusion that the system has a very small viscosity, close to the lower bound conjectured using AdS/CFT correspondence [2].
Recently, the ATLAS collaboration has presented results on the expansion of two-particle pseudorapidity correlations into Legendre polynomials [3]. The obtained coefficients contain important information on the fluctuations of the particle multiplicity in pseudorapidity.
In this letter we explore how two-particle pseudorapidity correlations from hydrodynamic calculations can give insight into the number of sources for particle production and their correlations, as well as the shear and bulk viscosity of the system. To achieve this, we introduce a simple initial state model that extends the conventional Monte Carlo (MC) Glauber model [4] to include fluctuations of valence quarks in rapidity and thus produces three dimensional fluctuating distributions of net baryon number and entropy density. We then use this model to generate the input for 3+1 dimensional viscous hydrodynamic calculations and compute rapidity distributions of charged hadrons and two-particle rapidity correlations. We then analyze the effect of i) the number of sources and ii) the shear and bulk viscosity of the system on the Legendre coefficients.

INITIAL STATE MODEL AND HYDRODYNAMIC EVOLUTION
Fluctuations in rapidity have been included in hydrodynamic calculations via UrQMD [5,6], EPOS [7], or AMPT [8] initial conditions. To study the effect of fluctuations in rapidity in addition to fluctuations in the transverse plane without being dependent on a compli-cated string with many parameters, we introduce a simple model that is a straight forward extension to the MC Glauber model to include longitudinal fluctuations. In particular, we will employ a Monte Carlo implementation of the Lexus model [9] to provide a simple parametrization of the rapidity distribution of wounded nucleons (or constituent quarks), which is based on experimental observations in proton-proton collisions.
3DMC-Glauber model For each nucleus we sample the three-dimensional spatial distribution of nucleons from a Woods-Saxon distribution. We then sample valence quark positions within each nucleon from an exponential distribution. After overlaying the two nuclei in the transverse plane, separated by the sampled impact parameter b, wounded quarks are determined using the quark-quark inelastic cross-section σ qq , which can be obtained geometrically by requiring that the experimentally determined nucleon-nucleon cross section is recovered. We do, however, treat the quark-quark cross section in heavy ion collisions as an independent, free parameter. When we choose σ qq = 3 mb at √ s = 200 GeV, we can reproduce the experimental multiplicity and pseudorapidity distribution of charged hadrons without including additional negative binomial fluctuations. We employ a Gaussian wounding where quarks are determined to be participants with a Gaussian probability [10,11].
As mentioned above, to obtain the longitudinal distribution of the initial baryon number we employ the Lexus model [9]. Here, the idea is that the rapidity distribution of nucleons in heavy ion collisions can be obtained by linear extrapolation from the distribution in p+p collisions. That distribution is parametrized and fit to experimental data. We use this model for valence quarks such that the probability for a quark with rapidity y P to end up with rapidity y after a collision with a quark with rapidity y T (from the other nucleus) is given by [9]: which corresponds to a flat distribution in longitudinal momentum p L . In the original work [9] λ is the fraction of nucleon-nucleon scatterings that result in a hard collision. Here, we treat λ as a free parameter, regulating the stopping power of the collision. It can be fixed by fitting the net baryon distribution to experimental data. Generally, we find a good fit to the net baryon rapidity distribution when σ qq λ ≈ 2 mb. The initial rapidities are distributed according to the quarks' x value, which is determined by the valence quark parton distribution functions (PDFs). 1 Initial rapidities are thus y = ±ln(x √ s/2m q ) , with the sign depending on whether it is a quark in the projectile (+) or the target (-). We employ a valence quark mass of m q = 0.34 GeV and assume zero transverse momentum initially.
To systematically organize the collisions of all quarks, we have quark pairs collide subsequently with increasing inter-quark distance ∆z. We then convert rapidity to space-time rapidity η s to obtain a three dimensional event-by-event distribution of quarks. To assign a baryon density, this distribution needs to be smeared and we do this by introducing Gaussians in the transverse plane with width σ T = 0.4 fm, and width in space-time rapidity of σ = 0.2.
Next we determine the entropy density distribution. We deposit a fixed entropy between the rapidities of each wounded quark and its last collision partner and smear the edges with half Gaussians of the same width as used for the baryon density distribution. This leads to the following form of the rapidity dependence of the entropy density per wounded quark pair s(y,y P , y T ) = N exp[−θ(−y + min(y P , y T )) where y T and y P are the rapidities of the colliding quarks, and N determines the normalization, which is the same for every "tube" and adjusted to fit the experimental data.
In the transverse plane, we smear the entropy density around the center of mass position of each pair by a Gaussian of width σ T = 0.4 fm.
This way of initializing the entropy density leads to the correct centrality dependence of the total multiplicity. We note that energy and momentum conservation is not explicitly fulfilled, however, we are not including any transverse momentum production or very high momentum partons, which will carry away some of the energy and momentum of the collision and are not part of the bulk medium we are describing.
A different initial state model using random rapidities for wounded nucleons was employed in [15].
Hydrodynamics We integrate the above initial condition into the 3+1 dimensional viscous relativistic fluid dynamic simulation Music [16][17][18][19]. In addition to numerically solving the equations for the conservation of energy and momentum ∂ µ T µν = 0 , and the baryon current ∂ µ J µ B = 0 , we solve the relaxation-type equations derived from kinetic theory [20,21] The transport coefficients τ Π , δ ΠΠ , λ Ππ , τ π , δ ππ , φ 7 , τ ππ , and λ πΠ are fixed using formulas derived from the Boltzmann equation near the conformal limit [21]. At zero baryon chemical potential the ratio η/s will be chosen to be constant in this work, and the temperature dependence of the ratio of bulk viscosity to entropy density ζ/s is parametrized as in [22,23], except that we gradually reduce the constant value at low temperature such that effects of the bulk δf corrections [24] are kept to a minimum. Because we include finite baryon chemical potential µ B > 0, we replace s in the above expressions by (ε + P )/T , motivated by the relativistic limit of the fluidity measure introduced in [25]. The equation of state, which includes finite baryon chemical potential, is constructed by interpolating the pressures of hadronic resonance gas and lattice QCD [26,27].
We leave a detailed description of the initial state model and the newly constructed equation of state for a longer paper in the future.
After Cooper-Frye freeze out at an energy density of 0.1 GeV/fm 3 , the calculation of thermal spectra including δf corrections [24] for shear and bulk viscosities, and resonance decays, we obtain the final hadron spectra as functions of transverse momentum p T and pseudorapidity η p .

TWO PARTICLE RAPIDITY CORRELATIONS
The p T integrated (p T > 0.5 GeV) event-by-event rapidity distributions dN/dη p are then used to determine the two-particle rapidity correlation function [28] C(η 1 , where N (η) = dN/dη p . We follow the ATLAS collaboration [3] and remove the effect of residual centrality dependence in the average shape N (η) by redefining the correlation function as [29] C N (η 1 , where Importantly, the resulting distribution is then normalized such that the average value of C N (η 1 , η 2 ) is one. Following [29,30] C N (η 1 , η 2 ) is then decomposed into Legendre polynomials. The Legendre coefficients are given by a n,m = C N (η 1 , η 2 ) where T n (η p ) = n + 1/2 P n (η p /Y ) and P n are the standard Legendre polynomials. The a n,m are related to the Legendre coefficients of the single particle distribution a n via a n,m = a n a m , where the a n are defined through N (η)/ N (η) = 1 + n a n T n (η) [3].

RESULTS
For the calculations presented in this work, the parameters of the initial state model are adjusted to fit the measured dN/dη p [31][32][33] and dN/d 2 p T [34] distributions of charged hadrons, as well as the net-baryon rapidity distribution [35,36], when using the T -dependent bulk viscosity and η/s = 0.12. In particular we use a quark-quark cross section of 3 mb and λ = 0.66. The hydrodynamic evolution starts at τ 0 = 0.38 fm/c. We present results for the |a n,m | of the fluctuating initial entropy density distribution (averaged over the transverse directions) and the final |a n,m | obtained from the produced charged hadrons with |η 1 |, |η 2 | < 2.4 for 20-25% central Pb+Pb collisions at 2760 GeV in Fig. 1. Apart from the a n,n = a 2 n , the next largest correlations are those for m = n + 2. For this combination of Legendre indices, the a n,m are negative, in qualitative agreement with the experimental data from ATLAS [3]. The a n,n+1 vanish in symmetric collision systems.
The final |a n,m | show a reduction relative to those obtained from the initial entropy density distribution. This reduction increases with n, m, showing that the hydrodynamic expansion smears out the shorter range correlations more efficiently. The effect of viscosity is to slightly increase the lower a n,m and decrease the higher ones. This can be understood as the effect of diffusion that destroys shorter scale structures in N (η p ), while the reduction of the longitudinal pressure keeps the pseudorapidity distribution more compact as shown previously [18,37], leading to smaller suppression of the lower a n,m that are sensitive to long range correlations. Generally, the n dependence of the calculated final state |a n,m | is stronger than observed in the experimental data [3]. The calculated a n,m agree with the experimental data for small n = m = 1, but are smaller than the experimental values for all other n, m.
The initial state values of |a n,m | agree fairly well with the experimental data, which is interesting, however, a direct comparison is difficult in this case.
The final calculated |a n,m | for large n, m are much smaller than the data for any choice of transport parameters. Since we are only interested in the effects of flow in this study we have not included any non-flow effects, such as short range correlations from resonance decays (see e.g. [38]). The large difference between the calculation and the experimental data for larger n, m indicates that these have a significant effect for n, m ≥ 3. We leave a quantitative study of this effect using statistical hadronization for future work.
The same effect is seen in 65-70% central collisions shown in Fig. 2. Here the coefficients extracted from the initial distribution already underestimate the experimental data and after hydrodynamic evolution most |a n,m | are largely underestimated.
The agreement of the lower a n,m with experimental data for 20-25% centrality, but disagreement for 65-70%, indicates that the simple initial state model does not describe correctly the centrality dependence of the number of sources, which dictates the degree of fluctuations. In Fig. 3 we show predictions for |a n,m | in Au+Au collisions at 200 GeV. We find an overall increase of all |a n,m | compared to the Pb+Pb result at 2760 GeV by approximately 50 %. We further study the effect of varying the shear viscosity to entropy density ratio η/s from 0.12 to 0.2, and the effect of bulk viscosity separately. As opposed to the change of η/s from 0 to 0.12 for Pb+Pb collisions (see Fig. 1), we generally find an increase (or no change) of all |a n,m | when increasing η/s. This indicates that the reduction of the longitudinal pressure dominates over the increased diffusion. Including bulk viscosity has a similar effect as increasing the shear viscosity, however, the effect is not large compared to the statistical errors of our results. Next we explicitly study the dependence of the a n,m on the number of sources in the initial state. Fig. 4 shows how the a n,m of the initial entropy density distribution are increased when reducing the number of sources by employing nucleons instead of valence quarks. As alluded to above, this effect is expected, because fewer sources lead to more fluctuations. Finally, in Fig. 5 we present a prediction for |a n,m | extracted from the net-baryon distribution. We see that the hydrodynamic evolution does not dampen the higher a n,m more than the lower ones as was the case for the charged hadron distribution. We have checked that the ideal results agree with the ones shown for η/s = 0.12 and (ζ/s)(T ) within errors. We argue that measuring these quantities e.g. at RHIC and comparing to various models could reveal important details on the mechanism of baryon stopping and transport in heavy ion collisions.

CONCLUSIONS AND OUTLOOK
We have presented the first calculation of two-particle pseudorapidity correlations in a viscous fluid dynamic framework and studied the dependence of our results on the initial state and the transport properties of the medium. We have shown that the number of particle producing sources affects the Legendre coefficients a n,m , as does the presence of hydrodynamic evolution. The shear and bulk viscosities of the medium also modify the a n,m coefficients. The effect of viscosity is two-fold. On the one hand, it reduces the longitudinal pressure, leading to an increase of the a n,m . On the other hand, the diffusion has the opposite effect. This effect becomes larger for shorter range correlations, i.e., larger n, m. The underestimation of the experimental |a n,m | coefficients in Based on our results we conclude that the study of a n,m together with the Fourier coefficients v n at various collision energies has the strong potential to constrain the initial state particle production mechanism as well as the transport properties of the medium.
We further propose to perform the same experimental analysis using net-baryon or net-proton distributions to constrain in detail the mechanism of baryon stopping and transport in heavy ion collisions.
It will also be very interesting to analyze the twoparticle pseudorapidity correlations using a more sophisticated initial state model, such as the IP-Glasma [19,39,40] in the future. This will require an extension of the boost-invariant Glasma picture to three dimensions.
Note: An independent study of the two-particle pseudorapidity correlations within a similar framework and leading to similar conclusions has appeared simultaneously with our work [41].