Higher harmonics of azimuthal anisotropy in relativistic heavy ion collisions in HYDJET++ model

The LHC data on azimuthal anisotropy harmonics from PbPb collisions at center-of-mass energy 2.76 TeV per nucleon pair are analyzed and interpreted in the framework of the HYDJET++ model. The cross-talk of elliptic $v_2$ and triangular $v_3$ flow in the model generates both even and odd harmonics of higher order. Comparison with the experimental data shows that this mechanism is able to reproduce the $p_{\rm T}$ and centrality dependencies of quadrangular flow $v_4$, and also the basic trends for pentagonal $v_5$ and hexagonal $v_6$ flows.


Introduction
The study of the fundamental theory of strong interactions (Quantum Chromodynamics, QCD) in the regimes of extreme densities and temperatures is ongoing via the measurement of the properties of hot and dense multiparton systems produced in high-energy nuclear collisions (see, e.g., reviews [1,2,3,4]). The started LHC heavy ion program makes it possible to probe the new frontiers of the high temperature QCD providing the valuable information on the dynamical behavior of quark-gluon matter (QGM), as predicted by lattice calculations. A number of interesting LHC results from PbPb runs at √ s NN = 2.76 TeV, have been published by ALICE, ATLAS and CMS collaborations (see [5] for the overview of the results from the first year of heavy ion physics at LHC).
One of the modern trends in heavy ion physics at high energies is a study of Fourier harmonics of azimuthal particle distribution, which is a powerful probe of bulk properties of the created high density matter. It is typically described by a Fourier series of the form: where ϕ is the azimuthal angle with respect to the reaction plane Ψ n , and v n are the Fourier coefficients. The second harmonic, v 2 , referred to as "elliptic flow", is the most extensively studied one, because it directly relates the anisotropic shape of the overlap of the colliding nuclei to the corresponding anisotropy of the outgoing momentum distribution. The momentum and centrality dependencies of the elliptic flow in PbPb collisions were measured at the LHC [6,7,8] in the first instance. Then, the results of measurements of the higher azimuthal harmonics [9,10,11] and the anisotropic flow of identified particles [12] were published. The higher order coefficients v n (n> 2) are smaller than v 2 . They also carry important information on the dynamics of the medium created, and complement v 2 in providing a more complete picture of its bulk properties. The two coefficients that have been closely studied are the quadrangular (or hexadecapole) flow v 4 [13,14] and triangular flow v 3 [15]. Although the pentagonal and hexagonal flows v 5 and v 6 are studied to a lesser extent, there exist some predictions from hydrodynamics on them also [16]. At relatively low transverse momenta, p T < 3 ÷ 4 GeV/c, the azimuthal anisotropy results from a pressuredriven anisotropic expansion of the created matter, with more particles emitted in the direction of the largest pressure gradients [17]. At higher p T , this anisotropy is understood to result from the path-length dependent energy loss of partonic jets as they traverse the matter, with more jet particles emitted in the direction of shortest pathlength [18].
In Ref. [19] the LHC data on multiplicity, charged hadron spectra, elliptic flow and femtoscopic correlations from PbPb collisions were analyzed in the frameworks of the HYDJET++ model [20]. Taking into account both hard and soft components and tuning input parameters allow HYDJET++ to reproduce these data. Another study [21] with HYDJET++ was dedicated to the influence of jet production mechanism on the ratio v 4 /v 2 2 and its role in violation of the number-of-constituent-quark (NCQ) scaling [22], predicted within the HYDJET++ in [23]. In the current paper, tuned HYDJET++ is applied to analyze the LHC data on momentum and centrality dependences of azimuthal anisotropy harmonics in PbPb collisions, and then to illuminate the mechanisms of the generation of Fourier coefficients v 2 ÷ v 6 . The detailed study of hexagonal flow v 6 is also the subject of our recent paper [24].
Note that the LHC data on higher-order azimuthal anisotropy harmonics (v 2 ÷v 4 ) were analyzed with a multiphase transport model (AMPT) in [25]. It was shown that AMPT describes LHC data on the anisotropic flow coefficients v n (n=2÷4) for semi-central PbPb collisions at p T < 3 GeV/c. It also reproduces reasonably well the centrality dependence of integral v n for all but most central collisions. Another approach [26] reproducing v n data in ultrarelativistic heavy ion collisions is the glasma flow with the subsequent relativistic viscous hydrodynamic evolution of matter through the quark-gluon plasma and hadron gas phases (IP-Glasma+MUSIC model). This model gives good agreement to p T -dependence of v n (n=2÷5) and event-byevent distributions of v 2 ÷ v 4 at RHIC and LHC.
The study of generation of higher flow harmonics within the HYDJET++ has several attractive features. Firstly, the presence of elliptic and triangular flow permits us to examine the interference of these harmonics and its contribution to all higher even and odd components of the anisotropic flow. If necessary, the original eccentricities of higher order can be easily incorporated in the model for the fine tuning of the distributions. Secondly, very rich table of resonances, which includes about 360 meson and baryon species, helps one to analyze all possible final state interactions. Thirdly, the interplay of ideal hydrodynamics with jets can unveil the role of hard processes in the formation of anisotropic flow of secondary hadrons. The basic features of the model are described in Sect. 2.

