Development of heavy-flavour flow-harmonics in high-energy nuclear collisions

We employ the POWLANG transport setup, developed over the last few years, to provide new predictions for several heavy-flavour observables in relativistic heavy-ion collisions from RHIC to LHC center-of-mass energies. In particular, we focus on the development of the flow-harmonics $v_2$ and $v_3$ arising from the initial geometric asymmetry in the initial conditions and its associated event-by-event fluctuations. Within the same transport framework, for the sake of consistency, we also compare the nuclear modification factor of the $p_T$ spectra of charm and beauty quarks, heavy hadrons and their decay electrons. We compare our findings to the most recent data from the experimental collaborations. We also study in detail the contribution to the flow harmonics from the quarks decoupling from the fireball during the various stages of its evolution: although not directly accessible to the experiments, this information can shed light on the major sources of the final measured effect.


Introduction
Experimental measurements of heavy-flavour production in relativistic heavy-ion collisions are a major tool to get information on the properties of the (deconfined, if a sufficiently high energy-density is achieved) medium formed in these events, in particular on its transport coefficients [1][2][3][4][5][6][7]. At high momentum the major effect of the interaction with the medium is a quenching of the heavy-quark momentum spectra due to parton energy-loss: this provides information on the medium opacity [2,4]. At low/intermediate momenta, on the other hand, if the transport coefficients were large enough, heavy quarks would even approach local thermal equilibrium with the rest of the medium, taking part in its collective expansion [8]. This would lead to clear signatures in the final observables: the radial and elliptic flow of the fireball arising from the heavy-ion collision would leave their fingerprints also in the heavy-flavour sector, boosting the heavy quarks from low to moderate momenta and giving rise to azimuthal anisotropies in their angular distributions [3,5]. Furthermore, since nowadays higher flow-harmonics (v 3 , v 4 , v 5 ...) of soft-hadron azimuthal distributions are measured (providing information on event-by-event fluctuations and granularity of the initial conditions), one would like to address this issue also in the heavy-flavour sector (for first experimental results, see Ref. [7]): this will be one of the major topics dealt with in this paper, based on the POWLANG transport setup developed by the authors over the last years [9][10][11][12]. A similar theoretical study was performed in [13] and, accounting only for the path-length dependence of parton energy-loss, in [14]. The long term goal would be to perform this kind of study on an event-by-event basis selecting, within the same centrality class, collisions characterized by different initial eccentricities or comparing events from different centrality classes but having a comparable initial eccentricity: we believe that this has the potential to further constrain the coupling of the heavy quarks with the background medium. However, the starting point is to check to be able to reproduce the trend of the experimental data in the case of the event-averaged results, which is the subject of the present paper.
Of course, non-trivial features in the heavy-flavour hadronic distributions experimentally measured can not be directly ascribed to the parent heavy quarks. As suggested by several studies, in the presence of a deconfined medium, rather then fragmenting like in the vacuum, heavy quarks may hadronize by recombining with the light thermal partons nearby to give rise to open charm/beauty hadrons. Belonging to the non-perturbative realm of QCD, there is no solid first-principle theory to describe hadronization, neither in the vacuum nor in the medium. The latter is modeled in several different ways in the literature -via coalescence [15][16][17], formation of color-singlet clusters/strings [11] or of resonances [18][19][20] -but the qualitative effect is always the same: the light thermal quark involved in the recombination process is part of a fluid cell sharing a common collective velocity and this provides an additional contribution to the (radial, elliptic and also triangular, as will be shown in the paper) flow of the final heavy-flavour hadron. Clearly, recombination with light partons from the medium, besides the kinematic distributions, can also affect the heavy-flavour hadrochemistry in nucleus-nucleus collisions, changing the relative yields of the various species with respect to the proton-proton case. This was modeled for instance in [20] and first experimental results are getting available [21][22][23], however we will not touch such an issue.
Concerning the flow acquired by charm and beauty quarks during the partonic phase it is of interest to disentangle the various sources of possible azimuthal anisotropies, in order to better understand how much heavy quarks really approach thermal equilibrium, tending to flow with the fireball, and how much of the final signal instead is simply due to trivial geometric effects. We will address this issue by studying the temporal development of the elliptic and triangular flow, disentangling the contribution to the final v 2 and v 3 from the heavy quarks decoupling at different times. A somehow similar analysis, referring to the bulk soft-particle production, was performed in [24] where the elliptic and triangular flow were studied within a transport model as a function of the number of collisions suffered by the partons; the authors found that the anisotropic escape probability of the partons, trivially arising from the initial geometry, provides a major contribution to the final signal, challenging the usual hydrodynamic interpretation of the data based on the formation of a strongly-interacting medium. Actually, our analysis deals only with the propagation of heavy quarks and is based on a more macroscopic approach, since the background medium is given a coarse-grained hydrodynamic description and the propagation of the heavy quarks throughout the fireball is not modeled through the individual collisions with the other partons, but just in terms -according to the Langevin equation -of an average squared-momentum exchange per unit time. Our findings will be presented and discussed in a devoted section; here we only anticipate that the final result come from a non-trivial interplay of contributions from the heavy quarks decoupling during all the stages of the fireball evolution.
Our paper is organized as follows. In Sec. 2 we present a detailed description of the transport equations implemented in the POWLANG setup. In Sec. 3 we describe how we model the initialization and the evolution of the background medium, in particular in the case of fluctuating initial conditions giving rise to a triangular flow. In Sec. 4 we study the temporal development of the heavy-quark v 2 and v 3 . In Sec. 5 POWLANG results for various heavy-flavour observables in nucleus-nucleus collisions are compared to recent experimental results obtained at RHIC and at the LHC. Finally, in Sec. 6 we summarize the main conclusions of our paper and outline possible future developments of our studies.

