The Collectivity of Heavy Mesons in Proton-Nucleus Collisions

Using a model based on the Color Glass Condensate framework and the dilute-dense factorization, we systematically study the azimuthal angular correlations between a heavy flavor meson and a light reference particle in proton-nucleus collisions. The obtained second harmonic coefficients (also known as the elliptic flows) for $J/\psi$ and $D^0$ agree with recent experimental data from the LHC. We also provide predictions for the elliptic flows of $\Upsilon$ and $B$ meson, which can be measured in the near future at the LHC. This work can shed light on the physics origin of the collectivity phenomenon in the collisions of small systems.


I. INTRODUCTION
In the last decade, the collectivity phenomenon in small collisional systems (such as proton-proton and proton-nucleus collisions) has been an extremely interesting topic in heavy-ion physics, with a lot of experimental evidences [1][2][3][4][5][6][7][8] observed in high multiplicity events. This particular phenomenon, which implies non-trivial angular correlations among many produced particles in small systems, can be conveniently described by the Fourier harmonics of the azimuthal angular distribution of the measured particles. If one chooses a reference particle and define ∆φ as the azimuthal angle difference between the measured particle and the reference, one can quantitatively extract these Fourier harmonics coefficients v n ≡ cos n∆φ in high-multiplicity events and find that the azimuthal angular distribution of particles is close to isotropic with a small anisotropy characterized by v n coefficients whose magnitudes are about a few percent.
More remarkably, recent measurements in pP b collisions by ALICE [9] and CMS [10][11][12] showed that heavy mesons such as J/ψ and D 0 have the elliptic flow v 2 comparable to the v 2 values of light hadrons. In the near future, the elliptic flow of heavier mesons such as B mesons and Υ may also be measured with sufficiently abundant statistics at the high-luminosity LHC.
This type of azimuthal angular correlation can have an apparent and intuitive classical interpretation in terms of pressure gradients in relativistic hydrodynamics. In the hydrodynamical approach, it is believed that a small droplet of quark-gluon plasma (QGP) is created even in the high-multiplicity events when two protons or a proton and a heavy nucleus collide. Small spatial anisotropies are generated by initial collisional geometries with small fluctuations in the high-multiplicity events in small collisional systems. The relativistic hydrodynamical evolution, which essentially conserves energy and momentum, is used to describe the space-time evolution of the QGP droplet in the final state after its initial production. Together with other final-state effects including the interaction between the probes and the QGP medium, it can convert the initial spatial anisotropy of the QGP droplet into the anisotropy in the momentum space of the measured particles in the final state. Quantitatively, hydrodynamical approaches [13][14][15][16][17][18][19][20][21][22][23][24] have been very successful in explaining the experimental results and making predictions for the collective behaviors involving light hadrons. Nevertheless, usually due to large masses, heavy flavor particles do not flow as much as light particles. A recent study [25] indicates that the final-state interactions between heavy quarks and the QGP medium are not sufficient to generate enough v 2 to explain the observed elliptic flow for J/ψ and D 0 mesons. Alternative mechanisms and additional sources of anisotropy besides the usual hydrodynamics approach are also proposed in Refs. [26][27][28][29].
In the meantime, the observed collectivity may also have an interesting underlying quantum interpretation from the perspective of the color glass condensate (CGC) framework . Due to the multiple interactions during the initial collision and the quantum evolution of the dense background gluon fields in high energy pP b collisions, the produced particles usually also have small correlations, which can also be interpreted as collective behavior. During the initial-state production process prior to the onset of hydrodynamical evolutions, uncorrelated active partons from the projectile proton strongly interact with the dense gluonic background fields inside the target heavy nucleus and can pick up non-trivial quan-tum correlations described by the so-called quadrupole correlators [61,62] in the CGC framework. In the language of Mueller's dipole model, the quadrupole configuration involved in the correlation arises by conversion from a two-color-dipole configuration due to the socalled inelastic multiple scatterings. Therefore, the angular correlation, which is 1/N 2 c suppressed, can be built up from two initially independent dipoles. In particular, the collective behavior of light hadrons produced in pA collisions can also be quantitatively explained within this framework [59]. Moreover, the CGC framework has been demonstrated to be important in understanding the heavy quarkonium productions [63][64][65][66] in pp and pP b collisions in the low transverse momentum region.
The objective of this paper is to extend our previous work [60] and compute the elliptic flow v 2 for both quarkonia and open heavy mesons. Our result shows that initial state effects due to the strong background gluon fields inside target nucleus can generate sufficient amount of collectivity for heavy quarkonia and open heavy flavor mesons. With reasonable choices of parameters and within the range of validity of our CGC model, we find that one can explain the large v 2 for both J/ψ and D 0 mesons measured by the CMS collaboration. We also make the prediction for the v 2 of Υ and B mesons which can be measured in the near future.
The paper is organized as follows. In Sec.II, we provide detailed derivations for the second anisotropy harmonics of heavy quarkoium, heavy quark or open heavy meson with respect to a light quark in pA collisions in a CGC model. In Sec.III, the numerical results are presented. The conclusion and further discussions are given in Sec.IV.

