On the power divergence in quasi gluon distribution function
Wei Wang Shuai Zhao
hep-ph hep-lat nucl-th






  INPAC, Shanghai Key Laboratory for Particle Physics
and Cosmology,

MOE Key Laboratory for Particle Physics, Astrophysics and Cosmology,

 School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai, 200240,   China







Recent perturbative calculation of quasi gluon distribution function at
one-loop level shows  the existence of extra  linear ultraviolet divergences in the cut-off scheme.   We employ the auxiliary field approach, and study the renormalization of  gluon operators.  The  non-local gluon operator can mix with new operators under renormalization, and the linear
divergences in quasi distribution function can be into the newly introduced operators. After including the mixing, we find the improved quasi gluon distribution functions  contain only logarithmic
divergences, and thus can be used to extract the gluon distribution in large momentum effective theory.




On the power divergence in quasi gluon distribution function
    Wei Wang[wei.wang@sjtu.edu.cn] and Shuai Zhao[shuai.zhao@sjtu.edu.cn] 
    December 30, 2023
==========================================================================







§ INTRODUCTION




The  high energy behaviour of  a hadron is encoded into  parton
distribution functions (PDFs). Because of the
non-perturbative nature, the PDFs can not be
calculated in QCD perturbation theory. Traditionally, PDFs are extracted using the   experimental data
for the processes in which the factorization theorem
is established <cit.>.
Another  popular non-perturbative
method is to use lattice QCD (LQCD). However, since   PDFs are defined as light-cone correlation matrix elements <cit.>, it is a formidable task to calculate the PDFs on the lattice. Instead only the lowest few moments are explored due to the technical difficulties.




A novel approach, named as large momentum effective theory (LaMET), has
been proposed to extract light-cone observables from lattice
calculation <cit.>. In this approach, one can calculate the corresponding quasi observables on the
lattice,  defined by equal-time correlation matrix
elements. The quasi   and light-cone observables share the same infrared (IR) structures. The difference between them is   IR independent,  and   can be calculated perturbatively. Factorization of the  quasi quark distribution has been examined at one-loop level  and proved to
hold to all orders <cit.>.
Therefore, if the matching
coefficient is calculated in perturbation theory, and the quasi
distribution function is evaluated on the lattice, one can extract
light-cone distribution functions. This approach has been
studied extensively
<cit.>, and some encouraging progresses on LQCD calculation are
reported <cit.>. Some other frameworks, e.g.,
the pseudo-PDFs <cit.> and the “lattice cross section” approach <cit.> are also developed.


The one-loop calculations of matching coefficient  in the cutoff scheme show that, the quasi quark PDF  contains linear power ultraviolet (UV)
divergences <cit.>.  These  divergences arise  from   Wilson-line's
self-energy.   Some suggestions have been made to solve the
problem, e.g., the non-dipolar Wilson line <cit.>,  a mass counter term of Wilson
line <cit.>, as well as the auxiliary field
formalism <cit.>.  All these methods are based on the fact
that the linear power divergences  arise from  Wilson lines' self-energy diagram. The renormalizability
of quasi quark distribution has also been proved to all orders  recently <cit.>. Based on the multiplicative renormalizability of quasi-PDF, one can perform  nonperturbative  renormalization of quasi-PDF in the regularization-invariant  momentum-subtraction (RI/MOM) scheme <cit.> , in which the UV divergence can be removed to all orders by a renormalization factor determined on the lattice <cit.>.



Among various PDFs, the gluon
distribution function is a key quantity especially for the  processes like the production of   Higgs particles. It is also crucial for understanding the spin structure of
hadron. However, compared to the quasi quark distribution, the quasi gluon distribution is  less explored. A first calculation of perturbative coefficient for quasi gluon distribution
is presented in Ref. <cit.>, in which the results are presented in both the   cut-off and dimensional regularization (DR)
schemes. It is shown that  linear divergences exist even in the diagrams without any Wilson line, and thus these linear divergences can not be absorbed into  the renormalization of Wilson line. Hence the  gluon quasi PDF is much more involved.
In this work we focus on   power divergences in quasi gluon
distribution. We will show that, in the formalism of auxiliary field, the linear
divergences can be renormalized  by considering the contribution from possible
operator mixing, as well as the mass
counter term of Wilson line.


The rest of this paper is organized as follows. In
Sec. <ref>, we introduce a slightly modified definition of
quasi gluon distribution. In Sec. <ref>, we present the
perturbative calculation of quasi gluon distribution at one-loop
level. In Sec. <ref>,
the renormalization of linear divergence will be analyzed with auxiliary
field method. Then we propose an improved quasi gluon distribution
and calculate the matching coefficient in Sec. <ref>.
Sec. <ref> is the summary and outlook of this work. In the appendix, we give some detailed calculations of the real diagrams with no Wilson lines involved.



§ DEFINITION OF QUASI GLUON DISTRIBUTION


The light-cone (or conventional) gluon distribution function is
defined by a non-local matrix element of two light-cone separated
gluon strength tensor

    f_g/H(x,μ) = ∫dξ^-/2π x P^+ e^-iξ^- xP^+⟨ P|G^+_ i(ξ^- n_+)
      W(ξ^- n_+,0; L_n_+)G^i +(0) |P⟩.

In  the light-cone
coordinate,  a four-vector a is expressed
as a=(a^+, a^-,a⃗_⊥), with a^+=(a^0+a^z)/√(2) and
a^-=(a^0-a^z)/√(2). n_+=(0, 1, 0, 0) is a light-cone vector.
Similarly, the quasi gluon distribution can be defined by non-local
matrix element which is equal-time separated
<cit.>

    f_g/H(x,P^z) = ∫dz/2π x P^z e^i z x P^z⟨ P|G^z_ i(z n_z)
      W(z n_z,0; L_n_z)G^i z(0) |P⟩.

In quasi PDF we still adopt the Cartesian coordinates, where n^μ_z=(0,0,0,1) is an unite vector along z direction.
In the above definitions, W(x,y;C) denotes a Wilson line along
contour C with two endpoints x and y, the gluon field is in
adjoint representation and L_n is a straight line with direction
vector n. The Wilson line can be
parametrized as

    W(z_1,z_2, L_n)=P exp[-ig∫^1_0 dλ ẋ(λ)· A(x(λ))],

where P denotes that the operator exponential is path ordered. Any point
on L_n can be expressed as
x(λ)=λ̅z_2 n+λ z_1 n, where
λ̅≡ 1-λ, x(0)=z_2 n and x(1)=z_1 n. ẋ^μ(λ)≡ d x^μ(λ)/dλ is the velocity
vector of the contour.

According to Refs. <cit.>, i sums over the transverse directions. In this paper, we adopt
a modified definition, where the Lorentz index of gluon field
strength tensor is summed over all the directions except “3”, i.e.,

    f_g/H(x,μ) = ∫d z/2π x P^z e^i z x P^z⟨ P|G^z_ μ(z n_z)
      W(z n_z,0;L_n_z)G^μ z(0) |P⟩ ,

for μ=0,1,2.

To see that Eq. (<ref>) is a proper definition of quasi gluon
distribution, we begin with the energy-momentum
tensor of gluon field

    T^μν= G^μλ G_λ^ ν-1/4 g^μν
      G^αβG_βα.

The averaged value of T^μν between  a hadron state |P⟩ can
be expressed as

    ⟨ P|T^μν|P⟩=a(μ)(P^μP^ν-1/4 M^2
      g^μν),

in which M is the hadron mass.
This tensor is symmetric and traceless with a scalar function
a(μ) characterizing the magnitude.  If we pick up the “++” component,
i.e. μ=ν=+, we have

    ⟨ P|T^++|P⟩=a(μ)(P^+)^2.

The term proportional to M^2 disappears because g^++=0. On the
other hand, if we take μ=ν=z, then

    ⟨ P|T^zz|P⟩= a(μ)[(P^z)^2-1/4
      g^zzM^2]=a(μ)[(P^z)^2+1/4
      M^2]=a(μ)(P^z)^2(1+O(M^2/(P^z)^2)).

