$b \bar b$ Kinematic Correlations in Cold Nuclear Matter

Background: The LHCb Collaboration has studied a number of kinematic correlations between $B$-hadron pairs through their subsequent decays to $J/\psi$ pairs at 7 and 8 TeV for four minimum values of the $J/\psi$ $p_T$. Purpose: In this work, these measurements are compared to calculations of $b \bar b$ pairs and their hadronization and inclusive decays to $J/\psi J/\psi$ are compared to the same observables. Potential cold matter effects on the $b \bar b$ pair observables are discussed to determine which are most likely to provide insights about the system and why. Methods: The calculations, employing the exclusive HVQMNR code, assume the same intrinsic $k_T$-broadening and fragmentation as in [R. Vogt, Phys. Rev. C {\bf 98} (2018) 034907]. The pair distributions presented by LHCb are calculated in this approach, both for the parent $b \bar b$ and the $J/\psi J/\psi$ pairs produced in their decay. The sensitivity of the results to the intrinsic $k_T$ broadening is shown. The theoretical uncertainties due to the $b$ quark mass and scale variations on both the initial $b \bar b$ pairs and the resulting $J/\psi$ pairs are also shown. Possible effects due to the presence of the nucleus are studied by increasing the size of the $k_T$ broadening and modification of the fragmentation parameter. Results: Good agreement with the LHCb data is found for all observables. The parent $b \bar b$ distributions are more sensitive to the $k_T$ broadening than are the final-state $J/\psi$ pairs. Conclusions: Next-to-leading order calculations with $k_T$ broadening, as in [R. Vogt, Phys. Rev. C {\bf 98} (2018) 034907], can describe all correlated observables. Multiple measurements of correlated observables are sensitive to different nuclear effects which can help distinguish between them.


I. INTRODUCTION
Heavy flavor pair production has long been of interest in elementary p + p and p + p collisions as a way to test perturbative QCD. Measurements of heavy flavor correlations provide information about how heavy quark pairs are produced in perturbative QCD, indeed much more information that can be gained from single inclusive heavy flavor production alone. In the case of bb production, measurements of pair observables can improve measurements of B 0 − B 0 mixing [1]. Finally, a good understanding of multiple correlated observables provides a better baseline for their production and modification in heavy-ion collisions. Correlated bb studies have been carried out at hadron colliders. The first measurements were in p + p collisions and carried out at the CERN SppS, UA1 ( √ s = 0.63 TeV) [1] and Fermilab Tevatron, D0 ( √ s = 1.8 TeV) [2], and CDF ( √ s = 1.8 TeV [3] and √ s = 1.96 TeV [4]). These measurements were primarily through studies of lepton pairs. The backgrounds for these measurements include cc decays, Drell-Yan dileptons and leptons from light meson decays. The light hadron decay leptons can be removed by like-sign subtraction. The Drell-Yan rate is much lower than the heavy flavor production rate because Drell-Yan is an electroweak process. In addition, if relatively high p T leptons are selected, the charm rate will be suppressed. CDF [3] and, more recently, at the LHC, ATLAS [5] studied B hadron pair production through J/ψ+lepton final states in p + p collisions at √ s = 1.8 TeV and p + p collisions at √ s = 8 TeV respectively.
More recent measurements of heavy flavor contribu-tions to low mass dilepton production in p + p collisions have been reported by PHENIX at the BNL relativistic heavy-ion collider ( √ s = 0.2 TeV) [6] and ALICE at the LHC ( √ s = 7 TeV [7] and 13 TeV [8]). In these measurements, the low mass kinematic region of interest makes the cc contributions larger so that both contributions have to be taken into account in the analysis.
Previous dilepton analyses [1][2][3]5] have generally focused on tests of NLO calculations, usually in conjunction with diagrams of different topologies in a leading order event generator such as ISAJET [9], HVQJET [10], HERWIG [11] or PYTHIA [12]. For example, UA1 favorably compared calculations at next-to-leading order made with the HVQMNR code [13] to ISAJET as long as the ISAJET "higher" order contributions (flavor excitation and gluon splitting) were included. They tried to separate different topologies but the results were rather inconclusive due to momentum requirements and low statistics. They did, however, determine that the NLO contribution was at least 40% of the measured cross section, depending on the muon p T [1]. The D0 Collaboration compared their data to calculations with HVQ-JET. They showed that their measured azimuthal separation, |∆φ|, between the decay leptons was compatible with HVQJET with higher order corrections and not with the leading order contributions alone [2]. The CDF measurement at √ s = 1.8 TeV [3] compared their final-state J/ψ+lepton data to HVQMNR calculations as well as to PYTHIA and HERWIG simulations. They found that ≈ 25% of the bb production as a function of azimuthal separation was found at |∆φ| < 90 • , suggesting the importance of higher-order corrections. CDF also studied the dependence of the |∆φ| distribution to the bottom quark mass, factorization and renormalization scales, and the intrinsic k T . They found that mass and scale variations could alter the magnitude of the cross section but not its shape: changing the shape of the distribution was only possible by modifying the k T [3]. These findings are in accord with the hadron-level studies of bb correlations in Ref. [14]. ATLAS compared their final-state J/ψ+lepton results with several event generators, finding good agreement between the simulations and the data [5]. The low mass dilepton measurements of PHENIX attempted to separate the dilepton data into cc and bb components, as well as separating the heavy flavor cross sections into their topological components, as if they were independent production mechanisms [6].
Most of the measurements mentioned so far have focused on the central rapidity region. The LHCb Collaboration [15] is the first to study bb correlations through J/ψJ/ψ final states as well as at forward rapidity. While such a measurement is less direct than reconstructed D or B mesons, as discussed in Ref. [14], along with comparison to DD pairs measured by CDF at √ s = 1.96 TeV [16] and by LHCb at √ s = 7 TeV [17] and B hadron-b jet pairs measured by CMS at √ s = 7 TeV [18], it allows a more straightforward comparison to calculations than the dilepton decay channel.
In this work, the model of QQ production developed in Ref. [14], with modified fragmentation function parameters and k T broadening, is employed to study bb → J/ψJ/ψ pair production. The measurement is discussed in more detail in Sec. II. In Sec. III, the model employed for bb production is briefly discussed. The pair observables are calculated both for the initial bb production and the final J/ψJ/ψ pairs, assuming the same minimum p T for both the parent B meson and the decay J/ψ in Sec. IV. The results are compared to the LHCb observables and their dependence on the experimental p T cut is shown. Section V shows how neglecting k T broadening affects the calculated observables. The mass and scale dependence of the results and how they change with p T cut is discussed in Sec. VI while the dependence of the results on the rapidity range of the measurement is shown in Sec. VII. Section VIII describes possible nuclear effects on the correlations. The work is then summarized in Sec. IX.

