Lattice Parton Collaboration ($\rm{\bf LPC}$)

# Impact of gauge fixing precision on the continuum limit
of non-local quark-bilinear lattice operators

Kuan Zhang  Corresponding author<EMAIL_ADDRESS>University of Chinese
Academy of Sciences, School of Physical Sciences, Beijing 100049, China  CAS
Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China  Yi-Kai Huo  Physics
Department, Columbia University, New York, NY 10027, USA  Xiangdong Ji
Department of Physics, University of Maryland, College Park, MD 20742, USA
Andreas Schäfer  Institut für Theoretische Physik, Universität Regensburg,
D-93040 Regensburg, Germany  Chun-Jiang Shi  Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  Peng Sun  Institute of
Modern Physics, Chinese Academy of Sciences, Lanzhou, 730000, China  Wei Wang
INPAC, Shanghai Key Laboratory for Particle Physics and Cosmology, Key
Laboratory for Particle Astrophysics and Cosmology (MOE), School of Physics
and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Southern
Center for Nuclear-Science Theory (SCNT), Institute of Modern Physics, Chinese
Academy of Sciences, Huizhou 516000, Guangdong Province, China  Yi-Bo Yang
Corresponding author<EMAIL_ADDRESS>University of Chinese Academy of
Sciences, School of Physical Sciences, Beijing 100049, China  CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  International Centre for Theoretical Physics Asia-Pacific,
Beijing/Hangzhou, China  Jian-Hui Zhang  School of Science and Engineering,
The Chinese University of Hong Kong, Shenzhen 518172, China

###### Abstract

We analyze the gauge fixing precision dependence of some non-local quark-
blinear lattice operators interesting in computing parton physics for several
measurements, using 5 lattice spacings ranging from 0.032 fm to 0.121 fm. Our
results show that gauge dependent non-local measurements are significantly
more sensitive to the precision of gauge fixing than anticipated. The impact
of imprecise gauge fixing is significant for fine lattices and long distances.
For instance, even with the typically defined precision of Landau gauge fixing
of $10^{-8}$, the deviation caused by imprecise gauge fixing can reach 12
percent, when calculating the trace of Wilson lines at 1.2 fm with a lattice
spacing of approximately 0.03 fm. Similar behavior has been observed in $\xi$
gauge and Coulomb gauge as well. For both quasi PDFs and quasi TMD-PDFs
operators renormalized using the RI/MOM scheme, convergence for different
lattice spacings at long distance is only observed when the precision of
Landau gauge fixing is sufficiently high. To describe these findings
quantitatively, we propose an empirical formula to estimate the required
precision.

##  I Introduction

In gauge field theories, the necessity of gauge fixing arises originally from
perturbation theory. If we just focus on non-perturbative lattice gauge theory
and gauge invariant quantities, we do not really need gauge fixing. Also, any
gauge dependent quantity will vanish if averaged over a sufficient number of
configurations. However, in certain non-perturbative renormalization schemes
Martinelli  _et al._ (1995), it is unavoidable to use gauge-dependent parton
matrix elements to renormalize composite operators in the Landau or other
gauges, and then match the results to the $\overline{\text{MS}}$ scheme
perturbatively. Furthermore, non-perturbative lattice gauge fixing is also
necessary to extract information from gauge dependent correlators Mandula
(1999).

In continuum gauge theory, the gauge potentials $A_{\mu}$ are the fundamental
degrees of freedom, and Wilson lines are defined through $A_{\mu}$ as

$\displaystyle U(x,y)\equiv{\cal P}{\rm
exp}\left[ig_{0}\int_{y}^{x}\textrm{d}z^{\mu}A^{m}_{\mu}(z)T^{m}\right].$ (1)

where $T^{m}$ are the $\mathrm{SU}(3)$ generators in the fundamental
representation, index $m$ is summed over and $g_{0}$ is the bare coupling
constant.

On the other hand, Wilson link $U_{\mu}(x)=U(x,x+a\hat{\mu})$ are treated as
the fundamental degrees of freedom in lattice gauge theory Wilson (1974),
where $a$ is the lattice spacing, and $\hat{\mu}$ is the unit vector in $\mu$
direction. The gauge potential $A_{\mu}$ can be extracted by short Wilson
lines up to discretization effects,

$\displaystyle A_{\mu}(x+\frac{a}{2}\hat{\mu})$
$\displaystyle\equiv\bigg{[}\frac{U_{\mu}(x)-U_{\mu}^{\dagger}(x)}{2ig_{0}a}\bigg{]}_{\text{Traceless}}.$
(2)

where $g_{0}$ is the bare coupling at lattice spacing $a$.

Based on the gauge transformation
$U^{G}_{\mu}(x)=G(x)U_{\mu}(x)G^{\dagger}(x+a\hat{\mu})$, the discretized
gauge condition for the Landau and Coulomb gauge can be defined as ($l=4$ for
Landau gauge and $l=3$ for Coulomb gauge),

$\displaystyle\Delta^{G}(x)\equiv$
$\displaystyle\sum^{l}_{\mu=1}\big{[}U^{G}_{\mu}(x)-U^{G}_{\mu}(x-a\hat{\mu})-(U^{G}_{\mu})^{\dagger}(x)$
$\displaystyle+(U^{G}_{\mu})^{\dagger}(x-a\hat{\mu})\big{]}/(2ig_{0})=0$ (3)

Above condition corresponds to continuum gauge fixing condition
$\sum_{\mu=1}^{l}\partial_{\mu}A^{G}_{\mu}=0$ with
$A^{G}_{\mu}=GA^{G}_{\mu}G^{\dagger}+\frac{1}{ig_{0}}G\partial_{\mu}G^{\dagger}$,
up to discretization errors.

The standard way to fix the Landau and Coulomb gauges in lattice QCD
Cabibbo:1982zn; Mandula and Ogilvie (1987); Katz  _et al._ (1988) is based on
the numerical minimization of the functional

$\displaystyle F_{U}[G]\equiv-\frac{1}{3lV}\text{Re}\text{
Tr}\sum_{x}\sum^{l}_{\mu=1}U_{\mu}^{G}(x),$ (4)

where $V$ is the lattice volume. The extremum of $F_{U}[G]$ corresponds to the
discretized gauge condition $\Delta^{G}(x)=0$ up to discretization errors
Giusti  _et al._ (2001a).

Numerical non-perturbative gauge fixing is never perfect. Therefore, it is
necessary to introduce a definition for the precision of gauge fixing. Since
the latter is unphysical, it allows for a wide range of choices. For an
appropriate definition of precision, as it approaches zero, the numerical
gauge fixing should approach perfect gauge fixing (here we ignore some subtle
issues like Gribov copies). And we expect that most of the gauge dependent
measurements we are interested in for our lattice simulation will behave
monotonically with this precision. Then we require a certain precision and
make sure that the impact of imprecise gauge fixing will not cause a
significant systematic deviation compared to the statistical uncertainty and
other systematic uncertainties. There are two popular criteria to estimate the
quality of gauge fixing. One is defined as,

$\displaystyle\theta^{G}\equiv\frac{1}{3V}\sum_{x}\text{Tr}[(\Delta^{G}(x))^{\dagger}(\Delta^{G}(x))].$
(5)

Obviously, when $\theta^{G}$=0, the discretized gauge condition
$\Delta^{G}(x)=0$ is met. The other is determined by the change in the
functional $F_{U}$ during the iterative gauge fixing procedure, given by the
equation:

$\displaystyle\delta^{F}(n)=F_{U}[G(n)]-F_{U}[G(n-1)],$ (6)

where $G(n)$ represents the gauge rotation at the $n$-th step. The gauge
fixing will stop at the $m$-th step once $\delta^{F}(m)$ is smaller than the
preassigned value $\delta^{F}$. When $\delta^{F}$=0, the functional $F_{U}[G]$
reaches its extremum, and the discretized gauge condition $\Delta^{G}(x)=0$ is
fulfilled. To the best of our knowledge, $\delta^{F}$ is the precision
typically defined in Landau gauge fixing codes. The magnitude of $\theta^{G}$
is proportional to that of $\delta^{F}$ as shown in the appendix, but larger
by a factor of $\sim$20\. In the following, we will just show the preassigned
value of $\delta^{F}$ as the precision of gauge fixing, except for the $\xi$
gauge. The relationship between them has also been explored in other works,
such as in Ref. Giusti  _et al._ (2001a).

For local quark bi-linear operators, $\delta^{F}\sim 10^{-7}$ is small enough
given the other uncertainty which are at the 0.1% level or larger Lytle  _et
al._ (2018). In this work, our interest is gauge-fixing precision study for
non-local quark bilinear operators. We find that non-local gauge dependent
measurements are much more sensitive to the precision of gauge fixing than
local ones. This has important implications for quasi-PDF and related matrix
elements renormalization in RI/MOM scheme or fixed-gauge quark bilinear non-
local operator calculations.

