Rapidity-dependent eccentricity scaling in relativistic heavy-ion collisions

There is a well-established relation between the spatial asymmetry in the initial stage of a heavy-ion collision and the final momentum anisotropy, which allows for a separation of effects from initial conditions vs.~later evolution and has proved exceptionally powerful. However, until recently it has only been studied in two dimensions --- either through boost-invariant simulations or studying only quantities at mid-rapidity. We explore an extension to 3 dimensions, in order to determine whether a similar understanding can be obtained for the rapidity dependence of the collision system. In particular, we introduce rapidity-dependent eccentricities and investigate a trivial extension of the 2D eccentricity scaling of elliptic and triangular flow, as well as a way to systematically improve these initial-state estimators. We then explore the dependence of the resulting response coefficients on shear viscosity and initial total energy.


I. INTRODUCTION
The evolution of a relativistic heavy-ion collision is a complicated process, well described by simulations that utilize relativistic fluid dynamics. However, despite the non-linear nature of these equations, simple relations can be found between the state of the system at the onset of hydrodynamic behavior, and the final spectrum of detected particles [1]. In the simplest case a linear relation between a measure of the initial spatial azimuthal asymmetry ε n and the final momentum anisotropy V n can be written as where all relevant information about the initial state is contained in a single quantity ε n , and all information about the subsequent dynamics of the system is contained in the response coefficient κ. This separation of effects is quite powerful, and has allowed for numerous insights [2][3][4][5][6]. However, these relations were developed in systems with two effective dimensions -first in systems with an assumed invariance under longitudinal boosts [2], and later in 3D simulations, but analysis restricted to midrapidity [1]. Only recently has an extension to a full 3 dimensions has been considered [7].
Here we investigate possible extensions of this mapping of the initial to final state, to include the dependence on rapidity.

II. INITIAL TO FINAL STATE MAPPING
The essential idea of the mapping from initial to final state is the ansatz that hydrodynamic evolution of the system is more sensitive to the large-scale structure of the initial conditions, as compared to features at small length scales.
This separation of scales is accomplished via Fourier transform of the initial energy density ρ, and large length scales isolated via a truncated Taylor series around transformed variable k = 0. That is, The truncated set of complex cumulants {W n,m } are then used to construct initial state estimators of the final distribution of particles in momentum space, which are traditionally separated into azimuthal rotation modes where φ is the azimuthal angle of a particle's momentum. The leading order relation is a linear proportionality between V n and the lowest cumulant W n,n . Typically, one constructs dimensionless "eccentricities"; for example, In a coordinate system such that W 1,1 = 1, this has the explicit form for n = 2, 3 with The resulting linear relation (1) between V n and ε n has been shown to be an excellent approximation for n = 2 and n = 3 [1]. Other harmonics can have significant contributions from non-linear combinations of cumulants, so for simplicity we focus here on these two harmonics. The extension to estimators with non-linear or higher order cumulant terms is straightforward.

III. EXTENDING TO 3D
All the above relations are implicitly 2-dimensional, representing dynamics in the transverse plane, and with the longitudinal dimension ignored. The justification for this is that, at the highest collision energies, the system should have an approximate invariance under longitudinal boosts. That is, the initial energy density ρ should be approximately independent of spacetime rapidity η ≡ arctan z t (8) and the final particle distribution (4) independent of rapidity such that in fact, boost invariant simulations have been quite successful at describing experimental data near mid-rapidity, Y 0 [8,9]. Nevertheless, numerous rapidity-dependent measurements have been made, which can be utilized to gain insight into QCD dynamics. Because of this, it is of great interest to extend the 2D relations to 3D, to give a better understanding of the evolution of the system and to provide guidance for how to separate, e.g., initial state and final state effects.
To that end, we must consider the dependence of V n on rapidity, and construct estimators from the initial state that contain information about spacetime rapidity η.
A straightforward way to do this is to simply calculate eccentricities independently at each value of rapidity η, to construct rapidity-dependent eccentricities: To first approximation, we might expect the results to be local in rapidity. That is Of course, we don't expect the results to be perfectly local in rapidity. Hydro evolution can propagate information from one rapidity η in the initial state to another rapidity at later times. In addition, particles in a fluid (e.g., at freeze-out) do not all have rapidity Y equal to the spacetime rapidity η of the fluid cell, but instead have a distribution. Similarly, unstable particles can decay to daughter particles that do not have exactly the same rapidity as the parent.
Nevertheless, we do expect a quasi-locality. All of these effects should be largest at small rapidity difference |η − Y |, with an exponential decrease such that properties of the initial state at far forward rapidity η have essentially no effect on the particle distribution at far backward rapidity Y .
This decreasing effect can be encoded in a gradient expansion: where the primes on the eccentricities indicate a derivative with respect to spacetime rapidity η and where the eccentricities are then evaluated at η = Y . Note that symmetry dictates that the hydro response functions must be symmetric under reflection, and so only even derivatives of η can appear.