Eqs. (<ref>) and (<ref>) are related to the first
moment of light-cone and quasi gluon distribution respectively,

    ∫ dx  x f_g/H(x,μ)   =1/(P^+)^2⟨
      P|T^++|P⟩=a(μ),
    ∫ dx  x  f_g/H(x,P^z,μ)   =1/(P^z)^2⟨
      P|T^zz+1/4 g^zzG^αβG_αβ|P⟩+ O(M^2/(P^z)^2)=a(μ)+ O(Λ^2_QCD/(P^z)^2, M^2/(P^z)^2).

Under the P^z→∞ limit,  one can obtain  ∫ dx  x
f_g/H(x,μ)=∫ dx  x  f_g/H(x,P^z,μ). This relation can
be generalized to high order moments, and consequently, to the
distribution functions f and f themselves. It indicates that
Eq. (<ref>) is a reasonable definition of quasi gluon
distribution. As a  comparison, if λ is only summed
over transverse directions in Eq. (<ref>), the tensor
will not be Lorentz covariant, and  then the above
arguments will no longer hold.


Actually, one can easily validate that both definitions in Eqs. (<ref>) and (<ref>) will give the same results for the matching at tree-level. In addition,
the non-local operator G^z_ μ(z n_z)  W(z n_z,0;L_n_z)G^μ z(0) can be decomposed into the ± components:

    G^z_ μ(z n_z)  W(z n_z,0;L_n_z)G^μ z(0) = 1/2 [G^+_ μ(z n_z) -G^-_ μ(z n_z)] W(z n_z,0;L_n_z)[G^μ +(0) -G^μ -(0) ].

In the large P^z limit, the - component is suppressed and thus the operator for quasi-PDF recovers the same Lorentz structure with the light-cone PDF. The factor 1/2 arises from the conversion of Euclidean and light-cone coordinates.




§ ONE-LOOP CORRECTIONS TO QUASI GLUON DISTRIBUTION



Based on the definition in Eq. (<ref>), we  will re-calculate the one-loop correction to quasi gluon
distribution. To perform  the matching, one can   replace the hadron state with a single parton state. In this  paper we are interested in the gluon-in-gluon case. The hadron state |P⟩ can be replaced by an on-shell gluon state |g(P)⟩, where P  is the gluon momentum. To regulate the UV divergence, we introduce an  UV cut-off Λ on the transverse momentum. This cut-off corresponds to 1/a in LQCD, with  a being the lattice spacing. The collinear divergence is regulated  by   a small gluon mass  m_g.  We will perform the perturbative calculation  in Feynman gauge.


The Feynman rules are listed in Fig. <ref>. The left two sub-diagrams are from the Abelian part of the gluon field strength tensor, while the right two are from the non-Abelian part.  In these Feynman rules, we have explicitly separated the non-Abelian terms.








We firstly make some remarks on the self-energy corrections of the external
gluon. It is known for a long time that the UV cut-off scheme (on all the
components of internal momentum) leads to a quadratic divergence in
QCD vacuum polarization diagram, which breaks the Ward identity and
the gauge invariance.   An exception is the lattice perturbation
theory that provides  an UV cut-off scheme which preserve
the gauge invariance.  For a detailed review on lattice perturbation
theory, one can refer  to Ref. <cit.>.
The quadratic divergence in lattice perturbation
theory will  be canceled by  additional tadpole-like
diagrams, as well as the counter term which originates from the measure of
the path integral. These diagrams are shown in
Fig. <ref>, (a) has the  continuum analog and diverge like 1/a^2. This 1/a^2 divergence corresponds to the quadratic divergence in continuum
theory with a naive UV cut-off. The (b), (c) and (d) only exist in lattice theory but are essential to restore the gauge symmetry.   One can set the same renormalization scales (or the
UV cut-offs) for quasi and light-cone gluon distributions, then the
gluon self energy diagram contributes a same factor to light-cone
and quasi distributions, and will not contribute  to the matching coefficient. Based on the above discussions, it is not necessary to include the
gluon's self energy at present.













