Heavy Flavor Azimuthal Correlations in Cold Nuclear Matter

Background: It has been proposed that the azimuthal distributions of heavy flavor quark-antiquark pairs may be modified in the medium of a heavy-ion collision. Purpose: This work tests this proposition through next-to-leading order (NLO) calculations of the azimuthal distribution, $d\sigma/d\phi$, including transverse momentum broadening, employing $$ and fragmentation in exclusive $Q \bar Q$ pair production. While these studies were done for $p+p$, $p + \bar p$ and $p+$Pb collisions, understanding azimuthal angle correlations between heavy quarks in these smaller, colder systems is important for their interpretation in heavy-ion collisions. Methods: First, single inclusive $p_T$ distributions calculated with the exclusive HVQMNR code are compared to those calculated in the fixed-order next-to-leading logarithm approach. Next the azimuthal distributions are calculated and sensitivities to $$, $p_T$ cut, and rapidity are studied at $\sqrt{s} = 7$ TeV. Finally, calculations are compared to $Q \bar Q$ data in elementary $p+p$ and $p + \bar p$ collisions at $\sqrt{s} = 7$ TeV and 1.96 TeV as well as to the nuclear modification factor $R_{p {\rm Pb}}(p_T)$ in $p+$Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV measured by ALICE. Results: The low $p_T$ ($p_T<10$ GeV) azimuthal distributions are very sensitive to the $k_T$ broadening and rather insensitive to the fragmentation function. The NLO contributions can result in an enhancement at $\phi \sim 0$ absent any other effects. Agreement with the data was found to be good. Conclusions: The NLO calculations, assuming collinear factorization and introducing $k_T$ broadening, result in significant modifications of the azimuthal distribution at low $p_T$ which must be taken into account in calculations of these distributions in heavy-ion collisions.


I. INTRODUCTION
Recently there has been interest in heavy flavor correlations and how they might be modified in heavy-ion collisions [1][2][3]. However, before drawing any conclusions, it is worth studying these correlations in p + p collisions and how they may be affected in cold nuclear matter. This paper focuses on heavy quark pair correlations in azimuthal angle.
Heavy flavor correlations were previously studied at high p T and high pair invariant mass where they contribute to the background for Z 0 boson decays [4]. They were also studied in more elementary collisions, including the effects of final-state radiation in event generators [5] and in the k T -factorization approach [6].
Correlated production, specifically of the azimuthal angle between heavy flavors, φ, either by direct reconstruction of both D mesons or of a D meson and the decay product of its partner, either a light hadron or a lepton [1], is a stronger test of QQ production than single inclusive distributions. Naively, at LO QQ pairs are produced back-to-back with a peak at φ = π. Higher order production will, however, result in a more isotropic distribution in φ due to light parton emission in the final state.
Single inclusive heavy flavor production is discussed and compared with data in Sec. II. Because it is not possible to study pair correlations with current approaches to single inclusive distributions such as FONLL [7] and GV-VFNS [8], the exclusive HVQMNR NLO code [9] is employed to calculate the azimuthal correlations. By modifying the fragmentation function and the intrinsic transverse momentum broadening in the HVQMNR code, it is possible to reproduce the shape of the sin-gle heavy flavor transverse momentum distributions from codes like FONLL. This gives confidence in the approach to the calculation of the pair distributions. In this section, methods of calculating exclusive QQ pair production are briefly introduced and a comparison of the nextto-leading order calculation with QQ production in leading event generators is discussed.
After demonstrating that the calculational approaches give reasonably equivalent results, the sensitivity of the azimuthal distributions to the heavy flavor fragmentation function and the intrinsic transverse momentum, k T , is explored. Sensitivities to the heavy quark transverse momentum cut, the rapidity range probed, and the renormalization and factorization scales are also discussed in Sec. III. Section IV describes the sensitivity of the azimuthal distributions to the size of k 2 T , independent of the fragmentation function.
The results are compared to QQ data from p + p and p + p collisions in Sec. V and to ALICE R pPb (p T ) single D meson data from p+Pb collisions at √ s N N = 5.02 TeV in Sec. VI. In Sec. VI, predictions for cold matter effects on the φ distribution are also shown.

II. HEAVY FLAVOR PRODUCTION
A. Single inclusive production approaches There are currently two approaches to heavy flavor production at colliders: collinear factorization and the k T -factorization approach, usually employed at low x.
There are two main methods of calculating the spectrum of single inclusive open heavy flavor production in perturbative QCD assuming collinear factorization. The underlying idea is similar but the technical approach differs. Both seek to cure the large logarithms of p T /m arising at all orders of the perturbative expansion which can spoil the convergence. The first terms in the expansion are the leading (LL) and next-to-leading logarithmic (NLL) terms, α 2 s [α s log(p T /m)] k and α 3 s [α s log(p T /m)] k respectively.
It is worth noting that the single inclusive heavy flavor p T distribution is finite at leading order (LO) as p T → 0 because of the finite quark mass scale. This is in contrast to light hadron and jet production where there is no mass scale to regulate the low p T cross section.
While large uncertainties can arise at low p T , these are due to the choice of factorization scale and are larger for charm than for bottom quark production [10,11]. A next-to-leading order (NLO) calculation that assumes production of massive quarks but neglects LL terms, a "massive" formalism, can result in large uncertainties at high p T [12][13][14][15]. (This massive formalism is sometimes referred to as a fixed-flavor-number (FFN) scheme.) If, instead, the heavy quark is treated as "massless" and the LL and NLL corrections are absorbed into the fragmentation functions, the approach breaks down as p T approaches m even though it improves the result at high p T . The massless formalism is sometimes referred to as the zero-mass, variable-flavor-number (ZM-VFN) scheme [16,17]. There are methods to interpolate smoothly between the FO/FFN scheme at low p T and the massless/ZM-VFN scheme at high p T .
The fixed-order next-to-leading logarithm (FONLL) approach is one such method. In FONLL, the fixed order and fragmentation function approaches are merged so that the mass effects are included in an exact calculation of the leading (α 2 s ) and next-to-leading (α 3 s ) order cross section while also including the LL and NLL terms [7]. The NLO fixed order (FO) result is combined with a calculation of the resumed (RS) cross section in the massless limit. The FO and RS approaches need to be calculated in the same renormalization scheme with the same number of light flavors. The FONLL result is then, schematically, where FOM0 is the fixed order result at zero mass. The interpolating function G(m, p T ) ∼ p 2 T /(p 2 T + (cm) 2 ) is arbitrary but must approach unity for m/p T → 0. Note that the number of light flavors is thus 4 for charm and 5 for bottom quark production in FONLL, in contrast the 3 for charm and 4 for bottom in the fixed-order calculation.
The second interpolation scheme is the generalizedmass variable-flavor-number (GM-VFN) scheme [8]. The large logarithms in charm production for p T ≫ m are absorbed in the charm parton distribution function and are thus included in the evolution equations for the parton distributions. The logarithmic terms can be incorporated into the hard cross section to achieve better accuracy for p T ≥ m. By adjusting the mass-dependent subtraction terms, no interpolating function is required [8].
In the k T -factorization approach, off-shell leading order matrix elements for g * g * → cc are used together with unintegrated gluon densities that depend on the transverse momentum of the gluon, k T , as well as the usual dependence on x and µ F . The motivation for choosing the k T -factorized approach is that, at sufficiently low x, collinear factorization should no longer hold.
The LHC data has been compared to calculations in both approaches and those assuming collinear factorization compare well with the LHC data. Recent ALICE data [18], at 0 < p T < 2 GeV supports collinear factorization. Their D 0 results at 7 TeV for |y| < 0.5, were compared to FONLL, GV-VFNS and LO k T -factorization calculations. Only the k T -factorized calculation is inconsistent with the shape of the data in the p T range where the calculation should best apply. The forward rapidity data of LHCb at 7 TeV [19] and 13 TeV [20], also agree well with the collinear factorization assumption. The recent 13 TeV data from LHCb is within the uncertainty bands of FONLL and POWHEG (discussed below) down to p T → 0, albeit with large uncertainty bands and near the upper limit of the calculated band. The GM-VFNS calculation is given only for p T > 3 GeV but agrees well with the data with small uncertainties [20]. Perhaps a NLO k T -factorized result could lead to improved agreement with the data but, so far, collinear factorization appears to still work well for low x, moderate p T charm production.
This paper will thus employ the collinear factorization approach in the calculations.

