Pseudorapidity dependent hydrodynamic response in heavy-ion collisions

We propose a differential hydrodynamic response relation, $V_2(\zeta)=\int d\xi G(\zeta-\xi) \mathcal{E}_2(\xi)$, to describe the formation of a pseudorapidity dependent elliptic flow in heavy-ion collisions, in response to a fluctuating three-dimensional initial density profile. By analyzing the medium expansion using event-by-event simulations of 3+1D MUSIC, with initial conditions generated via the AMPT model, the differential response relation is verified. Given the response relation, we are able to separate the two-point correlation of elliptic flow in pseudorapidity into fluid response and two-point correlation of initial eccentricity. The fluid response contains information of the speed of sound and shear viscosity of the medium. From the pseudorapidity dependent response relation, a finite radius of convergence of hydrodynamic gradient expansion is obtained with respect to realistic fluids in heavy-ion collisions.


I. INTRODUCTION
In high energy heavy-ion experiments, the fluid nature of quark-gluon plasma (QGP) has been well demonstrated through the measurements of flow harmonics in multi-particle correlations [1]. These flow harmonics can be understood as the fluid response to the decomposed azimuthal modes associated with the initial state geometrical deformations, giving rise to a series of response relations (cf. [2,3] for recent reviews). In particular, the eccentricity of the initial density profile E 2 , which characterizes the elliptic deformation of the QGP fireball at initial time, and elliptic flow V 2 [4], which characterizes the asymmetric emission of final state particles in azimuthal angles, are linearly correlated: Both V 2 and E 2 fluctuate from event to event, but this linear relation has been found valid for the mid-central collisions via event-by-event hydrodynamic simulations [5][6][7]. In Eq. (1), the response coefficient G 0 depends on the fluid properties of QGP, which is a real constant in one specified centrality class. A consistent suppression of the response coefficient has been observed when dissipative corrections in the QGP become larger. Eq. (1) has lead to many remarkable applications in heavy-ion collisions, such as the background subtraction of chiral magnetic effects [8] using a selected collision geometry [9,10]. Nevertheless, as one relation between the global azimuthal asymmetry of initial and final states, knowledge of the local structure of the QGP medium is absent. In this letter, we propose a generalization of Eq. (1) to a pseudorapidity dependent hydrodynamic response. Our generalization is partly motivated by the experimental developments in heavy-ion collisions, where * cliyan@fudan.edu.cn pseudorapidity dependent flow harmonics have been explored [11][12][13]. Besides, the pseudorapidity dependent hydrodynamic response is expected to play a more significant role in the beam energy scan program carried out at the Relativistic Heavy-Ion Collider [14], since at lower collision energies, approximations for the fluid dynamics along the longitudinal direction with respect to Bjorken symmetry are no longer valid in theoretical models [15].
A pseudorapidity dependent generalization of Eq. (1) reveals the local properties of the QGP medium, which is of extremely significance. First, a pseudorapidity dependent hydrodynamic response allows one to detect the fluid locally, so that the transport properties of the QGP medium can be examined differentially. One such example is the attempt to extract a temperature dependent ratio of shear viscosity over the entropy density, η/s(T ) [16]. Besides, if the local gradients in pseudorapidity are accessible, the convergence behavior of the hydrodynamic gradient expansion can be analyzed in realistic QGP medium with non-trivial 3+1 dimensional expansions, which would substantially extends the recent developments with respect to systems undergoing one dimensional Bjorken expansion (cf. [17][18][19]).

