RI/MOM renormalization of the quasi-PDF in lattice regularization

Parton Distribution Functions (PDFs) play a key role for most processes in high energy and hadron physics, ranging, e.g., from the search for New Physics to the detailed analysis of nucleon properties. Therefore, it is one of the most important tasks of lattice QCD to determine PDFs by first principle theory calculations. While the calculation of PDF Mellin moments, which can be expressed as matrix elements of local operators has a long, high-profile history, the calculation of their full $x$-dependence, which requires the evaluation of non-local correlators became only feasible in recent years. Among the different suggested approaches, the quasi-PDF method based on Large Momentum Effective Theory (LaMET) Ji (2013) has been widely used, because it is relatively direct and has no model dependence and proved to be quite successful. For example, iso-vector unpolarized PDFs have been calculated using different types of gauge and fermion actions down to the physical pion mass Alexandrou et al. (2018); Lin et al. (2018); Hua et al. (2020). A crucial topic which was not yet studied in great detail is the continuum extrapolation of quasi-PDFs. Controlling the continuum limit is already usually the main challenge for practicle lattice calculations of Mellin moments of PDFs. Because quasi-PDFs are far more delicate objects one has to expect that for them this task will be even more demanding, especially in combination with RI/MOM renormalization. At the same time, taking the continuum limit is indispensable for obtaining physical results, and therefore, intensifying such studies is the next logical step in the further development of the LaMET approach.

LaMET starts from the quasi-PDF operator $O_{\Gamma}(z)=\bar{\psi}(0)\Gamma U(0,z)\psi(z)$ with spatial gauge link $U(0,z)=\exp(-ig\int_{0}^{z}dz^{\prime}A_{z}(z^{\prime}))$, and connects its renormalized nucleon matrix element to PDFs through factorization theorems Ji (2013, 2014); Ma and Qiu (2018a, b); Izubuchi et al. (2018). However, the bare $O_{\Gamma}(z)$ in lattice regularizations contains an $1/a$ term

at 1-loop level Alexandrou et al. (2017) with lattice spacing $a$, and can deviate from its tree level value exponentially at either large $z$ or small $a$ when higher order effects are summed over. Thus, non-perturbative renormalization with accurate cancellation on the linear divergence is essential for the quasi-PDF, to ensure the existence of a finite continuum limit.

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 that the linear divergence just comes from the Wilson link and is independent of the external quark or hadron state. Therefore, 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). However, an accurate non-perturbative examination of this cancellation in lattice regularization is still missing.

To close this gap, we study the RI/MOM renormalized pion quasi-PDF matrix element in the rest frame, for 13 lattice spacings ranging from from 0.032 fm to 0.121 fm using different quark and gluon actions. Our results show that the cancelation deteriorates with decreasing lattice spacing, and that the RI/MOM method leaves a linearly divergent residue for quasi-PDFs. We also show that in the Landau gauge the interaction between the Wilson link and the external state results in a linear divergence which depends on the discretized fermion action.

When we consider the quasi-PDF nucleon matrix element in the moving frame, the gap $\delta m$ between the ground and first excited state decreases with increasing momentum, requireing a large source/sink separation $t_{\textrm{sep}}$ to eliminate excited state contaminations. At the same time, the signal to noise ratio also decreases exponentially with $t_{\textrm{sep}}$, reducing our ability to identify the lattice spacing dependence from data.

Using the pion rest frame can, therefore, be the better choice as it avoids this problem. It is known that the first excited state of the pion lies higher than 1GeV, and that the signal to noise ratio is almost independent of $t_{\textrm{sep}}$. Thus we can choose $t_{\textrm{sep}}$ to be half of the temporal lattice length $T$ ($T$ is mostly larger than 6 fm for the ensembles we use), and then obtain the ground state matrix element $h_{\pi}=\langle\pi|O_{\gamma_{t}}(z)|\pi\rangle$ with high accuracy:

where $O_{\pi}$ is the pion interpolating field. Note that the denominator includes both forward and backward two point functions. These are needed because of the possibility to loop around the lattice in the temporal direction for periodic boundary conditions.

In RI/MOM renormalization, we define the renormalization constant as Liu et al. (2020),

where $|q(p)\rangle$ is the off-shell quark state with external momentum $p$, and $Z_{q}$ is defined from the pion matrix element of a local vector current,

When the Landau gauge fixed volume ($V$) 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 error is much smaller than for $\tilde{Z}_{q}$ in lattice regularization as shown in Ref. Chang et al. (2021).

