Anisotropic flow fluctuations in hydro-inspired freeze-out model for relativistic heavy ion collisions

The LHC data on event-by-event harmonic flow coefficients measured in PbPb collisions at center-of-mass energy 2.76 TeV per nucleon pair are analyzed and interpreted within the HYDJET++ model. To compare the model results with the experimental data the unfolding procedure is employed. The essentially dynamical origin of the flow fluctuations in hydro-inspired freeze-out approach has been established. It is shown that the simple modification of the model via introducing the distribution over spatial anisotropy parameters permits HYDJET++ to reproduce both elliptic and triangular flow fluctuations and related to it eccentricity fluctuations of the initial state at the LHC energy.


Introduction
Azimuthal anisotropy of multi-particle production in relativistic heavy ion collisions is a powerful probe of collective properties of sub-nuclear matter created at extremely high densities and temperatures (see, e.g., recent reviews [1,2] and references therein). It is commonly described by the Fourier decomposition of the invariant cross section in the form (1) where p T is the transverse momentum, η is the pseudorapidity, ϕ is the azimuthal angle with respect to the reaction plane Ψ n , and v n are the Fourier coefficients. The observation of strong elliptic flow, which is the second harmonic, v 2 , in heavy ion collisions at RHIC, was argued as one of the main pieces of evidence for strongly interacting partonic matter a e-mail: igor@lav01.sinp.msu.ru ("quark-gluon fluid") formation [3][4][5][6]. At the LHC, a number of interesting measurements involving momentum and centrality dependencies of second-and higher-order harmonic coefficients in PbPb collisions at √ s NN = 2.76 TeV have been done by ALICE [7][8][9], ATLAS [10][11][12][13][14][15] and CMS [16][17][18][19][20] Collaborations. In particular, the event-by-event (EbyE) distributions of second, third, and fourth harmonics of the anisotropic flow have been obtained [12]. Other important observations are the azimuthal anisotropy of jet [21] and charmed meson [22,23] yields in PbPb collisions, and elliptic v 2 and triangular v 3 flow of inclusive [24][25][26] and identified [27] hadrons in pPb collisions.
In our previous study [28] the second-and higher-order harmonics of inclusive and identified charged hadrons in PbPb collisions at √ s NN = 2.76 TeV were analyzed in the framework of HYDJET++ model [29]. It has been shown that the cross-talk of elliptic v 2 and triangular v 3 flow in the model generates both even and odd harmonics of higher order. This mechanism is able to reproduce the p T and centrality dependencies of quadrangular flow v 4 , and also the basic trends for pentagonal v 5 and hexagonal v 6 flows. Moreover, it reproduces also specific angular dihadron correlations including the so-called "ridge effect" [30]. However, here we restricted ourselves to the analysis of the event-averaged harmonics v n ( p T ). In recent years, the study of anisotropic flow fluctuations has attracted much interest because of their direct connection with the geometry of the initial state of a relativistic heavy ion collision [31][32][33][34][35][36][37][38][39][40][41][42]. In the present paper, therefore, we analyze the EbyE distributions of the flow coefficients in PbPb collisions at the LHC within the HYDJET++ model.
The paper is organized as follows. The flow fluctuations intrinsic to the HYDJET++ are discussed in Sect. 2. Here the probability densities of both longitudinal and transverse flow components, as well as the flow modulus, obtained at different collision centralities are shown to be nicely fitted to Gaussian. The fluctuations can be enhanced by the EbyE Gaussian smearing of the spatial anisotropy parameters of the model. Section 3 describes the unfolding procedure proposed by the ATLAS Collaboration to get rid of the non-flow fluctuations. This procedure is utilized in Sect. 4 in HYDJET++ calculations to compare the model results with the experimental data on the same footing. The agreement with the data on fluctuations of both elliptic and triangular flow is quite good. Conclusions are drawn in Sect. 5.