IV. NUMERICAL TESTS
In order to determine whether these estimators describe the evolution of a heavy-ion collision system, we perform hydrodynamic simulations. For this we use the MUSIC [10-12] 3+1D hydrodynamic code. Specifically, all simulations use the equation of state is s95p-v1 and vanishing bulk viscosity. The initial viscous tensor is set to 0, and the initial fluid velocity is set to boost-invariant Bjorken flow. For simplicity, we calculate the spectrum of direct pions, without contribution of resonance decays.
We aim to study the hydrodynamic response to initial conditions. To that end, we construct a parameterized toy initial condition that can be varied by hand in order to probe the resulting response. Specifically, we choose a simple Gaussian energy density distribution, which can be deformed in the desired rotational harmonic n: The parameterε n controls the azimuthal asymmetry, such that the eccentricity can be directly controlled while the triangularity can also be easily controlled, having the approximate relation In general, we will allow the normalization A and asymmetry ε n to vary with rapidity η.

V. RESULTS AND DISCUSSIONS
First we perform ideal (zero viscosity) boost-invariant simulations, with all parameters independent of rapidity η. This allows us to verify the traditional linear relationship (1), but also to extract the response coefficients κ n , which should remain the same when the system is boost variant.
In the following, we fix ρ = 3f m and A = 50f m −4 , corresponding roughly to a typical collision at RHIC. By varying the asymmetry parameters ε n , we indeed find a linear behavior, from which we can extract values With these values in mind, we can finally explore rapidity-dependent systems. For example, we can choose a linear function ε n (η), and compare it to the resulting rapidity-dependent V n (Y ), as shown in Fig. 1. Using the same κ n as found in the boost invariant case, we find that the estimator (1) provides an excellent description of the results. Not only is the result a linear function of rapidity, but the proportionality constant between V n and ε n is the same as found previously. We found this to be to be true independent of the slope of the rapidity dependence.
Note that, in this case, all higher terms in Eq. (15) vanish, since all even derivatives of ε n (η) vanish. To test to what extent this gradient expansion describes the hydrodynamic response, we can include curvature to the initial eccentricity.
In Fig. 2, we show the corresponding results for quadratic dependence on rapidity, with various values for the curvature, compared to the first-order estimator (1).
The first-order estimator (1) continues to give a reasonable description of the results, but the accuracy gets poorer as the curvature increases. In fact, a close inspection shows that the difference between the anisotropic flow and the estimator is constant in rapidity, but grows proportionally to the curvature of the initial eccentricity. This is exactly what is expected from the proposed gradient expansion (15), which in this case has exactly 2 terms: From these results we can extract the values of κ n . We find The gradient corrections are small, as expected, but not necessarily negligible, depending on the curvature of the initial eccentricity. Next, in Fig. 3 we consider a Gaussian rapidity dependence. In this case, an infinite number of derivatives are present, and we can probe Eq (15) in more detail. Again, we see that the leading-order estimator (1) offers a reasonable description of the resulting anisotropic flow, but a better description can be obtained by adding a subleading term, Eq. (21). We find that the contribution from the next correction (from the fourth derivative) is suppressed even more than the second order term, such that it is negligible in this case.

