The Collectivity of Heavy Mesons in Proton-Nucleus Collisions
Cheng Zhang, Cyrille Marquet, Guang-You Qin, Yu Shi, Lei Wang, Shu-Yi Wei, Bo-Wen Xiao
TThe Collectivity of Heavy Mesons in Proton-Nucleus Collisions
Cheng Zhang,
1, 2
Cyrille Marquet, ∗ Guang-You Qin,
1, 4, † Yu Shi, Lei Wang, Shu-Yi Wei,
5, 3, ‡ and Bo-Wen Xiao
1, 3, § Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Key Laboratory of Particle Physics and Particle Irradiation (MOE),Institute of frontier and interdisciplinary science,Shandong University, QingDao, Shandong 266237, China CPHT, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91128 Palaiseau, France Nuclear Science Division Mailstop 70R0319, Lawrence Berkeley National Laboratory, Berkeley, California 94740, USA European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*) and Fondazione Bruno Kessler,Strada delle Tabarelle 286, I-38123 Villazzano (TN), Italy
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 alight reference particle in proton-nucleus collisions. The obtained second harmonic coefficients (alsoknown as the elliptic flows) for
J/ψ and D agree with recent experimental data from the LHC. Wealso provide predictions for the elliptic flows of Υ and B meson, which can be measured in the nearfuture at the LHC. This work can shed light on the physics origin of the collectivity phenomenonin the collisions of small systems. I. INTRODUCTION
In the last decade, the collectivity phenomenon insmall collisional systems (such as proton-proton andproton-nucleus collisions) has been an extremely interest-ing topic in heavy-ion physics, with a lot of experimentalevidences[1–8] observed in high multiplicity events. Thisparticular phenomenon, which implies non-trivial angu-lar correlations among many produced particles in smallsystems, can be conveniently described by the Fourierharmonics of the azimuthal angular distribution of themeasured particles. If one chooses a reference particleand define ∆ φ as the azimuthal angle difference betweenthe measured particle and the reference, one can quan-titatively extract these Fourier harmonics coefficients v n ≡ (cid:104) cos n ∆ φ (cid:105) in high-multiplicity events and find thatthe azimuthal angular distribution of particles is closeto isotropic with a small anisotropy characterized by v n coefficients whose magnitudes are about a few percent.More remarkably, recent measurements in pP b colli-sions by ALICE[9] and CMS[10–12] showed that heavymesons such as J/ψ and D have the elliptic flow v com-parable to the v values of light hadrons. In the near fu-ture, the elliptic flow of heavier mesons such as B mesonsand Υ may also be measured with sufficiently abundantstatistics at the high-luminosity LHC.This type of azimuthal angular correlation can have anapparent and intuitive classical interpretation in termsof pressure gradients in relativistic hydrodynamics. Inthe hydrodynamical approach, it is believed that a smalldroplet of quark-gluon plasma (QGP) is created even in ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] the high-multiplicity events when two protons or a protonand a heavy nucleus collide. Small spatial anisotropiesare generated by initial collisional geometries with smallfluctuations in the high-multiplicity events in small col-lisional systems. The relativistic hydrodynamical evolu-tion, which essentially conserves energy and momentum,is used to describe the space-time evolution of the QGPdroplet in the final state after its initial production. To-gether with other final-state effects including the inter-action between the probes and the QGP medium, it canconvert the initial spatial anisotropy of the QGP dropletinto the anisotropy in the momentum space of the mea-sured particles in the final state. Quantitatively, hydro-dynamical approaches[13–24] have been very successful inexplaining the experimental results and making predic-tions for the collective behaviors involving light hadrons.Nevertheless, usually due to large masses, heavy flavorparticles do not flow as much as light particles. A re-cent study[25] indicates that the final-state interactionsbetween heavy quarks and the QGP medium are not suffi-cient to generate enough v to explain the observed ellip-tic flow for J/ψ and D mesons. Alternative mechanismsand additional sources of anisotropy besides the usual hy-drodynamics approach are also proposed in Refs.[26–29].In the meantime, the observed collectivity may alsohave an interesting underlying quantum interpretationfrom the perspective of the color glass condensate (CGC)framework[30–60]. Due to the multiple interactions dur-ing the initial collision and the quantum evolution of thedense background gluon fields in high energy pP b colli-sions, the produced particles usually also have small cor-relations, which can also be interpreted as collective be-havior. During the initial-state production process priorto the onset of hydrodynamical evolutions, uncorrelatedactive partons from the projectile proton strongly inter-act with the dense gluonic background fields inside thetarget heavy nucleus and can pick up non-trivial quan- a r X i v : . [ h e p - ph ] M a r tum correlations described by the so-called quadrupolecorrelators[61, 62] in the CGC framework. In the lan-guage of Mueller’s dipole model, the quadrupole con-figuration involved in the correlation arises by conver-sion from a two-color-dipole configuration due to the so-called inelastic multiple scatterings. Therefore, the an-gular correlation, which is 1 /N c suppressed, can be builtup from two initially independent dipoles. In particu-lar, the collective behavior of light hadrons produced in pA collisions can also be quantitatively explained withinthis framework[59]. Moreover, the CGC framework hasbeen demonstrated to be important in understanding theheavy quarkonium productions[63–66] in pp and pP b col-lisions in the low transverse momentum region.The objective of this paper is to extend our previ-ous work[60] and compute the elliptic flow v for bothquarkonia and open heavy mesons. Our result shows thatinitial state effects due to the strong background gluonfields inside target nucleus can generate sufficient amountof collectivity for heavy quarkonia and open heavy fla-vor mesons. With reasonable choices of parameters andwithin the range of validity of our CGC model, we findthat one can explain the large v for both J/ψ and D mesons measured by the CMS collaboration. We alsomake the prediction for the v of Υ and B mesons whichcan be measured in the near future.The paper is organized as follows. In Sec.II, we pro-vide detailed derivations for the second anisotropy har-monics of heavy quarkoium, heavy quark or open heavymeson with respect to a light quark in pA collisions ina CGC model. In Sec.III, the numerical results are pre-sented. The conclusion and further discussions are givenin Sec.IV. II. CORRELATIONS IN THE CGCFORMALISM
To study the angular correlation of heavy mesons, letus first build a model for the production of a heavy-quarkpair and a reference quark in the CGC formalism in high-multiplicity events of pA collisions. In high-multiplicity pP b collisions, we assume that many active partons fromthe proton projectile participate in the strong interac-tion with the target nucleus and get produced in the fi-nal state. For the purpose of studying the flow of heavymesons, 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 refer-ence particle, in the presence of the strong gluonic back-ground field generated by the target nucleus. In ourmodel, we assume that the incoming gluon and the refer-ence quark are initially independent, and therefore theyhave little correlation in both rapidity and azimuthal an- gle before they interact with the target nucleus. Fur-thermore, the dominant contribution of this process doesnot generate any correlation, since the incoming gluon(or the split heavy-quark pair) and the reference particlecan interact with the background gluon fields indepen-dently, without knowing of the existence of the other.Nevertheless, angular correlations can be built up due tocolor interferences, if they start to interact with the samecolor charges in the nucleus simultaneously.As we show below, in order to obtain non-trivial cor-relations, we need to evaluate the expectation value ofproducts of dipole amplitudes, up to N c order, in thedense background gluon fields of target nucleus in thelanguage of the CGC framework. Furthermore, we canobtain the correlation between the final state quarkoniumand the spectator quark by keeping the total momentumof the heavy-quark pair fixed while integrating over theirrelative momentum. This is used to calculate the ellipticflow of J/ Ψ and Υ in pA collisions. Alternatively, onemay integrate over the momentum of either the heavy-quark or the anti-heavy-quark, then the collectivity ofopen heavy flavor mesons can be studied. As the commonpractice to conserve transverse momentum, we choose towork in the coordinate space which makes the summationof the multiple scattering with the dense nuclear targetstraightforward in the eikonal limit.
A. Differential rate of the p + A → Q ¯ Q + q + X process ~x Q ~x ¯ Q ~x q ~x g ~x ′Q ~x ′ ¯ Q ~x ′ q ~x ′ g ~r ⊥ ~r ′⊥ FIG. 1. An example of Feynman diagrams that contribute tothe differential spectrum of heavy quark pair production ac-companied by a spectator quark. The transverse coordinatesof partons are illustrated in the diagrams. This diagram showsthe interference contribution before and after the splitting ofthe heavy quark pair. The total contribution also includesthe squares of the above individual amplitudes.
In the CGC formalism, the differential spectrum forthe production of a heavy quark-antiquark pair plus aspectator quark can be obtained by employing techniquesdeveloped in [62, 67–73] and calculating diagrams as il-lustrated in Fig. 1. The resulting expression is given by dNd PS (cid:12)(cid:12)(cid:12)(cid:12) Q ¯ Q q = N (cid:90) d b d b d r d r d rd r (cid:48) (2 π ) e − i(cid:126)k Q⊥ · [ (cid:126)r +(1 − ζ )( (cid:126)r − (cid:126)r (cid:48) )] e − i(cid:126)k ¯ Q⊥ · [ (cid:126)r − ζ ( (cid:126)r − (cid:126)r (cid:48) )] e − i(cid:126)k q ⊥ · (cid:126)r × W ( x g , x q , b , b , r , r ) k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:104) DDD (cid:105) , (1)where N is the overall normalization factor and ζ = k + Q /k + g is the longitudinal momentum fraction carried by the heavyquark. The transverse coordinates used above are illustrated in Fig. 1 with (cid:126)r = (cid:126)x g − (cid:126)x (cid:48) g and (cid:126)r = (cid:126)x q − (cid:126)x (cid:48) q . Here d PS = dy Q d k Q⊥ dζd k ¯ Q⊥ dy q d k q ⊥ is the differential phase space volume of the final states. To model the double partondistribution from the proton projectile, we use W ( x g , x q , b , b , r , r ) = x g f g ( x g ) x q f q ( x q ) π B p exp[ − b + b B p − ∆ ( r + r )4 ]which specifies the momentum and spatial distributions of the initial state gluon-quark pair in the proton, with anoverall normalization factor being absorbed into N . The parameter B p controls the transverse size of the protonwhile ∆ represents the intrinsic transverse momentum distribution of a parton. In addition, f g ( x g ) and f q ( x q ) arethe 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 → Q ¯ Q splitting function is k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) = 8 π m Q [( ζ + (1 − ζ ) ) K ( m Q | (cid:126)r | ) K ( m ¯ Q | (cid:126)r (cid:48) | ) (cid:126)r · (cid:126)r (cid:48) | (cid:126)r || (cid:126)r (cid:48) | + K ( m Q | (cid:126)r | ) K ( m ¯ Q | (cid:126)r (cid:48) | )] . (2)Before taking the target expectation value, the scattering amplitudes between the incoming partons (or the final state Q ¯ Q pair) and the background gluon fields in the nuclear target (cid:104) DDD (cid:105) are given by [60]
DDD = [ D ( (cid:126)x Q , (cid:126)x (cid:48)Q ) D ( (cid:126)x (cid:48) ¯ Q , (cid:126)x ¯ Q ) + D ( (cid:126)x g , (cid:126)x (cid:48) g ) D ( (cid:126)x (cid:48) g , (cid:126)x g ) − D ( (cid:126)x Q , (cid:126)x (cid:48) g ) D ( (cid:126)x (cid:48) g , (cid:126)x ¯ Q ) − D ( (cid:126)x (cid:48) ¯ Q , (cid:126)x g ) D ( (cid:126)x g , (cid:126)x (cid:48)Q )] D ( (cid:126)x q , (cid:126)x (cid:48) q ) , (3)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 ofthree dipole amplitudes. In the McLerran-Venugopalan model [74, 75], the color average expectation of each correlatorin (cid:104) DDD (cid:105) can be evaluated order by order in the large N c expansion. Up to the N c order, the expectation value ofthe three-dipole amplitude can be written as follows (cid:104) D ( (cid:126)x Q , (cid:126)x (cid:48)Q ) D ( (cid:126)x (cid:48) ¯ Q , (cid:126)x ¯ Q ) D ( (cid:126)x q , (cid:126)x (cid:48) q ) (cid:105) = exp (cid:20) − Q s (cid:0) ( (cid:126)x Q − (cid:126)x (cid:48)Q ) + ( (cid:126)x (cid:48) ¯ Q − (cid:126)x ¯ Q ) + ( (cid:126)x q − (cid:126)x (cid:48) q ) (cid:1)(cid:21) (cid:40) (cid:90) dξ (cid:90) ξ dη (cid:34) (cid:18) Q s N c ( (cid:126)x Q − (cid:126)x (cid:48)Q ) · ( (cid:126)x ¯ Q − (cid:126)x (cid:48) ¯ Q ) (cid:19) exp (cid:18) ηQ s (cid:126)x Q − (cid:126)x (cid:48) ¯ Q ) · ( (cid:126)x ¯ Q − (cid:126)x (cid:48)Q ) (cid:19) + (cid:18) Q s N c ( (cid:126)x Q − (cid:126)x (cid:48)Q ) · ( (cid:126)x q − (cid:126)x (cid:48) q ) (cid:19) exp (cid:18) ηQ s (cid:126)x Q − (cid:126)x q ) · ( (cid:126)x (cid:48) q − (cid:126)x (cid:48)Q ) (cid:19) + (cid:18) Q s N c ( (cid:126)x ¯ Q − (cid:126)x (cid:48) ¯ Q ) · ( (cid:126)x q − (cid:126)x (cid:48) q ) (cid:19) exp (cid:18) ηQ s (cid:126)x (cid:48) ¯ Q − (cid:126)x q ) · ( (cid:126)x (cid:48) q − (cid:126)x ¯ Q ) (cid:19) (cid:35)(cid:41) , (4)where Q s is the saturation momentum of the large nucleus and the coordinates of the active partons are given by (cid:126)x Q = (cid:126)b + (cid:126)r − ζ ) (cid:126)r, (cid:126)x ¯ Q = (cid:126)b + (cid:126)r − ζ(cid:126)r, (5) (cid:126)x (cid:48)Q = (cid:126)b − (cid:126)r − ζ ) (cid:126)r (cid:48) , (cid:126)x (cid:48) ¯ Q = (cid:126)b − (cid:126)r − ζ(cid:126)r (cid:48) , (6) (cid:126)x q = (cid:126)b + (cid:126)r , (cid:126)x (cid:48) q = (cid:126)b − (cid:126)r , (7) (cid:126)x g = (cid:126)b + (cid:126)r , (cid:126)x (cid:48) g = (cid:126)b − (cid:126)r . (8)In order to arrive at the expression in Eq. (4), we have neglected the dipole-size logarithmic dependence of thesaturation 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 for ¯ Q and one for the referencequark) with the dense target gluon fields. This term does not yield any anisotropy since these three dipole amplitudesare uncorrelated in the coordinate space. The other three terms, which are suppressed by N c , arise when two colorsinglet dipoles out of those three are broken and converted into a color quadrupole during the interaction with thetarget nucleus. Due to this peculiar type of color interactions, correlations can be generated between the heavy quarkpair and the reference quark in the coordinate space. After the Fourier transform, it is then natural for us to obtainanisotropies for heavy mesons in the momentum space of the order N c . 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 levelis given by the convolution of the partonic spectrum and the corresponding probability density for the hadron ofinterest. In our study, the probability distribution of producing a heavy-quarkonium from a Q ¯ Q -pair is given by thecolor evaporation model (CEM). To generate an open heavy meson, we use the simple Peterson fragmentation function(FF) [76] for both D -meson and B -meson, or the KKKS FF for D meson [77, 78] and KKSS FF for B -meson [79].In the following, we will use the J/ψ and D 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 c ¯ c -pair can form a J/ψ only if its invariant mass satisfies M J/ψ < ( k c + k ¯ c ) < M D and thecorresponding probability is denoted as F c ¯ c → J/ψ . To simplify the multiple dimensional numerical calculation, we takethe approximation ζ = 1 / M J/ψ − M c < (cid:126)k ⊥ < M D − M c , (9)where, ∆ (cid:126)k ⊥ = ( (cid:126)k c ⊥ − (cid:126)k ¯ c ⊥ ) /
2. In our previous study[60], we analytically integrated ∆ (cid:126)k ⊥ from 0 to ∞ for the sake ofsimplicity. We have numerically checked that the above kinematic constraint does not play an important role in thefinal v as shown later in Fig. 2.The two-particle differential spectrum for the heavy-quarkonium and a reference quark production is dNd PS (cid:12)(cid:12)(cid:12)(cid:12) J/ψ + q = (cid:90) d k c ⊥ d k ¯ c ⊥ δ ( (cid:126)k J/ψ ⊥ − (cid:126)k c ⊥ − (cid:126)k ¯ c ⊥ ) θ ( | ∆ (cid:126)k ⊥ | − k min ) θ ( k max − | ∆ (cid:126)k ⊥ | ) (cid:90) dζδ ( ζ −
12 ) dNd PS (cid:12)(cid:12)(cid:12)(cid:12) c ¯ cq F c ¯ c → J/ψ = N J/ψ (cid:90) d b d b d r d r d rd r (cid:48) (2 π ) (cid:90) k max k min d ∆ k ⊥ ∆ k ⊥ J (∆ k ⊥ | (cid:126)r − (cid:126)r (cid:48) | ) e − i(cid:126)k J/ψ ⊥ · (cid:126)r e − i(cid:126)k q ⊥ · (cid:126)r × W ( x g , x q , b , b , r , r ) k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) ζ = (cid:104) DDD (cid:105)| ζ = , (10)where d PS = dy J/ψ d k J/ψ ⊥ dy q d k q ⊥ is the final state phase space for this process, k min = 0 .
98 GeV and k max = 2 . F c ¯ c → J/ψ as a constant which can be absorbed into thenew overall normalization factor N J/ψ . The overall normalization will disappear in the calculation of the elliptic flowdue to cancellation in ratios. The full phase space integration of d ∆ k ⊥ yields a delta function, δ ( (cid:126)r − (cid:126)r (cid:48) ), whichallows us to get rid of the d r (cid:48) 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 + q productionprocess can be obtained from Eq. (1) by integrating over the phase space of ¯ c quark and convoluting with the c → D fragmentation function, D ( z ). The full space integral of d k ¯ c ⊥ generates a delta function, (2 π ) δ ( (cid:126)r − ζ ( (cid:126)r − (cid:126)r (cid:48) )),which removes the d r (cid:48) integral. The final expression is then given by dNd PS (cid:12)(cid:12)(cid:12)(cid:12) D + q = (cid:90) dζd k ¯ c ⊥ (cid:90) dz D ( z ) z dNd PS (cid:12)(cid:12)(cid:12)(cid:12) c ¯ cq = N (cid:90) dζ (cid:90) dz D ( z ) z (cid:90) d b d b d r d r d r (2 π ) ζ e − i (cid:126)kD ⊥ z · (cid:126)r ζ e − i(cid:126)k q ⊥ · (cid:126)r W ( x g , x q , b , b , r , r ) × k + g (cid:88) αβγ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ ) (cid:104) DDD (cid:105)| (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ , (11)where the phase space for this process is d PS = dy D d k D ⊥ dy q d k q ⊥ and (cid:126)k c ⊥ = (cid:126)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 -thFourier harmonic from the differential spectrum of particle X plus the reference quark production as [81] dκ n dy X dk X ⊥ = k X ⊥ (cid:90) dφ X dy q d k q ⊥ e − in ( φ X − φ q ) dNd PS (cid:12)(cid:12)(cid:12)(cid:12) X + q . (12)Using the commonly used two-particle correlation method adopted by the CMS collaboration[10] for heavy mesonflows, the elliptic flow can be computed as follows v ( y X , k X ⊥ ) = dκ /dy X dk X ⊥ dκ /dy X dk X ⊥ v [ref] , (13)where v [ref] is the integrated elliptic flow of the reference quark given by v [ref] = (cid:112) κ [ref] /κ [ref]. In the following,we will present the expression for the second and zeroth harmonics for the cases of J/ψ + q and D + 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 calculatethe elliptic flow of heavy quarkonium. Since there are many dimensions of integrations, in order to obtain ellipticflows numerically, our strategy is to analytically perform as many integrations as possible and evaluate the rest of theintegrations numerically. The second harmonic is given by dκ dy J/ψ dk J/ψ ⊥ = k J/ψ ⊥ N J/ψ (cid:90) dy q k q ⊥ dk q ⊥ (cid:90) r dr d r d rd r (cid:48) (2 π ) (cid:90) dξ (cid:90) ξ dηJ ( | (cid:126)k J/ψ ⊥ || (cid:126)r | ) J ( | (cid:126)k q ⊥ || (cid:126)r | ) cos(2 φ r ) × k max J ( k max | (cid:126)r − (cid:126)r (cid:48) | ) − k min J ( k min | (cid:126)r − (cid:126)r (cid:48) | ) | (cid:126)r − (cid:126)r (cid:48) | x g f g ( x g ) x q f q ( x q ) 14 π exp (cid:20) − ∆ ( r + r )4 (cid:21)
11 + ηQ s B p × k + g (cid:88) αβγ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) ζ = Q s N c (cid:104) F ( r , r , r, r (cid:48) ) + F ( r , r ) − F ( r , r , r ) − F ( r , r , r (cid:48) ) (cid:105) , (14)where J i ’s are Bessel functions of the first kind and F i ’s are defined as F ( r , r , r, r (cid:48) ) = 2 (cid:20) ( (cid:126)r + (cid:126)r − (cid:126)r (cid:48) · (cid:126)r (cid:21) exp (cid:20) − ηQ s ( (cid:126)r + (cid:126)r (cid:48) ) ηQ s B p ) (cid:21) exp (cid:20) − Q s r + r + ( (cid:126)r − (cid:126)r (cid:48) ) (cid:21) × exp (cid:34) ηQ s ( (cid:126)r − (cid:126)r + (cid:126)r − (cid:126)r (cid:48) ) (cid:35) , (15) F ( r , r ) = 2( (cid:126)r · (cid:126)r ) exp (cid:20) − Q s r + r ) (cid:21) exp (cid:20) ηQ s (cid:126)r − (cid:126)r ) (cid:21) , (16) F ( r , r , r ) = 2 (cid:18) ( (cid:126)r + (cid:126)r · (cid:126)r (cid:19) exp (cid:20) − ηQ s r ηQ s B p ) (cid:21) exp (cid:20) − Q s r + r + r (cid:21) × exp (cid:20) ηQ s (cid:126)r − (cid:126)r + (cid:126)r (cid:21) , (17) F ( r , r , r (cid:48) ) = F ( r , r , r (cid:48) ) . (18)Since the dy q and d k q ⊥ integrals approximately factorize, one writes (cid:82) dy q x q f q ( x q ) = (cid:82) dx q f ( x q ) which is just totalnumber of the reference quark. To obtain the corresponding v , we have to rely on the numerical evaluation of therest of the 10-dimension integral.The zeroth Fourier harmonic on the other hand is relatively simple. The d r integral is eliminated by the deltafunction obtained from the d k q ⊥ integration. The integrals over d b d b d r can be carried out analytically. Forself-consistency in terms of the large N c expansion, we only need to keep the leading- N c contributions as well. In theend, we only need to numerically compute a three-dimentional integral, which is given by dκ dy J/ψ dk J/ψ ⊥ = k J/ψ ⊥ N J/ψ (cid:90) dy q rdrd r (cid:48) π k max J ( k max | (cid:126)r − (cid:126)r (cid:48) | ) − k min J ( k min | (cid:126)r − (cid:126)r (cid:48) | ) | (cid:126)r − (cid:126)r (cid:48) | π + 2 Q s × x g f g ( x g ) x q f q ( x q ) k + g (cid:88) αβγ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) ζ = exp[ − k J/ψ ⊥ ∆ + 2 Q s ] × (cid:26) exp (cid:20) − Q s ( (cid:126)r − (cid:126)r (cid:48) ) (cid:21) + 1 − exp (cid:20) − Q s r (cid:21) − exp (cid:20) − Q s r (cid:48) (cid:21)(cid:27) . (19)It is easy to note that the integral (cid:82) dy q x q f q ( x q ) gives an overall factor which cancels the identical one in κ , togetherwith the overall constant.For the open heavy meson production accompanied by a reference quark process, we can simplify the expressionby using several tricks which are provided in Appendix A. These procedures yield the expression which is less timeconsuming and more accurate in the numerical evaluation. The second harmonic of the differential spectrum of theopen heavy meson plus a reference quark production is given by dκ dy D dk D ⊥ = k D ⊥ N (cid:90) dy q k q ⊥ dk q ⊥ (cid:90) dζ (cid:90) r dr d r π ζ (cid:90) dz D ( z ) z (cid:90) dξ (cid:90) ξ dηJ ( | (cid:126)k D ⊥ || (cid:126)r | zζ ) × x g f g ( x g ) x q f q ( x q ) 14 π e − ∆2 r
11 + ηQ s B p k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ × Q s N c [ F D ( r , r ) + F D ( r ) − F D ( r , r ) − F D ( r , r )] , (20)where the F Di ’s are defined as F D ( r , r ) = exp (cid:20) − ηQ s ηQ s B p ) ((1 − ζ ) (cid:126)r − − ζ ζ (cid:126)r ) (cid:21) exp (cid:20) − (2 − η ) Q s r ζ (cid:21) × ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m (cid:18) r ζ (cid:19) m +2 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k q ⊥ a q ) , (21) F D ( r ) = exp (cid:20) − (4 − η ) Q s r (cid:21) ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m r m +21 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k q ⊥ a q ) , (22) F D ( r , r ) = exp (cid:20) − Q s r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s ηQ s B p ) (1 − ζ ) r (cid:21) × exp (cid:20) ηQ s (cid:126)r + (1 − ζ ) (cid:126)r ) (cid:21) cos(2 φ A ) × ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m | (cid:126)r + (1 − ζ ) (cid:126)r | m +2 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k q ⊥ a q )+ exp (cid:20) − Q s r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s ηQ s B p ) ζ r (cid:21) × exp (cid:20) ηQ s (cid:126)r − ζ(cid:126)r ) (cid:21) cos(2 φ B ) × ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m | (cid:126)r − ζ(cid:126)r | m +2 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k q ⊥ a q ) , (23) F D ( r , r ) = exp (cid:104) − Q s (cid:126)r ζ − (1 − ζ ) (cid:126)r ) + ζ r ] (cid:105) exp (cid:20) − ηQ s ηQ s B p ) ( (cid:126)r − ζ(cid:126)r ) (cid:21) exp (cid:20) ηQ s ζ r (cid:21) cos(2 φ r ) × ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m | ζr | m +2 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k s a q )+ exp (cid:104) − Q s (cid:126)r ζ − (1 − ζ ) (cid:126)r ) + ζ r ] (cid:105) exp (cid:20) − ηQ s (1 − ζ ) ηQ s B p ) ( (cid:126)r − (cid:126)r ζ ) (cid:21) × exp (cid:20) ηQ s (cid:126)r ζ − (1 − ζ ) (cid:126)r ) (cid:21) cos(2 φ C ) × ∞ (cid:88) m =0 (cid:20) ηQ s (cid:21) m | (cid:126)r ζ − (1 − ζ ) (cid:126)r | m +2 (2 m + 2)! k q ⊥ m +6 a m +3 q (2 m )! m ! F ( m + 3 , , − k q ⊥ a q ) , (24)with a q = [2∆ + (2 − η ) Q s ] / φ A the angle between (cid:126)r and (cid:126)r + (1 − ζ ) (cid:126)r , φ B the angle between (cid:126)r and (cid:126)r − ζ(cid:126)r and φ C the angle between (cid:126)r and (cid:126)r /ζ − (1 − ζ ) (cid:126)r . In fact, F D ( r , r ) and F D ( r , r ) give exactly the same contribution.This symmetry also indicates that the elliptic flow of the ¯ c quark or ¯ D meson are exactly the same with that of the c quark or D meson.Employing the same tricks presented in Appendix A, we can carry out the d b d b integration analytically for thezeroth harmonic contribution. It is then given by dκ dy D dk D ⊥ = k D ⊥ N (cid:90) dy q (cid:90) r dr d r π (cid:90) dζζ (cid:90) dz D ( z ) z J ( | (cid:126)k D ⊥ || (cid:126)r | zζ ) × x g f g ( x g ) x q f q ( x q ) 14 π e − ∆2 r k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ × (cid:26) exp (cid:20) − Q s r ζ (cid:21) + exp (cid:20) − Q s r (cid:21) − exp (cid:20) − Q s r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) − exp (cid:20) − Q s r ζ + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21)(cid:27) . (25)The above expressions allow us to compute the second harmonics for both heavy quarkonia and open heavy mesonsimilar 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 zerothharmonic of two-particle spectrum becomes the differential spectrum of single inclusive particle production. Therefore,we should expect that dκ /dydk ⊥ matches to the single inclusive particle cross section, if the overall normalizationfactor is properly recovered.The differential cross section for single inclusive J/ψ production can be found in Refs. [63–65]. As is pointed out inthose works, the CGC framework can only describe the low transverse momentum spectrum of the heavy quarkoniumproduction when p T ≤ Q s , while the high transverse momentum region is described by higher-order QCD calculationsin the collinear framework.In contrast, the situation of the open heavy meson is a bit different: the other heavy quark is unobserved and it canprovide sufficient amount of transverse momentum recoils without the need of additional gluon radiation. Usually asimple leading order CGC model calculation with properly chosen saturation momentum can describe the spectrumof open heavy quark or meson from low p T to high p T region. We find that the inclusive D meson cross section inthe dilute-dense factorization framework is [61, 66, 68, 82, 83] dσdy D d k D ⊥ = α s T R S ⊥ (2 π ) (cid:90) d r d r (2 π ) (cid:90) dz D ( z ) z (cid:90) dζζ e − i(cid:126)k D ⊥ · (cid:126)r /zζ x g f g ( x g ) k + g (cid:88) αβλ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) ) (cid:12)(cid:12) (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ × (cid:26) exp (cid:20) − Q s r ζ (cid:21) + exp (cid:20) − Q s r (cid:21) − exp (cid:20) − Q s r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) − exp (cid:20) − Q s r ζ + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21)(cid:27) . (26)where S ⊥ is the effective area of the target hadron. If we set N = α s T R S ⊥ π and remove the reference quark number (cid:82) dy q x q f q ( x q ), Eq. (25) and Eq. (26) only differ from each other in the Gaussian distribution e − ∆2 r , which isnumerically insignificant.As shown in Fig. 4, even with the Gaussian parameterization of the scattering amplitude, we can obtain a gooddescription of the B meson spectrum measured by CDF by using Eq. (26). This may imply that we could push our v results for open heavy flavor to a regime of higher p T . III. NUMERICAL RESULTS
In this section, we provide the results of the numeri-cal evaluation of the heavy meson v derived in previoussections. Together with the input of parton distribution functions for the incoming proton[84], we can describe thecollective behavior of heavy quarkonia and open heavyflavors with the same CGC model, which indicates therobustness of the anisotropy generated from the initialstate effects in the CGC formalism in small collision sys-tems. A. Elliptic flow of heavy quarkonia
In Ref. [60], we have presented the numerical resultsof elliptic flow of heavy quarkonia by analytically inte-grating over the relative momentum between the heavyquark pair Q ¯ Q from 0 to ∞ . This is an approximationwhich does not affect the resulting v as we show in thefollowing numerical calculations.In this work, we perform a more sophisticated calcu-lation with the proper kinematic constraints implied inthe CEM using Eqs. (14) and (19). As suggested in Ref.[60], we indeed find little difference between the two cal-culations as shown in Fig. 2. This indicates that thekinematic constraints, which reduce the yield of heavyquarkonia, have little impact on the angular correlationsuch as v . Similarly to the results in Ref. [60], a veryweak mass dependence of the heavy quarkonium v isfound in the numerical evaluation, mainly due to the factthat the mass dependent terms in κ and κ cancel eachother in the leading small r expansions where r is thecoordinate separation of the Q ¯ Q pair. . . . . . B p = 6 GeV − ∆ = 0 . Q s = 5 GeV k ⊥ (GeV) v ( k ⊥ ) CMS
J/ψJ/ψ , m c = 1 . R k max k min d∆ kJ/ψ , m c = 1 . R ∞ d∆ k Υ, m b = 4 . R k max k min d∆ k Υ, m b = 4 . R ∞ d∆ k FIG. 2. Elliptic flow of heavy quarkonium in g → Q ¯ Q + q channel compared with the experimental data from the CMScollaboration [10]. To compare with the experimental data which mea-sures the correlation between J/ Ψ and a charged hadronwhich serves as the reference particle, we need to takeinto account both g → Q ¯ Q + q and g → Q ¯ Q + g channelssince the charged hadron can be fragmented either froma quark or a gluon. When the reference quark is replacedby a reference gluon, we need to compute a slightly differ-ent set of diagrams. The elliptic flow of heavy quarkoniafrom the g → Q ¯ Q + g channel has also been studied andthe numerical results are shown in Fig. 3. Since thedifference between these two channels are negligible, wedo not expect a significant modification on the final re-sults. In addition, we have also found that the use ofthe charged-hadron fragmentation function [85] does notaffect the final results for heavy quarkonium flow, since we also need to integrate over the phase space of thecharged-hadron as indicated in the measurement. . . . . . B p = 6 GeV − ∆ = 0 . Q s = 5 GeV k ⊥ (GeV) v ( k ⊥ ) CMS
J/ψJ/ψ & g, m c = 1 . J/ψ & q, m c = 1 . m b = 4 . m b = 4 . FIG. 3. Elliptic flows of heavy quarkonia with a quark or agluon as the reference particle.
As pointed out in Ref. [60], our above calculation forheavy quarkonia is valid only in the low transverse mo-mentum region for several reasons. At leading order inthe CGC framework, the transverse momentum of finalstate heavy quarkonium receives transverse momentumcontribution only from the transverse momenta of the in-coming partons. In the large p T region, where a turnoverin v as a function of the transverse momentum shouldoccur, our simple parametrization of the Golec-Biernatand Wusthoff type Gaussian distribution[86] becomes in-sufficient to describe the heavy quarkonia p T spectrum.Instead, we would need to employ a more accurate andsophisticated implementation of small-x dipole ampli-tudes, such as the numerical solution to the small- x evo-lution equations. Furthermore, as shown in Ref. [63, 64],one needs to take the extra gluon radiation into accountin order to generate large momentum recoils in the high p T region. Interestingly, the situation changes for theopen heavy meson production as we discuss below. B. Elliptic flow of open heavy mesons
For the case of the open heavy flavor meson produc-tion, we only measure one heavy quark and integrate overthe phase space of the other one. The transverse momen-tum distribution of the final state open heavy meson atlarge transverse momentum is mostly controlled by thehard g → Q ¯ Q splitting, which can provide sufficientlylarge momentum recoils. Although the total transversemomentum of the Q ¯ Q pair is predominantly small, eachindividual heavy quark can have a large transverse mo-mentum due to the hard g → Q ¯ Q splitting. The situationis similar to the inclusive hadron or jet production in thecollinear factorization framework, in which the leading k ⊥ (GeV) d σ / d k ⊥ ( nb / G e V ) B , KKSS FF, Q s = 0 . , α s = 0 . B , KKSS FF, Q s = 1 GeV , α s = 0 . p ¯ p , √ S = 1 .
96 TeV, | y | < . FIG. 4. Differential cross section for B mesons computed inour model and compared with experimental data from theCDF collaboration [87] in p ¯ p collisions. Note that we usesmaller values of Q s due to the consideration of smaller targetsize and lower collision energy and S ⊥ = 12 mb for a protontarget. order calculation can already provide a good descriptionof the transverse momentum distribution. For example,the transverse momentum spectrum of B mesons in p ¯ p collision can be computed and the results are shown inFig. 4. Within our simplified CGC model, the shapeof the experimental data reported in Ref. [87] in bothlow and high p T regime can be described with a corre-sponding saturation momentum for the proton target. We expect that a better agreement could be reached ifone uses the numerical solution to the small- x evolutionequation instead of the Golec-Biernat and Wusthoff typeGaussian distribution. This implies that we may be ableto extend the region of validity of our calculation to hightransverse momentum for the open heavy flavor meson,as we show below.Using the same set of parameters as in the calculationfor heavy quarkonia, we numerically compute the v ofopen heavy mesons up to 8 GeV in g + q channel usingEqs. (20) and (25) and show the results in Figs. 5 and6. In our numerical evaluation, we have adopted the FFsprovided by the Peterson model [76] for both D -mesonand B -meson, and also the KKKS FF for D meson [77,78] and the KKSS FF for B -meson [79]. There is a clearshift from the v of c quark to that of the D mesonwhile it is less obvious for the v of the b quark and the B meson. This is mainly due to the b → B fragmentationfunction, which is strongly peaked at a larger value of z compared to the c → D one. As shown in Fig. 6, we Strictly speaking, our model is more applicable to pA collisions.Nevertheless, the comparison shown in Fig. 4 serves as a quanti-tative example in order to demonstrate our discussion regardingthe high transverse momentum region. . . . . k ⊥ (GeV) v ( k ⊥ ) CMS, prompt D mesonCMS, D from BD meson, KKKS FF c quark, m c = 1 . ∼ . B mesons, KKSS FF b quark m b = 4 . ∼ . FIG. 5. Elliptic flow of open heavy mesons in the g → Q ¯ Q + q channel compared with the experimental data from the CMScollaboration [11, 12]. In the numerical calculation, the sameset of values for parameters have been adopted than in thequarkonia case, namely, B p = 6 GeV − , ∆ = 0 . Q s = 5 GeV . find that the resulting v is insensitive to the choices ofthe FFs.As shown in Fig. 5, our calculation of the v for D meson production can describe the CMS data [11] rea-sonably well within the uncertainties of the experimentaldata. In addition, we can make prediction for the secondharmonic coefficients of b quarks and B mesons as well,which are strongly suppressed as compared to those co-efficients of D mesons and heavy quarkonia. In contrastto the heavy quarkonia case, we observed a strong massdependence for the open heavy meson v in our theo-retical and numerical calculations. For the v of heavyquarkonia, the mass dependent term coming from thesplitting function is only present in the d r and d r (cid:48) in-tegrals which can be factorized out from the integrationof azimuthal angle of the reference quark. Therefore, itcontributes little to the elliptic flow. Physically speak-ing, since we always require the heavy Q ¯ Q pair to beclose together in order to produce the quarkonium, thesplitting process does not really modify the correlationbetween the Q ¯ Q pair and the reference quark. For theFourier harmonics of open heavy mesons, as can be seenfrom Eqs. (23-24), the d r , d r (cid:48) and d r integrals areentangled together. In this case, we are studying the cor-relation between one final state heavy quark out of thesplitting and a reference parton. Since the distance be-tween the Q ¯ Q pair can be arbitrarily large, the mass de-pendence naturally comes in the correlation. In addition,we know that usually the mass of heavy quarks alwayscontributes as a suppression in the propagator, and wealso note the scale ordering m c < Q s < m b . Therefore,it is reasonable to expect that the D meson can have asizable v coefficient, while the correlation between the B pA collisions,the feature of the heavy flavor v shown in Fig. 5 is alsoqualitatively in agreement with the recent ATLAS mea-surement on the elliptic flow of muons from the decay ofcharm and bottom hadrons measured in pp collisions[88],which indicates that the bottom flow is suppressed ascompared to the charm flow. B p = 6 GeV − ∆ = 0 . Q s = 5 GeV . . . . k ⊥ (GeV) v ( k ⊥ ) CMS, prompt D CMS, D from BD , KKKS D , Peterson B , KKSS B , Peterson FIG. 6. Elliptic flows of open heavy mesons in the g → Q ¯ Q + q channel calculated with different fragmentation functions. . . . . B p = 6 GeV − ∆ = 0 . Q s = 5 GeV k ⊥ (GeV) v ( k ⊥ ) CMS, D meson v sub2 D meson & q, Peterson FF D meson & g, Peterson FF c quark & q, m c = 1 . c quark & g, m c = 1 . FIG. 7. Elliptic flows of open heavy mesons in the g → Q ¯ Q + q and g → Q ¯ Q + g channels. At last, we have also numerically evaluated the ellipticflow of open heavy mesons in the g → Q ¯ Q + g channelwith an incoming gluon from the proton projectile as thereference. As shown in Fig. 7, we find again that thereis little numerical difference between these two channels,which means that the elliptic flow coefficients are insen-sitive to the type of reference particle. This conclusionis the same as that for the ellipitic flow of heavy quarko- nium. The direct comparison with the experimental dataon the correlation should involve a proper combinationof two channels. However, we expect no significant de-viations from the contribution of each individual chan-nel. We have also tested that our numerical result doesnot depend strongly on the fragmentation function of thereference parton, again due to the cancellation in the v calculation. IV. CONCLUSION
In summary, we have studied the azimuthal angle cor-relation, and derived the analytic expressions of the sec-ond Fourier harmonic coefficients, for heavy quarkonium,heavy quarks or heavy mesons with respect to the refer-ence quark, in the dilute-dense factorization in pA col-lisions. Our calculations of the elliptic flow of heavymesons (
J/ψ and D mesons) in the CGC formalism areconsistent with the data from the CMS collaboration. Inaddition, 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 of B mesonsshould be significantly suppressed as compared to thatof the D meson due to the mass dependence in the openheavy flavor channels, although we find little mass de-pendence in the heavy quarkonium channels.To explore robustness of the collectivity of heavymesons due to initial state effects in the CGC formalism,we computed the corresponding v in several slightly dif-ferent setups of approximations and model inputs. First,we computed the heavy quarkonium v by integratingover the relative momentum of the heavy quark pair Q ¯ Q analytically from 0 to infinity as an approximation inRef. [60]. In the meantime, we can also calculate the cor-relation by numerically integrating the relative momen-tum of the heavy quark pair Q ¯ Q within the kinemati-cal range according to the CEM. These two calculationsyield similar numerical results for the correlations. Fur-thermore, we find that the resulting correlation remainsalmost the same if we use a gluon or a charged hadron asthe reference particle instead of a quark as we originallychose. In addition, little dependence on various types ofheavy meson fragmentation functions was found in ournumerical calculations.Future experimental and theoretical efforts along thisline may be able to help us explore the origin of collec-tivity in high-multiplicity events in high-energy pA colli-sions. The collectivity of heavy mesons could also provideus an interesting gateway to understand the propertiesof the initial-state dense gluon matter inside high-energyhadrons. Note added:
During the completion of this manuscript,we have made the prediction for the v of the B mesonavailable to the CMS collaboration. After convolutingwith the decay kinematics for the B → D decay by usingPYTHIA [89], CMS collaboration finds our calculationwith only initial state effects to be consistent with their1recent non-prompt D meson v data[12]. ACKNOWLEDGMENTS
We thank Zhenyu Chen, Wei Li, and Feng Yuanfor useful discussions and comments. This material is partly supported by the Natural Science Foundation ofChina (NSFC) under Grant Nos. 11575070, 11775095,11890711, 11935007, and by the China Scholarship Coun-cil (CSC) under Grant No. 201906775042. CM and SYWare supported by the Agence Nationale de la Rechercheunder the project ANR-16-CE31-0019-02.We have employed Jaxodraw [90, 91] to draw the Feyn-man diagrams in this paper.
Appendix A: Detailed derivation of the correlations between heavy mesons and the reference quark
By substituting Eq. (11) into Eq. (12), we get dκ dy D dk D ⊥ = k D ⊥ N (cid:90) dφ k D ⊥ (cid:90) dy q d k q ⊥ e − i φ kD ⊥ − φ kq ⊥ ) (cid:90) d b d b d r d r d r (2 π ) (cid:90) dζζ (cid:90) dz D ( z ) z e − i (cid:126)kD ⊥ z · (cid:126)r ζ e − i(cid:126)k q ⊥ · (cid:126)r × W ( x g , x q , b , b , r , r ) k + g λ (cid:88) αβ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ ) (cid:104) DDD (cid:105)| (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ . (A1)The leading- N c terms of the scattering amplitude, (cid:104) DDD (cid:105) , do not contribute to the second harmonic. Therefore, inthis subsection, we only keep the following non-zero contributions, (cid:104)
DDD (cid:105)| (cid:126)r (cid:48) = (cid:126)r − (cid:126)r /ζ = Q s N c (cid:90) dξ (cid:90) ξ dη (cid:104) (cid:104) DDD (cid:105) + (cid:104) DDD (cid:105) − (cid:104) DDD (cid:105) − (cid:104) DDD (cid:105) (cid:105) , (A2)where, (cid:104) DDD (cid:105) = exp (cid:20) − Q s r ζ + r ) (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b + (1 − ζ ) (cid:126)r − − ζ ζ (cid:126)r ) (cid:21) exp (cid:20) ηQ s (cid:126)r ζ − (cid:126)r ) (cid:21) ( (cid:126)r ζ · (cid:126)r ) , (A3) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r + r ) (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b ) − r + r (cid:21) exp (cid:20) − ηQ s (cid:126)r · (cid:126)r (cid:21) (cid:126)r · (cid:126)r ) , (A4) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r + r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b + 1 − ζ (cid:126)r ) − r − ( (cid:126)r + (1 − ζ ) (cid:126)r ) (cid:21) exp (cid:20) − ηQ s (cid:126)r · (cid:126)r + (1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) ( (cid:126)r · (cid:126)r + (1 − ζ ) (cid:126)r · (cid:126)r ) + exp (cid:20) − Q s r + r + ( ζ + (1 − ζ ) ) r ) + 2(1 − ζ ) (cid:126)r · (cid:126)r (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b − ζ (cid:126)r ) − r − ( (cid:126)r − ζ(cid:126)r ) (cid:21) exp (cid:20) ηQ s (cid:126)r · (cid:126)r − ζ(cid:126)r · (cid:126)r ) (cid:21) ( (cid:126)r · (cid:126)r − ζ(cid:126)r · (cid:126)r ) , (A5) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r ζ + r + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b + (cid:126)r − ζ (cid:126)r ) − r − ζ r (cid:21) exp (cid:20) ηQ s ζ(cid:126)r · (cid:126)r (cid:21) ζ ( (cid:126)r · (cid:126)r )
2+ exp (cid:20) − Q s r ζ + r + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s (cid:126)b − (cid:126)b − (1 − ζ ) (cid:126)r ζ + (1 − ζ ) (cid:126)r − r − ( (cid:126)r ζ − (1 − ζ ) (cid:126)r ] (cid:21) exp (cid:20) − ηQ s ζ (cid:126)r · (cid:126)r − (1 − ζ ) (cid:126)r · (cid:126)r (cid:21) ( 1 ζ (cid:126)r · (cid:126)r − (1 − ζ ) (cid:126)r · (cid:126)r ) . (A6)The b , dependece in Eqs. (A3-A3) always takes the form exp[ − ηQ s ( (cid:126)b − (cid:126)b − (cid:126)X ) ], where (cid:126)X could be anytwo-dimensional vector.Utilizing the following relation (cid:90) d b d b (2 π ) π B p e − b b Bp e − ηQ s ( (cid:126)b − (cid:126)b − (cid:126)X ) = 14 π
11 + ηQ s B p exp (cid:20) − ηQ s X ηQ s B p ) (cid:21) , (A7)the differential κ becomes, dκ dy D dk D ⊥ = k D ⊥ N (cid:90) dφ k D ⊥ (cid:90) dy q d k q ⊥ e − i φ kD ⊥ − φ kq ⊥ ) (cid:90) d b d b d r d r d r (2 π ) (cid:90) dζζ (cid:90) dz D ( z ) z e − i (cid:126)kD ⊥ z · (cid:126)r ζ e − i(cid:126)k q ⊥ · (cid:126)r × x g f g ( x g ) x q f q ( x q ) exp (cid:20) − ∆ ( r + r )4 (cid:21) k + g λ (cid:88) αβ ψ T λαβ ( (cid:126)r ) ψ T λ ∗ αβ ( (cid:126)r (cid:48) = (cid:126)r − (cid:126)r ζ ) × Q s N c (cid:90) dξ (cid:90) ξ dη (cid:20)(cid:90) (cid:104) DDD (cid:105) + (cid:90) (cid:104) DDD (cid:105) − (cid:90) (cid:104) DDD (cid:105) − (cid:90) (cid:104) DDD (cid:105) (cid:21) , (A8)where, (cid:82) (cid:104) DDD (cid:105) i is defined as (cid:90) (cid:104) DDD (cid:105) i = (cid:90) d b d b (2 π ) π B p e − b b Bp (cid:104) DDD (cid:105) i , (A9)with (cid:90) (cid:104) DDD (cid:105) = exp (cid:20) − ηQ s ηQ s B p ) ((1 − ζ ) (cid:126)r − − ζ ζ (cid:126)r ) (cid:21) exp (cid:20) − Q s r ζ + r ) (cid:21) exp (cid:20) ηQ s (cid:126)r ζ − (cid:126)r ) (cid:21) ( (cid:126)r ζ · (cid:126)r ) , (A10) (cid:90) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r + r ) (cid:21) exp (cid:20) ηQ s r + r ) (cid:21) exp (cid:20) − ηQ s (cid:126)r · (cid:126)r (cid:21) (cid:126)r · (cid:126)r ) , (A11) (cid:90) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r + r + ( ζ + (1 − ζ ) ) r + 2(1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s ηQ s B p ) (1 − ζ ) r (cid:21) exp (cid:20) ηQ s r + ( (cid:126)r + (1 − ζ ) (cid:126)r ) ] (cid:21) exp (cid:20) − ηQ s (cid:126)r · (cid:126)r + (1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) ( (cid:126)r · (cid:126)r + (1 − ζ ) (cid:126)r · (cid:126)r ) (A12)+ exp (cid:20) − Q s r + r + ( ζ + (1 − ζ ) ) r ) + 2(1 − ζ ) (cid:126)r · (cid:126)r (cid:21) exp (cid:20) − ηQ s ηQ s B p ) ζ r (cid:21) exp (cid:20) ηQ s r + ( (cid:126)r − ζ(cid:126)r ) ] (cid:21) exp (cid:20) ηQ s (cid:126)r · (cid:126)r − ζ(cid:126)r · (cid:126)r ) (cid:21) ( (cid:126)r · (cid:126)r − ζ(cid:126)r · (cid:126)r ) , (A13) (cid:90) (cid:104) DDD (cid:105) = exp (cid:20) − Q s r ζ + r + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s ηQ s B p ) ( (cid:126)r − ζ(cid:126)r ) (cid:21) exp (cid:20) ηQ s r + ζ r ] (cid:21) (cid:20) ηQ s ζ(cid:126)r · (cid:126)r (cid:21) ζ ( (cid:126)r · (cid:126)r ) + exp (cid:20) − Q s r ζ + r + ( ζ + (1 − ζ ) ) r − − ζ ) ζ (cid:126)r · (cid:126)r ) (cid:21) exp (cid:20) − ηQ s (1 − ζ ) ηQ s B p ) ( (cid:126)r − (cid:126)r ζ ) (cid:21) exp (cid:20) ηQ s r + ( (cid:126)r ζ − (1 − ζ ) (cid:126)r ) ] (cid:21) exp (cid:20) − ηQ s (cid:126)r ζ · (cid:126)r − (1 − ζ ) (cid:126)r · (cid:126)r ) (cid:21) ( 1 ζ (cid:126)r · (cid:126)r − (1 − ζ ) (cid:126)r · (cid:126)r ) . (A14)Furthermore, we can remove the integrals over the azimuthal angles of transverse momenta using the followingrelation (cid:90) dφ k D ⊥ dφ k q ⊥ (2 π ) e − i φ kD ⊥ − φ kq ⊥ ) e − i (cid:126)kD ⊥ z · (cid:126)r ζ e − i(cid:126)k q ⊥ · (cid:126)r = J ( | (cid:126)k D ⊥ z || (cid:126)r ζ | ) J ( | (cid:126)k q ⊥ || (cid:126)r | ) cos[2( φ r − φ r )] . (A15)To compute the d r integral analytically, we can perform a Taylor expansion of the exponentials that contain thefactor (cid:126)r · (cid:126)Y , where (cid:126)Y could be any vector. The | (cid:126)r | and φ r dependences are now separated. Using the followingrelations (cid:90) π dφ r π cos(2 φ r ) cos m +2 φ r = 12 m +2 (2 m + 2)! m !( m + 2)! , (A16) (cid:90) π dφ r π cos(2 φ r ) cos m +1 φ r = 0 , (A17) (cid:90) π dφ r π sin(2 φ r ) cos l +2 φ r = 0 , (A18) (cid:90) r dr J ( | (cid:126)k q ⊥ || (cid:126)r | ) exp (cid:2) − a q r (cid:3) r m +22 = ( m + 2)! k q ⊥
16 1 a m +3 q F ( m + 3 , , − k q ⊥ a q ) , (A19)we can obtain Eq. (20). [1] V. Khachatryan et al. [CMS Collaboration], JHEP ,091 (2010).[2] S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B , 795 (2013).[3] B. Abelev et al. [ALICE Collaboration], Phys. Lett. B , 29 (2013).[4] G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. , no. 18, 182302 (2013).[5] A. Adare et al. [PHENIX Collaboration], Phys. Rev.Lett. , no. 21, 212301 (2013).[6] A. Adare et al. [PHENIX Collaboration], Phys. Rev.Lett. , no. 19, 192301 (2015).[7] V. Khachatryan et al. [CMS Collaboration], Phys. Rev.Lett. , no. 1, 012301 (2015).[8] C. Aidala et al. [PHENIX Collaboration], Nature Phys. , no. 3, 214 (2019).[9] S. Acharya et al. [ALICE Collaboration], Phys. Lett. B , 7 (2018).[10] CMS Collaboration [CMS Collaboration], CMS-PAS-HIN-18-010.[11] A. M. Sirunyan et al. [CMS Collaboration], Phys. Rev.Lett. , no. 8, 082301 (2018).[12] CMS Collaboration [CMS Collaboration], CMS-PAS- HIN-19-009.[13] P. Bozek and W. Broniowski, Phys. Rev. C , no. 1,014903 (2013).[14] A. Bzdak, B. Schenke, P. Tribedy and R. Venugopalan,Phys. Rev. C , no. 6, 064906 (2013).[15] G. Y. Qin and B. Mller, Phys. Rev. C , no. 4, 044902(2014).[16] K. Werner, M. Bleicher, B. Guiot, I. Karpenko andT. Pierog, Phys. Rev. Lett. , no. 23, 232301 (2014).[17] P. Bozek, W. Broniowski and G. Torrieri, Phys. Rev.Lett. , 172303 (2013).[18] J. L. Nagle et al. , Phys. Rev. Lett. , no. 11, 112301(2014).[19] B. Schenke and R. Venugopalan, Phys. Rev. Lett. ,102301 (2014).[20] M. Habich, J. L. Nagle and P. Romatschke, Eur. Phys.J. C , no. 1, 15 (2015).[21] P. Bozek and W. Broniowski, Phys. Lett. B , 308(2014).[22] C. Shen, J. F. Paquet, G. S. Denicol, S. Jeon and C. Gale,Phys. Rev. C , no. 1, 014906 (2017).[23] R. D. Weller and P. Romatschke, Phys. Lett. B , 351(2017). [24] W. Zhao, Y. Zhou, H. Xu, W. Deng and H. Song, Phys.Lett. B , 495 (2018).[25] X. Du and R. Rapp, JHEP , 015 (2019).[26] Z. w. Lin and D. Molnar, Phys. Rev. C , 044901 (2003).[nucl-th/0304045].[27] A. Kurkela, U. A. Wiedemann and B. Wu, Phys. Lett. B , 274 (2018).[28] H. Li, Z. W. Lin and F. Wang, Phys. Rev. C , no. 4,044911 (2019).[29] A. Kurkela, U. A. Wiedemann and B. Wu, Eur. Phys. J.C , no. 9, 759 (2019).[30] N. Armesto, L. McLerran and C. Pajares, Nucl. Phys. A , 201 (2007).[31] A. Dumitru, F. Gelis, L. McLerran and R. Venugopalan,Nucl. Phys. A , 91 (2008).[32] S. Gavin, L. McLerran and G. Moschelli, Phys. Rev. C , 051902 (2009).[33] A. Dumitru and J. Jalilian-Marian, Phys. Rev. D ,094015 (2010).[34] A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian,T. Lappi and R. Venugopalan, Phys. Lett. B , 21(2011).[35] A. Kovner and M. Lublinsky, Phys. Rev. D , 034017(2011).[36] Y. V. Kovchegov and D. E. Wertepny, Nucl. Phys. A ,50 (2013).[37] K. Dusling and R. Venugopalan, Phys. Rev. Lett. ,262001 (2012).[38] Y. V. Kovchegov and D. E. Wertepny, Nucl. Phys. A ,254 (2014).[39] A. Dumitru and A. V. Giannini, Nucl. Phys. A , 212(2015).[40] A. Dumitru, L. McLerran and V. Skokov, Phys. Lett. B , 134 (2015).[41] A. Dumitru and V. Skokov, Phys. Rev. D , no. 7,074006 (2015).[42] T. Lappi, Phys. Lett. B , 315 (2015).[43] B. Schenke, S. Schlichting and R. Venugopalan, Phys.Lett. B , 76 (2015).[44] T. Lappi, B. Schenke, S. Schlichting and R. Venugopalan,JHEP , 061 (2016).[45] L. McLerran and V. Skokov, Nucl. Phys. A , 83(2017).[46] A. Kovner, M. Lublinsky and V. Skokov, Phys. Rev. D , no. 1, 016010 (2017).[47] E. Iancu and A. H. Rezaeian, Phys. Rev. D , no. 9,094003 (2017).[48] K. Dusling, M. Mace and R. Venugopalan, Phys. Rev.Lett. , no. 4, 042002 (2018).[49] K. Dusling, M. Mace and R. Venugopalan, Phys. Rev. D , no. 1, 016014 (2018).[50] K. Fukushima and Y. Hidaka, JHEP , 114 (2017).[51] Y. V. Kovchegov and V. V. Skokov, Phys. Rev. D ,no. 9, 094021 (2018).[52] D. Boer, T. Van Daal, P. J. Mulders and E. Petreska,JHEP , 140 (2018).[53] M. Mace, V. V. Skokov, P. Tribedy and R. Venugopalan,Phys. Rev. Lett. , no. 5, 052301 (2018) Erratum:[Phys. Rev. Lett. , no. 3, 039901 (2019)].[54] M. Mace, V. V. Skokov, P. Tribedy and R. Venugopalan,Phys. Lett. B , 161 (2019).[55] T. Altinoluk, N. Armesto, A. Kovner and M. Lublinsky,Eur. Phys. J. C , no. 9, 702 (2018).[56] A. Kovner and V. V. Skokov, Phys. Lett. B , 372 (2018).[57] A. Kovner and A. H. Rezaeian, Phys. Rev. D , no. 7,074018 (2017).[58] A. Kovner and A. H. Rezaeian, Phys. Rev. D , no. 7,074008 (2018).[59] M. K. Davy, C. Marquet, Y. Shi, B. W. Xiao andC. Zhang, Nucl. Phys. A , 293 (2019).[60] C. Zhang, C. Marquet, G. Y. Qin, S. Y. Wei andB. W. Xiao, Phys. Rev. Lett. , no. 17, 172302 (2019).[61] F. Dominguez, C. Marquet, B. W. Xiao and F. Yuan,Phys. Rev. D , 105005 (2011).[62] F. Dominguez, C. Marquet, A. M. Stasto and B. W. Xiao,Phys. Rev. D , 034007 (2013).[63] Y. Q. Ma and R. Venugopalan, Phys. Rev. Lett. , no.19, 192301 (2014).[64] Y. Q. Ma, R. Venugopalan and H. F. Zhang, Phys. Rev.D , 071901 (2015).[65] K. Watanabe and B. W. Xiao, Phys. Rev. D , no. 11,111502 (2015).[66] Y. Q. Ma, P. Tribedy, R. Venugopalan and K. Watanabe,Phys. Rev. D , no. 7, 074025 (2018).[67] F. Gelis and A. Peshier, Nucl. Phys. A , 879 (2002).[68] J. P. Blaizot, F. Gelis and R. Venugopalan, Nucl. Phys.A , 57 (2004).[69] J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D , 114017 (2004) Erratum: [Phys. Rev. D , 079901(2005)].[70] F. Dominguez, C. Marquet and B. Wu, Nucl. Phys. A , 99 (2009).[71] C. Marquet and H. Weigert, Nucl. Phys. A , 68(2010).[72] Y. Shi, C. Zhang and E. Wang, Phys. Rev. D , no. 11,116014 (2017).[73] C. Zhang, M. K. Davy, Y. Shi and E. Wang, Phys. Rev.D , no. 3, 034009 (2019).[74] L. D. McLerran and R. Venugopalan, Phys. Rev. D ,2233 (1994).[75] L. D. McLerran and R. Venugopalan, Phys. Rev. D ,3352 (1994).[76] C. Peterson, D. Schlatter, I. Schmitt and P. M. Zerwas,Phys. Rev. D , 105 (1983)..[77] T. Kneesch, B. A. Kniehl, G. Kramer and I. Schienbein,Nucl. Phys. B , 34 (2008).[78] B. A. Kniehl and G. Kramer, Phys. Rev. D , 037502(2006).[79] B. A. Kniehl, G. Kramer, I. Schienbein and H. Spies-berger, Phys. Rev. D , 014011 (2008).[80] J. W. Qiu, P. Sun, B. W. Xiao and F. Yuan, Phys. Rev.D , no. 3, 034007 (2014).[81] N. Borghini, P. M. Dinh and J. Y. Ollitrault, Phys. Rev.C , 054901 (2001).[82] H. Fujii, F. Gelis and R. Venugopalan, Nucl. Phys. A , 146 (2006)[83] H. Fujii and K. Watanabe, Nucl. Phys. A , 78 (2013).[84] A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt,Eur. Phys. J. C , 189 (2009).[85] S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys. B , 42 (2008).[86] K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D ,014017 (1998) [hep-ph/9807513].[87] D. Acosta et al. [CDF Collaboration], Phys. Rev. D ,032001 (2005).[88] G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. , no. 8, 082301 (2020). [89] T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys.Commun. , 852 (2008).[90] D. Binosi and L. Theussl, Comput. Phys. Commun. ,76 (2004). [91] D. Binosi, J. Collins, C. Kaufhold and L. Theussl, Com-put. Phys. Commun.180