HYDJET++ model
HYDJET++ (the successor of HYDJET [27]) is the model of relativistic heavy ion collisions, which includes two independent components: the soft state (hydro-type) and the hard state resulting from the in-medium multi-parton fragmentation. The details of the used physics model and simulation procedure can be found in the HYDJET++ manual [20]. Main features of the model are sketched below as follows.
The soft component of an event in HYDJET++ is the "thermal" hadronic state generated on the chemical and thermal freeze-out hypersurfaces obtained from the parametrization of relativistic hydrodynamics with preset freeze-out conditions (the adapted event generator FAST MC [28,29]). Hadron multiplicities are calculated using the effective thermal volume approximation and Poisson multiplicity distribution around its mean value, which is supposed to be proportional to a number of participating nucleons for a given impact parameter of a AA collision. To simulate the elliptic flow effect, the hydro-inspired parametrization is implemented for the momentum and spatial anisotropy of a soft hadron emission source [20,30].
The model used for the hard component in HYDJET++ is based on the PYQUEN partonic energy loss model [27].
The approach describing the multiple scattering of hard partons relies on accumulated energy loss via gluon radiation which is associated with each parton scattering in expanding quark-gluon fluid. It also includes the interference effect in gluon emission with a finite formation time using the modified radiation spectrum dE/dx as a function of the decreasing temperature T . The model takes into account radiative and collisional energy loss of hard partons in longitudinally expanding quark-gluon fluid, as well as the realistic nuclear geometry. The simulation of single hard nucleon-nucleon sub-collisions by PYQUEN is constructed as a modification of the jet event obtained with the generator of hadron-hadron interactions PYTHIA 6.4 [31]. Note, that Pro-Q20 tune was used for the present simulation. The number of PYQUEN jets is generated according to the binomial distribution. The mean number of jets produced in an AA event is calculated as a product of the number of binary NN sub-collisions at a given impact parameter per the integral cross section of the hard process in NN collisions with the minimum transverse momentum transfer p min T (the latter is an input parameter of the model). In HYDJET++, partons produced in (semi)hard processes with the momentum transfer lower than p min T , are considered as being "thermalized". So, their hadronization products are included "automatically" in the soft component of the event. In order to take into account the effect of nuclear shadowing on parton distribution functions, we use the impact parameter dependent parametrization [32] obtained in the framework of Glauber-Gribov theory.
The model has a number of input parameters for the soft and hard components. They are tuned from fitting to experimental data values for various physical observables, see [19] for details.
In order to simulate higher azimuthal anisotropy harmonics, the following modification has been implemented in the model. HYDJET++ does not contain the fireball evolution from the initial state to the freeze-out stage. Instead of application of computational relativistic hydrodynamics, which is extremely time consuming, HYDJET++ employs the simple and frequently used parametrizations of the freeze-out hypersurface [20]. Then, anisotropic elliptic shape of the initial overlap of the colliding nuclei results in a corresponding anisotropy of the outgoing momentum distribution. To describe the second harmonic v 2 the model utilizes coefficients δ(b) and ǫ(b) representing, respectively, the flow and the coordinate anisotropy of the fireball at the freeze-out stage as functions of the impact parameter b. These momentum and spatial anisotropy parameters δ(b) and ǫ(b) can either be treated independently for each centrality, or can be related to each other through the dependence on the initial ellipticity ǫ 0 where R A is the nucleus radius. The last option allows us to describe the elliptic flow coefficient v 2 for most centralities at the RHIC [20] and LHC [19] energies using only two independent on centrality parameters.
Non-elliptic shape of the initial overlap of the colliding nuclei, which can be characterized by the initial triangular coefficient ǫ 03 (b), results in an appearance of higher Fourier harmonics in the outgoing momentum distribution. Our Monte-Carlo (MC) procedure allows us to parametrize easily this anisotropy via the natural modulation of final freeze-out hypersurface, namely where φ is the spatial azimuthal angle of the fluid element relatively to the direction of the impact parameter. R(b, φ) is the fireball transverse radius in the given azimuthal direction φ with the scale R f (b), which is a model parameter. The phase Ψ RP 3 allows us to introduce the third harmonics possessing its own reaction plane, randomly distributed with respect to the direction of the impact parameter (Ψ RP 2 = 0). This new anisotropy parameter ǫ 3 (b) can again be treated independently for each centrality, or can be expressed through the initial ellipticity ǫ 0 (b) = b/2R A . Note, that such modulation does not affect the elliptic flow coefficient v 2 , which was fitted earlier with two parameters δ(b) and ǫ(b) [19,20]. Figure 1 illustrates second and third harmonics generation in HYDJET++ by representing particle densities in the transverse plane. One should be aware that the triangular deformation shown here is very strong. The actual deformations needed to describe triangular flow at LHC energies are typically order of magnitude weaker.
The modulation of the maximal transverse flow rapidity, first considered in Eq. (28) of Ref. [20] at the parametrization of 4-velocity u, also permits the introduction of higher azimuthal harmonics related, however, to the direction of the impact parameter (Ψ RP 2 = 0) only. In this case we get the modulation of the velocity profile in all freeze-out hypersurface, and can not "rotate" this modulation with independent phase. The new anisotropy parameters, ρ 3u (b) and ρ 4u (b), can again be treated independently for each centrality, or can be expressed through the initial ellipticity ǫ 0 For current simulations we have introduced the minimal modulation in HYDJET++ using just simple parame- being taken equal to zero. The corresponding proportionality factors were selected from the best fit of the data to v 3 (p T ) and v 4 (p T ).
Let us mark that the azimuthal anisotropy parameters ǫ(b), δ(b) and ǫ 3 (b) are fixed at given impact parameter b. Therefore they do not provide dynamical eventby-event flow fluctuations, and specify v n (b) accumulated over many events. The main source of flow fluctuations in HYDJET++ is fluctuations of particle momenta and multiplicity. Recall, that the momentum-coordinate correlations in HYDJET++ for soft component is governed by collective velocities of 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 collision impact parameter for each centrality class also increases such fluctuations. In the current paper we restrict ourselves to analysis of the eventaveraged v n (p T ). The detailed study of event-by-event flow fluctuations is the subject of our future investigation. The possible further modification of HYDJET++ to match experimental data on flow fluctuations would be smearing of parameters ǫ, δ and ǫ 3 at a given b.