This paper is organized as follows. In Sec. II, we will present the numerical
setup, investigate the impact of imprecise gauge fixing on longer Wilson
lines, and propose an empirical formula to describe the impact of imprecise
gauge fixing. In the next two sections, we will show this impact on LaMET Ji
(2013); Ji  _et al._ (2021) lattice simulations. Quasi PDFs and Coulomb gauge
quasi PDFs Gao  _et al._ (2023) are shown in Sec. III while quasi TMD-PDFs are
shown in Sec IV. The last section will provide a brief summary.

##  II Numerical setup and Wilson lines

In this calcualtion, We use the 2+1+1 flavors (degenerate up and down,
strange, and charm degrees of freedom) of highly improved staggered quarks
(HISQ) and one-loop Symanzik improved Hart  _et al._ (2009) gauge ensembles
from the MILC Collaboration Bazavov  _et al._ (2013) at five lattice spacings.
We tuned the valence light quark mass using the HYP smeared tadpole improved
clover action to be around that in the sea, and the information about the
ensembles and parameters we used are listed in Table 1.

Tag | $6/g_{0}^{2}$ | $L$ | $T$ | $a(\mathrm{fm})$ | $m^{\textrm{w}}_{q}a$ | $c_{\mathrm{sw}}$
---|---|---|---|---|---|---
a12m310 | 3.60 | 24 | 64 | 0.1213(9) | -0.0695 | 1.0509
a09m310 | 3.78 | 32 | 96 | 0.0882(7) | -0.0514 | 1.0424
a06m310 | 4.03 | 48 | 144 | 0.0574(5) | -0.0398 | 1.0349
a045m310 | 4.20 | 64 | 192 | 0.0425(5) | -0.0365 | 1.0314
a03m310 | 4.37 | 96 | 288 | 0.0318(5) | -0.0333 | 1.0287
Table 1: Setup of the ensembles, including the bare coupling constant $g_{0}$,
lattice size $L^{3}\times T$ and lattice spacing $a$. $m^{\textrm{w}}_{q}$ is
the bare quark masses. The pion masses in the sea are $\sim$310 MeV and the
valence pion mass are in the range of 280-320 MeV.

The gauge fixing procedure used in this work is proposed in Ref.
Cabibbo:1982zn and can be briefly described as follows Edwards and Joo (2005):

1\. Separate the lattice sites into the even
($\mathrm{mod}(\sum_{i=1}^{4}x_{i},2)=0$) and odd
($\mathrm{mod}(\sum_{i=1}^{4}x_{i},2)=1$) sectors;

2\. On each site $x$ of the even sector, calculate
$K(x)\equiv\sum_{\mu=1}^{l}\big{(}U_{\mu}(x)+U_{\mu}^{\dagger}(x-a\hat{\mu})\big{)}$,
obtain the projected SU(2) matrix $a(x)$ from a given 2x2 sub-matrix of
$K(x)$, and then apply $a^{\dagger}(x)$ on the left hand side of $U$ and $K$.
One can apply the over-relaxation algorithm Adler:1981sn; Giusti  _et al._
(2001a) by replacing $a\equiv a_{0}+\sigma_{i}a_{i}$ into

$\displaystyle
a^{\prime}=\mathrm{cos}(\omega\theta)+\sigma_{i}a_{i}*\mathrm{sin}(\omega\theta)/\mathrm{sin}(\theta),$
(7)

where $\theta=\mathrm{cos}^{-1}a_{0}$, and $1<\omega<2$ to accelerate the
convergence.

3\. Repeat the step 2 on all the other 2x2 sub-matrices of $K(x)$ and apply
them to $U$ and $K$.

4\. Repeat the steps 2-3 for the odd sector.

5\. Compute $F_{U}$, and iterate through steps 2-4 until the change in
$F_{U}$, denoted by $\delta^{F}$, is less than the required gauge fixing
precision.

To get a clear image of imprecise gauge fixing, we choose to look at the ratio
of the gauge dependent quantity $X$ under two gauge fixing with precisions
$\delta^{F}$ and $\delta^{F}_{0}$,

$\displaystyle r^{\delta^{F}}(X)\equiv X^{\delta^{F}}/X^{\delta^{F}_{0}}.$ (8)

where $\delta^{F}_{0}$ is chosen to be “sufficiently small” to ensure that
$X^{\delta^{F}_{0}}$ remains unchanged within the statistical uncertainty. In
this definition, $r^{\delta^{F}}(X)$ can deviate from 1 when
$\delta^{F}\gg\delta^{F}_{0}$.

Obtaining higher precision in numerical gauge fixing necessitates more
iteration steps, resulting in increased computational cost, especially for
larger lattices such as a03m310. Consequently, in practice, it is necessary to
select a specific precision that strikes a balance between systematic error
and computational cost. In this study, we opt for the gauge rotation $G_{0}$
to ensure that the gauge fixing can attain the precision $\delta^{F}=10^{-12}$
in most cases. Subsequent discussions will demonstrate the reasonableness of
this choice.

###  II.1 Straight Wilson line under Landau gauge

We start with the most straightforward gauge dependent quantities in lattice
simulations, the Wilson lines. They are defined in Eq. (1), and can also be
expressed as,

$\displaystyle
U(x,x+na\hat{\mu})\equiv\prod_{k=0}^{n-1}U_{\mu}(x+ka\hat{\mu}),n>0.$ (9)

They are referred to as long Wilson lines when $n>1$. They are gauge dependent
and transform under a gauge transformation similar to the short line,

$\displaystyle
U^{G}(x,x+na\hat{\mu})=G(x)U(x,x+na\hat{\mu})G^{\dagger}(x+na\hat{\mu}).$ (10)

As shown in the upper and middle panels of Fig. 1, at short distance, the
ratio $r^{\delta^{F}}(W(z))$ of the trace of Wilson line with given length $z$

$\displaystyle W(z)\equiv\frac{1}{3}\text{Tr}\left[\sum_{x}U(x,x+z)\right],$
(11)

depends only weakly on precision and lattice spacing. But at long distance,
the ratio $r^{\delta^{F}}(W(z))$ decays more rapidly for smaller lattice
spacings and lower precisions. For example, on a03m310, at $z=$ 1.2 fm, the
difference between $\delta^{F}=10^{-8}$ and $\delta^{F}=10^{-12}$ is 11.6(9)
percent. This indicates that even with a high precision of
$\delta^{F}=10^{-8}$, the Wilson lines remain very sensitive to the precision
of gauge fixing, especially for small lattice spacings at long distance. As
for HYP smearing, it enhances the signal but does not result in additional
systematic deviations as shown in the lower panels of Fig. 1 using the a06m310
ensemble at $a\sim 0.06~{}$fm, in comparison to the case without HYP smearing
shown in the middle panel of the same figure.

Figure 1: The ratio $r^{\delta^{F}}(W(z))$ defined in Eq. (8) of the trace of
Wilson lines defined in Eq. (11) for different precisions, with one-step HYP
smearing after gauge fixing (upper and middle panels) and without HYP smearing
(lower panel) on different ensembles. The result at $\delta^{F}_{0}=10^{-12}$
is chosen as the denominator. The impact of imprecise gauge fixing is
significant at long distance on fine lattice.

To describe and estimate the impact of imprecise gauge fixing, we propose the
following empirical formula

$\displaystyle X(\delta^{F})=X(0)e^{-c(X)(\delta^{F})^{n(X)}},$ (12)

where $c$ and $n$ are the parameters to be fitted and $X(0)$ corresponds to
the observable $X$ using precise enough $G_{0}(x)$ (for example,
$\delta^{F}=10^{-15}$).

Figure 2:  The parameters $n(W(z))$ (top panel) and $c(W(z))$ (middle panel)
determined from the fitting ansatz Eq. (12) for the trace of Wilson lines
defined in Eq. (11) at different $z$, and also the corresponding
$\chi^{2}/$d.o.f. (bottom panel) for a12m310 and a06m310 ensembles. We use
data from gauge fixings with
$\delta^{F}=10^{-4},10^{-5},10^{-6},10^{-7},10^{-8}$. One-step HYP smearing
has been applied after gauge fixing. The uncertainty is got from bootstrap
resampling. The two curves in the middle panel are the fits of $c(W(z))$ using
a polynomial form.

As shown in the upper and middle panels of Fig. 2, $n(W(z))$ is approximately
0.5 and does not show a very strong dependence on $z$ and lattice spacing.
Conversely, $c(W(z))$ is much more sensitive to $z$ and lattice spacing. The
lower panel indicates that the $\chi^{2}/$d.o.f. is sufficiently small. Thus,
this empirical formula adequately captures the impact of imprecise gauge
fixing in the Wilson line case.

