Rapidity dependence of heavy-flavour production in heavy-ion collisions within a full 3+1 transport approach: quenching, elliptic and directed flow

We extend our POWLANG transport setup for the modelling of heavy-flavour production in heavy-ion collisions to the case of full 3+1 simulations, dropping the approximation of longitudinal boost-invariance of the background medium. This enables us to provide predictions for observables for which the rapidity dependence is essential in order to obtain a non-vanishing signal, like the directed flow $v_1$, and to get reliable results also for kinematic distributions of heavy-flavour particles at forward rapidity. We compare our predictions with experimental data obtained in Au-Au and Pb-Pb collisions at RHIC and at the LHC.


Introduction
Heavy-flavour (HF) hadrons and their decay products (electrons and muons) arising from parent charm and beauty quarks are among the best probes of the medium formed in relativistic heavy-ion collisions. Due to their large mass heavy quarks (HQ's) arise from hard scattering processes occurring very early after the collision, before the formation of a thermalized plasma of gluons and light quarks (QGP) and their initial production can be calculated using pQCD. Hence, before hadronizing they cross the fireball produced in the collision probing all the stages of its evolution, performing a tomography of the latter. Furthermore, theoretical estimates indicate that the relaxation time of HQ's is of the same order as the lifetime of the QGP, so that one expects that charm and beauty quarks approach only partial thermal equilibrium with the fireball through which they propagate: deviations from full thermalization have thus the potential to put phenomenological constraints on the values of the HQ transport coefficients.
Most HF measurements concern the nuclear modification factor R AA [1,2] and the elliptic anisotropy v 2 [3,4] of the momentum and angular distributions of charmed and beauty hadrons and their decay products. This allows one to get access to the quenching of the distributions at high momentum due to in-medium parton energy-loss and to the collective (radial and elliptic) flow acquired by the heavy quarks in the expanding fireball. More recently HF studies have been extended to observables sensitive to event-by-event fluctuations of the initial geometry of the fireball, like the triangular flow v 3 [5][6][7] and the event-shape-engineering analysis of the D-meson v 2 [6,8], both satisfactorly reproduced by several transport calculations [9][10][11][12][13][14][15][16][17]. Another interesting development concerns the modification of HF hadrochemistry in nucleus-nucleus collisions (and, more generally, in high-multiplicity hadronic collisions), with a relative enhancement of the production of D s and B s mesons and Λ c baryons relative to non-strange mesons observed by the ALICE [18,19] and CMS [20,21] collaborations at the LHC and by the STAR [22][23][24] collaboration at RHIC.
So far most theoretical calculations developed for the study of HF observables are based on a (2+1)D modelling of the background medium, assuming its invariance for longitudinal boosts: most experimental measurements are in fact performed around midrapidity and this justifies such an approximation which allows one to save computing and storage resources. However, there are observables for which a full (3+1)D description of the medium is in order, either because measurements are performed at quite forward rapidity (e.g. the muons from HF decays measured by ALICE) or because their dependence on the particle rapidity is fundamental in order to get a non-vanishing signal (e.g. the directed flow v 1 , which would be zero if integrated over a symmetric interval around mid-rapidity). Hence in this paper we extend our POWLANG transport setup [25] to correctly deal with the rapidity dependence of the observables. This is done interfacing the Langevin simulation of the HQ propagation to the full (3+1)D output of the ECHO-QGP code [26] employed to model the evolution of the background medium, for which we take an initial condition with a non-trivial dependence on the space-time rapidity η s . We apply our improved setup to provide predictions for the directed flow v 1 of D mesons as a function of rapidity and for the nuclear modification factor R AA and elliptic flow v 2 of muons from semileptonic decays of charmed and beauty hadrons at forward (and central) rapidity. In particular the HF directed flow turns out to be an interesting probe for three main reasons. First of all, due to the mismatch between the initial location of the HQ's production points and the bulk matter distribution of the fireball, it allows one to perform a full three-dimensional tomography of the medium [27]: this by the way leads to a v 1 signal for D mesons much larger then the one of light hadrons, which are thermally produced at freeze-out when the expanding fluid gets converted into a system of free-streaming particles decoupling from the latter. Secondly, extending the theory-to-data comparison to the directed flow v 1 (y) of charm allows one to further constrain the HF transport coefficients. Finally, HQ'sproduced immediately after the heavy-ion collision -witness the initial huge magnetic field present in the medium (B ∼ 10 15 T). Any possible difference in the (electrically neutral) D 0 and D 0 meson distributions (in particular the directed flow v 1 (y), the magnetic field being orthogonal to the reaction plane) can only arise from the partonic stage, in which charm quarks and antiquarks carry opposite electric charge and hence feel an opposite force once embedded in an electromagnetic field. This possibility was suggested in a few seminal papers [28,29] and people looked for evidences of such an effect in the experimental data. So far a coherent picture of STAR [30] and ALICE [31] data at RHIC and at the LHC is missing. Before including the effect of the electromagnetic field, we believe it is important to correctly reproduce the size of the average (v D+D 1 ) signal, which is what we plan to do in our paper. In the next future we plan to go one step further, including the effect of the Lorentz force in the Langevin evolution of the HQ's.
Our paper is organized as follows. In Sec. 2 we discuss our modelling of the background medium through which the propagation of the HQ's takes place, focusing in particular on the initial conditions of the hydrodynamic evolution. In Sec. 3 we present the results of our transport setup for various HF observables: in Sec. 3.1 we condider the directed flow of D mesons, comparing our findings with recent results by the STAR and ALICE collaborations; in Sec. 3.2 we address the nuclear modification factor and the elliptic flow of HF particles in various rapidity intervals, comparing results in the forward and central rapidity regions. Finally, in Sec. 4 we discuss our major findings, providing an outlook of future improvements.  In order to set the initial conditions for our (3+1)D hydrodynamic equations following Ref. [32] we employ the optical Glauber model, assuming that the initial entropy deposition arises from a linear combination of the contributions from the participant ("wounded") nucleons and from the binary nucleon-nucleon collisions, with relative weight (1 − α h ) and α h , respectively. We assume that at the initial longitudinal proper time τ 0 the fireball has a finite extension in rapidity (η s ≡ 1 2 ln t+z t−z ) modelled by the function

The 3+1 hydrodynamic background
with a central flat plateau of extension 2η flat beyond which the density drops to zero with a Gaussian smearing σ η . Furthermore right (left)-moving wounded nucleons are assumed to produce more particles at forward (backward) rapidity, respectively. The effect is parametrized by the function (2. 2) The initial entropy density at τ 0 is then given by In the above equation the different modulation in rapidity of the contribution of participant nucleons from nuclei A and B is responsible for the tilting of the fireball, which eventually gives rise to the non-vanishing directed flow of the emitted hadrons. Furthermore, the tilted fireball is characterized by a sizable angular momentum and vorticity of the velocity field [33], which could explain the non-zero polarization of Λ baryons and vector mesons recently observed at RHIC [34] and at the LHC [35]. Concerning the initial condition for the fluid velocity, consistently with the pioneering study of Ref. [32], we assume at longitudinal proper time τ 0 a Bjorken-like flow with V z = z/t and V x = V y = 0, neglecting any initial transverse expansion. Hence, at the beginning of the hydrodynamic evolution, the fluid rapidity Y ≡ atanh(V z ) coincides with the spatial rapidity η s ≡ atanh(z/t). Later on the non-vanishing pressure gradient along the longitudinal axis provides an acceleration which makes Y > η s for τ > τ 0 . More complex initial flow profiles can be explored in the future.
In this paper we consider Au-Au and Pb-Pb collisions at √ s NN = 200 GeV and 5.02 TeV, respectively. In Table 1 we summarize the parameters employed to initialize the hydrodynamic evolution. In Fig. 1 we display the initial entropy density in the η s − x plane for non-central Au-Au and Pb-Pb collisions at RHIC and LHC center-of-mass energies. We notice the milder tilting of the fireball at LHC energy necessary in order to reproduce the experimental data. Employing the values given in Table 1 for the overall factor s 0 in Eq. (2.3), corresponding to the initial entropy density at the center of the fireball for a collision at zero impact parameter (b = 0), after integrating over the transverse plane at τ = τ 0 we obtain an entropy per unit rapidity dS/dη s ≈ 14600 for central (0 − 5%) Pb-Pb collisions at the LHC. This quantity, neglecting dissipative effects during the expansion, is directly related to the final rapidity density of charged hadrons measured by the detectors. Notice that, in case we assumed an exact scaling of s(x) with n coll (x) setting α h = 1 in Eq. (2.3) keeping s 0 untouched, we would get a lower value for the initial total entropy dS/dη s ≈ 12000. This is due to the fact that the distribution of binary nucleon-nucleon collisions in the transverse plane is narrower than the one of participant nucleons and hence the integration covers a smaller spatial domain. Au-Au coll. @ 200 GeV 0-5% centr. class. obtained with the ECHO-QGP code starting from the initial conditions in Table 1 are compared to experimental results obtained by the BRAHMS [36] and ALICE [37] collaborations at RHIC and at the LHC, respectively. ALICE data actually refer to the pseudo-rapidity density.   Table 1 are compared to experimental results obtained by the STAR [38] and ALICE [31] collaborations at RHIC and at the LHC, respectively. Before performing our transport simulations to describe HF observables, here we validate our hydrodynamic background against experimental data for soft particle production. We do not aim at performing a global precision fit, but simply to show that we are able to provide a reasonable description of the collective expansion of the fireball employing the parameter set in Table 1. For the kinetic decoupling temperature for soft hadrons we set the value T FO = 140 MeV, common to all centrality classes. In Fig. 2 Table 1 are compared to experimental results obtained by the PHOBOS [39] and ATLAS [40] collaborations at RHIC and at the LHC, respectively. PHOBOS data are obtained with two different methods depending on the rapidity interval: hits (h) and tracks (t). ATLAS data actually refer to √ s NN = 2.76 TeV.
the rapidity density of soft hadrons. We compare ECHO-QGP results for primary pions (i.e. excluding the feed-down contribution from resonance decay) to the experimental data obtained in central (0-5%) Au-Au and Pb-Pb collisions at RHIC and at the LHC: negative pions measured by the BRAHMS experiment (left panel, [36]) and charged hadrons measured by the ALICE collaboration (right panel, [37]). In Fig. 3 we focus on the directed flow of soft particles at RHIC and at the LHC. We compare the results for primary pions obtained with our hydrodynamic background to the directed flow of negative pions and charged hadrons measured in non-central Au-Au and Pb-Pb collisions by the STAR and ALICE collaborations. We notice that the directed flow of soft particles is larger at lower center-of-mass energy. In order to reproduce the latter, for the parameter describing the initial tilting of the fireball we have to fix the values η m = y beam − 2 = 3.36 at √ s NN = 200 GeV and η m = y beam = 8.58 at √ s NN = 5.02 TeV. Finally, in Fig. 4 we display the rapidity dependence of the elliptic flow of soft hadrons for various centrality classes. Hydrodynamic predictions for primary pions obtained with the ECHO-QGP code starting from the initial conditions summarized in Table 1 are compared to experimental data for the v 2 (η) (η being the pseudorapidity) of charged particles collected by the PHOBOS and ATLAS collaborations in Au-Au and Pb-Pb collisions at RHIC and LHC, respectively. Notice that ATLAS data refer to √ s NN = 2.76 TeV. Overall the agreement is quite satisfactory, even if we tend to overestimate the v 2 in the most peripheral collisions and to get a flatter behaviour for a larger rapidity range as compared to the experimental data, which (in particular at RHIC) tend to decrease linearly for increasing rapidity. We also slightly underestimate the elliptic flow in the most central collisions at the LHC, but this can be attributed to our optical-Glauber initial conditions, which miss event-by-event eccentricity fluctuations.
In summary, our modelling of the bulk medium produced in heavy-ion collisions allows us to get a quite satisfactory description of several soft observables and can be used as a reasonable background to study the propagation of HQ's. Our description of the fireball can be surely improved in the next future, exploring more realistic profiles for the initial longitudinal velocity of the fluid than the current Bjorken ansatz V z = z/t, adopting a MC-Glauber initialization of the density of the medium also in the study of the directed flow and including all soft hadrons and resonances in the validation against the experimental data.

Heavy-flavour transport simulations
We now address the issue of the propagation of HQ's in the tilted background described in the previous section. As discussed in detail in previous publications the transport of charm and beauty quarks in the QGP is described through the following relativistic Langevin equation [16,25] containing a deterministic friction force quantified by the drag coefficient η D and a random noise term specified by its temporal correlator Different values for the longitudinal/transverse momentum-diffusion coefficients κ and κ ⊥ are employed, arising either from a weak-coupling calculation with resummation of medium effects (HTL scheme) or from lattice-QCD simulations (l-QCD scheme). Both schemes display some shortcomings due to the mismatch between the kinematic and temperature conditions of experimental relevance and the ones in which theoretical calculations can provide results based on solid grounds. The drag coefficient η D is connected to the momentum-diffusion coefficients κ and κ ⊥ by a generalized Einstein relation, ensuring the asymptotic approach to kinetic equilibrium. Hadronization is modelled recombining with probability 1 the heavy quarks with light thermal partons from the same fluid-cell (an instantaneous decoupling at temperature T H = 155 MeV, with no further rescattering in the hadron-gas phase is assumed), forming excited Qq (or Qq) strings which are then fragmented in a 2 → 1 * → N process according to the Lund model implemented in PYTHIA 6.4. This turns out to have a major effect on the final charm and beauty hadron distributions. For further details we refer the reader to our past publications [16,25]. It is important to stress that, at variance with the background fireball, the spatial distribution of HQ's at the initial thermalization time of the bulk medium τ 0 does not display any distorsion in the η s − x plane. HQ's are in fact produced in hard scattering processes in binary nucleon-nucleon collisions, hence their initial location in the transverse plane simply follows the n coll (x) distribution (we do not consider any transverse expansion before τ 0 neither for the medium nor, for consistency, for the HQ's); concerning their longitudinal position, neglecting any interaction with the medium before τ 0 , one has simply z = v z t (v z being the HQ velocity) and hence the spatial rapidity of their fluid cell coincides with the rapidity with which they are produced, namely with no forward/backward asymmetry. Thus, in the case of HF particles, we have a second possible source of directed flow v 1 beside the anisotropic pressure gradients, i.e. the mismatch between the initial spatial distribution of the parent HQ's and the one of the bulk matter produced in the collision. In Sec. 3.1 we will investigate the quantitative consequences of such an occurrence, while in Sec. 3.2 we will profit from our (3+1)D setup to provide predictions for HF observables at forward rapidity, without the necessity of relying on the approximation of a longitudinal boost-invariant expansion. In this section we focus on the study of the directed flow v 1 of D mesons in non-central nucleus-nucleus collisions and on the information on the initial conditions and on the values of the transport coefficients that one can extract combining this measurement with the ones of radial and elliptic flow. In Fig. 5 we display the predictions for the D-meson v 1 obtained with our transport setup, comparing our results to the data collected by the STAR and ALICE collaborations in Au-Au and Pb-Pb collisions at √ s NN = 200 GeV and √ s NN = 5.02