Results
It was demonstrated in [19] that tuned HYDJET++ model can reproduce the LHC data on centrality and pseudorapidity dependence of inclusive charged particle multiplicity, p T -spectra and π ± π ± correlation radii in central PbPb collisions, and p T -and η-dependencies of the elliptic flow coefficient v 2 (up to p T ∼ 5 GeV/c and 40% centrality). However the reasonable treatment of higher and odd Fourier harmonics of particle azimuthal distribution v n (n > 2) needs the additional modifications of the model, which does not effect azimuthally-integrated physical observables (see previous section). We have compared the results of HYDJET++ simulations with the LHC data on v n for inclusive as well as for identified charged hadrons.

Anisotropy harmonics for inclusive charge hadrons
The standard way of measuring v n corresponds to the inclusive particle harmonics on the base of Eq. (1). Then v n is extracted using the special methods, such as the event plane v n {EP} [33], or m-particle cumulant v n {m} [34,35], or Lee-Yang zero methods v n {LYZ} [36,37]. In order to estimate the uncertainties related to the experimental definitions of flow harmonics, HYDJET++ results for different methods of v n extraction were compared with its "true" values, known from the event generator and determined relatively to Ψ RP 2 for even and Ψ RP 3 for odd harmonics, respectively. Figures 2-11 show anisotropic flow coefficients v n as a function of the hadron transverse momentum p T . Let us discuss first the results of HYDJET++ simulations. It can be separated in two groups: (i) results obtained with respect to the true reaction plane straight from the generator, i.e., v 2,4,6 (Ψ RP 2 ) and v 3,5 (Ψ RP 3 ), and (ii) those obtained by using the (sub)event plane method with rapidity gap |∆η| > 3. The last method provides us with v n {EP}. The main systematic uncertainties for the methods come from non-flow correlations and flow fluctuations. The last one (as it is kept in the model currently) almost does not affect mean v n values restored by the EP method, while the non-flow correlations can be effectively suppressed by applying η-gap in v n reconstruction. This gives us a good reconstruction precision for elliptic v 2 , triangular v 3 , and quadrangular v 4 flows up to p T ∼ 5 GeV/c. At higher transverse momenta some differences appear due to nonflow effects from jets. However, Figs. 8 and 9 show that pentagonal flow v 5 determined from the model w.r.t. Ψ RP 3 and v 5 restored w.r.t. the event plane of 5-th order Ψ EP 5 differ a lot. The reason is that although no intrinsic Ψ RP 5 is generated in HYDJET++, pentagonal flow v 5 emerges here as a result of the "interference" between v 2 and v 3 , each is determined with respect to its own reaction plane, , in line with the conclusions of Ref. [38]. Hexagonal flow v 6 is also very sensitive to the methods used due to nonlinear interplay of elliptic and triangular flows generating v 6 , see [24] for details. The results of HYDJET++ for v 6 {EP} are not shown on the plots because of too large statisitcal errors.
Note, that the experimental situation is even more complicated, and the dependence of measured v n on methods applied may be more crucial for all n due to apparently larger fluctuations in the data than in the model. For instance, it was shown in [39] that event-by-event fluctuations in the initial state may lead to characteristically different p T -dependencies for the anisotropic flow coefficients extracted by different experimental methods.
It is also worth mentioning here that the hump-like structure of the simulated v 2 (p T ) and v 3 (p T ) signals appears due to interplay of hydrodynamics and jets. At transverse momenta p T ≥ 3 GeV/c the spectrum of hadrons is dominated by jet particles which carry very weak flow. Thus, the elliptic and triangular flows in the model also drop at certain p T . Higher flow harmonics arise in the model solely due to the presence of the v 2 , v 3 and its interference. Therefore, transverse momentum distributions of these harmonics inherit the characteristic hump-like shapes. Now let us consider the ATLAS [10] and CMS [8,11] data plotted onto the model results in Figs. 2-11 for different centrality classes. The event plane for v n {EP} was defined experimentally with respect to n-th harmonics in all cases with the exeption of CMS data for v 6 {EP/Ψ 2 }, which was measured using second harmonics. One can see that HYDJET++ reproduces experimentally measured p Tdependences of v 2 , v 3 and v 4 {LYZ} up to p T ∼ 5 GeV/c. The centrality dependence of v 4 measured by event plane and two-particle cumulant methods is significantly weaker than that of v 4 measured by Lee-Yang zero method, presumably due to large non-flow contribution and increase of the flow fluctuations in more central events. Since the model is tuned to fit the p T −dependencies of v 4 {LYZ}, it underestimates the quadrangular flow, restored by the EP or two-particle cumulant methods, in (semi-)central collisions. Recall, that in ideal hydrodynamics (at the limit of small temperatures, large transverse momenta and absence of the flow fluctuations) v 4 {Ψ 2 }/v 2 2 = 0.5 [40]. The same trend is seen for p T -dependencies of the pentagonal flow. For central and semi-central topologies up to σ/σ geo ≈ 20% the v 5 {EP} in the model underestimates the experimentally measured v 5 {EP}, whereas for more peripheral collisions the agreement between the model and the data is good. Unfortunately, there are no data on pentagonal flow extracted by the LYZ method. As we have seen, for v 2 , v 3 and v 4 in central and semi-central collisions the LYZ method provides noticeably weaker flow compared to that obtained by the EP method. One may expect, therefore, that the pentagonal flow, v 5 {LYZ}, almost free from non-flow contributions, should be closer to the v 5 generated by the HYDJET++. If the future experimental data on v 5 will persist on stronger flow, this fact can be taken as indication of the possible presence of additional pentagonal eccentricity ǫ 5 (b) with the new phase Ψ RP 5 responsible for genuine v 5 . Both parameters can be easily inserted in Eq. (2) for the modulation of the final freeze-out hypersurface.
At last, p T -dependencies of the hexagonal flow in HY-DJET++ are similar to that seen in CMS data within the uncertainties related to methods used. However v 6 (Ψ RP 2 ) in the model visibly underestimates ATLAS data on v 6 {EP} for most central events. The latter fact may be explained by a siginificant v 3 contribution to v 6 {EP} in central collisions, which is not presented in v 6 (Ψ RP 2 ) component: On the other hand, the relative contribution to v 6 {EP} coming from v 2 is instantly increasing as the reaction becomes more peripheral [24], and starting from 20−30% centralities we already get v 6 {EP} ∼ v 6 (Ψ RP 2 ) ≫ v 6 (Ψ RP 3 ) with the approximate agreement between the model and the data.