The impact of imprecise gauge fixing increases as the Wilson line becomes
longer, and it appears that this is primarily due to $c(W(z))$. Therefore, we
fit $c(W(z))$ and attempt to model it using a quadratic function
$c(W(z))=c_{2}z^{2}+c_{1}z+c_{0}$ as shown in the middle panel in Fig. 2. For
a12m310, we obtain $c_{2}=42(3)$, $c_{1}=7.2(1.3)$, $c_{0}=-0.38(0.12)$, and
$\chi^{2}/$d.o.f.=0.06. For a06m310, the values are $c_{2}=166(8)$,
$c_{1}=27.6(1.9)$, $c_{0}=-0.97(0.09)$, and $\chi^{2}/$d.o.f.=0.51. This
result suggests that $c(W(z))$ will increase quadratically as the Wilson line
lengthens, allowing us to use a quadratic function to forecast $c(W(z))$ and
assess the effect of imprecise gauge fixing. Interestingly, it appears that
the coefficient $c_{2}$ is proportional to $1/a^{2}$, and this term will
become dominant when the Wilson line is sufficiently long.

The Wilson lines with gauge rotations $G_{1}(x)$ and $G_{2}(x)$ can be related
by the following expression,

$\displaystyle U^{G_{1}}(x,x+z)$
$\displaystyle=G_{1}(x)G_{2}^{-1}(x)U^{G_{2}}(x,x+z)G_{2}(x+z)G_{1}^{-1}(x+z)$
$\displaystyle=[G_{1}(x)G^{-1}_{2}(x)]U^{G_{2}}(x,x+z)[G_{1}(x+z)G^{-1}_{2}(x+z)]^{-1},$
(13)

where $U^{G}(x,x+z)\equiv G(x)U(x,x+z)G^{\dagger}(x+z)$ and $U^{G_{1}}$ and
$U^{G_{2}}$ just differ by the relative gauge rotations at both ends of the
Wilson line. We can then use the correlation between the relative gauge
rotation $\tilde{G}\equiv GG_{0}^{-1}$ between the gauge fixing with different
precision,

$\displaystyle
C^{G}(z)=\frac{1}{3V}\sum_{x}\text{Tr}\left[[G(x+z)G_{0}^{-1}(x+z)]^{\dagger}[G(x)G_{0}^{-1}(x)]\right],$
(14)

to gain insights into why $W(z)$ depends on the gauge fixing precision. $G$
and $G_{0}$ are the gauge rotations with certain gauge fixing precision
$\delta^{F}$ and $\delta_{0}^{F}=10^{-12}\ll\delta^{F}$.

Figure 3: The correlation of the relative gauge rotation $\tilde{G}\equiv
GG^{-1}_{0}$ defined in Eq. (14) as function of the spatial distance $z$ for
different gauge fixing precision $\delta^{F}$ and different ensembles. The
gauge rotations $G$ and $G_{0}$ correspond to gauge fixing with precision
$\delta^{F}$ and $\delta^{F}_{0}=10^{-12}\ll\delta^{F}$.

As shown in the upper panel of Fig. 3, the correlation $C^{G}(z)$ approaches 1
when the gauge fixing is precise enough, with $\delta^{F}=\delta^{F}_{0}$. On
the other hand, $C^{G}(z)$ deviates from 1 for larger $\delta^{F}$ and
approaches 0 in the limit of no gauge fixing. Such a deviation becomes more
significant at large z, while suffering from loop-around effects when z is
approximately L/2. These observations are consistent with the idea that gauge
fixing creates a correlation between different spatial positions, which will
be stronger with more precise gauge fixing.

At the same time, $C^{G}(z)$ with a given $\delta^{F}$ but different lattice
spacing $a$ also increases with $1/a$, as shown in the lower panel of Fig. 3.
This is also understandable since the correlation would depend on the distance
in lattice units instead of physical units.

###  II.2 Straight Wilson line under $\xi$ gauge

We also extend our study of Wilson lines from Landau gauge with $\xi=0$ to the
general covariant $\xi$ gauge. The gauge condition of $\xi$ gauge is defined
as Giusti  _et al._ (2001b); Zwanziger (1986); Giusti (1997),

$\displaystyle\Delta^{G}(x)\equiv\sum^{4}_{\mu}(A^{G}_{\mu}(x+\frac{a}{2}\hat{\mu})-A^{G}_{\mu}(x-\frac{a}{2}\hat{\mu}))=\Lambda(x),$
(15)

where random matrix $\Lambda(x)$ on each lattice site $x$ of each
configuration belonging to the $\mathrm{SU}(3)_{c}$ group are sampled
according to the weight $e^{-\frac{1}{2\xi}\text{Tr}\Lambda^{2}(x)}$. Note
that the Landau gauge is recovered when $\Lambda(x)=0$ at every $x$. The
precision of gauge fixing for the $\xi$ gauge is defined as,

$\displaystyle\theta^{G}_{\xi}\equiv\frac{1}{3V}\sum_{x}\text{Tr}[(\Delta^{G}(x)-\Lambda(x))^{\dagger}(\Delta^{G}(x)-\Lambda(x))].$
(16)

We use this expression in $\xi$ gauge instead of $\delta^{F}$ in this case.

Figure 4: The ratio $r^{\theta^{G}_{\xi}}(W(z))$ of the trace of Wilson lines
defined in Eq. (11) for different precision with one-step HYP smearing. The
result with the precision $\theta_{\xi,0}^{G}=2\times 10^{-7}$ is used as the
denominator in the ratio. The impact of imprecise gauge fixing is similar to
that in Landau gauge.

As we did for Landau gauge, we plot the ratio

$\displaystyle r^{\theta^{G}_{\xi}}(W)\equiv
W^{\theta^{G}_{\xi}}/W^{\theta^{G}_{\xi,0}}$ (17)

of the traces of Wilson lines defined in Eq. (11) for different precision of
gauge fixing, with $\theta^{G}_{\xi,0}=2\times 10^{-7}$ since it is very hard
to reach higher precision here. The numerical gauge fixing requires more
iteration steps to converge when $\xi$ becomes larger. Additionally, some
configurations can not converge within 100,000 iterations if the required
precision is too high, particularly for large $\xi$.

As shown in the upper panel of Fig. 4, it seems that the impact of imprecise
gauge fixing becomes slightly smaller as $\xi$ increases. The middle and lower
panels show that the effects of lower required precision are larger on finer
lattices, at longer distances and for lower average precision. This
observation is consistent with our findings in Landau gauge.

###  II.3 Staple shaped Wilson line under Landau gauge

Staple shaped Wilson lines ${\cal W}(b,z,L)$ are also widely used gauge
dependent measurement. For example, a quasi TMD-PDF operator Ji  _et al._
(2020) contains a staple shaped Wilson line. Staple shaped Wilson lines are
defined as,

$\displaystyle{\cal W}(b,z,L)$ $\displaystyle\equiv{\cal P}{\rm
exp}\left[ig_{0}\int_{-L}^{z}\textrm{d}s\ \hat{n}_{z}\cdot
A^{m}(b\hat{n}_{\perp}+s\hat{n}_{z})T^{m}\right]$ $\displaystyle\times{\cal
P}{\rm exp}\left[ig_{0}\int_{0}^{b}\textrm{d}s\ \hat{n}_{\perp}\cdot
A^{m}(s\hat{n}_{\perp}-L\hat{n}_{z})T^{m}\right]$ $\displaystyle\times{\cal
P}{\rm exp}\left[ig_{0}\int_{0}^{-L}\textrm{d}s\ \hat{n}_{z}\cdot
A^{m}(s\hat{n}_{z})T^{m}\right],$ (18)

where $b=|\vec{b}_{\perp}|$ and $z$ are the separations between the two
endpoints of the staple shaped Wilson line along the transverse direction
$\hat{n}_{\perp}$ and longitudinal direction $\hat{n}_{z}$, respectively.
Straight Wilson lines are recovered when $z=0,L=0$. The trace is denoted as
$W_{\rm st}(b,z,L)$ as in the Wilson line case.

Figure 5: The ratio $r^{\delta^{F}}(W_{\rm st}(b,z,L))$ Eq.(8) of the trace of
staple shaped Wilson lines defined in Eq. (II.3) for different precision with
one-step HYP smearing after gauge fixing. We use the $z=0$ case as an
illustration. The result for $\delta^{F}_{0}=10^{-12}$ is chosen as the
denominator. It is observed that the impact of imprecise gauge fixing depends
on $L$, but this dependence is much weaker than that on $b$.