A. Viscosity dependence
Next, we include a non-zero shear viscosity, in order to probe the behavior of the new response coefficients κ n . Here, we use a Gaussian rapidity dependence, exactly as in Fig. 3, in which case we can extract both κ n and κ n from a single simulation. Figure 4 shows the resulting  As expected, viscosity suppresses anisotropic flow, represented by a decrease in leading response coefficient κ n . This dependence is shown explicitly in the top panels of Fig. 5.
The bottom panels of Fig. 5 show an interesting new result. The subleading response coefficient κ 2 increases with viscosity, indicating an increased sensitivity to rapidity regions η = Y when viscosity is larger. In contrast, κ 3 , decreases with viscosity (although the relative contribution κ 3 /κ 3 still increases, as for the n = 2 case).
We have thus verified the ansatz (15), and investigated the dependence of new response coefficients to viscosity. However, in realistic simulations, not only the asymmetry depends on rapidity, but also the total energy and multiplicity. An important question is whether in this case we can still understand the hydrodynamic response to initial conditions in a simple way.

B. Rapidity dependence via energy dependence
Even in boost-invariant simulations, the anisotropic flow V n resulting from an initial condition with spatial anisotropy ε n will change if the total energy/entropy of the initial condition changes significantly. That is, the hydrodynamic response κ n depends on initial energy. The natural question is whether this effect alone can describe the hydrodynamic response of systems with rapidity-dependent initial energy.
To verify this, we first calculate κ n (E) in boostinvariant simulations, where in our initial condition the total energy per unit rapidity is given by with τ = 1f m the initial proper time. The result can be seen in Fig. 6. Then we verify whether the rapidity-dependent results are described by We choose a Gaussian profile for the initial energy as a function of rapidity and compute the resulting elliptic and triangular flows. In Fig. 7 we can see that the behavior is reasonably captured by Eq. (25), and locality in rapidity can indeed be maintained.

VI. CONCLUSIONS
We showed that the traditional eccentricity scaling of heavy-ion collisions can be extended to a rapiditydependent form, and systematically improved to include non-local influence in rapidity. We used hydrodynamic simulations to verify this framework, and to extract the viscosity and total energy dependence of newlyintroduced response coefficients.
We found that an increase in viscosity makes elliptic flow and triangular flow less dependent on the local value of eccentricity (i.e., the eccentricity at spacetime rapidity η equal to particle rapidity Y of the flow in question), and more sensitive to the eccentricity in the neighborhood surrounding rapidity Y .
An alternative ansatz was presented in Appendix A, but is unable to simultaneously capture the behavior of systems with large and small rapidity gradients.
In realistic systems, the total energy of the system will have a non-trivial dependence on rapidity as well. We find that locality in rapidity is maintained as long as the dependence of the response coefficients κ n on energy is taken into account.
These results should allow for new insight into the dynamics of heavy ion collisions, such as a more direct isolation of rapidity-dependent processes in the initial stages.

Appendix A: An alternative ansatz
Rather than assuming an expansion in gradients, (15), one might try the following ansatz: where the effects of eccentricity ε n at each pseudorapidity η influences the anisotropic flow V n at rapidity Y additively, according to the kernel W n (Y − η). A perfectly local relation, Eq. (1), corresponds to W n (Y −η) ∝ δ(Y − η), though in general one expects a function with finite width. In principle, one could find the function W n by setting ε n (η) = δ(η), in which case In practice, we can approach this limit by choosing a step function for ε(η) and reducing the width to the size of one hydro cell, in our case ∆η = 0.1.
Then, noting the discrete grid in η, we have We show the resulting W 2 and W 3 in Fig. 8, for various values of shear viscosity.
We can see that the influence of ε n (η) extends not much more than 1 unit of rapidity, with a smooth shape (that is not in general well described by a Gaussian).
If the ansatz is correct, we should be able to make contact with the previous results from Eq. (15) by expanding ε n (η) in a Taylor series around Y = η, where κ n and κ n are the previous response coefficients. FIG. 9. Resulting values of κn and κ n for each value of shear viscosity η/s Using the results of Fig. 8, we compute these response coefficients as a function of viscosity. The results are shown in Fig. 9. The coefficients κ n behave as expected, and agree reasonably well with the previous results of Fig.5. However, the higher coefficients κ n do not agree.
We find that it is difficult to accurately extract these coefficients with these large gradient simulations, and that the results do not accurately reflect the results of simulations with smaller gradients. While the ansatz A1 may still work in some regime of validity, it cannot simultaneously describe the results from simulations with initial eccentricity distributions spanning from very narrow to very wide, and we cannot therefore directly extract a universal weight function W (Y − η) as attempted here.