II. FORMULATION
We first present a step-by-step derivation of the generalization of Eq. (1). The only assumption we have is the fluidity of QGP, so that long wave-length (small wavenumber) modes dominate.
In heavy-ion collisions, from each collision event the single-particle distribution that describes the probability of particle emission with azimuthal angle (φ = atan(p y /p x )) and pseudorapidity (ζ = atanh(p z /|p|)) dependence can be decomposed into Fourier modes arXiv:1907.10854v3 [nucl-th] 30 Jan 2020 where the Fourier coefficient defines the complex harmonic flow V n (ζ) ≡ v n (ζ)e inΨn(ζ) . It is worth mentioning that the definition of flow harmonics in Eq. (2) differs from the usual one by a factor of normalized multiplicity distribution. Given the definition in Eq. (2), the integrated flow coefficient is then V n = ∞ −∞ dζ V n (ζ) . Throughout this letter, we shall focus on the second harmonic, namely, elliptic flow V 2 (ζ), while generalizations to other flow harmonics are straightforward.
Initial eccentricity E 2 is a complex quantity, which characterizes the elliptical asymmetry of the threedimensional initial density profile ρ( x ⊥ , ξ). Note that ξ = atanh(z/t) is the space-time rapidity, to be distinguished from the pseudorapidity ζ. In order to generalize the response relation along the longitudinal direction, we consider the definition for a ξ-dependent eccentricity as so that the azimuthal deformation of the density profile at a specified space-time rapidity ξ is quantified. Eq. (3) is not a unique definition that introduces dependence on space-time rapidity in eccentricity [20,21], but Eq. (3) satisfies the condition: , and more importantly, for m = 0, 1, 2, . . ., which will be found useful. Correspondingly, we generalize the linear response relation Eq. (1) as 1 where the real constant response coefficient is replaced by a response function G(ζ − ξ). This function is expected entirely determined by the fluid properties of the QGP medium, which does not fluctuate in one specified centrality class where the system multiplicity is a fixed constant. Eq. (5) is a non-local and differential response relation, relating the eccentricity at one space-time rapidity ξ at initial time to the flow observed at pseudorapidity ζ through the evolution of hydrodynamic modes. Note that the form of the response function implies a boost invariant background, on top of which the apparent breaking of the boost invariant symmetry in realistic heavy-ion collisions can be accounted for as perturbations introduced by the decomposed modes in E 2 (ξ).
To identify the response function, it is advantageous to work with a Fourier transformation, namely, for the elliptic flowṼ 2 (k) = ∞ −∞ dζ V 2 (ζ)e −ikζ and initial ec- In terms of the wave-number k, Eq. (5) becomes, In the long wave-length limit with |k| k * , corresponding to the hydrodynamic regime, one is allowed to expand the response function in series of k. For instance, up to the second order the expansion is which amounts to in the ζ-space. By construction, these expansion coefficients are real constants, to be determined later. Because of the boost-invariant background, reflection symmetry is implied so that odd order response coefficients must vanish.
Several comments are in order. First, the leading order term can be identified as the response of the global elliptic asymmetry, while higher order terms (O(k) and beyond) do not contribute after an integration over pseudorapidity. Therefore, Eq. (1) will be recovered from Eq. (8) after an integration over ζ. Secondly, these higher order terms in the series expansion correspond to the hydro gradients, which characterize modes satisfying a gapless dispersion relation. As is evident in the response relation, these hydrodynamic gradients are specified as derivatives in space-time rapidity, as a consequence of perturbations of the initial state density profile along ξ. Thirdly, validity of the hydrodynamic gradient expansion is determined by its convergence behavior. The radius of convergence k * , for instance, is finite if the hydrodynamic gradient expansion is convergent [22]. If, on the other hand, the gradient expansion is asymptotic, practical applications of the response formulation would rely on a Borel resummation method [17]. As will be clear in Appendix A, the existence of a finite k * results in an effective regularization of the eccentricity distribution along space-time rapidity, which smooths out local structures. More detailed discussion on the convergence behavior and k * will be given in Sec. III B.
To access these expansion coefficients G n 's, it is natural to define new sets of flow observables and initial eccentricity variables weighted with powers of ζ (or ξ), These new variables are defined according to the series expansion form Eq. (8), reflecting the nature of hydrodynamic mode evolution. Especially, E (n) 2 provides an alternative mode decomposition of the initial eccentricity distribution along the space-time rapidity, comparing to those developed on grounds of initial state geometrical fluctuations [23,24]. Using integration by parts repeatedly, the generalized linear response relation Eq. (5) can be written in an alternate form as Note that Eq. (4) has been taken into account implicitly to obtain Eq. (10). Note also that the leading order relation V is the linear response relation presented in Eq. (1), and the response coefficient can be calculated in event-by-event hydrodynamic simulations [7]: , where the angular brackets indicate event average. Given the leading order response coefficient G 0 , one realizes a linear response relation between E with the slope identified as G 1 . In a similar way, a set of linear relations can be realized between E , with G n calculated recursively as