When considering the staple shaped Wilson lines ${\cal W}(b,z,L)$ and its
trace $W_{\rm st}(b,z,L)$ for different gauge fixing, the only difference lies
in the gauge transformation $G(x)$ on the two endpoints. When $b$ and $z$
remain unchanged, the endpoints are fixed, but $L$ still has an influence on
the impact of imprecise gauge fixing, as depicted in Fig. 5. When $L$ becomes
longer, the impact of imprecise gauge fixing is enhanced as shown in the upper
panel of Fig. 5. Furthermore, longer $b$ and smaller lattice spacing make this
phenomena more pronounced, as shown in the middle and lower panels of Fig. 5.
As the gauge fixing precision dependence is not very sensitive on $\xi$, we
will not delve into the discussion of the staple-shaped Wilson line in the
$\xi$ gauge.

##  III Quasi PDF

A quasi-PDF operator is a quark bilinear operator with a Wilson line,

$\displaystyle O^{\text{GI}}_{\Gamma}(z)=\bar{\psi}(0)\Gamma U(0,z)\psi(z).$
(19)

Its renormalized nucleon matrix element is connected to PDFs through
factorization theorems Ji (2013, 2014); Ma and Qiu (2018a, b); Izubuchi  _et
al._ (2018). Its hadron matrix element is

$\displaystyle\tilde{h}^{\text{GI}}_{\chi,\Gamma}(z,P_{z};1/a)\equiv\langle\chi(P_{z})|O^{\text{GI}}_{\Gamma}(z)|\chi(P_{z})\rangle.$
(20)

The UV divergence is independent of the hadron state, as shown in Refs. Zhang
_et al._ (2021); Huo  _et al._ (2021). In the following we will concentrate on
the pion matrix element
$\tilde{h}^{\text{GI}}_{\pi}(z,P_{z};1/a)\equiv\langle\pi(P_{z})|O^{\text{GI}}_{\gamma_{t}}(z)|\pi(P_{z})\rangle$
of the quasi PDF operator $O^{\text{GI}}_{\gamma_{t}}(z)$ at different lattice
spacing $a$ as an illustrative example. The pion matrix element can be
extracted from the following ratio

$\displaystyle R_{\pi}(t_{1},z;a,t_{2})$ $\displaystyle\equiv\frac{\langle
O_{\pi}(t_{2})\sum_{\vec{x}}O^{\text{GI}}_{\Gamma}(z;(\vec{x},t_{1}))O_{\pi}^{\dagger}(0)\rangle}{\langle
O_{\pi}(t_{2})O_{\pi}^{\dagger}(0)\rangle}$
$\displaystyle=\langle\pi|O^{\text{GI}}_{\Gamma}(z)|\pi\rangle+{\cal
O}(e^{-\Delta mt_{1}},e^{-\Delta m(t_{2}-t_{1})},e^{-\Delta mt_{2}}),$ (21)

where $\Delta m$ is the mass gap between the pion and its first excited state
which is aournd 1 GeV. One-step HYP smearing is applied to the Wilson lines to
enhance the signal to noise ratio (SNR). The source/sink setup and also the
ground state matrix element extraction are similar to that in Refs. Zhang  _et
al._ (2021); Huo  _et al._ (2021); Zhang  _et al._ (2022).

###  III.1 Quasi PDF renormalized using the Wilson line

In lattice regularization, one encounters a linear divergence of the quasi PDF
operator $O^{\text{GI}}_{\gamma_{t}}(z)$ defined in Eq. (19). To ensure the
existence of a finite continuum limit, it is necessary to remove this linear
divergence. Studies in the continuum Ji and Zhang (2015); Ji  _et al._ (2018);
Ishikawa  _et al._ (2017); Green  _et al._ (2018) suggest that the quasi-PDF
operator is multiplicatively renormalizable, and the linear divergence just
comes from the Wilson line. We plot the ratio of the pion matrix element
$\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,0;1/a)$ and the trace of the Wilson
line $W(z)$ to see whether there is a residual linear divergence on the
lattice.

The central value of the pion matrix element with Coulomb wall source shows
negligible dependence on the precision of Coulomb gauge fixing, as
demonstrated in the appendix. So we will concentrate on the precision of
Landau gauge fixing used for the Wilson line rather than Coulomb gauge fixing
used for the pion matrix element with wall source.

Figure 6: Pion matrix element in lattice regularization
$\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,0;1/a)$ defined in Eq. (20) divided
by $W(z)$ defined in Eq. (11) for different lattice spacings and different
precision of gauge fixing. One-step HYP smearing has been applied after gauge
fixing.

As shown in Fig. 6, the results at $\delta^{F}=10^{-7}$ differ for different
lattice spacings. The discrepancy increases for finer lattices and longer
Wilson lines. When the precision of gauge fixing is increased, such as to
$\delta^{F}=10^{-8}$, the results show convergence. Moreover, for
$\delta^{F}=10^{-12}$, the results from different lattice spacings are
consistent within uncertainty even for the finest lattice and at long
distance, indicating the complete elimination of the linear divergence.
Although at short distance, the results for $\delta^{F}=10^{-7}$,
$\delta^{F}=10^{-8}$ and $\delta^{F}=10^{-12}$ are very close to each other,
the impact of imprecise gauge fixing becomes significant at long distance,
especially for smaller lattice spacing. This is due to the unphysical
contribution caused by imprecise gauge fixing of Wilson lines, as mentioned
previously.

Since the ratio has removed the linear divergence properly, one can do the
perturbative calculation of $W(z)$ to match the hadron matrix element
$\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}$ to the $\overline{\textrm{MS}}$
scheme. First define a ”subtracted” quasi PDF renormalized in the short-
distance ratio (SDR) scheme which divide the quasi PDF by the Wilson line and
normalized with a pion matrix element in the rest frame at a short distance
$z_{0}$,

$\displaystyle h^{{\rm
GI,SDR}}_{\chi,\gamma_{t}}(z,P_{z};1/z_{0})=\frac{\tilde{h}^{\text{GI}}_{\chi,\gamma_{t}}(z,P_{z};1/a)/W(z)}{\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z_{0},0;1/a)/W(z_{0})}.$
(22)

Then $h^{{\rm GI,SDR}}_{\chi,\gamma_{t}}$ can be matched to the
$\overline{\textrm{MS}}$ scheme using the perturbative calculation of $W(z)$,
$W(z_{0})$ and the on-shell $\tilde{h}^{\text{GI}}_{\chi,\gamma_{t}}$ which is
independent of $\chi$,

$\displaystyle h^{{\rm
GI,\overline{\textrm{MS}}}}_{\chi,\gamma_{t}}(z,P_{z};\mu)=$ $\displaystyle
C(z,z_{0},\mu)h^{{\rm GI,SDR}}_{\chi,\gamma_{t}}(z,P_{z};z_{0}),$ (23)

where

$\displaystyle
C(z,z_{0},\mu)=\frac{\tilde{h}^{\text{GI},\overline{\textrm{MS}}}_{{\rm
pert},\gamma_{t}}(z_{0},0;\mu)}{W^{\overline{\textrm{MS}}}_{\rm
pert}(z_{0},\mu)}W^{\overline{\textrm{MS}}}_{\rm pert}(z,\mu)$
$\displaystyle=\Big{\\{}1+\frac{\alpha_{s}C_{F}}{4\pi}\Big{[}\xi\text{ln}z_{0}^{2}\mu^{2}+5+\xi(2\gamma_{E}-\text{ln}4)\Big{]}\Big{\\}}$
$\displaystyle\quad\times\Big{\\{}1+\frac{\alpha_{s}C_{F}}{4\pi}(3-\xi)\Big{[}\text{ln}z^{2}\mu^{2}+2\gamma_{E}-\text{ln}4\Big{]}\Big{\\}}+{\cal
O}(\alpha_{s}^{2})$ $\ {}_{\overrightarrow{\xi\rightarrow 0}}\
\tilde{h}^{\text{GI},\overline{\textrm{MS}}}_{{\rm
pert},\gamma_{t}}(z,0;\mu)+{\cal O}(\alpha_{s}^{2}).$ (24)

Note that since the long Wilson line does contain non-perturbative physics at
large-$z$, the above approach will introduce some unknown non-perturbative
effects. The comparison with other renormalization schemes can be found in the
following subsection.

###  III.2 Quasi PDF renormalized in RI/MOM scheme and self renormalization

In RI/MOM renormalization, we need to calculate gauge dependent quark matrix
elements to renormalize hadron matrix elements to get results free of any
residual divergences. The definition of the renormalization constant is (e.g.,
Liu  _et al._ (2020)),

$\displaystyle Z^{\text{GI}}_{\Gamma}(z,P_{z};\mu,1/a)$
$\displaystyle=\frac{Z_{q}(P_{z};\mu,1/a)}{\textrm{Tr}[\Gamma\langle
q(p)|O^{\text{GI}}_{\Gamma}(z)|q(p)\rangle]_{p^{2}=-\mu^{2},p_{z}=p_{t}=0}},$
(25)

