$D$-meson production in high energy $pA$ collisions within the QCD color dipole transverse momentum representation

The $D$-meson production is investigated by considering the unintegrated gluon distribution within the dipole approach in the momentum representation. We analyze the $D$-meson spectrum accounting for the effects of nonlinear behavior of the QCD dynamics which can be accordingly addressed in the dipole framework. The unintegrated gluon distribution is obtained by using geometric scaling property and the results are compared to the Glauber-Gribov framework. The absolute transverse momentum spectra and the nuclear modification ratios are investigated. Predictions are compared with the experimental measurements by the ALICE and LHCb Collaborations in $pA$ collisions for different rapidity bins.


I. INTRODUCTION
The heavy flavor production has been extensively studied along the last years and robust theoretical formalisms have been developed, especially with the establishment of the high energy particle accelerators, such as the Large Hadron Collider (CERN-LHC). As a consequence, the precision of the measurements in conjunction with a wide window of center-ofmass energy, transverse momentum, and rapidity distributions offer interesting prospects to investigate heavy quarks and heavy meson productions. In particular, the D-meson production may be considered as an useful source for investigating the heavy quarks and their interactions [1]. The heavy quark mass is large enough to be taken as a hard scale, allowing to evaluate the production cross sections via perturbative methods, such as the framework of perturbative Quantum Chromodynamics (pQCD) [2][3][4][5]. The D-mesons in the final state are produced by the hadronization process of these heavy quarks. Hence, the charmed meson production may carry information about the heavy flavor fragmentation function [6] as well as the partonic distribution in the nucleon or nucleus [7][8][9][10]. The relatively low mass of D-mesons allows investigations based on parton saturation effects mostly in the low-p T region. At forward rapidities the nuclear saturation scale, Q s,A , is sufficiently high and should control the suppression in the nuclear modification factor, R D pA . From theoretical point of view, the D-meson production cross section is described within the collinear factorization [11] or the k T -factorization approach [12][13][14][15]. Examples of pQCD calculations are the general-mass variable-flavour-number scheme (GM-VFNS) [16,17] and the fixed order plus next-to-leading logarithms approach (FONLL) [4,18,19]. Investigations in the context of k T -factorization framework are found in Refs. [20][21][22][23][24][25][26][27]. On the experimental side, the ALICE Collaboration at the CERN-LHC has recently reported the measurements of the azimuthal-correlation function of prompt D-mesons with charged particles [28,29] and the measurements of prompt D-meson production as a function of multiplicity [30] in pp/pA collisions at √ s = 5.02 TeV. The measurement of angular correlations is a powerful tool to investigate collective effects in a complex system as heavy-ion collisions. The measurements provided by the ALICE experiment are obtained by non-central collision events where an anisotropy is introduced in the angular distribution. In particular, saturation effects may modified the azimuthal correlation and be a window to investigate the saturation physics.
Charmed meson hadroproduction in the k T -factorization approach is calculated by con-sidering the unintegrated gluon distribution (UGD) which includes the transverse momenta of the initial partons. The UGDs can be parametrized, where their rapidity dependence, Y = ln(1/x), and gluon transverse momentum, k ⊥ , vary depending on the underlying physical assumptions. Consequently, observables strongly dependent on gluon initiated processes are crucial to constrain the UGDs in nucleons [31]. On the other hand, cold nuclear matter effects are present in collisions involving a nucleus. These effects can be associated to the large density of initial-state partons in the nucleus. The proton-nucleus reactions allow to investigate the different QCD dynamics at low x and high gluon densities [32] and offer a baseline for the analyses in heavy-ion collisions. In the high energy regime, where the processes are dominated by gluons, the nucleus can be described in terms of the Color Glass Condensate (CGC) effective theory [33][34][35] as a saturated gluonic system. Then, heavy meson production may be useful to disentangle between the scenarios based on the distinct QCD descriptions of nonlinear saturation [36,37] or collinear factorization. Cold nuclear matter effects are related to the initial-state effects and constraints may be evaluated by investigating the D-meson production in pP b collisions.
In particular, the D-meson production at small-x can be described within the color dipole formalism [38]. In such an approach, the phenomenology is based on the universal dipole cross section fitted to DIS data. The corresponding phenomenology has been proven suitable to evaluate inclusive and exclusive processes in the high energy limit. Moreover, a scaling property related to the DIS process at small-x is naturally addressed in the parton saturation framework. Namely, the geometric scaling phenomenon is traced out on the scaling property of the dipole-target scattering amplitude, N dip (x, r) → N dip (rQ s (x)). The dipole cross section accounts for the nonlinear behavior and high-order corrections of QCD dynamics [39]. The process is described in terms of a projectile that emits a gluon, which fluctuates into a quark-antiquark color dipole with definite transverse separation that interacts with the color field of the target. Namely, the associated hard process is pictured in terms of qq dipole scattering off the target. The dipole amplitude is connected to the intrinsic dipole k T -distribution, i.e., the transverse momentum distribution (TMD). In the region of large gluon transverse momentum the dipole TMD corresponds approximately to the UGD.
In this work, we implement analytical expressions for the TMDs based on gluon saturation physics in order to obtain the double differential cross section. Predictions for the D-meson production in pA collisions at the CERN-LHC regime are provided and a wide range of transverse momenta and rapidity are covered by the present analysis. The present work is an extension of studies performed for proton-proton (pp) collisions in Ref. [40]. We compute the D 0 , D + , and D * + production cross sections, the ratios of σ(D + )/σ(D 0 ) and σ(D * + )/σ(D 0 ).
Moreover, the nuclear modification factor in pP b collisions are evaluated focusing on the kinematic range available at the CERN-LHC.
The paper is organized as follows. In Sec. II the theoretical approach is presented, including the main expressions used in the calculations for the D-meson production in the color dipole framework in transverse momentum representation. In Sec. III, the predictions are shown for applying distinct analytical models for the gluon TMD, which are compared to the measurements obtained at the CERN-LHC by the ALICE and LHCb Collaborations.
The last section summarizes our main conclusions and remarks.