B. Exclusive approaches to heavy flavor production
Recall that the FONLL and GM-VFNS calculations are for single inclusive production only and can thus not address QQ pair observables. There are NLO heavy flavor codes that, in addition to inclusive heavy flavor production, also calculate exclusive QQ pair production.
The HVQMNR code [9] uses negative weight events to cancel divergences numerically. Smearing the parton momentum through the introduction of intrinsic transverse momenta, k T , reduces the importance of the negative weight events at low p T . HVQMNR does not include any resummation. POWHEG-hvq [21] is a positive weight generator that includes leading-log resummation. The entire event is available since PYTHIA [22] and HERWIG [23] are employed to produce the complete event after production of the QQ pair.
The k T -factorization approach can also be employed to calculate correlated cc production since the unintegrated parton densities have a transverse component, giving finite p T and φ distributions even at LO [6].
The HVQMNR code is employed here to focus on the effects due to the NLO contribution alone, including k Tbroadening and fragmentation but excluding the parton showers which can further randomize the pair momenta 3 and thus affect the angular correlations.
C. Heavy flavor production in leading order event generators compared to NLO calculations In addition to these NLO codes, heavy flavor correlations can also be simulated employing LO event generators such as PYTHIA.
In an event generator like PYTHIA or HERWIG, heavy flavor production is divided into three different categories: flavor or pair creation, flavor excitation, and gluon splitting [24]. This classification depends on the number of heavy quarks in the final state of the hard process, defined as the process in the event with the highest virtuality. Pair creation is equivalent to the four leading order diagrams, shown in Fig. 1. In this case the final state has two heavy quarks, the Q and Q. In the case of flavor excitation, a heavy quark from the splitting g → QQ in an initial-state parton shower is put on mass shell by scattering with a parton, either a quark or gluon, from the other beam, qQ → qQ or gQ → gQ. There is thus one heavy quark in the final state of the hard scattering. Finally, gluon splitting is defined as having no heavy flavor in the hard scattering such as gg → gg. The QQ pair is produced in an initial-or final-state parton shower via g → QQ. Double counting of these processes is avoided by requiring that the hard scattering should be of greater virtuality than the parton shower [24].
There are parameters that can be tuned, depending on the generator employed, that can match the distributions from the LO generator to those of a NLO calculation. However, that does not mean that the mix of physics processes is identical. For example, all three of the QQ production processes in PYTHIA contribute to perturbative QCD production at next-to-leading order with 2 → 3 processes. Pair creation is realized at NLO by gluon emission from one of the final-state heavy quarks, as in Fig. 2(a), the NLO version of Fig. 1(c). Flavor excitation off a gluon from the opposite hadron is shown in Fig. 2(b) while gluon splitting is shown in Fig. 2(c).
The diagrams in Fig. 2(a)-(c) are an illustration of the possible NLO diagrams included in gg-initiated QQ production, gg → QQg. For example, virtual corrections and crossed diagrams are not shown in Fig. 2. In a NLO calculation, the amplitudes from these gg contributions are summed with weights determined by the spin and color factors of each diagram, not parton showers and individual virtualities. Instead, production processes are separated by the initial state. The gg initial state is dominant at collider energies over most of the rapidity range. The next largest contribution is due to qg or qg scattering, qg → qQQ, shown in Fig. 2(d), and not present at leading order. This diagram is considered to be flavor excitation in event generators. Finally, real and virtual corrections in the qq channel, qq → QQg, makes a small contribution to QQ production at high energies.
Due to the common use of event generators in simu- lations and analysis, there is a great temptation to try and analyze heavy flavor production in terms of separate production mechanisms, as characterized by creation, excitation and splitting. However, as described here, there is no distinction, all three are components of NLO production. Separating the diagrams in the initial gg process in such a manner will result in improper interferences in the calculation and an incorrect cross section. Thus it is preferable to do these types of analyses with an exclusive NLO code such as described in Sec. II B.
Examples of real contributions to next-to-leading order QQ production. Diagrams (a)-(c) illustrate contributions to gg → QQg while (d) shows an example of qg → qQQ production.

D. kT broadening and fragmentation
The transition from bare quark distributions to those of the final-state hadrons is accomplished by including a fragmentation function and intrinsic transverse momentum, k T , broadening. The implementation of these two effects are described here.

Intrinsic kT broadening
Results on open heavy flavors at fixed-target energies indicated that some level of transverse momentum broadening was needed to obtain agreement with the low p T data after fragmentation was applied [25]. Broadening is typically applied by including some intrinsic transverse momentum, k T , smearing to the initial-state parton densities. The implementation of intrinsic k T in HVQMNR is not handled in the same way as calculations of other hard processes due to the nature of the code. In HVQMNR, the cancellation of divergences is done numerically. Since adding further numerical Monte-Carlo integrations would slow the simulation of events, as well as require multiple runs in the same kinematics but with different intrinsic k T kicks, the kick is added in the final, rather than the initial, state. The Gaussian function g p (k T ) [25], multiplies the parton distribution functions for both hadrons, assuming the x and k T dependencies in the initial partons completely factorize. If factorization applies, it does not matter whether the k T dependence appears in the initial or final state as long as the kick is not too large. In Ref. [25], k 2 T = 1 GeV 2 was chosen to describe the p T dependence of fixed-target charm production.
In HVQMNR, the QQ system is boosted to the rest frame from its longitudinal center-of-mass frame. Intrinsic transverse momenta of the incoming partons, k T 1 and k T 2 , are chosen at random with k 2 T 1 and k 2 T 2 distributed according to Eq. (2). A second transverse boost out of the pair rest frame changes the initial transverse momentum of the QQ pair, p T , to p T + k T 1 + k T 2 . The initial k T of the partons could have alternatively been given to the entire final-state system, as is essentially done if applied in the initial state, instead of to the QQ pair. There is no difference if the calculation is LO but at NLO an additional light parton can also appear in the final state, making the correspondence inexact. In Ref. [25], the difference between the two implementations is claimed to be small if k 2 T ≤ 2 GeV 2 . The p T -integrated rapidity distribution is unaffected by the intrinsic k T .
The level of the intrinsic k T kick required in these calculations was determined by comparison with the shape of the quarkonium p T distributions. In this case, the Color Evaporation Model is employed to calculate production with a cut on the pair invariant mass in the HVQMNR code, giving an upper limit of 2m H where H = D, B for charm and bottom respectively. The CEM pair distributions require augmentation by k T broadening to make them finite at p T → 0, as do the QQ pair p T distributions. The k T broadening has the most important effect on the azimuthal correlation at low p T as well, as we will show.
The effect of the k T kick on the p T distribution can be expected to decrease as √ s increases because the average p T also increases with energy. However, the value of k 2 T is assumed to increase with √ s so that effect remains important for low p T production at higher energies. The energy dependence of k 2 T in Ref. [11] is Comparison with the RHIC J/ψ data found that n = 12 gave the best description of the J/ψ p T distribution both at central and forward rapidity [11]. The value of n was unchanged in the Improved Color Evaporation Model [26]. A smaller value of n and thus a larger k 2 T is required for the Υ p T distribution. For Υ, n = 3 is set by comparison to the Tevatron results at √ s = 1.8 TeV [27]. These same values are also used in the charm and bottom pair distributions.