Inherent flow fluctuations and eccentricity fluctuations in HYDJET++ model
The event generator HYDJET++ (the successor of HYD-JET [43]) is the Monte-Carlo model of relativistic heavy ion collisions, which incorporates two independent components: the soft hydro-type state with preset freeze-out conditions, and the hard state resulting from the in-medium multiparton fragmentation and taking into account the jet quenching effect. The details of this model can be found in the HYDJET++ manual [29]. Its input parameters have been tuned to reproduce the experimental LHC data on various physical observables measured in PbPb collisions [28,44], namely, centrality and pseudorapidity dependence of inclusive charged particle multiplicity, transverse momentum spectra, and π ± π ± correlation radii in central PbPb collisions, and momentum and centrality dependencies of elliptic and higher-order harmonic coefficients. In order to simulate higher azimuthal anisotropy harmonics, the following simple modification [28,45] has been implemented in the model. In the original HYDJET++ version [29] the direction and strength of the elliptic flow are governed by two parameters. The spatial anisotropy (b) represents the elliptic modulation of the final freeze-out hypersurface at a given impact parameter b, whereas the momentum anisotropy δ(b) deals with the modulation of flow velocity profile. Both δ(b) and (b) can be treated independently for each centrality, or (basic option of the model) can be related to each other through the dependence of the elliptic flow coefficient v 2 ( , δ) obtained in the hydrodynamical approach [46]: Then, due to the proportionality of v 2 (b) to the initial ellipticity 0 (b) = b/2R A , where R A is the nuclear radius, the relation between δ(b) and (b) takes the form [29] where the two parameters C and k are independent on centrality and should be obtained from the fit to the data. Compared to the former transverse radius of the fireball, which reproduces the elliptic deformation the altered radius of the freeze-out hyper-surface in the azimuthal plane takes into account triangular deformation as well: Here ϕ is the spatial azimuthal angle of the fluid element relatively to the direction of the impact parameter. R f 0 is the model parameter which determines the scale of the fireball transverse size at freeze-out, and the new parameter 3 (b) is responsible for the triangular spatial anisotropy. The event plane of the triangular flow, Ψ 3 , is randomly oriented with respect to the plane Ψ 2 , which is fixed to zero in the model calculations. This means that the elliptic and triangular flows are generated independently, in accordance with the experimental observations. Higher flow harmonics are not explicitly generated in the model, therefore these harmonics are absent if both v 2 and v 3 are absent.
It should be noted that although the azimuthal anisotropy parameters (b), δ(b), and 3 (b) are fixed at given impact parameter b, they define v n (b) only after the averaging over many events due to the inherent model fluctuations. The main source of the flow fluctuations in HYDJET++ is fluctuations of particle momenta and multiplicity. Recall that the momentum-coordinate correlations in HYDJET++ for soft component are governed by collective velocities of the fluid elements, and so the fluctuations in particle coordinates are reflected in their momenta. The fluctuations became stronger as resonance decays and (mini-)jet production are taken into account. An event distribution over the collision impact parameter for each centrality class also increases such fluctuations.
The detailed study of the EbyE flow fluctuations is the subject of our present investigation. The possible further modification of HYDJET++ to match the experimental data on the flow fluctuations would be smearing of all three parameters, , δ, and 3 , at a given b.
To get some notion for the inherent model fluctuations, we start with HYDJET++ simulations for the simplest case of central PbPb collisions (b = 0) in which all azimuthal anisotropy parameters (b), δ(b), and 3 (b) are equal to zero. Figure 1 shows the probability densities both for each component of the flow vector V n and for its modulus V n = |V n |, n = 2, 3, 4. 1 In this "fluctuation-only" scenario the proba- bility densities of V n are well described by two-dimensional (2D) Gaussian functions [47,48]: whereas the probability densities of V n have the forms of one-dimensional (1D) Gaussians: which are obtained from Eq. (6) by integration over the azimuthal angle. These distributions are characterized by a single parameter σ n only, which regulates both the modulus mean V n and the Footnote 1 continued components obtained by the averaging over all particles in an event and over all events in the data sample.
In HYDJET++ the value of this fitting parameter, the Gaussian width σ n , appears to be unique for all harmonics with a good enough accuracy: σ 2 σ 3 σ 4 ≈ 0.013. It depends on the number of model parameters which were already fixed. The main regulator of σ n is the mean multiplicity, which determines the variation of σ n (b) with centrality.
Recall that in all Bjorken-like models with cylindrical parameterization the azimuthal anisotropy of the freezeout surface transforms into the azimuthal anisotropy of particle momentum distribution proportionally to a term [29], arising in the scalar product of 4-vectors of particle momentum and flow velocity of the fluid element. Here φ is the azimuthal angle of the fluid element, ϕ is the particle azimuth, T is the freeze-out temperature and Y T is the transverse flow rapidity, respectively. The pre-factor before the cosine controls the azimuthal angle structure of particle spectrum and its inverse characterizes the fluctuation width squared. We have also verified numerically that at fixed mean multiplicity in selected p T window the width of σ n ( p T ) is approximately proportional to the factor is the maximal transverse flow rapidity.
The "true" direction of the flow vector V n in HYDJET++ for any azimuthal harmonic is pre-defined in each event. Therefore, V nL can be calculated as a longitudinal component  (6) and (7) we get here [47,48] where I 0 is the modified Bessel function of the first kind of zeroth order. Both the width σ V n = V 2 n − V n 2 and the modulus mean V n are controlled by the true value of the flow vector V nL and the width σ n (b), but they cannot be cast analytically as functions of V nL and σ n (b). Note also that V n is not equal to V nL exactly, and the azimuthal anisotropy parameters (b), δ(b), and 3 (b) have been tuned earlier at a given impact parameter b in such a way that the value of V nL extracted from the distribution (12) reproduces just the experimentally observed value of v n (b) entering in Eq. (1). Similar to the non-flow case, presented in Fig. 1, the widths σ n of the longitudinal and transverse distributions shown in Fig. 2 are approximately the same, σ 2 σ 3 σ 4 ≈ 0.02, but the distributions become broader. Also, the maxima of p(V nL ) and p(V n ) distributions are shifted toward zero with rising harmonic number n, indicating that v 2 > v 3 > v 4 at this centrality. Surely, it is interesting to compare our inherent model probability densities, obtained without any additional special parameters for the azimuthal fluctuations, with the experimental data.
We have also considered including of additional "eccentricity" fluctuations in HYDJET++ model. The simplest modification for this purpose is to introduce EbyE Gaussian smearing of the spatial anisotropy parameters (b) and ing experimental data is impossible. The EbyE distributions of anisotropic flow harmonics have been obtained by the ATLAS Collaboration [12] for the distribution of the modulus of the flow vector by application of the so-called "unfolding procedure". The goal of the unfolding procedure was to extract the "true" flow vector from the observed one by excluding the influence of non-flow effects, such as resonance decays and jet fragmentation, as well as the finite event multiplicity effect. Therefore, in the following we will show our results before and after the unfolding to be adequate. In order to employ the EbyE unfolding procedure for simulated events, the analysis method from [12] was utilized.
-The EbyE distributions of charged particles in PbPb collisions at √ s NN = 2.76 TeV with p T > 0.5 GeV/c and |η| < 2.5 are used as input distributions. Fourier decomposition of the azimuthal distribution is rewritten as where V obs n is the magnitude of the observed per-particle flow vector and Ψ obs n is the event plane angle. -The single-particle EbyE distributions are constructed and used in our unfolding procedure: V obs n,y = V obs n sin nΨ obs n = sin nϕ .
The averaging in the last two equations of (15) is performed over all hadrons in a single event. As was shown in [12], the distribution obtained after the unfolding procedure did not depend on the method applied to obtain V obs n . Thus, the single-particle method can be used for our study.
-The response function is constructed using the "two subevent method" (2SE), namely, the charged particles are divided into two sub-events with η < 0 and η > 0. The smearing effects are estimated by the difference of the flow vectors between the two sub-events, for which the flow signal cancels. This distribution is fitted to the Gaussian with the width δ 2SE determined mainly by the finite multiplicity effect and the non-flow contributions. It was shown in [12] that the response function constructed by such a procedure can be expressed as Here δ = δ 2SE /2 because we use the full-event V obs n distribution as an input, and δ 2SE is the width obtained from the difference between the EbyE per-particle flow vectors of the two sub-events.
-The constructed response function is used to obtain the unfolding matrix: where A ji is the response function between e j = V obs n ("effect") and c i = V n ("cause"). The true V n distribution (cause "c") is obtained from the measured V obs n distribution (effect "e") using an iterative algorithm. Then the Bayesian unfolding procedure is performed by means of the RooUnfold package [49].
The difference between the δ 2SE and σ n arises mostly due to the dynamical flow fluctuations. Therefore, the EbyE unfolding analysis excludes the effects related to δ 2SE , and leaves the genuine flow fluctuations.
The application of Bayesian unfolding method for the anisotropic flow analysis was checked with heavy ion event generators in [50]. It was shown that the restored density distributions were able to reproduce the input V n distributions. The non-flow effects were estimated by using the HIJING event generator [51] with and without the implementation of the flow signal, and were found to be of the order of statistical errors. In the case of the AMPT model [52], which includes both flow and non-flow fluctuations, the difference between the "generated" and "unfolded" flow harmonics in semiperipheral Au+Au collisions at RHIC was found to be small for elliptic flow and significantly increasing in the tails for triangular and quadrangular flows. Fluctuations originating from the finite multiplicity effect can be evaluated under the assumption of Gaussian multiplicity distribution [50]: where N and σ N are the mean value and the width of the distribution, respectively.  Fig. 3 The correlation between the event-averaged elliptic flow v 2 and triangular flow v 3 of charged hadrons at transverse momentum 0.5 < p T < 2 GeV/c and pseudorapidity |η| < 2.5 for ten 5 %centrality intervals in the centrality range 0-50 % of PbPb collisions at √ s NN = 2.76 TeV. Closed circles denote ATLAS data from [15], asterisks represent HYDJET++ events. The line is drawn to guide the eye event-averaged elliptic and triangular flow coefficients. The results of the model simulations are plotted onto the ATLAS data [15] in Fig. 3. One can see that the calculations and the data agree well within the 7 % accuracy limit. Now we consider the anisotropic flow fluctuations. The simulations and analysis were performed for three centrality intervals, namely, 5-10, 20-25, and 35-40 %, and for two settings of the HYDJET++ model: (i) without and (ii) with the additional smearing of spatial anisotropy parameters (b) and 3 (b); see Sect. 2. The results for p(V 2 ) and p(V 3 ) are presented in Figs. 4 and 5, respectively, and listed in Table 1.