Directed flow and HF transport coefficients
TeV, respectively. As one can see our results can qualitatively reproduce the trend of the data measured at RHIC (left panel), where the STAR experiment obtained evidence for a non-vanishing v 1 in Au-Au collision, with no significant difference between D 0 and D 0 . The situation is less clear for Pb-Pb collisions at the LHC, where we obtained a much weaker signal than at RHIC (consistent with the milder tilting of the fireball), but for which the ALICE collaboration found a significant ∆v 1 between D 0 and D 0 mesons, although the results are affected by large statistical and systematic uncertainties which prevent one from drawing firm conclusions. Since neither D 0 nor D 0 mesons carry electric charge and hence their interaction with a hot hadron gas should be the same even in the presence of a non-vanishing electromagnetic field, the only conceivable origin of a different behaviour should be looked for in the response of the parent charm quarks/antiquarks to the strong primordial magnetic field generated by the spectator protons and partially frozen in the fireball if the QGP is characterized by a large electric conductivity. Addressing such an issue and understanding the origin of the diffent size of ∆v 1 observed at RHIC and at the LHC would require coupling our transport setup to a full set of RMHD equations, providing a consistent description of the evolution of the matter and of the electromagnetic fields. We leave this for future work, waiting also for experimental data affected by smaller statistical and systematic uncertainties. In the following our results always refer to the average of particle and antiparticle contributions. As above discussed, the initial spatial distribution of the HQ's not matching the density of the bulk medium provides a second possible contribution to the D-meson v 1 , besides the directed flow of the background arising from the initial tilting of the fireball: this may lead to a larger v 1 for D mesons than for light hadrons. We perform such a comparison in Fig. 6. As one can see, the directed flow of charmed hadrons quantified by the slope |dv 1 /dy| turns out to be about one order of magnitude larger than the one of primary pions; both at RHIC and at LHC energies one gets a larger v 1 for D mesons with lattice-QCD than HTL transport coefficients.
Hadronization was found to significantly affect the final momentum and angular distributions of charmed and beauty hadrons. Recombination of HQ's with light thermal partons makes HF hadrons inherit part of the collective flow of the medium, which turns out to be important in order to get values of the v 2 and v 3 coefficients in agreement with the exper- imental data. In Fig. 7 we show how also the directed flow is affected by hadronization. For p T < ∼ 3 GeV/c the recombination of charm quarks with light thermal partons from the same fluid cell (which in POWLANG gives rise to the formation of colour-singlet strings to further fragment) enhances the v 1 (y) of the final D mesons. If one considers the very low value of the v 1 of charged hadrons this may appear surprising, since the overall momentum distribution of light quarks should accordingly display a very small dipole asymmetry. However, HQ's do not perform a uniform sampling of the hadronization hypersurface at T = T c and of the momenta of the thermal partons with which they recombine: there is a strong correlation between the momentum of the HQ's and the position and velocity of the fluid cell in which they are found at T = T c and hence the momentum of the light quarks with which recombination occurs. This explains why in any case the directed flow of D mesons is larger than the one of the parent charm quarks.
To what extent adding as an observable the directed flow v 1 (y) allows one to put tighter constraints on the value of the HF transport coefficients in the QGP? Clearly, with the same set of coefficients one must reproduce both the momentum distribution -and specifically its nuclear modification factor R AA -and the various azimuthal flow harmonics v 2 , v 3 (the latter not addressed here, since it would require dealing with event-by-event eccentricity fluctuations) and also v 1 . It is interesting to notice that the twofold origin of the directed flow of charm hadrons leads to a non-trivial behaviour. In fact, one expects that as the HQ momentum-diffusion coefficient κ increases HF particles should approach local thermal equilibrium with the fireball and hence their final momentum distribution should coincide with the one predicted by hydrodynamics and provided by a standard Cooper-Frye algorithm describing sudden particle decoupling from a freeze-out hypersurface. Starting   Figure 9. The same as in Fig. 8, but for charm quarks. For the elliptic flow transport results for very large values of κ are even closer to the hydrodynamic limit. from the lattice-QCD result for κ and multiplying the latter by a factor 5 and 10, in Figs. 8 and 9 we display the resulting R AA (left panels), v 2 (middle panels) and v 1 (right panels) of charm hadrons and quarks, respectively. We notice that for such large values of κ the results of transport calculations look very similar up to moderate values of p T . However, the situation is very different for elliptic and directed flow. In the case of v 2 transport results approach the Cooper-Frye hydrodynamic limit, also shown for comparison in Figs. 8 and 9 (in the case of charm quarks the matching is perfect up to p T ≈ 4 GeV/c). On the other hand, in the case of v 1 , transport calculations at strong coupling lead to a completely different behaviour, with a much larger directed flow than the one predicted by hydrodynamics. This may appear an inconsistency of the model, but it admits a quite natural explanation. As already mentioned, the initial spatial distribution of the QQ pairs is completely fixed by the position of the primordial nucleon-nucleon collisions in the transverse plane and by the momentum with which the HQ's are produced. In particular, at variance with the density of the background medium, the QQ distribution does not display any spatial tilting. Hence, at a given rapidity, there is a mismatch between the center-of-mass of the fireball and the one of the QQ distribution and this represents an important contribution to the final v 1 . Furthermore, one has for the spatial diffusion coefficient D s = 2T 2 /κ: the larger κ, the smaller D s . Since x 2 ∼ 6D s t, for very large κ the HQ's undergo a very limited spatial diffusion and the asymmetry between their distribution and the one of the bulk medium tends to survive for the entire hydrodynamic expansion of the fireball.   Finally, in Fig. 10 we address the p T -distribution of HF particles, both c quarks (left panel) and charm hadrons (right panel). Also in this case transport outcomes for very large values of κ tend to converge to the same result up to moderate p T , but for the D mesons the difference from the pure hydrodynamic prediction is slightly more pronounced than for charm quarks. The initial HQ distribution, following n coll (x), is narrower than the one of the background medium, whose initial density here follows a linear combination of n coll (x) and n part (x). Due to the small value of D s the HQ spatial diffusion is limited, this HQ-bulk asymmetry survives and at decoupling HF particles perform a non-uniform sampling of the freeze-out hypersurface, with an overpopulation of HQ's decoupling from fluid cells closer to the origin and hence recombining with thermal partons with smaller collective radial velocity.