Fragmentation
The default fragmentation function in HVQMNR is the Peterson function [28], where z represents the fraction of the parent heavy flavor quark momentum carried by the resulting heavy flavor hadron. In the original Peterson function, the nominal values of the fragmentation parameter ǫ P were 0.06 for charm and 0.006 for bottom. These values result in z = 0.671 for charm and 0.828 for bottom, a reduction in the heavy quark momentum of ∼ 33% and 17% respectively, as shown by the red curves in Fig. 3. Currently, the Peterson function with ǫ P = 0.06 is considered too strong for charm production. The FONLL fragmentation scheme for open heavy flavor is softer [29], as will be discussed.
The value of ǫ P in the Peterson fragmentation function must be modified to match the FONLL p T distribution for single D and B meson production since the standard ǫ P values for Eq. (4) result in much softer p T distributions than those of FONLL. Therefore, the value of ǫ P needs to be reduced to match the average ratio of the heavy quark momentum transferred to the final state hadron, z . The average z is calculated as Fragmentation functions for D and D * mesons have been calculated in an approach consistent with an FONLL calculation [30]. The FONLL charm fragmentation function was determined from Mellin moments of the distributions calculated consistently in the same framework. The D 0 fragmentation function includes both pseudoscalar and vector parts that account for the ground state, c → D 0 , and excited state contributions, c → D * 0 → D 0 and c → D * + → D 0 respectively. The fragmentation parameter for D 0 with a central charm quark mass of m = 1.5 GeV [29] was adjusted to the chosen central value of the charm quark mass, 1.27 GeV, used here. Assuming a linear dependence of ǫ P on charm mass from 1.2 to 1.7 GeV, the total FONLL charm fragmentation function is given by the blue dashed curve in Fig. 3(a). The vector and pseudoscalar contributions are shown in the cyan and magenta curves respectively. Individually, their average z values are z = 0.812 for the pseudoscalar channel and 0.843 for the more massive D * in the vector channel. Note that the vector channel is a small overall contribution to the combined fragmentation function, as is reflected in the average z for the total FONLL fragmentation function, 0.822. This z can be compared to the result obtained with the default value of ǫ P in the Peterson function, z = 0.671. The calculations in this work will use the reduced value, ǫ P = 0.008, shown in the black curve in Fig. 3(a). Using this ǫ P , z = 0.820, in good agreement with the average z for the combined FONLL charm fragmentation function. As will be shown in Sec. II E, the single inclusive D meson p T distribution calculated with HVQMNR for this ǫ P , combined with the value of k 2 T determined from RHIC with the energy dependence of Eq. (3), is in good agreement with the FONLL p T distribution.
The bottom quark fragmentation function in FONLL is of the form where ǫ b = 34 for a bottom quark mass of 4.75 GeV [29]. Again, assuming a linear dependence of ǫ b on b quark mass from 4.5 to 5 GeV, with m b = 4.65 GeV, ǫ b = 27.5. The FONLL fragmentation function with this value of ǫ b , is shown in the dashed blue curve of Fig. 3(b). In this case, z = 0.934. In contrast, the red curve shows the default Peterson function result for ǫ P = 0.006, with z = 0.828. If the same overall reduction of the Peterson function parameter for charm is used for bottom, from 0.006 to 0.0008, z = 0.911, insufficient to replicate the FONLL B meson p T dependence. Thus, ǫ P = 0.0004 is used here, giving z = 0.930, in good agreement with the average z from FONLL, as shown in the solid black curve of Fig. 3(b). The charm and bottom p T distributions employing these fragmentation functions will be compared to those from FONLL in the next section.

