Anisotropic flow, flow fluctuation and flow decorrelation in relativistic heavy-ion collisions: the roles of sub-nucleon structure and shear viscosity

One of the most significant discoveries in relativistic heavy-ion collisions performed at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) is the creation of the strongly-coupled Quark-Gluon Plasma (QGP). This novel hot and dense nuclear matter exhibits perfect fluid-like transport properties, as characterized by an extremely low ratio of shear viscosity to entropy density $({\eta_{v}/s})$ Adams:2003zg ; Aamodt:2010pa ; ATLAS:2011ah ; Chatrchyan:2012ta . Relativistic hydrodynamic approach Wu:2021fjf ; Pang:2018zzo ; Schenke:2010rr ; Schenke:2010nt ; Schafer:2021csj ; Karpenko:2013wva ; Shen:2014vra at high temperatures, combined with microscopic hadronic cascades SMASH:2016zqf ; Bass:1998ca ; Du:2019obx at low temperatures, have been very successful to describe the dynamical space-time evolution of heavy-ion collisions.

The strong (anisotropic) collective flow is one of the most important evidence for the formation of QGP, and has attracted significant attention over the past decades. It refers to the azimuthal anisotropy of the final produced particles, stemming from the initial state spatial anisotropy which converts to the final state momentum anisotropy as a result of the strong interaction among the QGP constituents Aguiar:2001ac ; Broniowski:2007ft ; Andrade:2008xh ; Hirano:2009ah ; Alver:2010gr ; Petersen:2010cw ; Qin:2010pf ; Staig:2010pn ; Teaney:2010vd . The final-state anisotropies are typically observed by performing the Fourier analysis of the azimuthal angle distribution of the final state particles relative to the event-plane angle $\Psi_{{n}}$, i.e., $\frac{dN}{p_{T}dp_{T}dyd\phi}\propto 1+2\sum_{{n}=1}^{\infty}v_{{n}}(p_{T},y)\cos\left[{n}\left(\varphi-\Psi_{{n}}(p_{T},y)\right)\right]$, where the Fourier coefficient $v_{n}$ is called the $n$-th order harmonic flow. Early research has mainly focused on the elliptic flow $v_{2}$ because of the underlying elliptic geometry of the collision zone. Comprehensive studies Song:2007ux ; Song:2007fn ; Schenke:2010rr ; Schenke:2010nt ; Ryu:2015vwa ; Ryu:2017qzn ; JETSCAPE:2020mzn ; JETSCAPE:2020shq have demonstrated that the magnitude of the elliptic flow $v_{2}$ is sensitive to the initial condition and transport properties of the QGP, such as specific shear viscosity ${(\eta_{v}/s)}$ and bulk viscosity ${(\zeta/s)}$. Recently, a Bayesian statistical method JETSCAPE:2020mzn ; JETSCAPE:2020shq ; Parkkila:2021yha ; Bernhard:2016tnd ; Heffernan:2023utr has been utilized to simultaneously constrain the initial conditions and transport coefficients by calibrating hydrodynamic calculations, e.g., the particle yield and flow magnitude, with the experimental data from the LHC and RHIC.

In recent years, much attention has been paid to the initial state fluctuations in relativistic heavy-ion collisions, such as the fluctuations of nucleon positions or sub-nucleon structures inside nuclei, which can not only induce the odd-order harmonics flow (e.g., $v_{3},v_{5}$) Alver:2010gr ; Alver:2010dn ; Qin:2010pf ; Teaney:2010vd , but also cause the fluctuations and correlations among flow magnitudes and flow angles from event to event in both transverse plane and longitudinal direction. Various flow observables have been defined to quantify flow fluctuations and correlations, such as event-plane correlations Aad:2014fla ; Bhalerao:2013ina ; Qin:2011uw ; Qiu:2012uy ; Magdy:2022jai , symmetric cumulants ALICE:2016kpq ; Zhu:2016puf ; Giacalone:2016afq ; Li:2021nas , non-linear response coefficients Yan:2015jma ; Qian:2016fpi ; Acharya:2017zfg ; Giacalone:2018wpp and $v_{n}$-$p_{T}$ correlation PhysRevC.93.044908 . Compared to the flow coefficients, the flow fluctuation and correlation observables can provide additional constraints on the evolution dynamics and transport properties of the QGP. For example, symmetric cumulants are sensitive to the temperature dependence of shear viscosity $\eta_{v}/s(T)$ Parkkila:2021yha ; ALICE:2016kpq ; ALICE:2017kwu , and $v_{n}$-$p_{T}$ correlation is sensitive to the initial nuclear structure Giacalone:2021clp ; Bally:2021qys ; Giacalone:2020awm ; Nijs:2020ors ; Schenke:2020uqq ; Fortier:2023xxy ; PhysRevC.93.044908 .