II. LHCB MEASUREMENT OF bb → J/ψJ/ψX
LHCb reconstructed two J/ψs from their decays to dimuons in the forward rapidity region, 2 < y < 4.5. The two J/ψs must be associated with the same primary vertex to ensure that they came from the same collision. The J/ψs were also required to be displaced from their primary vertex to be b-decay candidates. This requirement essentially eliminated prompt J/ψs from different collisions as well as events with two prompt J/ψs and associated J/ψ and b quark production.
They chose different minimum J/ψ transverse momenta, p T , to study the effect of an increasing minimum p T on the pair correlations. The data from proton-proton collisions at √ s = 7 and 8 TeV were combined for greater statistics. Because the shapes of the distributions at this energy are independent of √ s even if the integrated production cross sections differ, the results were presented as (1/σ)dσ/dX where X refers to the pair observable. This way of displaying the data makes it easier to compare the shapes of the distributions with different minimum p T .
LHCb presented results for six pair observables, |∆φ * |, the difference in azimuthal angle between the b and b mesons; |∆η * |, the difference in pseudorapidity between the b and b mesons; A T , the asymmetry between the transverse momenta of the J/ψs; and the mass, M , transverse momentum, p Tp , and rapidity, y p of the J/ψ pair. The first two observables, |∆φ * | and |∆η * |, were assumed to be directly related to the parent b mesons because φ * and η * were estimated from the direction of the vector from the primary vertex to the J/ψ decay vertex [15]. They also included, in an appendix, the distributions |∆φ|, |∆η|, and |∆y|,the differences in azimuthal angle, pseudorapidity, and rapidity respectively between the J/ψ mesons themselves. In this work, |∆y| is presented rather than |∆η| for the parent bb. All the pair observables studied by LHCb will be calculated for both the parent bb mesons and the subsequent J/ψJ/ψ decays.
In Ref. [15], the LHCb Collaboration compared their data to PYTHIA [19,20] and POWHEG [21] calculations as well as simulations of uncorrelated bb production [22,23] based on the transverse momentum and rapidity distributions for single b → J/ψX decays measured by LHCb. They noted that the pair distributions generated by both PYTHIA6 [19] and PYTHIA8 [20] were identical and thus the results from the two simulations were combined in their comparison to the data.
As in a number of the previous bb measurements analyzed via dilepton decays [1][2][3], LHCb looked for evidence of different topological contributions to heavy flavor production in their data: gluon splitting, flavor excitation and flavor creation.
As discussed in detail in Ref. [14], these artificial designations are not indicative of different production mechanisms but of distinct diagram topologies at leading order (flavor creation) and next-lo-leading order (including gluon splitting and flavor excitation). These processes are distinguished by having two (flavor creation), one (flavor excitation), or no (gluon splitting) heavy quarks in the hard scattering. Heavy quarks not involved in the hard scattering are produced in an initial-or final-state parton shower [24]. This separation is necessary because PYTHIA includes only LO diagrams. While all three topological contributions are part of the gg → QQX production process, when they are treated as individual components, not all NLO production diagrams are actually included (such as virtual corrections) and the interferences between diagrams are not accounted for. However, none of these contributions to bb production constitute a new production mechanism. The implementation of heavy flavor production in PYTHIA is more completely described in Ref. [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, see for example Ref. [25] for more detail. However, such tuning may change the relative contributions of distinct topologies from the same initial state in PYTHIA relative to a NLO calculation. Double counting of these processes is avoided by requiring that the hard scattering should be of greater virtuality than the parton shower [25]. The parton showers also effectively provide a leading-log resummation of light emissions while keeping the pair distributions finite over all phase space.
POWHEG, a NLO generator using PYTHIA for hadronization and decay [21], does not separate these topologies in the same way that PYTHIA does, all diagrams, with their interference terms, are included. In Ref. [15], they conclude that, because POWHEG and PYTHIA both describe the data, NLO effects on bb production are small. They also note that, aside from the |∆φ * | distributions, the data are consistent with uncorrelated bb production. They reach this conclusion by suggesting that gluon splitting is a small contribution to bb production. However, it is not feasible to separate this diagram from all other NLO contributions because it interferes with the amplitudes of other gg diagrams.
The conclusion the NLO contributions to bb production are small, reached by the LHCb Collaboration in Ref. [15], can be tested by comparison to a NLO calculation of both the bb and J/ψJ/ψ final states. This comparison, in Sec. IV, is at the center of this work.

III. MODEL DESCRIPTION
The calculations here, using the HVQMNR code [13] designed to calculate QQ pair production at NLO, follow those outlined in Ref. [14]. The bottom quark mass, m b , factorization scale, µ F , and renormalization scale µ R and their uncertainties were set by comparison to the bb total cross section data with m b = 4.65 ± 0.09 GeV, µ F /m = 1.40 +0.77 −0.49 , and µ R /m = 1.10 +0.22 −0.20 [14]. Hadronization was accomplished through the use of the Peterson fragmentation function [26] and k T broadening. The value of ǫ P , the Peterson function parameter, was set by comparison to the FONLL B meson p T distribution in Ref. [14] while k 2 T , the average broadening, was fixed previously by comparing the measured Υ p T distributions to calculations of Υ production in the color evaporation model. Here, for bb production, where the parameter ∆ was introduced to study the sensitivity of the azimuthal correlations to the amount of broadening. The value ∆ = 1 is the default value [14], resulting in k 2 T ≈ 3 GeV 2 for √ s = 7 TeV. For further details on the determination of ǫ P and the sensitivity of the QQ results to the magnitude of k 2 T , see Ref. [14]. Note that it is necessary to use a code such as the exclusive HVQMNR calculation because the QQ pair quantities are calculable in such an approach while only single inclusive quark distributions are so far available with the FONLL [27] and generalized mass -variable flavor number approaches [28,29].
The HVQMNR code [13] uses negative weight events to cancel divergences numerically. Without a k T kick there can be a mismatch in the cancellation, leading to a negative value of the cross section for pair p T at p Tp = 0 and azimuthal separation at φ = π, as can be seen in some of the bb distributions with k 2 T = 0. 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; the broadening plays this role in the code at low p T [30]. Open charm results at fixed-target energies required transverse momentum broadening to obtain agreement with the data after fragmentation was applied [31]. Such broadening was also used as a proxy for resummation in Drell-Yan production, see e.g. Refs. [30,32,33].
In HVQMNR, the kick is added in the final state using the Gaussian function g p (k T ) [31], which multiplies the parton distribution functions, assuming the x and k T dependencies factorize. As explained in Ref. [31], it does not matter whether the k T dependence is added in the initial or final state as long as the kick is not too large. The application of the k T is now described. The QQ system is first boosted to the rest frame from its longitudinal center-of-mass frame. The 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. (3). The quarks are then boosted out of the pair rest frame, changing the initial transverse momentum of the hard scattering from p T to p T + k T 1 + k T 2 . While the k T is here applied to the QQ pair, it could have alternatively been given to the entire final-state system, including the light parton in 2 → 3 processes, as if it were applied directly to the initial state. The two methods of introducing k T are equivalent if the calculation is LO but at NLO a light parton in the final state can make the correspondence inexact.
The effect of a k T kick on p T -related distributions (p Tp , M , A T ) should decrease as √ s increases because p T also increases with √ s. Because k 2 T is assumed to increase with √ s, see Eq.
(2), the effect is most important at low p T . The effect of a k T kick also decreases with increasing quark mass, as shown in Ref. [14], requiring a larger k 2 T for bottom quarks relative to charm quarks to have a nonnegligible effect on bottom production at higher energies.
Although LHCb suggested in Ref. [15] that the similarity of the PYTHIA and POWHEG simulations are indicative of a small NLO contribution, it is important to recall that gluon splitting is an explicit contribution to gg → QQX at O(α 3 s ) and thus a real NLO contribution, as is flavor excitation. As previously discussed, it is not feasible to separate individual diagrams since such a procedure no longer allows for interferences between diagrams. The LO flavor creation contributions, gg → QQ and qq → QQ, only produce back-to-back QQ pairs, a delta function for |∆φ * | = π, A T = 0 and p Tp = 0 without broadening. The NLO contributions are modeled in PYTHIA by flavor excitation and gluon splitting. These contributions have a significantly weaker ∆φ dependence in PYTHIA. Flavor excitation is weakly enhanced at |∆φ| ≈ π while, since gluon splitting generally produces collinear QQ pairs, it results in a weak enhancement at |∆φ| ≈ 0 [24]. These contributions and summed together with flavor creation at O(α 2 s ), without interference terms but with parton showers.
The introduction of k T broadening at NLO in HVQMNR softens and widens the peak at |∆φ| = π for bb with a finite tail as |∆φ| → 0. It does not produce a significant enhancement at |∆φ| → 0 as it does for charm pairs at similar values of p T because k 2 T < m 2 b while, for charm, k 2 T ≈ m 2 c [14]. The effect of k T broadening also depends strongly on the quark momentum, as discussed for the |∆φ| distributions in Ref. [14].
Observables related to the rapidity, either the rapidity difference or the pair rapidity, should be independent of the broadening. However, the other pair observables studied by LHCb should be affected by broadening, at least for the parent b quarks. Thus the calculations here compare results with and without k T broadening on both the initial bb pairs and the final state J/ψJ/ψ pairs.

IV. COMPARISON TO THE LHCB DATA
In this section, the pair quantities, |∆φ|, |∆y|, y p , A T , p Tp and M are calculated and compared to the LHCb data. The bb pair distribution include both fragmentation and k T broadening as described in the previous section. The J/ψJ/ψ pair distributions are calculated with the B → J/ψX decay with a 1.094% branching ratio [34].
All results are shown for the minimum p T cuts of 2, 3, 5 and 7 GeV on the J/ψ and the parent B meson. Note that the J/ψs from B decay would, of course, generally arise from parent B mesons with p T larger than those of the final-state J/ψ.
The LHCb data are also shown for comparison on each plot. All quantities are divided by the total cross section, (1/σ)(dσ/dX) ≡ d ln σ/dX where X denotes the observables on the y-axes of the plots, so that the 7 and 8 TeV LHCb measurements can be combined for improved statistics. Note that even though the LHCb data shown here are from the √ s = 7 and 8 TeV runs combined, the ≈ 15% difference in √ s between the two data sets gives only a 1-2% change in k 2 T based on Eq. (2). Given the small change in k 2 T for p+p collisions and the uniform shapes of d ln σ/dX, the calculations compared to the data in this section are all done for √ s = 7 TeV. It was verified that the normalized pair distributions calculated here remain unchanged at √ s = 7 and 8 TeV in p + p collisions.
Note that, in the calculations, the J/ψs from B decays have lower statistics than the parent B mesons, especially as the minimum p T increases. Thus red histograms are generally used to represent the J/ψ pair quantities while black curves are used for the bb pair distributions. The LHCb results for J/ψ pairs are rendered as red points while the reported bb quantities are given as black points.
A. |∆φ| and |∆φ * | LHCb presented |∆φ| distributions for both the initial B meson pair, reported as bb in Fig. 1 and in all the figures in this section, and the J/ψ pairs. Recall that LHCb estimated the azimuthal angle of each B meson from the direction of the vector from the primary vertex to the J/ψ decay vertex. They also determined the azimuthal angles for the J/ψs individually. As shown in Fig. 1, the |∆φ * | and |∆φ| distributions for bb and J/ψJ/ψ respectively are compatible with each other within the uncertainties.
The bb |∆φ * | distribution has a peak slightly below |∆φ * | ≈ π with a flatter distribution as |∆φ * | → 0 relative to that of the J/ψ pair. As the minimum p T grows, the peak near back-to-back (|∆φ * | ≈ π) grows higher and becomes narrower for the bb pairs. Likewise, the distribution at |∆φ * | ≈ 0 starts to increase from approximately flat at low |∆φ * | to a slight enhancement that becomes more pronounced with increasing minimum p T . This is because that, as the minimum p T grows from 2 to 7 GeV, the relative values of k 2 The larger m T allows the development of a double-peaked ∆φ * distribution, more closely connected to diagrams with a high p T bb pair balanced against a hard parton in the opposite direction, such as 'gluon splitting'.
This trend in the increase of (π/σ)dσ/d|∆φ| can especially be seen for the lighter mass J/ψ decay products. In this case, because m J/ψ /m b ≈ 2/3 and the minimum J/ψ p T is generally smaller than that of the parent B meson and the k T kick is applied to the parent meson, not the J/ψ decay product, the enhancement seen in the bb distributions sets in at lower p T for J/ψ pairs and is correspondingly larger. Here T . Because the k T kick is on the bottom quarks as they hadronize rather than on the J/ψ itself, the p T selected is larger relative to the primary B hadron so that the enhancement grows faster with minimum p T for J/ψ pairs, as shown in Fig. 1. If a higher minimum p T were chosen for B mesons, to more closely match the average p T of the B meson producing the minimum J/ψ p T , the enhancement at |∆φ| → 0 would grow larger, closer to matching the peak at |∆φ| ≈ π, as shown for cc correlations with p T > 10 GeV in Ref. [14].

B. |∆y| and yp
The difference in rapidity, ∆y, (or, in the case of the LHCb measurement, ∆η), was determined both for the initial bb pairs and the final-state J/ψ pairs. The pair rapidity, y p , was only determined for the J/ψ pairs. Given the acceptance of the LHCb spectrometer of 2 < y < 4.5, the limits on ∆y is constrained to be in the range 0 < ∆y < 2.5 while the pair rapidity reported by LHCb lies in the range 2 < y p < 4.5.
As is evident from Fig. 2, the |∆y| distribution for bb and J/ψJ/ψ are in good agreement. They decrease from a peak at |∆y| = 0 to 0 at |∆y| = 2.5. The shape of both distributions is more concave than linear but is approximately identical for bb and J/ψJ/ψ. The behavior is also relatively independent of the minimum p T . In the case of the bb pairs, the average |∆y| decreases from 0.75 for p T > 2 GeV to 0.73 for p T > 7 GeV, only a 2% difference. On the other hand, the average values of |∆y| for the J/ψ pairs decreases from 0.74 to 0.69 as the minimum p T increases from 2 to 7 GeV. At the highest minimum p T , the average |∆y| is reduced by 5% for J/ψ pairs relative to bb pairs. The differences, while not significant, are not zero.
The pair rapidity distributions, shown in Fig. 3, exhibit a similarly small decrease in the average y p with increasing minimum p T , a 2% decrease in the average for bb pairs between p T > 2 and > 7 GeV, from 3.07 to 3.00 respectively. There is a 5% decrease in average y p for the J/ψ pairs, from 3.07 with p T > 2 GeV to 2.93 with p T > 7 GeV. This small difference on average is sufficient for a small backward shift between the bb and J/ψJ/ψ curves for p T > 7 GeV, especially given the average p T for the parent B mesons of J/ψs with the same minimum p T .
FIG. 2: (Color online) The rapidity difference |∆y| between the b and b (black dashed curve) and the J/ψ's resulting from the bottom quark decays (red solid curve) are shown compared to the LHCb data [15] (black for bb, red circles for J/ψ pairs) for the pT cuts on the b quarks and the J/ψ of 2 (a), 3 (b), 5 (c) and 7 GeV (d).

C. AT
The p T asymmetry, Fig. 4, would be zero for bb pairs produced in a back-to-back configuration at leading order. At NLO, the pairs are no longer back-to-back and d ln σ/dA T decreases with increasing A T . The A T distribution for the J/ψ pairs is maximal at A T = 0, in accord with the maximum |∆φ| ≈ π. This same relation also results in a steeper A T distribution for higher minimum p T . The distribution goes to zero at A T = 1 if the p T of one of the b quarks or J/ψ mesons is very soft or the final states are in alignment. On the other hand, the bb distributions peak away from A T = 0 due to the inclusion of k T broadening, as discussed further in Sec. V. The peak of the A T distribution is at A T ≈ 0.25 for p T > 2 GeV. As the minimum p T is increased, the distribution for bb pairs becomes narrower with a higher peak, akin to the |∆φ * | distributions shown in Fig. 1. The average value of A T decreases from 0.025 at p T > 2 GeV to 0.17 at p T > 7 GeV.
As previously mentioned, the J/ψ pair A T distribution is maximum at A T = 0 instead of a finite A T , as for bb. At lower minimum p T , the distribution is narrower for the J/ψJ/ψ, with an average of A T ≈ 0.21 for p T > 2 GeV. By the highest minimum p T , p T > 7 GeV, the average is approximately the same for both, A T ≈ 0.16 for J/ψJ/ψ.
The trends for the calculated J/ψ pairs are in very good agreement with the data for all values of minimum p T studied. Note also that, above A T ≈ 0.4, the calculated bb and J/ψJ/ψ distributions are in agreement. tude, the J/ψ pair peak is shifted backward by about 1.7 GeV relative to the bb, as is evident in Fig. 5(a). At p T > 3 GeV, while the average p Tp is still about 1.3 GeV smaller for the J/ψ pairs, most of the difference is at p Tp < 5 GeV. Above this value, the distribution is significantly broader than for p T > 3 GeV.
As the minimum p T is increased, the average values of p Tp for the initial bb and the final J/ψJ/ψ become more similar. In addition, a feature develops at high p Tp , a shoulder in the distribution that appears at effectively twice the minimum p T , independent of whether the calculation is for the initial bb pairs or the decay J/ψ pairs although the statistics at high pair p T is significantly degraded for p T > 7 GeV. The rise of this shoulder appears to correspond to the rise of the peak at |∆φ| = 0 in Fig. 1 where the bb and J/ψ pairs are aligned. In all cases, the calculated J/ψJ/ψ pair distributions agree quite well with the LHCb data.
A similar trend seen for the pair mass distributions in Fig. 6. The minimum bb pair mass is 2m b = 9.3 GeV for m b = 4.65 GeV. Assuming that the p T of both of the individual mesons are equal, the square of the pair mass can be written as M 2 = 2m 2 T (1 + cosh(∆y)). Thus as the minimum single meson p T increases, m T also increases and the average pair mass moves to higher M . This estimate is accurate for 2 → 2 processes but is an underestimate for the 2 → 3 diagrams that dominate NLO bb production. Nonetheless, one can see a clear trend that the bb peak shifts to higher mass with an in- crease in minimum p T , with a residual enhancement at 2m b for the highest minimum p T .
When J/ψ pairs from b decays are considered, the pair mass does not have a specific threshold any longer. As shown for M > 2m b , the J/ψ pair mass is steeply decreasing for p T > 2 and 3 GeV while for p T > 5 and 7 GeV, a peak at higher M also develops. The average mass of the J/ψ pairs for the higher minimum p T is shifted backward by several GeV: compare the black curves and the red histograms in Fig. 6(c) and (d). The calculations of the J/ψ pairs follow the LHCb data very closely.
In general, the calculations presented here are in very good agreement for the pair observables obtained by LHCb for all values of the minimum p T chosen.
In this section, the sensitivity of observables to the presence of an intrinsic k 2 T is explored. While in Ref. [14], the sensitivity was studied by gradually dialing up k 2 T to its default value, here only the results with k 2 T = 0 and the default 3 GeV 2 (at √ s = 7 TeV) are compared. The comparison is made for both the bb pairs and the final J/ψJ/ψ.
Because |∆y| and y p are independent of k 2 T , the comparison is only shown for |∆φ|, M , p Tp and A T in Fig. 7. The left-hand side of the figure shows the results for p T > 2 GeV while calculations for p T > 7 GeV are shown on the right-hand side. The behavior of calculations with p T > 3 and 5 GeV follow similar trends. All results with k 2 T = 0 are given in black, curves for bb and histograms for J/ψJ/ψ.
It is clear that the bb distributions are most affected by the presence of broadening. The peaks at |∆φ * | ≈ π, low p Tp and low A T are enhanced. They are not delta functions without broadening, as they would be at leading order, but have finite tails indicative of a NLO process. There is no significant change in the pair mass distributions, independent of minimum p T .
Note that the bins at |∆φ * | ≈ π, p Tp → 0 and A T → 0 do not go directly to zero but show enhanced peaks due to the incomplete numerical cancellation of divergences with HVQMNR, as discussed in Sec. III. For example, with k 2 T = 0, the normalized A T distribution would still have a peak at finite A T but it would be closer to A T ≈ 0 and decrease faster with A T . The addition of k T broadening smears out this behavior. particularly at lower p T .
On the other hand, the J/ψ pair distributions are largely unaffected by k T broadening even though the parent b meson pairs are sensitive to the presence of k T broadening. This is because the decay randomizes the direction of the J/ψ relative to the b meson parent, independent of the choice of k T . Thus it is not possible to learn much about broadening in the initial state by studying the final-state decay products. It would be better to look at the B meson pair correlations themselves than studying pair observables through the J/ψ decay products.

VI. THEORETICAL UNCERTAINTIES
Finally, the mass and scale uncertainties on the bb distributions and their transition to the J/ψ pair decay products are discussed here. The results for all pair observables are shown for the lowest and highest minimum p T values, p T > 2 GeV in Fig. 8 and p T > 7 GeV in Fig. 9.
The mass and scale uncertainties are calculated based on results using the one standard deviation uncertainties on the quark mass and scale parameters. The uncertainty band can be obtained for the best fit sets [35,36] 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. Here X is the individual pair observable for a given minimum p T . 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. The central values of the bottom quark mass, µ F /m and µ R /m, as well as their one standard deviation uncertainties, can be found in Sec. III. Note that in the calculation of the uncertainites in the normalized ratios, all distributions are divided by the central value of the total cross section before calculating the uncertainty as in Eqs. (4) and (5). This is consistent with calculating the uncertainty on the distributions via these equations and then dividing by the central value of the integrated cross section and is more consistent with the uncertainties obtained on the distributions themselves [37]. If one instead divided by the total cross section for each mass and scale combination, the uncertainties would be underestimated [37].
The mass and scale variations do not significantly change the shapes of the distributions relative to the  [15] (black for bb, red circles for J/ψ pairs) are also shown.  [15] for pT > 2 GeV. The limits on the uncertainties are shown by dot-dashed curves in both cases. Results are given for the azimuthal difference (a); rapidity difference (b); pair rapidity (c); pT asymmetry (d); pair pT (e); and pair mass (f).
shape of the central distribution. In the case of bottom quark production, the mass uncertainties on the integrated cross section are smaller than those due to the scales by a few percent. The uncertainties on the integrated bb cross section are smaller for the higher p T cuts, decreasing by about a factor of 10 between p T > 2 GeV and 7 GeV.
The J/ψ pair cross sections generally reflect these trends. The J/ψ pair integrated cross sections are smaller and decrase faster with p T cut, resulting in a factor of ≈ 100 decrease between minimum p T of 2 and 7 GeV. This relative difference in cross section is due to the fact, as mentioned previously, that a J/ψ satisfying the same minimum p T originates from a higher p T bottom quark. In addition, due to the decay kinematics, some of the J/ψ's will no longer fall within the accep-tance and J/ψ pairs will not be measured. The mass and scale uncertainty bands are shown for the bb pairs (black dashed curves) and J/ψJ/ψ pairs (red solid curves) and compared to the LHCb data [15] for pT > 7 GeV. The limits on the uncertainties are shown by dot-dashed curves in both cases. Results are given for the azimuthal difference (a); rapidity difference (b); pair rapidity (c); pT asymmetry (d); pair pT (e); and pair mass (f).

VII. RAPIDITY DEPENDENCE
The rapidity dependence of the correlation is studied by calculating the same pair quantities considered by LHCb in the midrapidity region, |y| ≤ 0.8. The results for both the parent bb and the J/ψ pair decay productions are shown in Fig. 10 for p T > 2 and 7 GeV.
Generally, the shapes of the distributions at mid and forward rapidity are rather similar. While there are some differences between the shapes in the two rapidity regions, they are not large. Most of the differences are due to the narrower rapidity acceptance employed at midrapidity, 1.6 units rather than 2.5 units at forward rapidity. The main difference in the azimuthal distributions is the behavior at small |∆φ|. The peak at |∆φ| → 0 is higher at midrapidity, particularly for larger p T b mesons and J/ψs. Thus the narrower rapidity distribution covered at midrapidity seems to favor production topologies where the bb is produced at smaller angles, i.e. gg → bbg where the final-state gluon is hard and balanced against the bb pair. In addition, at |∆φ| → π, the back-toback peak is narrower. The ratios of the distributions at |∆φ| = π to that at |∆φ| = 0 are reduced at midrapidity, particularly for the J/ψ pairs where the results at |∆φ| = π and 0 are ≈ 1 : 1 for both p T values while it is ≈ 2 : 1 for the azimuthal separation between J/ψ pairs at forward rapidity. The differences between the peaks for bb pairs are less pronounced between central and forward rapidities but are still visible.
The rapidity gap distributions, |∆y|, are steeper in the chosen central rapidity region, simply because the ra-pidity range is ≈ 1 unit of rapidity narrower than the forward region studied so far. The pair rapidity distribution, y p , is now symmetric around y p = 0 at midrapidity where there is ample phase space for production. On the other hand, at forward rapidity, the average pair rapidity is not symmetric around the center of the rapidity bin ( y p = 3.25) but closer to the lower end of the range with the calculated y p = 3.08. Regardless of the rapidity region, the |∆y| and y p distributions are independent of p T minimum and whether bb or J/ψJ/ψ pairs are considered.
The average value of the p T asymmetry, A T , is somewhat larger for central rapidity. On the other hand, the average pair p T , p Tp , is higher for midrapidity. This is not unexpected because the average b hadron p T is reduced at forward rapidity relative to central. The average p Tp is ≈ 1 GeV higher at central rapidity for p T > 2 GeV and 3 GeV higher for p T > 7 GeV. Finally, the average pair mass does not vary significantly with rapidity. It increases slightly for forward rapidity, especially for p T > 7 GeV, likely because of the smaller ∆y for the midrapidity interval chosen.

VIII. COLD NUCLEAR MATTER EFFECTS
In Ref. [14], the effects of shadowing and enhanced k T broadening on cc production in cold nuclear matter, with the focus on 5.02 TeV p+Pb collisions, was studied. As discussed there, it has been suggested [38][39][40][41] that energy loss by heavy quarks in heavy-ion collisions could change the azimuthal correlations.
However, it must first be ascertained how the heavy flavor pair distributions are influenced by the presence of cold nuclear matter. For example, additional k T broadening may be present with a nuclear target due to multiple scattering with nucleons along the path of the initial proton (or nucleus). The strength of the effect depends on the impact parameter of the collision. Energy loss in matter, on the other hand, may result in a shift of the transverse momentum of the heavy quark, akin to a change in the fragmentation function. These effects would be in addition to the modification of the parton densities in the nucleus, referred to as shadowing or nPDF effects, calculated assuming collinear factorization.
Here the focus is on the parent bb correlations instead of their decays to J/ψ which do not convey as clear of an effect because the decay is isotropic in the rest frame of the b meson. For illustrative purposes, two particular pair variables are studied out of the six discussed previously: the pair rapidity and the azimuthal separation. Several different scenarios are studied: shadowing alone for both p+Pb and Pb+Pb collisions with k 2 T and ǫ P given in Eqs. (1) and (2); enhanced broadening in p+Pb collisions; and enhanced broadening with energy loss, represented by an increase in the Peterson function parameter, in Pb+Pb collisions. Shadowing is represented by the central EPS09 NLO shadowing ratio [42] for each of the LHCb p T cuts. As noted in Ref. [14] and in Figs. 8 and 9, if k 2 T is kept fixed, the mass and scale uncertainties do not substantially change the shapes of the distributions.
The central gluon modification of the latest nPDF set, EPPS16 [43], is not significantly different from EPS09. However, EPPS16 has five additional parameters relative to EPS09, resulting in larger uncertainty bands with an uncertainty of 25-30% due to shadowing [44]. Although the uncertainty due to shadowing is significant, it is smaller than the uncertainties due to the heavy quark mass and scale variations, particularly for charm quarks [37]. The larger bottom quark mass and comparably larger scales reduce both the overall uncertainty in the baseline p + p cross section and the shadowing effect in p+Pb and Pb+Pb collisions because of the larger parton momentum fraction accessed and the evolution of the shadowing due to the larger factorization scale. To avoid overlapping ratios in the following figures and better illustrate the effects, only results with the central nPDF set are shown.
In Ref. [14], the sensitivity to the magnitude of k T broadening was studied, varying the factor ∆ in Eq. (2) between 0 and 1. So far, in this work, ∆ = 1 has been used as a default. Here, to model broadening in medium, ∆ = 2 is used for p+Pb collisions and ∆ = 4 is used in the Pb+Pb calculations relative to p + p collisions with ∆ = 1. In the case with 'shadowing only', ∆ = 1 is still employed. In addition, energy loss in Pb+Pb collisions is modeled by changing the Peterson function parameter, ǫ P from the value used in these calculations heretofore, ǫ P = 0.0004 [14], to the default value used previously, ǫ P = 0.006 [26]. This change reduces the average z in the Peterson function from 0.93 to 0.83, a difference of about 10%. See Ref. [14] for the sensitivity of the single B meson p T distribution to ǫ P .
The calculations shown here are done at 8.16 TeV for p+Pb collisions and 5 TeV for Pb+Pb collisions. The p + p results used to calculate the nuclear modification factors, R pPb and R PbPb respectively, are calculated at the same energies. The results are calculated both at central (|y| ≤ 0.8) and forward (2 < y < 4.5) rapidity.
Note that there is a rapidity shift in p+Pb collisions due to the requirement of equal velocity beams at the LHC. The calculations shown assume the proton is moving in the positive y direction so that the parton momentum fraction, x, probed by the nucleus, is low. The change in the shadowing ratios is then small. If the beam directions were switched, the momentum fraction in the nucleus is in the antishadowing region. In Pb+Pb collisions, the parton from the forward-going nucleus is large (in the antishadowing region) while that in the backwardgoing direction is small (shadowing) and the collisions are again forward-backward symmetric, as in p + p.
First, the p + p distributions calculated at 5 and 8. 16 TeV were checked to see if the shapes of the pair distributions were modified at different energies. The shapes remain the same for all p T cuts, even at the lowest energy and highest minimum p T . Note that this agreement will eventually break down at lower energies, especially for higher p T , as one reaches the edge of available phase space, particularly at forward rapidity. Because the shapes of the distributions remain unchanged between 5 and 8.16 TeV, these results are not illustrated.
There will be some uncorrelated background to the correlated calculations shown here, particularly in ionion collisions. The background would be larger for cc pairs due to the larger production cross section. Scaling p+p production by the number of binary nucleon-nucleon collisions, several hundred cc pairs can be produced in a single Pb+Pb collision at the LHC [45]. This uncorrelated background would be reduced for bb production because of its relatively smaller production cross section: only a few bb pairs would be produced in a given Pb+Pb event. Even so, the correlated pair signals suggested here could be substantially washed out if they are not seen to be arising from a common decay vertex. Uncorrelated production may also arise in high multiplicity p + p and p+Pb events which could also affect the proposed correlation in these collisions. In lepton pair channels, uncorrelated background could be removed by like-sign pair subtraction [46] although, for bb production, correlated bb pair decays can also lead to like-sign lepton pairs [47].
Aside from independent uncorrelated production, two relatively independent QQ pairs can be produced in double parton scattering in all these collision systems. The double parton scattering contribution has been calculated in Ref. [48] for DD and DD production. The probability of such contributions should be reduced for bottom pair production due to the larger bottom quark mass and higher associated scales.
This section is divided into three subsections. To set the stage, first the single b meson modifications are shown as a function of p T for all four cases (p+Pb with shadowing alone; p+Pb with shadowing and ∆ = 2; Pb+Pb with shadowing only; and Pb+Pb with shadowing, ∆ = 4 and ǫ P = 0.006) at both forward and central rapidity. The pair results are then shown for the pair rapidity and the azimuthal separation between the heavy mesons. Here the nuclear modifications are shown for forward and central p+Pb and Pb+Pb collisions, both with shadowing alone and then including enhanced k T broadening, as well as modification of the fragmentation function in Pb+Pb collisions. However, now the results are shown for the same minimum p T cuts on the b mesons used by LHCb for their bb → J/ψJ/ψ analysis. All results will be presented as the modification factors, R pPb and R PbPb , calculated as the per nucleon cross section in p+Pb and Pb+Pb collisions respectively relative to the p + p result at the same energy.
It has already been noted that there is no modification of the p + p distributions, d ln σ/dX, as a function of center of mass energy. However, some modification of the distributions can be observed in p+Pb and Pb+Pb collisions relative to p + p, as will also be shown. Differences can arise with nuclear beams because the mo-mentum fraction probed changes with changing √ s. The ≈ 40% increase in √ s between 5 and 8.16 TeV reduces the x values correspondingly. Thus the shadowing effect could potentially lead to a modification, especially in Pb+Pb collisions where one of the lead nuclei is probed at relatively high momentum fraction, x ≈ 0.02, in the forward rapidity region. Choosing a higher p T cut also probes higher x and larger scales. Nonetheless, shadowing alone does not modify the shapes of the distributions at different energies. Observable differences only appear with enhanced broadening or modification of the fragmentation function, as is discussed in the remainder of this section.   The results for p+Pb collisions at forward rapidity are compared to LHCb data from nonprompt J/ψ [49] and direct B + mesons [50]. The calculations, using the cen-tral EPS09 NLO set only, with or without any additional k T kick, agree very well with the LHCb data, especially that for nonprompt J/ψs, shown in blue. While the direct B + data, shown in red, are within one standard deviation of the nonprompt J/ψ uncertainty, they are below the central EPS09 NLO calculation for p T < 7 GeV. If the full EPS09 NLO uncertainty band was shown, however, the B + data should be within the limits of the calculated band. The ratios including a higher k T kick for p+Pb collisions are very similar to those with shadowing only, similar to the midrapidity calculations shown for single D mesons at midrapidity in Ref. [14]. The maximum shadowing effect on b mesons at low p T is ≈ 15% at forward rapidity and ≈ 10% at midrapidity.
The single bottom R pPb (p T ) at forward, backward and mid-rapidity was also calculated for shadowing only at 8 TeV [51] employing a p T and rapidity dependent parameterization of the p + p cross section. Bands were shown for several nuclear parton densities and compared to calculations with k T broadening and energy loss [52] in Ref. [44]. The shadowing parameterizations employed in Ref. [44] ranged from ≈ 5 − 40% effects at p T ≈ 5 GeV at forward rapidity with a slightly weaker effect at central rapidity. On the other hand, for the same p T , a 0 − 8% enhancement was seen for the calculations with broadening and energy loss [44].
While not shown, the calculated ratio at backward rapidity is subject to 5−20% antishadowing at low p T . This is consistent with the p T -integrated modification factor R pPb shown as a function of rapidity in Ref. [44]: antishadowing at backward rapidity and increasing shadowing at central and forward rapidity. The LHCb B meson data at forward and backward rapidity are consistent with this trend [49,50].
The NLO results calculated for p+Pb collisions shown here with the central EPS09 NLO set at forward and central rapidity are in good agreement with the bands shown for the parameterizations with EPS09 NLO shadowing applied to the parameterization of the cross section in Refs. [44,51].
The results shown in Fig. 11(c) and (d) for Pb+Pb collisions are, perhaps, somewhat more surprising. First, for shadowing only, it is notable that, at forward rapidity, R PbPb (p T ) is approximately unity with a negligible p T dependence while the modification factor is a rather strong function of p T at midrapidity. These results are easily explained, however. The Pb+Pb modification factor at forward rapidity is the product of the p+Pb results at forward and backward rapidity, R AA (y; X) = R pA (y; X) × R pA (−y; X) [53] where X is another kinematic quantity such as p T or |∆φ|.
The combination of ≈ 15% suppression at forward rapidity with a ≈ 15% enhancement at backward rapidity is effectively unity. (Because Pb+Pb collisions are symmetric about midrapidity, one would see a similar modification factor for Pb+Pb collisions at backward rapidity.) On the other hand, the midrapidity shadowing results are effectively the R pPb (p T ) result in Fig. 11(b) squared.
Note that the analogy between R PbPb and the product of R pPb at forward and backward rapidity here is not exact because of the different energies of the two collisions: the lower energy Pb+Pb collisions would be at higher x than the corresponding rapidity in p+Pb collisions. However, the difference between the two should be small.
The Pb+Pb calculations shown in the dashed histograms, including increased broadening and a modified fragmentation function parameter, exhibit quite different behavior. As shown in Figs. 11(a) and (b), increasing the relative k T broadening in the lead nucleus does not strongly change the modification factor. Therefore the change in slope seen in these results is due to the change in the fragmentation parameter ǫ P . Increasing ǫ P changes the slope of the p T distribution, enhancing the low p T part of the spectrum and depleting the high p T contribution, see the b meson p T distributions in Ref. [14]. Thus, this behavior, while perhaps initially surprising, is easily understood.
In the following subsections, the pair results will be presented. In these, the p T cuts used by LHCb are applied. One must keep in mind that pair quantities are all integrated over p T from the minimum value and thus probe, on average, higher p T and, consequently, somewhat larger x than the single meson quantities shown here.

B. Modifications of yP
The modifications of the pair rapidity are shown in Figs. 12 and 13 for forward and central rapidities respectively. In p+Pb collisions at 8.16 TeV, for an average pair rapidity of 3 in the LHCb acceptance, x ≈ 10 −4 for b quark production for the minimum p T of 2 GeV. The minimum x remains of this order for all values of the minimum p T . Thus R pPb < 1 for all minimum p T values.
Given the x range, it is not surprising that R pPb is nearly independent of y p since the EPS09 NLO gluon nPDF ratio is approximately flat for x < 0.001 [42]. The factorization scale is also important for the nPDF ratio because the QCD scale evolution reduces the shadowing effect at higher p T as well. In addition, the average pair mass, which should be considered when calculating x instead of the transverse mass of a single b quark, is ≈ 15 GeV for p T > 2 GeV and ≈ 23 GeV for p T > 7 GeV. Thus one sees a mild tendency for R pPb (y p ) to increase slightly as the minimum p T increases, an effect more visible at central rapidity since the b meson R pPb (p T ) rises faster with p T at central than at forward rapidity.
When the average k T broadening is effectively doubled, as in Figs. 12(b) and 13(b), the ratios are still relatively independent of y p but the values of R pPb increase by a few percent relative to calculations with shadowing alone. A large effect is not expected because, even for a doubling of the k T kick, k 2 T = 5 GeV 2 in this case, m T is still larger than k 2 T 1/2 and, as seen in Ref. [14], changing the k T kick does not have a large effect on the shape of the single b meson p T distributions. The change in R pPb (p T ) is also minimal, see Fig. 11. A much larger effect was seen on the c and b quark p T distributions by modifying the fragmentation function [14]. Figure 12(c) shows the effect of shadowing alone on results at forward rapidity from Pb+Pb collisions at 5 TeV. As noted earlier, the x range probed is slightly higher even though the factorization scale remains the same, due to the lower energy. Now, however, the ratio R PbPb is no longer independent of y p but shows some structure due to the combination of nuclear effects from both nuclei. This is because one of the lead nuclei is now also probing the nPDFs at higher x, x ≥ 0.01, and moves through the antishadowing region as y p increases.
The increase in curvature with minimum p T is primarily due to the antishadowing contribution from backward rapidity. As p T and thus the input scale of the nPDF increases, the antishadowing peak both reduces its maximum and moves closer to midrapidity. At the same time, the modification in the shadowing region at forward rapidity is reduced but remains relatively independent of rapidity. Thus the results for shadowing only shows more change with y p for p T > 7 GeV than 2 GeV.
Perhaps the most intruiging result is seen in Fig. 12(d) where the average k T kick is again doubled over that employed in p+Pb collisions, to k 2 T ≈ 8.4 GeV 2 . (There are, of course, some small variations in k 2 T between 5 and 8.16 TeV due to the energy dependence assumed for k 2 T , given in Eq. (2).) If that was the only effect assumed for Pb+Pb collisions, one would have expected R PbPb to be similar to the p+Pb resultsin Fig. 12(b).
However, here the ratio with the highest minimum p T is now lowest at the largest y p . This is because, in addition to the k T broadening, an effective energy loss has been introduced by changing the Peterson fragmentation function parameter from the value determined in Ref. [14], to agree with the FONLL b meson p T distribution, to the e + e − default value, ǫ P = 0.006. As stated in Sec. VIII A, this is an effective reduction in the average fraction of momentum transfered from the quark to the meson, from ≈ 93% with ǫ P = 0.0004 to ≈ 83% for ǫ P = 0.006 [14]. As seen in Fig. 11(c), integration starting from p T > 2 GeV includes the peak of the shifted p T distribution where there is an enhancement while p T > 7 GeV includes a region of relative suppression compared to p+ p, resulting in stronger modification for the higher p T cut than the lower. This is an inversion of normally expected behavior for heavy flavor R PbPb . Note that this is in no way intended to replace a real energy loss calculation but is rather intended to illustrate the possible effect on R PbPb for correlated observables.
Results as a function of y p in the central rapidity region are shown in Fig. 13. The trends are quite similar for p+Pb collisions although the level of shadowing is reduced in both cases and a slightly larger separation of the results for different minimum p T values can be seen. However, for Pb+Pb collisions, the results are now also independent of y p . This is because, around the narrow midrapidity window, the x values probed do not change significantly and whatever change occurs is probed symmetrically around y p = 0. The stronger modification for higher p T with the increase of ǫ P is still evident here, albeit with less separation between results for different minimum p T . The y p distributions are shown for p + p collisions at √ s = 7 TeV (blue); p+Pb at √ s N N = 8.16 TeV (red); and Pb+Pb collisions at √ s N N = 5 TeV (black) for p T > 2 and 7 GeV in Fig. 14, including all cold matter effects. As one might expect from the discussion in Sec. IV, the ∆y and y p distributions are unaffected by k T broadening. Thus the y p distributions for p+p and p+Pb collisions are the same shape. However, the Pb+Pb distribution is clearly shifted backward to lower y p , a small but visible effect that increases with minimum p T . This is due to the change in ǫ P . The effect is stronger for p T > 7 GeV than p T > 2 GeV since the larger lower limit of p T integration is more sensitive to the fragmentation function. The steeper p T distributions with the higher value of ǫ P mean that fewer b quarks may be found at higher rapidity, especially for higher values of the minimum p T , reducing the average y p in Pb+Pb collisions relative to the other cases.  Figures 15 and 16 show R pPb and R PbPb as a function of the azimuthal separation between the b and b, at forward and central rapidities respectively. Note that for shadowing only at forward rapidity, the modification factors are rather independent of |∆φ|, with a mild decrease in the ratio as |∆φ| → π. At central rapidity, the shadowing only results show a somewhat stronger decrease in the modification factor as |∆φ| increases. A similar result was obtained for cc production at 5 TeV in Ref. [14] where the p T -and rapidity-integrated R pPb shadowing ratios were independent of |∆φ| with k 2 T = 0 but showed a slight decrease with increasing |∆φ| with k 2 T = 0. Note that in Figs. 15(a) and 16(a), (c), the modification factor decreases rather gradually with |∆φ| over most of |∆φ| with an increase in the slope as |∆φ| approaches π. (All the ratios are compatible with unity for Pb+Pb collisions with shadowing alone in Fig. 15(c), as might be expected from the result in Fig. 11(c).) When the minimum p T is increased, the ratios are independent of |∆φ| until they begin to decrease at larger |∆φ|. This can be attributed to the narrowing and sharpening of the peak in the |∆φ| distribution with increasing p T , seen in Fig. 1, while the enhancement at |∆φ| → 0 is increasing more slowly. The more striking effect is for p+Pb and, in particular, Pb+Pb collisions with enhanced k T broadening. As shown in Fig. 17, the effect of broadening on the azimuthal distributions in p+Pb and Pb+Pb collisions reduces and broadens the peak at |∆φ| ≈ π and enhances the distribution at |∆φ| ≈ 0. Recall that these distributions are the same shape in p + p collisions so that the differences seen in the figure arise primarily from enhanced broadening. Only the results for the lowest and highest minimum p T values are again shown to illustrate the effect.

C. Modifications of |∆φ|
There is an interesting change of behavior at |∆φ| ≈ 0 in p+Pb relative to Pb+Pb collisions for the two different p T cuts. At lower p T , where a change in broadening has a larger effect on the shape of the |∆φ| distribution: the Pb+Pb result is slightly enhanced over that of p+Pb since ∆ = 2 for p+Pb and 4 for Pb+Pb [14]. On the other hand, for the higher p T cut, the enhancement due to broadening is reduced and the change in the fragmentation function parameter suppresses the |∆φ| enhancement at |∆φ| ≈ 0 in Pb+Pb relative to p+Pb, even though the k T broadening is larger in Pb+Pb collisions, see Ref. [14].
The p+Pb ratios with enhanced k T broadening in both rapidity regions exhibit a kink that occurs at higher ∆φ for increasing minimum p T . This can be understood from the ratios of increasing k 2 T relative to the results with no broadening, k 2 T = 0. Reference [14] studied the turn on of the effect at k 2 T > 0, becoming increasingly isotropic as k 2 T increases. As shown in Ref. [14], the |∆φ| distributions peak more sharply at both |∆φ| → π and |∆φ| → 0. The effect at |∆φ| = 0 is reduced in bb production relative to cc since it requires a much harder gluon to balance a more massive bb pair than the lighter cc pair. This change in relative height of the peak for fixed k 2 T and increasing minimum p T causes the location of the kink in the ratio to increase from |∆φ| ≈ 1.7 to 2.5 radians as the minimum p T increases from 2 to 7 GeV.
The heirarchy is more clearly reversed for Pb+Pb collisions, shown in Figs. 15(d) and 16(d). The fragmentation function parameter ǫ P has almost no effect on the shape of the ∆φ distribution, as also shown in Ref. [14] when integrated over all p T . However, it will change the number of bb pairs with both quarks in the rapidity acceptance, as illustrated in Fig. 14, producing the inverted hierarchy of ratios seen here. Note that the larger k T kick assumed for Pb+Pb collisions also result in the kink in R pPb seen in Figs. 15(b) and 16(b), moving to lower ∆φ, now between 1.5 to 2.4 radians in Figs. 15(d) and 16(d).

IX. SUMMARY
The bb → J/ψJ/ψ pair observables measured by LHCb in p + p collisions were studied in detail in an exclusive NLO calculation with fragmentation and k T broadening, as first described in Ref. [14]. The calculations reproduced the data very well in all cases and for all p T cuts. The sensitivity of the results to the k T broadening is shown and, while the direct bb observables are indeed sensitive to the k T broadening, the resulting J/ψ pairs are not since the decays produce a decorrelation of the J/ψs relative to the parent b hadrons. The mass and scale dependence has also been studied and shown not to be large, as expected for bb production. The dependence of the results on rapidity were also shown.
Finally, the nuclear modification factors for enhanced k T broadening and fragmentation function modification in cold nuclear matter were presented. The potential cold nuclear matter effects calculated here for p+Pb and Pb+Pb collisions are not intended to be definitive but illustrative only. The calculations have demonstrated how effects like broadening and energy loss could be disentangled by specific correlated observables more sensitive to each. Although both observables discussed are affected by the two effects, the pair rapidity is more senstive to fragmentation while the azimuthal correlation depends most strongly on the k T broadening. While the effects were modeled in the context of cold nuclear matter, similar decorrelation, as produced by enhanced k T broadening, could be due to hot matter effects, as produced in the quark-gluon plasma [38]. A thermal medium also results in heavy quark energy loss, as modeled by the modified ǫ P . These calculations thus suggest that additional correlated observables are required to better quantify such effects, regardless of the medium.