We start with  the calculation  of one-loop diagrams with no Wilson line. The Feynman diagrams   are shown in Fig. <ref>.
In Fig. <ref>(a), the non-local vertex is from the
 Abelian term of field strength tensor. A direct calculation gives

    xf̃_g/g^(1)(x, P^z, Λ)|_Fig.<ref>(a)   =α_s C_A/2π{   (2 x^3-3 x^2+2 x-2) lnx-1/x+2 x^2-5
    x/2+3+Λ/P^z ,     x>1
    
       (2 x^3-3
    x^2+2 x-2)ln(1-x+x^2)m_g^2/4x(1-x)(P^z)^2
         +(2x-1)(4x^4-6x^3+10x^2-5x+2)/2(1-x+x^2)
    +Λ/P^z,     0<x<1
       -(2 x^3-3 x^2+2 x-2) lnx-1/x-2 x^2+5
    x/2-3+Λ/P^z .    x<0.


Fig. <ref>(b, c) and (d) correspond to the contributions from
the non-Abelian term in the field strength tensor. Note that in
Ref. <cit.>, the non-Abelian term is combined with the
Wilson line. In the present work, the contributions
from these two parts are separated. This is more convenient to find   the source of
linear divergence. These diagrams give

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(b)=α_s C_A/4π{   x(x+1)lnx-1/x+2x-1-Λ/P^z,       x>1
       x(x+1)ln(1-x+x^2)m^2_g/4x(1-x)(P^z)^2+2x(x-1)+1-Λ/P^z,       0<x<1
       -x(x+1)lnx-1/x-2x+1-Λ/P^z,       x<0
    .
       xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(c)=xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(b),

and

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(d)=α_s C_A/2π{   1-x+Λ/P^z,       x>1
       x-1+Λ/P^z,       0<x<1
       x-1+Λ/P^z.       x<0
    .


The four-gluon
coupling diagram,  depicted in Fig. <ref>(e), is
zero for light-cone distribution. However this diagram  contributes with a non-zero result
to the quasi distribution. This contribution is   linearly divergent:

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(e)=α_s C_A/2π{   x/2-Λ/P^z,       x>1
       x/2-Λ/P^z,       0<x<1
       -x/2-Λ/P^z.       x<0
    .


Fig. <ref>(f,g) are the virtual correction diagrams. Their contributions are proportional to the
tree-level result δ(1-x):

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(f)=δ(1-x)α_s C_A/4π∫ dy {   (y+1)lny-1/y,       y>1
       (y+1)ln(1-y+y^2)m^2_g/4y(1-y)(P^z)^2,       0<y<1
       -(y+1)lny-1/y,       y<0
    .
       xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(g)=xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(f).

Summing the contributions from Fig. <ref>(a)∼(g)
together, we get the total contribution from Fig. <ref>:

    xf̃_g/g^(1)(x,P^z,Λ)|_Fig.<ref>=   α_s C_A/2π{   (2x^3-2x^2+3x-2)lnx-1/x+2x^2-2x+3,       x>1
       (2x^3-2x^2+3x-2)ln(1-x+x^2)m^2_g/4x(1-x)(P^z)^2
         +8 x^5-12 x^4+21 x^3-15 x^2+8 x-2/2 (1-x+x^2),       0<x<1
       -(2x^3-2x^2+3x-2)lnx-1/x-2x^2+2x-3,       x<0
    .
       +δ(1-x)α_s C_A/2π∫ dy {   (y+1)lny-1/y,       y>1
       (y+1)ln(1-y+y^2)m^2_g/4y(1-y)(P^z)^2,       0<y<1
       -(y+1)lny-1/y.       y<0
    .

This result  is free of linear divergence, and  the linear divergences cancel with each other in Fig. <ref>(a)∼(e). In the appendix, we give more details for this cancellation.










We come to the
diagrams which involve one Wilson line shown in
Figs. <ref> and <ref>. In Fig. <ref>, only
one end of the internal gluon line is connected to the Wilson line,
while  both the ends of internal gluon line are
connected to the Wilson line in Fig. <ref>. A direct calculation shows that

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(a)=α_s C_A/4π{   [x(x+1)/1-xlnx/x-1+1-2x/1-x+Λ/(1-x)P^z]_,   x>1
       [-x(x+1)/1-xlnm^2_g(1-x+x^2)/4x(1-x)(P^z)^2-2x^2-2x+1/1-x+Λ/(1-x)P^z]_+,   0<x<1
       [-x(x+1)/1-xlnx/x-1-1-2x/1-x+Λ/(1-x)P^z]_,   x<0
    .
       xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(b)=xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(a),

and

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(c)=xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>(d)=α_s C_A/2π{   [-1-Λ/(1-x)P^z]_,       x>1
       [1-Λ/(1-x)P^z]_+,       0<x<1
       [1-Λ/(1-x)P^z]_.       x<0
    .

In the above expressions, three generalized functions are introduced to regulate the singular integral <cit.>:

    ∫^∞_1 dx  [f(x)]_ T(x)   =∫^∞_1 dx  f(x)[T(x)-T(1)],
    ∫^1_0 dx  [f(x)]_+ T(x)   =∫^1_0 dx  f(x)[T(x)-T(1)],
    ∫^0_-∞ dx  [f(x)]_ T(x)   =∫^0_-∞ dx  f(x)[T(x)-T(1)].

Here  T(x) is an arbitrary smooth   function. Then,
the total contribution from Fig. <ref> is given as

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>=α_s C_A/2π{   [x(1+x)/1-xlnx/x-1-1/1-x-Λ/(1-x)P^z]_,       x>1
       [-x(1+x)/1-xln(1-x+x^2)m^2_g/4x(1-x)(P^z)^2+1-2x^2/1-x-Λ/(1-x)P^z]_+,       0<x<1
       [-x(1+x)/1-xlnx/x-1+1/1-x-Λ/(1-x)P^z]_.       x<0
    .

One can find that the linear divergences  do not vanish.








At last we calculate the one-loop correction to the self energy of
Wilson line, which is presented by Fig. <ref>. The result is

    xf̃_g/g^(1)(x, P^z,Λ)|_Fig.<ref>=α_s C_A/2π{   [1/1-x+Λ/(1-x)^2P^z]_,       x>1
       [1/x-1+Λ/(1-x)^2P^z]_+,       0<x<1
       [1/x-1+Λ/(1-x)^2P^z]_.       x<0
    .

The linear divergence in Eq. (<ref>) has the similar
structure with the linear divergence in quasi quark distribution,
thus may  be renormalized in the same approach.


At last we note that there are also contributions from crossed  diagrams, which have  been discussed in Ref. <cit.>.
These contributions can be easily derived by using the property f(-x)=- f(x).







§ AUXILIARY FIELD METHOD AND RENORMALIZATION OF LINEAR
DIVERGENCES


The light-cone and quasi parton distribution functions are defined
by the gauge invariant quark bi-linear or gluonium operator. A
Wilson line,  a path ordered exponential of a line integral
over gluon field along a contour, is essential to retaining the
gauge invariance. It has been known for a long time that loop
corrections for Wilson line lead to power divergences
<cit.>. The power
divergences can be summed to an exponential factor exp(-π/a
L(c)) where 1/a is the UV cut-off and L(c) is the length of
Wilson line. Based on this, the linear divergence in quasi quark
distribution are
cured <cit.>, because
all the power divergences in quasi quark distribution are from the
Wilson line's self energy. In gluon quasi PDF, it has been shown in
Sec. <ref> that Wilson line's self energy is not the only
source of linear divergence. Since the Wilson line is a functional
of gluon field, the gluon operator and Wilson line might be
entangled in loop corrections.




 It is convenient to discuss gauge invariant non-local operator with the method of auxiliary
field <cit.>. In this method, the non-local operator
can be converted into local operators by introducing an auxiliary
field. The action of the theory is given by

    S=S_QCD+i∫^λ_0 𝒵̅(d/dλ'+i g ẋ^μ A^a_μ t^a
      )𝒵
      dλ'+ i𝒵̅(0) 𝒵(0),

where S_QCD is the action of QCD. Obviously the Wilson
line can be expressed by the two point Green's function of
𝒵- field
<cit.>

    W(z_1,z_2; C)=⟨𝒵(λ_1)𝒵̅(λ_2)⟩_𝒵,

⟨⋯⟩_𝒵 denotes the expectation value respect to
𝒵 only. Therefore, when 𝒵- field is
integrated out in Eq. (<ref>), all the internal
𝒵 line is just the Wilson line.

With the 𝒵 field technique, the gauge invariant
non-local quark bi-linear operator is then given by
<cit.>

    ψ̅_i(z_2) W_ij(z_2,z_1; C)ψ_j(z_1)=⟨ (ψ̅(z_1)𝒵(λ_1))(𝒵̅(λ_2)ψ(z_2))⟩_𝒵,

and the gauge invariant non-local gluonium operator is

    G^a_μν(z_1)W_ab(z_1,z_2; C)G^b_ρσ(z_2)=⟨(G^a_μν(z_1)𝒵_a(λ_1))(𝒵̅_b(λ_2)G^b_ρσ(z_2))⟩_𝒵.

Now the problem of renormalizing the non-local operator has been
converted into the problem of renormalizing the local operators. The
operator is divided into two local operators,  which can be considered individually.


The operator Ω^(1)_μν=G^a_μν𝒵_a, which is a gauge invariant field strength tensor,
may mix with other operators under renormalization. The allowed mixing
operators are constrained by the gauge symmetry, or more strictly,
the BRST symmetry, as well as the dimensional constraint. It has been demonstrated  that  the Ω^(1)_μν can
only  mix with two operators under renormalization
<cit.>

    Ω^(1)_μν   =G^a_μν𝒵_a,
    Ω^(2)_μν   =Ω^(1)_μαẋ_αẋ_ν/ẋ^2-Ω^(1)_ναẋ_αẋ_μ/ẋ^2,
    Ω^(3)_μν   =|ẋ|^-2(ẋ_μA^a_ν-ẋ_νA_μ^a)(D 𝒵)_a,

with

    (D 𝒵)_a=∂_λ𝒵_a- g f_abc ẋ· A_b 𝒵_c.

The Hermitian adjoint of these operators are

    Ω̅^(1)_μν   =𝒵̅_aG^a_μν,
    Ω̅^(2)_μν   =Ω̅^(1)_μαẋ_αẋ_ν/ẋ^2-Ω̅^(1)_ναẋ_αẋ_μ/ẋ^2,
    Ω̅^(3)_μν   =(-𝒵̅D)_a|ẋ|^-2(ẋ_μA^a_ν-ẋ_νA_μ^a),

with

    (𝒵̅D)_a=-∂_λ𝒵̅_a- g f_abc ẋ· A_b 𝒵̅_c.


After renormalization, the operator Ω^(1)_μν can be replaced by the operators Ω^(i)(r.)_μν(i=1,2,3) as

    Ω^(1)_μν→ c_1 Ω^(1)(r.)_μν+ c_2 Ω^(2)(r.)_μν+ c_3Ω^(3)(r.)_μν,

where the c_i (i=1,2,3) are the mixing parameters, and the superscript (r.) denotes that the 𝒵 field in the operators are renormalized.






However, the operator we are interested in is not the gluonium operator with full Lorentz indices. In the definition of quasi gluon
distribution, the operator should contract with n_z, which is the
unit tangent vector of the Wilson line contour. One can easily find that
Ω^(1)_zμ=Ω^(2)_zμ. Therefore, in the present work, Ω^(1)_μν and Ω^(2)_μν are not distinguishable.













To study the renormalization of operator we should calculate the
tree level matrix element from the operator Ω^(3)_μν.
Note that the operator we are interested in is
Ω^(3)_zμ. Since an arbitrary point on W(z, 0, L_n_z) can be expressed as x(λ)=λ z n_z, we have

    Ω^(3)_zμ   =|ẋ|^-2(ẋ_z
      A^a_μ-ẋ_μA^a_z)(D𝒵)_a=-z/|z|^2 A^a_μ(D 𝒵)_a.

Because the path in the Wilson line is along the z direction, the μ component (μ=0, 1, 2)
of ẋ is zero.



Before calculating the tree-level matrix elements we make some
remarks on the renormalization of 𝒵 field. Under
renormalization, the 𝒵 field gets a renormalized mass
δ m_z. Therefore, after renormalization, the equation of
motion for 𝒵 and 𝒵̅ should be

    ∂_λ𝒵_a- g f_abc ẋ· A_b 𝒵_c+δ m_z |z|𝒵^a=   0,
    ∂_λ𝒵̅_a+ g f_abc ẋ· A_b 𝒵̅_c-δ m_z |z|𝒵̅^a=   0.

Then, the contraction of Ω^(3)(r.)_zμ and external gluon state gives

    ⟨1g(p)|1Ω^(3)(r.)_zμ(z')=-ϵ^*_μ(p)e^i p · z'z/|z|^2
     (D𝒵)_a=δ m_z z/|z|ϵ^*_μ(p) e^i p · z'
      𝒵_a,

where the equation of motion for 𝒵 (Eq. (<ref>)) is used.

The contributions from Ω^(3)(r.)_zμ at tree level are shown in Fig. <ref>.
For Fig. <ref>(a), with the contraction rule presented in Eq. (<ref>), its contribution to quasi gluon distribution reads

    x  f_g/g(x,P^z,Λ)|_Fig.<ref>(a)=   ∫dz/2π P^z e^i x P^z z⟨ g(p)|⟨(G^z_ μ𝒵)(z n_z)
      Ω̅^(3)(r.)μ z(0)⟩_𝒵|g(p)⟩^(0)
    
    =   1/2(N^2_c-1)δ^a a'∑_iϵ_μ, i(p)ϵ^*μ_i(p)∫dz/2π
    P^z e^i x
    P^z z i P^z e^-i P^z z(δ m_z) -z/|z|⟨𝒵_a(1)𝒵̅_a'(0)⟩_𝒵
    
    =   i P^z δ m_z∫dz/2π P^z e^i (x-1) P^z zlim_ϵ→
    0z/|z| e^-ϵ|z|
    
    =   i P^zδ m_z ∫dz/2π P^z e^i(x-1)P^z
    zlim_ϵ→ 0∫^∞_0 dαz e^-αϵ^2-z^2/4α/2√(π)α^3/2
    
    =   i P^zδ m_z ∫dz/2π P^z e^i(x-1)P^z
    zlim_ϵ→ 0∫^∞_0
    dα∫^+∞_-∞ d k^z k^z/iπ
    e^-α((k^z)^2+ϵ^2)(1-e^-i k^z
    z)
    
    =   P^z δ m_z ∫dz/2π P^z e^i(x-1)P^z
    zlim_ϵ→ 0∫^+∞_-∞d k^z/πk^z/(k^z)^2+ϵ^2(1-e^-i k^z z)
    
    =   δ m_z ∫^+∞_-∞d
    k^z/πδ((x-1)P^z)-δ((x-1)P^z-k^z)/k^z
    
    =   δ m_z/π P^z[1/1-x-δ(1-x)∫
    dy1/1-y]
    
    =   δ m_z/π
    P^z{   [1/1-x]_,    x>1
       [1/1-x]_+,    0<x<1
       [1/1-x]_.     x<0
    .

Here the superscript (0) denotes the tree level matrix element.
In the above expression, ⟨𝒵_a𝒵̅_a'⟩_𝒵 is taken to be δ_a a' by recalling Eq. (<ref>) and also the fact that the matrix element is at tree level. It is interesting to notice that the divergence has the same structure   with Eq. (<ref>).
Fig. <ref>(b) gives the same contribution. Thus we have

    x  f_g/g(x, P^z, Λ)|_Fig.<ref>=2δ
      m_z/π{   [1/1-x]_,    x>1
       [1/1-x]_+,    0<x<1
       [1/1-x]_.    x<0.









Now we calculate the contribution from Fig. <ref>.



For a Wilson line in fundamental representation,
it has been shown
that <cit.>

    W(z_1,z_2; C)^(r.)=Z^-1exp(-δ m_z
      C(λ_1,λ_2))W(z_1,z_2; C),

where C denotes the length of the contour C. Z is a renormalization constant.

Similarly, for a Wilson line W(z,0,L_n^z) in the adjoint representation, we have

    W(z n_z,0; L_n_z)^(r.)   =e^-δ m_z |z|W(z n_z,0; L_n_z)
       =1-δ m_z |z|+W^(1)(z n_z,0; L_n_z)+⋯.

Therefore, with the same trick  in Ref. <cit.>,
the leading order contribution from the counter term reads

    ∫dz/2π P^ze^i x P^z z(-δ m_z |z|)⟨
      g(p)|G^z_ μ(z n_z)G^μ z(0)|g(p)⟩^(0)
    
    =   (-δ m_z)∫dz/2π P^ze^i x P^z z |z| (P^z)^2 e^-i
    z P^z
    
    =   (-δ m_z)P^z∫dz/2π e^i (x-1) P^z z|z|
    
    =   -δ m_z P^z∫d
    k^z/πδ((1-x)P^z)-δ(k^z-(1-x)P^z)/(k^z)^2
    
    =   -δ m_z/π P^z[-1/(1-x)^2+δ(1-x)∫
    dy1/(1-y)^2]
    
    =   δ m_z/π
    P^z{   [1/(1-x)^2]_,    x>1
       [1/(1-x)^2]_+,    0<x<1
       [1/(1-x)^2]_.    x<0.


Now we turn back to the one-loop corrections, shown by Figs. <ref> and <ref>. Fig.
<ref> denotes the self-energy of Wilson line, and  the linear divergence in this diagram
 should be canceled by the one in Fig. <ref>. From
Eqs. (<ref>) and (<ref>), one can
determine that

    δ m_z=-α_s C_A/2π(πΛ).


Fig. <ref> is the one-loop correction to the gluon-Wilson line vertex.
In the language of 𝒵 field formalism, it can also been
interpreted as the one-loop corrections to Ω^(1)_zμ.
The linear divergence in Fig. <ref> should be canceled by the one in
Fig. <ref>. From Eqs. (<ref>) and
(<ref>) one can find  that the structures of linear
divergences are identically  the same. It indicates that all these linear
divergences can be renormalized by operator Ω^(3)(r.)_zμ,
with only the mixing parameter c_3 undetermined.  From
Eq. (<ref>) and Eq. (<ref>), we have
c_3=-1/2. Then, at one-loop level, all the linear divergences in quasi gluon distribution can be absorbed by δ m_z, and
the matrix elements of  Ω^(3)(r.)_zμ.


The above calculation has been performed in the “naive” cut-off scheme, which may break  the gauge symmetry at first sight.
A complete analysis relies on the lattice perturbation theory. But  it is viable  to examine the gauge invariance from another viewpoint. In Ref. <cit.>, we have calculated
the matching coefficient in the dimensional regularization scheme. It is interesting to notice that
the differences between the cut-off and DR schemes   only reside in the linearly divergent terms, while the finite terms are the same. It is widely believed that gauge invariance is preserved in DR. In the present work, the linearly divergent terms have been absorbed into the matrix elements of the gauge invariant operators, leaving the finite terms consistent with the DR.
Therefore, it is   likely that   our results are  gauge invariant.






To further validate the above proposal, one needs  to investigate the cancelation of linear divergences for all loops. A complete proof is left for a future study, and here we outline the proof at
two-loop level.
Higher order corrections can be
generated by adding internal lines step by step from lower order
diagrams. The
contribution from a n-loop diagram with k propagators can be generally written as

    ∫ d^4 l_1 d^4 l_2⋯ d^4 l_nN_1 N_2 ⋯ N_j/D_1 D_2⋯
      D_kδ(∑_i=1^n c_i
      l^z_i-(1-x)P^z)
    
    =   ∫ d l^z_1 d l^z_2⋯ d l^z_nδ(∑_i=1^n c_i
      l^z_i-(1-x)P^z) [∫ d^3 l_1 d^3 l_2⋯ d^3 l_n N_1 N_2
    ⋯ N_j/D_1 D_2⋯
      D_k] ,

with d^3 l_i≡ d l^0_i d^2 l_i⊥.
Here N_i and D_i
denote the numerator and denominator
respectively, and  c_i are some constants. In the discussions of Sec. <ref> and Sec. <ref>, the z component of loop momentum is either constrained by the δ function, or left unintegrated.
It has been shown that the linear divergences come from the integration over the “0” and transverse components of loop momenta.
Therefore, we separate  the z-components of the loop momenta l_i and leave them
unintegrated. Then the loop integrals in the bracket [...] are all
three-dimensional (3D) integrals.
We define a few quantities for a general diagram as follows:

  
  * L: The number of loops.
  
  * F: The number of internal fermions.
  
  * G: The number of internal gluons and ghosts.
  
  * V_3g: The number of three-gluon vertices and gluon-ghost vertices.
  
  * V_G: The number of vertices corresponding to the Abelian part of field strength tensor. In addition, no external line should connects to these vertices.

With these quantities, we introduce the superficial degree of divergence D for the 3D integrals, with D=3L-F-2G+V_3g+V_G. The propagator of Wilson line only involves the z-component of momentum, hence will not affect the 3D superficial degree of divergence.

There are three generic ways to generate  two-loop diagrams from the
one-loop diagram.


  *  None of the added internal line touches the Wilson line, which are shown in Fig. <ref>.
This is the same with ordinary Feynman diagrams in QCD. For Fig. <ref>(a, b), adding one gluon line will
introduce two three-gluon vertices and three gluon propagators, i.e., G→ G+3, V_3g→ V_3g+2 and L→ L+1. In this case D will be decreased by 1. A linear divergence will reduce to a logarithm divergence, so there is no need to further subtract the linear divergence. Fig. <ref>(c) contains gluon's self energy, although D is increased by 1, the power divergence must be canceled by additional diagrams due to the gauge invariance, which has been discussed at the beginning of Sec. <ref>.












  *  Only one endpoint of the added gluon line is attached  to
Wilson line. This is shown in Fig. <ref>.  For Fig. <ref>(a), there will be two more
gluon propagators, one three-gluon vertex, and one Wilson line
propagator, i.e., G→ G+2, V_3g→ V_3g+1, and L→ L+1,
then D will not be modified; For Fig. <ref>(b),
there will be two more
gluon propagators, one three-gluon vertex, and one Wilson line
propagator as well, but V_G will be modified from 1 to 2,
then D will be increased by 1. The power divergence in this diagram can be
viewed as the product of two linear divergences, which are from the
one-loop corrections of the gluon-Wilson line (or gluon-𝒵) vertices.
In this case,
the power
divergence can be removed by adding contributions from
Ω^(3)(r.)_μν. For Fig. <ref>(a) just one counter
term is needed, while for (b) we need two counter terms.










  *  Both the endpoints of the added gluon line touch the Wilson
line, which are shown in Fig. <ref>. In this case there will be one more
gluon propagator, i.e., G→ G+1, and two more Wilson line propagators, then D will be increased by 1. This case can be explained as Wilson line's self
energy. The power divergence for this case should be cured by adding
mass counter term δ m_z of the Wilson line, as well as the corresponding counter term for
the linear divergence from the lower order diagram.








The above discussions on two-loop corrections are conjectured to
hold to all orders, though we are lack of a complete all-order proof.

After discussing the linear divergence we add some remarks on the
logarithm divergence.
 From Sec. <ref>  one can find
that there is no logarithm divergence in the one-loop real
corrections.  By employing the auxiliary field approach the
renormalization of gauge invariant quark bi-linear and gluonium
operators has been discussed in DR scheme for a long
time <cit.>. Note that in DR
the linear divergence is not regularized. A remarkable
simplification is achieved when the gluonium operator is contracted
with the tangents of the contour, in this case, the operator will
not be renormalized <cit.>, i.e.,

    ẋ_μ G^μνρσẋ_σ=ẋ_μ (G^μνρσ)^(r.)ẋ_σ

with

    G^μνρσ=G^μν(z_1)W(z_1,z_2;C)G^ρσ(z_2).

This  can be
traced back to the property of functional derivation of Wilson line.
The quasi gluon distribution belongs to this case.

It is also shown that the quasi distributions can   be defined
from other components of the correlators. For example, the quasi
quark distribution can   be defined by the matrix element of
ψ̅(z)γ^0ψ(0) <cit.>. For gluon, one
can pick up any indices of μ and ν in the matrix elements
of G^μρ(z)W(z,0, L_n_z)G_ρ^ ν(0), in present work we choose
μ=ν=z. Eq. (<ref>) indicates that, if one
defines quasi gluon distribution by specifying μ and ν
components of gluonium operator, the renormalization of logarithm divergence
will be greatly simplified when the tangents of the contour are collinear to the
projection vector corresponding to these components; On the other
hand, if one modifies the contour of the Wilson line (e.g., in the
non-dipolar Wilson line approach <cit.>),   it is more convenient to modify the components of
gluon strength tensor correspondingly.







§ MATCHING COEFFICIENT FOR IMPROVED GLUON DISTRIBUTION


With the discussions in Sec. <ref>, we propose the subtracted quasi gluon distribution
as

    x  f^(r.)_g/H(x,P^z,Λ)=∫d z/2π P^ze^i x z P^z⟨ P|⟨(Ω^(1)(r.)z_       μ+c_3 Ω^(3)(r.)z_       μ)(z n_z)(Ω̅^(1)(r.)μ z+c_3 Ω̅^(3)(r.)μ z)(0)⟩_𝒵|P⟩,

where c_3=-1/2+O(α_s) is the mixing parameter determined by perturbation theory. δ m_z and the matrix elements can be evaluated
non-perturbatively in LQCD. According to LaMET, one can expect a
factorization formula as

    x f^(r.)_g/H(x, P^z ,Λ )=∫ ^1_0 dy/yZ_g i(x/y,
      Λ/P^z) y f_i/H(y, Λ).



In perturbation theory the Z_gg can be expanded as

    Z_gg(ξ, Λ/P^z)=∑_i=0^∞(α_s/2π)^n
      Z^(n)_gg(ξ,Λ/P^z),

with Z^(0)_gg(ξ)=δ(1-ξ), ξ≡ x/y.

The one-loop correction to light-cone gluon distribution has already
been calculated in Ref. <cit.>, with result

    x f_g/g^(1)(x, Λ)|_total.=   α_s
     C_A/2π{[(2x^3-3x^2+2x-2)ln(1-x+x^2)m^2_g/Λ^2-3x/2(1-x+x^2)+2x^2(x-1)+7/2 x-2]
       -x^2(1+x)/1-xln(1-x+x^2)m^2_g/Λ^2-δ(1-x)∫^1_0 dy
     1+y/1-ylnΛ^2/(1-y+y^2)m^2_g}
    
     =   α_s C_A/2π{(2x^3-2x^2+4x-1)ln(1-x+x^2)
    m^2_g/Λ^2+[x+1/x-1ln(1-x+x^2)
    m^2_g/Λ^2]_+
       +2x^2(x-1)+7/2x-2-3x/2(1-x+x^2)}

for 0<x<1. The result is  0 for x>1 and x<0. On the other hand, the
one-loop correction to   gluon-in-gluon quasi distribution reads

    x  f_g/g^(r.)(1)(x, P^z,Λ)|_total.=α_s C_A/2π{   (2x^3-2x^2+4x-1)lnx-1/x+[x+1/x-1lnx-1/x]_
       +2x^2-2x+3,       x>1
       (2x^3-2x^2+4x-1)ln(1-x+x^2)m^2_g/4x(1-x)(P^z)^2
       +[x+1/x-1ln(1-x+x^2)m^2_g/4x(1-x)(P^z)^2+2x^2/x-1]_+
       +8 x^5-12 x^4+21 x^3-15 x^2+8
    x-2/2 (1-x+x^2),       0<x<1
       -(2x^3-2x^2+4x-1)lnx-1/x-[x+1/x-1lnx-1/x]_
       -2x^2+2x-3.       x<0
    .


From Eqs. (<ref>) and (<ref>), one can determine the
matching coefficient at one-loop level,

    Z^(1)_gg(ξ, Λ/P^z)=C_A{   (2ξ^3-2ξ^2+4ξ-1)lnξ-1/ξ+[ξ+1/ξ-1lnξ-1/ξ]_+2ξ^2-2ξ+3,       ξ>1
       (2ξ^3-2ξ^2+4ξ-1)lnΛ^2/4ξ(1-ξ)(P^z)^2+[ξ+1/ξ-1lnΛ^2/4ξ(1-ξ)(P^z)^2]_+
       +[2ξ^2/ξ-1]_++2ξ^3+ξ+1,       0<ξ<1
       -(2ξ^3-2ξ^2+4ξ-1)lnξ-1/ξ-[ξ+1/ξ-1lnξ-1/ξ]_-2ξ^2+2ξ-3.       ξ<0
    .

This coefficient has only the logarithmic dependent on the cut-off Λ. It is
free of IR divergence, which indicates that the LaMET factorization
holds for gluon distribution at one-loop level.




§ SUMMARY AND OUTLOOK




In this work we have  discussed   power divergences  in the quasi gluon
distribution.  The gluon PDF is parametrized by two gluon field strength tensors contracted with each other.  From  the viewpoint of covariance, we proposed in the contraction to sum all components  of Lorentz indices instead of only transverse ones. With this definition, we have re-calculated  the gluon-in-gluon quasi distribution function at one-loop level, in which the Abelian and non-Abelian contributions are separated.

We found that power divergences cancel in diagrams without Wilson line,  and not all linear divergences  can be attributed into the renormalization of Wilson line.
Employing auxiliary field method,  we express the Wilson line as the expectation value of the product of two auxiliary fields. We found that a new operator can mix with the gauge invariant gluon field strength tensor, and  pointed out that the power divergence can be absorbed into the matrix elements of this newly introduced operator at one-loop level.  This cancelation of power divergence is conjectured to
hold to all orders.

We note that  since one can identify the 𝒵 field as a Wilson line 𝒵(z)=W(z,∞,L_n_z),  this new operator is a product of  covariant derivatives of Wilson line.
We then have used an improved quasi gluon distribution, and determined the matching
coefficient  at   one-loop level. This matching  coefficient   is IR
finite and free of UV power divergence. Therefore, the present work provides the possibility of extracting gluon distributions from
LaMET and LQCD calculation.

We have to note that it is still a challenge of evaluating gluon PDF on the lattice and the present work leaves substantial room for improvement.


  * In this work, we only consider the quasi gluon PDF in continuum QCD, the non-local gluon operators are constructed by the gluon field strength tensor and Wilson lines are under adjoint representation.  To simulate quasi gluon PDF on the lattice, one should define quasi PDF and perform the matching calculation in lattice gauge theory, in which the gauge field is expressed by the Wilson loops under fundamental representation. For a practical lattice evaluation, one needs to define the gauge invariant non-local gluon operators by the renormalized Wilson loops under fundamental representation. The mixing of the gauge invariant gluon operators should be considered as well.

  * Since the quasi PDF should be evaluated on the lattice while the light-cone PDF is defined in continuum, one has to match the light-cone PDF in continuum to the quasi PDF calculated with the lattice regularization. Furthermore, to arrive at a reliable result, renormalizing the matrix elements on the lattice is also a crucial step. Two approaches have been employed to renormalize the quasi quark PDF, one is the lattice perturbation theory <cit.>, another one is the nonperturbative renormalization approach, e.g., the RI/MOM scheme <cit.>. The renormalized matrix element in RI/MOM is independent of UV regularization, thus one can match the lattice evaluated quasi PDF to the MS renormalized light-cone PDF, and work out the matching coefficient in DR, instead of the naive cut-off. However, applying RI/MOM to gluon operators is not straightforward for quasi gluon PDF, the reason is that the multiplicative renormalizability is not clear for non-local gluon operator, during to the operator mixing discussed in this paper. At the moment, a perturbative renormalization with lattice perturbation theory will be useful in practical extraction of gluon PDF.

These issues will be addressed in the forthcoming work.





§ ACKNOWLEDGMENTS


We  are grateful to  J. W. Chen, Y. Jia, J. P. Ma, and  R. L. Zhu for inspiring discussions, and T. Ishikawa, X. D. Ji,  H. N. Li, J. Xu, Y. B. Yang and J. H. Zhang
for valuable suggestions. W.W. thanks Cai-Dian Lü and Qiang Zhao for their hospitality when this work is finalized.
This work is supported  in part   by National  Natural
Science Foundation of China under Grant
 No.11575110, 11655002, 11735010,  Natural  Science Foundation of Shanghai under Grant  No. 15DZ2272100 and No. 15ZR1423100, Shanghai Key Laboratory for Particle Physics and Cosmology, and  by  MOE  Key Laboratory for Particle Physics, Astrophysics and Cosmology.







§ CANCELATION OF THE LINEAR DIVERGENCES IN NO-WILSON LINE DIAGRAMS


In Sec. <ref>, we have concluded that the linear divergences cancel between the diagrams which with no Wilson line (Fig. <ref>). In the following discussion, we will show that how the cancelation works.

The UV divergence appears when loop momenta go to infinity. Since the z component of loop momentum in one-loop real correction is constrained by the δ function, the UV divergence is from the region k∼ (Λ, Λ, Λ, P^z), where k is the loop momentum and Λ≫ P^z is the largest scale. In this region, one can neglect k^z as well as the external momentum P. Then, under this approximation, k is replaced by k=(k^0, k^x, k^y, 0). To compare with the original definition in Refs. <cit.> and the new definition in Eq. (<ref>), we calculate the correlation matrix element

    ℳ^μν=∫d z/2π P^ze^i x z P^z⟨ g(P)|G^zμ(z n_z)G^ν z(0)|g(P)⟩.

Since we are interested in the no-Wilson diagrams, the Wilson line has been eliminated.

In the region k∼ (Λ, Λ, Λ, P^z), for Fig. <ref>, we have

    ℳ^μν|_Fig.<ref>(a)   ∼ -2 i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^μk^νk^2_⊥/(k^2)^3+O(Λ^0),
    ℳ^μν|_Fig.<ref>(b)   ∼  i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^μ_⊥k^ν/(k^2)^2+O(Λ^0),
    ℳ^μν|_Fig.<ref>(c)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^μk^ν_⊥/(k^2)^2+O(Λ^0),
    ℳ^μν|_Fig.<ref>(d)   ∼ -i g^2 C_A n^2/2P^z∫d^3 k/(2π)^4g^μν_⊥/k^2+O(Λ^0),
    ℳ^μν|_Fig.<ref>(e)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^μk^ν/(k^2)^2+O(Λ^0).


By contracting ℳ^μν with g_μν we arrive at the definition of quasi gluon PDF adopted in this work. Then, from the above results, we have

    xf(x, P^z, Λ) |_Fig.<ref>(a)   ∼ -2 i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0),
    
      xf(x, P^z, Λ)|_Fig.<ref>(b)   ∼  i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0),
    
      xf(x, P^z, Λ)|_Fig.<ref>(c)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0),
    
      xf(x, P^z, Λ)|_Fig.<ref>(d)   ∼ -i g^2 C_A n^2/P^z∫d^3 k/(2π)^41/k^2+O(Λ^0),
    
      xf(x, P^z, Λ)|_Fig.<ref>(e)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^41/k^2+O(Λ^0).

One can find that in the region k∼ (Λ, Λ, Λ, P^z),

    xf(x, P^z, Λ) |_Fig.<ref>(a)+ xf(x, P^z, Λ) |_Fig.<ref>(b)+ xf(x, P^z, Λ) |_Fig.<ref>(c)∼ O(Λ^0),
       xf(x, P^z, Λ) |_Fig.<ref>(d)+ xf(x, P^z, Λ) |_Fig.<ref>(e)∼ O(Λ^0).

Therefore, the linear divergences cancel between the no-Wilson line diagrams under the definition proposed in this work.

On the other hand, if ℳ^μν is contracted with g_⊥μν, we return to the definition proposed by Refs. <cit.>. In the region k∼(Λ, Λ, Λ, P^z), we have

    g_⊥μνℳ^μν|_Fig.<ref>(a)   ∼ -2 i g^2 C_A n^2/P^z∫d^3 k/(2π)^4(k^2_⊥)^2/(k^2)^3+O(Λ^0),
    
      g_⊥μνℳ^μν|_Fig.<ref>(b)   ∼  i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0),
    
      g_⊥μνℳ^μν|_Fig.<ref>(c)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0),
    
      g_⊥μνℳ^μν|_Fig.<ref>(d)   ∼ -i g^2 C_A n^2/P^z∫d^3 k/(2π)^41/k^2+O(Λ^0),
    
      g_⊥μνℳ^μν|_Fig.<ref>(e)   ∼ i g^2 C_A n^2/P^z∫d^3 k/(2π)^4k^2_⊥/(k^2)^2+O(Λ^0).

Totally,

    g_⊥μνℳ^μν|_Fig.<ref>∼ i g^2 C_An^2/P^z∫d^3 k/(2π)^4k^2_0 k^2_⊥-k^2_0 k^2_0/(k^2)^3+O(Λ^0)∼ O(Λ^1),

which indicates that the linear divergence will not cancel between diagrams with no Wilson line under the original definition.

99

Butterworth:2015oua
  J. Butterworth et al.,
  
  J. Phys. G 43, 023001 (2016)
  doi:10.1088/0954-3899/43/2/023001
  [arXiv:1510.03865 [hep-ph]].
  
  




Hou:2016nqm
  T. J. Hou et al.,
  
  Phys. Rev. D 95, no. 3, 034003 (2017)
  doi:10.1103/PhysRevD.95.034003
  [arXiv:1609.07968 [hep-ph]].
  
  




Hou:2017khm
  T. J. Hou et al.,
  
  arXiv:1707.00657 [hep-ph].
  
  




Gao:2017yyd
  J. Gao, L. Harland-Lang and J. Rojo,
  
  arXiv:1709.04922 [hep-ph].
  
  



Collins:2011zzd
  J. Collins,
  
  Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol.  32, 1 (2011).
  
  


Ji:2013dva
  X. Ji,
  
  Phys. Rev. Lett.  110, 262002 (2013)
  doi:10.1103/PhysRevLett.110.262002
  [arXiv:1305.1539 [hep-ph]].
  
  




Ji:2014gla
  X. Ji,
  
  Sci. China Phys. Mech. Astron.  57, 1407 (2014)
  doi:10.1007/s11433-014-5492-3
  [arXiv:1404.6680 [hep-ph]].
  
  




Xiong:2013bka
  X. Xiong, X. Ji, J. H. Zhang and Y. Zhao,
  
  Phys. Rev. D 90, no. 1, 014051 (2014)
  doi:10.1103/PhysRevD.90.014051
  [arXiv:1310.7471 [hep-ph]].
  
  




Ma:2014jla
  Y. Q. Ma and J. W. Qiu,
  
  arXiv:1404.6860 [hep-ph].
  
  




Gamberg:2014zwa
  L. Gamberg, Z. B. Kang, I. Vitev and H. Xing,
  
  Phys. Lett. B 743, 112 (2015)
  doi:10.1016/j.physletb.2015.02.021
  [arXiv:1412.3401 [hep-ph]].
  
  




Ji:2014hxa
  X. Ji, P. Sun, X. Xiong and F. Yuan,
  
  Phys. Rev. D 91, 074009 (2015)
  doi:10.1103/PhysRevD.91.074009
  [arXiv:1405.7640 [hep-ph]].
  
  




Jia:2015pxx
  Y. Jia and X. Xiong,
  
  Phys. Rev. D 94, no. 9, 094005 (2016)
  doi:10.1103/PhysRevD.94.094005
  [arXiv:1511.04430 [hep-ph]].
  
  




Ji:2015qla
  X. Ji, A. Schäfer, X. Xiong and J. H. Zhang,
  
  Phys. Rev. D 92, 014039 (2015)
  doi:10.1103/PhysRevD.92.014039
  [arXiv:1506.00248 [hep-ph]].
  
  




Ji:2015jwa
  X. Ji and J. H. Zhang,
  
  Phys. Rev. D 92, 034006 (2015)
  doi:10.1103/PhysRevD.92.034006
  [arXiv:1505.07699 [hep-ph]].
  
  




Monahan:2015lha
  C. Monahan and K. Orginos,
  
  Phys. Rev. D 91, no. 7, 074513 (2015)
  doi:10.1103/PhysRevD.91.074513
  [arXiv:1501.05348 [hep-lat]].
  
  




Xiong:2015nua
  X. Xiong and J. H. Zhang,
  
  Phys. Rev. D 92, no. 5, 054037 (2015)
  doi:10.1103/PhysRevD.92.054037
  [arXiv:1509.08016 [hep-ph]].
  
  




Gamberg:2015opc
  I. Vitev, L. Gamberg, Z. Kang and H. Xing,
  
  PoS QCDEV 2015, 045 (2015)
  [arXiv:1511.05242 [hep-ph]].
  
  




Ishikawa:2016znu
  T. Ishikawa, Y. Q. Ma, J. W. Qiu and S. Yoshida,
  
  arXiv:1609.02018 [hep-lat].
  
  




Monahan:2016bvm
  C. Monahan and K. Orginos,
  
  JHEP 1703, 116 (2017)
  doi:10.1007/JHEP03(2017)116
  [arXiv:1612.01584 [hep-lat]].
  
  




Bacchetta:2016zjm
  A. Bacchetta, M. Radici, B. Pasquini and X. Xiong,
  
  Phys. Rev. D 95, no. 1, 014036 (2017)
  doi:10.1103/PhysRevD.95.014036
  [arXiv:1608.07638 [hep-ph]].
  
  




Ji:2017rah
  X. Ji, J. H. Zhang and Y. Zhao,
  
  Nucl. Phys. B 924, 366 (2017)
  doi:10.1016/j.nuclphysb.2017.09.001
  [arXiv:1706.07416 [hep-ph]].
  
  




Nam:2017gzm
  S. i. Nam,
  
  Mod. Phys. Lett. A 32, 1750218 (2017)
  doi:10.1142/S0217732317502182
  [arXiv:1704.03824 [hep-ph]].
  
  




Rossi:2017muf
  G. C. Rossi and M. Testa,
  
  Phys. Rev. D 96, no. 1, 014507 (2017)
  doi:10.1103/PhysRevD.96.014507
  [arXiv:1706.04428 [hep-lat]].
  
  




Stewart:2017tvs
  I. W. Stewart and Y. Zhao,
  
  arXiv:1709.04933 [hep-ph].
  
  




Xiong:2017jtn
  X. Xiong, T. Luu and U. G. Meißner,
  
  arXiv:1705.00246 [hep-ph].
  
  




Hobbs:2017xtq
  T. J. Hobbs,
  
  arXiv:1708.05463 [hep-ph].
  
  




Broniowski:2017wbr
  W. Broniowski and E. Ruiz Arriola,
  
  Phys. Lett. B 773, 385 (2017)
  doi:10.1016/j.physletb.2017.08.055
  [arXiv:1707.09588 [hep-ph]].
  
  




Carlson:2017gpk
  C. E. Carlson and M. Freid,
  
  Phys. Rev. D 95, no. 9, 094504 (2017)
  doi:10.1103/PhysRevD.95.094504
  [arXiv:1702.05775 [hep-ph]].
  
  




Chen:2017mie
  J. W. Chen, T. Ishikawa, L. Jin, H. W. Lin, Y. B. Yang, J. H. Zhang and Y. Zhao,
  
  arXiv:1710.01089 [hep-lat].
  
  




Monahan:2017hpu
  C. Monahan,
  
  arXiv:1710.04607 [hep-lat].
  
  




Lin:2014zya
  H. W. Lin, J. W. Chen, S. D. Cohen and X. Ji,
  
  Phys. Rev. D 91, 054510 (2015)
  doi:10.1103/PhysRevD.91.054510
  [arXiv:1402.1462 [hep-ph]].
  
  




Alexandrou:2015rja
  C. Alexandrou, K. Cichy, V. Drach, E. Garcia-Ramos, K. Hadjiyiannakou, K. Jansen, F. Steffens and C. Wiese,
  
  Phys. Rev. D 92, 014502 (2015)
  doi:10.1103/PhysRevD.92.014502
  [arXiv:1504.07455 [hep-lat]].
  
  




Chen:2016utp
  J. W. Chen, S. D. Cohen, X. Ji, H. W. Lin and J. H. Zhang,
  
  Nucl. Phys. B 911, 246 (2016)
  doi:10.1016/j.nuclphysb.2016.07.033
  [arXiv:1603.06664 [hep-ph]].
  
  




Alexandrou:2017huk
  C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos and F. Steffens,
  
  Nucl. Phys. B 923, 394 (2017)
  doi:10.1016/j.nuclphysb.2017.08.012
  [arXiv:1706.00265 [hep-lat]].
  
  




Lin:2017ani
  H. W. Lin, J. W. Chen, T. Ishikawa and J. H. Zhang,
  
  arXiv:1708.05301 [hep-lat].
  
  




Briceno:2017cpo
  R. A. Briceño, M. T. Hansen and C. J. Monahan,
  
  Phys. Rev. D 96, no. 1, 014502 (2017)
  doi:10.1103/PhysRevD.96.014502
  [arXiv:1703.06072 [hep-lat]].
  
  




Chen:2017mzz
  J. W. Chen, T. Ishikawa, L. Jin, H. W. Lin, Y. B. Yang, J. H. Zhang and Y. Zhao,
  
  arXiv:1706.01295 [hep-lat].
  
  




Constantinou:2017sej
  M. Constantinou and H. Panagopoulos,
  
  Phys. Rev. D 96, no. 5, 054506 (2017)
  doi:10.1103/PhysRevD.96.054506
  [arXiv:1705.11193 [hep-lat]].
  
  




Zhang:2017bzy
  J. H. Zhang, J. W. Chen, X. Ji, L. Jin and H. W. Lin,
  
  Phys. Rev. D 95, no. 9, 094514 (2017)
  doi:10.1103/PhysRevD.95.094514
  [arXiv:1702.00008 [hep-lat]].
  
  




Radyushkin:2016hsy
  A. Radyushkin,
  
  Phys. Lett. B 767, 314 (2017)
  doi:10.1016/j.physletb.2017.02.019
  [arXiv:1612.05170 [hep-ph]].
  
  




Radyushkin:2017cyf
  A. V. Radyushkin,
  
  Phys. Rev. D 96, no. 3, 034025 (2017)
  doi:10.1103/PhysRevD.96.034025
  [arXiv:1705.01488 [hep-ph]].
  
  




Radyushkin:2017ffo
  A. Radyushkin,
  
  Phys. Lett. B 770, 514 (2017)
  doi:10.1016/j.physletb.2017.05.024
  [arXiv:1702.01726 [hep-ph]].
  
  




Radyushkin:2017gjd
  A. V. Radyushkin,
  
  Phys. Rev. D 95, no. 5, 056020 (2017)
  doi:10.1103/PhysRevD.95.056020
  [arXiv:1701.02688 [hep-ph]].
  
  




Orginos:2017kos
  K. Orginos, A. Radyushkin, J. Karpie and S. Zafeiropoulos,
  
  Phys. Rev. D 96, no. 9, 094503 (2017)
  doi:10.1103/PhysRevD.96.094503
  [arXiv:1706.05373 [hep-ph]].
  
  




Karpie:2017bzm
  J. Karpie, K. Orginos, A. Radyushkin and S. Zafeiropoulos,
  
  arXiv:1710.08288 [hep-lat].
  
  




Radyushkin:2017lvu
  A. V. Radyushkin,
  
  arXiv:1710.08813 [hep-ph].
  
  




Ma:2014jga
  Y. Q. Ma and J. W. Qiu,
  
  Int. J. Mod. Phys. Conf. Ser.  37, 1560041 (2015)
  doi:10.1142/S2010194515600411
  [arXiv:1412.2688 [hep-ph]].
  
  




Ma:2017pxb
  Y. Q. Ma and J. W. Qiu,
  
  arXiv:1709.03018 [hep-ph].
  
  




Li:2016amo
  H. n. Li,
  
  Phys. Rev. D 94, no. 7, 074036 (2016)
  doi:10.1103/PhysRevD.94.074036
  [arXiv:1602.07575 [hep-ph]].
  
  




Li:2017udt
  H. n. Li,
  
  JPS Conf. Proc.  13, 020055 (2017).
  doi:10.7566/JPSCP.13.020055
  
  




Chen:2016fxx
  J. W. Chen, X. Ji and J. H. Zhang,
  
  Nucl. Phys. B 915, 1 (2017)
  doi:10.1016/j.nuclphysb.2016.12.004
  [arXiv:1609.08102 [hep-ph]].
  
  




Green:2017xeu
  J. Green, K. Jansen and F. Steffens,
  
  arXiv:1707.07152 [hep-lat].
  
  




Ji:2017oey
  X. Ji, J. H. Zhang and Y. Zhao,
  
  arXiv:1706.08962 [hep-ph].
  
  




Ishikawa:2017faj
  T. Ishikawa, Y. Q. Ma, J. W. Qiu and S. Yoshida,
  
  Phys. Rev. D 96, no. 9, 094019 (2017)
  doi:10.1103/PhysRevD.96.094019
  [arXiv:1707.03107 [hep-ph]].
  
  




Martinelli:1994ty
  G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas,
  
  Nucl. Phys. B 445, 81 (1995)
  doi:10.1016/0550-3213(95)00126-D
  [hep-lat/9411010].
  
  




Wang:2017qyg
  W. Wang, S. Zhao and R. Zhu,
  
  arXiv:1708.02458 [hep-ph].
  
  




Capitani:2002mp
  S. Capitani,
  
  Phys. Rept.  382, 113 (2003)
  doi:10.1016/S0370-1573(03)00211-4
  [hep-lat/0211036].
  
  




Polyakov:1980ca
  A. M. Polyakov,
  
  Nucl. Phys. B 164, 171 (1980).
  doi:10.1016/0550-3213(80)90507-6
  
  




Dotsenko:1979wb
  V. S. Dotsenko and S. N. Vergeles,
  
  Nucl. Phys. B 169, 527 (1980).
  doi:10.1016/0550-3213(80)90103-0
  
  




Brandt:1981kf
  R. A. Brandt, F. Neri and M. a. Sato,
  
  Phys. Rev. D 24, 879 (1981).
  doi:10.1103/PhysRevD.24.879
  
  




Gervais:1979fv
  J. L. Gervais and A. Neveu,
  
  Nucl. Phys. B 163, 189 (1980).
  doi:10.1016/0550-3213(80)90397-1
  
  




Samuel:1978iy
  S. Samuel,
  
  Nucl. Phys. B 149, 517 (1979).
  doi:10.1016/0550-3213(79)90005-1
  
  




Arefeva:1980zd
  I. Y. Arefeva,
  
  Phys. Lett.  93B, 347 (1980).
  doi:10.1016/0370-2693(80)90529-8
  
  




Dorn:1986dt
  H. Dorn,
  
  Fortsch. Phys.  34, 11 (1986).
  doi:10.1002/prop.19860340104
  
  




Dorn:1981wa
  H. Dorn, D. Robaschik and E. Wieczorek,
  
  Annalen Phys.  40, 166 (1983).
  
  




Craigie:1980qs
  N. S. Craigie and H. Dorn,
  
  Nucl. Phys. B 185, 204 (1981).
  doi:10.1016/0550-3213(81)90372-2
  
  