Another important set of flow observables is the factorization breaking (or the decorrelation effect) of the anisotropic flow Pang:2018zzo ; Pang:2015zrq ; Pang:2014pxa ; Bozek:2015bna ; ATLAS:2020sgl ; Bozek:2017qir ; Wu:2018cpc ; Sakai:2020pjw ; Zhao:2017yhj ; CMS:2015xmx . A common method for studying the factorization breaking effect is to construct two-particle correlations using flow vectors. Such method takes into account the combined contributions from flow magnitude and flow angle ATLAS:2017rij ; Bozek:2017qir ; Bozek:2021mov ; Wu:2018cpc ; Jia:2017kdq , The factorization breaking effect has been confirmed in Au+Au collisions at RHIC, as well as in Xe+Xe and Pb+Pb collisions at the LHC ATLAS:2020sgl ; CMS:2017xnj ; ATLAS:2017rij ; CMS:2015xmx ; ALICE:2022dtx . Recently, the ALICE Collaboration has developed new observables based on four-particle correlations, which can separately measure the contributions from flow magnitude and flow angle in the transverse plane ALICE:2022dtx . Similar observables have also been defined by the ATLAS Collaboration to separate the contributions from flow magnitude and flow angle along the longitudinal direction ATLAS:2017rij . Distinguishing the contributions from flow magnitude and flow angle offers a unique opportunity to understand and constrain the fundamental properties of the QGP medium.

Our paper is organized as follows. In Sec. II, we introduce the setup of the event-by-event (3+1)-dimensional CLVisc viscous hydrodynamics framework. In Sec. III, the $p_{T}$ differential flow observables are constructed based on two-particle and four-particle correlations. In Sec. IV, we present our numerical results for the elliptic flow fluctuation for identified particles, and the angle and magnitude decorrelations of elliptic flow for charged particles. The summary will be presented in Sec. IV.

In this work, we utilize the event-by-event (3+1)-dimensional CLVisc viscous hydrodynamic framework Pang:2012he ; Pang:2018zzo ; Wu:2021fjf , which includes the fluctuating TRENTO initial condition Moreland:2014oya ; Soeder:2023vdn ; Ke:2016jrd , the second-order dissipative hydrodynamic evolution and Cooper-Frye hadronic sampler, to simulate the space-time evolution of the QGP medium and hadron production in Pb+Pb collisions at $\sqrt{s_{NN}}$=5.02 TeV at the LHC.

The theoretical descriptions of the initial states of the QGP created in relativistic heavy-ion collisions from the first principle remain a significant challenge. In this study, we employ the parameterized initial condition model, TRENTO Moreland:2014oya ; Soeder:2023vdn ; Ke:2016jrd , to phenomenally simulate the initial space-time distribution of the QGP fireball. In the TRENTO initial condition, the positions of nucleons are sampled following the Woods-Saxon distributions. Compared with the traditional Glauber model Miller:2007ri , the collisional probability of nucleons is determined by the effective parton cross-sections, which are tuned to fit the inelastic nucleon cross-section. Recently, the sub-nucleon degree of freedom has been incorporated into the TRENTO model, which provides a possible approach to investigate the substructure of nucleons and extra sources of flow fluctuations Giacalone:2021clp ; JETSCAPE:2020shq ; Nijs:2021clz ; Schenke:2020unx ; Shen:2016zpp ; Loizides:2016djv ; Noronha-Hostler:2015coa ; Mantysaari:2016jaz ; Mantysaari:2016ykx . Once the positions of wounded nucleons are determined, one can obtain the thickness functions of the projectile and target nuclei as follows:

$$T_{A/B}=\sum_{i=1}^{N_{A/B}}wT_{p}\left(x,y;x_{i},y_{i}\right).$$

