Nonperturbative renormalization of nonlocal quark bilinears for quasi-PDFs on the lattice using an auxiliary field

Quasi-PDFs provide a path toward an ab initio calculation of parton distribution functions (PDFs) using lattice QCD. One of the problems faced in calculations of quasi-PDFs is the renormalization of a nonlocal operator. By introducing an auxiliary field, we can replace the nonlocal operator with a pair of local operators in an extended theory. On the lattice, this is closely related to the static quark theory. In this approach, we show how to understand the pattern of mixing that is allowed by chiral symmetry breaking, and obtain a master formula for renormalizing the nonlocal operator that depends on three parameters. We present an approach for nonperturbatively determining these parameters and use perturbation theory to convert to the MS-bar scheme. Renormalization parameters are obtained for two lattice spacings using Wilson twisted mass fermions and for different discretizations of the Wilson line in the nonlocal operator. Using these parameters we show the effect of renormalization on nucleon matrix elements with pion mass approximately 370 MeV, and compare renormalized results for the two lattice spacings. The renormalized matrix elements are consistent among the different Wilson line discretizations and lattice spacings.

Parton distribution functions (PDFs) describe the distribution of quarks and gluons inside a proton with respect to its longitudinal momentum. They are universal: the same PDFs appear in many different scattering processes, and they are phenomenologically determined from global fits to collider data. Except for their lowest Mellin moments, PDFs have resisted ab initio calculation. Lattice QCD can be used to calculate many properties of protons, but it is an inherently Euclidean space method, whereas PDFs are defined in Minkowski space via the matrix elements of operators with quark 1 creation and annihilation separated along the light cone. A possible solution to this problem was proposed by Ji [1]: compute quasi-PDFs using matrix elements of equal-time operators O Γ (x, ξ, n) ≡ψ(x + ξn)ΓW (x + ξn, x)ψ(x), where the ψ andψ are spatially separated by distance ξ in direction n and connected by a straight Wilson line W . From quasi-PDFs we can obtain PDFs via a matching formula [2,3], in the limit where the proton's momentum component p · n goes to infinity. One of the challenges in lattice calculations of quasi-PDFs is renormalization of O Γ , which is a nonlocal operator and is known to be power-law divergent [4,5]. This operator also appears in the related approach of pseudo-PDFs [6,7], and similar nonlocal operators are used to study transverse momentum-dependent PDFs [8][9][10]. The initial lattice studies of quasi-PDFs [3,[11][12][13] did not include a complete renormalization, but they did build the Wilson line using smeared gauge links, which 1 Here we focus on quark and antiquark PDFs rather than gluon PDFs.
has been shown in perturbation theory to reduce the power divergence [4]. Renormalization of O Γ was studied in one-loop lattice perturbation theory in Ref. [14], where it was found that chiral symmetry breaking allows O Γ to mix with O {/ n,Γ} . Numerical evidence for this mixing was also found in Ref. [10]. The study of nonperturbative renormalization was pioneered in Refs. [15,16], which used the Rome-Southampton method [17] to obtain a complex ξdependent renormalization factor (or matrix, when there is mixing) Z(ξ) in the RI-MOM scheme. This was then supplemented by a perturbative conversion [14] to the MS scheme used in phenomenology. The problems with this method are that a whole function, rather than a handful of parameters, must be determined, and that conversion at large ξ may occur outside the regime where perturbation theory is valid. Since the intermediate scheme fixes more than the minimal number of parameters, this means that nonperturbative information in correlation functions is lost, only to be recovered perturbatively in the conversion to MS.
In this work, we study the use of an auxiliary scalar, color triplet field ζ(ξ) defined only on the line x + ξn to simplify the renormalization of O Γ . In this approach, we replace correlation functions in QCD involving O Γ with correlation functions in the extended theory QCD+ζ involving the local color singlet bilinear φ ≡ζψ. This approach has been used long ago in the continuum [18,19].
On the lattice, for now we restrict n to point along one of the axes, n = ±μ, and use the action where ∇ and ∇ * are the forward and backward lattice covariant derivatives and a is the lattice spacing. This yields the bare propagator in fixed gauge background where m = a −1 log(1+am 0 ) and W is the simple product of lattice gauge links connecting x and x + ξn. Smeared gauge links can be used to construct W by using the same gauge links to define ∇ n . The mass term can not be forbidden by any symmetry, and corresponds to an O(a −1 ) counterterm. Using this propagator, we obtain for m = 0 and ξ > 0 the relation To obtain O Γ for ξ < 0, we reverse the direction in which ζ propagates, and use In addition to determining the counterterm m 0 , we must also renormalize the local composite operator φ. The lattice quark action may break chiral symmetry, in which case φ can mix with / nφ. We thus obtain the renormalization pattern We can also use projectors to form operators . This pattern leads to the form of the renormalized O Γ , for ξ = 0: O Γ can therefore be renormalized by determining three parameters: the linearly divergent m, the log-divergent Z φ , and the finite r mix . Since r mix is O(g 2 ), this is the same pattern at one-loop order as in Ref. [14]. At ξ = 0, O Γ is a local quark bilinear with different divergence structure, and should have a separate renormalization factor [20], which can be computed in the usual way; we use results from Ref. [21]. We also note that since the local operator φ is not flavor singlet, there is no mixing between quark and gluon quasi-PDFs even when O Γ is flavor singlet. In the latter case mixing will occur in the matching to PDFs. If we choose n =t and give ζ spin degrees of freedom (which do not couple in the action), then S ζ becomes the action for a static quark on the lattice [22,23]. In particular, φ is related to the static-light bilinears and we find relations between renormalization factors: Z stat The static quark theory also tells us how to remove O(a) lattice artifacts [24], which are present even if chiral symmetry is preserved on the lattice. In the continuum, the relation between O Γ and the static quark theory was previously discussed in Ref. [25].
We determine the renormalization parameters nonperturbatively using a variant of the Rome-Southampton method. In Landau gauge on N f = 4 twisted mass lattice ensembles [21], we compute the position-space ζ propagator, the momentum-space quark propagator S ψ (p), and the mixed-space Green's function for φ ± : The renormalization parameters, as well as the ζ and ψ field renormalizations can be determined by the conditions (11) and the RI -MOM/RI-SMOM condition for S ψ [17,26]. Eq. (9) is sensitive to m, whereas the others are constructed to eliminate dependence on it. These conditions define a two-parameter family of renormalization schemes at scale µ 2 = p 2 0 , that depend on the dimensionless quantities y ≡ |p 0 |ξ 0 and z ≡ p 0 · n/|p 0 |. We call this family of schemes RI-xMOM. Restricting to the kinematics p 0 ∝ n (i.e. z = 1), we have computed the conversion to the MS scheme at one-loop order, using dimensionally regularized perturbation theory [27]. We obtain in Landau gauge: where γ E is the Euler-Mascheroni constant and Ci is the cosine integral function, Ci(z) ≡ − ∞ z cos(t) t dt. For converting m to the MS scheme, we use the three-loop results for the static quark propagator from Refs. [28,29].
In Fig. 1, we show the quantity 2 in Eq. (9), for two different lattice spacings: a = 0.082 fm (β = 1.95) and a = 0.064 fm (β = 2.10). This is renormalized by adding m. Without smearing there is a significant difference between the two lattice spacings due to the linear divergence, but this is greatly reduced by applying one or five steps of HYP smearing [31], which also reduces the statistical uncertainty at large ξ. At small ξ/a, smearing distorts the shape and therefore we choose to impose our condition at ξ 0 ≈ 0.6 fm, where the shapes are similar. We then convert to the MS scheme at short distance by matching the result using unsmeared links on the fine ensemble to the perturbative result. From Eq. (11) an estimator for r mix can be isolated. In our data we see indications that this suffers from significant O(a) lattice artifacts. These could be remedied through Symanzik improvement [24], but instead we choose to focus on the helicity quasi-PDF, which is unaffected by mixing 3 and depends only on r 2 mix . We find that the latter is not much greater than 1% and can be neglected at the current level of precision.
The remaining factor Z φ depends on the kinematics that define our scheme, but the ratio between Z φ for different discretizations is scheme independent. We evaluate this ratio in the plateau region at large ξ and at small p (Fig. 2) to reduce lattice artifacts. Finally, we determine Z φ using unsmeared gauge links, converting to the MS scheme and evolving to the scale 2 GeV, as shown in Fig. 3. We find that the one-loop conversion factor is effective at removing much of the dependence on the scheme parameter |p|ξ, and the two-loop evolu- The green open squares are given in our family of schemes, the blue filled diamonds show the result from conversion to MS at scale |p| using Eq. (12), and the orange filled triangles with black outlines show the MS results evolved to the scale 2 GeV, using the two-loop anomalous dimension of the staticlight current [28,33,34].
tion removes most of the dependence on the scale |p|.
We apply these renormalization parameters to nucleon matrix elements computed on one N f = 2 + 1 + 1 twisted mass ensemble at each lattice spacing. The physical parameters are matched on the two ensembles: m π ≈ 370 MeV and p·n ≈ 1.85 GeV. The coarser ensemble was previously used in Refs. [3,13] and our methodology is similar to Ref. [13], including the use of momentum smearing [35] in the nucleon interpolating operator to obtain a good signal at large momentum. Figure 4 shows the effect of renormalization on the isovector helicity matrix element ∆h u−d , p, λ O / nγ5 p, λ ≡ū(p, λ )/ nγ 5 u(p, λ)∆h(ξ, p · n), (13) computed on the fine ensemble, for different link dis- cretizations. Without renormalization there is a significant disagreement between the different link types, and renormalization brings them into good agreement. In the renormalized matrix elements, we also see the benefit of smearing: without it, the statistical errors grow rapidly at large ξ and there is no useful signal for ξ > ∼ 10a. With smearing, we are able to see that the matrix elements return toward zero at large ξ. Five steps of HYP smearing also yields more precise data than one step, at large ξ. In Fig. 5, we compare the renormalized data with five steps of HYP smearing on the two ensembles. They are in excellent agreement, which suggests that the linear divergence is under control and discretization effects are not large.
Finally, we examine the helicity quasi-PDF, which is given by a Fourier transform of the matrix element: ∆q(x, p · n) ≡ p · n 2π dξ e −ixξp·n ∆h q (ξ, p · n). (14) Without renormalization, unsmeared links lead to a much broader distribution, as shown in Fig. 4 of Ref. [3]. We show our renormalized results on the fine ensemble in Fig. 6. Because the data become noisy at large ξ, we restrict the integral to |ξ| ≤ 16a (|ξ| ≤ 10a for the unsmeared case). When this restriction removes part of the signal, it leads to oscillations, which are clearly visible in the unsmeared case. In future studies improved results could be obtained by replacing the hard cutoff with a model for the large-ξ behavior of the matrix elements, or by applying one of the methods proposed in Refs. [36,37] for suppressing the contributions at large ξ. Ignoring the oscillations, we see that renormalization brings the data with different link discretizations into reasonably good agreement.
In this work, we have shown that the nonlocal problem of renormalizing lattice quasi-PDFs can be turned into a local problem by introducing an auxiliary field. The auxiliary field approach can also be applied to operators with staple-shaped gauge connection used for lattice studies of transverse momentum-dependent PDFs [8][9][10], where the mixing pattern will be different. Because this approach is closely connected with the static quark theory, it is possi-ble to make use of existing results from that theory such as the three-loop continuum calculation in Ref. [28]. We have contributed to the evidence that link smearing is a very beneficial technique for these calculations, as it leads to greatly reduced uncertainty at large ξ; this was also explained some time ago in the static quark theory [38]. There are several steps from the results shown here to a full calculation of PDFs: matching from quasi-PDFs to PDFs, control over the p · n → ∞ limit, use of a physical pion mass, as well as control over finite-volume and excited-state effects.
Note: after one of us presented a preliminary version of this work at a conference [39], we were made aware of an independent parallel effort based on the auxiliary field approach [40]. We also note another very recent article discussing renormalizability of quasi-PDFs, without using the auxiliary field approach [41].
We thank our colleagues in ETMC for a pleasant collaboration. Gauge fixing was performed using Fourieraccelerated conjugate gradient [42], implemented in GLU [43]. Calculations were performed using the Grid library [44] and the DD-αAMG solver [45] with twisted mass support [46]. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA [47] at Jülich Supercomputing Centre (JSC).