Multiple scattering effects on heavy meson production in p+A collisions at backward rapidity

We study the incoherent multiple scattering effects on heavy meson production in the backward rapidity region of p+A collisions within the generalized high-twist factorization formalism. We calculate explicitly the double scattering contributions to the heavy meson differential cross sections by taking into account both initial-state and final-state interactions, and find that these corrections are positive. We further evaluate the nuclear modification factor for muons that come form the semi-leptonic decays of heavy flavor mesons. Phenomenological applications in d+Au collisions at a center-of-mass energy $\sqrt{s}=200$ GeV at RHIC and in p+Pb collisions at $\sqrt{s}=5.02$ TeV at the LHC are presented. We find that incoherent multiple scattering can describe rather well the observed nuclear enhancement in the intermediate $p_T$ region for such reactions.


I. INTRODUCTION
In recent years, the experimental study and theoretical understanding of nuclear effects that affect hadronic observables in proton-nucleus (p+A) collisions has attracted significant attentions [1][2][3][4][5][6][7]. On one hand, quantifying the differences between p+A and p+p collisions can provide a solid baseline for extracting the properties of the quarkgluon plasma (QGP) created in ultra-relativistic nucleus-nucleus (A+A) collisions, where both the initial-state cold nuclear matter effects and final-state hot dense medium effects modify the final-state observables [8], including heavy flavor [9][10][11][12]. On the other hand, the physics in p+A collisions is interesting in its own right in that it helps illuminate the QCD dynamics of multiple parton interactions [13], the transport properties of cold nuclear matter [14,15], the dense gluon structure of the nucleus [16], and the multi-parton correlations probed by a propagating parton in p+A collisions [17].
So far, most of the theoretical efforts have been devoted to the study of the nontrivial QCD dynamics in the forward rapidity region, where the parton momentum fraction x in the nucleus is small and the external probe interacts with the partons inside the nucleus coherently. Furthermore, the parton momentum fraction in the proton is large and the effects of energy loss from the external probe are amplified. The resulting nuclear suppression of inclusive particle production cross sections relative to the binary collisions scaled p+p baseline [18,19] has been addressed in the framework of several theoretical formalisms [20][21][22][23][24][25]. In a previous paper we explored a different regime -the backward rapidity region [26]. Specifically, we studied the single inclusive light hadron production in p+A collisions and demonstrated explicitly that in such a regime all interference Feynman diagrams drop out due to the lack of nuclear-size A 1/3 enhancement. Thus, only incoherent multiple scatterings are relevant. We adopted a generalized high-twist factorization formalism to study these incoherent multiple scattering effects. Within such a formalism, multiple parton interactions manifest themselves as power-suppressed corrections to the differential cross section and the contributions can be written in terms of high-twist multi-parton correlation functions. For recent progress on the next-to-leading order corrections within this formalism and QCD evolution of the associated high-twist correlation functions, see Refs. [27][28][29]. By taking into account both initial-state and final-state interactions, we derived the incoherent double scattering contributions to the differential cross section for single inclusive light hadron production in p+A collisions. We showed that these contributions are positive and lead to nuclear enhancement in the backward rapidity region.
In the current paper we generalize our earlier study to open heavy flavor production in p+A collisions, and use our results to understand the nuclear enhancement observed in the backward rapidity region at both RHIC and the LHC. Other approaches in understanding the backward rapidity region include the use of a universal nuclear parton distribution functions (nPDFs), see for example [30], which appear to give a somewhat unsatisfactory description of the experimental data on heavy flavor decay muons at RHIC [31]. One of the main differences in our approach is that the calculated double scattering contributions are process-dependent, that is non-universal. We find that, because of the heavy quark mass, the double scattering corrections can no longer be written in the same simple compact form as for light hadron production. However, as we will show below, such incoherent double scattering still gives a positive contribution to the heavy meson differential cross section in the backward rapidity region. We calculate the nuclear modification factor for single muons coming from heavy flavor meson decays, and find that the incoherent double scattering effects can describe rather well the corresponding RHIC and LHC data. We expect that our results will shed light on the origin of cold nuclear matter effects in the backward rapidity region.
The rest of our paper is organized as follows: in Sec. II we first review the single scattering contribution to single inclusive heavy meson production in p+A collisions. We then extend our previous calculation of the double scattering contribution for light hadron production to heavy meson production. In Sec. III, based on our analytical results, we evaluate the nuclear modification factor for muons coming form the semi-leptonic heavy flavor decays in d+Au collisions at RHIC and p+Pb collisions at the LHC. We also present comparison to the experimental data. A summary of our paper is given in Sec. IV. In this section we consider single inclusive heavy meson production in p+A collisions, where H represents the observed charm or beauty meson with momentum P h and mass m h , P ′ is the momentum of the incoming proton, and P is the momentum per nucleon in the nucleus. In general, the differential cross section for single inclusive particle production in proton-nucleus collisions can be expanded in terms of single scattering, double scattering, and even larger number of scatterings [17]: where the superscript "(S)" and "(D)" represent the contributions of single scattering and double scattering, respectively. The single scattering contribution for heavy meson production at large transverse momentum can be derived within the usual leading-twist perturbative QCD factorization formalism [32], and at leading order in the strong coupling α s it has the following form: where a,b represents the sum over all parton flavors and the center-of-mass energy squared is s = (P ′ + P ) 2 . f a/p (x ′ ) and f b/A (x) are the usual leading-twist parton distribution functions, and D c→H (z) is the fragmentation function for a parton c fragmenting into a heavy flavor meson H. H U ab→c is a short-distance hard-part function for two partons of flavor a and b to produce a parton c. It is important to realize that both heavy flavor and light flavor partons (heavy quark, light quark, and gluon) can fragment into the heavy flavor meson. Following Refs. [33,34], we take all these contributions into consideration.
The hard-part functions H U ab→c for the light flavor contribution (i.e., parton c is a light quark q or gluon g, which then fragments into the heavy meson) are well-known and can be found in Ref. [35]. On the other hand, the hard-part functions H U ab→c for heavy flavor contribution (i.e., parton c is a heavy quark Q) gets contributions from both the light quark-antiquark annihilation qq → QQ and gluon-gluon fusion gg → QQ subprocesses, as illustrated in Fig. 1, and are given by [13,36] where N c = 3 is the number of colors, and m c is the mass of the heavy quark c that fragments into the heavy flavor meson H. The slightly modified Mandelstam variablesŝ,t,û are defined at the partonic level aŝ where p c is the momentum of the heavy quark c, andŝ,t, andû satisfiesŝ +t +û = 0 as indicated in Eq. (3).

