Prepared for submission to JHEP An Effective Theory of Quarkonia in QCD Matter

For heavy quarkonia of moderate energy, we generalize the relevant successful theory, non-relativistic Quantum Chromodynamics (NRQCD), to include interactions in nuclear matter. The new resulting theory, NRQCD with Glauber gluons, provides for the first time a universal microscopic description of the interaction of heavy quarkonia with a strongly interacting medium, consistently applicable to a range of phases, such as cold nuclear matter, dense hadron gas, and quark-gluon plasma. The effective field theory we present in this work is derived from first principles and is an important step forward in understanding the common trends in proton-nucleus and nucleus-nucleus data on quarkonium suppression. ar X iv :1 90 6. 04 18 6v 1 [ he pph ] 1 0 Ju n 20 19


Introduction
It is widely believed today that novel phases of nuclear matter, such as the quark-gluon plasma (QGP) and a hot, dense gas of hadrons, are integral and important parts of the evolution of the early universe. These extreme environments are inaccessible to direct observation, but can be recreated in the laboratory by colliding heavy nuclei at relativistic energies. One of the main goals of nuclear physics is to accurately determine the properties of these new states of matter [1]. Since their lifetimes are very short, of order 10 −23 s, one must use the produced particles themselves to probe the QGP and the hadron gas. Quarkonia have emerged as premier diagnostics of the QGP. It was predicted that, when immersed in the plasma characterized by very high temperature, the color interaction between the heavy quarks will be screened and quarkonia will dissociate [2]. Excited, weakly-bound states are expected to melt away first, ground tightly-bound states are expected to melt away last, provide a way to determine the plasma temperature [3].
In the past decade phenomenological studies of quarkonia have evolved significantly to include effects that range from heavy quark recombination to dissociation through collisional interactions of J/ψ and Υ states propagating through the QGP [4][5][6][7][8]. The physics input in such calculations comes from the hard thermal loop calculations of the real and imaginary parts of the heavy quark-antiquark potentials [9,10], lattice QCD calculations [11], a T −matrix approach [12] to obtain interaction and decay rates of thermal states, and lightcone wavefunction approach to obtain the dissociation rate of quarkonia from collisional and thermal effects [13]. The evolution of the quarkonium system has been described by rate equations [13,14], stochastic equations [15][16][17] such as the Lindblad equation, and the Boltzmann equation [18]. Those studies has focused almost exclusively on quarkonia in a thermal QGP medium.
In spite of the advances described above, a fully coherent theoretical picture of quarkonium production at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) has not yet emerged. In proton-nucleus (p+A) collisions, where QGP is much less likely to be formed, attenuation similar to the one seen in nucleus-nucleus (A+A) reactions is still observed, albeit of smaller magnitude. Even in high multiplicity proton-proton (p+p) collisions there is evidence for Υ(2S) disappearance as a function of the hadronic activity (N tracks ) in the event. Specifically, the relative suppression of the excited versus ground bottomonium states Υ(2S)/Υ(1S) as a function of the number of charged particle tracks, shows the same dissociation trend for high-multiplicity proton-proton, proton-lead, and lead-lead reactions at the LHC [19]. This experimental finding has not yet found satisfactory theoretical expectations. It was argued very recently that quarkonium dissociation by co-movers might be responsible for those trends [20]. Differential ψ , χ c and Υ suppression was also established at RHIC [21,22] in d+Au reactions. Upcoming experimental detector upgrades at RHIC and luminosity upgrades at the LHC will allow extensive studies of J/ψ and Υ states with improved precision in high-multiplicity hadronic and nuclear collisions. There is an opportunity to further develop microscopic QCD approaches that describe this quarkonium physics in nuclear matter and that will facilitate the quantitative determination of the transport properties of the QGP and the hadron gas.
With this motivation, we first notice that calculations of heavy quarkonium production encounter hierarchies of momentum and mass scales, which is precisely where effective filed theories (EFTs) excel in reducing theoretical uncertainties and improving computational accuracy [23]. Usually the scales one encounters are p T , m Q , m Q λ, m Q λ 2 , and Λ QCD , where p T is the quarkonium transverse momentum, m Q the heavy quark mass, and λ the heavy quark-antiquark pair relative velocity in the quarkonium rest frame. For moderate and high transverse momentum p T 2m Q the established and most successful theory that describes quarkonium production and decays is non-relativistic QCD (NRQCD) [24]. Many recent theoretical studies take full advantage of the EFT capabilities to significantly boost the theoretical precision of J/ψ and Υ analyses and propose modern observables [25] that can probe the quarkonium production mechanisms. Most of those studies focus their efforts on quarkonium states in the high energy (E m QQ ) region, where theoretical advances are now possible based upon NRQCD, SCET [26][27][28][29], and the picture of parton fragmentation [30,31].
The challenge that we face is to develop a microscopic theory of quarkonia applicable to different phases of nuclear matter in p+A and A+A reactions. We approach this challenge from the effective field theory point of view. The distinct advantage of an EFT approach is that it can provide a model-independent description of the universal physics of energetic particle production in the background of a QCD medium. This universal description can be applied equally well to the QGP or to a hadron gas, with model dependence entering only in the choice of the medium. In the past several years there were important developments in applying an EFT approach to describe particle production in the presence of strongly interacting matter. Particularly relevant to this work is the formulation and application of an effective theory of QCD, soft collinear effective theory with Glauber gluons (SCET G ) [32,33] for light particles (π 0,± , K ± , · · · ). It was also demonstrated that rigorous treatment of heavy flavor in matter is possible by constructing the necessary extension of SCET G to nonzero quark masses, giving us the applicable theory for energetic mesons containing a single heavy quark [34]. SCET G allowed us for the first time, to overcome known limitations of traditional phenomenological approaches, use the same computational techniques in high energy and heavy ion physics, and increase the accuracy and quantify the theoretical uncertainties in the calculations of light particle [35,36] and heavy meson [34] production in A+A reactions.
As is the case in the vacuum, production of quarkonia in nuclear matter remains a multiscale problem. For this reason, we identify the EFT approach the correct way to attack it. In this paper we demonstrate how one can generalize NRQCD to incorporate interactions of the non-relativistic heavy quarks with the medium. This is achieved through incorporating the Glauber and Coulomb gluon exchanges of the heavy quarks with three different sources: collinear, soft, and static. We believe this version of NRQCD will facilitate a much more robust and accurate theoretical analysis of the wealth of quarkonium measurements in dense QCD matter.
The outline of this paper is as follows: In Section 2, after a brief overview of NRQCD, we explore the applicability of the well-established energy loss approach to quarkonia. We take the leading power factorization limit, where a quarkonium state is produced thought the fragmentation process from a parton that undergoes energy loss in matter and demonstrate that the predicted magnitude and hierarchy of suppression for ground and excited charmonium states is not compatible with the experimental data. With this in mind, we, consider the propagation of the quarkonium state itself in QCD matter in Section 3. The possible off-shell gluon exchanges between the heavy quark/antiquark and the medium are discussed for several sources of scattering and we identify two relevant modes that mediate the interaction: Coulomb and Glauber gluons. In the following Section 4, we give the Lagrangian and derive the Feynman rules for such exchanges. Finally, we conclude in Section 5. We discuss how a self-consistent background field approach to quarkonium propagation in matter can be formulated in Appendix A.

Energy loss approach within the NRQCD formalism
Before we proceed to the formulation of a generic effective theory of quarkonium production in matter, we have to explore whether medium-induced radiative processes might contribute significantly to the modification of quarkonium cross sections in reactions with nuclei. It was suggested [37,38] that such effects can reduce the cross section of high transverse momentum J/ψ production at the LHC [39,40].
After we give a brief review of vNRQCD we proceed by describing the leading power factorization of NRQCD for quarkonium production and introduce the quarkonium fragmentation functions within the NRQCD framework. We then apply energy loss effects to obtain quarkonium production rates in medium.

Non-relativistic QCD: a brief overview
In the quarkonium rest frame, the heavy quark and antiquark have small relative velocity, (λ 2 ∼ 0.1 for bottomonium and λ 2 ∼ 0.3 for charmonium). Therefore, NRQCD, which is an effective field theory that describes Quantum Chromodynamics in the non-relativistic limit, provides the correct theoretical framework for studying their interactions.
There are three important scales that appear when studying the dynamics of nonrelativistic heavy quarks: the mass of the heavy quark, m, the size of their momentum in the quarkonium rest frame, mλ, and their kinetic energy, mλ 2 . The distance r ∼ 1/(mλ) gives an estimate on the size of the quarkonium state and the separation between the heavy quark-antiquark pair. The non-relativistic kinetic energy ∆E ∼ mλ 2 is of the same order as the energy splittings of radial excitations. We refer to mλ and mλ 2 as the soft and ultra-soft scales respectively. Correspondingly, gluons that have all of their four-momentum components scaling as mλ and mλ 2 are called soft and ultra-soft gluons. While the ultra-soft scale is well within the non-perturbative regime the soft scale is about 1.5 GeV for both bottomonium and charmonium.
The effective theory of vNRQCD is a version non-relativistic QCD introduced in Ref. [41] and recently formulated in a manifestly gauge invariant form in Ref. [42]. What we find appealing about this version of NRQCD is the clear distinction of soft and ultra-soft degrees of freedom and the use of label-momentum notation. Both of those aspects are essential for the purposes of our work. We work in the limit where the measurement is sensitive to the kinematics of the heavy quark-antiquark pair (in the quarkonium rest frame) and therefore is critical we can separate the various infrared degrees of freedom. Using the four-vector v µ = (1, 0), the four-momenta of the heavy quark, p, can be written as follows, where r 0 is the kinetic energy and r is the three momentum of the heavy quark. Since the heavy quarks we consider are on-shell, i.e. p 2 = m 2 , then in the non-relativistic limit where the three momentum is small compared to the mass, |r| ∼ λm, with λ 1 we have which has solution only if r µ ∼ (λ 2 , λ). In the presence of both soft and ultra-soft modes, it is important to decompose the small momentum component in its soft (label) and ultra-soft (residual) parts, where r µ us ∼ (λ 2 , λ 2 , λ 2 , λ 2 ), and r µ s ∼ (0, λ 1 , λ 1 , λ 1 ). Then the connection with the convention in Eq. (2.1) can be made with the replacement, The QCD heavy quark field (Ψ) can then be decomposed in the vNRQCD heavy quark field (ψ (x)) as follows, where are the label components of the heavy quark momentum and x is the coordinate space conjugate of the residual components. The soft (A µ ) and ultra-soft (A µ us ) gluon fields have momenta which scale (all four components) as soft (∼ mλ) or ultra-soft (∼ mλ 2 ) respectively.
The Lagrangian of the EFT can then be written in terms of those fields in the following form [41,42], where ψ denotes the heavy quark field and χ the corresponding antiquark. The Lagrangian terms L (2) are higher order terms, L s is the soft gluon and ghost part of the Lagrangian, and L V contains the potential terms which have the following generic structure, Double soft gluon emissions: Interactions with soft fermions: Heavy quark-antiquark potential: where U µ,ν , Z µ , and V are functions of the momenta of the field included in the corresponding interactions. The soft fermion fields,φ , acting on the vacuum creates a light quark with soft momenta, µ ∼ (λ, λ, λ, λ), and similarly φ for the antiquark. The Lagrangian that describes the interaction of soft fermions with soft gluons is identical to QCD, see Ref [42]. The label momentum operator [28], P µ = (P 0 , −P), is defined such that it projects only onto the label momentum space, and the covariant derivative is iD µ ≡ i∂ µ − gA µ us (x). In collider physics, quarkonium production is studied within the NRQCD factorization conjecture, based on which the cross section is written as a sum of products of short distance matching coefficients and the corresponding long distance matrix elements (LDMEs) (2.8) The short distance coefficients (SDCs), dσ ij→QQ[n]+X , describe the production of the QQ[n] pair in a particular angular momentum and color configuration, J . In the case of hadronic initial states, SDCs are expressed as a convolution of the partonic cross section and the collinear PDFs. The partonic cross section is then calculated in the matching of NRQCD onto QCD as an expansion in the strong coupling constant [43][44][45][46][47][48][49][50]. In contrast, the LDMEs, O Q (n) , describe the decay of the QQ[n] pair into the final color-singlet quarkonium state, Q, through soft and ultra-soft gluon emissions. LDMEs are universal and fundamentally non-perturbative objects, and need to be extracted from experiment [50][51][52][53][54]. Although in principle all possible intermediate QQ[n] configurations contribute to the final quarkonium state, LDMEs scale with powers of λ, thus, we can truncate the sum up to the desired accuracy.

Quarkonium fragmentation functions
In order to envision energy loss processes as contributors to the modification of quarkonium cross sections in QCD matter two conditions must be satisfied. First, quarkonium production must be expressed as fragmentation of partons into the various J/ψ and Υ states. The energy of the hard parton is then reduced through inelastic processes in matter prior to fragmentation. Second, the process of fragmentation of quarkonia must happen at time scales larger than the size of the QCD medium, τ f orm ≥ L. This condition must also be investigated phenomenologically in reactions with nuclei, as the simpler hadronic collisions do not give relevant constraints.
Fortunately, in the last decade a leading power (LP) factorization of NRQCD has been established [55][56][57][58][59][60] and is expected to hold at high transverse momenta (p T m Q ). In the large transverse momentum limit the NRQCD short distance coefficients suffer from logarithmic enhancements of the form α m s ln n (p T /2m Q ). These terms could spoil the perturbative expansion and, thus, resummation is necessary in order to make meaningful predictions. This is achieved through the LP factorization of NRQCD, where the cross section is now factorized into short distance matching coefficients (that describe the production and propagation of a parton k) and the so called NRQCD fragmentation functions, The dependence on the factorization scale, µ, of the factorized terms is exactly what allows for the resummation of large logarithms through the use of renormalization group techniques and, particularly, the DGLAP evolution for the fragmentation functions. Comparison of the above equation with Eq. (2.8) immediately gives that the NRQCD fragmentation functions can be written in terms of the same LDMEs that appear in the fixed order factorization and perturbatively calculable matching coefficients, ) pair. The LP factorization is expected to hold for p T m Q but the precise p T region of validity cannot be be determined analytically. However, phenomenological applications to charmonia have shown that it may hold to transverse momenta as low as p T = 10 GeV [50]. U E c x T g C M q C F Z s p g n C g u p f I Z 4 g g b D S W V X q v 8 / k y 2 M V P u b R 2 M t B r J L u e c N u N q y 7 i 1 r r t A i p D I 7 B C T g D N r g E L X A D 2 q A D M E j A M 3 g B r 8 a T 8 W a 8 G x + L 1 p J R z B y B P z A + f w C b v J x p < / l a t e x i t > 3 S R B g c 9 I 1 5 / e Z H 7 3 k Q h J Q / 6 g Z h F x A z T m d E Q x U l r y z H K S D p J 6 e u / Z g 8 R p u K l n V q y a N Q d c J 3 Z O K i B H y z N / + s M Q x w H h C j M k p W N b k X I T J B T F j K S l f i x J h P A U j Y m j K U c B k W 4 y f z y F V a 0 M 4 S g U u r i C c 3 V 5 I k G B l L P A 1 5 0 B U h O 5 6 m X i f 5 4 T q 1 H D T S i P Y k U 4 X h w a x Q y q E G Y p w C E V B C s 2 0 w R h Q f W v E E + Q Q F j p r E r V 5 T P Z 8 k g F T 1 k 0 9 m o Q 6 6 R z W b P r N e v u q t I 8 z 0 M q g l N w B i 6 A D a 5 B E 9 y C F m g D D G L w A l 7 B m / F s v B s f x u e i t W D k M y f g D 4 y v X 4 N G n F o = < / l a t e x i t > 3 S T h J C R c Y Y a k d G w r V m 6 K h K K Y k a w y S C S J E Z 6 i M X E 0 5 S g k 0 k 3 n j 2 e w r p U R D C K h i y s 4 V 3 9 P p C i U c h b 6 u j N E a i K X v V z 8 z 3 M S F V y 5 K e V x o g j H i 0 N B w q C K Y J 4 C H F F B s G I z T R A W V P 8 K 8 Q Q J h J X O q l L / f S Z f H q v w M Y / G X g 5 i l X T P G 3 a z Y d 1 e 1 F q n R U h l c A x O w B m w w S V o g R v Q B h 2 A Q Q K e w Q t 4 N Z 6 M N + P d + F i 0 l o x i 5 g j 8 g f H 5 A 3 g C n F M = < / l a t e x i t > 3 P L e y + B / n h U p r 2 r H 1 A 8 j R X w 8 W e R F D K o A Z u n A H h U E K z Z K B c K C p m + F e I A E w i r N s F C a X Z N d H i r + m E V j z g e x K J p n Z b N S N m 7 P i 7 W j a U h 5 c A A O w Q k w w Q W o g W t Q B w 2 A w R N 4 A a / g T X v W 3 r U P 7 X P S m t O m M / v g T 2 n f v 6 Z K p d w = < / l a t e x i t > g < l a t e x i t s h a 1 _ b a s e 6 4 = " G 4 + 2 j T g h m x 2 t W + 9 C C n R t B 2 5 H w K o = " > A A A B / n i c b V D L S g M x F M 3 4 r P V V d e k m W C q u y o w K u i y 4 c d m C f U A 7 l E x 6 p w 1 N M k O S E c t Q c O 9 W f 8 G d u P V X / A M / w 0 w 7 i 9 p 6 I H A 4 5 7 5 y g p g z b V z 3 2 1 l b 3 9 j c 2 i 7 s F H f 3 9 g 8 O S 0 f H L R 0 l i k K T R j x S n Y B o 4 E x C 0 z D D o R M r I C L g 0 A 7 G d 5 n f f g S l W S Q f z C Q G X 5 C h Z C G j x F i p M e y X y m 7 V n Q G v E i 8 n Z Z S j 3 i / 9 9 A Y R T Q R I Q z n R u u u 5 s f F T o g y j H K b F X q I h J n R M h t C 1 V B I B 2 k 9 n h 0 5 x x S o D H E b K P m n w T F 3 s S I n Q e i I C W y m I G e l l L x P / 8 7 q J C W / 9 l M k 4 M S D p f F G Y c G w i n P 0 a D 5 g C a v j E E k I V s 7 d i O i K K U G O z K V Y W 1 2 T D Y y O e p j Y a b z m I V d K 6 r H p X V b d x X a 6 d 5 y E V 0 C k 6 Q x f I Q z e o h u 5 R H T U R R Y B e 0 C t 6 c 5 6 d d + f D + Z y X r j l 5 z w n 6 A + f r F 2 t 0 l l o = < / l a t e x i t > q < l a t e x i t s h a 1 _ b a s e 6 4 = " U L s L f k i n w F x u + c Y I f x P 9 c 1 L A S K 0 = " H a r 7 g x 4 l X g 5 K a M c 9 X 7 p p z e I a C K Z A i q I M V 3 P j c F P i Q Z O B Z s W e 4 l h M a F j M m R d S x W R z P j p 7 N A p r l h l g M N I 2 6 c A z 9 T F j p R I Y y Y y s J W S w M g s e 5 n 4 n 9 d N I L z 1 U o j p q I I o Z e 0 C t 6 c 5 6 d d + f D + Z y X r j l 5 z w n 6 A + f r F 0 g 4 l k Q = < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " u z t T e 6 u I R U g g b G 1 G p u r g m G y 4 N f 5 r a a P z l I F Z J 6 6 L m X 9 a 8 + 6 t K / S w P q Q g n c A r n 4 M M 1 1 O E O G t A E D A x e 4 B X e n G f n 3 f l w P u e l B S f v O Y Y / c L 5 + A f Y Q m W 0 = < / l a t e x i t > 2 s < l a t e x i t s h a 1 _ b a s e 6 4 = " r U v D Y p 8 z J i i M r l 4 W p y 8 J I U Z H r I g = " > A A A C B 3 i c b V D L T g I x F L 3 j E / G F u n T T S D C u y A y a 6 J L E j U t M 5 B F h J J 3 S g Y a 2 M 2 k 7 R k L 4 A P d u 9 R f c G b d + h n / g Z 9 i B W S B 4 k i Y n 5 9 x X T x B z p o 3 r f j s r q 2 v r G 5 u 5 r f z 2 z u 7 e f u H g s K G j R B F a J x G P V C v A m n I m a d 0 w w 2 k r V h S L g N N m M L x O / e Y j V Z p F 8 s 6 M Y u o L 3 J c s Z A Q b K 9 1 3 M I 8 H u K s f K t 1 C 0 S 2 7 U 6 B l 4 m W k C B l q 3 c J P p x e R R F B p C M d a t z 0 3 N v 4 Y K 8 M I p 5 N 8 J 9 E 0 x m S I + 7 R t q c S C a n 8 8 v X i C S l b p o T B S 9 k m D p u p 8 x x g L r U c i s J U C m 4 F e 9 F L x P 6 + d m P n Q / n c 1 a 6 4 m Q 9 R / A H z t c v L a a a E Q = = < / l a t e x i t > 3 s < l a t e x i t s h a 1 _ b a s e 6 4 = " R g j X + x F + R a 4 q B w / n 4 7 S 1 h X X f w j I = " < l a t e x i t s h a 1 _ b a s e 6 4 = " r U v D Y p 8 z J i i M r l 4 W p y 8 J I U Z H r I g = " > A A A C B 3 i c b V D L T g I x F L 3 j E / G F u n T T S D C u y A y a 6 J L E j U t M 5 B F h J J 3 S g Y a 2 M 2 k 7 R k L 4 A P d u 9 R f c G b d + h n / g Z 9 i B W S B 4 k i Y n 5 9 x X T x B z p o 3 r f j s r q 2 v r G 5 u 5 r f z 2 z u 7 e f u H g s K G j R B F a J x G P V C v A m n I m a d 0 w w 2 k r V h S L g N N m M L x O / e Y j V Z p F 8 s 6 M Y u o L 3 J c s Z A Q b K 9 1 3 M I 8 H u K s f K t 1 C 0 S 2 7 U 6 B l 4 m W k C B l q 3 c J P p x e R R F B p C M d a t z 0 3 N v 4 Y K 8 M I p 5 N 8 J 9 E 0 x m S I + 7 R t q c S C a n 8 8 v X i C S l b p o T B S 9 k m D p u p 8 x x g L r U c i s J U C m 4 F e 9 F L x P 6 + d m P n Q / n c 1 a 6 4 m Q 9 R / A H z t c v L a a a E Q = = < / l a t e x i t > 2 s < l a t e x i t s h a 1 _ b a s e 6 4 = " r U v D Y p 8 z J i i M r l 4 W p y 8 J I U Z H r I g = " > A A A C B 3 i c b V D L T g I x F L 3 j E / G F u n T T S D C u y A y a 6 J L E j U t M 5 B F h J J 3 S g Y a 2 M 2 k 7 R k L 4 A P d u 9 R f c G b d + h n / g Z 9 i B W S B 4 k i Y n 5 9 x X T x B z p o 3 r f j s r q 2 v r G 5 u 5 r f z 2 z u 7 e f u H g s K G j R B F a J x G P V C v A m n I m a d 0 w w 2 k r V h S L g N N m M L x O / e Y j V Z p F 8 s 6 M Y u o L 3 J c s Z A Q b K 9 1 3 M I 8 H u K s f K t 1 C 0 S 2 7 U 6 B l 4 m W k C B l q 3 c J P p x e R R F B p C M d a t z 0 3 N v 4 Y K 8 M I p 5 N 8 J 9 E 0 x m S I + 7 R t q c S C a n 8 8 v X i C S l b p o T B S 9 k m D p u p 8 x x g L r U c i s J U C m 4 F e 9 F L x P 6 + d m P n Q / n c 1 a 6 4 m Q 9 R / A H z t c v L a a a E Q = = < / l a t e x i t > 2 s < l a t e x i t s h a 1 _ b a s e 6 4 = " r U v D Y p 8 z J i i M r l 4 W p y 8 J I U Z H r I g = " > A A A C B 3 i c b V D L T g I x F L 3 j E / G F u n T T S D C u y A y a 6 J L E j U t M 5 B F h J J 3 S g Y a 2 M 2 k 7 R k L 4 A P d u 9 R f c G b d + h n / g Z 9 i B W S B 4 k i Y n 5 9 x X T x B z p o 3 r f j s r q 2 v r G 5 u 5 r f z 2 z u 7 e f u H g s K G j R B F a J x G P V C v A m n I m a d 0 w w 2 k r V h S L g N N m M L x O / e Y j V Z p F 8 s 6 M Y u o L 3 J c s Z A Q b K 9 1 3 M I 8 H u K s f K t 1 C 0 S 2 7 U 6 B l 4 m W k C B l q 3 c J P p x e R R F B p C M d a t z 0 3 N v 4 Y K 8 M I p 5 N 8 J 9 E 0 x m S I + 7 R t q c S C a n 8 8 v X i C S l b p o T B S 9 k m D p u p 8 x x g L r U c i s J U C m 4 F e 9 F L x P 6 + d m P n Q / n c 1 a 6 4 m Q 9 R / A H z t c v L a a a E Q = = < / l a t e x i t > 3 s < l a t e x i t s h a 1 _ b a s e 6 4 = " R g j X + x F + R a 4 q B w / n 4 7 S 1 h X X f w j I = " < l a t e x i t s h a 1 _ b a s e 6 4 = " r U v D Y p 8 z J i i M r l 4 W p y 8 J I U Z H r I g = " > A A A C B 3 i c b V D L T g I x F L 3 j E / G F u n T T S D C u y A y a 6 J L E j U t M 5 B F h J J 3 S g Y a 2 M 2 k 7 R k L 4 A P d u 9 R f c G b d + h n / g Z 9 i B W S B 4 k i Y n 5 9 x X T x B z p o 3 r f j s r q 2 v r G 5 u 5 r f z 2 z u 7 e f u H g s K G j R B F a J x G P V C v A m n I m a d 0 w w 2 k r V h S L g N N m M L x O / e Y j V Z p F 8 s 6 M Y u o L 3 J c s Z A Q b K 9 1 3 M I 8 H u K s f K t 1 C 0 S 2 7 U 6 B l 4 m W k C B l q 3 c J P p x e R R F B p C M d a t z 0 3 N v 4 Y K 8 M I p 5 N 8 J 9 E 0 x m S I + 7 R t q c S C a n 8 8 v X i C S l b p o T B S 9 k m D p u p 8 x x g L r U c i s J U C m 4 F e 9 F L x P 6 + d m P n Q / n c 1 a 6 4 m Q 9 R / A H z t c v L a a a E Q = = < / l a t e x i t > 3 s < l a t e x i t s h a 1 _ b a s e 6 4 = " R g j X + x F + R a 4 q B w / n 4 7 S 1 h X X f w j I = " In this work we consider both the direct production and the feed-down from decays of excited quarkonium states. For J/ψ the following feed-down contributions are implemented, For the direct fragmentation of a parton to J/ψ and ψ(2S) we consider the following intermediate QQ states: 3 S 1 channel, for each other channel we only conciser the leading in α s contribution. As a result, the various channels will be evaluated at different order in the perturbative expansion. For the case 3 S where the leading mechanism is the heavy quark fragmentation, in addition we include the gluon channel due to the abundance of gluons in hadronic collisions. These contributions are summarized in Figure 1.
The dominant production channels for the χ cJ come from the intermediate QQ[n] → χ cJ states for which n ∈ { 3 P 1 }. For these mechanisms, we identify the gluon and heavy quark initiating processes to be the most relevant, see Figure 1. Therefore, the fragmentation functions we need for our analysis are: i/Q to an arbitrary scale µ > 2m c we use the standard DGLAP evolution [61][62][63] at leading logarithmic (LL) accuracy. From Ref. [46] we have, For the same channel, the heavy quark short distance coefficients are given by: whereD J (z, 2m c ) are given in Eq. (3.3) of Ref. [47]. For the octet production mechanism, 3 S 1 , also present in the case of ψ(nS), we have (see Refs. [48,49]): Our analysis for the direct production of J/ψ and ψ(2S) follows Ref. [30]. All relevant fragmentation functions and the corresponding Mellin transforms are collected in the Appendix of Ref. [30]. A comprehensive analysis and extraction of the non-perturbative LDMEs, consistent with LP factorization, is given by Ref. [50]. Throughout this paper we use their results for the values of the LDMEs.

Medium-induced energy loss
Let us now turn to the application of energy loss to quarkonium production. If a parton c loses momentum fraction during its propagation in the medium to escape with momentum p med Tc , in the short distance hard process its momentum is given by p Tc = p med Tc /(1 − ). This also gives rise to an additional Jacobian factor |d 2 p med Tc /d 2 p Tc | = (1 − ) 2 , similar to the z 2 factor in the factorization formula for hadron production. The cross section for hadron production and quarkonium production per elementary nucleon-nucleon (N N ) collision in the leading power limit is then written down as (2.16) In Eq. (2.16) we have omitted the renormalization and factorization scale dependences for brevity. P ( ) is the probability distribution for the hard parton c to lose energy due to multiple gluon emission, dσ c (p T ) dyd 2 p Tc is the hard partonic cross section, and N coll. is the average number of binary nucleon-nucleon colliions.  In the approximation that the fluctuations of the average number of medium-induced gluons are uncorrelated [64,65], the spectrum of the total radiative energy loss fraction due to multiple gluon emissions, = i ω i /E, can be expressed via a Poisson expansion P ( , E) = ∞ n=0 P n ( , E), with P 1 ( , E) = e − N g ρ( , E). We note that in our notation ρ(x, E) is the medium-induced gluon spectrum where x = ω/E is the fraction of the energy of the parent parton taken by an individual gluon and x 0 = Λ QCD /2E. We keep explicitly the dependence on the parent parton energy but remark that medium-induced gluon radiation also depends on the parton's flavor and mass. The terms of the Poisson series are generated iteratively as follows We note that in the presence of a medium radiation is attenuated at the typical Debye screening scale and the number of medium-induced gluons is finite. Therefore, we have explicitly a finite n = 0 no radiation contribution P 0 ( , E) = e − N g (E) δ( ). The normalized Poisson distribution that enters Eq. (2.16) then gives Several formalisms have been developed in the literature to evaluate medium-induced gluon radiation [66][67][68][69][70][71]. In this work, we use the soft gluon emission limit of the full in-medium splitting kernels [34,72,73] and evaluate them in a viscous 2+1 dimensional hydrodynamic model of the background medium [74].
We now turn to the evaluation of the prompt J/ψ and ψ(2S) suppression in lead-lead (Pb+Pb) collisions at the LHC. We calculate the partonic cross sections as in Ref. [36]. We chose the values of the coupling between the hard partons and the QCD medium that they propagate in to be in the range g = 1.7 − 1.9. These values are slightly smaller than the ones used in [36] and the difference can be traced to the different hydrodynamic models of the medium. Earlier works used ideal Bjorken expanding medium with purely gluonic degrees of freedom. As we will show below, the suppression of quarkonia, especially the J/ψ, obtained in the energy loss framework is too large when compared to experimental measurements. Thus, if there is an uncertainty in the choice of the coupling constant g, we must err on the side of smaller couplings. A larger coupling constant will produce an even larger discrepancy. Results are presented as the ratio of the cross sections in nucleus-nucleus (AA) collisions to the ones in nucleon-nucleon collisions scaled with the number of binary nucleon-nucleon interactions (2.20) In Figure 2 we first show the transverse momentum dependence of the of J/ψ (yellow band) and ψ(2S) (cyan band) suppression. We use minimum bias Pb+Pb collisions at √ s N N = 5.02 TeV for illustration and the suppression is calculated as a sum over centrality classes i corresponding to mean impact parameters b i with weights W i [7] R min. bias We find that the theoretical calculation produces a rather flat transverse momentum dependence of the quarkonium suppression factor R AA . The magnitude of this suppression is large, a factor 3 to 5, and is very similar between the J/ψ and ψ(2S) states. This is easy to understand, as in the parton energy loss picture the nuclear modification depends on the flavor and mass of the propagating parton, the fragmentation functions and the steepness of paticle spectra. The ground and excited J/ψ states have very similar partonic origin and fragmentation functions. The ψ(2S) spectra are slightly harder than the ones for the J/ψ and this accounts for the slightly smaller suppression. Comparison of theoretical calculations to ATLAS experimental data on the transverse momentum dependence of J/Ψ attenuation from √ s N N = 5.02 TeV Pb+Pb collisions at the LHC [40] is presented in Figure 3. The top panel shows results for 0-10% central collisions. As can be seen from the figure, the data is not described by the theoretical predictions. Energy loss calculations overpredict the suppression of J/ψ even in the lowest transverse momentum bin around p T ∼ 10 GeV. At higher transverse momenta the discrepancy is as large as a factor of 3. The bottom panel of Figure 3 shows similar comparison but for minimum bias collisions (ATLAS measurements cover 0-80% centrality). The same conclusion can be reached, i.e. the theoretical calculation predicts significantly the nuclear modification in comparison to the one measured measured by the experiment. Next, we address the relative medium-induced suppression of ψ(2S) to J/ψ in matter in Figure 4. The purple bands correspond to variation of the coupling between the parton and the medium of g = 1.7 − 1.9. Since these are double ratios, the sensitivity to the variation of g is significantly reduced. The upper panel of Figure 4 shows the double nuclear modification ratio as a function of p T compared to CMS data [39]. Theory and experimental measurements are for minimum bias collisions and are clearly very different. The energy loss model predicts slightly smaller suppression for the ψ(2S) state when compared to J/ψ and the double ratio is 10-20% above unity. In contrast, experimental results show that the suppression of the weakly bound ψ(2S) is 2 to 3 times larger than that of J/ψ. It is clear that the energy loss model is incompatible with the hierarchy of excited to ground state suppression of quarkonia in matter. The bottom panel of Figure 4 shows the same ratio as a function of the number of participants N part. and the transverse momenta are integrated in the range of 9-40 GeV. Similar conclusion about the tension between data and the theoretical model calculations can be reached, which is inherent to the model and cannot be resolved by varying the coupling between the partons that fragment into quarkonia and the medium.
In summary, in this section we demonstrated that in the currently accessible transverse momentum range of up to ∼ 50 GeV for quarkonium measurements in heavy ion collisions, the energy loss approach combined with leading power factorization is not compatible with existing experimental data from the LHC. The tensions are both in the overall magnitude of J/ψ suppression and in the relative suppression of the ψ(2S) to the ground J/ψ. This implies that the quarkonium states coexist with the medium and motivates us to pursue the formulation a general theory for quarkonium interactions with nuclear matter.  Figure 4. The double ratio of ψ(2S) to J/ψ suppression (purple bands) as a measure of the relative significance of QCD matter effects on ground and excited states is compared to energy loss model calculations. Upper panel: comparison between theory and CMS data [39] as a function of transverse momentum p T for minimum bias collisions. Lower panel: comparison between theory and ATLAS data [40] as a function of centrality integrated in the p T region of 9-40 GeV.

Toward a formulation of NRQCD G : the Glauber and Coulomb regions
The main goal of this work is to devise a framework where quarkonia propagate in a variety of strongly-interacting media, such as cold nuclear matter, QGP, or a hadron gas. We are interested in the regime where matter itself might be non-perturbative, but the interaction with its quasiparticles is mediated by gluon fields and can be described by perturbation theory. Such approach has proven to be extremely successful in constructing theories of light flavor, heavy flavor, and jet production in heavy ion collisions.
When an energetic particle propagates in matter, the interaction with the quasiparticles of the medium is typically mediated by t−channel exchanges of off-shell gluons, called Glauber gluons. We will, thus, call the new effective theory NRQCD with Glauber gluons, or NRQCD G . We have noticed in the past [33] that when the sources of interaction do not have large momentum component, the exchange gluon field's momentum can scale as soft. Here, we call them Coulomb gluons and treat this limit explicitly. The Lagrangian of NRQCD G is constructed by adding to the vNRQCD Lagrangian the additional terms that include the interactions with quark and gluon sources through (virtual) Glauber/Coulomb gluons exchanges. We may then write, where the effective fields A µ,a G/C incorporate the information about the source fields. In order to extract the form and perform the power-counting of the terms in L Q−G/C (ψ, A µ,a G/C ) we will follow three different approaches: 1. Perform a shift in the gluon field in the NRQCD Lagrangian (A µ us → A µ us + A µ G/C ) and then perform the power-counting established in Table 1 to keep the leading contributions. This approach is also known as the background field method.

A hybrid method, where from the full QCD diagrams for single effective Glauber/Coulomb
gluon insertion, and after performing the corresponding power-counting, one can read the Feynman rules for the relevant interactions.

3.
A matching method where we expand in the power-counting parameter, λ, the full QCD diagrams describing the interactions of an incoming heavy quark and a light quark or a gluon. To get the NRQCD G Lagrangian, we then keep the leading and subleading contributions and focus on the dominant contributions in forward scattering limit. In contrast to the hybrid method, here we also derive the tree level expressions of the effective fields in terms of the QCD ingredients.
The first two methods do not directly involve the source fields, since this information is compressed in the effective fields, A µ,a G/C . We show that the background field method, naively applied in the vNRQCD Lagrangian, yields an ambiguous result. In Appendix A we discuss how to properly implement this method in agreement with the other two methods. The fact that all three approaches then give the same Lagrangian is a non-trivial test of our derivation.
We now consider the scaling of the gluon momenta, q µ G/C , for the Glauber and Coulomb regions and the corresponding scaling of the effective gluon fields, A µ G/C . This is done for three types of sources: collinear, soft, and static. We will use the four-component notation (p 0 , p 1 , p 2 , p 3 ) rather than the light-cone coordinates, (p + , p − , p ⊥ ), since this more compatible with the NRQCD formalism. We use n = (0, 0, 1) as the direction of motion of the collinear source.
Note that, for any gluon interacting with the vNRQCD heavy quark, we require q 0 G/C ∼ λ 2 and q i G/C λ such that the heavy quark momenta, both on the left and right of the insertion, scale as (λ 2 , λ), as illustrated in Figure 5. If all of the three-momenta components scale as λ, i.e. q µ C ∼ (λ 2 , λ) then this corresponds to Coulomb (or potential) gluons. The exchange of such modes between the heavy quarks and soft particle has already been investigated up to next-to-next-leading order in the non-relativistic limit in vNRQCD [41,42]. We compare our derivations with theirs in Section 4.4. On the other hand, collinear particles cannot interact with the heavy quarks through the exchange of Coulomb gluons since this will push the collinear particles away from their canonical angular scaling. The relevant mode here is the Glauber gluons, which scale as q µ G ∼ (λ 2 , λ, λ, λ 2 ). We will, therefore, consider Coulomb gluons for the interaction of the heavy quarks with soft and static modes and Glauber gluons for the interactions with collinear modes: for static and soft sources: q µ C ∼ (λ 2 , λ 1 , λ 1 , λ 1 ) , for collinear sources: q µ G ∼ (λ 2 , λ 1 , λ 1 , λ 2 ) .  We now follow the discussion in Sec. 4.1 of Ref. [33] and [32] to establish the scaling of the gluon fields A µ G and A µ C for the three sources of the virtual gluons. Using Eqs. (4.2) and (4.3) along with the first row of Table 1 in Ref. [33], we establish the scaling shown in Table 1 of this paper. These scalings corresponds to the maximum allowed components for each source. For example Glauber scaling for soft and static sources is also kinematically allowed but the Lagrangian terms resulting from such scaling are power suppressed due to the phase-space integration for the sources.

The background field method
We now proceed with the calculation of the Glauber/Coulomb and heavy quark interactions within the naive background filed method. Here, we shift the ultra-soft gluon fields in the vN-RQCD Lagrangian in Eq.(2.6): A µ,a us → A µ,a us + A µ,a G/C . After this shift, we read the interaction Lagrangian, L Q−G/C , from the leading expansion in λ linear in A µ,a G/C . As mentioned above, this approach is problematic and yields the wrong results. Nonetheless, we proceed with this exercise since it will help us set up the goals of the following section and, in addition, it demonstrates the dangers of not carefully consider the distinction of soft and ultra-soft scales.
We only consider the heavy quark sector, i.e. L Q−G/C , since the antiquark can follow trivially. We will organize the result by powers of λ, where if L For the sub-leading Lagrangian we have contributions only from the collinear and soft sources: where A n = n · A and n is the collinear direction (in our convention n = (0, 0, 1)). Note that, for both L (0) and L (1) , the creation and annihilation of the heavy quark (or antiquark) are not evaluated at the same momenta, i.e. p = p , since momentum is shifted by the Glauber/Coulomb gluon. This suggests that the naive shift of the fields might not yield the correct result due to the ambiguity in the choice of p and p in the Lagrangian L (1) . Indeed, the correct L (1) can be calculated in the non-relativistic limit of QCD with the hybrid and matching methods which we will discuss in the following section. In Appendix A we include a detailed discussion on how to properly implement the background field approach consistent with the power counting procedure. This, then will give results in agreement with the nonrelativistic limit of QCD.

Non-relativistic limit of QCD (NRQCD)
To approach more systematically the inclusion of Glauber/Coulomb gluons in the NRQCD Lagrangian, we begin with some definitions and establishing the notation and conventions we will be using in the rest of this section. We then continue with an exercise to establish some of the terms of the known vNRQCD Lagrangian. This will help us to smoothly transition into the main goal of this analysis, which is introducing the Glauber and Coulomb gluon interaction with the heavy quarks. We will consider the leading and sub-leading corrections to the NRQCD Lagrangian from Glauber and Coulomb gluon exchanges and start with fermonic sources (collinear, static, and soft). We will work in the chiral representation of Dirac matrices.
Then the Dirac spinors in this representation take the following form: 2) The non-relativistic limit of those (|p| p 0 ) is given by where the ellipsis denotes terms of higher order in |p|/p 0 . The normalized rest frame spinors u (0) and v (0) are given by and satisfy the equations of motion with v µ = (1, 0).

Interactions with ultra-soft gluons
In this subsection we will show how one can reconstruct the tree-level NRQCD Lagrangian involving single ultra soft gluon interactions with the heavy quarks. In this exercise we will build the formalism and all ingredients necessary to introduce the Glauber and Coulomb gluon interactions. We do that by studying the non-relativistic limit of the expectation value of the QCD operator O 1 , We will consider the single particle expectation value of the operator O 1 for extracting the kinematic terms in the NRQCD Lagrangian, where we interpret the RHS of the above diagrammatic equation as the corresponding terms generated by the non-relativistic version of O 1 . Similarly, for the interaction terms we then consider an expectation value where the initial state contains an additional gluon. This corresponds to, QCD(λ 1) = . (4.8) In principle, in the above equation we need to consider insertions from the QCD Lagrangian in the LHS and the corresponding NRQCD contributions in the RHS. Its easy to demonstrate that, including those terms and after some simplifications, the result reduced to the same equation as above.
We start with the kinematic terms in Eq. (4.7).
The RHS of Eq. (4.9) vanishes from the equation of motion (EoM), but instead of applying EoM, we will first take the non-relativistic limit which will give the corresponding EoM for the non-relativistic heavy quark (i.e. Schrödinger's equation for free particles). To better understand this statement, imagine a function f (λ) that depends on a small parameter λ.
If the function vanishes for all values of 1 > λ > 0, then if we expand in powers of λ the coefficients have to vanish independently. In the context of NRQCD, λ is the velocity of the heavy quark and we are interested in the leading non-trivial coefficient. Since all coefficients vanish, by non-trivial we mean that an additional condition needs to be imposed for them to vanish. We then interpret this condition as the equation of motion for the non-relativistic theory. Alternatively, one may add a small offshellness to the momenta p and p using r 0 →r 0 and r 0 →r 0 . Then the first non-vanishing term is what we are after. In Eq. (4.8) we have not yet specified the scaling of the vector field or its momenta. For constructing the vNRQCD Lagrangian we will take this gluon to be ultra-soft, where g(q) is an ultra-soft gluon with momenta q ∼ (λ 2 , λ 2 , λ 2 , λ 2 ). We take the nonrelativistic limit of Eq. (4.9) by expanding up-to the leading correction the spinors, and up-to the subleading propagator. For this, we use the Eqs. (4.3), (4.4), and (2.1). We explicitly show all steps.
• O(λ 0 ): At leading power (LP), we expand all relevant elements only in the leading velocity terms, that is the absolute non-relativistic limit where the heavy quark is at rest: Each of the three terms in curly brackets comes from expanding at leading order one of the following:ū(p ), ( / p − m), and u(p). All three terms vanish independently. We will see later that this is a consequence of what we will define as the equation of odd gammas.
• O(λ 2 ): For the next-to-next-to-leading power (NNLP) expansion we need the O(r 2 /m 2 ) from each ofū(p ), ( / p − m), and u(p) but also contributions from mixed NLP expansion:

(4.13)
To simplify this result we note that the first and last term in the curly brackets of the first line, vanish from application of Eq. (4.5). To simplify the second line we use: (4.14) With these modifications the result significantly simplifies to give a familiar expression, 2Q need to vanish, then r 0 = r 2 /2m, which is exactly the well-known non-relativistic relation between the kinetic energy and the three-momenta.
• O(λ 3 ): All terms that contribute to this order can easily be shown to have one or three γ i squeezed between the (u (0) ) † and u (0) . This means that all of them vanish. This statement can be generalized to any odd power, n, of γ i : For future reference we will refer to the above equation as the equation of odd gammas. Thus: In order to account for the O(λ 3 ) terms that come for the decomposition of soft and ultra-soft (see in Eq. (2.3)), we need to make replacements as described in Eq. (2.4). This will give for the leading and subleading contributions, We can now write the Lagrangian that would generate such term, We kept the term proportional to ∂ 2 even though is of higher order (O(λ 4 )) than what we are considering here. This will later help us write the final Lagrangian in a gauge invariant form. In the above equation, ψ p (x) is the two-component Pauli spinor that satisfy the twocomponent Schrödinger's equation: We now turn to the V 2Q,A . Since A µ U is an ultra-soft gluon we have, and thus our expansion of V 2Q,A starts from O(λ 2 ), compared to V 2Q .
• O(λ 2 ): This result, we can trivially get from the LP expansion ofū(p) and u(p).
Then from the equation of odd gammas we have (4.23) • O(λ 3 ): We would like to utilize the result we get in this section later, when we extent to Glauber and Coulomb regions instead of ultra-soft. For this reason, we work with generic three-momenta and we will implement the momentum conservation delta function at the end, Again, from the equation of odd gammas only the µ = k = {1, 2, 3} will contribute to this result Using the momentum conservation delta function and expanding r in its soft and ultra-soft components we get We now have all the ingredients to construct the interaction Lagrangian of NRQCD up-to corrections of O(λ 3 ). Adding the two terms together (4.27) The term 2r us +q is of O(λ 4 ) but we keep it anyway because will help to write the Lagrangian in a gauge invariant form. We, thus, have Therefore, for the total Lagrangian we obtain where we have introduced an O(λ 4 ) term, quadratic in the vector field A, such that we can write the Lagrangian in a gauge invariant form. The interaction terms we constructed here involve only a single gluon vertex. Larger number of gluons contribute only at O(λ 4 ) and higher. For example, from conservation of momentum the difference of the three momentum of the in and out heavy quark is simply the ultra-soft component of the gluon. Of course, up-to the order we are working here this contribution is not relevant, but if we have kept this term we would have, This corresponds to a term in the Lagrangian of the form which is the abelian part of the chomomagnetic operator B i = ijk G jk /2. The complete chromo-magnetic operator contains also a non-abelian part with two gluon fields which we do not reproduce here, but they can be introduced through gauge completion. Alternatively, one can explicitly calculate the contribution of the terms quadratic in the vector field by evaluating the following: where is understood that in the RHS the first term is to be evaluated in the non-relativistic limit. The subtraction of the NRQCD diagram is necessary to avoid double counting. We will no further pursue this analysis here.

Introducing the Glauber and Coulomb interactions
Here we introduce the Glauber/Coulomb interactions by repeating the analysis of expanding in λ the O 1 expectation value V 2Q,A , but this time assuming Glauber/Coulomb gluon scaling instead of ultra-soft. This approach we refer to as hybrid method. The relevant scalings that control the power-counting expansion are then given by Eq. (3.2) and Table 1. To simplify the discussion we will utilize many of the results from the last subsection.
• L (0) : We can use the results from Eqs. (4.22) and (4.23) and directly get: 2Q,A from the last line of Eq. (4.26) and, performing the proper power-counting for q, we have: Since the components A i G/C for i = 1, 2, 3 have different scaling for each source, in order to continue we need to specify the source of the Glauber/Coulomb gluon. Collinear: Static: 2Q,A C = 0 , (4.36) Soft: We are now ready to write the leading and subleading correction to the NRQCD G Lagrangian in the heavy quark sector from virtual (Glauber/Coulomb) gluon insertions, i.e. L Q−G , : and Q−C (ψ, A µ,a C ) = 0 (static) , where we use squared brackets in order to denote the region in which the label momentum operator, P µ , acts. Eqs.
Q−G/C , we find that for the cases of collinear and soft sources there are additional terms that appeared in the hybrid method. We further discuss the origin of the discrepancy in Appendix A.

Matching from QCD including source fields
Here, we will reproduce the results in Eqs. (4.38) and (4.39) by considering the non-relativistic limit of the t-channel diagram for a particular source. We consider both quark and gluon sources. This will give the fields A G and A C , appearing in Eqs. (4.38) and (4.39), as a function of the source currents. We begin with the collinear quark source where p n and p n are the momenta of the incoming and outgoing collinear quarks, respectively, and p and p are the momenta of the corresponding heavy quarks. Taking the collinear limit for the spinor u(p n ) and the non-relativistic limit for u(p) we get We then interpret this term as a Feynman diagram generated by the following Lagrangian: (4.42) In the above equation n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1). This is exactly the result we obtained in Eq. (4.38), but now we have an expression for the background Glauber gluon as a function of the source fields. For the next order result, t coll. , we will keep the expansion of the collinear sector up-to the leading accuracy end expand the heavy quark spinors one order higher in the non-relativistic limit. For that we can utilize the result of Eq. (4.25) to write: Tū n (p n )(gT a ) / n 2 u n (p n ) .
(4.44) This is the result we get using he Lagrangian terms L (1) Q−G given in Eq.(4.39), with A µ,a G given by Eq.(4.42). Since the non-relativistic expansion of the heavy spinors is independent of the sources, it is easy to extent this result for soft and static sources by simply performing the following replacements: With these substitutions, and using the expansion in Eq. (4.43), we find for the t-channel diagram with soft fermion source: (4.46) and with static fermion source, Is easy now to see how these terms for t (0) and t (1) are reproducing exactly the Lagrangian terms in Eqs. (4.38) and (4.39) with for a static source and for a soft source, where soft fermion fields φ are the same that appear in the vNRQCD Lagrangian in Eq.(2.6), and h v, are the heavy fermion field and its properties are governed by the HQET Lagrangian [75,76]. Next, we consider gluon field sources. In this case, in addition to the t-channel diagram we have additional two diagrams that contribute to the same process. These two diagrams correspond to absorbing and re-emitting a collinear (or soft) gluon and are necessary to establish a full gauge invariant result when considering all polarizations of the propagating gluons. As before, we begin with the analysis of collinear sources, Using the following power counting for the light-cone components (along the n µ direction) of the collinear fields, A a,µ n = (A +,a n , A −,a n , A a n⊥ ) ∼ (λ 2 , 1, λ) , (4.51) we expanding the spinors and the heavy quark propagators in the power-counting parameter λ to get for the leading contribution: (4.53) The gluon building block B n⊥ is only the leading term in the strong coupling expansion of the gauge invariant operator Written in terms of the effective Lagrangian, we have where Note that the form of the Lagrangian in terms of the effective Glauber field, A a,µ G , remains the same as in Eqs. (4.42) and (4.38). In the next-to-leading power expansion for the sum of all three diagrams we get This gives where the Glauber field, A a,µ G is given by Eq. (4.56). Comparing with the results for collinear quark sources we find that the Lagrangian in terms of the effective field A µ,a G is identical whichever collinear source (quark vs gluons) we are considering.
Repeating the same exercise for soft gluons, where we replace: p n → p s and p n → p s in Eq.(4.50), we find In the forward scattering limit (q → 0) this result can be further simplified and the corresponding Lagrangian, L Q−C (ψ, A µ,a C ), in terms of the Coulomb field, A µ,a C , can be written in the form of Eqs. (4.38) and (4.39) where the effective Coulomb field in terms of the source soft gluon can be written as follows,