We can apply $Z_{\gamma_{t}}$ to the bare pion matrix element $h_{\pi,\gamma_{t}}(z)$ to obtain the non-perturbatively renormalized and normalized matrix element at a given RI/MOM scale $\mu$,

In order to check if the linear divergence is related to the fermion action, we compare results for two kinds of fermion actions in this work: the clover (CL) and overlap (OV) actions. The clover action is computationally relatively cheap and widely used in quasi-pdf calculations, while the overlap action is much more expensive but conserves chiral symmetry.

The clover fermion action is

with an additional “clover” term,

where $\hat{n}_{\mu}$ is the unit vector along the $\mu$ direction, $m^{\textrm{w}}_{q}-m^{\textrm{cri}}$ is the multiplicatively renormalizable bare quark mass, and $m^{\textrm{cri}}$ is the value of the bare quark mass for which the pion mass vanishes. $m^{\textrm{cri}}$ is ${\cal O}(\frac{\alpha_{s}}{a})$ at leading order for the Wilson action and always negative. It can be reducsed to ${\cal O}(\frac{\alpha^{2}_{s}}{a})$ (but usually still negative) with a clover coefficient $c_{\textrm{sw}}=1+{\cal O}(\alpha_{s})$. It can be further suppressed by applying gauge link smearing and/or fine-tuning of $c_{\textrm{sw}}$.

To eliminate $m^{\textrm{cri}}$ exactly one can use chiral fermions which satisfies the Ginsburg-Wilson relation $D_{\textrm{ov}}\gamma_{5}+\gamma_{5}D_{\textrm{ov}}=\frac{a}{\rho}D_{\textrm{ov}}\gamma_{5}D_{\textrm{ov}}$ Ginsparg and Wilson (1982). For example, Refs. Narayanan and Neuberger (1995); Neuberger (1998); Chiu (1999); Liu (2005) define overlap fermion by

where $-\rho$ should be smaller than $m^{\textrm{cri}}$, to make $D_{\textrm{ov}}$ to be the same as the standard Dirac operator in the continuum limit. The chiral fermion propagator is defined by $D_{\textrm{ov}}$,

where $D_{c}$ satisfies the relation $D_{c}\gamma_{5}+\gamma_{5}D_{c}=0$. Then, $m_{q}\rightarrow 0$ make the pion mass vanish without any fine-tuning.

We use three sets of ensembles:

(i) $N_{f}=2+1+1$ highly improved staggered fermions (HISQ) and one-loop Symanzik improved gauge fields from the MILC Collaboration Bazavov et al. (2013) at five lattice spacings (MILC12–MILC03).

(ii) $N_{f}=2+1$ domain wall fermions (DWF) and Iwasaki gauge fields from the RBC/UKQCD collaboration Blum et al. (2016) at three lattice spacings (RBC11–RBC06).

(iii) $N_{f}=2+1$ clover fermions and Luescher-Weisz (equivalent to Symanzik) gauge fields from the CLS collaboration Bruno et al. (2015) at five lattice spacings (CLS10–CLS04).

For the MILC/RBC gauge configurations we apply 1-step HYP smearing Hasenfratz and Knechtli (2001). For the overlap fermions we use $\rho=1.5/a$ and we use two clover coefficients: one is the tadpole improved tree level coefficient $c_{\mathrm{sw}}$ which is very close to one after the configuration is HYP smeared (CL), and the other is $c^{\prime}_{\mathrm{sw}}\sim 1.3$ which gives a critical quark mass $m^{\textrm{cri}}_{q}$ around zero (tuned-CL). The bare overlap quark mass $m_{q}^{\textrm{ov}}a$ is roughly proportional to the lattice spacing if we require the pion mass to be around 310 MeV, but $m^{\textrm{w}}_{q}a$ is always negative with $c_{\mathrm{sw}}\sim 1$ unless we enlarge $c_{\mathrm{sw}}$ to $\sim$ 1.3. For the CLS configuration, we simply use the unitary setup (without any smearing of the action) in the clover fermion calculations. Since most of the CLS ensembles use the open boundary condition except CLS10, we just use the 1/3 of the time slices in the middle of the lattice to avoid lattice artifacts for the volume source used for the RI/MOM renormalization constants.

We start with the pion matrix element and the MILC ensemble using the overlap action, to demonstrate the elimination of excited state contaminations. As shown in Fig. 1, the bare ratio $R_{\pi}(T/2,t,z\simeq 0.36$ fm$;~{}a)$ defined in Eq. (II) is independent of $t$ in the region $t-T/4\in$[-1,1] fm. Thus we average $R_{\pi}(T/2,t,z;a)$ in the region $t\in[T/8,3T/8]$ to get a precise estimate of the ground state matrix elements. We can also see that the ratio becomes exponentially smaller at smaller lattice spacings. Non-perturbative renormalization is essential to recover a reasonable continuum limit. A similar figure for the clover case can be found in our related work Huo et al. (2021) which uses the same data sets. It is worths emphasizing that if we take the ratio between the pion matrix elements $h_{\pi,\gamma_{t}}(z)$ obtained using the clover fermion and overlap fermion, then the ratio is consistent with 1 within 2 $\sigma$, or in other words with at most 3% difference for all $z<$ 1.0 fm as shown in Fig. 2. This agrees with the previous impression that the linear divergence is independent of the discretized fermion action.

Based on a symmetry analysis and 1-loop calculation, we can decompose the amputated Green’s function $\Lambda_{\gamma_{t}}(z,p)\equiv\langle q(p)|O_{\gamma_{t}}(z)|q(p)\rangle$ into the following Lorentz structures:

This is equivalent to what was proposed in Ref. Liu et al. (2020) but the form Eq. (10) is more natural to avoid $p_{z}$ to appear in a denominator. When $z=0$, the $\tilde{F}_{z}$ term vanishes and Eq. (10) becomes a matrix element of the local operator $\bar{q}\gamma_{\mu}q$. Both, the $\tilde{F}_{z}$ and $\tilde{F}_{p}$ terms vanish when $p_{t}=0$.

At the same time, setting $p_{z}=0$ leads to a vanishing imaginary part of $F_{t}$, as shown in the 1-loop calculations Constantinou and Panagopoulos (2017); Liu et al. (2020). Since $h_{\pi,\gamma_{t}}(z)$ is real in the rest frame, using the RI/MOM momentum with $p_{z}=p_{t}=0$ is the optimal choice to avoid complications from a complex renormalization constant as well as operator mixings. Based on the numerical calculation, we verified that the renormalization constant $Z_{\gamma_{t}}$ is real and the mixing between $O_{\gamma_{t}}(z)$ and $O_{\Gamma\neq\gamma_{t}}(z)$ is consistent with zero within very small uncertainties. Fig. 3 shows the real part of the operator mixing coefficients using clover fermions on the HISQ sea ensemble with $a$=0.121 fm. They are consistent with zero within 0.1% relative statistical uncertainties in all cases.

It is also popular to require all momentum components to be non-zero to suppress discretization errors, and use the p-slash projection as proposed in Ref. Stewart and Zhao (2018) when operator mixing is unavoidable. Since the conclusion is still unchanged with the p-slash projection, we place the related discussion in the appendix for the interested readers.

As for the RI/MOM renormalization constant in this work, we use the momentum $p=2\pi(5,5,0,0)/L$ on the HISQ sea ensembles and tune $p_{x,y}$ on the DWF sea ensembles to make $\mu=\sqrt{p^{2}}$ to be the same on all ensembles within 6%. As shown in the upper panel of Fig. 4 for the ratio of the renormalization constants,

with $\mu^{\prime}=\sqrt{-(p^{\prime})^{2}}=1.8$ GeV ($p^{\prime}=2\pi(3,3,0,0)/L$) on the HISQ sea ensembles using different valence fermion actions and lattice spacings, the ratio just changes by about 2% at $z$=1 fm regardless of the actions and lattice spacings used. Since the 6% fluctuation of $\mu$ on different ensembles is much smaller than the difference between $\mu^{\prime}$ and $\mu$, the systematic uncertainty due to this fluctuation will be much smaller than 2% at $z$=1 fm and thus negligible. At the same time, Fig. 4 also shows that $R_{\rm RI/MOM}$ converges at small lattice spacing and is independent of the discretized fermion actions (clover or overlap).

However, the residual $z$ dependence indicates that the $\mu$ and $z$ dependences of $Z(z,\mu)$ are not independent. Thus it is natural to be curious about its behavior when $\mu$ decreases towards the on-shell limit. Thus we also consider the ratio with $\mu^{\prime\prime}=\sqrt{-(p^{\prime\prime})^{2}}=0$ GeV using the clover valence fermion action on MILC ensembles. In the lower panel of Fig. 4, the $z$ dependence is 10 times stronger than for $\mu^{\prime}\simeq 1.8$ GeV, but the ratios at different lattice spacings are still consistent with each other. We conclude from the results in the two panels of Fig. 4 that the UV divergence of the RI/MOM renormalization constant is independent of the external off-shell momenta.

Another thing we need to mention here are the normalization conditions Eqs. (3–4) of the vector current. Based on the definition of the RI/MOM scheme, the quark self energy is defined from the charge conservation condition which requires the renormalized local vector charge to be unity, $h^{r}_{\pi,\gamma_{t}}=1$. In Fig. 5, we show the normalization factors $Z_{\gamma_{t}}(0,\mu)$ we obtained with overlap fermions on HISQ sea ensembles (red boxes) and DWF sea ensembles (blue dots), and also for clover fermions on HISQ sea ensembles (purple boxes). After performing simple fits for the $\alpha_{s}$ and $a^{2}$ corrections Bazavov et al. (2018), we can also show the lattice spacing dependence of the normalization factor with the ${\cal O}(a^{2})$ correction subtracted, using the red line for the overlap case and purple line for the clover case. It is clear that the normalization we used with the $\alpha_{s}$ correction approaches 1 smoothly and does not introduce any additional linear divergence.

In Ref. Huo et al. (2021), the Lattice Parton Collaboration parameterized the linear divergence as

where $M$ can be either the RI/MOM renormalization constant $Z$ or the pion matrix element $h$, and used our data to extract the linear divergence coefficients $k$ in different cases. As in Table IV-V of Ref. Huo et al. (2021) for the HYP smeared Wilson link, the $k$ is in the range k$\in[0.49,0.52]$ for the bare pion matrix elements with different valence and sea actions using $\Lambda_{QCD}=0.39$ GeV, but the $k$ in $Z$ is 0.55 for overlap fermions and $\sim$0.63 for clover fermions, on HISQ sea ensembles. The value of k in the RI/MOM renormalized pion matrix element is $\sim$ 0.05 and 0.13 for the overlap and clover fermion cases respectively, and much smaller than that in the bare pion matrix elements.

According to the lattice perturbative theory calculation Constantinou and Panagopoulos (2017), the linear divergence just comes from the Wilson link and is independent of the external quark state at the 1-loop level. This suggests that if there is any external state (like a quark state using different actions or the pion state) dependence in the quasi-PDF matrix element, it can only come from 2-loop or even higher level contributions, which is consistent with our observation as the residual linear divergence left after the cancelation is only at the 20% level.

We guess that the Landau gauge fixing we used for the quark matrix element introduces an additional linear divergence. As the Wilson link can have a gauge dependent logarithmic divergence at 1-loop level Constantinou and Panagopoulos (2017), we would expect it to introduce sub-leading linear divergences at higher loop levels. Thus we considered the ratio between the Wilson link in Coulomb and Landau gauge, and show the results in Fig. 10. The curves seem to converge at small lattice spacings, which suggests that the linear divergence of the Wilson link itself is unlikely to be gauge dependent, at least for Coulomb and Landau gauge fixing.

We also checked the correlation between the Wilson line and external states, by considering the following correlation

where $S(p)=\sum_{x}S(0,x)e^{ipx}$ is the momentum projected quark propagator. In the practical calculation we averaged over the reference points of all three expectation values in the right hand side of Eq. (13) independently to improve statistics. If the linear divergence only comes from the Wilson link and is independent of the external state, then $C_{US}(z)$ has to be free of any linear divergence, as confirmed by the 1-loop lattice perturbative theory calculation Constantinou and Panagopoulos (2017). However, as shown in Fig. 11 for the clover (upper panel) and overlap fermion (lower panel) on HISQ sea ensembles with 1-step of the HYP smearing of the link, we find the linear divergences to exist in both, the clover and overlap case. It is, however, much smaller in the latter. This is similar to what we saw for the RI/MOM renormalized pion quasi-PDF matrix element. However, since a similar test cannot be done for the Coulomb gauge due to the residual gauge degree of freedom along the temporal direction, we cannot yet verify whether the residual linear divergence we saw in $h^{r}_{\pi,\gamma_{t}}(z)$ and $C_{US}(z)$ are gauge dependent or not.

Eventually we consider the ratio

with the Wilson link in the Coulomb gauge, to test whether the linear divergence for the hadron matrix element is exactly the same as that for the Wilson link. As seen in Fig. 12, $h^{r,U}_{\pi,\gamma_{t}}(z)$ also has a similar lattice spacing and $z$ dependence pattern as we observed for $h^{r}_{\pi,\gamma_{t}}(z)$ and $C_{US}(z)$. This suggests that the linear divergence for the pion quasi-PDF matrix element is also not exactly the same as that of the Wilson line itself. Since we have shown that the linear divergence in the pion matrix element is independent of the fermion action (Fig. 2), and that in the Wilson line it is insensitive to whether one chooses Landau or Coulomb gauge fixing (Fig. 10), we will skip the discussion on the other combinations.