Anisotropy harmonics for identified charge hadrons
Finally, let us consider distributions for some hadronic species measured in PbPb collisions at the LHC. Before addressing to azimuthal anisotropy harmonics of identified hadrons, the comparision of HYDJET++ results with AL-ICE data [41] on p T -spectra of negatively charged pions, kaons and anti-protons in PbPb collisions is displayed in Fig. 12. One can see that HYDJET++ reproduces well the measured transverse momentum spectra of identified hadrons within the whole range of accessible p T . Figure 13 presents the comparision of HYDJET++ results and the ALICE data [42] for the elliptic and triangular flow of pions, kaons and anti-protons at 10-20% and 40-50% centrality of PbPb collisions. The agreement between the model and the data for kaons and anti-protons looks fair. For pions the model underestimates the data a bit. The discrepancy is more pronounced for more central collisions indicating, perhaps, presence of strong non-flow correlations in the data.

Conclusion
Azimuthal anisotropy harmonics of inclusive and identified charged hadrons in PbPb collisions at √ s NN = 2.76 TeV have been analyzed in the framework of HYDJET++ model. The effects of possible non-elliptic shape of the initial overlap of the colliding nuclei are implemented in HY-DJET++ by the modulation of the final freeze-out hypersurface with the appropriate fitting triangular coefficient.
This modulation is not correlated with the direction of the impact parameter, and two independent "strong" lower azimuthal harmonics, v 2 and v 3 , being obtained as a result. They are of different physical origin, coded partly in the different centrality dependence. Interference between v 2 and v 3 generates as "overtones" both even and odd higher azimuthal harmonics, v 4 , v 5 , v 6 , etc. This mechanism allows HYDJET++ to reproduce the LHC data on p T -and centrality dependencies of the anisotropic flow coefficients v n (n=2÷4) up to p T ∼ 5 GeV/c and 40% centrality, and also the basic trends for pentagonal v 5 and hexagonal v 6 flows. Some discrepancy between the model results and the data on the pentagonal flow in central events requires further study of additional sources of the non-flow correlations and flow fluctuations, which may be absent in the model. Although the introduction of internal higher harmonics is also possible in the HYD-JET++, there is no clear evidence in the data to do so at present. Obtained results show that higher harmonics of the azimuthal flow get very significant contributions from the lower harmonics, v 2 and v 3 . This circumstance makes it difficult to consider the higher harmonics as independent characteristics of the early phase of ultrarelativistic heavy ion collisions.   (d) 40-50% Fig. 12. Negatively charged pion, kaon and anti-proton transverse momentum spectra at pseudo-rapidity |η| < 0.5 for different centralities of PbPb collisions at √ s NN = 2.76 TeV. The points are ALICE data [41], histograms are the simulated HYDJET++