B. Double scattering contributions
Let us now study the multiple scattering contributions to charm and beauty meson production. In particular, we focus on the backward rapidity region where the parton momentum fraction in the nucleus is outside the smallx regime, and incoherent multiple interactions are important [26]. Within the high-twist collinear factorization formalism [37], the first non-trivial multiple scattering (double scattering) contributions are attributed to the twist-4 power-suppressed corrections to the differential cross section. This formalism has been used to study cold nuclear matter effects in lepton+nucleus and hadron+nucleus collisions, including parton energy loss [38][39][40][41], transverse momentum broadening [21,27,[42][43][44], and dynamical shadowing [13,20,45,46]. The detailed techniques are well explained in our previous paper [26], where we computed the incoherent double scattering contributions to the single inclusive light hadron production in p+A collisions.
As we have emphasized above, both light and heavy flavor partons can fragment into a heavy meson. For the case when the light flavor parton fragments into a heavy meson, since multiple scattering occurs at the partonic level, we can immediately obtain the double scattering contributions to the differential cross section by replacing the fragmentation function of light hadron with heavy meson, D c→h (z) → D c→H (z), in the result from our previous paper [26]. For later convenience in the phenomenological study we write down explicitly the contributions from light flavor fragmentation at twist-4: where the subscript "light" denotes the light flavor fragmentation contributions, i=I,F represents the sum over the initial-state and final-state double scattering as explained in [26] and below, and the factors c i are given by while the hard-scattering functions H i ab→cd (ŝ,t,û) have the following form: The relevant initial-state correlation functions are the so-called quark-gluon correlation function T (I) q/A (x) and gluongluon correlation function T (I) g/A (x), and they have the following definitions [26]: On the other hand, the relevant final-state correlation functions T , except for the θ-functions that are replaced as follows [21,44,47] One interesting feature of the light flavor contributions in Eq. (7) is that the second derivative, the first derivative, and the non-derivative T I,F q,g/A (x) terms share the common hard-part function, and have a very compact simple form. We will see below that such a feature will no longer hold for the heavy flavor contributions because of the mass terms.
We now turn to the double scattering contribution for the prompt heavy quark final states, which is the main new result of our analytical calculations. In other words, we study both the initial-state and final-state double scattering corrections to the partonic processes qq → QQ and gg → QQ. The relevant Feynman diagrams are shown in Figs. 2 and 3, respectively. Here, the initial-state multiple scattering represents the situations where the incoming parton from the proton undergoes multiple interactions with the soft partons inside the nucleus before the hard collisions, while final-state multiple scattering stands for the situations where the leading outgoing parton undergoes multiple interactions in the large nucleus after the hard collisions.
FIG. 2. The central-cut diagrams for initial-state (left) and final-state (right) double scatterings in quark-antiquark annihilation process. The "H"-blobs represent the hard qq → QQ processes as shown in Fig. 1(a).
Since all techniques for computing the double scattering corrections are the same as those for light hadron production, we skip all the details of the calculations and make only a few remarks about our derivation. First, besides those Feynman diagrams presented in Figs. 2 and 3 in which there is only one additional rescattered gluon in each side of the unitarity cut line in a classical double scattering picture, we also include the diagrams in which two rescattered gluons are on the same side of the unitarity cut line, representing the interferences between single and triple scatterings. These interference diagrams give the so-called "contact" contributions (i.e., have no nuclear-size A 1/3 enhancement) when combined with the central-cut diagrams in the large-x regime and will be neglected. For more details, see [26]. Because of this, our final results only depend on the four-parton correlation functions associated with the central-cut diagrams and, thus, represent the classical incoherent scattering regime. Second, the Feynman diagrams in which the rescattering happens between the unobserved outgoing parton d and the nuclear medium also leads to "contact" contributions and is also neglected. We finally obtain the following result for the double scattering contributions to the qq → QQ process, where the double scattering hard-part functions are proportional to the single-scattering one H U qq→QQ with the prefactor coefficients c qi 0,1,2 given by where we have the color factor C F for both initial-state and final-state double scattering, representing the color interaction strength between the incoming light quark q (or outgoing heavy quark Q) and the soft partons in the nucleus. We also recall that for prompt heavy quark final states the definition of the Mandelstam variables is given in Eq. (6). It is instructive to notice that because the mass generates extra terms in the above pre-factor coefficients, we do not have the simple compact form as in the light flavor fragmentation. As a consistency check, if one takes heavy quark mass m c → 0 limit, we recover the same result for qq → q ′q′ [26].
FIG. 3. The central-cut diagrams for initial-state (left) and final-state (right) double scattering in the gluon-gluon fusion process. The "H"-blobs represent the hard gg → QQ processes as shown in Fig. 1(b).
On the other hand, the double scattering contributions to the gg → QQ fusion process are given in Fig. 3. The calculations are similar, and the final result reads where again the double scattering hard-part functions are proportional to the single scattering one H U gg→QQ with the following pre-factor coefficients: Here, we have the color factor C A (C F ) for initial-state (final-state) double scattering, representing the color interaction strength between the incoming gluon g (outgoing heavy quark Q) and the soft partons in the nucleus. Again, we can not express the final results in the most compact possible form due to the finite mass correction terms. If we take m c → 0, however, we recover the same result for gg → qq [26].
Combining the cross sections in Eq. (7), (15), and (22) with the corresponding hard-part functions, we have the final result for the double scattering contribution to heavy meson production in p+A collisions as which will be used in our phenomenological studies in the next section. It is important to emphasize again that only the central-cut Feynman diagrams contribute in the backward region (outside small-x). In other words, the incoherent double scattering contributions in Eq. (29) has no interference effects and thus are expected to give positive contributions to the heavy meson cross sections in the backward rapidity region in p+A collisions, as we will show in the next section.

III. PHENOMENOLOGY
In this section we present phenomenological applications of our analytic results. We study heavy meson production in p+A collisions at both RHIC and LHC energies.
We first show numerical results for heavy meson production in p+p collisions, in which only single scattering contributions are relevant and we start from the cross section in Eq. (3). We use CTEQ6L1 parton distribution functions [52] and the KKKS parametrization for heavy meson fragmentation functions [53]. We choose the factorization and renormalization scales to be equal throughout the numerical study and set µ ∼ m T = m 2 c + p 2 h⊥ with m c = 1.5 GeV for D meson production. In order to account for higher order QCD contributions, we include a phenomenological K-factor. We find that with a same value of the K-factor, K = 2, our leading order formalism describes reasonably well the data from both RHIC at center-of-mass energy √ s = 200 GeV and LHC at √ s = 2.76, 7 TeV for D-meson productions. In the left and middle panels of Fig. 4 we compare our calculations to LHC data for D * + meson production at √ s = 2.76 TeV [48] and √ s = 7 TeV [49], respectively. We also consider RHIC D * + meson production at √ s = 200 GeV, shown in the right panel. Similar agreements are also found for D 0 and D + data. This provides us with a reasonable p+p baseline for our study of heavy meson production in p+A collisions. To this order, any remaining discrepancy in the overall normalization of the cross section will cancel out in the nuclear modification ratio discussed below.
Keeping in mind that the Cronin-like enhancement can be considerable for both open heavy flavor and quarkonia [9,11], we turn our attention to the study of nuclear effects on heavy meson production in the backward rapidity region in p+A collisions. As we have emphasized in the last section, this is the region where the incoherent multiple scattering are relevant and nuclear enhancement (because of the positive contributions from incoherent double scatterings) 76 TeV with rapidity |y| < 0.5 [48]; Middle: LHC case with √ s = 7 TeV with rapidity |y| < 0.5 [49]; Right: RICH case at √ s = 200 GeV with rapidity |y| < 1 [50,51]. We set µ = mT and use a same K-factor, K = 2 for both RHIC and LHC energies.
should be expected. Our numerical simulations below confirm that this is indeed the case. The nuclear effect is usually quantified by the nuclear modification factor R pA defined as follows: where N coll is the average number of binary collisions. R pA is defined such that the deviation from unity reveals the presence of non-trivial nuclear effects in p+A collisions. In our formalism, the denominator in Eq. (30) represents the heavy meson production cross section in p+p collisions, as given by Eq. (3). On the other hand, the numerator represents the heavy meson cross section in p+A collisions, which receives the contributions from both single and double scatterings, i.e., given by the sum of Eq. (3) and (29). To numerically evaluate the double scattering contributions, the only unknown ingredients in our formalism are the twist-4 quark-gluon and gluon-gluon correlation functions. These functions represent non-perturbative properties of the nuclear medium, and should in principle be extracted from the experimental data. In Refs. [21,45], they have been parametrized as where ξ 2 is a universal quantity, representing a characteristic scale and the strength of parton multiple scattering. At tree level, ξ 2 can be treated as a fixed number. High order corrections may lead to residual energy and factorization scale dependence of ξ 2 after the leading dependence captured by f q,g/A (x) is taken into account, see [27]. ξ 2 = 0.09 − 0.12 GeV 2 was extracted from deep inelastic scattering data in e+A collisions [45], which has been used to describe successfully the nuclear suppression of single inclusive hadron production [20] and the di-hadron transverse momentum imbalance and correlations in d+Au collisions at forward rapidities at RHIC √ s = 200 GeV [21]. For the purpose of numerical study below, we will use the same value ξ 2 = 0.09 − 0.12 GeV 2 .
In Fig. 5 we plot the nuclear modification factors R pA for muons coming form the semi-leptonic open heavy flavor decays in the backward rapidity region in minimum bias d+Au (or p+Pb) collisions as a function of the muon transverse momentum p T . We compare to experimental data from both RHIC [31] (left) and LHC (right) [54] energies. As in Ref. [9], we use the PYTHIA event generator [55] to simulate the full kinematics of Dalitz decays of heavy mesons in both p+p and p+A collisions. RHIC data is from the PHENIX collaboration at √ s = 200 GeV and rapidity −2 < y < −1.4 [31]. LHC data is from the ALICE collaboration at √ s = 5.02 TeV and rapidity −4 < y < −2.96 [54]. As anticipated, the incoherent double scatterings give positive contributions to the heavy meson cross section, which lead to a Cronin-like enhancement in the intermediate p T region. Such an enhancement disappears at large p T , simply because of the high-twist nature of the double scatterings in this framework, i.e. the twist-4 contributions is power suppressed ∼ 1/p 2 T . On the other hand, it increases at the low p T region (even slightly above the RHIC data), because of the same reason, as well as due to the theoretical uncertainty of the p + p baseline in the choice of factorization scale. Our numerical calculations give a reasonable description of both RHIC and LHC data, though slightly below the LHC data. Such a slight difference could be further studied in the future. On the  [31]; Right: LHC at √ s = 5.02 TeV with −4 < y < −2.96 in p+Pb collisions [54]. We set µ = mT , and the band corresponds to ξ 2 = 0.09 − 0.12 GeV 2 .
experimental side there is no baseline measurements for heavy meson production in p+p collisions at the same centerof-mass energy, √ s = 5.02 TeV, and this can introduce uncertainties in the determined nuclear modification factor R pA . On the theory side, one could study the energy/scale dependence of the parameter ξ 2 , as well as other nuclear matter effects. As pointed out in the introduction, other approaches, such as the use of universal nPDFs, are also considered in the literature. In Fig. 6 we present a comparison between our approach and the nPDFs approach (using EPS09 parametrization [30]). As can be seen from Fig. 6 (left), the result from EPS09 cannot describe the experimental data on heavy flavor decay muons at RHIC energy [31]. For LHC energy, the results from high-twist and EPS09 are consistent at large p T region, and disagree in the low p T region. As we have clarified already, the main difference between our approach and the one that uses nPDFs is that the nuclear modification in high-twist formalism comes from the double scattering contributions, which are process-dependent, that is non-universal. However, the latter comes from the universal nPDFs. Therefore, the agreement between our results and experimental data could indicate the non-universality (process dependence) of cold nuclear matter effects.
We would also like to point out that our calculations are restricted to LO, while higher order contributions are accounted for by implementing a phenomenological K-factor. This introduces uncertainty in the theoretical calculation of the differential cross sections. As one can see from the uncertainty band in Fig. 6, the effect of scale variation m T < µ < 2m T in the nuclear modification ratio is greatly reduced and smaller than the uncertainty due to the choice of ξ 2 .

IV. SUMMARY AND DISCUSSION
In this paper, we studied the effect of multiple scattering on heavy flavor meson production in p+A collisions. We concentrated on the backward rapidity regime, where the parton momentum fraction x in the nucleus is relatively large and the multiple scattering between the probe and the nuclear medium is incoherent. Within the high-twist factorization formalism, we evaluated the initial-state and final-state interactions relevant to heavy meson production in p+A collisions. We found that our final results depend on both the twist-4 quark-gluon and gluon-gluon correlation functions. Using the existing parametrization for these correlation functions, we calculated numerically the double scattering contribution to the differential cross section of heavy mesons in the kinematic region relevant to both RHIC and LHC experiments. Our simulations describe quite well the nuclear modification factor for muons coming from the semi-leptonic decays of heavy flavor mesons in the backward rapidity region in p+A collisions, where a Cronin-like enhancement is observed experimentally. This feature is understood as the incoherent multiple scattering of hard partons in the large nucleus. We conclude that the backward rapidity measurement provide a unique opportunity to investigate the perturbative QCD dynamics in a region that has so far not received adequate attention and to help further constrain the properties of cold nuclear matter.
It is important to emphasize that the incoherent multiple scatterings for heavy meson production in p+A collisions involve both initial-state and final-state interactions. To disentangle these two multiple scattering effects, it is instructive to consider processes involving a photon either in the initial state or final state of the process, as shown in Ref. [26]. With the advent of a future electron-ion collider, it will be interesting to study heavy meson production in e+A collisions. This process can provide us with a clean channel to investigate purely the final-state multiple scattering effect. Such a study could also help clarify the non-universality of nuclear effects in different processes due to the process dependent hard part coefficient [44], and test the predictive power of the high twist formalism due to the universality of the twist-4 parton-parton correlation functions. We leave this study for future work.