Longitudinal eccentricity decorrelations in heavy ion collisions

In heavy-ion collisions, the harmonic flow $V_n$ of final-state particles are driven by the eccentricity vector ${\mathcal{E}}_n$ that describe the shape of the initial fireball projected in the transverse plane. It is realized recently that the structure and shape of the fireball, and consequently the ${\mathcal{E}}_n$, fluctuate in pseudorapidity $\eta$ in a single event, ${\mathcal{E}}_n(\eta)$. This leads to eccentricity decorrelation between different $\eta$, driving the longitudinal flow decorrelations observed in the experiments. Using a Glauber model with a paramerterized longitudinal structure, we have estimated the eccentricity decorrelations and related them to the measured flow decorrelation coefficients for elliptic flow $n=2$ and triangular flow $n=3$. We investigated the dependence of eccentricity decorrelations on the choice of collision system in terms of the size, asymmetry and deformation of the nuclei. We found that these nuclear geometry effects lead to significant and characteristic patterns on the eccentricity decorrelations, which describe the measured ratios of the flow decorrelations between Xe+Xe and Pb+Pb collisions. These patterns can be searched for using existing experimental data at RHIC and the LHC, and if confirmed, they will provide a mean to improve our understanding of the initial state of the heavy-ion collisions.


I. INTRODUCTION
Heavy ion collisions produce a quark-gluon plasma (QGP) [1,2] whose space-time evolution is well described by relativistic viscous hydrodynamics [2][3][4]. The QGP expansion converts the initial-state spatial anisotropies into finalstate momentum anisotropies. These are characterized by Fourier expansion of azimuthal distribution of particle density, dN dφ ∝ 1 + 2 ∑ ∞ n=1 v n cos n(φ − Φ n ), where v n and Φ n represent the amplitude and phase of the n th -order flow vector V n = v n e inΦn . The V n reflects the hydrodynamic response of the produced medium to the n th -order initial-state eccentricity vector [5,6], denoted by E n = ε n e inΦ ε n . Due to event-by-event (EbyE) density fluctuations in the initial state, the E n and consequently the V n also fluctuate event to event. However, model calculations show that an approximate linear relation V n ∝ E n is valid for n = 2 (elliptic flow) and 3 (triangular flow) within a fixed centrality class, and the proportionality constant depends on the transport properties of the QGP [6][7][8][9][10][11].
Most previous efforts assumed that E n and V n are boost invariant within a single event. But recent studies [12,13] show significant fluctuations of harmonic flow along the longitudinal direction within the same event. This so called "flow decorrelations" appear as differences in flow magnitude (v n (η 1 ) ≠ v n (η 2 )) and its phase (Φ n (η 1 ) ≠ Φ n (η 2 )) along pseudorapidity (η). The origin can be attributed to the fact that the number of particle production sources and their transverse distribution fluctuates along η, which leads to longitudinal decorrelation of eccentricity vector in configuration space in a single event. For example, the number of forward-going and backward-going nucleon participants, N F part and N B part , are not the same in a given event [14,15], and the corresponding eccentricity vectors E F n and E B n would also be different. Since these participants contribute differently to the final-state particles in the forward and backward rapidity, the eccentricity vector is closer to E F n (E B n ) in the forward (backward) rapidity [16]. Hydrodynamic model simulations [15,[17][18][19][20][21] show that the flow decorrelations reflect mainly the longitudinal structure of the initial state, and is insensitive to the viscosity of the QGP. Therefore, flow decorrelations serve as a unique probe for the early time dynamics of the heavy-ion collisions.
The first measurement of flow decorrelations was performed by the CMS Collaboration [12], followed by a more detailed study by the ATLAS Collaboration [13] in Pb+Pb collisions. Preliminary results have also been obtained at RHIC energies as well [22]. These results were described reasonably by several hydrodynamic model simulations with a 3D initial condition based on the lund-string picture [18,19]. Very recently, ATLAS also measured the flow decorrelations in the Xe+Xe system [23]. Compared with the Pb+Pb system, the decorrelation signal is observed to be larger for v 2 , but smaller for v 3 . Current hydrodynamic models [24,25] reproduce the v n in both systems but fail to describe simultaneously the centrality dependence of the v n decorrelations, which implies that the hydrodynamic models tuned to describe the transverse dynamics may not have the correct initial-state geometry in the longitudinal direction. However, in order to pin down exactly how to improve the description of 3D initial-state geometry, further systematic measurements and model studies in different collision systems are required.
In this paper, we explore longitudinal decorrelations of the initial-stage geometry using Glauber model simulations for different collision systems. We study the qualitative trends of the system size dependence in symmetric collision system, as well as the effects of the nuclear deformation and asymmetric collision system. We found that the decorrelations are sensitive to all these variations.

II. SETUP
The longitudinal flow decorrelations are studied with a factorization ratio proposed by the CMS collaboration [12], where η r is a reference pseudorapidity range common to both the numerator and the denominator, and the average is done over events in a given centrality interval. Measurements show that r n (η) is an approximately linear function close to unity, and the slope parameter F n characterizes the strength of the decorrelation.
Since flow vector and eccentricity vector are linearly correlated, V n ∝ E n , the r n (η) can be directly related to the initial eccentricity in spatial rapidity defined analogously to Eq. 1: Hydrodynamic model calculations show that r n (η) ≈ r s n (η) [19], nearly independent of the value of shear viscosity in the final state.
Following our previous work [16], the η dependence of the eccentricity is estimated from the eccentricities of the forward-going and backward-going quark participants E F n and E B n , where f n (η) is an odd function that controls the relative mixture of the eccentricity vectors for the forward and backward going quark participants: f n (∞) = 1 and f n (−∞) = −1, and E n+ ≈ E n is the eccentricity calculated using all participants 1 . Note that the E n± fluctuate event to event but are constants within an event. Assuming f n (η) in each event is a slowly varying function near mid-rapidity, Ref. [16] shows that where a n is a constant that encodes information about the f n (η r ), and A n controls the strength of the eccentricity decorrelations.
In the linear response picture, the flow harmonics are driven by the overall eccentricity: where we use the fact that harmonic flow can only be measured via the two-particle correlation method which corresponds to ⟨v 2 n ⟩. The response coefficient κ n captures the effects of the viscous damping and depends mainly on the overall size of the system (N part or number of quark participant N qp ). With similar argument, we hypothesize that flow decorrelations should be driven by eccentricity decorrelations, The coefficients κ ′ n ≈ a n is controlled by the mixing function f n (η), whose dependences on centrality is currently unknown. Furthermore, although the influence of κ n to A n is expected to largely cancel between the numerator and denominator, some residual dependence could remain since the κ n for ε F n and ε B n can be different if N F part ≠ N B part in a given event.
In studying the system-size dependence, it is useful to consider ratios of flow harmonics or flow decorrelations as a function of The N part is a proxy for absolute system size while N part 2A can be considered as a measure for scaled system size. When plotted as a function of N part , the κ n is expected to cancel in the v n -ratio and v A+A In contrast, for the same N part 2A, the longitudinal structure of the initial state is expected to have similar F-B asymmetry in the number of sources and similar f n (η) 2 , and therefore the F n -ratio is expected to approximately scales with A n -ratio, i.e. F A+A The Eq. 6 and above arguments are the main assumptions used in this paper for our predictions of the system-size dependence of the eccentricity decorrelations.
The eccentricity and its decorrelations are calculated using a standard quark Glauber model from Ref. [26]. Three quark constituents are generated for each nucleon according to the "mod" configuration [27], which ensures that the radial distribution of the three constituents after recentering follows the proton form factor ρ proton (r) = e −r r0 with r 0 = 0.234 fm [28]. The nucleons are assumed to have a hard-core of 0.4 fm in radii, their density distribution is given by the Woods-Saxon profile, where, ρ 0 is the nucleon density, R 0 is the nuclear radius and a = 0.55 is the skin depth. The value of quark-quark cross-section is chosen to be σ qq = 18 mb, which corresponds to nucleon-nucleon inelastic cross-section σ nn = 68 mb for √ s NN = 5.02 TeV. For the study of system size dependence, six spherical nuclei are considered, see Table I. The effect of the deformation are considered for Xenon and Uranium, denoted by Xe d and U d , according to, ρ(r, θ) = ρ 0 1 + e (r−R0(1+β2Y20(θ)+β4Y40(θ))) a where Y 20 and Y 40 are Legendre polynomials and β 2 and β 4 are deformation parameters. The deformation parameters are chosen as β 2 = 0.28 and β 4 = 0.093 for U d [29] and β 2 = 0.162 and β 4 = −0.003 for Xe d [30,31].  The Glauber simulation is performed for various collision systems to generate the positions of the participant nucleons and quark constituents, which are used to calculate eccentricity ε n and eccentricity decorrelations A n . The eccentricity vector is calculated using the transverse positions of quarks as E n = − ⟨r n e inφ ⟩ ⟨r n ⟩. Similarly, the E F n and E B n are calculated using only the forward-going and backward-going quarks, respectively, which are then used to obtain the A n .

III. RESULT
The top panels of Fig. 1 show the ε 2 calculated in different collision systems. A clear hierarchy is observed when ε 2 is plotted as a function of N part . However, when plotted as a function of N part 2A, the ε 2 values for different systems nearly collapse on a common curve that simply reflects the centrality-dependent shape of the elliptic geometry of the overlap region. The bottom panels of Fig. 1 show the results for ε 3 . The ε 3 values from different systems overlap at small N part region, but deviate from each other at larger N part values. This behavior suggests that although the ε 3 is driven by the random fluctuations of quark constituents, it results in common ε 3 values only when N part is not too large. In the large N part region, the ε 3 also depends on the size and the ε 2 of the overlap region (for example due to the anti-correlation between ε 2 and ε 3 [32]). The left panels of Fig. 2 show the results of eccentricity decorrelations A 2 and A 3 as a function of N part . The A 2 values are larger for small systems, while the opposite trend is observed for the A 3 . This opposite system-size dependence trend between A 2 and A 3 is much more obvious when they are plotted as a function of N part 2A in the right panels. Recently, the ATLAS Collaboration has performed the first measurement of the system-size dependence of flow decorrelations for v 2 and v 3 [23]. We can check how well the Glauber model describes the change between Xe+Xe and Pb+Pb observed in the ATLAS data. The left panels of Fig. 3 show the N part dependence of ε n -ratios and A n -ratios, and they are compared with the v n -ratios and F n -ratios respectively from the data and hydrodynamic model predictions. The ε n -ratios agree with the v n -ratio very well and this is because the response coefficient κ n depend only on the overall size of the overlap region described by N part , and therefore cancel in the ratios. On the other hand, the A n -ratios show qualitatively similar trends as the F n -ratios, but are quantitatively different especially for n = 2. Note that the hydrodynamic model predictions reproduce the v n -ratios but fail to describe the F n -ratios, implying the model does not have the correct initial-state condition in the longitudinal direction. They are compared also with hydrodynamic model predictions (lines) for vn-ratio [31] and Fn-ratio [24,25].
The right panel of Fig. 3 show the same ratios calculated as a function of N part 2A. It is clear that ε n -ratios do not describe the v n -ratios due to the fact that the κ n do not cancel, which lead to about 10% difference between ε n -ratio and v n -ratio for n = 2 and 10-20% difference for n = 3. However, the A n -ratios, which are expected to be relatively insensitive to κ n , show an overall good agreement with the F n -ratios. This agreement implies that the f n (η) function controlling the mixing between forward-going and backward-going sources is mostly a function of centrality percentile or N part 2A between different systems. This result also supports the opposite hierarchy between the system-size dependence of A 2 and the system-size dependence of A 3 in Fig. 2.
The deformation of colliding nuclei is known to influence the N part dependence of ε n and v n [33][34][35]. An interesting question is whether the eccentricity decorrelations are also affected. Figure 4 shows the ε n -ratios and A n -ratios for Xe (top panels) and U (bottom panels) with and without deformation. In the case of Xe, the deformation influences the ε 2 and A 2 in central collisions but in the opposite direction, i.e. deformation increases the ε 2 but reduces the A 2 . The deformation has very little influences on the ε 3 and A 3 . In the case of U, the deformation increases ε 2 over a broader centrality range. The influence on A 2 is a bit non-trivial: the deformation has little effect in the mid-central and peripheral collisions, but decreases the A 2 in the ultra-central collisions. The deformation increases the values of ε 3 but decreases the values of A 3 . These features can be searched for in the experimental analyses, for example by comparing the Ru+Ru and Zr+Zr at RHIC which have the same atomic number but different amount of deformation [36]. Figure 5 shows our prediction of the eccentricity decorrelations in asymmetric collision system Cu+Au, for which the intrinsic asymmetry between the N F part and N B part should also influence the behavior of ε n and A n . Since the overall system size for Cu+Au is in between Zr+Zr and Xe+Xe, we compare the ε n and A n among these three systems. The N part dependences of A n in central Cu+Au region are distinctly different from those in the Zr+Zr and Xe+Xe systems. This is because over a wide range in the central collisions region, all the nucleons from Cu participate in the collisions, while the nucleon participants in Au still increases, resulting in a weak dependence of both ε n and A n on the N part .

IV. SUMMARY
We discussed the dependence of elliptic flow v 2 and triangular flow v 3 and their longitudinal decorrelation coefficients F 2 and F 3 on the choice of collision systems in terms of the size, deformation and asymmetry of the nuclei with atomic number A. Hydrodynamic model simulation shows that the harmonic flow are driven by the initial-state eccentricity, ε n v n ∝ ε n , and the flow decorrelations are directly determined by the eccentricity decorrelations in the longitudinal direction A n , F n ≈ A n . We estimate the values of ε n and A n in various collision systems using a Monte-Carlo quark Glauber model, which assumes three constituent quarks for each nucleon in determining the initial state, and the results are presented as a function of number of nucleon participants N part or that normalized by the total number of nucleons of the collision systems N part 2A.
We found that the A 2 is larger for smaller collision systems, while the opposite ordering is observed for the A 3 . The ratios of ε n or A n between Xe+Xe and Pb+Pb are compared with the ratios of v n or F n measured by the ATLAS Collaboration as a function of both N part and N part 2A. The ε n -ratios approximately agree with v n -ratios as a function of N part , while the A n -ratios agree with F n -ratios as a function of N part 2A. This behavior is consistent with our understanding that the flow response coefficient κ n = v n ε n depend on the overall system size described by the N part , while the coefficient for flow decorrelations F n A n might depend only on the overall shape of the overlap region controlled by the N part 2A. Current hydrodynamic models fail to describe simultanously the flow decorrelations in Xe+Xe to Pb+Pb, and this failure implies that their initial longitudinal structure based on lund-string picture of