Comparison of HYDJET++ simulations with LHC data
Both figures indicate that the original version of HYDJET++ [without the smearing of (b) and 3 (b)] already includes some dynamical fluctuations due to radial flow. Moreover, the mean values of V n and widths of the distributions σ V n in the default version of HYDJET++ are quite close to the measured ones in collisions with centralities up to 25 % for the triangular flow and up to 45 % for the elliptic flow; see Table 1. Although the agreement can be further improved by rescaling of V n to match the data, such approach would be completely misleading. The implementation of the unfolding procedure clearly demonstrates in Figs. 4 and 5 that the initial distributions become more narrow. Thus, the intrinsic fluctuations appear to be too weak to match the experimental data [12]. On the other hand, simple modification of the model via introducing the normal distribution over the spatial anisotropy parameters allows us to reproduce the measured EbyE fluctuations of elliptic and triangular flow, the distribution widths and the event-averaged values of V 2 and V 3 . The distributions p(V n ) obtained with the set of the smeared out parameters are broader in the tails compared with the experiment. Unfolding makes it narrower. Here the difference between the "initial" and "unfolded" spectra are not so dramatic although still noticeable in contrast to that at RHIC energy, calculated in [50] within the AMPT model. The two additional parameters of the model appearing in this case are the coefficients of proportionality between the Gaussian widths of the distributions p( ) and p( 3 ) and their "unsmeared" values. These two coefficients are fixed to fit the  The same as Fig. 4 but for the triangular flow V 3 in three centrality intervals data on p(V 2 ) and p(V 3 ), respectively, for only one arbitrary centrality, whereas for other centralities they are the same. It is worth noting that such a simple modification of the model also increases EbyE fluctuations for higher-order harmonics v n (n > 3), which arise in HYDJET++ due to the presence of elliptic v 2 and triangular v 3 flows, and its interference. However, significant sensitivity of high harmonic values on their extraction methods makes the direct comparison of our simulations with the data even more tricky than for v 2 and v 3 . For example, the centrality dependence of quadrangular flow v 4 measured by event plane and twoparticle cumulant methods is significantly weaker than that of v 4 measured by Lee-Yang zero method due to large nonflow contribution and increase of the flow fluctuations in more central events. Since HYDJET++ was tuned to fit the p T -dependence of v 4 {LY Z}, it underestimates v 4 extracted by the event plane or two-particle cumulant methods in (semi)central collisions [28]. We plan to study the EbyE fluctuations of higher-order flow harmonics in the future.
A few important issues should be clarified still. We utilized normal smearing of the parameters (b), δ(b), and 3 (b) in the modified version of the event generator. What will happen if the parameters are smeared out with respect to a non-Gaussian distribution? To check this possibility we opted for a uniform distribution of the key parameters within the interval ±σ to ensure the same mean and width values. The distribution p(V 2 ) is displayed in Fig. 6 for centrality bin 20-25 %, where the signal heavily dominates over the fluctuations. Interestingly enough, the generated distribution is very close to the ATLAS unfolded curve everywhere, but in the low-V 2 range. The unfolded HYDJET++ spectrum, how-ever, is narrow. It resembles the implemented rectangularshaped V 2 -distribution with some rounding of the shoulders because of the intrinsic model fluctuations.
Next question is the Gaussian-like behavior of the obtained spectra after the unfolding. ATLAS Collaboration reported some deviations from the Gaussians observed for p(V 2 ) distributions in peripheral collisions [15]. The last centrality bin in ATLAS analysis is 60-65 %. At this centrality the default version of HYDJET++, which works reasonably well up to 40-45 % [29], needs further fine tuning in line with other semi-phenomenological models. Simply, the linear dependence for v 2 (b) becomes too crude here. The new tuned values of the parameters are = 0.14 and δ = 0.25 (cf. with = 0.16 and δ = 0.38 in the default version). It is worth noting that this is the only modification of the model, whereas the ratio /σ is kept constant for all centralities in question. Figure 7 shows the observed and the unfolded distributions of p(V 2 ) obtained in HYDJET++ for centrality 60-65 % in comparison with the ATLAS data. One can see that the model calculations agree well with the data. The unfolded distribution provided by HYDJET++ was also fitted to complex Bessel-Gauss (BG) product given by Eq. (13). Results are plotted onto the simulated spectra in Fig. 7 as well. At this centrality the BG fit clearly deviates from the data, indicating that the model possesses some intrinsic fluctuations which cause the distortion of the initial Gaussians. This interesting problem definitely deserves further investigation.
Finally, our results can be used for the analysis of fluctuations of the initial anisotropy ε n . This approach relies on the assumption of a linear response of the flow coefficient V n to the corresponding initial eccentricity, Table 1 Mean values V n , widths σ Vn , and ratios σ Vn / V n (n = 2, 3) for HYDJET++ simulations and ATLAS data [12] Centrality where k n is the response coefficient. For both elliptic and triangular flow this assumption works very well, as was confirmed by hydrodynamic model calculations [53,54]. The probability distribution of the flow is connected to the initial anisotropy distribution via [55] p(V n ) = dε n dV n p(ε n ).
Inserting Eq. (19) into Eq. (20) we get Then, in [37], the elliptic power distribution was proposed to parametrize the eccentricity distributions, Here the parameter ε 0 is approximately the mean reaction plane eccentricity and α describes the eccentricity fluctuations. In [55] the authors fitted the ATLAS data on the elliptic flow to Eq. (22). It is very tempting, therefore, to  fit the HYDJET++ generated distributions to Eq. (22) and compare the extracted parameters, α, ε 0 , and k 2 . The fitted curves are plotted onto the model calculations (with smearing and unfolding procedure) of p(ε 2 ) distributions for two centralities, 20-25 and 35-40 %, in Fig. 8. The extracted fit parameters are compared in Table 2 with those obtained in [55]. The agreement between the two sets is good, indicating that HYDJET++ quantitatively reproduces the anisotropic flow fluctuations.