The transport setup
Different approaches are adopted in the literature to model the heavy-flavour transport throughout the plasma of light quarks and gluons expected to be produced in heavyion collisions. The POWLANG setup is based on the relativistic Langevin equation; for the latter different implementations can be found and this can be sometimes a source of confusion. Hence, here we briefly summarize the essential points of our transport scheme.
The starting point of any transport calculation is the Boltzmann equation for the evolution of the heavy-quark phase-space distribution The direct solution of the Boltzmann integro-differential equation is numerically demanding (for a detailed description of the approach see for instance [25]); however, as long as q ≪ p (in a relativistic gauge plasma q is typically of order gT , g being the QCD coupling and T the temperature), one can expand the collision integrand in powers of the momentum exchange. Truncating the expansion to second order corresponds to the Fokker-Planck (FP) approximation, which, for a homogeneous and isotropic system, gives The study of the heavy-quark propagation in the medium is then reduced to the evaluation of three transport coefficients expressing the friction -A(p) -and the momentumbroadening along the transverse and longitudinal directions -B 0/1 (p) -suffered in the plasma. Actually, since one must enforce the asymptotic approach to thermal equilibrium with the medium, the above coefficients (in principle all derived from the scattering matrix) cannot be taken as independent, but are related by the Einstein fluctuation-dissipation relation which establishes a link between the momentum broadening and the friction force felt by the heavy quark (d being the number of spatial dimensions). Our choice, in the POWLANG setup, is to evaluate B 0 (p) and B 1 (p) from first principles and to get A(p) from Eq. (2.3).
In order to embed the study of the heavy-quark transport into a setup including the simulation of the initial QQ production through a pQCD event-generator and the modeling of the evolution of the background medium through a hydrodynamic calculation, it is more convenient to rephrase the FP equation in the form of a discretized Langevin equation: (2.4) One no longer deals with the time evolution of a phase-space distribution but rather with the one of a (large) sample of relativistic particles. Eq. (2.4) provides a recipe to update the heavy quark momentum in the time-step ∆t through the sum of a deterministic friction force and a random noise term specified by its temporal correlator It can be shown that there is a one-to-one correspondence between the transport coefficients entering into the Langevin equation and the FP ones: κ ⊥ (p) = 2B 0 (p) and κ (p) = 2B 1 (p). Concerning the friction term, the momentum-dependence of the noise-noise correlator (multiplicative noise) requires to consider carefully the discretization of the equation. In the pre-point Ito scheme, in updating the heavy-quark position and momentum during the time-step t → t + ∆t, the transport coefficients are evaluated at time t and one can show that, in this case, the friction term coincides with the FP one, η D (p) = A(p), given in Eq. (2.3). Other schemes are sometimes employed in the literature: for an overview we refer the reader to Ref. [26]. Since the heavy quarks propagate throughout an expanding fireball, the evaluation of the transport coefficients and the update of their momentum has to be performed at each time-step in the local rest-frame of the fluid, eventually boosting back the result to the laboratory frame [9], in which the medium flows with four-velocity u µ . If the coupling with the background medium were sufficiently strong, the heavy quarks would tend to approach kinetic equilibrium with the plasma in its local rest frame and hence, after boosting to the laboratory frame, to share its collective hydrodynamic flow. Within the Langevin setup the interaction between the heavy quark and the medium is summarized (thanks to the Einstein relation) in just two transport coefficients, κ ⊥ and κ , which reduces to a single one, κ, in the non relativistic limit. Theoretical calculations for κ in hot-QCD exist in the M → ∞ static-quark limit. Lattice-QCD calculations for the case of a gluon plasma were performed, for various temperatures, in [27] and recently first continuum-extrapolated results have become available [28], although extracting real-time information from simulations in a Euclidean spacetime introduces large systematic uncertainties. Furthermore, NLO analytic weak-coupling calculations for κ were performed in [29], introducing large positive corrections with respect to the tree-level result. Unfortunately the kinematic range in which the most solid theoretical results for κ are (so far) available is not the one of relevance for describing (or extracting information from) the experimental data, referring mainly to heavy-flavour particles in a relativistic regime. Hence, in our simulations with weak-coupling transport coefficients, we have to account for their full momentum dependence, evaluating κ ⊥ (p) and κ (p) within a tree-level calculation with Hard Thermal Loop (HTL) resummation of medium effects in the case of interactions mediated by the exchange of soft gluons. In the case of lattice-QCD transport coefficients, on the other hand, no information on the momentum dependence is available and we simply take κ from the static calculation of Ref. [27], which covers the largest range of temperatures. Actually, an alternative strategy could consist in exploiting the experimental data on heavy-flavour observables (e.g. the nuclear modification factor and the elliptic flow) to estimate a posteriori the most probable value of the heavy quark transport coefficients. This was done for charm through a Bayesian analysis in [30], obtaining results compatible with lattice-QCD calculations.

Modeling of the background medium
In order to simulate the heavy quark transport in the fireball produced in heavy-ion collisions one needs to model the initial conditions and the subsequent hydrodynamic expansion of the background medium. The initial state is simply taken from the Glauber model, either in its optical or Monte Carlo implementation. As in our past publications [9][10][11][12], the system is initialized via the entropy-density at the longitudinal proper-time τ 0 ranging, depending on the center-of-mass energy of the collision, from GeV to τ 0 = 0.5 fm/c at √ s NN = 5.02 TeV. The hydrodynamic equations describing its evolution are solved through the ECHO-QGP code [31] in 2+1 dimensions, assuming longitudinal boost-invariance, which is a reasonable approximation to describe observables around mid-rapidity 1 .
s(x,y) (fm -3 ) 0-10% Pb-Pb coll. In the case of a smooth optical-Glauber initialization, as described in [9][10][11], the entropy density at τ 0 is taken as proportional to the local density of binary nucleon-nucleon s(x,y) (fm -3 ) 0-10% Pb-Pb coll.  For observables like the nuclear modification factor and the elliptic flow in non-central nucleus-nucleus collisions the optical-Glauber model is sufficient to capture the relevant features of the initial conditions driving the medium evolution. On the other hand, for the study of observables arising from event-by-event fluctuations of the initial geometry like the triangular flow, this is not enough: smooth initial conditions would lead to v 3 = 0 for any impact parameter of the colliding nuclei and only the granularity of the initial condition can give rise to a non-vanishing triangular flow. Here, as done in [12], we assume that the above lumpiness arises mainly from event-by-event fluctuations in the positions of the nucleons inside the colliding nuclei. We proceed as follows, generalizing to the nucleus-nucleus case the Monte Carlo approach adopted in [12] for proton(deuteron)-nucleus collisions. We generate several thousands (∼ 6000) of Pb-Pb collisions at random impact parameter and we organize them in centrality classes according to the number of binary nucleon-nucleon collisions. For a given event each nucleon-nucleon collision is taken as a source of entropy production, so that, employing a Gaussian smearing (with σ = 0.2 fm), we have for the initial entropy density in the transverse plane For each event the above entropy density can be used as a weight to define complex eccentricities, which characterize the initial state (i.e. both the amount of anisotropy and its orientation in the transverse plane) and are mapped into the final hadron distributions by the hydrodynamic evolution [32]: Modulus and orientation of the various azimuthal harmonics are given by: Exploiting the fact that on an event-by-event basis one has v m ∼ ǫ m for the lowest-order harmonics m = 2, 3, one can consider an average background obtained through an average of all the events of a given centrality class, each one properly rotated to have the reference angle ψ m (with m depending on the harmonic being considered) aligned along the x-axis and weighted by the number of binary nucleon-nucleon collisions (QQ production scales according to N coll ). We applied the above procedure to model the initial conditions of Pb-Pb collisions at √ s NN = 5.02 TeV, with the purpose of studying within a consistent setup both the elliptic and the triangular flow of heavy-flavour particles after their propagation throughout the medium. As in Ref. [12] the contribution to entropy production by each nucleon-nucleon collision was fixed via a matching to an optical-Glauber calculation at the same center-of-mass energy, obtaining Kτ 0 = 6.37 with an initialization time τ 0 = 0.5 fm/c. The resulting initial entropy-density profiles in the transverse plane are displayed, for different centrality classes, in Figs. 1 and 2. Notice that, being the angles ψ 2 and ψ 3 essentially uncorrelated, one gets average initial conditions displaying an almost perfect elliptic/triangular eccentricity.

Development of the heavy-quark flow
The mass-dependent flattening of the hadron p T -spectra observed in relativistic heavy-ion collisions as well as the azimuthal anisotropy of their angular distributions, parametrized in terms of various harmonic coefficients (v 2 , v 3 , v 4 ...), have been interpreted for long as signatures of the formation of a strongly interacting medium undergoing a hydrodynamic expansion which, via pressure gradients, translates the initial spatial anisotropy of the system into the final momentum distribution of the particles decoupling from the fireball (for a recent review, see e.g. [33]). More and more observables have been analyzed which can be accommodated within a hydrodynamic description like higher flow-harmonics [34], event-by-event flow fluctuations [35] and non-linear effects like interference between different flow-harmonics [36,37]. Notice that, within a kinetic description, in order for a system to behave as a fluid the mean-free-path of its constituents has to be much smaller than the system size, λ mfp ≪ L. The above condition is only marginally satisfied with perturbative partonic cross-sections and hence the idea of the formation of a strongly-interacting QGP was proposed. A further surprise came in the last few years from the observation of analogous effects (mass-dependent radial, elliptic and triangular flow) also in small systems, like the ones produced in high-multiplicity deuteron-nucleus, proton-nucleus and even proton-proton collisions [38][39][40][41]: in light of the small size of the medium this makes the hydrodynamic interpretation of the experimental measurements in these events quite challenging and alternative explanations have been proposed (see e.g. [42][43][44]). Recently some authors proposed a different paradigm to interpret the above experimental observations. Employing a transport setup with relatively mild partonic cross section of a few mb, they identified the major source of elliptic and triangular flow in their model in the anisotropic escape probability of the partons which decouple from the medium with no or very few interactions, getting a non-vanishing v 2 even in the case of small medium opacity [24]. Similar analytic estimates based on kinetic theory were performed in the past in order to explain the elliptic flow in peripheral nucleus-nucleus collisions in which one expects to produce a less dense medium [45,46]. In [24], however, the authors aim at delivering a much stronger message, suggesting that the above mechanism can account for most of the observed effect and questioning the picture of the formation of a strongly-interacting medium, with a collective flow arising from multiple collisions.  Although the above considerations mainly refer to the bulk particle production, dominated by soft, light hadrons, it is of interest to perform a similar analysis with our heavyflavour transport model, studying the decoupling of the charm quarks from the fireball (schematically assumed to occur at a temperature T dec = 155 MeV) during the various stages of its evolution and how they separately contribute to the anisotropies (elliptic and triangular) of their final (time-integrated) angular distribution. For an independent and somehow similar analysis, focused mostly on the different time-development of the heavy-flavour R AA and v 2 , see Ref. [47]. At variance with a kinetic calculation based on the Boltzmann equation, in which it is possible to keep track of the collisions suffered by each particle, in the Langevin setup the picture is more coarse-grained: in each time-step ∆t, the particle is given a random momentum kick, depending on the local value of the transport coefficients. However, it is possible to isolate the contribution to the anisotropy from the quarks decoupling at various values of the longitudinal proper-time τ ≡ √ t 2 − z 2 . The study is performed for Pb-Pb collisions at √ s NN = 5.02 TeV and Glauber-MC initial conditions, properly averaged depending on the considered flow-harmonic, as discussed in Sec. 3. In Fig. 3 we display the distribution of the decoupling time of the charm quarks in the 10-30% centrality class. Notice (as can be seen from the time-integrated red curves) that half of the quarks escape from the fireball only after a quite long time τ > ∼ 7 fm/c and this holds for both choices of transport coefficients, which give rise to very similar curves. Only a small fraction of about 10% of quarks spend in the medium a time < ∼ 4 fm/c. Hence, we expect that the interaction with the medium, in light of the average long time spent in the latter by the heavy quarks, provides a non negligible effect in determining the final angular distribution of their momenta.
Pb-Pb @ 5.02 TeV 10-30% centr. class  In Figs. 4 and 5 we display the differential contribution to the elliptic and triangular flow of charm quarks from particles decoupling at various values of the longitudinal proper-time τ . In both cases the pattern is quite similar. Quarks decoupling very early (τ < ∼ 2−3 fm/c) provide a positive contribution, interpreted as arising from the previously discussed anisotropic escape probability. For larger values of the decoupling time the situation changes, and the Fourier coefficients start to decrease with increasing τ FO , getting even negative until reaching a minimum around the time τ ≈ 7 fm/c at which most of the quarks decouple. One has then a sudden increase of the v 2 and v 3 of the heavy quarks decoupling during the latest stage, which makes the integrated final result positive. Interestingly, the picture depends only mildly on the transport coefficients (weak-coupling HTL  or non-perturbative l-QCD) employed. The peculiar behaviour of flow development can be interpreted also in light of the freeze-out τ − r correlation plotted in the right panel of Fig. 3, in which one can clearly identify two bands -corresponding to heavy quarks decoupling along the x and y-axis respectively -which at a certain value of τ cross each other: at the very latest times only quarks moving along the x-axis have still to decouple and this gives rise to the very large contribution to the v 2 seen in the left panel of Fig. 4. As it can be seen, in our framework in which one considers the heavy-quark propagation throughout a background medium undergoing an hydrodynamic expansion, the final flow signal (both for v 2 and v 3 ) is not dominated by the few particles escaping very early as in the study performed within a pure transport setup in [24], but arises from the non-trivial interplay of opposite-sign contributions from all the different decoupling times. Interestingly, as can be seen from the green curves in the right panels of Figs. 4 and 5, the trend of the v 2 and v 3 of the quarks looks in qualitative agreement with the collective elliptic and triangular flow (v fluid 2/3 ) of the fluid cells from which they decouple. In Figs. 6 and 7 we show the findings of a p T -differential study of the time-development of the elliptic flow, considering both the time-differential and integrated results, respectively. As can be seen, at early times, when the signal is dominated by the anisotropic escape probability of the partons, the p T -dependence of the effect is negligible, whereas it gets important during the later stages, where the interaction with the medium affects differently the propagation of heavy quarks of different momenta. Similar considerations hold for the dependence on the transport coefficients, HTL and l-QCD results being significantly different only at late times.

Model results versus experimental data
Here we display some new results obtained with our setup, pushing its prediction from RHIC ( √ s NN = 200 GeV) to top LHC energies ( √ s NN = 5.02 TeV) and focusing mainly on the elliptic and (for the first time) triangular flow of charm quarks and hadrons. In Fig. 8 we display the POWLANG predictions for the elliptic flow of charm quarks and hadrons in Au-Au collisions at √ s NN = 200 GeV in the 0-80% centrality class, comparing our results to STAR measurements of the D 0 meson v 2 [6]. In a previous publication [11] we already showed how the flow inherited from the light thermal partons at hadronization via in-medium recombination was able to boost the spectra of charmed hadrons towards slightly larger values of p T , leading to an enhancement of the D 0 R AA at intermediate p T not observed in the results obtained with independent vacuum fragmentation functions. results with HTL and l-QCD transport coefficients are compared to recent STAR data [6]. The right panel illustrates the effect of in-medium hadronization, which enhances the anisotropy due to the additional flow inherited from the light thermal partons.
As can be seen in Fig. 8, a similar effect occurs for the elliptic flow. In POWLANG the elliptic flow of charm quarks at the end of the partonic phase is non-negligible, but not sufficient to describe the sizable D 0 v 2 measured in the experiment. Notice that, at the quark level, results obtained with weak-coupling HTL and non-perturbative l-QCD transport coefficients differ substantially, the l-QCD curve displaying a much larger v 2 at low p T due to the larger value of the momentum diffusion coefficient (this can be also appreciated comparing the left and right panels of Fig. 7), the HTL curve saturating instead at a larger value of v 2 at high p T , simply reflecting the different amount of parton energy-loss in-plane versus out-of-plane, larger in the HTL case due to the steep rise of κ (p). On the other hand, after hadronization via recombination with light thermal partons feeling the collective expansion of the medium, the v 2 of charmed hadrons turns out to increase at low-moderate p T and looks in better agreement with the experimental data. In Figs. 9 and 10 we consider POWLANG predictions for the elliptic and triangular flow of charmed hadrons in Pb-Pb collisions at √ s NN = 5.02 TeV. As in the previous case, differences between HTL and l-QCD results are more evident at high p T , due to the different momentum dependence of the transport coefficients, in particular of κ . For the background medium we employ Glauber-MC initial conditions, taking, as explained in detail in Sec. 3, a proper weighted average of hundreds of collisions belonging to the same centrality class. The agreement with the D-meson v 2 and v 3 values measured by the ALICE [48] and CMS [7] collaborations is quite good. As in the case of light hadrons, the triangular flow does not arise from the finite impact parameter of the collisions (in fact the signal does not change so much in the different centrality classes) but is due to event-by-event fluctuations: in our study we limited ourselves to geometric fluctuations in the nucleon positions.
Also in this case we can disentangle in the model the effect of the heavy-quark trans- TeV, for various centrality classes. POWLANG predictions with HTL and l-QCD transport coefficients are compared to ALICE [48] and CMS data [7]. TeV, for various centrality classes. The signal arises from event-by-event fluctuations in the initial conditions. POWLANG predictions with HTL and l-QCD transport coefficients are compared to CMS data [7]. port through the deconfined plasma and of the in-medium hadronization, both for the elliptic and the triangular flow; the results are shown in Fig. 11. Similar considerations to what found at √ s NN = 200 GeV hold. The additional flow acquired at hadronization via recombination with light thermal partons enhances the D-meson anisotropies at low and intermediate p T , moving the curves closer to the experimental data. Also in this case, at the partonic level, the flow signal at low p T obtained with l-QCD transport coefficients is larger due to the stronger coupling with the medium; notice however that in this kinematic range hadronization tends to wash out the differences arising from the transport coefficients. On the other hand, at high p T larger v 2 and v 3 values are obtained in the HTL case, reflecting the larger energy-loss at high momentum.
In Fig. 12 we address the elliptic flow of electrons from heavy-flavour semi-leptonic decays in the case of non-central Pb-Pb collisions at √ s NN = 5.02 TeV. Curves with HTL and lattice-QCD transport coefficients look very similar at low momentum, while they display a quite different behaviour at higher p T arising mainly from the one of the parent As usual, it is important to check that the same transport setup provides a consistent description not only of the azimuthal anisotropies of heavy-flavour hadron distributions, but also of the medium modifications of their p T -spectra, reflecting, depending on the kinematic region, either the radial flow (dominant at low-moderate p T ) or the energy loss (the relevant effect at high-p T ) acquired/suffered by the heavy particles. Hence, in Figs. 13 and 14 we display the POWLANG predictions for the R AA of charmed hadrons (and parent quarks) in Pb-Pb collisions at √ s NN = 5.02 TeV for various centrality classes, from central to peripheral ones. In the 0-10% centrality class our transport results are compared to experimental measurements of the nuclear modification factor of D 0 mesons performed by the CMS collaboration [49]. Transport results are characterized by a pronounced peak (supported also by the available experimental data) around p T ≈ 3 GeV/c, which we interpret as due to the radial flow, acquired in part crossing the deconfined medium (whose collective motion tend to boost the heavy quarks), in part at hadronization. This looks evident from the right panel of Fig. 13, in which the bump in the charm hadron R AA looks shifted to larger p T with respect to the corresponding partonic one. In POWLANG, as already discussed, hadronization is modeled through the formation of color-singlet strings/clusters obtained via recombination of the heavy quarks with light thermal partons flowing with the medium: hence, this provides a further boost to spectrum, which causes the bump to move to larger p T . Transport results at high momenta display a strong sensitity to the choice of the transport coefficients. Weak-coupling HTL results, due to the steep rise of the longitudinal momentum-broadening, tend to overpredict the amount of energy-loss. On the other hand, information on the momentum dependence of the non-perturbative lattice-QCD result for κ is missing and keeping it as a constant leads to a too small friction force acting on the heavy quarks at high momentum and hence to underestimate the energy loss. Experimental data suggest that reality sits perhaps in between these two scenarios. Finally, in Fig. 15 we extend our predictions to the electrons from semi-leptonic decays of charm and beauty hadrons, considering their nuclear modification factor in central Pb-Pb collisions at √ s NN = 5.02 TeV. We notice that up to p T ≈ 4 GeV/c the curves with weak-coupling HTL and lattice-QCD transport coefficients are very similar whereas, as usual, they tend to diverge for higher p T , due to the different momentum dependence of κ ⊥/ in the two cases. This behaviour is mainly driven by the one of the parent charm hadrons, displaying very similar momentum distributions at low p T and a quite different quenching at higher p T , as shown in the right panel of Fig. 15.

