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

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,

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,

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$

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.

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

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}$).

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,

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,

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}$.

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.

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),

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,

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

As we did for Landau gauge, we plot the ratio

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.

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,

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.

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.

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.

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}$,

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$,

where

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.

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)),

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,

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$,

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).

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.

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.

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$:

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

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.

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

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,

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.

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.