Conclusions
The phenomenological analysis of the EbyE distributions of anisotropic flow harmonics measured in lead-lead collisions at the center-of-mass energy 2.76 TeV per nucleon pair has been performed within the two-component HYDJET++ model. The unfolding procedure was applied to examine the simulated events, thus allowing for the direct comparison of model calculations with the experimental data. To the best of our knowledge, this is for the first time that the model-generated spectra are filtered by means of the unfolding procedure and then compared directly with the LHC data obtained by the same method. This procedure removes the non-flow effects, originating, e.g., from the decays of resonances and fragmentation of jets, as well as the finite event multiplicity effect. The essentially dynamical origin of the flow fluctuations in hydro-inspired freeze-out approach has been established. The effect is traced to the correlation between the momenta and coordinates of final particles and the velocities of hadronic fluid elements. The simple modification of the model via introducing the distribution over spatial anisotropy parameters permits HYDJET++ to reproduce both elliptic and triangular flow fluctuations in heavy ion collisions at the LHC energy. In contrast, an attempt to utilize the uniform non-Gaussian smearing with the same /σ ratio failed shortly. The unfolding procedure is sensitive, therefore, to the initial distributions of the parameters. It should be implemented in any model in the case of comparison with the unfolded data.
For the peripheral topologies the model calculations deviate from the Bessel-Gauss fit to Eq. (13) thus hinting to some intrinsic fluctuations in HYDJET++ which cause the distortion of initial Gaussians. This interesting problem deserves to be studied in the future.
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .