Heavy-flavour production in high-energy d-Au and p-Pb collisions

Soft-hadron measurements in high-energy collisions of small systems like p-Pb and d-Au show peculiar qualitative features (long-range rapidity correlations, flattening of the $p_T$-spectra with increasing hadron mass and centrality, non-vanishing Fourier harmonics in the azimuthal particle distributions) suggestive of the formation of a strongly-interacting medium displaying a collective behaviour, with a hydrodynamic flow as a response to the pressure gradients in the initial conditions. Hard observables (high-$p_T$ jet and hadron spectra) on the other hand, within the current large systematic uncertainties, appear only midly modified with the respect to the benchmark case of minimum-bias p-p collisions. What should one expect for heavy-flavour particles, initially produced in hard processes but tending, in the nucleus-nucleus case, to approach kinetic equilibrium with the rest of the medium? This is the issue we address in the present study, showing how the current experimental findings are compatible with a picture in which the formation of a hot medium even in proton-nucleus collisions modifies the propagation and hadronization of heavy-flavour particles.


Introduction
One of the most surprising findings in the experimental search for the Quark-Gluon Plasma is certainly the signature of possible collective effects, suggestive of the formation of a hot strongly-interacting medium, recently observed in the collisions of small systems like p-Pb at the LHC and d-Au (and now also 3 He-Au) at RHIC, in particular when selecting events characterized by a high multiplicity of produced particles. Various observables support the above picture: the structure of two-particle correlations in the ∆η−∆φ plane (double ridge), suggestive of a boost-invariant initial condition, with azimuthal spatial asymmetries mapped by the strong interactions into the final particle spectra [1][2][3][4]; the non-vanishing values of the elliptic, triangular and higher flow-harmonics [5,6], obtained also through the study of higher-order cumulants [7], which seem to indicate a common correlation of all the particles with the same symmetry-plane; the hardening of the p T -spectra moving towards more central events [8,9], which can be described as the effect of the collective radial flow of an expanding medium. The above effects display also a characteristic dependence on the particle species (mass ordering), still in agreement with the expectations of a hydrodynamic description [10][11][12]. Theory predictions based on such a picture can be found e.g. in Refs. [13][14][15][16][17].
The possible formation of a medium featuring a collective behaviour in proton-nucleus collisions, on the other hand, does not seem to significantly affect the yield and momentum distribution of hard observables like jets and high-p T hadrons: minor changes with respect to the p-p case (after accounting for the proper scaling with the number of binary nucleonnucleon collisions) can be attributed to initial-state effects, like the nuclear modification of the Parton Distribution Functions (PDF's). The nuclear modifications factor R pPb of jets [18,19] and charged hadrons [20] at the LHC was found to be compatible with unity within the experimental error bars, although these findings have to be taken with a grain of salt, due to the absence of a p-p benchmark at the same center-of-mass energy which introduces large systematic uncertainties: alternative interpolations of p-p data taken at lower and higher energies can lead to different results [21]. These findings are not necessarily in contradiction with what observed in the soft sector, since in-medium parton energy-loss -playing a major role in heavy-ion (A-A) collisions -has a strong dependence on the the temperature and the size of the medium (for the average energy-loss due to coherent gluon radiation one has ∆E ∼qL 2 ∼ T 3 L 2 [22]), while on the contrary for the radial velocity of the medium at time t one finds the approximate behaviour v x/y ∼ c 2 s t/σ x/y [23]: in this case the speed of sound c s has only a mild temperature dependence and the smaller transverse size σ x/y of the system in p-A with respect to A-A leads to larger pressure gradients and hence to a larger radial flow.
In light of the above findings it is clearly of interest what happens to heavy flavour (HF) particles in such small systems: in fact, due to their large mass, the initial charm and beauty production occurs in hard pQCD processes on a very short time-scale, like all other hard particles; at the same time, however, experimental data show that in A-A collisions they tend to acquire at least part of the elliptic and radial flow of the medium. Is the hot medium possibly formed in p-A collisions, although of small size and short lifetime, able to leave its signatures in the final hadronic observables also in the heavy-flavour sector? For several years the paradigm was that such collisions were only useful to point out initial-state effects, like nuclear modifications of the PDF's. First experimental data, due to their large systematic uncertainties, couldn't rule out such an interpretation arising from the above theoretical prejudice. However, as more experimental data are getting accessible, it is clearly of interest to use them to discriminate among the different theoretical scenarios, including or not final-state effects. The current experimental situations is the following. Results for the nuclear modification factor of the spectra of various heavy-flavour particles have been obtained for different colliding systems and center-of-mass energies: non-photonic electrons (NPE's) from charm and beauty-hadron decays in d-Au collisions at √ s NN = 200 GeV by PHENIX [24], D mesons and HF electrons by ALICE [25,26], highp T B mesons and b-jets by CMS [27,28] and J/ψ's from B decays by LHCb [29] in p-Pb collisions at √ s NN = 2.76 TeV. Overall, one can claim that -within the large systematic uncertainties -the measured R pPb /R dAu is compatible with unity, i.e. no final-state effects, although PHENIX data lie systematically above unity. Besides the p T -spectra, it is clearly of interest to study how HF particles are distributed in the azimuthal plane and how they are correlated among themselves and with the other particles produced in the collision, with a double purpose: both for looking for possible medium-modifications of the original Q−Q correlations from the initial hard production and for checking whether correlations with the other hadrons reflect a common correlation with the same symmetry plane characterizing the initial condition, i.e. whether also HF particles tend to follow the flow of the medium (in case the latter is created). Due to the small branching ratio (∼4% for D 0 → Kπ) a direct measurement of D−D correlations is currently out of reach. Preliminary results from ALICE for D-h and e-h correlations in p-Pb collisions are available [30,31]. However, if the purpose is to display the possible angular decorrelation of the original QQ pairs from the hard event, the comparison with theoretical calculations may result difficult, since part of the away-side hadrons can come from the fragmentation of light jets from NLO contributions to heavy-quark production, like flavour-excitation or gluon-splitting processes. A more direct link with the parent heavy-quarks is provided by PHENIX measurements of e-µ correlations [32], both leptons (after background subtraction) having a heavy quark as an ancestor; comparing p-p and d-Au measurements one observes a suppression of the awayside yields, suggesting a decorrelation due to the interaction with a medium. On the other hand, as above mentioned, correlations of heavy-flavour particles with the other charged hadrons from a given event may be studied in order to check whether they display the same long-range structure and azimuthal modulation observed in the soft sector and interpreted as arising from the elliptic (and possible higher harmonics) flow of a strongly-interacting medium. Besides the above mentioned preliminary study of e-h correlations [31], the ALICE collaboration has recently presented results for forward-central µ-h correlations in p-Pb collisions [33]. In central events, if the jet-like contribution estimated from peripheral collisions is subtracted, a double ridge structure reminiscent of what found for light hadrons appears: this can be interpreted as a signal of elliptic flow of the muons, part of which coming from charm and beauty decays.
To summarize: even if strong conclusions cannot be drawn, there are hints (although not an evidence) from recent experimental data that heavy-flavour particles produced in the collisions of small systems may be characterized by a finite elliptic flow, signature of the rescattering with a strongly-interacting medium. Can this lead one to reconsider the interpretation of the results for the nuclear modification factor? Can the statement that the latter looks compatible with unity arise from the theoretical prejudice that no medium can be produced in p-A or d-A events and, in case it were produced, it would lead to a quenching of the spectra? An interesting analysis was carried out in Ref. [34], where blastwave spectra for D and B mesons -with parameters fixed using light hadron spectraturned out to be able to explain the R dAu < ∼ 1.5 of the HF decay electrons in d-Au collisions at √ s NN = 200 GeV, the value larger than unity at moderate p T being interpreted as due to the radial flow acquired in the medium. Although the above interpretation looks suggestive, it it nevertheless based on a quite simplified picture, assuming that heavy mesons follow the flow of the other hadrons, without wondering whether the latter is a realistic assumption.
Here we wish to provide a more solid theoretical framework to study HF observables in the collisions of small systems, focusing on the d-Au and p-Pb cases at RHIC and LHC.
Here, a proper hydrodynamic background for the heavy-quark propagation is developed (some preliminary results were shown in [35,36]). Initial-state nuclear effects (nPDF's and transverse-momentum broadening) are included in the initial hard production of the QQ pairs, which are then distributed in the transverse plane according to the local density of binary collisions. Assuming that, within a quite short time interval τ 0 , a thermalized medium is formed, which will live for about 2-3 fm/c in the deconfined phase, one can follow the evolution of charm and beauty quarks in such a hydrodynamic background by solving the relativistic Langevin equation developed in Refs. [37][38][39]. Finally, heavy quarks are hadronized according to a mechanism [40] which involves their recombination with light partons from the medium to form colour-singlet objects (strings), eventually fragmented to produce the final hadrons. The results of the above setup, accounting for the heavy-quark propagation and hadronization in the presence of a hot deconfined medium, is that final HF particles (D and B mesons and their decay electrons) are characterized by a non-vanishing radial and elliptic flow, providing a new paradigm to interpret current experimental data.

The setup
In A-A collisions the main source of initial eccentricity (at least for what concerns elliptical deformations) is represented by the finite impact parameter of the two nuclei. In this case, for the initialization of the hydrodynamic evolution of the medium, smooth average initial conditions based on the optical Glauber model are sufficient to describe the gross features of the system, which tend to develop an elliptic-flow as a response to the initial spatial anisotropy. Hence, in our previous studies focused on Au-Au and Pb-Pb collisions [38][39][40], we could rely on such a picture. The situation gets different already when one starts looking at different observables, like elliptic flow in ultra-central events or triangular flow, which can only arise from event-by-event fluctuations, not captured by the optical Glauber model. It is reasonable to assume that the major source (although not the only one) of fluctuations in the A-A case is represented by the random positions of the nucleons inside the colliding nuclei in each event. Such fluctuations can be easily simulated by standard Monte-Carlo (MC) implementations of the Glauber model [41], in which randomly distributed nucleons of the two different nuclei collide if their transverse distance d is such that d < σ in NN /π (σ in NN being the inelastic nucleon-nucleon cross-section). If the necessity of a Glauber-MC approach in the A-A case was realized only when addressing observables ignored in the first experimental studies, in the collisions of small systems (d-Au, p-Pb and also, recently, 3 He-Au) it is absolutely mandatory: initial-state anisotropies have little to do with the value of the impact parameter and are instead dominated by event-by-event fluctuations in the nucleon positions and possibly, at a more microscopic scale, in the colour-fields generated by the valence partons of the colliding hadrons. Although there are attempts in the literature to account for these sub-nucleonic fluctuations in the initialization of hydrodynamic equations [42], we neglected them, adopting a much simpler approach. As above mentioned, for the simulation of the initial conditions of the hydrodynamic evolution of the medium -which will represent the background in which the heavy quarks will be eventually distributed and made propagate -we relied on a Glauber-MC model. Each binary nucleon-nucleon collision was assumed to deposit some entropy in the transverse plane, described by a Gaussian distribution centered around the scattering position and depending on the smearing parameter σ smear . This leads to the initial entropy-density profile (2.1) In the above K is a constant, fixed linking Kτ 0 (τ 0 being the thermalization time at which the hydrodynamic evolution starts) to the final rapidity density dN ch /dη of charged hadrons (taken as a proxy of the initial entropy), according to the procedure described in Appendix A. The full set of parameters employed in the estimate of the initial entropy distribution is given in Table 1. Due to the random location of the participant nucleons and hence of the binary collisions, the above initial condition is characterized by an elliptic deformation quantified by the eccentricity ǫ 2 (a positive real number) and by its azimuthal orientation Ψ 2 [43]: where the curly brackets denote an average over the transverse plane weighted by the entropy density in Eq. (2.1). The eccentricity ǫ 2 is then given by and the event-plane angle 1 Ψ 2 ∈ [−π/2, π/2] by The eccentricity ǫ 2 coincides with the usual reaction-plane eccentricity if one rotates the nucleon positions by the angle Ψ 2 , so that the minor axis of the ellipse is aligned along the x-axis: To initialize the hydrodynamic evolution we proceeded then as follows. A few thousands nuclear configurations were generated, distributing the nucleons randomly according to the Fermi distribution with parameters given in Table 2. In the deuteron case, the relative position of the two nucleons was taken from a Hulten wave-function. A random impact parameter was then extracted from a distribution dP (b) ∼ b db. In order to increase the statistics without storing too much information, a fixed number of trials was made for each configuration and the event was kept if at least one binary nucleon-nucleon collision occurred. At the end we generated 8520 and 8260 minimum bias collisions in the d-Au and p-Pb cases, respectively, each one characterized by the entropy distribution in the transverse plane given by Eq. (2.1). According to the number of binary collisions events were divided in centrality classes. Notice that in reality, depending on the their impact parameter and on the fluctuations of colour sources, each nucleon-nucleon collision can provide a different contribution (often modeled in the literature through a negative-binomial distribution) to the final particle multiplicity N ch : experimentally, events are classified according to N ch , rather than N coll . For this work, we treated all nucleon-nucleon collisions on equal footing, neglecting such a further source of fluctuations.
Full event-by-event hydro+transport simulations of HF production would require huge storage and computing resources. For this first study we decided to adopted a simplified approach. For each centrality class considered in our analysis defined as a fraction (percentile) of the total cross section (e.g. 0-20%, 0-100%) we took the average of all the events belonging that class, each one weighted by its value of N coll (since we wish to have an average backround for HF propagation, whose production is a hard process scaling with the number of binary collisions) and rotated by Ψ 2 . Fig. 1 shows the result of such a procedure, comparing the initial profile of a single central p-Pb collision -with a quite irregular shape -to the average one for the 0-20% centrality class, the latter being characterized by a clear elliptic eccentricity. The smooth entropy-density profile thus obtained was used as the initial condition of the hydrodynamic evolution, calculated with the ECHO-QGP code [44]. We performed viscous runs with η/s = 0.08, corresponding to the universal lower bound predicted by the gauge-gravity duality. Concerning the longitudinal direction, in order to employ a reasonable amount of computing resources for this first analysis, we took a rapidity-flat profile with a Bjorken-like flow v z = z/t, reducing the calculation of the hydrodynamic background to a (2+1)D problem. The resulting temperature evolution of the medium is displayed in Fig. 2. TeV from the hydrodynamic evolution of our average initial condition, with different smearing parameters. Also shown, for comparison, are ALICE [10] and CMS [7] data obtained with 2 and 4-particle correlations.
In spite of the above simplifications, we believe to have a sufficiently realistic background for our purposes and, in particular, to be able to provide predictions for the possible elliptic flow acquired by heavy-flavour particles. This might be of relevance in order to attempt an interpretation of the recent electron-hadron and muon-hadron correlations measured in p-Pb collisions. As a validation of our hydrodynamic background, in Fig. 3 we compare the outcomes for the pion elliptic-flow in central p-Pb collisions arising from our calculations (with different smearing parameters) with the data obtained by the ALICE [10] and CMS [7] collaborations (the 4-particle cumulant analysis by CMS should remove nonflow effects). Notice also that the two measurements, based on two and four-particle correlations, differ in the sensitivity to flow fluctuations, tending to overestimate/underestimate the actual magnitude of the average v 2 , respectively. The size of the effect is approximately reproduced; The decrease of v 2 with the increase of the smearing parameter can be understood as a consequence of the lower initial eccentricity obtained with larger values of σ smear . We postpone a full event-by-event analysis (clearly desirable, due to the importance of fluctuations in such small systems), with a (3+1)D hydrodynamic evolution of the medium, to a future publication.
Having settled the background, both in the case of p-Pb and d-Au collisions, we can now address the initial QQ production, which is simulated through the POWHEG-BOX package [45]. As in our previous studies, EPS09 [46] nuclear modifications of the PDF's are adopted for the Au and Pb nuclei (only the central value is employed); on the contrary, no correction is used for the deuteron projectile. Within the POWHEG-BOX framework the calculation of the hard event is interfaced to PYTHIA, which takes care of the simulation of other processes such as initial and final-state radiation and intrinsic-k T . In the p-Pb and d-Au cases heavy quarks are assigned a further transverse-momentum broadening proportional to the average number of binary nucleon-nucleon collisions N coll in the considered centrality class. For such an additional k T -broadening in p-A events one has, ∆ 2 being the average squared-momentum acquired by the incoming partons of the proton in each individual nucleon-nucleon collision [47,48], since, on average, one-half of the collisions will occur before the hard scattering and the transverse momentum broadening acquired by the initial-state parton will be shared by the quark and the antiquark of the heavy pair; in the d-A case the above estimate is further divided by a factor 2, since the collisions involve with equal probability both nucleons of the deuteron. Heavy quarks are then distributed in the transverse plane according to the local entropy density s(τ 0 , x). Heavy-quark propagation in the QGP (assuming, as a working hypothesis that a hot deconfined medium is formed) is then simulated through the Langevin equation described at length in our previous works [38][39][40]. The latter requires the knowledge of the transport coefficients of the heavy quarks in the medium. As in our previous studies, we choose the ones provided by weak-coupling [38,39] and lattice-QCD calculations [49][50][51][52] and compare the predictions obtained in the two different scenarios. Weak coupling calculations are performed by separating hard (treated in pQCD) and soft (including HTL resummation of medium effects) collisions. Lattice-QCD results, on the other hand, refer to a non perturbative setup, in which however the quark is treated as a static colour source; hence the kinematic range in which we can rely on a solid first-principle calculation in a non-perturbative domain is quite limited, being restricted to small quark velocities. The heavy-quark stochastic dynamics is followed until they reach a fluid-cell below a decoupling temperature T d (in this study set to 155 MeV), where they are made hadronize. The kinematics of the charm and beauty hadrons is defined at this stage and their possible rescatterings in the hadronic phase are neglected.
To describe hadronization in the presence of a hot medium we adopt the model developed in a previous study by us [40], to which we refer the reader for its detailed description; here, we just summarize its main features. In the fluid-cell reached by the heavy-quark Q one extracts a light antiquark q light (up, down or strange, with relative thermal abundances dictated by the ratio m/T dec ) from a thermal momentum distribution corresponding to the temperature T dec in the Local Rest Frame (LRF) of the fluid; information on the local fluid four-velocity u µ fluid provided by hydrodynamics allows one to boost the momentum of q light from the LRF to the laboratory frame. A string is then constructed joining the endpoints given by Q and q light and is then passed to PYTHIA 6.4 [53] to simulate its fragmentation into hadrons (and their final decays). In agreement with PYTHIA, in evaluating their momentum distribution, light quarks are taken as "dressed" particles with the effective masses m u/d = 0.33 GeV and m s = 0.5 GeV.
In A-A collisions the above hadronization mechanism turns out to provide a better agreement between the results of our model and the experimental data, with respect to the employment of standard in-vacuum fragmentation functions [40]. The collective motion inherited from the light thermal parton increases the radial and elliptic flow of the final charmed and beauty hadrons, whose momentum and angular distributions display features closer to the actual experimental findings, such as the bump in the D-meson R AA at lowp T (p T ≈ 1.5 GeV/c) measured by STAR in Au-Au collisions at RHIC [54] and the sizable v 2 observed by ALICE in Pb-Pb collisions at the LHC [55]. A few words are in order concerning the relationship between the model here adopted and the standard coalescence calculations. Both of them involve some mechanism of recombination with light medium particles, which transfer their collective motion to the heavy-flavour hadrons, giving rise overall to a more satisfactory agreement with the experimental findings compared to vacuum fragmentation. Coalescence is a 2 → 1 process, occurring with high probability when the wave-function of the final D meson (if one considers charm) and the wave-packets of the two partons (the c quark and a light antiquark from the medium) display a sizable overlap; this happens when the initial partons are sufficiently close in space and have comparable velocities. The effect of heavy-quark coalescence with light partons on the final particle spectra was studied in detail for instance in [56], where it turned out to provide a better description of the R AA and v 2 of heavy-flavour decay electrons at RHIC. Its relevance for possible modifications of the heavy-flavour hadrochemistry, such as an enhanced production of D s mesons, was also pointed out in [57]. Our hadronization model, on the other hand, is based on a multistep 2 → 1 → N mechanism. One first combines a Q with a thermal q, independently of its kinematics, giving rise to a string of invariant mass M , which eventually decays into N hadrons through excitation of qq pairs from the vacuum.
In the following section we check whether the above setup, which assumes the formation of a hot deconfined medium (albeit of small size) affecting the propagation and hadronization of heavy quarks even in the collision of small systems such as d-Au and p-Pb, is able to provide a consistent description of the current experimental data.

Results: D/B-mesons and HF electrons
In this section we display the results obtained with our transport setup (referred to as POWLANG, as for the A-A case) for the production of HF hadrons and decay-electrons    in d-Au and p-Pb collisions at RHIC and LHC center-of-mass energies, respectively. We will compare them to the currently available experimental data obtained by the PHENIX and ALICE collaborations. Actually, also CMS [27] has recently obtained results for the B-meson R pPb , but so far limited to too high-p T to make a comparison with an approach based on the Langevin equation meaningful. The first experimental hints of possible final-state effects affecting heavy-flavour production in small systems were provided by the PHENIX results on the nuclear modification factor of non-photonic electrons in central d-Au collisions at √ s NN = 200 GeV [24], with central values around R dAu ≈ 1.4 over a quite extended p T -range, from 1 to 5 GeV/c. Indeed, since -in the light of A-A results -one tended to associate a possible in-medium interaction with a quenching of the spectrum (i.e. R dAu < 1) and since the systematic uncertainties from the background subtraction were large, people did not give at the beginning the proper importance to these results, considering them compatible with no mediumeffect. A study was actually carried out [34], showing that, employing blast-wave spectra for D and B-mesons with parameters fixed by light hadron spectra (i.e. assuming that they share a common flow with an expanding medium), one would have been able to explain the enhanced production of NPE's in the moderate-p T domain analyzed. Such a picture, albeit interesting, lacks a microscopic dynamical justification. Here we wish then to display the findings of the transport setup presented in the previous section, checking whether the combined effect of Langevin dynamics in the QGP and in-medium hadronization can lead to results in agreement with the experimental data.
We start considering d-Au collisions at RHIC, focusing on the 20% most central events, for which experimental data are available. In Fig. 4 we show how the formation of a hot deconfined medium in the collision affects the HF quark and hadron spectra, by modifying their propagation and subsequent hadronization. The nuclear modification factor is characterized by a bump at intermediate p T (smaller for the quarks and more pronounced for the hadrons) that we attribute to radial flow; in particular, the larger effect at the hadron level is due to the additional flow inherited from the light quarks (the recombination and subsequent string fragmentation can of course slightly smear also the final hadron rapidity with respect to the one of the parent heavy quark). The curves display the results of calculations performed with weak-coupling heavy-quark transport coefficients, characterized by a steep increase with the particle momentum, hence the sizable quenching of charm at high p T , shown in the figure up to a kinematic region out of the domain of validity of a Langevin picture. We also show the curves containing only cold nuclear matter(CNM) effects (nPDF's and k T -broadening). Notice how the effect of the Cronin broadening is quite large at intermediate p T , but, in our calculations, is washed out by the subsequent energy-loss in the deconfined medium. Therefore, the bump in the nuclear modification factor at the hadron level comes from the interplay of several effects, besides the CNM ones: low-p T charm quarks are pushed to higher momenta by their scatterings in the expanding medium, high-p T quarks tend to lose part of their energy and, finally, at hadronization the light partons transfers to the charmed hadrons part of their radial flow.
In Fig. 5 we display the predictions of our transport setup at the hadron level, both for charm and beauty, exploring different choices of the smearing of the initial entropy-density in the transverse plane and of the transport coefficients, from weak-coupling and lattice-QCD calculations. In the last case no velocity-dependence of the heavy-quark momentumdiffusion coefficient (evaluated on the lattice for the case of a static quark) was assumed, which explains why, at large p T , results stay close to unity. The important phenomenological outcome, however, is that, independently of the choice of the transport coefficients, all the curves display a bump around the same p T values, attributed -within our setup -to radial flow. GeV. Results of Langevin simulation with weak-coupling and lattice-QCD transport coefficients and different values of the smearing parameter in the Glauber-MC initialization are shown and compared to PHENIX data [24]. For comparison, we also display, in grey, the curve including only CNM effects (nPDF's and k T -broadening).
Finally, we let HF hadrons decay semi-leptonically and we study the resulting electron spectra. Results for the nuclear modification factor R dAu of HF decay electrons, corresponding to different transport coefficients and Glauber-MC initialization, are collected in Fig. 6. Due to the large systematic uncertainties all curves look compatible with the data, which are not able to discriminate among the various choices of parameters. However, they are sufficient to see that the enhancement in the R dAu of NPE's at moderate p T can be accommodated by models which, on top of cold nuclear matter effects, include also a stage of partonic transport in the hot plasma accompanied by in-medium hadronization. Therefore, within such a framework, the enhancement of the R dAu of HF decay electrons reflects the radial flow acquired by the parent D and B mesons. We also show in grey the result obtained including only CNM effects (nuclear modification of the PDF's and k T -broadening) followed by in-vacuum independent fragmentation, which tends to slightly undershoot the data in the intermediate p T region. Notice that the fact that difference with the POWLANG curves is not dramatic is not due to the absence of final-state medium effects in the last case, but to the combined effect of parton energy-loss (tending to quench the spectra at high p T ) and of in-medium hadronization (responsible for most of the final HF radial-flow): this could be already inferred from the curves displayed in Fig. 4.
We now move to p-Pb collisions at the LHC, at √ s NN = 5.02 TeV, where first experimental data for D-mesons, high-p T B-mesons, and HF decay electrons are available. Before addressing a systematic comparison with the accessible observables we wish to study sep- arately, also in this case, the effects of the transport in the QGP phase and of in-medium hadronization. In Figs. 7 and 8 we show our model predictions for the nuclear modification factor (integrated over centrality, i.e. 0-100%) and elliptic flow (in the 0-20% centrality class) of charmed quarks and hadrons. Results refer to weak-coupling transport coefficients. As one can see, part of the final radial and elliptic flow of the charmed hadrons is actually acquired at hadronization, via the recombination with light thermal partons which participate in the collective motion of the medium. The effect is particularly evident in the case of v 2 , which is small at the level of charm quarks at decoupling and reaches values around 5% for D mesons; interestingly, in this case the p T dependence looks similar to the one from a Cooper-Frye decoupling of thermalized particles.   [25]. For comparison, we also display, in grey, the curve including only CNM effects (nPDF's and k Tbroadening).
We now want to compare the outcomes of our transport calculations with the currently available experimental data. We start from the case of charmed hadrons, whose nuclear modification factor is shown in Fig. 9 after averaging over all D-meson species, since our model is not conceived to address changes in the HF hadrochemistry, which might affect exclusive measurements. POWLANG results with HTL and l-QCD transport coefficients and different values of the smearing in the Glauber-MC initial conditions are compared to ALICE data for the D-meson R pPb in the 0-100% centrality class, resulting from the average of D 0 , D + and D * + . Predictions obtained with different transport coefficients differ sizably at high-p T , due to their different dependence on the quark momentum; notice that, in particular for charm, this is a region out of the domain of validity of a Langevin approach. However, all the explored cases display some qualitative similarities: model results are characterized by a bump around p T ≈ 3 GeV/c (interpreted, as usual, as a signature of radial flow), accompained by a depletion at low-p T , partially arising from the nPDF's (shadowing) and partially from the conservation of the total number of charmed particles, moved by the rescatterings to larger p T . Also shown, in grey, is the case including only CNM effects (nuclear PDF's and k T -broadening), characterized by a monotonic rise of the nuclear modification factor with increasing p T .   Finally, in Fig. 10 we display the outcomes of our transport calculations for the spectra of electrons from charm and beauty decays e c + e b (left panel), together with the corresponding results for the isolated e b contribution (right panel). As for the previous cases, theory calculations are characterized by a systematic uncertainty from the different possible choices of initialization and transport coefficients. However, overall, they provide indications for an R pPb > ∼ 1 over a quite extended p T -range, in qualitative agreement with the current experimental data [26], also affected by large systematic error bars.

Conclusions and perspectives
In this paper we have tried to answer the question whether the possible formation of a hot deconfined medium in the collisions of small systems, like d-Au or p-Pb, suggested by several soft observables involving light hadrons, may affect also the momentum and angular distributions of heavy-flavour particles. We have shown that transport calculations, accounting for medium modifications of both the propagation and the hadronization of heavy quarks, lead to results in qualitative and quantitative agreement with the current experimental data, some of which, on the contrary, display some tension when compared to approaches including simply initial-state effects, like nPDF's. The results we got lead one to look at HF observables in small systems from a conceptually new point of view with respect to the traditional idea of considering p(d)-A collisions just as a tool to estimate cold nuclear-matter effects.
Medium effects on heavy-flavour production in small systems have been recently studied also by other authors. The approach followed in [58] is close to the one adopted for our analysis, namely a Langevin-like dynamics for the heavy quarks on top of a proper hydrodynamic evolution of the medium. On the other hand the perspective employed in [59] is quite different: modifications of the HF p T -spectra in p(d)-A collisions are attributed, in the A-going direction, to higher-twist effects such as incoherent multiple (double) parton scattering in the nuclear medium. Furthermore, ALICE data for the D-meson R pPb at mid-rapidity turned out to be well described, within the current large systematic errorbars, also by the model of Ref. [60], based on momentum broadening and energy-loss in cold nuclear matter.
The growing experimental and theoretical interest on the subject and the first results we obtained encourage us to perform more detailed studies, requiring a much heavier numerical effort, like full hydro+transport event-by-event simulations, with a more realistic (3+1)D hydrodynamic background, in order to access a wider set of observables in a larger rapidity-range. This will also allow one to discriminate among our predictions and the ones of the other models, which don't assume the formation of a hot medium, but involve only cold nuclear matter effects [59][60][61][62]. We leave such an issue for future publications.

A Initial conditions: setting the parameters
The parameters necessary to specify the Glauber-MC initial conditions for the hydrodynamic evolution in the d-Au and p-Pb cases, in particular the constant K (having dimensions of an inverse length) providing the entropy released by each nucleon-nucleon inelastic interaction, are fixed by matching the results of the optical and Monte-Carlo versions of the Glauber model in A-A collisions.
We proceed as follows. We assume the initial entropy density in a nucleus-nucleus collision at impact parameter b to be given, within the optical Glauber model, by where s 0 (in fm −3 ) sets its maximum value at the origin. Other choices are of course legitimate; we took this n coll scaling for consistency with our previous A-A studies, based on the Glauber-model initialization for the background medium of Ref. [63,64]. The entropy per unit (space) rapidity (τ 0 being the thermalization time) in the centrality-class C 1 − C 2 is obtained integrating over the transverse plane and averaging over the impact parameter: In the optical-Glauber approach N coll opt C 1 −C 2 is given by: .

(A.3)
For the sake of completeness we remind the meaning of the various terms in the above (nuclear thickness and overlap functions are assumed normalized to 1): • Probability of having at least one inelastic interaction Centrality classes are defined as percentiles of the total geometric cross-section: In the Glauber-MC approach on the other hand the entropy-density in a given event is distributed in the transverse plane according to Eq. (2.1). Its integral over the transverse plane is directly proportional to the number of binary collisions of the considered event.
The average entropy per unit rapidity in a given centrality class is now given by:  Table 1. At the LHC the situation is different, since so far we don't have Pb-Pb and p-Pb collisions at the same center-of-mass energy. We have then to extrapolate the parameters fixed at √ s NN = 2.76 TeV to the 5.02 TeV case. We proceed then as follows. The initial entropy deposited by the colliding nuclei is assumed to be proportional to the rapidity density of the final charged particles produced in the collision. Considering for instance the percentile of the most central events one has dN ch dη centr ∼ dS dη S centr = s 0 τ 0 N coll opt centr n coll (x = 0, b = 0) . (A.7) For the product s 0 τ 0 in Pb-Pb collisions at √ s NN = 2.76 TeV in our past studies [39,40] we employed the value s 0 τ 0 = 166.8 fm −2 . In order to fix the parameter K to employ in the Glauber-MC simulations at √ s NN = 5.02 TeV through Eq. (A.6) we need to estimate dN ch /dη, and then s 0 , at such a center-of-mass energy for which so far we don't have experimental data in nucleus-nucleus collisions. Hence we must rely on extrapolation of experimental data at lower values of √ s NN . In the case of A-A collisions the following empirical formula was found [65] 1 which nicely describe the evolution from RHIC to LHC, from dN ch /dη ≈ 660 in the 0-6% most central Au-Au events at √ s NN = 200 GeV [66,67] to dN ch /dη ≈ 1600 in the 0-5% centrality class for Pb-Pb collisions at √ s NN = 2.76 TeV [68]. In extrapolating the initial entropy-density from √ s 1 (2.76 TeV) to √ s 2 (5.02 TeV) we will then employ the following which leads to the estimate Kτ 0 = 6.37 quoted in Table 1.