where $|q(p)\rangle$ is the off-shell quark state with external momentum $p$,
and we have assumed a lattice regularization and included the $1/a$ dependence
in the matrix elements. The setup in the calculation of quark matrix elements
is the same as that in Ref. Zhang  _et al._ (2021), except for the precision
of gauge fixing. $Z_{q}$ is defined from the hadron $\chi$ matrix element of a
local vector current,

$\displaystyle Z_{q}(P_{z};\mu,1/a)$
$\displaystyle=\frac{\textrm{Tr}[\Gamma\langle
q(p)|O^{\text{GI}}_{\Gamma}(0)|q(p)\rangle]_{p^{2}=-\mu^{2},p_{z}=p_{t}=0}}{\tilde{h}^{\text{GI}}_{\chi,\Gamma}(0,P_{z};1/a)}.$
(26)

When the Landau gauge fixed volume source Chen  _et al._ (2018) is used in the
calculation, the statistical uncertainty is suppressed by a factor
$1/\sqrt{V}$ compared to the point source case. Note that the $Z_{q}$
definition used here is exactly the same as
$\tilde{Z}_{q}\equiv\textrm{Tr}[p\\!\\!\\!/S^{-1}]/p^{2}$ from the quark
propagator in dimensional regularization Gracey (2003), but the discretization
effect is much smaller than that for $\tilde{Z}_{q}$ in lattice regularization
as shown in Ref. Chang  _et al._ (2021).

We multiply the pion matrix element
$\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,Pz;1/a)$ with $Z^{\rm
GI}_{\gamma_{t}}$ to obtain the non-perturbatively renormalized and normalized
matrix element at a given RI/MOM scale $\mu$,

$\displaystyle
h^{{\text{GI}},\text{RI/MOM}}_{\pi,\gamma_{t}}(z,P_{z};\mu)=Z^{\text{GI}}_{\gamma_{t}}(z,P_{z};\mu,1/a)\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,P_{z};1/a).$
(27)

Since the quasi-PDF operator is multiplicatively renormalizable and the linear
divergence just comes from the Wilson line and is independent of the external
quark or hadron state, RI/MOM renormalization Martinelli  _et al._ (1995) is
believed to be a good candidate to remove the linear divergence Constantinou
and Panagopoulos (2017); Green  _et al._ (2018); Chen  _et al._ (2018);
Stewart and Zhao (2018); Alexandrou  _et al._ (2017, 2020); Lin  _et al._
(2020).

Figure 7: The ratio
$r^{\delta^{F}}(Z^{\text{GI}}_{\gamma_{t}}(z,P_{z};\mu,1/a))$ of the
renormalization factor defined in Eq. (III.2) in RI/MOM scheme for two values
of the required precision, with one-step HYP smearing. The case
$\delta^{F}_{0}=10^{-12}$ is used as denominator in the ratios.

If the linear divergence just comes from the Wilson line and is independent of
the external quark or hadron state, then RI/MOM renormalization Martinelli
_et al._ (1995) can be a good candidate to remove the linear divergence
Constantinou and Panagopoulos (2017); Green  _et al._ (2018); Chen  _et al._
(2018); Stewart and Zhao (2018); Alexandrou  _et al._ (2017, 2020); Lin  _et
al._ (2020). But similar to Wilson line, the quasi PDF quark matrix element in
RI/MOM scheme is a non-local gauge dependent quantity, and thus it is also
highly sensitive to the precision of Landau gauge fixing at long distance on
fine lattices as shown in Fig. 7.

Compared to Fig. 1, we observe that the impact of imprecise gauge fixing on
the quark matrix element is weaker than that on the Wilson line. This suggests
that the quark external state mitigates the impact of imprecise gauge fixing.
For example, the difference of the quark matrix elements at
$\delta^{F}=10^{-8}$ and $\delta^{F}=10^{-12}$ on a06m310 at 1.2 fm is
0.021(0.012) whereas this difference is 0.032(0.003) for the Wilson line.
Similarly, on a03m310, the difference is 0.06(0.03) and 0.116(0.09) for the
quark matrix element and the Wilson line, respectively.

Figure 8: The fitted results for $n(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$
and $c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$ defined in Eq. (III.2) and
Eq. (12), along with the corresponding $\chi^{2}/$d.o.f. on different
ensembles. We use data for gauge fixing with
$\delta^{F}=10^{-4},10^{-5},10^{-6},10^{-7},10^{-8}$. one-step HYP smearing
has been applied after gauge fixing. The uncertainty is got by bootstrap
resampling. The two curves in the middle panel are the fitted results for
$c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$.

It is important to check whether the empirical formula defined in Eq. (12) can
also describe the impact of imprecise gauge fixing on quark matrix elements.
The fitted results are presented in Fig. 8. Here,
$n(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$ is about 0.5 and
$c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$ increases for longer $z$ and
smaller lattice spacing. The lower panel indicates that $\chi^{2}/$d.o.f. is
not consistently small and can reach around 3 on a06m310. Nevertheless, we
maintain that this empirical formula performs reasonably well. Ultimately, the
precision of gauge fixing is unphysical, and there is some flexibility in its
definition. Since $\delta^{F}$ is simply the trace of the Wilson line, it is
more suitable for the Wilson line case than for the quark matrix element case
when investigating the impact of imprecise gauge fixing. We believe that this
empirical formula is sufficient to help us estimate the impact of imprecise
gauge fixing for the quark matrix element. Additionally, we find that
$c(Z_{q})$ is negative on both a12m310 and a06m310, as
$n(Z_{q}(P_{z};\mu,1/a))=0.48(0.02)$, $c(Z_{q}(P_{z};\mu,1/a))=-1.1(0.2)$ and
$\chi^{2}/$d.o.f.=0.30 on a12m310, $n(Z_{q}(P_{z};\mu,1/a))=0.507(0.09)$,
$c(Z_{q}(P_{z};\mu,1/a))=-6.0(0.5)$ and $\chi^{2}/$d.o.f.=0.36 on a06m310. The
negative value of $c(Z_{q})\propto 1/a^{2}$ implies that $Z_{q}$ defined in
Eq. (III.2) will increase as the precision of Landau gauge fixing decreases,
and that its impact can increase significantly at small lattice spacing.

We can perform a similar fit for $c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$
using a quadratic function
$c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))=c_{2}z^{2}+c_{1}z+c_{0}$ as shown
in the middle panel in Fig. 8. On a12m310, we obtain $c_{2}=39(4)$,
$c_{1}=3.5(1.4)$, $c_{0}=-0.20(0.12)$, with $\chi^{2}/$d.o.f.=0.15. For
a06m310, the values are $c_{2}=136(10)$, $c_{1}=19(2)$, $c_{0}=-0.80(0.10)$,
and $\chi^{2}/$d.o.f.=1.40. This result indicates that
$c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$ will grow quadratically as the
Wilson line becomes longer. Consequently, we can employ a quadratic function
to predict $c(Z^{\text{GI}}_{\gamma_{t}}(z,0;\mu,1/a))$ and estimate the
impact of gauge fixing.

Figure 9: Quasi PDF matrix element in RI/MOM scheme defined in Eq. (11) at
different lattice spacings and for different precision of Landau gauge fixing.

Thus the imprecise gauge fixing has an significant impact on quasi PDF in
RI/MOM scheme. As shown in Fig. 9, the outcome at $\delta^{F}=10^{-7}$
exhibits inconsistencies across various lattice spacings, with a more
pronounced discrepancy observed for finer lattices. Furthermore, the
discrepancy amplifies as the line lengthens. However, with an increased
precision of gauge fixing, such as $\delta^{F}=10^{-8}$, the results behave
better. Additionally, for $\delta^{F}=10^{-12}$, the results for different
lattice spacings agree within the uncertainties, even on the finest lattice
and at long distances, indicating the complete elimination of the linear
divergence. If the precision is not high enough, the divergent behavior may be
mistaken for residual linear divergence, as what we did in Ref. Zhang  _et
al._ (2021).

Even though the RI/MOM renormalization with precise enough gauge fixing can
eliminate the linear divergence, the non-perturbative effect at large $z$
remains. Thus one can also extract divergent renormalization factors of quasi-
PDF by parameterizing the results in the rest frame on the lattice, which is
called self renormalization Huo  _et al._ (2021).

The fit strategy used in Ref. Huo  _et al._ (2021) can be improved into a
joint fit of the results at different lattice spacing $a$ and Wilson line
length $z$:

