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 nth-order flow vector V n = v n e in n . The V n reflects the hydrodynamic response of the produced medium to the nth-order initial-state eccentricity vector [5,6], denoted by E n = ε n e in ε n . Due to event-by-event 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. These so-called Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

"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 the eccentricity vector in configuration space. For example, the number of forward-going and backwardgoing 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 particle multiplicity in the forward (backward) rapidity is more correlated with N F part (N B part ) (see Ref. [16] for a recent detailed study on this), and similarly the eccentricity vector in the forward (backward) rapidity is closer to E F n (E B n ) [17]. Hydrodynamic model simulations [15,[18][19][20][21][22] show that the flow decorrelations reflect mainly the longitudinal structure of the initial state, and are 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 the Relativistic Heavy-Ion Collider (RHIC) energies as well [23]. These results were described reasonably by several hydrodynamic model simulations with a three-dimensional (3D) initial condition based on the Lund-string picture [19,20]. Very recently, ATLAS also measured the flow decorrelations in the Xe+Xe system [24]. 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 [25,26] 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 a 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 the 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 (η) [20], nearly independent of the value of shear viscosity in the final state.

Following our previous work [17], the η dependence of the eccentricity is estimated from the eccentricities of the forward-going and backward-going quark participants E F n and

where f n (η) is an odd function that controls the relative mixture of the eccentricity vectors for the forward-and backwardgoing quark participants:

, and E n+ ≈ E n is the eccentricity calculated using all participants. 1Note 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 midrapidity, Ref. [17] 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 participants N qp ). With a similar argument, we hypothesize that flow decorrelations should be driven by eccentricity decorrelations,

The coefficient κ ′ n ≈ a n is controlled by the mixing function f n (η), whose dependence 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 N part or N part /2A where A is the atomic number,

The N part is a proxy for absolute system size while N part /2A can be considered as a measure for scaled system size or centrality. When plotted as a function of N part , the κ n is expected to cancel in the v n ratio and

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 scale with the A n ratio, i.e., (6) and the 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. [27]. Three quark constituents are generated for each nucleon according to the "mod" configuration [28], which ensures that the radial mass locations for forward-going and backward-going nucleons (see Ref. [15]). 2 In the limit of many sources per nucleon or optical Glauber, the F-B asymmetry should be a universal function of N part /2A. distribution of the three constituents after recentering follows the proton form factor ρ proton (r) = e -r/r 0 with r 0 = 0.234 fm [29]. 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 fm is the skin depth. The value of the quark-quark cross section is chosen to be σ qq = 18 mb, which corresponds to the nucleon-nucleon inelastic cross-section σ nn = 68 mb at √ s NN = 5.02 TeV. For the study of system size dependence, six spherical nuclei are considered (see Table I). The effects of the deformation are considered for xenon and uranium, denoted by Xe d and U d , according to

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 [30] and β 2 = 0.162 and β 4 = -0.003 for Xe d [31,32]. 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 the 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 (e.g., due to the anticorrelation between ε 2 and ε 3 [33]).

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 [24]. 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, . They are compared also with the hydrodynamic model predictions (lines) for the v n ratio [32] and the F n ratio [25,26].

from the data and hydrodynamic model predictions. The ε n ratios agree with the v n ratio very well (within 5%-10%), and this is because the response coefficient κ n depends only on the overall size of the overlap region described by N part , and therefore cancels 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.

The right panel of Fig. 3 shows 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 leads to about a 10%-15% difference between the ε n ratio and the v n ratio for n = 2 and a 10%-25% 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 the centrality percentile or the N part /2A between different systems. This result also supports the opposite hierarchy between the systemsize 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 [34][35][36]. 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 defor- mation. 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 influence 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 nontrivial: the deformation has little effect in the midcentral and peripheral collisions, but decreases the A 2 in the ultracentral 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 amounts of deformation [37].

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 dependencies of A n in the 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 increase, 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 is 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 the 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 the v n ratios as a function of N part , while the A n ratios agree with the 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 depends 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 simultaneously the flow decorrelations in Xe+Xe to Pb+Pb. This failure implies that future investigations are needed on the shortcomings of the initial longitudinal structure based on the Lund-string picture of a multi-phase transport (AMPT) model, as well as on the role of final-state longitudinal expansion for the flow decorrelations.

We further compared the ε n and A n for Xe and U nuclei with and without the effects of nuclear deformation. For the modest deformation parameter β 2 = 0.162 of Xe, the deformation influences the ε 2 and A 2 only in ultracentral collisions. For the large deformation parameter β 2 = 0.28 of U, the nuclear deformation influences both the n = 2 and n = 3 of ε n and A n over a broad centrality range. The deformation always increases the value of ε n while it decreases the value of A n . These features might be searchable using the Zr+Zr and Ru+Ru isobar data from the STAR Collaboration since Zr and Ru are expected to have slightly different β 2 values [37]. We also considered the Cu+Au asymmetric collision system, which shows a N part dependence of ε n and A n different from symmetric systems in the central region.

This paper serves as an exploratory study of the possible influence of various nuclear geometry effects on the harmonic flow and longitudinal flow decorrelation in heavy-ion collisions. More quantitative predictions would require coupling the 3D initial condition with state-of-the-art hydrodynamic model simulation. The final-state effects are expected to significantly change the relation between the v n and ε n in a system-dependent manner, but they may not change too much the relationship between the F n and A n as predicted by our calculations. In addition, it would be important to extend this study to small collision systems such as p+A and pp collisions, where the eccentricity decorrelations should be sensitive to the nature of subnucleonic fluctuations, which is not yet modeled realistically in our framework.