In this contribution we study systematically the continuum limit for quasi-PDFs with RI/MOM renormalization in the LaMET approach. We compare results for a variety of valence and sea quark fermion actions and lattice set-ups. We find that the traditional RI/MOM method cannot eliminate the characteristic linear divergences of quasi-PDFs completely, and that the remnants of this linear divergence blow up with decreasing lattice spacing. Obviously this greatly complicates control of the continuum limit.

We make a number of observations which can prove useful to clarify the origin of these remnants and develop strategies to cope with them:

The RI/MOM renormalized pion matrix elements in the rest frame include residual linear divergences from certain higher loop effects, so, higher loop calculation in lattice perturbation theory are needed.

Simulations with chiral fermions (overlap or domain-wall) are much less affected than those with clover fermions, so, chiral symmetry seems to play an important role.

Since the renormalization of operators should be frame independent, we believe that this conclusion should also apply to all quasi-PDF calculations in a moving frame. For a quasi-PDF calculation using an ensemble with $a=0.09$ fm and clover fermions, we estimate the systematic bias to be $\sim 5$% based on the correction at $z\sim$ 1.0 fm which is smaller than kinds of the systematic uncertainties in the practical quasi-PDF calculation; but it should increase to more than 50% at $a=0.03$ fm. Thus, removal of the residual linear divergence is crucial to obtain a meaningful continuum limit. For overlap fermions the corresponding bias could be 10% at 0.03 fm which might still be acceptable compared to other sources of uncertainty, but the fundamental problem is the same.

We believe that similar checks should be done for all LaMET calculations using bi-linear operators with Wilson link, e.g., quasi-PDFs, quasi-DAs, quasi-TMDs and so on. Comparison to calculations with large momentum for the nucleon suggests that the proposed check just requires the calculation of ${\cal O}(200)$ propagators on each ensemble, while providing essential information to avoid possible misinterpretations of the results obtained, which is obviously a very well justified investment. For the pseudo-PDFs for which a ratio of hadron matrix elements is used for renormalization, it is also essential to verify that the linear divergence is indeed completely canceled, independent of all details of the specific simulation, e.g., of the momenta used.

The outcome of many such tests should allow to better understand the origin of these residues and should provide clues for how to remove them systematically Huo et al. (2021). Let us stress that these residues become only sizable at small lattice spacings and that it is highly probable that the results obtained for not so fine lattices extrapolate smoothly to the continuum. So, existing results stay valid. However, the ultimate goal must be full control of the continuum limit and, therefore, this problem has to get solved.

The perturbative and non-perturbative study in lattice regularization beyond the 1-loop level is also essential, as a residual linear divergence is forbidden at the 1-loop level. We also confirmed the following expectations:

1) The linear divergences in the pion quasi-PDF matrix element with different fermion actions (overlap and clover) are the same. Combining with the consistency of the pion and nucleon matrix element shown in Ref. Huo et al. (2021), we conclude that the linear divergences in different hadrons are the same, regardless of the fermions we use.

2) The linear divergence in the Wilson link is indeed gauge independent up to certain lattice artifacts, based on our calculation in the Coulomb and Landau gauges.

3) There is a residual linear divergence due to the gluon exchange between the external state and the Wilson line, which is absent at 1-loop level. Such a linear divergence is action-dependent for a quark state in the Landau gauge, while it is action-independent in the hadron state.

We thank the CLS, MILC and RBC/UKQCD collaborations for providing us their gauge configurations, and Long-Cheng Gui, Xiangdong Ji, Keh-Fei Liu, Wei Wang, Jian-Hui Zhang and Yong Zhao for useful information and discussion. 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) and GWU-code Alexandru et al. (2012, 2011) through HIP programming model Bi et al. (2020). The numerical calculation has majorly been done on CAS Xiaodao-1 computing environment, and supported by Strategic Priority Research Program of Chinese Academy of Sciences, Grant No. XDC01040100, HPC Cluster of ITP-CAS, and Jiangsu Key Lab for NSLSCS. P. Sun is supported by Natural Science Foundation of China under grant No. 11975127, as well as Jiangsu Specially Appointed Professor Program. Y. Yang is supported by Strategic Priority Research Program of Chinese Academy of Sciences, Grant No. XDC01040100, XDB34030303, and XDPB15. A. Schäfer, P. Sun and Y. Yang are also supported by a NSFC-DFG joint grant under grant No. 12061131006 and SCHA 458/22.