$\displaystyle\text{ln}\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,0;1/a)=\frac{kz}{a\text{ln}[a\Lambda_{\text{QCD}}]}+f(z)a+m_{0}z$
$\displaystyle\quad+\frac{3C_{f}}{b_{0}}\text{ln}\bigg{[}\frac{\text{ln}(1/a\Lambda_{\text{QCD}})}{\text{ln}(\mu/\Lambda_{\text{QCD}})}\bigg{]}+\frac{1}{2}\text{ln}\bigg{[}1+\frac{d}{\text{ln}[a\Lambda_{\text{QCD}}]}\bigg{]}^{2}$
$\displaystyle\quad+\begin{cases}\text{ln}\big{[}Z_{\overline{{\text{MS}}}}(z,\mu,\Lambda_{\overline{{\text{MS}}}})\big{]}&\text{if
}z_{0}\leq z\leq z_{1}\\\ g(z)&\text{if }z_{1}<z\end{cases},$ (28)

where the linear divergence factor $k$, the effective QCD scale
$\Lambda_{\text{QCD}}$, the discretization effect coefficient $f(z)$, the
renormalon factor $m_{0}$, the overall $\alpha_{s}$ correction coefficient
$d$, and the renormalized matrix element $g(z)$ are the parameters to be
fitted, and $Z_{\overline{{\text{MS}}}}(z,\mu,\Lambda_{\overline{\text{MS}}})$
represents the quasi-PDF in the $\overline{\text{MS}}$ scheme at 1-loop level.

It is important for $z_{0}$ and $z_{1}$ to be sufficiently large to mitigate
discretization effects, yet small enough to ensure the validity of
perturbative theory. Additionally, $z_{1}-z_{0}$ should not be excessively
small. In this work, we choose $z_{0}=0.06~{}\text{fm}$ and
$z_{1}=0.15~{}\text{fm}$, and interpolate the data at different lattice
spacing with intervals of 0.03 fm since both $f(z)$ and $g(z)$ depend on $z$.

Based on this strategy, the renormalized quasi-PDF in the rest frame in the
$\overline{\text{MS}}$ scheme will be

$\displaystyle\tilde{h}^{\text{GI}}_{\pi,\gamma_{t}}(z,0;1/a)=\begin{cases}Z_{\overline{{\text{MS}}}}(z,\mu,\Lambda_{\overline{{\text{MS}}}})&\text{if
}z_{0}\leq z\leq z_{1}\\\ \mathrm{Exp}[g(z)]&\text{if }z_{1}<z\end{cases}.$
(29)

Figure 10: Quasi PDF matrix element in $\overline{\text{MS}}$ scheme with
scale $\mu=2\text{GeV}$. The results through the intermediate RI/MOM and SDR
schemes are calculated on a06m310 and converted to the $\overline{\text{MS}}$
scheme by perturbative matching at 1-loop level. The result from self
renormalization is extracted from joint fits with 5 lattice spacings, and $z$
is interpolated to multiples of 0.03 fm. The curve is the perturbative result
at 1-loop level in the $\overline{\text{MS}}$ scheme.

The conversion factor between RI/MOM and $\overline{\text{MS}}$ at 1-loop
level can be found in Ref. Constantinou and Panagopoulos (2017). As shown in
Fig. 10, the $\overline{\text{MS}}$ results using the self renormalization,
the intermediate RI/MOM or SDR schemes are consistent with each other and
close to the perturbative result at short distance, but they deviate from each
other at long distance due to the non-perturbative effect from the intrinsic
scale $1/z$ introduced in RI/MOM and SDR schemes.

###  III.3 Coulomb gauge quasi PDF

The Coulomb gauge quasi PDF Gao  _et al._ (2023) (CG qPDF) operator is

$\displaystyle O^{\text{CG}}_{\Gamma}(z)=\bar{\psi}(0)\Gamma\psi(z).$ (30)

At next-to leading order (NLO), its LaMET matching to PDFs has been
demonstrated in Ref. Gao  _et al._ (2023). This operator is non-local and its
hadron matrix element $\tilde{h}^{\text{CG}}_{\chi,\Gamma}(z,P_{z};1/a)$ is
gauge dependent because the quark fields in the operator are non-local and not
linked by Wilson lines. If it is the free of the linear divergence, the
residual log divergence can be removed by a short distance matrix element,

$\displaystyle
h^{\text{CG}}_{\chi,\Gamma}(z,P_{z};1/a)\equiv\tilde{h}^{\text{CG}}_{\chi,\Gamma}(z,P_{z};1/a)/\tilde{h}^{\text{CG}}_{\chi,\Gamma}(z_{0},P_{z};1/a).$
(31)

The hadron matrix element is non-local and gauge dependent, similar to the
quark matrix elements of quasi PDFs and quasi TMD-PDFs. Therefore, it is
expected to be sensitive to the precision of gauge fixing.

Figure 11: Renormalized Coulomb gauge quasi PDF
$h^{\text{CG}}_{\pi,\gamma_{t}}(z,0;1/a)$ for different required precision of
gauge fixing and on different ensembles. We studied the unpolarized case
($\Gamma=\gamma_{t}$) for the pion in the rest frame as an illustration.
$z_{0}$ is interpolated to 0.3 fm in this work.

We performed calculations on the MILC ensembles a06m310 and a045m310, with
lattice spacings close to 0.06 fm and 0.04 fm, respectively, as used in Ref.
Gao  _et al._ (2023). As shown in Fig. 11, the results are not consistent with
each other when the precision of Coulomb gauge fixing is set to
$\delta^{F}=10^{-6}$. If we improve the precision to $\delta^{F}=10^{-7}$, the
results are more consistent. This suggests absence of linear divergences or
residual logarithm divergences in normalized CG qPDFs.

##  IV Quasi TMD-PDF

As a natural generalization of the collinear PDFs, the transverse-momentum-
dependent (TMD) PDFs provide a useful description of the transverse structure
of hadrons. In LaMET, the calculation of TMD-PDFs starts from the unsubtracted
quasi TMD-PDF Ji  _et al._ (2020); Ebert  _et al._ (2019); Ji  _et al._
(2021). The quasi TMD-PDF operator is a quark bilinear operator with staple-
shaped Wilson line,

$\displaystyle O^{\text{TMD}}_{\Gamma}(b,z,L)$
$\displaystyle\equiv\bar{\psi}(\vec{0}_{\perp},0)\Gamma{\cal
W}(b,z,L)\psi(\vec{b}_{\perp},z).$ (32)

$L$ should be large enough to approximate the infinitely long Wilson line in
the continuum. Similar as for the quasi PDF case, we focus on the pion matrix
element and the unpolarized case ($\Gamma=\gamma_{t}$). Details in the RI/MOM
scheme can be found in Ref. Zhang  _et al._ (2022). Similar to quasi PDFs,
quasi TMD-PDF hadron matrix elements are denoted as
$\tilde{h}^{\text{TMD}}_{\chi,\Gamma}(b,z,P_{z};1/a)$, and renormalized matrix
elements in the RI/MOM scheme are denoted as
$h^{{\text{TMD}},\text{RI/MOM}}_{\chi,\Gamma}(b,z,P_{z};\mu)$.

Figure 12: Quasi TMD-PDF matrix element in RI/MOM scheme. The operator
$O^{\text{TMD}}_{\Gamma}(b,z,L)$ is defined in Eq. (32). The upper panel shows
the $L$ dependence for different $b$ on a06m310. The lower panel shows the $b$
dependence on five ensembles. All the quark matrix elements are calculated at
$\delta^{F}=10^{-8}$.

As shown in the upper panel of Fig. 12, the results from different $b$ all
saturate to a plateau. This indicates that the pinch pole singularity has been
eliminated by the quark matrix element on the lattice. For a well-defined
quasi-TMD operator, $L$ should be larger than the spatial parton distribution
and independent of $L$ after proper subtraction. However, the statistical
uncertainty will be larger at longer $L$. In practice, we can choose a proper
$L$ to balance the statistical uncertainty and systematic uncertainty, as we
did in Ref. Zhang  _et al._ (2022).

As shown in the lower panel of Fig. 12, the results at $\delta^{F}=10^{-8}$
from different lattice spacings exhibit consistency, similar to the behavior
observed in the quasi PDF case with $\delta^{F}=10^{-12}$. This suggests that
the elimination of linear divergence, logarithmic divergence and cusp
divergence have been achieved, and the precision $\delta^{F}=10^{-8}$ is
adequate in this case. If the precision of Landau gauge fixing is
$\delta^{F}=10^{-6}$ as as in Ref. Zhang  _et al._ (2022), convergent behavior
cannot be observed in the RI/MOM scheme.

The conversion factor between RI/MOM and $\overline{\text{MS}}$ can be found
in Ref. Constantinou  _et al._ (2019); Ebert  _et al._ (2020). The reason why
their perturbative results do not converge is that there is a pinch pole
singularity in quasi TMD-PDFs. As $L$ increases, the conversion factor does
not converge, resulting in a huge value. One method that can be used to remove
the pinch pole singularity Ji  _et al._ (2021) is dividing quasi TMD-PDF
matrix elements by the square root of a rectangular Wilson loop. The
calculation on a lattice will not be affected because the contributions of a
Wilson loop in the quark matrix element and the hadron matrix element cancel
each other out. This procedure has also been mentioned in other work, such as
Ref. Alexandrou  _et al._ (2023); Spanoudes  _et al._ (2024). What we really
need to do is adding a Wilson loop term to the conversion factor.