Comparison with the literature
The interaction of heavy quarks with soft fermions and gluons was studied in the framework of vNRQCD in Refs. [41,42]. Here, we are interested in the case where the fields are sourced by partons originating from a quark-gluon plasma (or some other medium), but the formalism (non-relativistic expansion) up-to the effective coupling remains the same. Therefore, we test our approach be comparing our result in Eq. (4.46) for soft fermion sources with those of Eqs. (2.9), (2.10), and (3.11) of Ref. [42] and find that the two agree. Note the overall i factor from expanding the action, also in our notation q = r s − r s . For interactions of the heavy quarks with soft gluons, one should then compare our Eq. (4.59) with Eqs. (3.6), (3.7), and (3.11) of Ref. [42]. Again, the two results are in agreement and we note also the factor of 1/2 introduced at the level of the Lagrangian for the symmetry of exchanging the two soft gluons.
The interactions of heavy quarks with collinear partons were studied in the context of SCET G in Ref. [33], where only the leading Lagrangian, L (0) Q−G , was investigated. For interactions with collinear quarks our result in Eq. (4.41) agrees with the equivalent result in Eq. (4.14) of Ref. [33]. In contrast, for interactions with collinear gluons our results in Eq. (4.52) disagree with the corresponding of Ref. [33]. The disagreement originates from the fact that in [33] the authors consider only the first of the three diagrams and assume the replacement A µ → B µ n⊥ . For forward scattering processes on the medium quasiparticles to lowest non-trivial order, this is the dominant diagram and the gauge invariance of the splitting kernels was checked explicitly by comparing three different gauges: covariant, lightcone, and hybrid. For the general cause, however, we expect that this will not be true. Here, we establish gauge invariance most generally at the level of the matching procedure. Furthermore, to our knowledge the results for L

Conclusions
In recent years, different phenomenological approaches have been proposed to describe the modification of the production cross sections of moderate and high transverse momentum quarkonia. Theoretical guidance on the relative significance of the various nuclear effects in the currently accessible transverse momentum range can be very useful. In this paper we used the leading power factorization limit of NRQCD, along with recent extractions of the LDMEs, to implement the energy loss approach to quarkonium production. We calculated the J/ψ and ψ(2S) suppression in the p T = 10 − 40 GeV range and compared the theoretical predictions to experimental measurements from ATLAS and CMS collaborations at √ s = 5.02 GeV for Pb-Pb collisions. We found that theoretical predictions overestimate of the J/ψ suppression for both 0-10% and 0-80% central collisions and the discrepancies persist even after taking the effective coupling to be smaller than traditionally used for in-medium jet propagation. Most importantly, comparing the double radio R AA [ψ(2S)]/R AA [J/ψ] to data, we also find a disagreement that cannot be resolved within the energy loss model. Wwhile the data show that suppression of exited states is clearly larger by more than a factor of two, the theoretical prediction yields a distinctly opposite trend, suppression of the J/ψ is slightly larger. The strong tension between experimental data and theoretical predictions suggests that the energy loss assumption for production and propagation of quarkonium states in medium needs to be revisited. As a formal step in that direction, we introduced a modified theory of non-relativistic QCD that accounts for the interactions of heavy quarks and antiquarks with the medium through soft-virtual gluon exchanges. We refer to the resulting effective theory as NRQCD G and considered three types of medium sources for the virtual gluons: static, soft, and collinear. For static and soft sources we identified the Coulomb region, q µ C ∼ (λ 2 , λ, λ, λ), to be the most relevant. On the other hand, for collinear sources the leading contributions come from the Glauber region, q µ G ∼ (λ 2 , λ, λ, λ 2 ). We derived the NRQCD G leading and sub-leading Lagrangians for a single virtual gluon exchange. To accomplish this task, we used three different approaches: i) the background field method, ii) a matching (with QCD) procedure, and iii) a hybrid method. Although we found that applying the background field method requires caution in the order of shifting the fields and applying power-counting (as discussed in Section 3.1 and Appendix A), all three methods give the same Lagrangian which serves as a non-trivial test of our derivation. A natural extension of this work will be to also extract the double virtual gluon interactions. This can be achieved with minimal effort in the background field method, as described in Appendix A, but a consistency check through one of the other two approaches is advisable. We have outlined the process of such derivation in the hybrid model below Eq. (4.32).
As we focused on the formal aspects of of NRQCD G , phenomenological applications to various topics of interest are left for the future. In particular, would be interesting to investigate using the EFT derived in this work the modification of the heavy quark-antiquark potential due to medium interactions, which in the vacuum is Coulomb-like. In addition, interactions with the medium could induce radial excitations which will likely cause transitions from one quarkonium state to another. Medium-induced transitions from and to exited states might modify the observed relative suppression rates. Moreover, it is interesting to entertain the possibility of using the terms from the matching procedure to investigate the effect of Glauber gluons in quarkonium production and decay factorization theorems in the vacuum.
following replacements where it is understood that after the replacement the partial derivatives act only on the conjugate of ultra-soft momenta. The four-momentum version of the label momentum operator is defined as P µ = (0, −P). In order to perform the analysis in an organized manner is important to establish the power-counting of the various operator that appear in the Lagrangian. We will conciser each source separately. We start with the collinear source.
We now have all the ingredients to expand the Lagrangian up to O(λ 3 ). 1 Collecting all the terms that do not involve the field A G will give us the heavy quark part of the vNRQCD Lagrangian. For L Q−G we need to collect all the terms that contain at least one power of A G . We, thus, get: +O(λ 4 ) .

(A.4)
Collecting all the terms that involve the field A C we get: Again, this is exactly what we can derive by summing the leading and sub-leading terms, i.e. L We are now ready for the final derivation of this appendix. We implement the same analysis as above for a soft source. Then we have Collecting all the terms that involve the field A C we get: (A.7) The sum of the leading and sub-leading terms, i.e. L As mentioned earlier these terms can be reproduced in the hybrid method or within the matching procedure by evaluating Eq. (4.32). We do not pursue this task here.
It is important to mention that in this section we only study the tree-level result for the NRQCD G Lagrangian. The coefficients for the various terms in the Lagrangian take loop corrections and the coefficients can be written as an expansion in the strong coupling. Logarithmic enchantments in the perturbative expansion of the coefficients need to be resumed through renormalization group equations (RGEs). An important question is if the perturbative expansion for these coefficients remains the same as in NRQCD. Using the background filed approach the coefficients, by construction, do not change. This is, obviously, a very nontrivial statement if using a matching approach. Further studies of this issue might be an important task for the future.