(1)

Here $T_{p}$ is the thickness function originating from each wounded nucleon with Gaussian form, and $w$ is a random weight taking from a Gamma distribution to account for the multiplicity fluctuations in p+p collisions. Naturally, one may define the reduced thickness function by considering the generalized average of the thickness functions of two nuclei,

$$T_{R}\left(p;T_{A},T_{B}\right)\equiv\left(\frac{T_{A}^{p}+T_{B}^{p}}{2}\right)^{1/p},$$

(2)

where the parameter $p$ characterizes the different approximations for the entropy deposition. In the new version of the TRENTO initial condition model, the sub-nucleon structure is included. The positions of constituent partons are assumed to follow a Gaussian distribution within the nucleons, and the same entropy deposition function is used as the nucleons. In this work, we choose an IP-Glasma-like entropy deposition, i.e., $p=0$ JETSCAPE:2020mzn ; JETSCAPE:2020shq . The number of constituent partons is set to 3. The Gaussian widths are set as 0.5 fm for nucleons and 0.2 fm for partons. By utilizing the scale factor $K$, the initial transverse entropy profile is assumed to be,

$$\left.\frac{dS}{dy}\right|_{\tau=\tau_{0}}=K\cdot T_{R}\left(p;T_{A},T_{B}\right),$$

(3)

where the factor $K$ can be tuned to fit the final particle yield ${dN/d{\eta}}$ in the most central collisions. Additionally, we assume the initial proper time $\tau_{0}$ as 0.6 fm. For the longitudinal direction, we use the following plateau-like envelope function,

$$H(\eta)=\exp\left[-\frac{\left(\eta-\eta_{\text{flat }}\right)^{2}}{2\eta_{{\rm gw}}{}^{2}}\theta\left(|\eta|-\eta_{\text{flat }}\right)\right].$$

(4)

In this work, we take the parameters ${\eta_{\rm flat}}=1.7$ and ${\eta_{\rm gw}}=2.0$.

After the initial entropy density distribution has been constructed, the (3+1)-dimensional CLVisc viscous hydrodynamics model is used to simulate the space-time evolution of the QGP fireball. The following energy-momentum conservation equation is solved:

$$\partial_{\mu}T^{\mu\nu}=0,$$

(5)

where $T^{\mu\nu}$ is the energy-momentum tensor of the QGP system. Since the QGP medium is not fully in local thermal equilibrium, we also solve the evolution of shear tensor and bulk pressure using the second-order Israel-Stewart equation Denicol:2018wdp ,

$$\begin{gathered}\tau_{\Pi}D\Pi+\Pi=-\zeta\theta-\delta_{\Pi\Pi}\Pi\theta+\lambda_{\Pi\pi}\pi^{\mu\nu}\sigma_{\mu\nu}\\ \tau_{\pi}\Delta_{\alpha\beta}^{\mu\nu}D\pi^{\alpha\beta}+\pi^{\mu\nu}=\eta_{v}\sigma^{\mu\nu}-\delta_{\pi\pi}\pi^{\mu\nu}\theta+\tau_{\pi\pi}\pi^{\lambda\langle\mu}\sigma_{\lambda}^{\nu\rangle}\\ +\varphi_{1}\pi_{\alpha}^{\langle\mu}\pi^{\nu\rangle\alpha}.\end{gathered}$$

(6)

Here $D$ is the covariant derivative, and $\theta=\partial_{\mu}u^{\mu}$ is the expansion rate. The traceless symmetric tensor is defined as $\pi^{\langle\mu\nu\rangle}=\Delta_{\alpha\beta}^{\mu\nu}\pi^{\alpha\beta}$ with second-order symmetric projection operator $\Delta_{\alpha\beta}^{\mu\nu}=\frac{1}{2}\left(\Delta_{\alpha}^{\mu}\Delta_{\beta}^{\nu}+\Delta_{\alpha}^{\nu}\Delta^{\mu}{}_{\beta}\right)-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$. $\tau_{\pi}$ and $\tau_{\Pi}$ are relaxation times of shear tensor and bulk pressure. $\delta_{\Pi\Pi},\lambda_{\Pi\pi},\delta_{\pi\pi},\tau_{\pi\pi},\varphi_{1}$ are second-order transport coefficients. The values of relaxation time and transport coefficients are summarized in Table 1.

In this work, the specific shear viscosity ${\eta_{v}}/s$ is taken as 0.16, the temperature dependence of bulk viscosity ${\zeta/s(T)}$ is taken from Ref. Bernhard:2016tnd , and the equation of state is taken from the HotQCD Collaboration HotQCD:2014kol . When the temperature of QGP fluid cell drops down to the freezeout temperature ($T_{\mathrm{frz}}$=150 MeV), the momentum of the thermal hadrons is obtained using the Cooper-Frye prescription,

$$\frac{dN_{i}}{dYp_{T}dp_{T}d\phi}=\frac{g_{i}}{(2\pi)^{3}}\int p^{\mu}d\Sigma_{\mu}f_{\mathrm{eq}}(1+\delta f),$$

(7)

where the thermal equilibrium term $f_{\rm{eq}}$ and out-of-equilibrium correction $\delta f$ are given as,

	$\displaystyle f_{\mathrm{eq}}$	$\displaystyle=\frac{1}{\exp\left[\left(p\cdot u-\mu_{i}\right)/T_{\mathrm{frz}}\right]\pm 1},$		(8)
	$\displaystyle\delta f$	$\displaystyle=\left(1\mp f_{\mathrm{eq}}\right)\frac{p_{\mu}p_{\nu}\pi^{\mu\nu}}{2T_{\mathrm{frz}}^{2}(\varepsilon+P)}.$

In this work, we neglect the correction from bulk pressure due to its theoretical uncertainty. After the thermal hadrons are produced, they undergo resonance decay procedures. The centrality classes are determined from the total entropy at midrapidity at initial time $\tau_{0}$ within the TRENTO initial condition. In order to improve the statistics, we run 2000 hydrodynamics events and oversample 3000 times for each hydrodynamics event at each centrality bin.

Under the assumptions that the flow vector follows a 2-dimensional Gaussian distribution and the collision geometry gives the dominant effect compared to the flow fluctuations, the two-particle and four-particle correlations can be written as follows Voloshin:2008dg ,

		$\displaystyle v_{{n}}^{2}\{2\}=\left\langle v_{{n}}\right\rangle^{2}+\sigma_{{v}_{{n}}}^{2},$		(9)
		$\displaystyle v_{{n}}^{2}\{4\}\approx\left\langle v_{{n}}\right\rangle^{2}-\sigma_{{v}_{{n}}}^{2},$

where $\left\langle v_{{n}}\right\rangle$ and $\sigma_{{v}_{{n}}}$ are the mean and fluctuation of the anisotropic flow coefficient $v_{n}$, respectively. Using the above formula, the mean flow $\left\langle v_{{n}}\right\rangle$ and flow fluctuation $\sigma_{{v}_{{n}}}$ can be estimated via $v_{{n}}\{2\}$ and $v_{{n}}\{4\}$ as:

		$\displaystyle\left\langle v_{{n}}\right\rangle\approx\sqrt{\frac{v_{{n}}^{2}\{2\}+v_{{n}}^{2}\{4\}}{2}},$		(10)
		$\displaystyle\sigma_{{v}_{{n}}}\approx\sqrt{\frac{v_{{n}}^{2}\{2\}-v_{{n}}^{2}\{4\}}{2}}.$

The relative flow fluctuation $F\left(v_{{n}}\right)$ can be calculated as ALICE:2022zks :

$$F\left(v_{{n}}\right)=\frac{\sigma_{{v}_{{n}}}}{\left\langle v_{{n}}\right\rangle}.$$

(11)

Note that the relative flow fluctuation can also be quantified by the ratio $v_{2}\{4\}/v_{2}\{2\}$.

Based on the standard multiple cumulant method Bilandzic:2010jr , the $p_{T}$ differential flow coefficient can be defined via two-particle correlation as follows,

	$\displaystyle{v_{{n}}\{2\}(p_{T})}$	$\displaystyle{=\frac{\langle\langle e^{in(\phi_{1}^{a}-\phi_{{2}}}\rangle\rangle}{\sqrt{\langle\langle e^{in(\phi_{1}-\phi_{{2}})}\rangle\rangle}}=\frac{\langle\vec{V}_{{n}}(p_{T})\vec{V}_{{n}}^{*}\rangle}{\sqrt{\langle\vec{V}_{{n}}\vec{V}_{{n}}^{*}\rangle}}}$		(12)
		$\displaystyle{=\frac{\langle v_{{n}}(p_{T})v_{{n}}\cos n(\Psi_{{n}}(p_{T})-\Psi_{{n}})\rangle}{\sqrt{\langle v_{{n}}^{2}\rangle}}},$

where the single bracket denotes the event averages, and the double angle brackets denote the average over both particles and events. $\phi_{1}^{a}$ is the azimuthal angles of final particles from some specific $p_{T}$ range, ${\Psi_{{n}}(p_{T})}$ is the corresponding $p_{T}$ differential event-plane angle, $\phi_{1}$ and $\phi_{2}$ are the azimuthal angles of final particles (without $p_{T}$ constraint), and $\Psi_{{n}}$ stands for $p_{T}$-integrated event plane angle. Due to event-by-event fluctuation, $p_{T}$ differential event plane ${\Psi_{{n}}(p_{T})}$ does not always align with $p_{T}$ integrated $\Psi_{{n}}$. Therefore, the above differential flow coefficient includes the decorrelation effect. Following the standard procedure, the $\eta$-gap technology is utilized in this work to eliminate the short-range correlation.

One may also define the $p_{T}$ differential flow coefficient using the four-particle correlation method as follows,

$$v_{n}\{4\}(p_{T})=-\frac{d_{n}\{4\}(p_{T})}{\left(-c_{n}\{4\}\right)^{3/4}},$$

(13)

where $d_{n}\{4\}(p_{T})$ is related to the fourth-order $p_{T}$ differential cumulant and $c_{n}\{4\}$ corresponds to the fourth-order $p_{T}$ integrated cumulant, as described in Ref. Bilandzic:2010jr .

In order to separate the flow angle and flow magnitude fluctuations in the ${p_{T}}$-dependent flow vector, ALICE collaboration ALICE:2022dtx has recently proposed a new set of four-particle correlation functions,

	$\displaystyle{A_{{n}}^{f}}$	$\displaystyle{=\frac{\langle\langle\cos n[\phi_{1}^{a}+\phi_{{2}}^{a}-\phi_{3}-\phi_{4}]\rangle\rangle}{\langle\langle\cos n[\phi_{1}^{a}+\phi_{{2}}-\phi_{3}^{a}-\phi_{4}]\rangle\rangle}}$		(14)
		$\displaystyle{=\frac{\langle v_{{n}}^{2}(p_{T}^{a})v_{{n}}^{2}\cos 2n[\Psi_{{n}}(p_{T}^{a})-\Psi_{{n}}]\rangle}{\langle v_{{n}}^{2}(p_{T}^{a})v_{{n}}^{2}\rangle}}$
		$\displaystyle{\approx\langle\cos 2n[\Psi_{{n}}(p_{T}^{a})-\Psi_{{n}}]\rangle},$

and

	$\displaystyle{M_{{n}}^{f}}$	$\displaystyle{=\frac{\langle\langle\cos n[\phi_{1}^{a}+\phi_{{2}}-\phi_{3}^{a}-\phi_{4}]\rangle\rangle}{(\langle\langle\cos n[\phi_{1}^{a}-\phi_{3}^{a}]\rangle\rangle\langle\langle\cos n[\phi_{{2}}-\phi_{4}]\rangle\rangle)}}$		(15)
		$\displaystyle{\qquad\bigg{/}\frac{\langle\langle\cos n[\phi_{1}+\phi_{{2}}-\phi_{3}-\phi_{4}]\rangle\rangle}{\langle\langle\cos n[\phi_{1}-\phi_{{2}}]\rangle\rangle^{2}}}$
		$\displaystyle{=\frac{\langle v_{{n}}^{2}(p_{T}^{a})v_{{n}}^{2}\rangle/(\langle v_{{n}}^{2}(p_{T}^{a})\rangle\langle v_{{n}}^{2}\rangle)}{\langle v_{{n}}^{4}\rangle/\langle v_{{n}}^{2}\rangle^{2}}},$

where $A_{n}^{f}$ ($M_{n}^{f}$) denotes the decorrelation between $p_{T}$ differential flow angle (magnitude) and $p_{T}$ integrated flow. If $p_{T}$ differential flow is the same as $p_{T}$ integrated flow, then the $A_{n}^{f}$ and $M_{n}^{f}$ should be the unity. It should be noted that the flow angle decorrelation $A_{n}^{f}$ should be less than or equal to unity due to the cosine function in the definition. However, $M_{n}^{f}$ does not have such upper limit, due to the fact that it not only involves the covariance between the $p_{T}$ differential flow and the $p_{T}$ integrated flow, but also their relative fluctuations.

In this section, we present the numerical results for the $p_{T}$ differential flow coefficients, flow fluctuation, flow angle and magnitude decorrelations in Pb+Pb collisions at $\sqrt{s_{NN}}$=5.02 TeV. These results are based on the (3+1)-dimensional CLVisc hydrodynamics framework with fluctuating TRENTO and AMPT initial conditions Lin:2004en . We first focus on the $p_{T}$ differential flow and flow fluctuations of identified particles in various centrality intervals and study the effects of shear viscosity. Later, we discuss the granularity of initial state fluctuations by considering the nucleon and sub-nucleon structures in the initial condition, respectively. Finally, we demonstrate the hydrodynamic description of flow angle and flow magnitude fluctuations, and their dependence on initial condition model, shear viscosity and sub-nucleon structure.

In this work, we have performed a systematic study on elliptic flow coefficients, flow fluctuation and flow decorrelation in Pb+Pb collisions at $\sqrt{s_{NN}}$=5.02 TeV based on the (3+1)-Dimensional CLVisc hydrodynamics framework with the TRENTO (and AMPT) initial conditions. In particular, we have studied the effect of shear viscosity and sub-nucleon structure on these flow observables.

Our numerical results with shear viscosity $\eta_{v}/s=0.16$ can reasonably describe the experimental data on two-particle $v_{2}\{2\}$ and four-particle $v_{2}\{4\}$, mean flow ${\langle v_{2}\rangle}$ and flow fluctuation ${\sigma_{v_{2}}}$ for $\pi^{\pm}$ and $K^{\pm}$ particles. For protons, our result overestimates the elliptic flow and flow fluctuation due to the absence of hadronic afterburner in the simulation. Shear viscosity tends to suppress the elliptic flow coefficient and flow fluctuation, due to its smearing effect on the local density fluctuation. In addition, the flow coefficients appear to be insensitive to the sub-nucleon structure, while the flow fluctuation does depend on the sub-nucleon structure: it tends to be suppressed by the sub-nucleon structure in central collisions but enhanced in peripheral collisions. Our result including the sub-nucleon structure can better describe the experimental data on the relative flow fluctuation observables, such as ${v_{2}\{4\}/v_{2}\{2\}}$ and $F({v_{2}})$, measured by the ALICE Collaboration at the LHC.

We have further studied the effects of the shear viscosity, the sub-nucleon structure and the initial condition on flow angle and flow magnitude decorrelations (${A_{2}^{f}}$, ${M_{2}^{f}}$) using four-particle correlation method as proposed by the ALICE Collaboration. We found that the decorrelation effect for both flow angle and flow magnitudes is largest in most central collisions. Also, the flow angle decorrelation effect increases with the transverse momentum. In addition, the flow angle decorrelation is insensitive to shear viscosity and only slightly depends on initial condition in most central collisions, whereas the flow magnitude decorrelation shows sensitivity to both shear viscosity and initial condition. In particular, the AMPT initial condition gives very different behavior for the flow magnitude decorrelation as compared to the TRENTO initial condition.

In the future, it would be interesting to include the hadronic afterburner in the simulation, which will allow us to investigate additional sources of flow fluctuation and decorrelation. Additionally, performing a Bayesian analysis based on calibrations including flow fluctuation and decorrelation observables can further help to constrain the transport properties and initial structure of the QGP produced in relativistic heavy-ion collisions.

This work is supported in part by Natural Science Foundation of China (NSFC) under Grant Nos. 12225503, 11890710, 11890711 and 11935007. Some of the calculations were performed in the Nuclear Science Computing Center at Central China Normal University (NSC³), Wuhan, Hubei, China.