II. THEORETICAL FORMALISM
Considering the k T -factorization framework, the connection between the gluon UGD, , and the usual collinear gluon distribution, F g (x, µ 2 F ), is as follows, In the scenario where the momentum of the gluon in the target is particularly large, satisfying , is approximately equivalent to the UGD function times α s [41][42][43][44]. This indicates that a connection between the k ⊥ -factorization and the dipole framework is feasible, i.e., This approximated relation can be safely employed on the evaluation of the D-meson production. In general, the gluon UGD is not well determined at small k ⊥ , however there is a correspondence concerning the intrinsic dipole TMD and the color dipole cross section σ qq (see, e.g. Refs. [42,45]). Therefore, for a given dipole cross section model is possible to determine the respective TMD by taking a specific Fourier transform.
The QCD dipole formalism assumes that the production process is described via a color dipole that interacts with the color field of the nucleon/nucleus considering the target rest frame. The D-meson production is determined by the cross section of the process g + N → QQ+X, where the QQ produced in singlet and color-octet states comes from a virtual gluon fluctuation. Consequently, the cross section associated to the hadronic collision pp → QQX is given by where there is a convolution between the gp → QQX cross section and the projectile gluon UGD. Along with this, the form as the cross section is obtained in Eq. (3) is similar to that of the k T -factorization scenario. Moreover, in Eq. (3) y is the rapidity and p T is the transverse momentum of the heavy quark, while α (ᾱ = 1−α) stands for the fractional gluon momentum exchanged with the heavy quark (antiquark). In the momentum representation the heavy quark TMD is written in connection with the dipole TMD [46], with α s (µ 2 F ) being the running coupling in the one-loop approximation evaluated at the scale µ 2 (m Q is the heavy quark mass). Furthermore, x 1 and x 2 correspond to the fractional lightcone momentum of the projectile and target given by x 1,2 = M QQ √ s e ±y , with √ s denoting the collision center-of-mass energy. Additionally, the Eq. (4) contains the auxiliary quantities I i (α,ᾱ, p T ) (i = 0, 1, 2, 3) -we quote Refs. [40,46] for their corresponding expressions.
Since that the UGD is not obtained from the first principles, it requires modeling, and distinct parameterizations for the UGDs are found in the literature. In this work, we will consider the analytical parametrizations from Refs. [41,47] for the UGD in protons. The following UGD parameterization is resulting from the Golec-Biernat and Wüsthoff (GBW) dipole cross section [41,48] based on the gluon saturation assumption, assuming the form with α s = 0.2 and Q 2 s (x) = (x 0 /x) λ GeV 2 being the saturation scale in the proton. The parameters σ 0 , x 0 , and λ are extracted from a fit to the proton structure function F 2 , which was recently done in Ref. [49].
On the other hand, an approach that accounts for the geometric scaling present in charged hadron production in pp collisions combined with a Tsallis-like distribution extracted from the measured hadron spectrum is proposed in Ref. [47] (hereafter MPM model). The corresponding UGD is expressed as: with α s = 0.2 and Q 2 s (x) = (x 0 /x) 0.33 GeV 2 . The powerlike behavior of the gluons produced at high momentum spectrum is determined via the function δn = aτ b , where τ is the scaling variable defined as τ = k 2 T /Q 2 s . Moreover, the set of parameters σ 0 , x 0 , a, and b are fitted from DIS data available at small-x. Thus, we consider the parameters from Fit A in Ref. [47] for our calculuations.
An essential ingredient to obtain the differential distribution of D-meson is to account for the hadronization process of the heavy quarks. Hence, the differential distribution of open heavy mesons is given by the convolution of the heavy quark cross section and the fragmentation function, where the heavy quark light-cone momentum exchanged with the D-meson is denoted by z. In addition, D Q/D (z, µ 2 F ) represents the fragmentation function. In the calculation the KKKS parametrization [50] will be considered. Furthermore, the mass, rapidity, and transverse momentum of the D-meson are given, namely m D , Y = y , and P T , respectively [51]. The quark and charmed hadron transverse momenta are related by p T = P T /z.
The lower limits regarding the z and α integration are defined by z min = (m ⊥ / √ s)e Y and α min = (z min /z) (m 2 Q z 2 + P 2 T )/m 2 ⊥ , respectively. The D-meson transverse mass is given by m ⊥ = m 2 D + P 2 T . As shown in Ref. [40] an approximate expression for the p T -spectrum in pp collisions can be evaluated by means of the GBW parameterization. The kinematic domain established here implies that the hard scale µ F achieves higher values than the saturation scale, It can be shown that the D-meson spectra is given approximately by [40]: where z c is the average momentum fraction [50], being B c the branching fraction c → D and x cut = 0.1 [50]. We will assume z ≡ z c (2m c ).
Here we focus on the evaluation of the D-meson production spectrum in proton-nucleus (pA) collisions. Considering a heavy target colliding at high-energy regime, the nuclear QCD effects are present and, in particular, those associated to multiple parton scattering and nonlinear gluon saturation. Within the color dipole approach, such nuclear effects can be embedded onto the dipole-nucleus amplitude N A by means of the geometric scaling (GS) property derived from parton saturation models [52]. The geometric scaling (GS) is based on the assumption that the nuclear effects are absorbed into the saturation scale and on the transverse area of the colliding nucleus, establishing an A-dependence in the scattering cross section. Thus, the nuclear effects are embedded into the saturation scale and on the nucleus transverse area, S A = πR 2 A , in correlation to the proton case, S p = σ 0 /2 = πR 2 p , where the nucleus radius is given by R A ≃ 1.12A 1/3 . Therefore, it is necessary to replace the proton saturation scale by the corresponding nuclear saturation scale, consequently, where the parameters ∆ = 1/δ and S p = πR 2 p are adjusted by data producing δ = 0.79 and πR 2 p = 1.55 fm 2 [52]. Hence, the assumptions encoded in the GS are converted into the D-meson cross section which is appropriately rescaled in the following form, The approximation above has been tested against the experimental measurements for prompt photon production in pA and AA collisions in Refs. [53][54][55].
Using the approximate analytical expression for open heavy meson production off protons, Eq. (8), and the GS arguments presented above, the following parametrization for the nuclear modification factor is obtained in the limit P T → ∞: where the small-x data on γ * A collisions support an increase of Q 2 s,A stronger than A 1/3 since ∆ ≃ 1.27. Hence, at sufficiently large p T one expects R pA ≃ 1. 26. In case ∆ = 1 which corresponds to Q 2 s,A ∝ A 1/3 one has R pA = 1. The low-P T limit of the nuclear modification factor has to be determined numerically. However, the qualitative behavior should be similar to the ratio for gluon production in pA in the context of CGC approach.
Alternatively, the D-meson production in pA collisions can be computed by using the nuclear version for the unintegrated gluon distribution in Eq. (4). Namely, the UGD of the proton is substituted for the nucleus one, F A . Here we will apply the model for F A provided in Ref. [58], which is based in a Glauber-Gribov expression for the dipole-nucleus cross section, σ dA (x, r, b) at a given impact parameter b. Considering a given b and a Bessel-Fourier transform, one can associate the F A with the dipole-nucleus cross section as follows [59,60]: and by applying the technique described in Ref. [60] the UGD for the nucleus (using the GBW parametrization for σ qq (x, r)) can be written as: where B = AT A (b)σ 0 /2 and T A (b) represents the nuclear profile function. For large nucleus the series is fastly convergent and in our calculations we take n = 7. Hereafter we will refer to this model by UGDnuc.
Based on the expressions introduced before to calculate the D-meson P T and Y distributions in pA collisions, in the next section we employ the referred UGD parameterizations in order to obtain the corresponding predictions and comparing them with the experimental measurements reported by the LHC collaborations.