Figure 13: Quasi TMD-PDF matrix element on a06m310 in $\overline{\text{MS}}$
scheme, $\mu=\text{2GeV}$. The result from the RI/MOM scheme is converted to
the $\overline{\text{MS}}$ scheme by perturbative matching at 1-loop level.
The curve is the perturbative result at 1-loop level in the
$\overline{\text{MS}}$ scheme.

On the other hand, Ref. Zhang  _et al._ (2022) proposed the short distance
ratio scheme which is gauge invariant and can also be converted to the
$\overline{\text{MS}}$ scheme through perturbative matching. To get rid of the
pinch pole singularity and the linear divergence from the Wilson link self-
energy, a more convenient “subtracted” quasi-TMD-PDF is then formed Ji  _et
al._ (2020),

$\displaystyle h^{\text{TMD}}_{\chi,\gamma_{t}}(b,z,P_{z};1/a)$
$\displaystyle=\lim_{L\to\infty}\frac{\tilde{h}^{\text{TMD}}_{\chi,\gamma_{t}}(b,z,L,P_{z};1/a)}{\sqrt{Z_{E}(b,2L+z;1/a)}},$
(33)

We write the fully renormalized quasi-TMD-PDF as

$\displaystyle h^{\text{TMD},{\rm
SDR}}_{\chi,\gamma_{t}}(b,z,P_{z};\frac{1}{b_{0}})=\frac{h^{\text{TMD}}_{\chi,\gamma_{t}}(b,z,P_{z};1/a)}{h^{\text{TMD}}_{\pi,\gamma_{t}}(b_{0},z_{0}=0,0,1/a)},$
(34)

where we also chosen the pion matrix element in the denominator (denoted by
the subscript $\pi$) likes the quasi-PDF case. For simplicity, we have also
chosen $z_{0}=0$. The singular dependence on $a$ on the r.h.s. of Eq. (34) has
been cancelled, leaving a dependence on the perturbative short scale $b_{0}$.
To perform the renormalization, the pion matrix element
$h_{\pi,\gamma_{t}}(b_{0},0,0,1/a)$ needs to be calculated non-perturbatively
on the lattice. In order to match the renormalized quasi-TMD-PDF to the
standard TMD-PDF, we also need the perturbative results of the above matrix
elements, for which we can choose on-shell quark external states. The
numerator of the r.h.s. of Eq. (34) has been calculated previously in
dimensional regularization (DR) and $\overline{\rm MS}$ scheme in Ref. Ebert
_et al._ (2019), whereas the denominator is independent of $\chi$ and reads

$\displaystyle h^{\text{TMD},\overline{\rm MS}}_{{\rm
pert},\gamma_{t}}(b_{0},z_{0},0;\mu)=1+\frac{\alpha_{s}C_{F}}{2\pi}\bigg{\\{}\frac{1}{2}+3\gamma_{E}-3\text{ln}2$
$\displaystyle\quad\quad+\frac{3}{2}\text{ln}[\mu^{2}(b_{0}^{2}+z_{0}^{2})]-2\frac{z_{0}}{b_{0}}\text{arctan}\frac{z_{0}}{b_{0}}\bigg{\\}}+{\cal
O}(\alpha_{s}^{2}).$ (35)

With this, we can also convert the SDR result to the $\overline{\rm MS}$
scheme via

$\displaystyle h^{\text{TMD},\overline{\rm
MS}}_{\chi,\gamma_{t}}(b,z,P_{z};\mu)=$ $\displaystyle
h^{\text{TMD},\overline{\rm MS}}_{{\rm
pert},\gamma_{t}}(b_{0},0,0;\mu)h^{\text{TMD},{\rm
SDR}}_{\chi,\gamma_{t}}(b,z,P_{z};\frac{1}{b_{0}}),$ (36)

where the $b_{0}$ dependence cancels.

As shown in Fig. 13, both results are consistent with each other and are close
to the perturbative result at short distance, and they deviate from each other
at long distance due to the non-perturbative effects from the intrinsic scales
$1/z$ and/or $1/b$ introduced in RI/MOM scheme.

##  V Summary

We have investigated the gauge fixing precision dependence for several
measurements, using 5 lattice spacings, in Landau gauge, Coulomb gauge and
$\xi$ gauge. Our results indicate that the linear divergences in quasi PDF
operators and quasi TMD-PDFs all comes from Wilson lines, and are independent
of the external state when the gauge fixing is precise enough. There is also
no linear divergence in Coulomb gauge quasi PDFs renormalized by short
distance matrix elements.

Thus, we verified that self-renormalization, which fits the hadron matrix
elements at different lattice spacings, RI/MOM renormalization, which uses the
quark matrix element and precise enough Landau gauge fixing, and the short-
distance renormalization, which divides by the corresponding Wilson line or
square root of the Wilson loop, can fully eliminate the UV divergence of the
quasi-PDF and quasi-TMD matrix elements, and make the renormalized ones finite
in the continuum limit. Obviously the combinations of RI/MOM and SDR schemes
(e.g. the short distance RI/MOM scheme Alexandrou  _et al._ (2023)) should
also have a finite continuum limit.

At the same time, even though all the approaches are consistent with the
perturbative results at short distances up to the discretization error, they
can be significantly different at long distances due to introduction of
spurious infrared physics in RI/MOM and similar schemes. Therefore, the proper
removal of the non-perturbative effects, such as self-renormalization, should
be crucial to obtain reliable matrix elements at long distances.

When the gauge fixing precision is restricted, the dependence of the results
on the gauge fixing residual $\delta^{F}$ can be roughly approximated by
$e^{-c\sqrt{\delta^{F}}}$, where $c$ grows quadratically with the Wilson line
length $z/a$ for both Wilson lines and quark matrix elements. Since the
accurate gauge fixing with $\delta^{F}$ can be extremely costly on large
lattices, this empirical formula can aid in estimating the systematic
uncertainty stemming from imprecise gauge fixing. For example, if we requires
the systematic uncertainty from the imprecise gauge fixing to be 1% at
$a=0.05$ fm and $z\sim 1$ fm $=20a$, then the $\delta^{F}$ should be around
$(10^{-2}/(z/a)^{2})^{2}\sim 6\times 10^{-10}$.

## Acknowledgement

We thank the MILC collaborations for providing us their gauge configurations,
and Xiang Gao, Xiangyu Jiang, Peter Petreczky, and Yong Zhao for useful
information and discussions. The calculations were performed using the Chroma
software suite Edwards and Joo (2005) with QUDA Clark  _et al._ (2010); Babich
_et al._ (2011); Clark  _et al._ (2016) through HIP programming model Bi  _et
al._ (2020). The numerical calculation were carried out on the ORISE
Supercomputer, and HPC Cluster of ITP-CAS. This work is supported in part by
NSFC grants No. 12293060, 12293062, 12293065 and 12047503, the science and
education integration young faculty project of University of Chinese Academy
of Sciences, the Strategic Priority Research Program of Chinese Academy of
Sciences, Grant No. XDB34030303 and YSBR-101, and also a NSFC-DFG joint grant
under Grant No. 12061131006 and SCHA 458/22.