E. Single Inclusive Heavy Flavor Distributions
In the calculations reported in this paper, the same values of the charm quark mass and scale parameters as in Ref. [11] are employed here, (m, µ F /m T , µ R /m T ) = (1.27 ± 0.09 GeV, 2.1 +2.55 −0.85 , 1.6 +0.11 −0.12 ) where µ F is the factorization scale and µ R is the renormalization scale. In the case of bottom production, (m, µ F /m T , µ R /m T ) = (4.65 ± 0.09 GeV, 1.4 +0.77 −0.49 , 1.1 +0.22 −0.20 ) is used [27]. The  CT10 proton parton densities [31] are employed in the calculations. The scale factors, µ F and µ R , are defined relative to the transverse mass of the pair, ). At LO in the total cross section, the QQ pair p T is zero. Thus, while the calculation is NLO in the total cross section, it is LO in the pair distributions. In the exclusive NLO calculation [9] both the Q and Q variables are retained to obtain the pair distributions. Unless otherwise noted, the calculations employ the central values of the heavy quark mass and scale factors. First, the relative importance of k T broadening and fragmentation on the single inclusive heavy quark p T distribution is studied. Undertaking a study of the azimuthal correlations between heavy quarks requires finding the values of k 2 T and ǫ P appropriate for a reasonably faithful reproduction of the FONLL single inclusive heavy quark distributions in the same kinematics employing the HVQMNR code. These same values of k 2 T and ǫ P are then used in the calculation of the azimuthal distributions. The values of k 2 T obtained for J/ψ and Υ production are assumed as a default. Figure 4 shows how k 2 T and ǫ P affect the low p T part of the spectrum at √ s = 7 TeV. As examples, the charm distributions are shown at forward rapidity, 2.5 < y < 5, the LHCb acceptance [19], while the bottom distributions are shown at central rapidity, |y| < 2.4, the rapidity region of b hadron distributions reported by ATLAS [32]. Note that the lower p T range is emphasized in Fig. 4 because this region is most affected by k T broadening, as discussed in more detail later.
The blue solid curves in Fig. 4 are the bare quark distributions, with k 2 T = 0, ǫ P = 0. This distribution is the hardest of all the cases shown and thus has the highest average p T . (See Table I for the values of p T with each k 2 T , ǫ P combination.) The dashed curves, with the lowest average p T , include no broadening but employ the default Peterson fragmentation function parameter from e + e − data [33], ǫ P = 0.06 for charm and 0.006 for bottom, corresponding to the red curves in Fig. 3. The effect on the charm quark distributions is particularly strong because this value of ǫ P results in a ∼ 33% decrease in p T relative to the bare quark, as already discussed. At fixed-target energies, such a strong reduction in momentum could only be made compatible with data, which agreed rather well with the low p T bare charm quark distribution, by setting k 2 T = 1 GeV 2 [25]. For fixed-target energies, √ s ∼ 20 GeV, this was sufficient because the average p 2 T of the charm quark was low since m c / √ s ∼ 0.065. Therefore, k 2 T ∼ p 2 T . However, at collider energies, m c / √ s is small, ∼ 1.8 × 10 −4 at √ s = 7 TeV. Thus an average k 2 T on the order of 1 − 3 GeV 2 has a rather small effect on the shape of the p T distribution, particularly at high p T . On the other hand, ǫ P reduces the quark p T uniformly, independent of the center-of-mass energy.
The effect of fragmentation on the bottom quark distribution is less striking because ǫ P is an order of magnitude smaller. Nonetheless, as will be shown in Fig. 5, it is clear that even this milder effect may be too strong to be compatible with the b-hadron data.
While k 2 T ∼ 1 GeV 2 was sufficient to mitigate the effect of fragmentation on charm production at √ s = 20 GeV [25], it is clear that any value of k 2 T large enough to do so at √ s = 7 TeV would be too large to be physical. Thus, a reduction of ǫ P by a factor of ∼ 7.5, to 0.008 for charm and 0.0008 for bottom, was checked (dotdashed curves in Fig. 4). This reduction is adequate for charm production, giving a result intermediate to the free quark distribution and that with the standard Peterson fragmentation parameter ǫ P . However, it was found necessary to reduce ǫ P by an additional factor of two, to 0.0004, for bottom (red solid curve in Fig. 4(b)). This additional reduction is sufficient to produce a b-hadron distribution similar to that of FONLL when broadening is included.
Next, the effect of k 2 T broadening is introduced in addition to fragmentation. As mentioned previously, the default values of k 2 T for charm and bottom are assumed to be those found to agree with low p T quarkonium production in the Improved Color Evaporation Model [26]. The energy dependence of k 2 T is given in Eq. (3). At √ s = 7 TeV, Eq. (3) results in k 2 T ∼ 1.5 GeV 2 for charm and ∼ 3 GeV 2 for bottom. It is clear that even these relatively large values, although still less than p 2 T at √ s = 7 TeV, have a rather small effect on the single inclusive heavy quark p T distribution. Indeed, doubling k 2 T for charm (magenta curve in Fig. 4(a)) results in only a small change in the p T distribution in the range shown. No further increase of k 2 T is shown for bottom quarks because such values, k 2 T > 3 GeV 2 , would be too large for the assumption of the equivalence of adding k 2 T in the initial or final state [25]. In Fig. 5, the final values of k 2 T and ǫ P , ∼ 1.5 GeV 2 , 0.008 for charm and 3 GeV 2 , 0.0004 for bottom, are used to calculate the uncertainty bands on the HVQMNR results and compare them to those of FONLL, with its default fragmentation functions. The same mass and scale parameters are employed in the HVQMNR and FONLL calculations. The results are compared to LHCb D 0 meson data in Fig. 5(a), and ATLAS b-hadron data in  The uncertainty band can be obtained for the best fit sets [11,27] by adding the uncertainties from the mass and scale variations in quadrature. The envelope contained by the resulting curves, defines the uncertainty on the cross section. The kinematic observable X can be e.g. p T or azimuthal angular separation φ. In the calculation labeled "cent", the central values of m, µ F and µ R are used while in the calculations with subscript µ, the mass is fixed to the central value while the scales are varied and in the calculations with subscript m, the mass is varied while the scales are held fixed.
As an example, Fig. 5(a) compares the calculations in the rapidity interval 2.5 < y < 3 to the LHCb D 0 data [19]. The LHCb data in the range 2.5 < y < 5 was divided into five rapidity bins, each of 0.5 units of rapidity in width. The lowest of these bins is shown as an example. The agreement between the two calculations is generally very good. Some disagreement becomes apparent for p T > 10 GeV where the upper limit of the band calculated using the exclusive HVQMNR code becomes slightly harder than the FONLL result, as can be expected because FONLL was designed to cure large logarithms of p T /m and here p T /m > 10. This region, as will be shown, is also relative insensitive to effects of k 2 T on the azimuthal distributions.
There is a somewhat larger difference in the calculations of the b-hadron distributions in Fig. 5(b). In this case, the results are shown, along with the ATLAS data [32], at midrapidity and up to considerably higher p T . It is, however, worth noting that, at p T ∼ 50 GeV, p T /m for bottom is similar to that for charm at p T ∼ 10 GeV. Here the one standard deviation fit to the bb cross section gives a narrower range on m b and µ F /m T , µ R /m T . To improve the agreement between the calculations still further it would be necessary to employ ǫ P = 0 in the HVQMNR calculation since the FONLL b-hadron p T distribution is equivalent to the HVQMNR bare b-quark distribution.