Discussion and perspectives
In this paper, besides extending to higher center-of-mass energies the predictions of our POWLANG transport model, we tried to study in more detail the development (as a function of the decoupling time from the medium) of the anisotropies in the angular distribution of heavy-flavour particles in relativistic heavy-ion collisions. We considered both the second and third harmonics of the Fourier expansion of the heavy-quark azimuthal distribution at the time of their decoupling, which may occur -depending on the initial production point -during the whole lifetime of the fireball arising from the collision of the two nuclei. The second harmonic v 2 arises mainly (except in the case of ultra-central events) from the finite impact parameter of the collision, while non-zero values of the third Fourier coefficient v 3 are entirely due to event-by-event fluctuations -in the nucleon positions, but possibly also (not addressed here) in the nucleon structure itself -giving rise to lumpy initial conditions with a non-vanishing triangular deformation. Within the hydrodynamic paradigm, an initial spatial deformation is transferred via the resulting anisotropic pressure gradients to the momentum and angular distribution of the final particles. There is quite a strong consensus in the literature that peculiar features displayed by the soft-hadrons distributions in heavyion collisions (flattening of the p T -spectra, baryon-over-meson ratios, non-vanishing value of the various azimuthal harmonics v 2 , v 3 , v 4 ...) reflect the underlying collective flow of the medium from which they decouple. Concerning heavy-flavour particles, initially produced off-equilibrium in hard pQCD processes, the measurement of non-vanishing Fourier harmonics of their azimuthal distribution may indicate that the interaction with the crossed medium was sufficiently strong to make them (at least partially) thermalize and take part in the collective flow of the latter. However, before drawing firm conclusions, it is necessary to examine whether other more trivial effects related to the collision geometry may give rise to the same kind of signals. Hence, we performed a systematic analysis of the development of the azimuthal anisotropy of particle distributions through the study of the heavy quarks decoupling from the fireball during the various stages of its expansion; both in the case of v 2 and v 3 we found a non trivial trend, the final signal arising from the interplay of very different contributions. Heavy quarks were then hadronized, through the fragmentation of color-singlet strings/clusters obtained joining them with thermal partons picked-up from the medium and hence sharing its collective flow. The final results for the elliptic and triangular flow of charmed hadrons in Au-Au and Pb-Pb collisions at RHIC and LHC energies, for various centrality classes, look in quite good agreement with recent experimental data. Actually, as already found in previous studies, the contribution provided by in-medium hadronization turns out to be quite important in moving the results at the partonic level closer to the experimental data. We also checked that our transport calculations, within the same consistent setup, provide reasonable results for the nuclear modification factor of charmed hadrons and electrons from heavy-flavour decays. Several items would deserve further investigation. Within the same centrality class (in our case identified via the minimum/maximum number of binary nucleon-nucleon collisions) one could examine the effect of eccentricity fluctuations and how they can affect the flow not only of light, but also of heavy-flavour hadrons, by selecting events characterized by a large/small eccentricity (with so-called event-shape engineering techniques). At the same time one could select events characterized by similar eccentricity, but belonging to different centrality classes. Secondly, we plan to perform transport calculations based on a full (3+1)D modeling of the background medium, dropping the assumption of longitudinal boost-invariance. This, although requiring greater storage and computing resources, will allow us to provide predictions for observables at forward/backward rapidity, so far neglected in our analysis. In particular, this will certainly provide a more realistic description of the background medium in proton-nucleus collisions, in which -due to the asymmetry of the system -the assumption of longitudinal boost-invariance is too drastic. The question of the possible hot-medium effects in small-system, also for what concerns heavy flavour, remains in fact open. In a previous publication we showed how our transport setup -with initial-state nuclear effects, partonic transport in a small QGP droplet and in-medium hadronization -provides results compatible with the experimental data, within their large systematic error-bars. Nowadays, experimental analysis with larger statistics are in progress and hopefully will provide more differential results for a wider set of observables, allowing one to put tighter constraints on theoretical models and to rule out scenarios not supported by the experimental data. We plan to address the above important items in forthcoming publications.