## References

  * Martinelli  _et al._ (1995) G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, and A. Vladikas, Nucl. Phys. B445, 81 (1995), arXiv:hep-lat/9411010 [hep-lat] .
  * Mandula (1999) J. E. Mandula, Phys. Rept. 315, 273 (1999), arXiv:hep-lat/9907020 .
  * Wilson (1974) K. G. Wilson, Phys. Rev. D 10, 2445 (1974).
  * Mandula and Ogilvie (1987) J. E. Mandula and M. Ogilvie, Phys. Lett. B 185, 127 (1987).
  * Katz  _et al._ (1988) G. Katz, G. Batrouni, C. Davies, A. S. Kronfeld, P. Lepage, P. Rossi, B. Svetitsky, and K. Wilson, Phys. Rev. D 37, 1589 (1988).
  * Giusti  _et al._ (2001a) L. Giusti, M. L. Paciello, C. Parrinello, S. Petrarca, and B. Taglienti, Int. J. Mod. Phys. A 16, 3487 (2001a), arXiv:hep-lat/0104012 .
  * Lytle  _et al._ (2018) A. T. Lytle, C. T. H. Davies, D. Hatton, G. P. Lepage, and C. Sturm (HPQCD), Phys. Rev. D 98, 014513 (2018), arXiv:1805.06225 [hep-lat] .
  * Ji (2013) X. Ji, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Ji  _et al._ (2021) X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Rev. Mod. Phys. 93, 035005 (2021), arXiv:2004.03543 [hep-ph] .
  * Gao  _et al._ (2023) X. Gao, W.-Y. Liu, and Y. Zhao, (2023), arXiv:2306.14960 [hep-ph] .
  * Hart  _et al._ (2009) A. Hart, G. M. von Hippel, and R. R. Horgan (HPQCD), Phys. Rev. D 79, 074008 (2009), arXiv:0812.0503 [hep-lat] .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), Phys. Rev. D87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 832 (2005), [,832(2004)], arXiv:hep-lat/0409003 [hep-lat] .
  * Giusti  _et al._ (2001b) L. Giusti, M. L. Paciello, S. Petrarca, and B. Taglienti, Phys. Rev. D 63, 014501 (2001b), arXiv:hep-lat/9911038 .
  * Zwanziger (1986) D. Zwanziger, NATO Sci. Ser. B 141, 345 (1986).
  * Giusti (1997) L. Giusti, Nucl. Phys. B 498, 331 (1997), arXiv:hep-lat/9605032 .
  * Ji  _et al._ (2020) X. Ji, Y. Liu, and Y.-S. Liu, Nucl. Phys. B 955, 115054 (2020), arXiv:1910.11415 [hep-ph] .
  * Ji (2014) X. Ji, Sci. China Phys. Mech. Astron. 57, 1407 (2014), arXiv:1404.6680 [hep-ph] .
  * Ma and Qiu (2018a) Y.-Q. Ma and J.-W. Qiu, Phys. Rev. D98, 074021 (2018a), arXiv:1404.6860 [hep-ph] .
  * Ma and Qiu (2018b) Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. 120, 022003 (2018b), arXiv:1709.03018 [hep-ph] .
  * Izubuchi  _et al._ (2018) T. Izubuchi, X. Ji, L. Jin, I. W. Stewart, and Y. Zhao, Phys. Rev. D98, 056004 (2018), arXiv:1801.03917 [hep-ph] .
  * Zhang  _et al._ (2021) K. Zhang, Y.-Y. Li, Y.-K. Huo, A. Schäfer, P. Sun, and Y.-B. Yang ($\chi$QCD), Phys. Rev. D 104, 074501 (2021), arXiv:2012.05448 [hep-lat] .
  * Huo  _et al._ (2021) Y.-K. Huo _et al._ (Lattice Parton Collaboration (LPC)), Nucl. Phys. B 969, 115443 (2021), arXiv:2103.02965 [hep-lat] .
  * Zhang  _et al._ (2022) K. Zhang, X. Ji, Y.-B. Yang, F. Yao, and J.-H. Zhang ([Lattice Parton Collaboration (LPC)]), Phys. Rev. Lett. 129, 082002 (2022), arXiv:2205.13402 [hep-lat] .
  * Ji and Zhang (2015) X. Ji and J.-H. Zhang, Phys. Rev. D92, 034006 (2015), arXiv:1505.07699 [hep-ph] .
  * Ji  _et al._ (2018) X. Ji, J.-H. Zhang, and Y. Zhao, Phys. Rev. Lett. 120, 112001 (2018), arXiv:1706.08962 [hep-ph] .
  * Ishikawa  _et al._ (2017) T. Ishikawa, Y.-Q. Ma, J.-W. Qiu, and S. Yoshida, Phys. Rev. D96, 094019 (2017), arXiv:1707.03107 [hep-ph] .
  * Green  _et al._ (2018) J. Green, K. Jansen, and F. Steffens, Phys. Rev. Lett. 121, 022004 (2018), arXiv:1707.07152 [hep-lat] .
  * Liu  _et al._ (2020) Y.-S. Liu _et al._ (Lattice Parton), Phys. Rev. D101, 034020 (2020), arXiv:1807.06566 [hep-lat] .
  * Chen  _et al._ (2018) J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Phys. Rev. D97, 014505 (2018), arXiv:1706.01295 [hep-lat] .
  * Gracey (2003) J. A. Gracey, Nucl. Phys. B662, 247 (2003), arXiv:hep-ph/0304113 [hep-ph] .
  * Chang  _et al._ (2021) L. Chang, Y.-B. Liu, K. Raya, J. Rodríguez-Quintero, and Y.-B. Yang, (2021), arXiv:2105.06596 [hep-lat] .
  * Constantinou and Panagopoulos (2017) M. Constantinou and H. Panagopoulos, Phys. Rev. D 96, 054506 (2017), arXiv:1705.11193 [hep-lat] .
  * Stewart and Zhao (2018) I. W. Stewart and Y. Zhao, Phys. Rev. D97, 054512 (2018), arXiv:1709.04933 [hep-ph] .
  * Alexandrou  _et al._ (2017) C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and F. Steffens, Nucl. Phys. B923, 394 (2017), arXiv:1706.00265 [hep-lat] .
  * Alexandrou  _et al._ (2020) C. Alexandrou, K. Cichy, M. Constantinou, J. R. Green, K. Hadjiyiannakou, K. Jansen, F. Manigrasso, A. Scapellato, and F. Steffens, (2020), arXiv:2011.00964 [hep-lat] .
  * Lin  _et al._ (2020) H.-W. Lin, J.-W. Chen, and R. Zhang, (2020), arXiv:2011.14971 [hep-lat] .
  * Ebert  _et al._ (2019) M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP 09, 037 (2019), arXiv:1901.03685 [hep-ph] .
  * Constantinou  _et al._ (2019) M. Constantinou, H. Panagopoulos, and G. Spanoudes, Phys. Rev. D 99, 074508 (2019), arXiv:1901.03862 [hep-lat] .
  * Ebert  _et al._ (2020) M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP 03, 099 (2020), arXiv:1910.08569 [hep-ph] .
  * Alexandrou  _et al._ (2023) C. Alexandrou _et al._ , Phys. Rev. D 108, 114503 (2023), arXiv:2305.11824 [hep-lat] .
  * Spanoudes  _et al._ (2024) G. Spanoudes, M. Constantinou, and H. Panagopoulos, (2024), arXiv:2401.01182 [hep-lat] .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis Seattle, Washington, November 12-18, 2011_ (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, B. Jo, A. Strelchenko, M. Cheng, A. Gambhir, and R. Brower, (2016), arXiv:1612.07873 [hep-lat] .
  * Bi  _et al._ (2020) Y.-J. Bi, Y. Xiao, M. Gong, W.-Y. Guo, P. Sun, S. Xu, and Y.-B. Yang,  _Proceedings, 37th International Symposium on Lattice Field Theory (Lattice 2019): Wuhan, China, June 16-22 2019_ , PoS LATTICE2019, 286 (2020), arXiv:2001.05706 [hep-lat] .

## Appendix

###  V.1 The distribution of $\Delta^{G}(x)$

The quantity $\Delta^{G}(x)$, defined in Eq. (I) in Landau gauge, does not
reach exactly zero at every lattice site due to imperfect gauge fixing. As
illustrated in Fig. 14, the distribution is symmetrical and exhibits a long
tail. It is important to note that the scale of the abscissa is smaller for
higher precision. With a precision that is 100 times greater, the tail becomes
approximately 10 times shorter. However, the tail does not display a strong
dependence on lattice spacing. Additionally, we compute $\theta^{G}$, as
defined in Eq. (5), on the same configuration for comparison. The magnitude of
$\theta^{G}$ is directly proportional to that of $\delta^{F}$ in the part of
the ensemble that we are investigating.

Figure 14: $\Delta(x)=\Delta^{m}(x)T^{m}(x)$ defined in Eq. (I) on different
ensembles and with different required precision. Without loss of generality,
we just show the first component $\Delta^{1}(x)$ as illustration. The shape of
the distribution is similar on different ensembles and for different
precision.

###  V.2 Coulomb wall source

As our focus is solely on the UV divergence of the quasi-PDF operator, we
utilize the pion in the rest frame as the external state to enhance the
signal-to-noise ratio. Additionally, we have employed a Coulomb wall source to
further improve the signal. It is well-known that wall source calculations
include a gauge variation part, but with a sufficient number of
configurations, the gauge variant part will vanish. And one finds that
applying Coulomb gauge fixing improves the signal-to-noise ratio, but does not
lead to an obvious systematic deviation. In this work, we aim to verify
whether this holds true for non-local operator hadron matrix elements, and use
the unpolarized quasi PDF as illustration.

Figure 15: The ratio of quasi PDF pion matrix element in Coulomb gauge for
$\delta^{F}=10^{-6}$ and $\delta^{F}=10^{-7}$. They are calculated in the rest
frame with awall source.

As depicted in Fig. 15, the ratio remains close to one for various required
precisions, as long as the Wilson lines are not excessively long.
Consequently, we can conclude that the imprecise Coulomb gauge fixing in the
wall source calculation in this study does not change the central values.