III. QQ PAIR AZIMUTHAL DISTRIBUTIONS
The single inclusive heavy quark distributions have been calculated both in the exclusive QQ code, HVQMNR, and the inclusive-only FONLL approach and have been shown to agree well with each other. The next step is to calculate the azimuthal angular separation between the two heavy quarks. In Fig. 6, the distributions, dσ/dφ, are calculated to next-to-leading order for the same choices of k 2 T and ǫ P , with the central mass and scale values, as in Fig. 4, using HVQMNR only since now an exclusive calculation is required to obtain dσ/dφ. First, some topological details relevant to the azimuthal distributions of QQ pair production are discussed.
The leading order perturbative QCD 2 → 2 diagrams for heavy flavor production, shown in Fig. 1, produce back-to-back QQ pairs with an azimuthal separation of φ = π. At leading order, the p T of the pair of heavy quarks is zero and dσ/dφ is represented by a delta function, δ(φ − π). Note also that, already at leading order, the single inclusive heavy quark p T distribution is maximally hard so that higher order corrections do not significantly change the shape of the p T distribution [34,35].
As discussed in Sec. II C, at next-to-leading order, both virtual and real conntributions arise. The virtual corrections are typically the exchange of soft gluons at the vertices while real corrections give rise to 2 → 3 processes, see Fig. 2. These corrections smear out the azimuthal separation in the leading order contribution so that while there is still a peak at φ = π, the pairs are no longer strictly back-to-back. Instead, there is a tail toward φ ∼ 0, even with k 2 T = 0 and ǫ P = 0. The NLO contribution can also lead to different azimuthal topologies, as in Fig. 2(b)-(d), due to 'flavor excitation' and 'gluon splitting'. These diagrams do not harden the single inclusive heavy flavor p T distribution. However, they do tend to make the azimuthal correlation more isotropic [24]. The 2 → 3 processes produce a non-singular azimuthal distribution at φ = π before any broadening or fragmentation is taken into account. A final-state light parton emitted from one of the outgoing heavy quarks is likely to be collinear to the heavy quark that emitted it. Thus the heavy quark and anti-quark remain close to a back-to-back configuration with a peak at φ ∼ π. However, if the light parton is from the initial state, as in qg scattering, this parton can be hard and its momentum can be balanced against that of the QQ pair so that φ ∼ 0. In addition, the 'gluon splitting' diagram of Fig. 2(c), will result in a pair with φ ∼ 0 against a hard parton. Thus, the low p T heavy quarks azimuthal distribution is likely to be more isotropic, while high p T heavy quarks will more likely result in a doubly-peaked φ distribution, with peaks at φ ∼ 0 and π.
Fragmentation results in a hadron with momentum nearly collinear to the original heavy quark. Thus fragmentation does not change the general direction of the parent heavy parton. The resulting heavy hadron will, however, be somewhat out of alignment with the initial heavy quark. The larger the value of ǫ P , then, the greater the energy lost to hadronization and the more likely that the Q and Q are out of alignment, resulting in a slight enhancement at φ ∼ 0 relative to ǫ P = 0. Lower values of ǫ P reduce any enhancement at φ ∼ 0, as is evident in the dashed and dot-dashed curves in Fig. 6(a) for ǫ P = 0.06 and 0.008 respectively.
On the other hand, including the randomization of the final-state parton momentum by k T broadening has a substantial efect on dσ/dφ. Indeed, the level of k T broadening required to describe J/ψ production at low p T completely changes the shape of dσ/dφ for cc, resulting in a peak at φ ∼ 0 and none at φ ∼ π. The larger value of k 2 T , ∼ 3 GeV 2 , increases the enhancement at φ ∼ 0, as shown in Fig. 6(a).
In the case of bb pairs, introducing k T broadening still gives a peak near φ ∼ π, even though the k 2 T for bottom shown in Fig. 6(b) is as large as the upper value of k 2 T for charm in Fig. 6(a). This is because k 2 T ∼ 1.5 GeV 2 is on the order of m 2 c , ∼ 1.7 GeV 2 , while k 2 T = 3 GeV 2 for botttom is considerably less than m 2 b ∼ 21.6 GeV 2 . Thus k 2 T /m 2 c ∼ 1 results in a much more isotropic φ distribution while k 2 T /m 2 b ∼ 0.14 results in a lower momentum kick imparted to a bottom quark.
To illustrate the results more clearly, the ratio of all calculations are shown relative to the results for bare quarks, with k 2 T = 0, ǫ P = 0 in Fig. 7. In both cases, the results with only finite ǫ P and k 2 T = 0 are relatively independent of φ with the largest values of (dσ( k 2 T , ǫ P )/dφ)/(dσ( k 2 T = 0, ǫ P = 0)/dφ) for the default values of ǫ P , 0.06 for charm and 0.006 for bottom. For charm production, the ratio is ∼ 1.05 − 1.15 at φ ∼ 0, increasing to ∼ 1.15 − 1.25 near φ ∼ π. For bottom, with its smaller ǫ P , the ratio is only a few percent above unity. Once finite k 2 T is included, there is an almost linear decrease with increasing φ for charm while, for bottom, there is an enhancement until φ ∼ π/2 with a steep falloff as φ → π. We will examine the sensitivity of the φ distributions to k 2 T in the next section. To determine the influence of the heavy quark mass and scale values on dσ/dφ for finite k 2 T , the full uncertainty bands are shown in Fig. 8 for the same values of k 2 T and ǫ P employed in Fig. 5. In this case, in addition to the central value and the upper and lower limits of the band, each mass and scale combination is also shown in the light dot-dashed curves. The results are given as histograms in both cases even though dσ/dφ is relatively smooth, especially for bb pairs. The differences in the shapes of the distributions are larger for charm. Thus the level of the enhancement at φ ∼ 0 seen for the central value varies considerably in this case.
The upper and lower limits of the band are obtained according to Eqs. (7) and (8) with X = φ. On a bin-bybin basis, the maximum and minimum mass and scale combination may vary. For example, changing the factorization scale changes the slope of the p T distribution so that, at some value of p T , the maximum of dσ/dp T may come from a different scale choice. The mass and scale choices affect the shape of dσ/dφ most strongly for charm. As φ → π, the maximum of dσ/dφ approaches the central value while the minimum of the band is essentially determined by a single mass and scale combination.
To end this section, the effect of a p T cut on dσ/dφ for the values of k 2 T and ǫ P used in the calculations of the uncertainty band in Figs. 5 and 8 is shown in Fig. 9. Results are given starting from low p T , p T < 10 GeV, and increasing with p T cuts above 10, 20, 30, 40, 50 and 75 GeV for charm with an additional cut of p T > 100 GeV for bottom. (Note that a minimum p T cut is chosen instead of a p T range because the azimuthal distribution is dominated by the behavior at the lowest p T .) The low p T results, with p T < 10 GeV, are equivalent to those integrated over all p T , shown in Fig. 6. Notably, already at p T > 10 GeV, the shape of dσ/dφ changes with peaks at both φ = 0 and φ ∼ π and a deepening dip developing between the two peaks. This is because that, as the minimum p T increases, the light parton in the 2 → 3 process is more likely to either be low momentum and aligned with one of the two heavy quarks and opposite the other (emission from a final-state heavy quark, Fig. 2(a), φ ∼ π) or high momentum and balancing its momentum against that of the QQ pair (φ ∼ 0, Fig. 2(b)-(d)).
Note that no scale factors are applied to separate the distributions, the decrease of dσ/dp T is sufficient to separate them without introducing an additional scale factor. However, in the case of charm production, while no scale factor is applied to the p T < 10 GeV distribution, the remaining distributions are multiplied by a factor of 10 3 to be visible on the plot due to the steeply-falling p T distribution of the charm quarks.
The cc distributions are shown as histograms with per bin uncertainties because the charm quark p T distributions decrease with p T faster than those for bottom quarks so that, by p T > 50 GeV the statistics for cc pairs are poor. There are few cc events produced for p T > 75 GeV, as seen in Fig. 9(a), and these are mostly at φ ∼ 0. Note that the cc results are shown for the forward rapidity region, 2.5 < y < 5, in the LHCb acceptance. In this rapidity range, the p T distributions are more steeply falling than at central rapidity because the edge of phase space for production is reached at lower p T at forward rapidity. However, as will be explained, the shape of dσ/dφ is not strongly dependent on the rapidity range.
The bb distributions in Fig. 9(b) show the same trends but the harder p T distributions for bottom quarks result in smoother azimuthal distributions for bb pairs even for p T > 100 GeV. The higher the p T , the more pronounced the peak at φ ∼ 0 becomes while the enhancement at φ ∼ π does not disappear. (Note that there is no enhancement at φ ∼ 0 for the lowest p T cut, p T < 10 GeV, as is the case for charm.) While the calculations shown in Fig. 9 are for the central values of quark mass and scale factors, the trends would remain the same for all the mass and scale combi-nations, especially for bb, see Fig. 8.
In this section, the sensitivity of dσ/dφ to the size of the intrinsic k 2 T is explored to determine whether or not there is any obvious onset of the change in dσ/dφ. Thus Eq. (3) is modified to introduce a parameter ∆ that reduces the average k 2 T , starting with k 2 T = 0. (Note that ∆ = 1 in Eq. (3).) For charm pair production with n = 12, ∆ = −3/2, −1, −1/2, 0, 1/2, and 1, effectively changing k 2 T by ∼ 0.25 GeV 2 at √ s = 7 TeV as ∆ is increased by 1/2. For bottom pair production, only ∆ = −1/2, 0, 1/2 and 1 result in k 2 T > 0 because, with n = 3, each change in ∆ increases k 2 T by ∼ 1 GeV 2 at √ s = 7 TeV. Here ∆ = −1/2 results in k 2 T = 0.02 GeV 2 . The broad p T cuts, p T < 10 GeV and p T > 10 GeV, are studied both at midrapidity (|y| < 2.4) and forward rapidity (2.5 < y < 5). Only the results for midrapidity are shown here, however, since trends were found to be independent of rapidity. The figures in this section include the cc and bb azimuthal distributions, dσ/dφ, as well as the ratios, (dσ( k 2 T )/dφ)/(dσ( k 2 T = 0)/dφ), to highlight the effects when differences in the distributions themselves may be small. Note that the choice of 10 GeV is somewhat arbitrary and not indicative of any threshold behavior. Comparison with data in the next section will be more indicative of the importance of k 2 T in a given p T range. Figure 10(a) shows the cc azimuthal distributions for p T < 10 GeV at central rapidity. Note that increasing k 2 T from 0 to 0.25 GeV 2 already shows a significant change in dσ/dφ for cc. The peak at φ ∼ π is largely erased and only a weak variation with φ can be seen on the log scale for φ < π/2. As k 2 T increases, the enhancement at φ ∼ 0 increases while the peak near φ ∼ π decreases and disappears. A similar but slower evolution is seen for bb in Fig. 10(b). Note that k 2 T = 0.02 GeV 2 still largely follows the k 2 T = 0 results, as may be seen in the dotted curve. Thus sufficiently small steps in k 2 T would likely reveal a gradual change in dσ/dφ.
The relative changes in the distributions care seen more clearly in Fig. 11 where the ratios of the azimuthal distributions relative to that with k 2 T = 0 are shown. The charm quark ratios go from relatively flat at φ < π/2 to an almost linear decrease with the maximum k 2 T , with ∆ = 1 corresponding to the red curves in Figs. 4(a) and 6(a). The ratio of the distributions pivot around φ ∼ 1 with a smooth evolution of the ratios with ∆.
On the other hand, the bb ratios pivot around φ ∼ 2, from a slight modification at φ ∼ π for k 2 T = 0.02 GeV 2 to a flat ratio for φ < 2 followed by a steep descent for φ > 2.
The difference in the location of the pivot point of the distributions is due to the heavy quark mass. The lighter mass of the charm quarks makes it more likely that their relative momentum will become more isotropic at low p T as ∆ increases. Thus the point around which the slope of the ratio is changing is closer to φ ∼ 0. The bottom quarks, more than three times more massive than the charm quarks, are less likely to have their azimuthal distributions become completely isotropic at low p T since the k T kicks are less effective on heavier quarks. Thus there is a reduced peak at φ ∼ π for bb and the slopes of the ratios pivot at an angle closer to φ ∼ π. Figures 12 and 13 show dσ/dφ for p T > 10 GeV at midrapidity and the ratios relative to dσ/dφ with k 2 T = 0. In this case, an enhancement is seen at φ ∼ 0 for both charm and bottom although the enhacement is larger for charm (note that the scales on the y-axes of Fig. 13 are the same). The shapes are consistent with the solid red distributions in Fig. 9.
It is not possible to distinguish between different values of ∆ for cc, even by examining the ratios in Fig. 13(a). Although the ratios increase with ∆ at φ ∼ 0, the curves are almost indistinguishable at low φ. Larger differences in the cc ratios as φ → π can be observed since the ratios at φ > 2.5 increase at slightly lower φ for higher values of k 2 T . This is because k 2 T 1/2 < p T for all p T , resulting in a negligible effect for charm. There is also only a small difference between the bb ratios in Fig. 13(b). However, the bb ratios are visibly separated at φ ∼ 0 because k 2 T 1/2 < p T remains true for the heavier bottom quarks, k 2 T 1/2 /p T ∼ 0.1 − 0.2 for bottom, rather than k 2 T 1/2 /p T ∼ 0.05 − 0.1 for charm. In addition, the harder p T distribution for bottom means that a smaller fraction of the total bb cross section is contained the region p T < 10 GeV than for charm where ∼ 98% of the charm p T distribution is contained in the range p T < 10 GeV.
Finally, it is worth noting that the azimuthal angle distributions and their ratios at low p T (p T < 10 GeV) and higher p T (p T > 10 GeV) at forward rapidity are effectively identical to those at midrapidity. Thus, these results are not shown. Results for p T cuts higher than p T > 10 GeV are not shown because the effect of k T broadening on higher p T cuts would be negligible.