III. EVENT-BY-EVENT HYDRODYNAMIC SIMULATIONS
To test the analytical formulations of the longitudinal hydrodynamics response, we carry out realistic event-byevent hydrodynamic simulations, for the Pb-Pb collision system with √ s N N = 2.76 TeV at the LHC, for events in centrality class of 30-40%. This is realized by the 3+1D MUSIC [25][26][27], with respect to random 3D initial conditions generated by the AMPT model [28,29]. The initial density profile is obtained in a similar strategy as in Ref. [30], with an overall constant factor adjusted to reconcile the mismatch of viscosity in AMPT and viscous hydrodynamics [20,31]. Note that non-trivial longitudinal density distribution with fluctuations in the AMPT model have been implemented. In Fig. 1, the distribution of eccentricity magnitude along space-time rapidity in one collision event is shown as the green dash-dotted line, where bumps appear as a consequence of longitudinal fluctuations. MUSIC solves the second order viscous hydrodynamics with respect to a 3+1D expanding system, in which a cross-over from the QGP to hadron gas is implied in an equation of state obtained from lattice QCD. We calculate elliptic flow V 2 from thermal pions. Throughout our analyses, statistical errors from simulations are estimated using the Jackknife resampling method. In the current study, we only consider shear viscous corrections with a constant ratio of shear viscosity to the entropy density, η/s, while second order transport coefficients are fixed accordingly. More details on the numerical simulations of MUSIC, such as the equation of state, the freeze-out prescription, etc., can be found in Ref. [25][26][27] and references therein.
Given the results from hydrodynamic simulations, the constant response coefficients can be calculated, as well as the linear relations between E (0) 2 . Scatter plots in Fig. 2 present these linear relations respectively from the real part (red points) and the imaginary part (yellow points), from event-by-event hydro simulations of approximately 5000 events with a constant η/s = 0.08. For even orders up to n = 16, apparent linear relations are shown, with slopes compatible with those calculated with respect to Eq. (11). 2 We notice that slopes of odd orders are consistent with zero within errors, as anticipated.
, with coefficients Gn's solved recursively with respect to Eq. (11) (slopes of solid lines), using event-by-event 3+1D hydro simulations with η/s = 0.08. Red and yellow points are for the real and imaginary parts in the linear relation respectively. Slopes of odd orders (not shown) are consistent with zero.

A. Two-point correlation of V2
The pseudorapidity dependent response relation Eq. (5) is non-local, which implies that the generation of V 2 at one pseudorapidity receives contributions from other space-time rapidities. But the response along longitudinal direction is limited by the speed of sound in fluid systems. This effect can be shown in the analysis of two-point correlations. We define to characterize the length of the two-point correlation measured via elliptic flow at different pseudo-rapidities.
With the response relation derived in Eq. (9) and Eq. (10), in particular, considering the fact that G 1 = 0, it can be proved that The length of the initial state eccentricity two-point correlation (∆ξ) 2 is defined according to Eq. (12) through E 2 (ξ). Eq. (13) states that the increase of the two-point correlation length in elliptic flow comparing to that in the initial eccentricity is purely an effect of fluid dynamics. This can be understood as a direct consequence of sound propagation, reflected in the ratio G 2 /G 0 , that initial perturbations propagate in speed of sound c s along the longitudinal direction [32]. etc., in dissipative fluids the two-point correlation at the final stage is reduced. Fig. 3 shows the numerical results of about 1000 events from hydrodynamic simulations with different values of η/s. The obtained length of the two-point correlation of V 2 is plotted as a function of η/s, which agrees with expectation from Eq. (13) within statistical errors. When η/s → 0, the length of two-point correlation approaches an upper bound determined by the sound horizon of QGP medium. Reduction of the correlation length, or the ratio G 2 /G 0 , due to shear viscous correction is clear and systematic. Subject to boundary corrections due to the experimental acceptance along pseudorapidity, (∆ζ) 2 is a measurable in heavy-ion experiments, providing a novel probe of η/s in the QGP medium.