III. RESULTS AND DISCUSSIONS
Before investigating the D-meson production in nuclear collisions, we show the theoretical predictions that can be found in Ref. [40] compared with the current setup with LHC data in nucleon-nucleon collisions [61,62]. Comparing the LHC pp data on the P T spectrum of The results with the UGDs GBW and MPM take into account the GS phenomenon following the prescription in Eq. (12). The results will be referred hereafter as GS (GBW) and GS (MPM). Moreover, we provided results by employing the nuclear UGD presented in Eq. (14) denoted as UGDnuc. Finally, the approximated expression given in Eq. (8) is labeled as GS (APPROX.).
Let us first present the predictions for the double differential cross section for D 0 production including the charge conjugated states in pP b collisions at √ s = 5 TeV. In Fig. 2 the respective results are compared to the LHCb data [63] for different forward and backward rapidity bins. Considering both forward and backward rapidities, we observe that the results with GS (MPM) describe the data except for the rapidity intervals, 3.5 < Y * < 4 and −5 < Y * < −4.5. However, GS (GBW) and GS (APPROX.) approaches reproduce similar predictions in almost the entire P T -distribution. A difference between them occurs in the spectrum at P T > 6 GeV. Also a reasonable agreement with the experimental data is observed in the very forward (3.5 < Y * < 4) and very backward (−5 < Y * < −4.5) rapidity bins. In these configurations, we can conclude that the predictions from the approximate expression mimic the estimates obtained with the GS (GBW). In particular the GBW parameterization presents a Gaussian shape that takes place in Eqs. (2) and (5), which results in a suppression that underestimates the data. Concerning the UGDnuc results in forward/backward rapidity bins the model provides a reasonable description of the experimental points in a narrow P T interval (2 < P T < 3 GeV). On the other hand, a better adjustment to the experimental measurements is found at the rapidity regions 3.5 < Y * < 4 and −5 < Y * < −4.5. It should be noticed that the results with the UGDnuc, where the parameter B < 3, are reliable [see Eq. (15)] due to the fast series convergence.
The cross section for D 0 +D 0 production in pP b collisions at 8.16 TeV is shown in Fig. 3 for forward/backward rapidity bins. Comparison is done with the preliminary measurements reported by the LHCb Collaboration [64]. Here the results indicate the same pattern found at 5 TeV, however one has a wider P T spectrum at √ s = 8.16 TeV. Also, the difference between the GS (GBW) and GS (APPROX.) results become more apparent for P T > 6 GeV   In what follows we present the results concerning the production of D 0 , D + , and D * + mesons assuming the rapidity interval −0.96 < Y * < 0.04. The theoretical predictions are shown in Fig. 4 and compared against the experimental data from the ALICE Collaboration Consequently, those ratios do not allow us to extract particular information concerning the charm quark fragmentation functions from the D 0 , D + , and D * + meson decays. Moreover, the D 0 production rate is clear to be higher than the D + and D * + ones because the respective ratios are smaller than unity.
As a final investigation, we perform an analysis considering the nuclear modification factor computed by The corresponding results for D 0 , D + , and D * + production are shown in Fig. 6 and compared to the measurements provided by the ALICE Collaboration [65]. Considering the experimental uncertainties, we notice that the results obtained with the UGDnuc e GS (MPM) approaches provide a better description than the GS (GBW) and GS (APPROX.) ones except for the D + case, where all models are in agreement with the data. However, the GS (APPROX.) predictions reproduce the ratio weakly dependent on P T as discussed in previous section. Apparently, the experimental measurements for the nuclear modification factor tend to unity, R pP b ≈ 1, suggesting that the nuclear effects have no important impact for D 0 , D + , and D * + production at midrapidity range.
It is timely to discuss the x 2 -values probed in the measured P T spectrum of D-meson production considering the kinematic range accessible at the LHCb and ALICE experiments, mainly for very forward/backward rapidity bins (see Tab. I). We can verify that the values are within the validity region of the color dipole formalism, namely x 2 ≤ 10 −2 and intermediate P T . This allows us to make feasible predictions applying such approach. One exception is the configuration of very backward rapidity bin at 5 TeV in the LHCb data, although the x 2 -value is near the validity limit.
Let us compare our calculations with other approaches in the literature. In Ref. [66] the D 0 production in pP b collisions has been addressed using a next-to-leading order pQCD calculation and a comparison between the input nuclear parton distribution functions (nPDFs) is done. The ratios R pA (y, P T ) are in good agreement with recent data using both the EPPS16 and nCTEQ15 nPDFs. Nuclear gluon shadowing at small-x is predicted and it is UGDnuc models are directly compared to the experimental data from the ALICE Collaboration [65].
argued that the description of a pure collinear approach is robust even at P T → 0. A study performing the reweighting of the two referred nPDFs using the LHC data for D, J/ψ, and Υ(1S) production in pP b collisions at LHC was done in Ref. [67]. The main conclusions in work of Ref. [67] are similar to the ones achieved in Ref. [66]. The role played by the nuclear effects driven by the fully coherent energy loss (FCEL) in cold nuclear matter on the open heavy-flavour production has been investigated in Ref. [68]. It has been demonstrated that the FCEL effects on D and B production is quite relevant and similar to those quantified in quarkonium and light hadron production. It is argued that the effect corresponds to about half of the nuclear suppression measured at the LHC at forward rapidities and low P T [68].  [69,71] concerning the forward configuration.
The experimental data are reported by the LHCb Collaboration [63].
After inspection, our results at midrapidity are compatible with those studies.
In the context of the CGC framework, Ref. [37] presents the nuclear modification factors R J/ψ pA (y, P T ) and R D pA (y, P T ) at forward rapidities. The calculations make use of the Glauber model to obtain the dipole-nucleus cross section, σ dA (x, r, b). Predictions are in good agreement with the existing data at the time and the formalism had been introduced in Refs. [69,70] for quarkonium production. A distinct procedure is employed in Ref. [71], where transverse momentum dependent multi-point Wilson line correlators are used to describe the target nucleus and proton projectile. The corresponding UGDs are obtained from  numerical results are consistent with data and there is a relevant dependence on the initial nuclear saturation scale, Q 2 s0,A , for the amount of nuclear shadowing appearing in R D pA as a function of P T (see also Ref. [72] for a compilation of results within the same framework for quarkonium and charged hadron production). In Fig. 7 we show our results for R pP b compared to those obtained by the CGC approach provided in Ref. [69] (denoted as CGC ) and in Ref. [71] [referred as CGC (FW)] taking the forward rapidity bin (2.  [71]. In our case, Q 2 s,A ≈ 3Q 2 s,p . The backward rapidity the QCD color dipole approach can be employ as well. This is associated to the x 2 probed in this particular kinematic region and the validity imposed by the approach. Considering the mean values associated to the nuclear modification factor at backward rapidity: P T = 5 GeV, Y = − 3.75 and 5 GeV of CM energy imply that x 2 ≈ 4 × 10 −2 . This is a small value for x 2 and is within the validity region of the color dipole formalism making such approach suitable to perform predictions. Still within the CGC effective theory, the absolute spectra for D 0 , D + , and D * + in pp and pA collisions have been presented in the comprehensive study of Ref. [73] (see Ref. [74] for similar studies addressing nuclear modification factors in quarkonium production). The focus in Ref. [73] is on the so-called event engineering, which is related to the spatial and momentum structure of rare parton configurations in high-energy collisions realized by changes in the system size, multiplicity, and energy. The potential of event engineered heavy flavor measurements to reveal the dynamics of these unparalleled configurations has been demonstrated.
Finally, we present the results concerning the rapidity dependence of the cross section and the nuclear modification factor of D-meson integrated over P T , 0 < P T < 10 GeV. and UGDnuc models and the data reported by the LHCb Collaboration [63] is performed.
The predictions for D 0 production at 5 TeV that are compared to the experimental points provided by the LHCb Collaboration are shown in Fig. 8. Apparently the models delivered a better description of the cross section taking the negative rapidity bin. In particular, the GBW model provides a better agreement concerning the data description at positive rapidity bin. In general, considering the rapidity distribution, the results overshoot the data. This is expected as we can see the results taking the P T spectrum, 0 < P T < 2 GeV, in Fig. 2, where the corresponding theoretical predictions overestimate the experimental data. On the other hand, the results for the nuclear modification factor show that R pP b < 1 with a decreasing towards to positive rapidity, i.e., forward configuration. This can be associated to the higher rate of D 0 production in pp collisions as we can observe from Fig. 1 and it directly affects the values of R pP b , see Eq. 16. The percentile value of the rate of D 0 production that exceeds the experimental measurements in pp collisions is approximately 30%.
For the sake of completeness, we show the rapidity dependence of the UGD parametrizations. Indeed, the x 2 dependence of the UGD allows us to investigate the rapidity dependence of the D-meson production. The results considering three values of the transverse momentum in terms of the nuclear saturation scale, k 2 T = 0.5 Q 2 s,A , k 2 T = Q 2 s,A , and k 2 T = 5 Q 2 s,A , are presented in Fig. 9. The predictions present a slightly difference starting from the central rapidity (Y = 0) and this difference becomes more pronounced as the rapidity evolution increases. In particular, the MPM result gives a larger deviation in comparison to the others parametrizations, except from k 2 T = 5 Q 2 s,A , where the deviation between the UGD predictions turn less apparent. Consequently, this difference observed with GBW, MPM, and UGDnuc parametrizations is explicitly converted into the results predict for the cross section and the nuclear modification factor of the D-meson production.

IV. SUMMARY
In this work we perform an analysis concerning the D 0 , D + , and D * + production in pP b collisions in the high-energy limit considering the kinematic region achieved at the CERN-LHC. The predictions are computed by applying the color dipole approach in transverse momentum representation with the GBW and MPM unintegrated gluon distribution including the geometric scaling property. Results for the nuclear gluon distribution based on Glauber-Gribov theory and an approximated expression valid at P T > Q s,A , Eq. (8), are presented as well. We have found that the predictions obtained with the MPM parameterization in conjunction with the geometric scaling are very consistent with measurements reported by the ALICE and LHCb Collaborations over a wide P T range in forward and backward rapidity bins. Nonetheless, we can not distinguish between the approaches by means of the ratios for D 0 , D + , and D * + production. The magnitude of nuclear effects associated to the production of D-mesons in pP b collisions seems to be small at not so small-P T at midrapidities, since the nuclear modification factor measured is consistent with unity given the experimental uncertainties.
The color dipole in transverse momentum representation offers a suitable framework to evaluate the D-meson production and effective in forthcoming investigations of D-meson measurements in heavy-ion programmes. Moreover, more data from future experimental measurements on nuclear modification factor at forward and backward rapidities are needed to further constrain the different approaches. Thus, one will be able to refine the associated phenomenology as well as the assumptions encoded in the UGDs. We restrict our investigations to analytical expressions for the UGD in protons/nuclei which parametrize the parton saturation effects. Future studies considering the numeric solutions of the nonlinear evolution equations, as the running coupling BK equation, would be valuable.