V. COMPARISON TO DATA
There are data on heavy quark pair correlations to test these calculations. There are charm pair correlation data from ALICE [36] and LHCb [37] and bottom quarkbottom jet correlations from CMS [38], both from p + p collisions at √ s = 7 TeV. In addition, there are D 0 D * − and D + D * − data from p + p collisions at √ s = 1.96 TeV from CDF in Tevatron Run II [39].

11
A. LHCb CC Data LHCb measured cc, cc, and (c + c)J/ψ correlations in p + p collisions at 7 TeV for 2 < y < 4 and 3 < p T < 12 GeV. The discussion here will focus on comparison to some of the cc correlations, in particular on the D 0 D 0 , D 0 D − and D + D − combinations where fragmentation functions are known. Thus final states that include a D s or Λ c are excluded. In addition, the finalstate cc data measured by LHCb are more likely to arise from production of a single cc pair. Figure 14 shows the azimuthal angle distributions for the same values of ∆ (from Eq. (9)) employed in the previous section. The results for φ < 2 are in good agreement with the data and with the results of Fig. 12(a) for p T > 10 GeV. However, as φ → π, there is a strong effect on the peak value as a function of ∆. The peak in this region gets broader and lower as ∆ increases. The larger values of ∆ are closer to the data although the data show only a slight trend upward in this region rather than a second peak at φ ∼ π. The calculations at large φ differ strongly with those at p T > 10 GeV in Fig. 12(a), likely because at p T ∼ 3 GeV a significant contribution from k T broadening still remains, see Fig. 4. The sharp peak near φ ∼ π is somewhat artificial since, at k 2 T = 0, all the peak is essentially contained in a single ∆φ bin. When k T broadening is included, the peak region is wider and is shared among several bins although, for p T values of a few GeV, not including p T < 2 p T ∼ 2m D , it is not completely washed out as it is for the results with p T < 10 GeV, dominated by p T ∼ 0, in Fig. 10. Thus the overall agreement of the calculations with the LHCb φ distribution data can be considered quite good.
Calculations of the DD pair invariant mass in the LHCb acceptance with the data are shown in Fig. 15. The agreement of the calculations with the data are very good. A similar minimum is seen in both for M ∼ 2 m 2 D + p 2 T ∼ 7 GeV, using the lower limit of the p T range for the detected D mesons. Below this value of M , the data have a higher cross section than the calculations. The gap in the calculated mass at M ∼ 7 GeV for k 2 T = 0 fills in as ∆ increases. In addition, at higher masses, the calculated distribution falls off somewhat faster for larger ∆ so that ∆ = 1 agrees best with the data.
Finally, Fig. 16 compares the data with the calcualated |∆y| distribution. In this case, the calculated results are independent of the k T broadening. The data fall off somewhat faster than the calculations in all cases.
Overall, the calculated results are in good agreement with the assumption that the D and D mesons both arise from the production of a single cc pair. Discrepancies between the calculations and the data may appear because the calculations are based on exclusive QQ production and include only real corrections at next-toleading order in addition to vertex corrections. No soft radiation from parton showers, as in event generators, which might smooth out the ∆φ distribution at φ ∼ π, are included in these calculations, only fragmentation. In addition, HVQMNR is a negative weight Monte Carlo so that additional weight is placed on the last φ bin, φ = π. A positive weight next-to-leading order exclusive Monte Carlo, such as POWHEG, can be coupled to PYTHIA or HERWIG and include soft initial or finalstate radiation. If the multiplicities of the p+p events are sufficiently large, heavy quark scattering with produced particles may also weaken the azimuthal correlation. See Ref. [6] for a comparison of these data to calculations in the k T -factorization approach.
Note also that, while most of the events are likely to arise from production of a single cc pair, the cc and (c+c)J/ψ yields constitute about 10% of the total charm hadron pair yields [37]. Because only the cc pairs can be produced in a single hard scattering, these other events are more consistent with production through double parton scattering or independent production of two cc pairs. The azimuthal and rapidity distributions for the cc and cJ/ψ events are consistent with isotropic production [37], as might be expected from two hard scatterings.
In addition the ALICE Collaboration made an analysis of azimuthal correlations between reconstructed D mesons and a light hadron trigger particle in p + p collisions at 7 TeV and p+Pb collisions at 5.02 TeV [36]. The light hadrons were primary particles, emitted from the collision points. These particles include those from heavy flavor decays, such as the unreconstructed partner D meson. The data were binned according to the transverse momentum of both the D meson and the light hadron. The minimum light hadron p T was soft, p T > 0.3 GeV. These data were further subdivided into two p T ranges, 0.3 < p T < 1 GeV and p T > 1 GeV. The D meson p T was considerably higher: 3 < p T < 5 GeV, 5 < p T < 8 GeV, and 8 < p T < 16 GeV. To improve statistics, the "D meson" is an average over the D 0 , D + and D * + . The ALICE measurements cover the central region, |y| < 0.5 for the D and |∆η| < 1 for the light hadron. The general behavior is, however, the same as the LHCb D 0 D 0 pairs, a peak at ∆φ = 0 and a smaller enhancement at ∆φ = π. The peak at ∆φ = 0 increases with increasing trigger particle p T and also with increasing D meson p T consistent with the trends of these calculations. The data were compared to simulations with various PYTHIA tunes and also POWHEG+PYTHIA. All simulations reproduced the trends of the data [36], consistent with what one might expect from the results shown here. New AL-ICE data on D-hadron correlations in p + p collisions at 13 TeV were presented recently [40]. In these higher energy collisions, the same models describe these data as well.