B. The radius of convergence
With the expansion coefficients calculated from eventby-event hydrodynamic simulations, we are able to explore the convergence behavior of the hydrodynamic gradient expansion in Eq. (7). Especially, the gradient expansion corresponds to the fluid properties in a realistic, 3+1D expanding QGP medium obeying lattice equation of state.
Owing to the reflection symmetry, we ignore contributions from odd orders of the expansion. For n > 2, we notice that the sign of G n of even orders flips. In Fig. 4, the calculated even order coefficients are shown, from hydrodynamic simulations of roughly 1000 events, using η/s = 0.001 (ideal) and η/s = 0.2, respectively. For n > 4, these G n 's grow exponentially: log |G n | ∝ αn, which implies a finite radius of convergence. Therefore, the hydrodynamic gradient expansion in this current case  is convergent. The slope in Fig. 4 determines the radius of convergence. Alternatively, one may also apply a diagonal Padé approximation for the gradient expansion. Convergence behavior can thus be revealed from the singularity structure of the Padé approximation. The resulted poles in the complex k 2 -plane are shown in Fig. 5. When n > 70, the singularity structure of the Padé approximation stabilizes. In particular, in addition to those poles lying around a circle, there exists one pole closest to the origin on the positive real axis, which determines the radius of convergence.
From both ways, the obtained radii of convergence are finite. Numerical values of the radii are consistent, which we identify as the upper bound of the gradient expansion k * . From ideal hydrodynamic simulations, we find k * = 1.209 (slope) and 1.209 (Padé), while for the case of η/s = 0.2 we find k * = 1.192 (slope) and 1.189 (Padé). When shear viscosity decreases, there seems to be a trend that k * grows [22,33]. However, due to the effect of expansion, we find that k * saturates to a constant in an ideal fluid even though the mean free path approaches zero. For a medium system away from local equilibrium, it is known that the gradient expansion consisting of higher order dissipative corrections is asymptotic [17][18][19]. Origin of such divergence can be intuitively attribute to the fact that the number of higher order viscous terms has a factorial increase. However, it is different in the current case concerning perturbations on top of the medium system in or close to local equilibrium, where gradients are taken into account with respect to the dispersion relation [22,33,34].
In the current study, the existence of a finite radius of convergence of the gradient expansion has an important physical interpretation. The upper bound k * of the expansion implies a length scale limit in the space-time rapidity of the initial eccentricity profile E 2 (ξ). This is the scale that determines the resolution of fluid response to initial state local structures. Although for a static medium this effect can be generically understood as the fact that the fluid behavior of a medium system must be visible at a scale greater than the mean free path and one expects k * ∼ 1/l mfp , in the medium with non-trivial 3+1D expansion it becomes more complicated. As a consequence of the hydro response and finite k * , a regularization scheme for the initial eccentricity profile can be introduced, leading to a regularized eccentricity profile E Reg 2 (ξ). Derivation of the regularization scheme is given in Appendix A. In Fig. 1, one finds in the resulted E Reg 2 the local perturbations are smoothed out, giving rise to a consistent shape comparing with the distribution of flow where bumpy structures are merely seen.

IV. SUMMARY AND DISCUSSIONS
We have derived the generalized formulation of pseudorapidity dependent hydrodynamic response. The generalized form is a non-local relation that interprets the formation of elliptic flow at different pseudo-rapdities as a consequence of longitudinal fluid response, which provides a novel picture of longitudinal flow generation that differs from some previous studies (cf. Ref. [35,36]). The formulation is applicable to other flow harmonics as well, when complexities involving nonlinearities do not arise [7,37].
The longitudinal fluid response can be quantified by a set of constant coefficients G n 's, and a new set of flow observables in Eq. (9). This is the main result of this letter, which has been shown valid through realistic event-byevent hydrodynamic simulations for a 3+1D expanding system in heavy-ion collisions. The second order expansion coefficient G 2 is found to be sensitive to the sound mode propagation, contributing to the increase of twopoint correlation length of the elliptic flow. Shear viscous correction in the fluid reduces two-point correlation, as expected from the damping of sound. Therefore, twopoint correlation of flow provides a new measurable to investigate the shear viscosity of the QGP medium in experiments.
Within the formulation of pseudorapidity dependent hydrodynamic response, the expansion coefficients G n 's can be calculated with respect to a realistic expanding QGP medium. Given these expansion coefficients, hydrodynamic gradient expansion can be analyzed. In this current study, we solved 3+1D MUSIC with respect to fluctuating initial conditions, which represents the fluid response to perturbations on top of a fluid medium that obeys second order viscous hydrodynamics. Generation of the elliptic flow from gradients along the longitudinal direction thus reflects the evolution of sound and shear diffusive modes, as in the hydrodynamic dispersion relations. A finite radius of convergence is observed, which confirms the finding that hydrodynamic gradient expansion associated with the dispersion relation is convergent [22,33,34].