Heavy-flavour at forward rapidity
Working with a full (3+1)D hydrodynamic background opens the possibility of addressing the study of HF observables at forward rapidity, relying on a more realistic description of the hot QCD medium than the one based on the assumption of longitudinal boost-invariance.
Measurements of open HF at forward rapidity were carried out so far only indirectly, via the detection of muons from the decay of charm and beauty hadrons through the muon spectrometer of the ALICE experiment. Results for the nuclear modification factor and the elliptic flow of HF decay muons were presented in Refs. [41][42][43] and compared to the predictions of some transport calculations [11,44,45]. Unfortunately, so far it is not possible to separate the contributions from charm and beauty semileptonic decays, which will become accessible in Run 3 at the LHC thanks to the new Muon Forward Tracker [46]. Important opportunities in Run 3 will arise also from the LHCb experiment, which after the upgrade of the detector will have the potential to measure D and B mesons at formard rapidity. The second and third azimuthal harmonics of HF decay muons were also studied by the ATLAS experiment in a quite broad window |η| < 2 around mid-rapidity [7] and compared to various model predictions [47,48]. The quite broad rapidity range covered by ATLAS measurements makes of interest having at disposal a theoretical setup non relying on the simplifying assumption of longitudinal boost-invariance. 60-80% |y|<1 2<y<4 Figure 11. The nuclear modification factor R AA (p T ) of charmed hadrons in Pb-Pb collisions for various centrality classes at mid (|y| < 1) and forward (2 < |y| < 4) rapidity. Both for HTL (left panel) and l-QCD (right panel) transport coefficients POWLANG results display only a mild dependence on the particle rapidity.
We start considering the nuclear modification factor of charm and beauty hadrons in Pb-Pb collisions at √ s NN = 5.02 TeV , shown in Fig. 11 and 12 for different centrality classes. As one can see, results at central |y| < 1 and forward 2 < y < 4 rapidity are very similar. On the contrary, our predictions are quite sensitive to the choice of transport coefficients, in particular for the case of charm. The momentum dependence included in the HTL weak coupling transport coefficients leads to a strong quenching of the spectra at high p T . On the contrary the R AA obtained with l-QCD transport coefficients display a milder quenching of the production of high-momentum HF particles. In both cases the R AA is characterized by an enhancement at moderate p T , reflecting the radial flow of HF hadrons, partially inherited from the recombination with light partons at hadronization. We now move to the elliptic flow. Also in this case, both for charm and for beauty,  60-80% |y|<1 2<y<4 Figure 12. The same as in Fig. 11 but for B-mesons. 60-80% |y|<1 2<y<4 Figure 13. The elliptic flow v 2 (p T ) of charmed hadrons in Pb-Pb collisions for various centrality classes at mid (|y| < 1) and forward (2 < |y| < 4) rapidity. Both for HTL (left panel) and l-QCD (right panel) transport coefficients POWLANG results display only a mild dependence on the particle rapidity.
the rapidity dependence of our results is quite mild for all centralities in the considered rapidity intervals |y| < 1 and 2 < |y| < 4. This is consistent with what found for our hydrodynamic background and displayed in the right panel of Fig. 4. We expect that at more forward rapidity the v 2 coefficient should actually decrease. Concerning charm hadrons, shown in Fig. 13, although the magnitude of the v 2 is similar both for HTL and l-QCD transport coefficients, in the last case the v 2 displays a more rapid decrease with  10%  |y|<1  2<y<4   10-30%  |y|<1  2<y<4   20-40%  |y|<1  2<y<4   30-50%  |y|<1  2<y<4 60-80% |y|<1 2<y<4 Figure 14. The same as in Fig. 13 but for B-mesons.
p T . This is consistent with what found for the R AA . Actually, at high p T a positive v 2 does not reflect the collective flow of the medium, but simply the different energy-loss suffered by the hard partons moving along or orthogonally to the reaction plane. Hence a milder energy-loss gives rise to a milder azimuthal asymmetry of high-p T particle production.
For what concerns beauty, the larger mass of the quarks leads to a smaller azimuthal asymmetry as compared to charm. The elliptic-flow v 2 turns out to be larger in the l-QCD case, reflecting the larger value of the momentum diffusion coefficient κ entering into the Langevin equation. Notice that in the limited kinematic region considered in the figure, the rise of the HTL transport coefficients for beauty as a function of the HQ momentum is mild and not sufficient to compensate their lower value at rest with respect to the nonperturbative l-QCD estimate. 60-80% Figure 15. The nuclear modification factor of muons from charm and beauty decay in Pb-Pb collisions at forward rapidity provided by our simulations with HTL (short-dashed) and l-QCD (long-dashed) transport coefficients. Numerical findings are compared to ALICE data [41,43] .
(GeV/c)  Figure 16. The same as in Fig. 15 but for the elliptic-flow coefficient v 2 . Results are compared to the so far available ALICE data at forward rapidity at lower center-of-mass energy [42].
We finally show our predictions for the muons arising from the semileptonic decays of charm and beauty hadrons, so far the only available observable to get access to open-HF production at forward rapidity in nucleus-nucleus collisions. In Figs. 15 and 16 we show our inclusive (µ ← c + µ ← b) results for their nuclear modification factor and elliptic flow in Pb-Pb collisions at √ s NN = 5.02 TeV at various centralities. As usual, the results obtained with HTL and l-QCD transport coefficients are displayed and compared to experimental data obtained by the ALICE collaboration [41][42][43]. Medium effects look milder in the l-QCD case, but this has to be interpreted as arising from our lack of information from nonperturbative calculations on the momentum dependence of the HQ transport coefficients, which is then neglected. This is particularly relevant for high-p T leptons, since they arise from the decays of parent HF hadrons of even larger momentum. So far the ALICE experiment does not have the possibity of distinguishing the muons from charm and beauty decays. Such a more differential analysis for the azimuthal asymmetry of HF-muon production in Pb-Pb collisions was instead performed by the ATLAS collaboration [7], although limited to central rapidity, |y| < 2. Hence, in Figs. 17  we display the separate findings for the elliptic-flow coefficient v 2 of muons from charm and beauty decays obtained with our transport calculations. Their combined results are shown in Fig. 19. Since ATLAS measurements are carried out in a quite broad range around mid-rapidity, we also show the predictions of (2+1)D transport simulations, assuming longitudinal boost-invariance for the hydrodynamic background and starting from the Monte-Carlo Glauber initial conditions discussed in Ref. [16,17]. As expected, the event-by-event fluctuations accounted for by the MC-Glauber initialization make the initial eccentricity of central collisions larger than the one provided by the optical Glauber model. Furthermore in the MC-Glauber initialization employed in Refs. [16,17] the initial entropy deposition in the transverse plane is assumed to arise entirely from the binary nucleon-nucleon collisions (i.e. α h = 1) and this also contributes to make the initial eccentricity larger. The comparison of the results for the average elliptic deformation of the fireball ǫ 2 ≡ y 2 − x 2 / y 2 + x 2 at the beginning of its evolution provided by the two different initializations is shown in Table 2. In the MC-Glauber case ǫ 2 is measured in a rotated frame in which the event-plane angle Ψ 2 lies along the x-axis. The effect of the event-by-event fluctuations is particularly strong in the 0-10% centrality class. Accordingly, for a given set of transport coefficients, the elliptic flow of HF decay muons in central (0 − 10%) Pb-Pb collisions obtained with our (2+1)D transport simulations is larger than what found with the (3+1)D calculations, so far limited to an initial geometry given by the optical Glauber model. Notice that ATLAS data for HF decay muons extend to much larger transverse momenta than ALICE results. As already discussed, for such large values of p T of the muons -which arise from the decays of parent HF hadrons of even higher momentum -the non-vanishing v 2 does not reflect the collective flow of the background fireball but rather the different amount of energy loss suffered by HQ's of very high energy propagating along or orthogonally to the reaction plane. In this kinematic range -in which however also other processes like medium-induced gluon radiation should be consideredthe very different behaviour of the HTL and l-QCD transport coefficients explains the larger v 2 obtained in the HTL case.