B. CDF DD Data
The CDF Collaboration studied charm hadron correlations in p+ p collisions at √ s = 1.96 TeV [39]. CDF combined their single inclusive charm hadron data to look for fully reconstructed charm quark pairs. Only D * ± were considered as candidates for the second charm hadron, in the decay chain D * ± → (D 0 /D 0 )π ± → (Kπ)π ± , so that the mass difference between the D * ± and the first charm hadron, D 0 /D 0 or D ± , ∆m = m(Kππ) − m(Kπ), can suppress the combinatorial background. Several thousand D 0 D * − and D + D * − pairs were constructed in this way.
Their goal was to try to better understand the underlying QQ production process [39]. As discussed in Sec. II C, commonly employed high energy event generators such as PYTHIA implement prompt QQ production by three different leading-order porcesses: 'pair creation', as in the LO diagrams in Fig. 1; 'flavor excitation', and 'gluon splitting', corresponding to the NLO processes shown in Fig. 2 (b) and (c) respectively.
The weighting of these different contributions in PYTHIA can be tuned to match the data to determine which processes are most important to heavy flavor production in the phase space covered by an experiment. According to this interpretation, the final-state QQ configurations are treated separately by virtue of their topology with φ = π, corresponding to pair creation, and φ = 0, corresponding to gluon splitting. This apporach stands in contrast to the NLO calculations shown in this paper where all NLO diagrams are summed coherently with weighting according to color and spin with interference terms. As described in Sec. II C, at NLO, flavor excitation and gluon splitting are both part of production by gg interactions in the initial state. They are not separate production processes and are not treated as such.
The data on each final state, in the rapidity interval |y| < 1 and p T ranges 5.5 < p D 0 ,D * − T < 20 GeV, 7 < p D + T < 20 GeV, are shown in Fig. 17. In the calculations, the y interval was the same but the lower limit on the D meson p T was taken to be 5.5 GeV in all cases. The calculations agree rather well with the data. The results are independent of ∆ except when the pairs are completely back-to-back (|∆φ| = π), as might be expected for the given p T range. In Ref. [39], the comparison with the separate LO 'production mechanisms' in PYTHIA showed that, at least for the PYTHIA tune used by the CDF Collaboration, the 'gluon splitting' or collinear contribution (∆φ = 0) to production was underestimated. In the calculation, with the NLO contributions added according to initial-state parton channel (gg, (q/q)g and qq), there is no significant discrepancy.