II. CORRELATIONS IN THE CGC FORMALISM
To study the angular correlation of heavy mesons, let us first build a model for the production of a heavy-quark pair and a reference quark in the CGC formalism in highmultiplicity events of pA collisions. In high-multiplicity pP b collisions, we assume that many active partons from the proton projectile participate in the strong interaction with the target nucleus and get produced in the final state. For the purpose of studying the flow of heavy mesons, we consider the production of an incoming gluon, which splits into a pair of heavy quark and anti-quark, together with an incoming quark, which serves as a reference particle, in the presence of the strong gluonic background field generated by the target nucleus. In our model, we assume that the incoming gluon and the reference quark are initially independent, and therefore they have little correlation in both rapidity and azimuthal an-gle before they interact with the target nucleus. Furthermore, the dominant contribution of this process does not generate any correlation, since the incoming gluon (or the split heavy-quark pair) and the reference particle can interact with the background gluon fields independently, without knowing of the existence of the other. Nevertheless, angular correlations can be built up due to color interferences, if they start to interact with the same color charges in the nucleus simultaneously.
As we show below, in order to obtain non-trivial correlations, we need to evaluate the expectation value of products of dipole amplitudes, up to 1 N 2 c order, in the dense background gluon fields of target nucleus in the language of the CGC framework. Furthermore, we can obtain the correlation between the final state quarkonium and the spectator quark by keeping the total momentum of the heavy-quark pair fixed while integrating over their relative momentum. This is used to calculate the elliptic flow of J/Ψ and Υ in pA collisions. Alternatively, one may integrate over the momentum of either the heavyquark or the anti-heavy-quark, then the collectivity of open heavy flavor mesons can be studied. As the common practice to conserve transverse momentum, we choose to work in the coordinate space which makes the summation of the multiple scattering with the dense nuclear target straightforward in the eikonal limit.
A. Differential rate of the p + A → QQ + q + X process 1. An example of Feynman diagrams that contribute to the differential spectrum of heavy quark pair production accompanied by a spectator quark. The transverse coordinates of partons are illustrated in the diagrams. This diagram shows the interference contribution before and after the splitting of the heavy quark pair. The total contribution also includes the squares of the above individual amplitudes.
In the CGC formalism, the differential spectrum for the production of a heavy quark-antiquark pair plus a spectator quark can be obtained by employing techniques developed in [62,[67][68][69][70][71][72][73] and calculating diagrams as illustrated in Fig. 1. The resulting expression is given by where N is the overall normalization factor and ζ = k + Q /k + g is the longitudinal momentum fraction carried by the heavy quark. The transverse coordinates used above are illustrated in Fig. 1 with r 1 = x g − x g and r 2 = x q − x q . Here dPS = dy Q d 2 k Q⊥ dζd 2 kQ ⊥ dy q d 2 k q⊥ is the differential phase space volume of the final states. To model the double parton distribution from the proton projectile, we use which specifies the momentum and spatial distributions of the initial state gluon-quark pair in the proton, with an overall normalization factor being absorbed into N . The parameter B p controls the transverse size of the proton while ∆ represents the intrinsic transverse momentum distribution of a parton. In addition, f g (x g ) and f q (x q ) are the collinear parton distributions with the parton momentum fractions x g and x q for the incoming gluon and quark, respectively. The expression of the g → QQ splitting function is Before taking the target expectation value, the scattering amplitudes between the incoming partons (or the final state QQ pair) and the background gluon fields in the nuclear target DDD are given by [60] where the above four terms correspond to four possible diagrams illustrated in Fig. 1. In the color dipole language, the multiple scatterings between the measured incoming partons and the target can be represented by the product of three dipole amplitudes. In the McLerran-Venugopalan model [74,75], the color average expectation of each correlator in DDD can be evaluated order by order in the large N c expansion. Up to the 1 N 2 c order, the expectation value of the three-dipole amplitude can be written as follows where Q s is the saturation momentum of the large nucleus and the coordinates of the active partons are given by In order to arrive at the expression in Eq. (4), we have neglected the dipole-size logarithmic dependence of the saturation scale, hence obtained the parametrizations valid in a Golec-Biernat-Wusthoff-like approximation. Let us also comment on the physical interpretation of each term in Eq. (4). The first term, which is the leading N c contribution, comes from the interaction of three independent dipoles (one for Q, one forQ and one for the reference quark) with the dense target gluon fields. This term does not yield any anisotropy since these three dipole amplitudes are uncorrelated in the coordinate space. The other three terms, which are suppressed by 1 N 2 c , arise when two color singlet dipoles out of those three are broken and converted into a color quadrupole during the interaction with the target nucleus. Due to this peculiar type of color interactions, correlations can be generated between the heavy quark pair and the reference quark in the coordinate space. After the Fourier transform, it is then natural for us to obtain anisotropies for heavy mesons in the momentum space of the order 1 B. Differential spectra of heavy quarkonium and open heavy meson together with a reference quark Let us first consider the production of the heavy quarkonium. The differential spectrum at the hadronic level is given by the convolution of the partonic spectrum and the corresponding probability density for the hadron of interest. In our study, the probability distribution of producing a heavy-quarkonium from a QQ-pair is given by the color evaporation model (CEM). To generate an open heavy meson, we use the simple Peterson fragmentation function (FF) [76] for both D 0 -meson and B-meson, or the KKKS FF for D 0 meson [77,78] and KKSS FF for B-meson [79]. In the following, we will use the J/ψ and D 0 mesons as examples. One can easily get the corresponding formulas for Υ or B meson production by replacing the c quark with a b quark.
In the CEM, a cc-pair can form a J/ψ only if its invariant mass satisfies M 2 J/ψ < (k c + kc) 2 < 4M 2 D and the corresponding probability is denoted as F cc→J/ψ . To simplify the multiple dimensional numerical calculation, we take the approximation ζ = 1/2 [80]. The kinematic constraint is then given by where, ∆ k ⊥ = ( k c⊥ − kc ⊥ )/2. In our previous study [60], we analytically integrated ∆ k ⊥ from 0 to ∞ for the sake of simplicity. We have numerically checked that the above kinematic constraint does not play an important role in the final v 2 as shown later in Fig. 2.
The two-particle differential spectrum for the heavy-quarkonium and a reference quark production is where dPS = dy J/ψ d 2 k J/ψ⊥ dy q d 2 k q⊥ is the final state phase space for this process, k min = 0.98 GeV and k max = 2.95 GeV are computed from Eq. (9). In our calculation, we take F cc→J/ψ as a constant which can be absorbed into the new overall normalization factor N J/ψ . The overall normalization will disappear in the calculation of the elliptic flow due to cancellation in ratios. The full phase space integration of d 2 ∆k ⊥ yields a delta function, δ 2 ( r − r ), which allows us to get rid of the d 2 r integral. In this paper, we carry out a more sophisticated calculation by using Eq. (10).
As to the case of the open heavy flavor meson, such as the D meson, the differential spectrum for D 0 + q production process can be obtained from Eq. (1) by integrating over the phase space ofc quark and convoluting with the c → D 0 fragmentation function, D(z). The full space integral of d 2 kc ⊥ generates a delta function, (2π) 2 δ 2 ( r 1 − ζ( r − r )), which removes the d 2 r integral. The final expression is then given by where the phase space for this process is dPS = dy D d 2 k D⊥ dy q d 2 k q⊥ and k c⊥ = k D⊥ /z.

C. Elliptic flow of heavy-quarkonium and open heavy-meson
Following the experimental setup and the usual convention, we define the transverse momentum dependent n-th Fourier harmonic from the differential spectrum of particle X plus the reference quark production as [81] Using the commonly used two-particle correlation method adopted by the CMS collaboration[10] for heavy meson flows, the elliptic flow can be computed as follows where In the following, we will present the expression for the second and zeroth harmonics for the cases of J/ψ + q and D 0 + q productions. Then, the corresponding elliptic flow can be easily obtained from Eq. (13). Substituting Eq. (10) into Eq. (12), we can obtain the second and zeroth harmonics that are required to calculate the elliptic flow of heavy quarkonium. Since there are many dimensions of integrations, in order to obtain elliptic flows numerically, our strategy is to analytically perform as many integrations as possible and evaluate the rest of the integrations numerically. The second harmonic is given by where J i 's are Bessel functions of the first kind and F i 's are defined as Since the dy q and d 2 k q⊥ integrals approximately factorize, one writes dy q x q f q (x q ) = dx q f (x q ) which is just total number of the reference quark. To obtain the corresponding v 2 , we have to rely on the numerical evaluation of the rest of the 10-dimension integral. The zeroth Fourier harmonic on the other hand is relatively simple. The d 2 r 2 integral is eliminated by the delta function obtained from the d 2 k q⊥ integration. The integrals over d 2 b 1 d 2 b 2 d 2 r 1 can be carried out analytically. For self-consistency in terms of the large N c expansion, we only need to keep the leading-N c contributions as well. In the end, we only need to numerically compute a three-dimentional integral, which is given by It is easy to note that the integral dy q x q f q (x q ) gives an overall factor which cancels the identical one in κ 2 , together with the overall constant. For the open heavy meson production accompanied by a reference quark process, we can simplify the expression by using several tricks which are provided in Appendix A. These procedures yield the expression which is less time consuming and more accurate in the numerical evaluation. The second harmonic of the differential spectrum of the open heavy meson plus a reference quark production is given by where the F D i 's are defined as with a q = [2∆ 2 + (2 − η)Q 2 s ]/8, φ A the angle between r 1 and r 1 + (1 − ζ) r, φ B the angle between r 1 and r 1 − ζ r and φ C the angle between r 1 and r 1 /ζ − (1 − ζ) r. In fact, F D 3 (r 1 , r) and F D 4 (r 1 , r) give exactly the same contribution.
This symmetry also indicates that the elliptic flow of thec quark orD 0 meson are exactly the same with that of the c quark or D 0 meson. Employing the same tricks presented in Appendix A, we can carry out the d 2 b 1 d 2 b 2 integration analytically for the zeroth harmonic contribution. It is then given by The above expressions allow us to compute the second harmonics for both heavy quarkonia and open heavy meson similar to the measurement conducted at the LHC.

D. Matching to the cross section of single inclusive particle production
In addition, if one integrates over the phase space of the reference quark, it is natural to find that the zeroth harmonic of two-particle spectrum becomes the differential spectrum of single inclusive particle production. Therefore, we should expect that dκ 0 /dydk ⊥ matches to the single inclusive particle cross section, if the overall normalization factor is properly recovered.
The differential cross section for single inclusive J/ψ production can be found in Refs. [63][64][65]. As is pointed out in those works, the CGC framework can only describe the low transverse momentum spectrum of the heavy quarkonium production when p T ≤ Q s , while the high transverse momentum region is described by higher-order QCD calculations in the collinear framework.
In contrast, the situation of the open heavy meson is a bit different: the other heavy quark is unobserved and it can provide sufficient amount of transverse momentum recoils without the need of additional gluon radiation. Usually a simple leading order CGC model calculation with properly chosen saturation momentum can describe the spectrum of open heavy quark or meson from low p T to high p T region. We find that the inclusive D 0 meson cross section in the dilute-dense factorization framework is [61,66,68,82,83] dσ where S ⊥ is the effective area of the target hadron. If we set N = αsT R S ⊥ 2π and remove the reference quark number , which is numerically insignificant.
As shown in Fig. 4, even with the Gaussian parameterization of the scattering amplitude, we can obtain a good description of the B meson spectrum measured by CDF by using Eq. (26). This may imply that we could push our v 2 results for open heavy flavor to a regime of higher p T .

III. NUMERICAL RESULTS
In this section, we provide the results of the numerical evaluation of the heavy meson v 2 derived in previous sections. Together with the input of parton distribution functions for the incoming proton [84], we can describe the collective behavior of heavy quarkonia and open heavy flavors with the same CGC model, which indicates the robustness of the anisotropy generated from the initial state effects in the CGC formalism in small collision systems.

A. Elliptic flow of heavy quarkonia
In Ref. [60], we have presented the numerical results of elliptic flow of heavy quarkonia by analytically integrating over the relative momentum between the heavy quark pair QQ from 0 to ∞. This is an approximation which does not affect the resulting v 2 as we show in the following numerical calculations.
In this work, we perform a more sophisticated calculation with the proper kinematic constraints implied in the CEM using Eqs. (14) and (19). As suggested in Ref. [60], we indeed find little difference between the two calculations as shown in Fig. 2. This indicates that the kinematic constraints, which reduce the yield of heavy quarkonia, have little impact on the angular correlation such as v 2 . Similarly to the results in Ref. [60], a very weak mass dependence of the heavy quarkonium v 2 is found in the numerical evaluation, mainly due to the fact that the mass dependent terms in κ 2 and κ 0 cancel each other in the leading small r 2 expansions where r is the coordinate separation of the QQ pair.  To compare with the experimental data which measures the correlation between J/Ψ and a charged hadron which serves as the reference particle, we need to take into account both g → QQ + q and g → QQ + g channels since the charged hadron can be fragmented either from a quark or a gluon. When the reference quark is replaced by a reference gluon, we need to compute a slightly different set of diagrams. The elliptic flow of heavy quarkonia from the g → QQ + g channel has also been studied and the numerical results are shown in Fig. 3. Since the difference between these two channels are negligible, we do not expect a significant modification on the final results. In addition, we have also found that the use of the charged-hadron fragmentation function [85] does not affect the final results for heavy quarkonium flow, since we also need to integrate over the phase space of the charged-hadron as indicated in the measurement.
3. Elliptic flows of heavy quarkonia with a quark or a gluon as the reference particle.
As pointed out in Ref. [60], our above calculation for heavy quarkonia is valid only in the low transverse momentum region for several reasons. At leading order in the CGC framework, the transverse momentum of final state heavy quarkonium receives transverse momentum contribution only from the transverse momenta of the incoming partons. In the large p T region, where a turnover in v 2 as a function of the transverse momentum should occur, our simple parametrization of the Golec-Biernat and Wusthoff type Gaussian distribution [86] becomes insufficient to describe the heavy quarkonia p T spectrum. Instead, we would need to employ a more accurate and sophisticated implementation of small-x dipole amplitudes, such as the numerical solution to the small-x evolution equations. Furthermore, as shown in Ref. [63,64], one needs to take the extra gluon radiation into account in order to generate large momentum recoils in the high p T region. Interestingly, the situation changes for the open heavy meson production as we discuss below.

B. Elliptic flow of open heavy mesons
For the case of the open heavy flavor meson production, we only measure one heavy quark and integrate over the phase space of the other one. The transverse momentum distribution of the final state open heavy meson at large transverse momentum is mostly controlled by the hard g → QQ splitting, which can provide sufficiently large momentum recoils. Although the total transverse momentum of the QQ pair is predominantly small, each individual heavy quark can have a large transverse momentum due to the hard g → QQ splitting. The situation is similar to the inclusive hadron or jet production in the collinear factorization framework, in which the leading order calculation can already provide a good description of the transverse momentum distribution. For example, the transverse momentum spectrum of B mesons in pp collision can be computed and the results are shown in Fig. 4. Within our simplified CGC model, the shape of the experimental data reported in Ref. [87] in both low and high p T regime can be described with a corresponding saturation momentum for the proton target. 1 We expect that a better agreement could be reached if one uses the numerical solution to the small-x evolution equation instead of the Golec-Biernat and Wusthoff type Gaussian distribution. This implies that we may be able to extend the region of validity of our calculation to high transverse momentum for the open heavy flavor meson, as we show below.
Using the same set of parameters as in the calculation for heavy quarkonia, we numerically compute the v 2 of open heavy mesons up to 8 GeV in g + q channel using Eqs. (20) and (25) and show the results in Figs. 5 and 6. In our numerical evaluation, we have adopted the FFs provided by the Peterson model [76] for both D 0 -meson and B-meson, and also the KKKS FF for D 0 meson [77,78] and the KKSS FF for B-meson [79]. There is a clear shift from the v 2 of c quark to that of the D 0 meson while it is less obvious for the v 2 of the b quark and the B meson. This is mainly due to the b → B fragmentation function, which is strongly peaked at a larger value of z compared to the c → D 0 one. As shown in Fig. 6, we 1 Strictly speaking, our model is more applicable to pA collisions.
Nevertheless, the comparison shown in Fig. 4 serves as a quantitative example in order to demonstrate our discussion regarding the high transverse momentum region.  [11,12]. In the numerical calculation, the same set of values for parameters have been adopted than in the quarkonia case, namely, Bp = 6 GeV −2 , ∆ = 0.5 GeV and Q 2 s = 5 GeV 2 .
find that the resulting v 2 is insensitive to the choices of the FFs.
As shown in Fig. 5, our calculation of the v 2 for D 0 meson production can describe the CMS data [11] reasonably well within the uncertainties of the experimental data. In addition, we can make prediction for the second harmonic coefficients of b quarks and B mesons as well, which are strongly suppressed as compared to those coefficients of D mesons and heavy quarkonia. In contrast to the heavy quarkonia case, we observed a strong mass dependence for the open heavy meson v 2 in our theoretical and numerical calculations. For the v 2 of heavy quarkonia, the mass dependent term coming from the splitting function is only present in the d 2 r and d 2 r integrals which can be factorized out from the integration of azimuthal angle of the reference quark. Therefore, it contributes little to the elliptic flow. Physically speaking, since we always require the heavy QQ pair to be close together in order to produce the quarkonium, the splitting process does not really modify the correlation between the QQ pair and the reference quark. For the Fourier harmonics of open heavy mesons, as can be seen from Eqs. (23)(24), the d 2 r, d 2 r and d 2 r 2 integrals are entangled together. In this case, we are studying the correlation between one final state heavy quark out of the splitting and a reference parton. Since the distance between the QQ pair can be arbitrarily large, the mass dependence naturally comes in the correlation. In addition, we know that usually the mass of heavy quarks always contributes as a suppression in the propagator, and we also note the scale ordering m c < Q s < m b . Therefore, it is reasonable to expect that the D meson can have a sizable v 2 coefficient, while the correlation between the B meson and the reference particle should be suppressed.
Although our model is only applicable to pA collisions, the feature of the heavy flavor v 2 shown in Fig. 5 is also qualitatively in agreement with the recent ATLAS measurement on the elliptic flow of muons from the decay of charm and bottom hadrons measured in pp collisions [88], which indicates that the bottom flow is suppressed as compared to the charm flow.  At last, we have also numerically evaluated the elliptic flow of open heavy mesons in the g → QQ + g channel with an incoming gluon from the proton projectile as the reference. As shown in Fig. 7, we find again that there is little numerical difference between these two channels, which means that the elliptic flow coefficients are insensitive to the type of reference particle. This conclusion is the same as that for the ellipitic flow of heavy quarko-nium. The direct comparison with the experimental data on the correlation should involve a proper combination of two channels. However, we expect no significant deviations from the contribution of each individual channel. We have also tested that our numerical result does not depend strongly on the fragmentation function of the reference parton, again due to the cancellation in the v 2 calculation.

IV. CONCLUSION
In summary, we have studied the azimuthal angle correlation, and derived the analytic expressions of the second Fourier harmonic coefficients, for heavy quarkonium, heavy quarks or heavy mesons with respect to the reference quark, in the dilute-dense factorization in pA collisions. Our calculations of the elliptic flow of heavy mesons (J/ψ and D 0 mesons) in the CGC formalism are consistent with the data from the CMS collaboration. In addition, we made predictions for the elliptic flow of Υ and B mesons, which could be measured in the future. As explained above, we predict that the v 2 of B mesons should be significantly suppressed as compared to that of the D meson due to the mass dependence in the open heavy flavor channels, although we find little mass dependence in the heavy quarkonium channels.
To explore robustness of the collectivity of heavy mesons due to initial state effects in the CGC formalism, we computed the corresponding v 2 in several slightly different setups of approximations and model inputs. First, we computed the heavy quarkonium v 2 by integrating over the relative momentum of the heavy quark pair QQ analytically from 0 to infinity as an approximation in Ref. [60]. In the meantime, we can also calculate the correlation by numerically integrating the relative momentum of the heavy quark pair QQ within the kinematical range according to the CEM. These two calculations yield similar numerical results for the correlations. Furthermore, we find that the resulting correlation remains almost the same if we use a gluon or a charged hadron as the reference particle instead of a quark as we originally chose. In addition, little dependence on various types of heavy meson fragmentation functions was found in our numerical calculations.
Future experimental and theoretical efforts along this line may be able to help us explore the origin of collectivity in high-multiplicity events in high-energy pA collisions. The collectivity of heavy mesons could also provide us an interesting gateway to understand the properties of the initial-state dense gluon matter inside high-energy hadrons.
Note added: During the completion of this manuscript, we have made the prediction for the v 2 of the B meson available to the CMS collaboration. After convoluting with the decay kinematics for the B → D 0 decay by using PYTHIA [89], CMS collaboration finds our calculation with only initial state effects to be consistent with their recent non-prompt D meson v 2 data[12]. We have employed Jaxodraw [90,91] to draw the Feynman diagrams in this paper. + exp − Q 2 s 4 ( r 2 1 ζ 2 + r 2 2 + (ζ 2 + (1 − ζ) 2 )r 2 − 2(1 − ζ) ζ r 1 · r) The b 1,2 dependece in Eqs. (A3-A3) always takes the form exp[− ηQ 2 , where X could be any two-dimensional vector.