Event shape Engineering analysis of D meson in ultrarelativistic heavy ion collisions

We describe the propagation of charm quarks in the quark-gluon plasma (QGP) by means of an event-by-event transport approach. In our calculations the non-perturbative interaction between heavy quarks and light quarks has been taken into account through a quasi-particle approach with thermal light quark masses tuned to reproduce lQCD thermodynamics. We found that the flow observables $v_2$ and $v_3$ of D mesons are comparable with the experimental measurements for Pb+Pb collisions at 5.02 ATeV in different ranges of centrality selections. The results are analyzed with Event-Shape Engineering technique. The comparison of the anisotropic flow coefficients $v_n$ with experimental data show a quite well agreement with experimental data for different flow vector $q_2$ selections, which confirms the strong coupling between charm quarks and light quarks in the QCD matter. Furthermore, we present here a novel study of the event-by-event correlations between flow harmonics of $D$ mesons and soft hadrons at LHC energy with the Event-Shape Engineering technique that can put further constraints on heavy quark transport coefficients toward a solid comparison between the phenomenological determination and the lattice QCD calculations.


I. INTRODUCTION
Heavy quarks (HQs), mainly charm and bottom quarks, are produced in the early stages of a ultrarelativistic Heavy Ion collision (uRHIC) by hard processes described by perturbative QCD calculations. They have been identified as one among the few probes which may allow for a direct study of the QGP properties. This is manly due to the following reasons related to their large masses: i) they have a very short formation time (τ ∼ 0.1 fm/c), ii) HQs are not expected to reach a full thermalization. Therefore, they can probe the entire space-time evolution of the matter created in these collisions. Furthermore, due to their interaction with the medium produced (quarks and gluons) they can keep this information in their final state as constituent of heavy hadrons mainly D, B mesons and Λ c , Λ b baryons. Among the several observables studied in uRHICs two of them have been long investigated for HF hadrons: the heavy mesons nuclear modification factor R AA (p T ) [1][2][3], and the so called elliptic flow, v 2 (p T ). The first gives an information about the change of the spectrum in AA collision with respect to a simple pp superposition. The second gives a measure of the anisotropy in the particle angular distribution and allows to investigate the coupling of the HQs with the medium and their participation in the collective motion. Both these quantity have been extensively used as probe of the Quark Gluon Plasma [4,5]. From a theoretical point of view, in recent years, several studies have done efforts in the direction of studying both of these observables to understand the heavy quark dynamics in QGP * Electronic address: salvatore.plumari@dfa.unict.it [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. More recent measurements on the D mesons and Λ c production in nuclear collisions unfolded the possibility to get access to richer information. In particular, there have been studies that have opened the way to start to address the issue of medium-modification of heavy-flavour hadronization making possible to get information about heavy-quark hadronization and possibly to validate models based on the recombination [23], and recently these studies have been extended to the case of small collision systems like proton-nucleus and proton-proton collision providing insight on the heavy flavor hadronization mechanism and the possibility that also in such small systems the QGP can be formed [24]. More information about the interactions of heavy quarks with the medium can be obtained through measurements of the azimuthal distributions of heavy-flavour hadrons produced in a HIC. Because of the interaction of heavy quarks with the medium, the spatial anisotropy in the early stages of nucleus-nucleus collisions is converted into an azimuthally anisotropic distribution in momentum space for the final particles produced from the QGP. This anisotropy is characterised in terms of the Fourier coefficients v n of the final distribution respect to the symmetry-plane angles Ψ n (for the harmonic of order n). The second coefficient of this expansion is the elliptic flow v 2 which is the dominant contribution in non-central collisions. The higher order harmonics usually studied are the triangular flow v 3 , and the v 4 ; the measurement of odd harmonics of heavy-flavour hadrons or their eventby-event fluctuations coming from the random nucleon positions can provide a richer information on the initial conditions. Recently, the triangular flow v 3 has been investigated in studies based on Langevin or Boltzmann approaches on an event-by-event analysis [25][26][27][28][29]. One of the new tools that can be used to further investigate the dynamics of heavy quarks in the medium is the Event Shape Engineering (ESE) technique [30]. With this procedure it is possible to select events in heavy-ion collisions based on the geometry of the fireball, enabling to select sub-samples of events characterized by an high or low eccentricity, furnishing the chance to study the system evolution dynamics and the role of the initial conditions. This technique was first used by the ALICE collaboration in the light-flavour sector to study the interplay between the initial geometry of the nucleus-nucleus collisions and the subsequent evolution of the system [31] and recently it was extended to investigate the sensitivity of the D meson v 2 to the event-shape selection [32]. These measurements have shown a sizeable sensitiveness to the event-shape selection, confirming a correlation between the D-meson azimuthal anisotropy and the underlying collective expansion of the bulk matter. In this paper we use an event-by-event transport approach for the evolution of heavy quarks and we perform an analysis of the heavy-quark production using the ESE technique, in order to study the correlation of heavy-quark anisotropic flows v n − v m correlations.
The paper is organized as follows. In section II we discuss briefly the Boltzmann transport approach used to describe the HQs evolution and the hadronization model used, furthermore we will describe the initial conditions employed in the simulation. In sections III we discuss the comparison between the simulation results and the experimental data for the anisotropic flows coefficients v n at different centralities. In section IV we introduce the ESE technique and we show and discuss the outcome for the anisotropic flows coefficients and their correlations. Finally, section V contains a summary and some concluding remarks.

II. HEAVY QUARKS TRANSPORT APPROACH
We use a transport code developed to perform studies of the dynamics of relativistic heavy-ion collisions at both RHIC and LHC energies [33][34][35][36][37][38][39][40][41][42][43]. The space-time evolution of gluons and light quarks as well as of charm quarks distribution functions is described by mean of the Relativistic Boltzmann Transport (RBT) equations given by where f i (x, p) is the on-shell phase space one-body distribution function for the i−th parton and in the right-hand side C[f q , f g , f Q ](x, p) is the relativistic Boltzmann-like collision integral accounting for 2 → 2 scattering processes. For the quarks and gluons distribution functions that compound the QGP medium, the evolution is described by Eqs. (1) and (2) in this paper. We assume that these Eqs. are independent of f Q , which correspond to discard collisions between charm quarks and the bulk in the determination of light quarks and gluons distribution functions, which is by far a solid approximation. In both collision integrals the total cross section is determined in order to keep the ratio η/s = 1/(4π) fixed during the evolution of the QGP, see Refs [36,40,42,44] for more details. Furthermore, in the present paper we have employed a bulk with thermal massive quarks and gluons according to a Quasi-Particle Model (QPM) such that the system has approximately the lattice QCD equation of state, with a softening of the equation of state consistent with a decreasing speed of sound approaching the cross-over region [45]. For the heavy quark interacting with the bulk medium we consider a QPM whose main feature is that the resulting coupling is significantly stronger than the one coming from pQCD running coupling, particularly as T → T c (see details in Refs [17,41]). Numerically we solve the RBT equation using the test particle method, where the collision integral is solved by Monte Carlo method based on stochastic interpretation of transition amplitude [33,46,47]. In our calculations the viscosity is fixed by determining the total isotropic cross section according to the Chapman-Enskog approximation: where z = m/T while the function f (z) is defined by the following expression f (z) = 15 16 where K n (z) are the modified Bessel functions. As shown in Ref [33] the above expression for the shear viscosity η is in quite good agreement with the Green-Kubo formula. At the end of the evolution at freeze-out temperature, we use the hybrid hadronization approach via coalescence plus fragmentation as shown in Refs [23].
In our simulation for parton initial conditions, we have employed a modified Monte Carlo Glauber model, used in [43] and inspired by wounded quark model. In this model each nucleon is composed by three constituent quarks that are randomly distributed within each nucleon according to the distribution dN/dr = r 2 r 3 0 e −r/r0 with r 0 = 0.3 f m and the final center of mass of the three quarks is translated to the position of the nucleon. The positions of nucleons in P b nuclei are instead distributed according to the standard Woods-Saxon distribution. In this way we generate the wounded quark profile using Monte-Carlo Glauber model and we decide whether each quark pair from target and projectile can collide or not with a probability p = e −πr 2 /σ qq with σ qq = 13.6 mb in 5.02 ATeV P b+P b collision [48]. Finally, the total initial parton distribution is given by: where N part is the number of participant quarks while x i and n i are the transverse position of each participant quark and the number of partons generated by each participant respectively. The profile in the transverse position ρ ⊥ (x ⊥ ) is expressed as follows: While the profile for the distribution of the spatial rapidity ρ(η) is given by longitudinal boost invariance in pseudo-rapidity. In Eq.6 the numbers n i fix the number of partons produced in each event. This number is given by n 0 N where N is sampled according to a negative binomial distribution given by where n 0 is fixed in order to have that the final charged particle multiplicity is same as that measured in experiments. This model is inspired by the distribution of the number of particles produced in pp collisions that fluctuates according to a negative binomial distribution. In order to reproduce the distribution of charged particles measured by ALICE Collaboration the parameters have been fixed to k = 0.224,n = 1.621 and n 0 ≈ 1.88.
In momentum space we have considered a mixture of a soft and an hard component, with the former consisting of quarks and gluons with an initial thermal equilibrium distribution, and the latter being equivalent to the products of initial binary pQCD collision consisting of 'minijets'. The soft partons are distributed according to the following (10) where g = 2×8+3×2×6 = 52 is the degree of freedom of partons (three flavor quarks and gluons), τ 0 = 0.4f m/c is the termalization time and m T = m 2 + p 2 T is the transverse mass of light partons. The transverse momentum distributions of minijets at mid-rapidity are taken from pQCD calculations and we have used the spectra of gluons and quarks from CUJET Collaboration [49] for pp collisions at 5.02 TeV.
In coordinate space we initialize the charm quark distribution in accordance with the number of binary nucleon-nucleon collisions, N coll , from the Monte Carlo Glauber model. Finally, for the charm quark initial distribution in momentum space we use the production calculated at Fixed Order + Next-to-Leading Log (FONLL) [50] which describes the D-meson spectra in proton-proton collisions after fragmentation. For detail we refer to our earlier work in Ref. [41].
In our simulation when the temperature of the QGP phase goes below the quark-hadron transition temperature, T c = 155 MeV, the charm quark can hadronize into D-meson. For the heavy quark hadronization process, we consider a hybrid approach of hadronization by coalescence plus fragmentation. The fragmentation is implemented using the Peterson fragmentation function [51]: is the momentum fraction of the D-meson fragmented from the charm quark and c is a free parameter to fix so that one can describe D-meson production in proton-proton collisions [52]. For the coalescence model, it is based on a Wigner formalism. In our calculation the Wigner function is assumed to be Gaussian where the widths are fixed by the mean square charged radius of the D meson, for detail we refer to our earlier work [23].

III. HEAVY FLAVOUR ANISOTROPIC FLOWS vn(pT )
In a non-central nucleus-nucleus collision, the interaction of the QGP constituents, convert the initial eccentricity of the overlap region into a final anisotropy in momentum space, characterised in terms of the Fourier coefficients v n (p T ). In our model, event-by-event fluctuations generate a profile in the transverse plane ρ ⊥ (x ⊥ ) that change with the event and it determines in each event the initial anisotropy in coordinate space. This is quantified in terms of different order spatial asymmetries named n and expressed by: Ψ n = 1 n arctan r n ⊥ sin(nφ) r n ⊥ cos(nφ) (12) where r ⊥ = x 2 + y 2 and φ = arctan(y/x). The results shown in this paper have been obtained using the two particle correlation method to calculate elliptic (v 2 ) and triangular flow (v 3 ) and v 4 [53,54]. We have found that the number of test particles N test ∼ 350 and with the lattice discretization ∆x = ∆y = 0.3 fm in the simulations guarantee a convergence for v 2,3 (p T ) up to p T ∼8 GeV/c.
In Fig. 1, we show the D mesons anisotropic flow coefficients v 2 (p T ) and v 3 (p T ) at mid-rapidity for P b + P b collisions at 5.02 AT eV for central (0 − 10%) in the upper panel and semi-peripheral(30 − 50%) collisions in the bottom panel. The red line corresponds to v 2 (p T ) and blue line to v 3 (p T ) both for charm quarks (dashed line) and D mesons(solid line). The charm quarks v 2 and v 3 are different from zero in both centrality classes suggesting that charm quarks take part in the collective motion. We observe that the v 2 is sensitive to the centrality of the collisions while the v 3 has a weak dependence with the centrality and it is comparable for both centralities. This suggests that the elliptic flow is mainly generated by the geometry of overlapping region and it is larger for larger centrality collision while triangular flow is mainly driven by the fluctuation of the triangularity 3 of overlapping region , and similar results observed also for light quarks [40,56]. The hadronization mechanism plays a role for the v 2 to describe D meson anisotropic flow at intermediate p T range. The effect of coalescence corresponds to an increase of the collective flows at p T > 2 GeV due to the fact that D mesons are formed by the recombination of a light and a heavy quarks. This effect is mainly due to the elliptic flow of light quarks which is larger than the one of charm quarks. The process of the couple recombination gives as a result a final hadron elliptic flow larger respect to the simple anisotropy of the initial charm, because it carries information from both quarks anisotropies. On the other hand, D mesons coming from fragmentation have always a smaller momentum with respect to the original charm quark, as a result the fragmentation contribution to the elliptic flow is similar to the one at charm level but shifted to lower p T . As shown in Fig. 1, our results are in good agreement with the ALICE experimental data [55] and capture the centrality dependence observed experimentally.

IV. EVENT-SHAPE ENGINEERING TECHNIQUE
Further investigation into the dynamics of heavy quarks in the medium can be performed with the eventshape engineering (ESE) technique [30]. This technique is based on the observation that the event-by-event variation of the anisotropic flow coefficient v n at fixed centrality is very large [57]. This technique allows for selection of events with the same centrality but different initial geometry on the basis of the magnitude of the average bulk flow q 2 . In this section, we present the main features of the ESE selection and our results with this approach focusing on the prediction for anisotropic flows correlations v n −v m in the heavy flavor sector. In order to apply the ESE technique, in our calculations the events in each centrality class are divided according to the magnitude of the second-order harmonic reduced flow vector q 2 which is defined as: where Q 2 is a vector built from the azimuthal angles φ k of the considered particles and given by while M is the number of charged particles in the specific range of transverse momentum p T and pseudorapidity η. In Fig. 2, it is shown q 2 distribution of charged particles with |η| < 0.8 and in the transverse momentum range 0.2 < p T < 5 GeV /c for 30 − 50% centrality class in P bP b collisions at 5.02 AT eV . The selection of the events according to the average elliptic flow was performed by defining q 2 percentiles in 1%-wide centrality intervals to make our results comparable with the experimental data [55]. Our results shown by the solid line are in in good agreement with the ALICE measurements [58] (dashed line). This kind of study can be performed with the implementation of event-by-event initial state fluctuations, and we are performing one of the first calculation of this type concerning the dynamics of heavy quark. A first study on the heavy-light flavor correlations of anisotropic flows at LHC energies with an event-by-event transport approach has been performed by means of a full Boltzmann transport approach in Ref. [29]. In this paper we extend that work and evaluate the correlations between different order anisotropic flows v n − v m both in the light and heavy flavor sectors with the ESE technique. In Fig. 3, we have shown different spatial bulk eccentricities n as a function of reduce flow vector q 2 for P b+P b collisions at √ s N N = 2.76AT eV in 30−50% centrality class. The results show that large q 2 corresponds to events where the fireball has a large initial eccentricity 2 . The triangularity 3 show a weak anti-correlation while 4 receives a non linear contribution with respect to q 2 . The correlation between the spatial eccentricities nm that is present in the initial geometry gives rise to correlation between different flow harmonics v n -v m . Hydrodynamic and transport calculations have shown that the response of the system to the initial spatial anisotropy is essentially linear for the second and third harmonics, meaning that the final state v 2 and v 3 are very well correlated with the second and third order eccentricities in the initial state, for small values of η/s [40,[59][60][61]. In Fig. 4 the correlations between v 2 and v 3 (solid black line), v 2 and v 4 (red solid line) for charged particles in P bP b collisions at √ s N N = 2.76AT eV in 15 − 20% centrality class are shown in comparison with the experimental data from ALICE collaboration [62]. In these calculations, for an established centrality class, every point is determined evaluating the v n varying the q 2 . The trend of the v n − v 2 correlations is related to the initial correlation (anti-correlation) between n and that having, for a determined centrality, a specific number of participating nucleons can be ideally arranged to have a large 2 , the geometric consequence of this config- uration leads to a smaller 3 . These results are in agreement with other measurement performed by the ATLAS collaboration [63]. As shown in Fig. 4 the v 4 has a non linear correlation with the v 2 . This behaviour can be attributed to the fact that v 4 receives a non-linear contribution from the initial 2 that goes like ( 2 ) 2 [56]. On the other hand, the v 4 can be parametrized as the quadrature sum of two components: one proportional to v 2 2 (as v 2 ∝ 2 ), and represents the non-linear component of v 4 driven by ( 2 ) 2 , and the other contribution that in principle could be a linear component that should be independent of v 2 [64]. Our results for charged particles are in good agreement with the experimental data suggesting that our approach may capture the initial spatial fluctuations of the bulk matter. In the same scheme in the ones of the bulk.
In order to perform a more systematic study we can now introduce a new type of observable for anisotropic flow analyses, the so-called Symmetric Cumulant (SC(m, n)), proposed in [65]. This observable is particularly useful for systems in which flow harmonics fluctuate in magnitude event-by-event, and it is defined as follows, starting from the 4-particle correlation where the φ i are the particle azimuthal angles, while the double averaging is a mean over all particles and over all events for a given centrality class [65]. In the left panel the experimental data taken from ALICE collaboration [66] we observe that SC(4, 2) and SC(3, 2) decrease considerably in magnitude from peripheral to central events. This large decrease does not imply that the correlations between these harmonics becomes negligible for centrality class smaller then (0 − 10)%. This decreasing in the magnitude of the correlation occurs since the SC ∼ v 4 n and for very central events the v n (especially v 2 ) become quite small. The experimental data are taken from ALICE measurements [55].
In Fig.7 we show the results for the D mesons v 2 at mid-rapidity for two centralities and for different q 2 selections: 20% large-q 2 (red solid and dashed lines) and 20% small-q 2 (blue solid and dashed lines) events. As shown in Fig. 7, the v 2 of 20% large-q 2 (red solid line) events is quite larger than the v 2 of 20% small-q 2 events (blue solid line) for most central 0 − 10% and semi-peripheral 30 − 50% centrality classes. This result suggests that the selection in q 2 distinguish between the initial eccentricities, large (small) q 2 means large (small) eccentricity, that consequently develop in large or small transverse flow translated into the final azimuthal anisotropy in momentum space for both bulk and D-meson. This difference increases with the collision centrality. Our results are compared to the ALICE measurements [55], where the good agreement shows that our model can essentially capture the strong coupling between the bulk medium and heavy quark. Finally, in Fig. 8 we show the ratio as a function of the transverse momentum between the D meson v 2 in events characterized by large or small q 2 divided by the v 2 in unbiased events. We find, for these ratios, a difference of about 50% between the q 2 -selected and the unbiased events in both centrality class, with no dependence on the transverse momentum, this behaviour might suggest that the ESE selection is related to a global property of the events.

V. CONCLUSION
In this paper we have studied the propagation of charm quarks in the quark-gluon plasma (QGP) by means of an event-by-event transport approach. In particular we have solved the relativistic Boltzmann transport equation with an event-by-event fluctuating initial condition which allow us to access the study of high order anisotropic flow coefficients v n for charmed mesons. We have studied the D meson v n at LHC energies in different centrality class selection. In particular, we have evaluated for the first time the heavy flavor v n −v m correlations with the Event-Shape Engineering (ESE) technique consisting in selecting events in the same centrality class but characterized by different average elliptic anisotropy. With our approach we have a good agreement with the available experimental data from ALICE collaborations for the bulk. Furthermore, we have provided predictions for charmed hadrons v n − v m correlations and we predict a weaker but similar correlations behaviour respect to the light flavours reflected in a smaller Symmetric Cumulant Correlator SC(m, n) for D mesons. Like for the light flavor we observe that SC(4, 2) is positive, indicating that for a given centrality interval, events with larger v 2 on average have larger v 4 while the SC(3, 2) is negative indicating an anti-correlation between v 2 and v 3 at fixed centrality. For both cases we observe that SC(4, 2) and SC(3, 2) decrease from peripheral to central events. Finally, we have found that the selection based on the event-shape leads to a visible effect on the final particle elliptic flow. We have shown that the v 2 of 20% large-q 2 events is larger than the v 2 of 20% small-q 2 events for both most central 0 − 10% and semi-peripheral 30 − 50% centrality class, suggesting a correlation between the D-meson azimuthal anisotropy and the collective expansion of the bulk matter. We notice that the sensitivity with the event-shape selection increase with the collision centrality. Our results are in good agreement with the experimental data from the ALICE coll. [55].