C. CMS Bottom Quark-Bottom Jet Data
The CMS Collaboration measured bb correlations through the use of their jet trigger [38]. The analysis employed the single jet trigger with calorimeter energy above trigger thresholds of 15, 30 and 45 GeV. The energy scale at which these three triggers are greater than 99% efficient corresponds to transverse momentum of the leading jet, p jet T , of 56, 84 and 120 GeV respectively. The leading jet is associated with one of the b quarks. The triggered events are required to have one reconstructed jet with the minimum p T defined as above and a reconstructed primary vertex. The leading jet is requred to be within the pseudorapidity interval |η jet | < 3. The jet along with identified B hadrons with p B T > 15 GeV and in the interval |η B | < 2, is considered to originate from a single bb pair. Thus the events also have two reconstructed secondary vertices. The angular correlations are calculated using the B hadron and b-jet flight directions.
The CMS results were presented as a function of the azimuthal angle difference, ∆φ, and a separation variable, ∆R, that includes the difference in polar angles between the B hadron and jet, given in terms of the difference in pseudorapidity, ∆η, ∆R = (∆η) 2 + (∆φ) 2 . The CMS predictions for the behavior of ∆R are based on the notion that gluon splitting dominates at low ∆R while flavor creation is most important at large ∆R. This is similar to what one expects from event generators for ∆φ, as discussed in Secs. II C and V B.
The calculations shown here are based on the assumption that a bb pair is produced with one final state B or B in the kinematic region p B T > 15 GeV and rapidity |y B | < 2. The jet is assumed to originate from the second B hadron of the pair with the minimum p T of the B equal to the p T of the jet, greater than 56, 84 and 120 GeV respectively, while η jet is assumed to be equivalent to the rapidity of the b quark. Since the pseudorapidity difference is replaced by a rapidity difference in the calculation of ∆R and the equivalence of the b quark and the jet p T might not be exact, one might expect that the calculation of the ∆φ distribution might be better reproduced than that of ∆R.
The comparisons of the calculations with the data are shown in Figs. 18 and 19. While the calculations are done with all values of ∆ in Eq. (9), no difference between them is visible except at large ∆R. As in Ref. [38], the data and calculations for the different minimum jet p T are scaled by factors of 4, 2 and 1 for p jet T greater than 56, 84 and 120 GeV respectively to more clearly separate the results.
The comparison to the ∆φ distributions is very good, as may be expected. The comparison to ∆R is not as good at low ∆R where the CMS Collaboration suggests gluon splitting is dominant. However, given the agreement with the low ∆φ results, where this mechanism is also expected to be dominant in event generators, the underestimate may be due to the calculation of ∆R with ∆y instead of ∆η. Note that the large ∆R results are best reproduced with ∆ = 1 in Eq. (9).
VI. COLD NUCLEAR MATTER: p+PB COLLISIONS AT 5.02 TEV Finally, the interaction of cc pairs in cold nuclear matter, such as p+Pb collisions at the LHC, are discussed. It has been suggested [3] that energy loss by heavy quarks in heavy-ion collisions could change the azimuthal correlations. First, it must be determined how the charm distributions and their azimuthal separation are influenced by the presence of cold nuclear matter. For example, the effect of cold matter energy loss, which could be manifested as additional k T broadening by multiple scattering in the nucleus, could result in changes in the p T distributions of heavy quarks in p+Pb relative p + p collisions. This would be in addition to modification of the parton densities in the nucleus, referred to as shadowing.
Here the results of shadowing alone on the single charm meson p T distributions and cc pair azimuthal distributions are compared to shadowing and k T broadening, both with k 2 T ∼ 1.5 GeV 2 alone, as in p + p collisions, and with an additional k T kick due to the presence of the nuclear medium in p+Pb collisions.
The calculations of the nuclear modification factor R pPb (p T ) are compared to data from ALICE [41]. AL-ICE measured prompt production of D 0 , D + , D * + , and D − s as well as their charge conjugates at central rapidity, −0.96 < y cms < 0.04 in p+Pb collisions at √ s = 5.02 TeV. The yields were compared to those in p + p collisions to form R pPb (p T ). Note that, for the analysis in Ref. [41], the p + p baseline at 5.02 TeV was interpolated between the 2.76 TeV and 7 TeV data since it was completed before p + p data were available at 5.02 TeV.
In Ref. [41], the average of the D 0 , D + and D * + were compared to several calculations, including nextto-leading order calculations with HVQMNR including only shadowing effects with no momentum braodening nor fragmentation. The uncertainty band gives suppression below p T ∼ 5 GeV, rising to be equivalent with unity at higher p T . See Ref. [41] for details of the other comparisons which included a leading order calculation of shadowing due to power corrections, k T broadening and cold matter energy loss [42], and a calculation assuming color glass condensate in the initial state [43].
To go beyond p + p collisions, in the NLO calculations, the proton parton densities must be replaced by those of the nucleus. If A is a nucleus, the nuclear parton densities, f A i (x 2 , µ 2 ), can be assumed to factorize into the nucleon parton density, f p i (x 2 , µ 2 F ), independent of A; and a shadowing ratio, S i (A, x 2 , µ 2 F ) that parameterizes the modifications of the nucleon parton densities in the nucleus. Here x 2 is the fraction of the nucleon momentum carried by the interacting parton in the nucleus. While the shadowing effect may also depend on the impact parameter, b, of the parton from the proton with the lead nucleus, only minimum bias results, independent of impact parameter, are shown here.
The calculations in Fig. 20 employ the EPS09 NLO [44] parameterization for S i , assuming collinear factorization. The EPS09 sets include 15 parameters, giving an error set of 30 additional parameterizations created by varying each parameter within one standard deviation. The EPS09 uncertainty band is obtained by calculating the deviations from the central value for the 15 parameter variations on either side of the central set and adding them in quadrature. There is a more recent shadowing parameterization by Eskola et al., EPPS16 [45], that includes some of the LHC data from the first p+Pb run at 5.02 TeV in the global analysis of the nuclear parton densities. The central result for the gluon distribution is very similar to that of EPS09 so that the central result with EPPS16 should be very similar to that shown here. However, it employs 5 additional parameters, thus resulting in a larger uncertainty band than the one in Fig. 20.
The calculations given as a function of p T in Fig. 20 are in the same rapidity interval as the ALICE data. Three bands are shown. The first, with k 2 T = 0 and ǫ P = 0, in blue, corresponds to the EPS09 calculation in Ref. [41] albeit the curves in that paper likely employed a larger quark mass, m ∼ 1.5 GeV, which would result in a reduced shadowing effect at p T ∼ 0 relative to the calculation shown here. The red histogram, employing the preferred parameters found for D mesons in this paper, and assuming the same intrinsic k T in p+Pb and p + p collisions, gives a nearly identical result to the blue histogram, calculated without either effect.
However, one might expect that a higher intrinsic k T is required in a nuclear medium relative to that in p + p due to multiple scattering in the nucleus, known as the Cronin effect [46]. If one assumes that k 2 T is doubled in the cold nuclear medium, then the magenta curves are obtained. The larger k 2 T results in an increase of the band at intermediate p T , with the upper limit of the band increasing above unity at p T ∼ 5 GeV. The statistics on the averaged D meson data are not sufficiently significant to distinguish between the results: all are consistent with the data.
The results for the nuclear modification factor of the azimuthal angular correlations between the heavy quarks are shown in Fig. 21 for the same three sets of calculations as in Fig. 20. In this case there are no cuts made on p T or rapidity.
When no intrinsic k T is included, R pPb (φ) is independent of φ. There is a slight φ dependence introduced when the p + p value of k 2 T is employed. It is, however, a small deviation relative to the calculation with no k T effect. There is a more significant effect, as expected, when the intrinsic k T effect is doubled in p+Pb relative to p + p collisions. Instead of being independent of φ, R pPb (φ) now decreases as φ increases, a result of the increased k 2 T , as in Fig. 7(a).
In heavy-ion collisions, in particular at heavy-ion colliders, multiple cc pairs are produced in a single event. Therefore, for an analysis of heavy-flavor azimuthal correlations in hot matter it is necessary to find the pair vertex to ensure that the c and c originate from the same hard scattering.

VII. SUMMARY
The similarity of the results for p T > 10 GeV, even for k 2 T = 0, shows that the enhancement at φ ∼ 0, with the QQ pair aligned opposite a hard parton, is independent of k 2 T and arises instead from the QQ production mechanism at high p T . Thus the high p T behavior of dσ/dφ is indicative of the contribution of next-to-leading order production while the low p T behavior of dσ/dφ is extremely sensitive to the chosen k 2 T and essentially independent of fragmentation.
It appears that, so far, the open heavy flavor results presented by the LHC collaborations, both single inclusive production and charm pair correlations, are in agreement with calculations based on collinear factorization with single hard scatterings. The exception, the cc and cJ/ψ events at LHCb, are consistent with double parton scattering.
Thus, for p T cuts on the order of a few GeV, the calculations of the azimuthal angle distributions are rather insensitive to fragmentation and k T broadening which affect the correlations at low p T . Thus, hot nuclear matter effects on these correlations should be rather robust for p T ≥ 3 − 5 GeV. sions are compared to data from LHCb [19] at 2.5 < y < 3 and ATLAS [32] at |η| < 2.5 respectively. The curves show the extent of the uncertainty bands calculated employing Eqs. (7) and (8). The HVQMNR code (blue dashed curves) utilizes ( k 2 T (GeV 2 ), ǫP ) = (1.5, 0.008) for charm and (3,0.0004) for bottom. The corresponding FONLL uncertainty band (red curves) is also shown. The same quark mass and scale parameters are used in both calculations.     [39]. The data are compared to calculations in the same acceptance with k 2 T = 0 and for values of ∆ from −3/2 to 1 in Eq. (9).