Proof of NRQCD Factorization at All Order in Coupling Constant in Heavy Quarkonium Production

Recently the proof of factorization in heavy quarkonium production in NRQCD color octet mechanism is given at next-to-next-to-leading order (NNLO) in coupling constant by using diagrammatic method of QCD. In this paper we prove factorization in heavy quarkonium production in NRQCD color octet mechanism at all order in coupling constant by using path integral method of QCD. Our proof is valid to all powers in the heavy quark relative velocity. We find that the gauge invariance and the factorization at all order in coupling constant require gauge-completed non-perturbative NRQCD matrix elements that were introduced previously to prove factorization at NNLO.


I. INTRODUCTION
In the last two decades, the NRQCD color octet mechanism [1] for heavy quarkonium production has been very successful in explaining experimental data at high energy colliders such as at Tevatron [2] and at LHC [3]. In its original formulation [1] the proof of factorization in heavy quarkonium production in NRQCD color octet mechanism was lacking. The proof of factorization is an essential requirement to study heavy quarkonium production at high energy colliders. Factorization refers to separation of short-distance effects from long-distance effects in quantum field theory.
Recently the proof of factorization in heavy quarkonium production in NRQCD color octet mechanism is given at next-to-next to leading order (NLLO) in coupling constant by using diagrammatic method of QCD [4][5][6]. However, the proof of factorization in heavy quarkonium production in NRQCD color octet mechanism at all order in coupling constant is still missing. In this paper we will prove factorization in heavy quarkonium production in NRQCD color octet mechanism at all order in coupling constant by using path integral method of QCD.
The typical non-perturbative NRQCD matrix element in heavy quarkonium production is given by [1] < 0|O n |0 >=< 0|χ † (0)K n ξ(0)(a † H a H )ξ † (0)K ′ n χ(0)|0 > (1) where ξ is the two component Dirac spinor field that annihilates a heavy quark, χ is the two component Dirac spinor field that creates a heavy quark, a † H is the operator that creates the heavy quarkonium H in the out state. The factors K n and K ′ n are products of a color matrix (either a unit matrix or T a ), a spin matrix (either a unit matrix or σ i ) and a polynomial of covariant derivative D. The color and spin indices on the fields χ and ξ have been suppressed.
The production cross section for heavy quarkonium H at transverse momentum P T in NRQCD factorizes into a sum of perturbative functions times universal matrix elements, where each NRQCD non-perturbative matrix element < O n > represents the probability of a heavy quark-antiquark pair in state [n], such as color singlet or color octet etc., to produce heavy quarkonium state H.
The fragmentation function for parton i to evolve into a heavy quarkonium at large P T is factorized according to [7] D H/i (z, m c , µ) = n d i→QQ [n] (z, m c , µ) < O n > (3) in terms of same NRQCD non-perturbative matrix elements, along with perturbative functions d i→QQ [n] (z, m c , µ) that describe the evolution of an off-shell parton into a heavy quarkantiquark pair in state [n], such as color singlet or color octet etc..
At a first glance it can be easily seen that the non-perturbative NRQCD matrix element in eq. (1) is not gauge invariant unless it is a color singlet S-wave non-perturbative matrix element. Hence one expects that any non-canceling infrared divergences in the perturbative Feynman diagrams of heavy quark-antiquark production short-distance coefficient can not be factorized in to the definition of the non-perturbative NRQCD matrix element in eq. (1) to study heavy quarkonium production at high energy colliders in the NRQCD color octet mechanism.
This is explicitly shown in [4][5][6] where the NNLO coupling constant calculation shows that the above non-perturbative NRQCD matrix element in eq. (1) is not consistent with factorization of infrared divergences unless it is a color singlet S-wave non-perturbative matrix element. By using the calculation at NNLO in coupling constant and to all powers in the heavy quark relative velocity it was shown in [4][5][6] that the octet S-wave non-perturbative NRQCD matrix element which is gauge invariant and is consistent with the factorization of infrared divergences is given by where Φ (A) is the gauge link or the non-abelian phase in the adjoint representation of SU(3), A µa (x) is the gluon field, P is the path ordering and l µ is the light-like four-velocity.
Note that a necessary condition for NRQCD factorization is that the long-distance behavior of the non-perturbative NRQCD matrix element must be independent of the light-like vector l µ . Such a dependence would be inconsistent with NRQCD factorization because the infrared divergences of < O n > must match those of cross sections, in which there is no information on l µ . In [4][5][6] we have verified the l µ independence of the infrared pole at NNLO in coupling constant and to all powers in heavy quark relative velocity.
Since the NRQCD matrix element in eq. (1) is a non-perturbative quantity in QCD it can not be calculated by using perturbative QCD methods. It is well known that a nonperturbative function can not be studied by using perturbative methods, no matter how many orders of perturbation theory is used. Hence path integral formulation (as opposed to diagrammatic method using perturbation theory) is necessary to study properties of nonperturbative quantities in QCD at all order in coupling constant. The only path integral formulation to study factorization of soft and collinear divergences at all order in coupling constant in quantum field theory is given by R. Tucci in [8] which is exact for QED but is not exact for QCD. We have extended this path integral approach to QCD to prove factorization in QCD at all order in coupling constant in [9]. In this paper we will extend this path integral approach to prove NRQCD factorization at all order in coupling constant in heavy quarkonium production. We will prove that the long-distance behavior of the nonperturbative NRQCD matrix element is independent of the light-like vector l µ at all order in coupling constant.
The paper is organized as follows. In section II we briefly describe the lagrangian density in NRQCD and in QCD. In section III we discuss infrared divergences in NRQCD and in QCD. In section IV we include heavy quark in the path integral formulation of QCD.
In section V we describe infrared divergences in NRQCD and the light-like Wilson line in QCD. In section VI we show that the eikonal current of the light-like charge generates pure gauge field in quantum field theory. In section VII we show how the pure gauge field in quantum field theory can be used to describe soft (infrared) divergences. In section VIII we study heavy quark-antiquark non-perturbative matrix element in the presence of light-like Wilson line in QCD. In section IX we prove factorization in heavy quarkonium production in NRQCD color octet mechanism at all order in coupling constant and to all powers in the heavy quark relative velocity. In section X we show that the factorization theorem is a key ingredient in calculation of NRQCD heavy quarkonium production cross section at collider experiments. Section XI contains conclusions.

II. LAGRANGIAN DENSITY IN NRQCD AND IN QCD
The lagrangian density in QCD including heavy quarks is given by [10] where F µνa is the full non-abelian gluon field tensor, ψ l is the Dirac field of the light quark (l = u, d, s), Ψ is the Dirac field of the heavy quark, D µ is the covariant derivative, γ µ is the Dirac matrix, m l is the mass of the light quark and M is the mass of the heavy quark.
In NRQCD an ultraviolet cutoff Λ ∼ M is introduced. The lagrangian density in NRQCD is given by [1] L NRQCD = L light + L heavy + δL (7) where and where D t and D are the time and space components of the covariant derivative D µ and E and B are electric and magnetic components of the gluon field tensor and σ is the Pauli spin matrix. The dimensionless coefficients c 1 , c 2 , c 3 , c 4 etc. in eq. (10) are obtained by matching NRQCD with QCD [1].

III. INFRARED BEHAVIOR IN NRQCD AND IN QCD
Note that in order for the factorization formula to hold in eqs. (2) and (3) the perturbative functions have to be infrared-safe by definition because infrared limit corresponds to longdistance regime [1]. However, as found in [4,5] the NNLO infrared pole contribution to order v 2 is given by which is not zero where v is the relative velocity of the heavy quark-antiquark pair. Eq. (11) is in the rest frame of the heavy quarkonium ( P = 0) where P µ is the four-momentum of the heavy quarkniuom and l µ is the four-velocity of the light-like Wilson line which is fixed to be l µ = δ µ− along the minus light cone direction in [4][5][6]. The presence of non-zero infrared pole in eq. (11) implies that infrared poles will appear in perturbative functions at NNLO and beyond when the factorization is carried out with octet non-perturbative NRQCD matrix element < χ † K n ξ(a † H a H )ξ † K ′ n χ > in the conventional manner as given by eq. (1) in eqs. (2) and (3). On the other hand, when defined according to its gauge-completed form as given by eq. (4) each octet non-perturbative NRQCD matrix element itself generates precisely the same pole terms given in eq. (11) above. This conclusion is valid to all powers in v at NNLO in coupling constant [6]. Thus NRQCD can accommodate these corrections. Hence our main aim in this paper is to prove that eq. (4) is valid at all order in coupling constant.
Note that in NRQCD an ultraviolet cutoff Λ ∼ M is introduced [1]. Hence the ultraviolet (UV) behavior of QCD and NRQCD differ. However, the infrared (IR) behavior of QCD and NRQCD remains same [11]. Hence the infrared behavior in NRQCD can be obtained by studying the corresponding infrared behavior in QCD. Since the matrix element of the type < χ † K n ξ(a † H a H )ξ † K ′ n χ > is the non-perturbative NRQCD matrix element, it is natural to study its infrared behavior at all order in coupling constant by using path integral method.
Hence we will use path integral method of QCD in this paper.

IV. HEAVY QUARKS AND THE PATH INTEGRAL FORMULATION OF QCD
The generating functional in QCD including the heavy quark is given by [10,12] where Q µa is the quantum gluon field, the symbols l = 1, 2, 3 = u, d, s stand for three light quarks u, d, s and the symbol h stands for heavy quark and In eq. (12) theη u ,η d ,η s are external sources for u, d, s quark fields respectively andη h is the external source for the heavy quark field and the term δ(∂µQ µa ) is the derivative of the gauge fixing term under an infinitesimal gauge transformation [10,12] δQ µa = gf abc Q µb ω c + ∂ µ ω a .
Note that the determinant det( δ(∂µQ µa ) ) in eq. (12) can be expressed in terms of path integration over the ghost fields [10]. However, we will directly work with the determinant ) in eq. (12).
For the heavy quark Dirac field Ψ(x), the non-perturbative matrix element of the type if the factors O n and O ′ n are independent of quantum fields where the suppression of the normalization factor Z[0] is understood as it will cancel in the final result (see eq. (133)).

IN QCD
The gauge transformation of the quark field in QCD is given by Hence one finds that the issue of gauge invariance and factorization of infrared divergences in QCD can be simultaneously explained if ω a (x) can be related to the gluon field A µa (x).
Before proceeding to the issue of gauge invariance and the factorization of infrared divergences in QCD let us first discuss the corresponding situation in QED. The gauge transformation of the Dirac field of the electron in QED is given by Hence we can expect to address the issue of gauge invariance and factorization of infrared divergences in QED simultaneously if we can relate the ω(x) to the photon field A µ (x).
In QED the infrared (or soft) divergence arises only from the emission of a photon for which all components of the four-momentum are small. The Eikonal propagator times the Eikonal vertex for a soft photon with momentum k interacting with a light-like electron moving with four momentum p µ is given by [8,9,[13][14][15][16][17][18][19][20][21] e p µ p · k + iǫ = e l µ l · k + iǫ (18) where l µ is the four-velocity of the light-like electron. Note that when we say the "light-like electron" we mean the electron that is traveling at its highest speed which is arbitrarily close to the speed of light (| l| ∼ 1) as it can not travel exactly at speed of light (| l| = 1) because it has finite mass even if the mass of the electron is very small. From eq. (18) we find where the photon field A µ (x) and its Fourier transform A µ (k) are related by From eq. (19) we find where the eikonal current density J µ (x) for the light-like charge e is given by Now consider the corresponding Feynman diagram for the infrared divergences in QED due to exchange of two soft-photons of four-momenta k µ 1 and k µ 2 . The corresponding Eikonal contribution due to two soft-photons exchange is analogously given by Extending this calculation up to infinite number of soft-photons we find that the Eikonal contribution for the infrared divergences due to soft photons exchange with the light-like electron in QED is given by the exponential where l µ is the light-like four velocity of the electron. The Wilson line in QED is given by When A µ (x) = A µ (lλ) as in eq. (24) then one finds from eq. (25) that the light-like Wilson line in QED for infrared divergences is given by [22] e ie x 0 dx µ Aµ(x) = e −ie ∞ 0 dλl·A(x+lλ) e ie ∞ 0 dλl·A(lλ) .
Note that a light-like electron traveling with light-like four-velocity l µ produces U(1) pure gauge potential A µ (x) at all the time-space position x µ except at the position x perpendicular to the direction of motion of the electron ( l· x = 0) at the time of closest approach [14,23,24].
When A µ (x) = A µ (λl) as in eq. (24) we find l · x = λ l · l = λ = 0 which implies that the light-like Wilson line finds the photon field A µ (x) in eq. (24) as the U(1) pure gauge. The

U(1) pure gauge is given by
which gives from eq. (26) the light-like Wilson line in QED for infrared divergences which depends only on end points 0 and x µ but is independent of the path. The path independence can also be found from Stokes theorem because for pure gauge which gives from Stokes theorem where C is a closed path and S is the surface enclosing C. Now considering two different paths L and M with common end points 0 and x µ we find which implies that depends only on end points 0 and x µ but is independent of path which can also be seen from eq. (28). Hence from eq. (28) we find that the abelian phase or the gauge link in QED is given by From eqs. (17) and (33) one expects that the gauge invariance and factorization of infrared divergences in QED can be explained simultaneously.
One can recall that the gauge invariant greens function in QED in the presence of background field A µ (x) was formulated by Schwinger long time ago [25].
When this background field A µ (x) is replaced by the U(1) pure gauge background field as given by eq. (27) then one finds by using the path integral method of QED that [8] which proves the gauge invariance and factorization of infrared divergences in QED simultaneously. In eq. (35) the < 0|ψ(x 2 ) ψ(x 1 )|0 > is the full Green's function in QED and < 0|ψ(x 2 ) ψ(x 1 )|0 > A is the corresponding Green's function in the background field method of QED. This path integral technique is also used in [18] to prove factorization of infrared divergences in non-equilibrium QED.
Hence we find that the gauge invariance and factorization of infrared divergences in QED can be studied by using path integral method of QED in the presence of U(1) pure gauge background field. Therefore one expects that the gauge invariance and factorization of infrared divergences in QCD can be studied by using path integral method of QCD in the presence of SU(3) pure gauge background field. Now let us proceed to QCD. In QCD the infrared (or soft) divergence arises only from the emission of a gluon for which all components of the four-momentum are small. The Eikonal propagator times the Eikonal vertex for a soft gluon with momentum k interacting with a light-like quark moving with four momentum p µ is given by [8,9,[13][14][15][16][17][18][19]21] where l µ is the four-velocity of the light-like quark. Note that when we say the "light-like quark" we mean the quark that is traveling at its highest speed which is arbitrarily close to the speed of light (| l| ∼ 1) as it can not travel exactly at speed of light (| l| = 1) because it has finite mass even if the mass of the light quark is very small. On the other hand the gluon is massless and hence it always travels at speed of light and is exactly light-like. From eq. (36) we find where the gluon field A µa (x) and its Fourier transform A µa (k) are related by Note that a path ordering in QCD is required which can be seen as follows, see also [20].
The Eikonal contribution for the infrared divergence in QCD arising from a single softgluon exchange in Feynman diagram is given by eq. (37). Now consider the corresponding Feynman diagram for the infrared divergences in QCD due to exchange of two soft-gluons of four-momenta k µ 1 and k µ 2 . The corresponding Eikonal contribution due to two soft-gluons exchange is analogously given by where P is the path ordering. Extending this calculation up to infinite number of softgluons we find that the Eikonal contribution for the infrared divergences due to soft gluons exchange with the light-like quark in QCD is given by the path ordered exponential where l µ is the light-like four velocity of the quark. The Wilson line in QCD is given by which is the solution of the equation [26] with initial condition When A µa (x) = A µa (lλ) as in eq. (40) we find from eq. (41) that the light-like Wilson line in QCD for infrared divergences is given by [22] Pe ig x A light-like quark traveling with light-like four-velocity l µ produces SU(3) pure gauge potential A µa (x) at all the time-space position x µ except at the position x perpendicular to the direction of motion of the quark ( l · x = 0) at the time of closest approach [14,23,24].
When A µa (x) = A µa (λl) as in eq. (40) we find l · x = λ l · l = λ = 0 which implies that the light-like Wilson line finds the gluon field A µa (x) in eq. (40) as the SU(3) pure gauge. The SU(3) pure gauge is given by which gives Hence when A µa (x) = A µa (λl) as in eq. (40) we find from eqs. (44) and (46) that the light-like Wilson line in QCD for infrared divergences is given by which depends only on end points 0 and x µ but is independent of the path. The path independence can also be found from the non-abelian Stokes theorem which can be seen as follows. The SU(3) pure gauge in eq. (45) gives Note that from eq. (48) we find the vanishing physical gauge invariant field strength square is the SU(3) pure gauge as given by eq. (45). Hence in classical mechanics the SU(3) pure gauge potential does not have an effect on color charged particle and one expects the effect of exchange of soft gluons to simply vanish. However, in quantum mechanics the situation is a little more complicated, because the gauge potential does have an effect on color charged particle even if it is SU(3) pure gauge potential and hence one should not expect the effect of exchange of soft gluons to simply vanish [14]. This can be verified by studying the non-perturbative matrix element in QCD such as <Ψ( in the presence of SU(3) pure gauge background field.
Using eq. (48) in the non-abelian Stokes theorem [27] we find where C is a closed path and S is the surface enclosing C. Now considering two different paths L and M with common end points 0 and x µ we find from eq. (49) which implies that the light-like Wilson line in QCD depends only on the end points 0 and x µ but is independent of the path which can also be seen from eq. (47). Hence from eq. (47) we find that the non-abelian phase or the gauge link in QCD is given by In the adjoint representation of SU(3) the corresponding path ordered exponential is given by To summarize this, we find that the infrared divergences in the perturbative Feynman It can be mentioned here that in soft collinear effective theory (SCET) [28] it is also necessary to use the idea of background fields [12] to give well defined meaning to several distinct gluon fields [15].
As mentioned earlier, in NRQCD an ultraviolet cutoff Λ ∼ M is introduced [1]. Hence the ultraviolet (UV) behavior of QCD and NRQCD differ. However, the infrared (IR) behavior of QCD and NRQCD remains same [11]. Hence the infrared behavior in NRQCD can be studied by studying the corresponding infrared behavior in QCD. Hence we find that the infrared behavior of the non-perturbative NRQCD matrix element < 0|χ † K n ξ(a † H a H )ξ † K ′ n χ|0 > in eq. (1) can be obtained by studying the infrared behavior of the non-perturbative matrix n are appropriate factors which identify the state of the heavy quark-antiquark system such as the color singlet state or color octet state etc..
Note that a massive color source traveling at speed much less than speed of light can not produce SU(3) pure gauge field [14,23,24]. Hence when one replaces light-like Wilson line with massive Wilson line one expects the factorization of infrared divergences to break down. This is in confirmation with the finding in [29] which used the diagrammatic method of QCD. In case of massive Wilson line in QCD the color transfer occurs and the factorization breaks down.

VI. EIKONAL CURRENT OF THE LIGHT-LIKE CHARGE GENERATES PURE GAUGE FIELD IN QUANTUM FIELD THEORY
In order to study factorization of infrared divergences by using the background field method of QED, the soft photon cloud traversed by the electron is represented by the pure gauge background field A µ (x) [8] due to the presence of the light-like Wilson line, where one represents the quantum photon field by Q µ (x). As mentioned above, in classical mechanics the assertion that the gauge field that is produced by a highly relativistic (light-like) particle is a pure gauge [14,23,24]. One may ask a question if this assertion is correct in quantum field theory. In this section we will show that this assertion is correct in quantum field theory. We will use path integral formulation of the quantum field theory for this purpose.
The generating functional for the gauge field in the quantum field theory in the presence of external source J µ (x) in the path integral formulation is given by where Q µ (x) is the quantum photon field and The effective action S ef f [J] is given by [34] < where where D µν (x − x ′ ) is the photon propagator.
The photon propagator in the coordinate space is given by Using eq. (58) in (57) we find From the continuity equation we have Using eq. (60) in (59) we find First of all, by using the path integral formulation of the quantum field theory we will derive Coulomb's law for static charge. Note that the derivation of the Coulomb's law by using path integral formulation of the quantum field theory is not necessary to prove factorization theorem. We have included it here only to demonstrate the correctness of the prediction of the path integral formulation in quantum field theory which we will use (see below) to show that the eikonal current of the light-like charge generates pure gauge field in quantum field theory.
In order to derive Coulomb's law by using path integral formulation of the quantum field theory we consider two static charges at positions X and X ′ respectively. The current density for this two static charges is given by Using eq. (62) in (61) and neglecting the self energies we find in the time interval t that which gives the (effective) potential energy V ef f [J] of the interaction between two static charges to be which reproduces the Coulomb's law. Hence we have shown that the assertion that a charge at rest generates a Coulomb gauge field is correct in quantum field theory.

B. Effective Lagrangian Density of Light-Like Eikonal Current in Quantum Field
Theory Similarly using the above procedure in quantum field theory we will show that the assertion that a light-like charge generates pure gauge field is correct in quantum field theory.
This can be shown as follows.
The eikonal current density of the charge e with light-like four-velocity l µ is given by eq. (22). By using the path integral formulation of the quantum field theory we find by using eq. (22) in (61) that for light-like eikonal current the effective lagrangian density is given by For light-like four-velocity we have Hence from (65) and (66) we find that for light-like eikonal current the effective lagrangian density is given by at all the time-space position x µ except at the spatial position perpendicular to the motion of the charge ( l · x = 0) at the time of closest approach (x 0 = 0).

Light-Like Eikonal Current in Quantum Field Theory
Similarly by using the above path integral formulation calculation we find from eq. (A5) that the interaction between the (light-like or non-light-like) non-eikonal current and the gauge field generated by the light-like eikonal current gives the effective (interaction) lagrangian density where q µ is the (light-like or non-light-like) four-momentum of non-eikonal current of charge e and l µ is the light-like four-velocity of the eikonal current of charge e.
For light-like eikonal current we find from eqs. (66) and (68) that effective (interaction) lagrangian density due to the interaction between the (light-like or non-light-like) non-eikonal current of four-momentum q µ and the gauge field generated by the light-like eikonal current of four-velocity l µ is given by This is also obvious from eq. (81).

D. Pure Gauge Field Generated By Eikonal Current of Light-Like Charge in Quantum Field Theory
Hence from eqs. (67) and (69) we find that the eikonal current for light-like charge generates pure gauge field in quantum field theory. From eqs. (67) and (69) we find that the assertion that a light-like charge generates a pure gauge field is correct in quantum field theory which is consistent with the corresponding result in classical mechanics [14,23,24].

VII. PURE GAUGE FIELD IN QUANTUM FIELD THEORY DESCRIBES SOFT (INFRARED) DIVERGENCE
In this section we will show how the pure gauge field is used in quantum field theory to describe soft (infrared) divergences. Consider an incoming electron of four momentum q µ and mass m emitting a real photon of four momentum k µ . The corresponding Feynman diagram contribution is given by [30] where we write and From eq. (3.2) of [30] we write the gauge field as is the physical gauge field [corresponding to transverse polarization of the gauge field] and is the pure gauge field [corresponding to longitudinal polarization of the gauge field]. Now using eq. (73) in eq. (70) we find that the total contribution of the Feynman diagram is given by and Hence in the soft photon limit (k 0 , k 1 , k 2 , k 3 ) → 0 we find from the eqs. (70) and (77) that  Hence we find that we do not need to calculate the finite value of the cross section (or the full cross section) [which will require the non-eikonal-line part of the diagram as given by eq.
(80)] to study the relevant infrared divergence behavior. The relevant infrared divergence behavior can be calculated by using eikonal approximation as given by eq. (79).
From eq. (80) we find that We are interested in the infrared divergence behavior due to the presence of the light-like

MENT IN THE PRESENCE OF LIGHT-LIKE WILSON LINE IN QCD
We have seen in section V that the infrared behavior of the non-perturbative NRQCD matrix element < 0|χ † K n ξ(a † H a H )ξ † K ′ n χ|0 > in eq. (1) can be obtained by studying the infrared behavior of the non-perturbative matrix element in QCD of the type Background field method of QCD was originally formulated by 't Hooft [31] and later extended by Klueberg-Stern and Zuber [32,33] and by Abbott [12]. This is an elegant formalism which can be useful to construct gauge invariant non-perturbative green's functions in QCD. This formalism is also useful to study quark and gluon production from classical chromo field [34] via Schwinger mechanism [35], to compute β function in QCD [36], to perform calculations in lattice gauge theories [37] and to study evolution of QCD coupling constant in the presence of chromofield [38].
In the background field method of QCD the generating functional is given by [12,31,32] where Q µa (x) is the quantum gluon field and the gauge fixing term is given by which depends on the background field A µa and We have followed the notations of [12,31,32] and accordingly we have denoted the quantum gluon field by Q µa and the background field by A µa . The determinant det( δG a (Q) δω b ) in eq. (82) can be expressed in terms of path integration over the ghost fields [10,32]. However, we will directly work with the determinant det( δG a (Q) δω b ) in eq. (82). Note that the gauge fixing term 1 2α (G a (Q)) 2 in eq. (82) [where G a (Q) is given by eq. (83)] is invariant for gauge transformation of A a µ : provided one also performs a homogeneous transformation of Q a µ [12,32]: The gauge transformation of background field A a µ as given by eq. (85) along with the homogeneous transformation of Q a µ in eq. (86) gives  [12,32]: gives eq. (87) which leaves − 1 4 F a2 µν [A + Q] invariant in eq. (82). Extending eq. (82) to include heavy quark [by using the lagrangian density from eq. (6)] we find that the generating functional in the background field method of QCD is given by Note that in the absence of external sources a pure gauge can be gauged away from the generating functional. However, in the presence of external sources a pure gauge can not be gauged away from the generating functional. It is useful to remember that, unlike QED [8], field is not easy. The main difficulty is due to the gauge fixing terms which are different in both the cases. While the Lorentz (covariant) gauge fixing term − 1 2α (∂ µ Q µa ) 2 in eq. (12) in QCD is independent of the background field A µa (x), the background field gauge fixing term − 1 2α (G a (Q)) 2 in eq. (90) in the background field method of QCD depends on the background field A µa (x) where G a (Q) is given by eq. (83) [12,31,32]. Hence in order to study non-perturbative matrix element in the background field method of QCD in the presence of SU(3) pure gauge background field we proceed as follows.
Hence we obtain eqs. (91), (92) and (93) whether we use the type I transformation or type II transformation. Hence we find that we will obtain the same eq. (120) whether we use the type I transformation or type II transformation.

Note that
in eq. (85) is valid for infinitesimal transformation (ω << 1) which is obtained from the finite equation Simplifying infinite numbers of non-commuting terms we find Hence from eqs. (95), (96) and [23] we find that where M ab (x) is given by eq. (97). Similarly, the equation in eq. (93) is valid for infinitesimal transformation (ω << 1) which is obtained from the finite equation which gives where M ab (x) is given by eq. (97).
Changing the variables of integration from unprimed to primed variables in eq. (91) we This is because a change of variables from unprimed to primed variables does not change the value of the integration.
(45) we find [23] A where M ab (x) is given by eq. (97). By using eqs. (101) and (109) in eq. (108) we find which gives From eq. (111) we find From [23] we find which in the adjoint representation of SU(3) gives (by using T a bc = −if abc ) where M ab (x) is given by eq. (97). Using eq. (115) in (113) we find which gives Since for n × n matrices A and B we have we find from eq. (116) that From eqs. (109) and (101) we find where M ab (x) is given by eq. (97).
Hence we find that eq. (127) is the relation between the generating functional Z[J, η u ,η u , η d ,η d , η s ,η s , η h ,η h ] in QCD and the generating functional Note that in QED the corresponding result is [8,18] when the background field A µ (x) is the U(1) pure gauge field given by A µ (x) = ∂ µ ω(x). Eq. For the heavy quark Dirac field Ψ(x), the non-perturbative matrix element of the type where the suppression of the normalization factor Z[0] is understood as it will cancel in the final result (see eq. (133)).
When the background field A µa (x) is the SU(3) pure gauge as given by eq. (45) we find from eqs. (15), (129), (127), (124) and (122) that if the factors O n and O ′ n are independent of quantum fields where, see eq. (52), Note that the creation operator a † q and annihilation operator a q of the quark are related to the quark field via the equation [39] ψ where color indices are suppressed. Hence one finds that the quark field ψ(x) or Ψ(x) depends on the the creation (annihilation) operator a † q (a q ) of the quark but is independent of the creation (annihilation) operator a † H (a H ) of the hadron. Similarly the gluon field Q µa (x) is independent of the creation (annihilation) operator a † H (a H ) of the hadron. Since a † H a H is independent of ψ(x), Ψ(x), Q µa (x) one can perform exactly the similar steps of the path integral calculation as above to find from eq. (130) that where Φ(x) is given by eq. (131).
Under non-abelian gauge transformation as given by eq. (95) the Wilson line in QCD transforms as Pe ig x f From eqs. (47) and (134) we find which gives from eq. (131) Hence we find that < 0|Ψ( (133) is obtained from the exact generating functional in QCD as given by eq. (12), see eq. (15). Similarly, the non-perturbative matrix element < Note that in [4,5] the proof of factorization is presented at NNLO in coupling constant and to v 2 order in relative velocity of the heavy quark-antiquark pair. This is done by restricting the result to v 2 order by using where p µ 1 is the momentum of the heavy quark, p µ 2 is the momentum of the heavy antiquark, P µ is the total momentum of the heavy quark-antiquark pair and p µ r is the relative momentum of the heavy quark-antiquark pair. In the rest frame of heavy quark-antiquark pair p r = M v where M is the mass of the heavy quark.
Similarly in [6] the proof of factorization is presented at NNLO in coupling constant and to all powers in relative velocity v of the heavy quark-antiquark pair. This is done by obtaining the result for arbitrary p µ 1 and p µ 2 without restricting to order p 2 r . Hence in order to be consistent with the proof of factorization of [6] to all powers of relative velocity v it is necessary to present the final result for arbitrary p µ 1 and p µ 2 . It can be seen that the non-perturbative matrix element < is obtained from the exact generating functional in QCD as given by eq. (12) without putting any restrictions on heavy quark and antiquark momenta, see eq. (15). Similarly, the non-perturbative matrix Note that the non-perturbative matrix element < 0|Ψ( (133) is independent of l µ . Hence all the l µ dependence in Φ(x) defined by eq. (131) in the non-perturbative matrix element < To summarize this, we find that eq. (141), which is found by using path integral method of QCD, is valid at all order in coupling constant and to all powers in heavy quark relative velocity. We have also shown that the long-distance behavior of the non-perturbative NRQCD matrix element is independent of the light-like vector l µ at all order in coupling constant and to all powers in heavy quark relative velocity. The eq. (4), which is found by using diagrammatic method of QCD at NNLO in coupling constant and to all powers in heavy quark relative velocity shows that long-distance behavior of the non-perturbative NRQCD matrix element is independent of the light-like vector l µ at NNLO in coupling constant and to all powers in heavy quark relative velocity. This implies that the gauge invariance and the factorization at all order in coupling constant require gauge-completed octet S-wave non-perturbative NRQCD matrix element that was introduced previously to prove factorization at NNLO.
Hence we find that eq. (4) is valid at all order in coupling constant and to all powers in the heavy quark relative velocity.

X. FACTORIZATION THEOREM IS A KEY INGREDIENT IN CALCULATION OF NRQCD HEAVY QUARKONIUM PRODUCTION CROSS SECTION
As mentioned earlier the definition of the NRQCD heavy quarkonium production matrix element from heavy quark-antiquark pair is a non-perturbative quantity which can not be calculated by using perturbation theory no matter how many orders of perturbation theory is used. From this point of view the path integral formulation (as opposed to diagrammatic methods in perturbation theory) is useful to study the properties of the NRQCD non-perturbative matrix element of heavy quarkonium production at all order in coupling constant. As mentioned earlier the only path integral formulation to study factorization of soft and collinear divergences at all order in coupling constant in quantum field theory available is by R. Tucci [8]. However, the calculation of R. Tucci [8] was exact for QED but was not exact for QCD. We have extended the exact path integral calculation of proof of factorization of R. Tucci in QED [8] to proof of factorization in QCD at all order in coupling constant in [9] and to proof of factorization in NRQCD heavy quarkonium production at all order in coupling constant in the previous section.
In this section we will show how the factorization theorem as given by eqs. (130) and (139) is actually a key ingredient in calculation of NRQCD heavy quarkonium production by using eqs. (2) and (3) where the NRQCD non-perturbative matrix element of heavy quarkonium production in color octet mechanism is given by eq. (141).
Let us prove how the eqs. (130) and (139) are key ingredients to prove eqs. (2) and (3) to calculate the NRQCD heavy quarkonium production in color octet mechanism. Suppose we calculate the cross section of heavy quark-antiquark production in color octet state in the presence of light-like quark (or gluon). Then from eq. (130) we find Hence from eq. (143) we find that the A dependence which arises due to the soft gluon exchanges with light-like quark (or gluon) is factorized and only appears in the gauge-links Φ(x) in the right hand side where Φ(x) is given by eq. (131). Eq. (143) implies that in the cross section for QQ production in color octet state at all order in coupling constant the infrared divergences due to the presence of light-like quark (or gluon) are factorized only to the gauge links Φ(x).
Eq. (143) is the exact extension of eq. (1.6) of [8] of the factorization in QED.
Hence from eq. (143) we find that the non-perturbative matrix element of NRQCD heavy quarkonium production in color octet mechanism which cancels these infrared divergences and is consistent with the factorization theorem is obtained from eq. (139) and is given by eq.
(141). This proves that the factorization theorem as given by eqs. (130) and (139) is actually a key ingredient to prove eqs. (2) and (3) to calculate the NRQCD heavy quarkonium production cross section in color octet mechanism at all order in coupling constant.

XI. CONCLUSIONS
Recently the proof of factorization in heavy quarkonium production in NRQCD color octet mechanism is given at next-to-next-to-leading order (NNLO) in coupling constant by using diagrammatic method of QCD. In this paper we have proved factorization in heavy quarkonium production in NRQCD color octet mechanism at all order in coupling constant by using path integral method of QCD. Our proof is valid to all powers in the heavy quark relative velocity. We have found that the gauge invariance and the factorization at all order in coupling constant require gauge-completed non-perturbative NRQCD matrix elements that were introduced previously to prove factorization at NNLO.