Discussion and prospects
A full (3+1)D realistic modelling of the medium is necessary to describe several observables in relativistic heavy-ion collisions, whose dependence on rapidity has at the same time the potential to get access to a richer information on the initial state of the system and to provide additional constraints to the transport coefficients of the medium. First of all, the produced fireball has a finite longitudinal extension, with an upper bound set by the rapidity of the colliding nuclei, but with a non-negligble energy density deposited only in a more limited range in space-time rapidity around η s = 0; hence, medium effects on the final particle distributions can be different in the forward and central regions. Secondly, the fact that participant nucleons tend to deposit a fraction of their energy mostly along their direction of motion rather then backwards leads to a tilting of the initial geometry of the medium in the η s − x plane. In this connection the asymmetry of the pressure gradients in the η s − x plane leads to the development of a non-vanishing directed flow v 1 of light hadrons, which is however quite small. The non-trivial longitudinal shape of the initial conditions affects also HF observables, which are the subject of this work, in which transport equations are interfaced to the output of full (3+1)D hydrodynamic calculations. The different production mechanism of HF particles with respect to the one of the light thermal partons of the fireball has the potential to provide further insight into several features of the medium, in particular concerning its initial stage, due to the very short formation time of HQ's. Our major finding is that also D mesons develop a non-vanishing directed-flow coefficient v 1 , which turns out to be at least one order of magnitude larger than the one of light hadrons. This occurrence was found also in previous independent studies [27,28] and cannot be interpreted as simply due to the collective flow acquired through the interaction with the medium. Its origin rather lies in the mismatch between the initial spatial distribution of heavy quarks and the energy density of the medium, the latter characterized by the above mentioned tilting in the η s − x plane. Remarkably, the stronger the coupling of the heavy quarks with the medium the smaller their spatial diffusion coefficient D s is, hence such an asymmetry tends to persist throughout the evolution of the medium, since in this case each heavy quark tends to move following its original fluid cell. As a consequence, while in the limit of very strong coupling the momentum distribution and the elliptic flow of the heavy quarks tend to the Cooper-Frye result expected for particles at local thermal equilibrium, the directed-flow coefficient v 1 displays a completely different behaviour, tending to a much larger value. This is due to the fact that, although each heavy quark flows with the fluid, as a consequence of the initial mismatch between the HF and bulk-medium spatial distribution, the different cells belonging to the decoupling hypersurface are characterized by an over-population or under-population of charm quarks with respect to the case of full thermodynamic equilibrium.
Although the major motivation of our work was the study of the D-meson directed flow, we also applied our setup to model medium effects on the HF particle distributions at forward rapidity, finding, overall, only small differences with respect to the results obtained in the central region. Concerning the comparison with the experimental data, the measurements available to date are based on the detection of forward muons from charm and beautys decay performed by the ALICE experiment. The current impossibility of disentangling the two contributions will be overcome through the upgrade of the experimental apparatus. Also future measurements by the LHCb collaboration at forward rapidity will provide important information in this kinematic domain.
A possible relevant issue to investigate in the near future is the effect of the initial strong magnetic fields present in non-central heavy-ion collisions on the final distributions of heavy-quarks, produced in the very early stage in hard scattering processes. Past studies suggest [28,29] that one should look for differences in the v 1 coefficient between D 0 and D 0 mesons, which can only arise from the partonic stage of the fireball evolution. The experimental situation is not yet completely settled [30,31], with some hints of a possible splitting found by the ALICE collaboration, but with the need of reducing the current experimental uncertainties and of exploring the effect in different centrality classes. Having shown that, with our transport model, we are able to account for the rapidity dependence of several HF observables and, in particular, to reproduce the size of the average D-meson directed-flow coefficient observed by the STAR collaboration, as a next step of our study we plan to insert also the electromagnetic fields into our transport equations, coupling them to a full RMHD code that some of us contributed to develop in the past [49,50] and which provides a consistent solution of the in-medium Maxwell equations. This will allow one to get access to the initial magnetic field and, extending the study to the resistive MHD case, to the electric conductivity of the QGP. We leave the above interesting topic for a forthcoming independent publication.
A Validation of the 3+1 transport simulations   Full 3+1 transport simulations require storing a huge amount of information. Hence we performed our hydrodynamic calculations employing a grid of 101 × 101 × 101 cells in the x, y and η s directions, ranging from -16.1 to 16.1 fm in the transverse plane and from -12.1 to 12.1 in space-time rapidity. A time-step ∆τ = 0.1 fm/c in longitudinal proper time was used to store the information on the temperature and velocity of the fluid. This corresponds to a broader grid with respect to the 161 × 161 × 1 one, with time-step ∆τ = 0.04 fm/c, usually employed in our past 2D simulations. Here we verify that this does not affect our findings for observables around mid-rapidity when results of 2D and 3D simulations with the same initial conditions in the transverse plane are compared. This is clearly shown in Figs. 20 and 21 for the nuclear modification factor and elliptic flow of charm in non-central Pb-Pb collisions, respectively. The figures show also that, as long as observables integrated over a not too large symmetric interval around mid-rapidity are concerned, the non-vanishing longitudinal pressure gradients of the fireball in the 3+1 case do not significantly affect the results.