Perturbative Renormalization of quasi-PDFs

TL;DR

This paper calculates the one-loop renormalization of gauge-invariant nonlocal fermion operators with Wilson lines in lattice QCD, addressing divergences, operator mixing, and proposing non-perturbative methods for divergence subtraction.

Contribution

It provides the first detailed one-loop lattice perturbation theory analysis of nonlocal fermion operators with Wilson lines, including mixing and divergence removal techniques.

Findings

Calculated finite mixing coefficients for operator subsets.

Identified linear and logarithmic divergences in renormalization.

Proposed non-perturbative methods to subtract divergences.

Abstract

In this paper we present results for the renormalization of gauge invariant nonlocal fermion operators which contain a Wilson line, to one-loop level in lattice perturbation theory. Our calculations have been performed for Wilson/clover fermions and a wide class of Symanzik improved gluon actions. The extended nature of such long-link operators results in a nontrivial renormalization, including contributions which diverge linearly as well as logarithmically with the lattice spacing, along with additional finite factors. On the lattice there is also mixing among certain subsets of these nonlocal operators; we calculate the corresponding finite mixing coefficients, which are necessary in order to disentangle individual matrix elements for each operator from lattice simulation data. Finally, extending our perturbative setup, we present non-perturbative prescriptions to extract the linearly divergent contributions.

Tables

Table 1. Table 1: Input parameters c 0 subscript 𝑐 0 c_{0} , c 1 subscript 𝑐 1 c_{1} , c 3 subscript 𝑐 3 c_{3} for selected gluon actions presented in this paper.

Action	$c_{0_{}}^{}$		$c_1$		$c_3$
Plaquette	1		0		0
Symanzik	5	/3	-1	/12	0
Iwasaki	3	.648	-0	.331	0

Table 2. Table 2: Numerical values of the coefficients e 1 subscript 𝑒 1 e_{1} - e 4 subscript 𝑒 4 e_{4} of the multiplicative renormalization functions and e 5 subscript 𝑒 5 e_{5} - e 6 subscript 𝑒 6 e_{6} of the mixing coefficients for Wilson, tree-level (TL) Symanzik and Iwasaki gluon actions.

Action	$e_1$	$e_2$	$e_3$	$e_4$	$e_5$	$e_6$
Wilson	24.3063	-19.9548	-2.24887	-1.39727	14.4499	-8.28467
TL Symanzik	19.8442	-17.2937	-2.01543	-1.24220	12.7558	-7.67356
Iwasaki	12.5576	-12.9781	-1.60101	-0.97321	9.93653	-6.52764

Equations

$$Z=e^{-cL/a},$$

$$S_F=$$

$${\hat{F}}_{\mu \nu }\equiv \frac{1}{8}(Q_{\mu \nu }-Q_{\nu \mu }),$$

$$Q_{\mu \nu }$$

$$\displaystyle S_{G}=\frac{2}{g_{0}^{2}}\Bigl{[}$$

$$+$$

$$c_0+8c_1+16c_2+8c_3=1.$$

$$C_0\equiv c_0+8c_1+16c_2+8c_3=1,$$

$$C_1\equiv c_2+c_3,$$

$$C_2\equiv c_1-c_2-c_3.$$

$${\mathcal{O}}_{\Gamma }\equiv \overline{\psi }(x)\Gamma \mathcal{P}{e}^{ig{\int }_0^z{A}_{\mu }(x+\zeta \hat{\mu })d\zeta }\psi (x+z\hat{\mu }),$$

$$\Gamma =\hat{1},{\gamma }^5,{\gamma }^{\nu },{\gamma }^5{\gamma }^{\nu },{\gamma }^5{\sigma }^{\nu \rho },{\sigma }^{\nu \rho }.$$

$$\begin{array}[]{lclrcl}{S}&\equiv&\mathcal{O}_{\hat{1}}&&&\\[2.15277pt] {P}&\equiv&\mathcal{O}_{\gamma^{5}}&&&\\[2.15277pt] {V}_{1}&\equiv&\mathcal{O}_{\gamma^{1}}&&&\\[2.15277pt] {V}_{\nu}&\equiv&\mathcal{O}_{\gamma^{\nu}}\,,&\nu&=&2,3,4\\[2.15277pt] {A}_{1}&\equiv&\mathcal{O}_{\gamma^{5}\gamma^{1}}&&&\\[2.15277pt] {A}_{\nu}&\equiv&\mathcal{O}_{\gamma^{5}\gamma^{\nu}}\,,&\nu&=&2,3,4\\[2.15277pt] {T}_{1\nu}&\equiv&\mathcal{O}_{\sigma^{1\nu}}\,,&\nu&=&2,3,4\\[2.15277pt] {T}_{\nu\rho}&\equiv&\mathcal{O}_{\sigma^{\nu\rho}}\,,&\nu,\rho&=&2,3,4\,.\end{array}$$

$$\left(\genfrac{}{}{0ex}{}{{\mathcal{O}}_{{\Gamma }_1}^R}{{\mathcal{O}}_{{\Gamma }_2}^R}\right)={\left(\begin{array}{cc}{Z}_{11}& Z_{12}\\ Z_{21}& Z_{22}\end{array}\right)}^{-1}\left(\genfrac{}{}{0ex}{}{{\mathcal{O}}_{{\Gamma }_1}}{{\mathcal{O}}_{{\Gamma }_2}}\right).$$

$$⟨{\psi }^R{\mathcal{O}}_{{\Gamma }_i}^R{\bar{\psi }}^R{⟩}_{\mathrm{amp}}=Z_{\psi }\sum _{j=1}^2(Z^{-1}{)}_{ij}⟨\psi {\mathcal{O}}_{{\Gamma }_j}\bar{\psi }{⟩}_{\mathrm{amp}},\psi =Z_{\psi }^{1/2}{\psi }^R,$$

$$Z_{ij}={\delta }_{ij}+g^2{z}_{ij}+\mathcal{O}(g^4),Z_{\psi }=1+g^2{z}_{\psi }+\mathcal{O}(g^4).$$

$$⟨{\psi }^R{\mathcal{O}}_{{\Gamma }_i}^R{\bar{\psi }}^R{⟩}^{DR,\overline{\mathrm{MS}}}={⟨{\psi }^R{\mathcal{O}}_{{\Gamma }_i}^R{\bar{\psi }}^R{⟩}^{LR,\overline{\mathrm{MS}}}|}_{a\to 0}.$$

$$\langle\psi^{R}\,{\cal O}^{R}_{\Gamma_{1}}\,\bar{\psi}^{R}\rangle^{DR,\,{\overline{\rm MS}}}-\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{LR}=g^{2}\,\Big{(}z_{\psi}^{LR,\,{\overline{\rm MS}}}-z_{11}^{LR,\,{\overline{\rm MS}}}\Big{)}\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{\rm tree}-g^{2}\,z^{LR,\,{\overline{\rm MS}}}_{12}\langle\psi\,{\cal O}_{\Gamma_{2}}\,\bar{\psi}\rangle^{\rm tree}+\mathcal{O}(g^{4}).$$

$$⟨{\psi }^R{\mathcal{O}}_{{\Gamma }_1}^R{\bar{\psi }}^R{⟩}^{DR,\overline{\mathrm{MS}}}=Z_{\psi }^{LR,\overline{\mathrm{MS}}}(Z_{11}^{LR,\overline{\mathrm{MS}}}{)}^{-1}⟨\psi {\mathcal{O}}_{{\Gamma }_1}\bar{\psi }{⟩}^{LR},$$

$$\langle\psi^{R}\,{\cal O}^{R}_{\Gamma_{1}}\,\bar{\psi}^{R}\rangle^{DR,\,{\overline{\rm MS}}}-\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{LR}=g^{2}\Big{(}z_{\psi}^{LR,\,{\overline{\rm MS}}}-z_{11}^{LR,\,{\overline{\rm MS}}}\Big{)}\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{\rm tree}+\mathcal{O}(g^{4}).$$

$${\rm Tr}\Big{[}\Lambda_{\Gamma_{i}}^{{{\rm RI}^{\prime}}}\,(\Lambda_{\Gamma_{j}}^{\rm tree})^{\dagger}\Big{]}_{q_{\nu}={\bar{q}}_{\nu}}={\rm Tr}\Big{[}\Lambda_{\Gamma_{i}}^{\rm tree}\,(\Lambda_{\Gamma_{j}}^{\rm tree})^{\dagger}\Big{]}\qquad\left(\ =12\delta_{ij}\ \right).$$

$${\Lambda }_{\Gamma }^{\mathrm{tree}}=\Gamma e^{i{q}_{\mu }z}.$$

$$\frac{1}{12}\,Z_{\psi}^{LR,{{\rm RI}^{\prime}}}\,\sum_{k=1}^{2}{(Z^{LR,{{\rm RI}^{\prime}}})^{-1}}_{ik}\,{\rm Tr}\Big{[}\Lambda_{\Gamma_{k}}\,(\Lambda_{\Gamma_{j}}^{\rm tree})^{\dagger}\Big{]}_{q_{\nu}=\bar{q}_{\nu}}=\delta_{ij}\,,$$

$$Z_{\psi}^{LR,{{\rm RI}^{\prime}}}=\frac{1}{12}\,{\rm Tr}\Big{[}S\,\left(S^{\rm tree}\right)^{-1}\Big{]}_{q^{2}={\bar{q}}^{2}}\,,$$

$$\left(\genfrac{}{}{0ex}{}{{\mathcal{O}}_{{\Gamma }_1}^{{\mathrm{RI}}^{\prime }}}{{\mathcal{O}}_{{\Gamma }_2}^{{\mathrm{RI}}^{\prime }}}\right)$$

$$\left(\genfrac{}{}{0ex}{}{{\mathcal{O}}_{{\Gamma }_1}^{\overline{\mathrm{MS}}}}{{\mathcal{O}}_{{\Gamma }_2}^{\overline{\mathrm{MS}}}}\right)$$

$${\mathcal{C}}^{\overline{\mathrm{MS}},{\mathrm{RI}}^{\prime }}=(Z^{LR,\overline{\mathrm{MS}}}{)}^{-1}\cdot (Z^{LR,{\mathrm{RI}}^{\prime }})=(Z^{DR,\overline{\mathrm{MS}}}{)}^{-1}\cdot (Z^{DR,{\mathrm{RI}}^{\prime }}).$$

$$\left(\genfrac{}{}{0ex}{}{{\Lambda }_{{\Gamma }_1}^{\overline{\mathrm{MS}}}}{{\Lambda }_{{\Gamma }_2}^{\overline{\mathrm{MS}}}}\right)=\frac{Z_{\psi }^{LR,\overline{\mathrm{MS}}}}{Z_{\psi }^{LR,{\mathrm{RI}}^{\prime }}}(Z^{LR,\overline{\mathrm{MS}}}{)}^{-1}\cdot (Z^{LR,{\mathrm{RI}}^{\prime }})\cdot (\genfrac{}{}{0ex}{}{{\Lambda }_{{\Gamma }_1}^{{\mathrm{RI}}^{\prime }}}{{\Lambda }_{{\Gamma }_2}^{{\mathrm{RI}}^{\prime }}})=\frac{1}{{\mathcal{C}}_{\psi }^{\overline{\mathrm{MS}},{\mathrm{RI}}^{\prime }}}({\mathcal{C}}^{\overline{\mathrm{MS}},{\mathrm{RI}}^{\prime }})\cdot (\genfrac{}{}{0ex}{}{{\Lambda }_{{\Gamma }_1}^{{\mathrm{RI}}^{\prime }}}{{\Lambda }_{{\Gamma }_2}^{{\mathrm{RI}}^{\prime }}}).$$

$${\int }_0^{z}d\zeta \int \frac{d{p}^D}{(2\pi {)}^D}\frac{e^{-i\zeta p_{\mu }}{p}_{\mu }}{p^2(p+q{)}^2}$$

$$\displaystyle\displaystyle\Lambda^{d1}_{\Gamma}\Big{|}_{1/\epsilon}$$

$$\displaystyle\Lambda^{d2+d3}_{\Gamma}\Big{|}_{1/\epsilon}$$

$$\displaystyle\Lambda^{d4}_{\Gamma}\Big{|}_{1/\epsilon}$$

$$\langle\psi\,{\cal O}_{\Gamma}\,\bar{\psi}\rangle^{DR}\Big{|}_{1/\epsilon}=g^{2}\,\lambda_{\Gamma}\,\Lambda_{\Gamma}^{\rm tree},\qquad\lambda_{\Gamma}=\frac{C_{f}}{16\,\pi^{2}}\,\frac{1}{\epsilon}\left(4-\beta\right)\qquad[C_{f}\equiv(N^{2}-1)/(2N)].$$

Full text

The extended nature of such ‘long-link’ operators results in a nontrivial renormalization, including contributions which diverge linearly as well as logarithmically with the lattice spacing, along with additional finite factors.

On the lattice there is also mixing among certain subsets of these nonlocal operators; we calculate the corresponding finite mixing coefficients, which are necessary in order to disentangle individual matrix elements for each operator from lattice simulation data. Finally, extending our perturbative setup, we present non-perturbative prescriptions to extract the linearly divergent contributions.

1 Introduction

Parton distribution functions (PDFs) provide important information on the quark and gluon structure of hadrons; at leading twist, they give the probability of finding a specific parton in the hadron carrying certain momentum and spin, in the infinite momentum frame. Due to the fact that PDFs are light-cone correlation functions, they cannot be computed directly on a Euclidean lattice. Nevertheless, there is an alternative approach, proposed by X. Ji [1], involving the computation of quasi-distribution functions, which are accessible in Lattice QCD. This formalism provides a promising means of studying quark distribution functions in nucleons, given that, for large momenta, one can establish connection with the physical PDFs through a matching procedure. Exploratory studies of the quasi-PDFs reveal promising results for the non-singlet operators for the unpolarized, helicity and transversity cases [2, 3, 4, 5].

A standard way of extracting quasi-distribution functions in lattice simulations involves computing hadronic matrix elements of certain gauge-invariant nonlocal operators; the latter are made up of a product of an anti-quark field at position $x$, possibly some Dirac gamma matrices, a path-ordered exponential of the gauge field (Wilson line) along a path joining points $x$ and $y$, and a quark field at position $y$ (particular cases are defined in Eq. (8) of the following Section). Given the extended nature of such operators, an endless variety of them, with different quantum numbers, can be defined and studied in lattice simulations and in phenomenological models. In pure gauge theories, prototype nonlocal operators are path-ordered exponentials along closed contours (Wilson loops); the contours may be smooth, but they may also contain angular points (cusps) and self-intersections.

The history of investigations of nonlocal operators in gauge theories goes back a long time, including seminal work of Mandelstam [6] and Polyakov [7], and encompassing several complementary viewpoints [8, 9]. In particular, the renormalization of Wilson loops was studied perturbatively, in dimensional regularization (DR) using a $D$-dimensional spacetime, for smooth contours [10] as well as for contours containing singular points [11]. Using arguments valid to all orders in perturbation theory, it was shown that smooth Wilson loops in DR are finite functions of the renormalized coupling, while the presence of cusps and self-intersections introduces logarithmically divergent multiplicative renormalization factors; at the same time, it was shown that other regularization schemes are expected to lead to further renormalization factors $Z$ which are linearly divergent with the dimensionful ultraviolet cutoff $a$:

[TABLE]

where $c$ is a dimensionless quantity and $L$ is the loop length.

Certain nonlocal operators have also been studied extensively via lattice simulations in the past; typical examples are products of open plaquettes at points $x$ and $y$, with varying orientations, joined with two Wilson lines running in both directions between $x$ and $y$ [12, 13]. Lattice results on matrix elements of such operators have found extensive use in the description of chromoelectric and chromomagnetic field correlations, in phenomenological models of the strong interactions (see, e.g., [14], [15], [16], etc.); however, the renormalization properties of these operators on the lattice need to be further explored.

At present, several aspects of Ji’s approach are being investigated: The matching between quasi-PDFs and physical PDFs [17, 18, 19, 20, 21], the relation with transverse momentum-dependent parton distributions (TMDs) [22, 23, 24, 25] and the extraction of the linear divergence through studies of the static quark potential [26, 27, 20]. Furthermore, there has been recent progress towards the investigation of the logarithmic divergences [27, 28, 29], the demonstration that the quasi-PDF extracted from a Euclidean correlation function is the same matrix element as that determined from the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula in Minkowski spacetime [30], as well as the quark-in-quark quasi-PDF in lattice perturbation theory [31]. There are several obstacles which need to be overcome before a transparent picture of PDFs can emerge via this approach; one such obstacle is clearly the intricate renormalization behavior, which is the object of our present study.

The paper is organized as follows: In Section 2 we formulate the problem, providing the definitions for the lattice action and for the operators which we set out to renormalize, along with the renormalization prescription. Section 3 contains our calculations, performed both in dimensional regularization and on the lattice; we address in detail new features appearing on the lattice, such as contributions which diverge linearly and logarithmically with the lattice spacing, and finite mixing effects allowed by hypercubic symmetry. We also provide a prescription for estimating the linear divergence using non-perturbative data and following arguments from one-loop perturbation theory. In Section 4 we summarize our results, and point out some open questions for future investigations.

2 Formulation

2.1 Lattice Actions

In the calculation we make use of the clover (Sheikholeslami-Wohlert) fermion action [32]; we allow the clover parameter, $c_{\rm SW}$, to be free throughout the calculation, in order to ensure wider applicability of our results. Using standard notation, this action reads:

[TABLE]

where $r$ is the Wilson parameter (henceforth set to 1), $f$ is a flavor index, ${\sigma_{\mu\nu}}=[\gamma_{\mu},\,\gamma_{\nu}]/2$ and $F_{\mu\nu}$ is the standard clover discretization of the gluon field tensor, defined through:

[TABLE]

where $Q_{\mu\nu}$ is given by the sum of the plaquette loops:

[TABLE]

We are interested in mass-independent renormalization schemes, and therefore we set the Lagrangian masses for each flavor, $m^{f}_{0}$ , to their critical value; for a one-loop calculation this corresponds to $m^{f}_{0}{=}0$. Such a choice simplifies the algebraic expressions but requires special treatment of potential IR singularities.

In the gluon sector we employ a 3-parameter family of Symanzik improved actions involving Wilson loops with 4 and 6 links, defined as [33]:

[TABLE]

The only restriction for the coefficients $c_{i}$ is a normalization condition, which ensures the correct classical continuum limit of the action:

[TABLE]

In this work we employ several values for the coefficients $c_{i}$, but for simplicity, we present numerical results for three choices widely used in numerical simulations. These are the Plaquette (Wilson), tree-level Symanzik-improved and Iwasaki actions; the corresponding values of the coefficients are shown in Table 1. Note that the one-loop Feynman diagrams which appear in our calculation do not involve pure gluon vertices and, thus, no vertex depends on $c_{i}$ ; furthermore, the gluon propagator depends only on three combinations of the Symanzik coefficients:

[TABLE]

Therefore, with no loss of generality we set $c_{2}=0$.

2.2 Definition of Operators

To establish notation, let us first write the operators we study in this work, which have the general form:

[TABLE]

with a Wilson line of length $z$ inserted between the fermion fields in order to ensure gauge invariance. In the limit $z\to 0$, Eq. (8) reduces to the standard ultra-local fermion bilinear operators. However, the calculation of the Green’s functions for $\mathcal{O}^{\mu}_{\Gamma}$ is for strictly $z\neq 0$ : the appearance of contact terms beyond tree level renders the limit $z\to 0$ nonanalytic.

We consider only cases where the Wilson line is a straight line along any one of the four perpendicular directions, which will be called $\mu$. Without loss of generality we choose 222Thus $\mu=1$ should be identified with the $z$ direction, which is conventionally chosen for the Wilson line in numerical studies. $\mu=1$. We perform our calculation for all independent combinations of Dirac matrices, $\Gamma$, that is:

[TABLE]

In the above, $\rho\neq\mu$ and we distinguish between the cases in which the index $\nu$ is in the same direction as the Wilson line ($\nu=\mu$), or perpendicular to the Wilson line ($\nu\neq\mu$). For convenience, the 16 possible choices of $\Gamma$ are separated into 8 subgroups, defined as follows:

[TABLE]

One might also consider an alternative definition for the “tensor” operators: $T^{\prime}_{\kappa\lambda}\equiv\mathcal{O}_{\gamma^{5}\sigma^{\kappa\lambda}}$ . Such a definition is clearly redundant if one employs a 4-dimensional regularization, such as the lattice, since the $T^{\prime}$ operators are just a renaming of the $T$ operators, and they will thus renormalize identically; as it turns out, even in dimensional regularization, where different treatments of $\gamma_{5}$ amount to different renormalization prescriptions (see Subsection 3.1.2), the renormalization of the above two sets of tensor operators remains identical.

2.3 Renormalization Prescription

We perform the calculation in both the dimensional (DR) and lattice (LR) regularizations, which allows one to extract the renormalization functions directly in the continuum ${\overline{\rm MS}}$-scheme. The setup of this process is extensively described in Refs. [34, 35] and is briefly outlined below. As is common practice, we will consider mass-independent renormalization schemes, so that fermion renormalized masses will be vanishing; for the one-loop lattice calculations this implies that the Lagrangian masses must be set to zero.

In the LR calculation we encounter finite mixing for some pairs of operators (see Subsection 3.2), and, thus, here we provide the renormalization prescription in the presence of mixing between two structures, $\Gamma_{1}$ and $\Gamma_{2}$, where one has, a $2\times 2$ mixing matrix 333All renormalization functions, generically labeled $Z$, depend on the regularization $X$ ($X$ = DR, LR, etc.) and on the renormalization scheme $Y$ ($Y$ = ${\overline{\rm MS}}$, RI*′*, etc.) and should thus properly be denoted as: $Z^{X,Y}$, unless this is clear from the context. ($Z$). More precisely, we find mixing within each of the pairs: $\{S,\,V_{1}\}$, $\{A_{2},\,T_{34}\}$, $\{A_{3},\,T_{42}\}$, $\{A_{4},\,T_{23}\}$, in the lattice regularization. In these cases, the renormalization of the operators is then given by a set of 2 equations:

[TABLE]

Once the mixing matrix $Z_{ij}$ is obtained through the perturbative calculation of certain Green’s functions, as shown below, it can be applied to non-perturbative bare Green’s functions derived from lattice simulation data, in order to deduce the renormalized, disentangled Green’s functions for each of the two operators separately. Such an application will be presented in a follow up publication using Twisted Mass fermions [36].

The one-loop renormalized Green’s function of operator $\mathcal{O}_{\Gamma_{i}}$ can be obtained from the one-loop bare Green’s function of $\mathcal{O}_{\Gamma_{i}}$ and the tree-level Green’s function of $\mathcal{O}_{\Gamma_{j}}$ ($j\neq i$); this can be seen starting from the general expression:

[TABLE]

where the renormalization matrix $Z$ and the fermion field renormalization $Z_{\psi}$ have the following perturbative expansion:

[TABLE]

Throughout this work, $g$ denotes the gauge coupling; the distinction between bare and renormalized coupling is immaterial for the one-loop perturbative calculation. All Green’s functions are intended to be amputated, and thus the label “${\rm amp}$” will be dropped from this point on.

Once the ${\overline{\rm MS}}$ renormalized Green’s functions have been computed in DR (see Section 3.1), the condition for extracting $Z^{LR,\,{\overline{\rm MS}}}_{11}$ and $Z^{LR,\,{\overline{\rm MS}}}_{12}$ is simply the requirement that renormalized Green’s functions be regularization independent:

[TABLE]

Substituting the right-hand side of the above relation by the expression in Eq. (12), there follows:

[TABLE]

The Green’s functions on the left-hand side of Eq. (15) are the main results of this work, where $\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{DR,\,{\overline{\rm MS}}}$ is the renormalized Green’s function for ${\cal O}_{\Gamma_{1}}$ which has been computed in dimensional regularization (Eqs. (72)-(76)) and renormalized using the ${\overline{\rm MS}}$-scheme, while $\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{LR}$ is the bare Green’s function of ${\cal O}_{\Gamma_{1}}$ in LR. The difference of the aforementioned Green’s functions is polynomial in the external momentum (of degree 0, in our case, since no lower-dimensional operators mix); in fact, verification of this property constitutes a highly nontrivial check of our calculations. Thus, Eq. (15) is an appropriate definition of the momentum-independent renormalization functions, $Z_{11}$ and $Z_{12}$. Note that in the absence of mixing ($Z_{12}=Z_{21}=0$), Eqs. (12), (15) reduce to:

[TABLE]

[TABLE]

Non-perturbative evaluations of the renormalization functions cannot be obtained directly in the ${\overline{\rm MS}}$ scheme; rather, one may calculate them in some appropriately defined variant of the RI*′* (“modified regularization-invariant”) scheme, and then introduce the corresponding conversion factors between RI*′* and ${\overline{\rm MS}}$. Here we propose a convenient RI*′* scheme which can be applied non-perturbatively, similar to the case of the ultra-local fermion composite operators, with due attention to mixing. Defining, for brevity: $\Lambda_{\Gamma_{i}}=\langle\psi\,{\mathcal{O}}_{\Gamma_{i}}\,\bar{\psi}\rangle$ , and denoting the corresponding renormalized Green’s functions by $\Lambda_{\Gamma_{i}}^{{{\rm RI}^{\prime}}}$ , we require:

[TABLE]

The factor of 12 above originates from the fact that the trace acts on both Dirac and color indices. The momentum of the external fermion fields is denoted by $q_{\nu}$, and the four-vector ${\bar{q}}_{\nu}$ denotes the RI*′* renormalization scale. We note that the magnitude of $\bar{q}$ alone is not sufficient to specify completely the renormalization prescription: Different directions in $\bar{q}$ amount to different renormalization schemes, which are related among themselves via finite renormalization factors. In what follows we will select RI*′* renormalization scale 4-vector to point along the direction $\mu=1$ of the Wilson line: $(\bar{q},0,0,0)$. The tree-level functions $\Lambda_{\Gamma_{j}}^{\rm tree}$ are given by:

[TABLE]

Using Eq. (12) we express Eq. (18) in terms of bare Green’s functions, obtaining:

[TABLE]

[TABLE]

where $S$ is the bare quark propagator and $S^{\rm tree}$ is its tree-level value; the one-loop computation of $S$ can be found, e.g., in Ref. [37].

Eq. (20) amounts to four conditions for the four elements of the matrix $Z^{LR,{{\rm RI}^{\prime}}}$. As it was intended, it lends itself to a non-perturbative evaluation of $Z^{LR,{{\rm RI}^{\prime}}}$, using simulation data for $\Lambda_{\Gamma_{k}}$ .

Converting the non-perturbative, RI*′*-renormalized Green’s functions $\Lambda_{\Gamma_{i}}^{{{\rm RI}^{\prime}}}$ to the ${\overline{\rm MS}}$ scheme relies necessarily on perturbation theory, given that the very definition of ${\overline{\rm MS}}$ is perturbative in nature. We write:

[TABLE]

The conversion factor $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$ is a $2{\times}2$ matrix in this case; it is constant ($q$-independent) and stays finite as the regulator is sent to its limit ($a\to 0$ for LR, $D\to 4$ for DR). Most importantly, its value is independent of the regularization:

[TABLE]

Thus, the evaluation of $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$ can be performed in DR, where evaluation beyond one loop is far easier than in LR; this, in a nutshell, is the advantage of using the RI*′* scheme as an intermediary. A further simplification originates from the fact that the DR mixing matrices $Z^{DR,{\overline{\rm MS}}}$ and $Z^{DR,{{\rm RI}^{\prime}}}$ are both diagonal; as a result, $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$ turns out to be diagonal as well. We stress that in the case of ultra-local operators, the conversion factor depends on the renormalized coupling and the ratio of the ${\overline{\rm MS}}$ over the RI*′* renormalization scales ($\bar{\mu}/\bar{q}$); for the Wilson line operators on the other hand, the conversion may (and, in general, will) depend on the length of the Wilson line and on the individual components of the RI*′* renormalization-scale four-vector, through the dimensionless quantities $z{\bar{q}}_{\nu}$ . Finally, let us also point out that, just as in the case of ultra-local operators, the non-perturbative evaluation of RI*′* renormalization functions is performed at nonzero renormalized quark masses $m_{q}$ ; this would, in principle, necessitate that the perturbative evaluation of $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$ also be performed at the same nonzero values of $m_{q}$ . However, given the near-critical values for light quark masses employed in present-day simulations, $m_{q}/\bar{\mu}\ll 1$, and given the smooth dependence of $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$ on $m_{q}$ , only imperceptible changes are expected by setting $m_{q}\to 0$.

Once the conversion factor is computed, the Green’s functions $\Lambda_{\Gamma_{i}}^{{{\rm RI}^{\prime}}}$ can be directly converted to the ${\overline{\rm MS}}$ scheme through:

[TABLE]

The fermion field conversion factor: $\mathcal{C}_{\psi}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}\equiv Z_{\psi}^{LR,{{\rm RI}^{\prime}}}/Z_{\psi}^{LR,{\overline{\rm MS}}}=Z_{\psi}^{DR,{{\rm RI}^{\prime}}}/Z_{\psi}^{DR,{\overline{\rm MS}}}$ is a finite function of the renormalized coupling constant, and its value is known well beyond one loop [38].

In non-perturbative studies of Green’s functions with physical nucleon states, performed via lattice simulations, the external states are normalized in a way which does not involve the quark field renormalization $Z_{\psi}$ ; thus, the only conversion factor necessary in this case is $\mathcal{C}^{{\overline{\rm MS}},{{\rm RI}^{\prime}}}$.

3 Calculation - Results

The Feynman diagrams that enter our one-loop calculations are shown in Fig. 1, where the filled rectangle represents the insertion of any one of the nonlocal operators $\mathcal{O}_{\Gamma}$ with a Wilson line of length $z$. Diagram d1 contains the 0-gluon vertex of the operator ($\mathcal{O}(g^{0})$ of Eq. (8)), whereas diagrams d2-d3 (d4) contain the corresponding 1-gluon (2-gluon) vertex. These diagrams will appear in our calculations in both LR and DR, since all vertices are present in both regularizations, and since even the “tadpole” diagram d4 does not vanish in DR, by virtue of the nonlocal nature of $\mathcal{O}_{\Gamma}$ . However, the LR calculation is much more challenging: The vertices of $\mathcal{O}_{\Gamma}$ are more complicated, and extracting the singular parts of the Green’s functions is a more lengthy and subtle procedure.

3.1 Dimensional Regularization

Let us first recall some of the essentials of dimensional regularization and of the ${\overline{\rm MS}}$ scheme. The computation is performed in $D$ Euclidean spacetime dimensions, where $D=4-2\epsilon$ and $\epsilon$ is the regularizing parameter. Bare $n$-loop Green’s functions in DR will be Laurent series in $\epsilon$, of the form $\sum_{i=-n}^{\infty}c_{i}\,\epsilon^{i}$ ; renormalization must eliminate all poles (negative powers) in $\epsilon$, before the limit $D\to 4$ can be taken. The simplest scheme for the elimination of divergences is the modified minimal subtraction (${\overline{\rm MS}}$) scheme according to which the renormalization functions are defined to only remove poles in $\epsilon$, without any finite parts 444Together with this operation one must also express every occurence of the dimensionful scale $\mu$ in terms of the ${\overline{\rm MS}}$ renormalization scale $\bar{\mu}\equiv\mu\left(4\pi/e^{\gamma_{E}}\right)^{1/2}$, where $\gamma_{E}$ is the Euler constant and $\mu$ appears in bare Green’s functions by virtue of the relation between the $D$-dimensional bare coupling $g$ and the renormalized coupling $g^{R}$: $g=\mu^{\epsilon}\,Z_{g}\,g^{R}$.. ${\overline{\rm MS}}$ has become a reference scheme, and although other schemes, such as RI*′, are more suitable for non-perturbative calculations, an appropriate conversion factor is applied to reach ${\overline{\rm MS}}$. In this Section we present our results for the renormalization functions in the ${\overline{\rm MS}}$ scheme, and provide the conversion factor between RI′*, as defined in Eq. (18), and ${\overline{\rm MS}}$.

For the extraction of UV divergences and determination of the poles in $\epsilon$, we follow the standard procedure of introducing Feynman parameters which allow us to perform the $D$-dimensional loop integrals [39]. However, unlike ordinary massless 2-point Green’s functions, in which the dependence on the external momentum $q$ is dictated purely on dimensional grounds, the results of integration are now considerably more complicated; this is due to the appearance of both $q$ and the length $z$ in the integrands, along with a nontrivial dependence on the preferred direction of the Wilson line. Just to give an example, we present one of the simpler integrals which appear in the one-loop calculation:

[TABLE]

where $p$ is the loop momentum, $q$ is the external momentum, $x$ is a Feynman parameter to be integrated over, and $K_{n}$ is the modified Bessel function of the second kind. In the first line of Eq. (25), integration was performed over the $D{-}1$ components of $p$ perpendicular to $p_{\mu}$ ; subsequently, a Laurent expansion was made, keeping terms up to ${\cal O}(\epsilon^{0})$. We see that the $1/\epsilon$ pole part of this expression has a constant coefficient; even though a potential $z$-dependence would not come into conflict with renormalizability, it cannot arise, since neither $q$ nor any renormalization scale can appear in the pole part and, therefore, no dimensionless combination can be formed out of $z$. The finite part, on the other hand, clearly exhibits a nontrivial dependence on the dimensionless quantities: $zq$, $zq_{\mu}$, in addition to the standard logarithmic dependence on $q/\bar{\mu}$.

In a similar way we manipulate all integrals appearing in the computation in dimensional regularization. In Subsection 3.2 we will discuss the evaluation of the lattice counterpart of the above integral and describe the highly intricate process for the proper extraction of the UV divergences.

3.1.1 Renormalization Functions

The perturbative calculation has been performed in an arbitrary covariant gauge in order to see first-hand the gauge invariance of the renormalization functions; this serves as a consistency check of our calculation. The gauge fixing parameter, $\beta$, is defined such that $\beta{=}0(1)$ corresponds to the Feynman (Landau) gauge. We find that $1/\epsilon$ terms arise from diagrams d2, d3 and d4, and thus d1 does not contribute to the ${\overline{\rm MS}}$ renormalization function of the operators under study. However, d1 contributes to the renormalized Green’s functions and to the conversion factors from RI*′* to ${\overline{\rm MS}}$. Below we present the ${\cal O}(1/\epsilon)$ contributions to $\Lambda_{\Gamma}=\langle\psi\,{\mathcal{O}}_{\Gamma}\,\bar{\psi}\rangle$ from each of the diagrams, including all combinatorial factors:

[TABLE]

and thus:

[TABLE]

Note that to one-loop level in the DR calculation the pole parts are multiples of the tree-level values, which indicates no mixing between operators of equal or lower dimension. Also, diagrams d2 and d3, the so-called ‘sail’ diagrams, are symmetric and give the same contribution to the total Green’s function. Another important characteristic of the ${\cal O}(g^{2})/\epsilon$ contributions is that they are operator independent, in terms of both the Dirac structure and the length of the Wilson line, $z$. Using Eq. (29) and $Z_{\psi}$ in DR [40]:

[TABLE]

the ${\overline{\rm MS}}$ condition to one loop reads:

[TABLE]

As expected from gauge invariance, we find a gauge independent renormalization function for the operators of Eq. (8):

[TABLE]

in agreement with Refs.  [41, 42]. Gauge invariance is not guaranteed in all schemes, since the Green’s functions on which the renormalization condition is imposed typically contain gauge variant renormalized external fields. Since the ${\overline{\rm MS}}$ scheme removes only the divergences, which are universal for all (gauge dependent and independent) Green’s functions, it has to be gauge invariant.

While the independence of $Z^{DR,\,{\overline{\rm MS}}}_{\Gamma}$ from the Dirac matrix insertion $\Gamma$ is a feature valid to one-loop level, its independence from the length of the Wilson line $z$ is expected to hold to all orders in perturbation theory; this, in essence, is due to the fact that the most dominant pole at every loop can depend neither on the external momenta nor on the renormalization scale, thus there is no dimensionless $z$-dependent factor that could appear in the pole part.

3.1.2 Conversion Factors

The Green’s functions in DR are also very useful for the computation of the conversion factors between different renormalization schemes, and here we are interested in the RI*′* scheme defined in Eq. (20). To one-loop level, the conversion factor is simplified to:

[TABLE]

which we have computed for all operators shown in Eqs. (10). Note that our one-loop calculations do not depend on the prescription which one adopts for extending $\gamma_{5}$ to $D$ dimensions (see, e.g., Refs. [43, 44, 45, 46, 47, 48] for a discussion of four relevant prescriptions and some conversion factors among them); this is because (anti-)commutation relations among $\gamma_{5}$ and $\gamma_{\nu}$ appear only in contributions which are finite as $\epsilon\to 0$. Also, the conversion factor is the same for each of the following pairs of operators: Scalar and pseudoscalar, vector and axial “parallel” ($V_{1},\ A_{1}$), vector and axial “perpendicular” ($V_{2,3,4},\ A_{2,3,4}$), as well as for the tensor with and without $\gamma^{5}$ ; furthermore, all components of $T_{\nu\rho}$ will have the same conversion factor, regardless of whether $\nu$ or $\rho$ are parallel or perpendicular to the Wilson line. The general expression for $\mathcal{C}_{\Gamma}$ are shown in Eqs. (34)-(37) for general gauge fixing parameter. They are expressed compactly in terms of the quantities $F_{1}(\bar{q},z)$ - $F_{5}(\bar{q},z)$ and $G_{1}(\bar{q},z)$ - $G_{5}(\bar{q},z)$, which are integrals over modified Bessel functions of the second kind, $K_{n}$. These integrals are presented in Eqs. (62)-(71) of Appendix A. Just as in the example of Eq. (25), the conversion factors depend on the dimensionless quantities $z\bar{q}$ and $\bar{q}/\bar{\mu}$. The RI*′* and ${\overline{\rm MS}}$ renormalization scales ($\bar{q}$ and $\bar{\mu}$, respectively) have been left independent.

The Green’s functions of operators with a Wilson line are complex, a property which is also valid for the non-perturbative matrix elements between nucleon states (see, e.g. Ref. [5]); this is also propagated to the conversion factor, as can be seen in the following equations.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $q{\equiv}\sqrt{q^{2}}$.

In Appendix B we present the ${\overline{\rm MS}}$-renormalized Green’s functions for each Wilson line operator. Starting from these Green’s functions, the conversion factors Eqs. (34)-(37) can be derived in a straightforward manner, but have been included for ease of reference.

Note that for a scale of the form $(\bar{q},0,0,0)$ the one-loop Green’s functions are a multiple of the tree-level value of the operator under consideration. The conversion factors take the form:

[TABLE]

where $F_{\Gamma}(\bar{q}z)$ are defined in Eqs. (50) - (54).

We remind the reader that the conversion factors as defined in Eq. (23) are to be multiplied by the $Z^{{\overline{\rm MS}}}$ to give $Z^{{\rm RI}^{\prime}}$. Alternatively, one may obtain $Z^{{\overline{\rm MS}}}$ by multiplying $Z^{{\rm RI}^{\prime}}$ with $\mathcal{C}_{\Gamma}$ provided that $g^{2}\to-g^{2}$, which is valid to one-loop level.

It is interesting to plot the conversion factors for the cases used in simulations, that is, $\mathcal{C}_{V_{1}(A_{1})}$ and $\mathcal{C}_{T}$. For convenience we choose the coupling constant and the RI*′* momentum scale to match the ensemble of twisted mass fermions employed in Ref. [5]: $g^{2}{=}3.077$, $a{=}0.082$fm, lattice size: $32^{3}\times 64$ and $a\bar{q}{=}\frac{2\pi}{32}(\frac{n_{t}}{2}{+}\frac{1}{4},0,0,n_{z})$, for $n_{t}{=}8$ and $n_{z}{=}4$ (the nucleon is boosted in the $z$ direction). The ${\overline{\rm MS}}$ scale is set to $\bar{\mu}=2{\rm GeV}$. The conversion factors are gauge dependent and we choose the Landau gauge which is mostly used in non-perturbative renormalization. Since the RI*′* scale is given in lattice units, $a\bar{q}$, we also rescale the length of the Wilson line with the lattice spacing, that is $z/a$. In Fig. 2 we plot the real (left panel) and imaginary (right panel) parts of $\mathcal{C}_{V_{1}(A_{1})}$ and $\mathcal{C}_{T}$, as a function of $z/a$. We remind the reader that the results of Eq. (38) are only valid for $z{\neq}0$. Thus, the points shown in the plot at $z{=}0$ (open symbols) have been extracted from Ref. [49], and is a real function. In particular, it is exactly one for the scale and scheme independent vector and axial operators.

In the plot we allow $z$ to take all possible values up to half the lattice size, for both forward and backward directions of the Wilson line. One observes that the real part is symmetric with respect to $z{=}0$, while the imaginary part is antisymmetric. Additionally, for large values of $z$ the dependence of the conversion factor on the choice of operator becomes milder. However, this behavior is not granted at higher loops, where a more pronounced dependence on the operator may arise.

3.2 Lattice Regularization

We now turn to the evaluation of the lattice-regularized bare Green’s functions $\langle\psi\,{\cal O}_{\Gamma_{1}}\,\bar{\psi}\rangle^{LR}$ ; this is a far more complicated calculation, as compared to dimensional regularization, because the extraction of the divergences is more delicate. The main task is to write lattice expressions in terms of continuum integrals, such as those in Appendix A, plus additional terms which are lattice integrals, independent of external momentum $q$; however, in contrast to the case of local operators, these additional terms are expected to have a nontrivial dependence on $z$. Thus, the renormalized Green’s functions stemming from the lattice can be made to coincide with the continuum ones, shown in Appendix B, once an appropriate $q$-independent (but $z$-dependent!) renormalization is applied.

The one-loop Feynman diagrams in LR are the same as in DR (Fig. 1). Further, diagram 1 gives exactly the same contribution as in DR; this was to be expected by the fact that the latter is finite as $\epsilon\to 0$, and thus the limit $a\to 0$ can be taken right from the start, with no lattice corrections.

Once the loop momentum $p$ is rescaled to fit the boundaries of the Brillouin zone, $p\to p/a$, lattice divergences manifest themselves as IR divergences in the external momentum $q\to 0$. We turn next to this issue, before giving our results in LR.

3.2.1 Isolation of IR divergences

One of the most laborious tasks in the lattice computation is the extraction of IR divergences, a process much more complex than in dimensional regularization. For demonstration purposes we present the lattice counterpart of the integral discussed in DR (Eq. (25)), in order to point to the necessity of applying novel techniques for its evaluation. In particular we developed a procedure for the isolation of IR divergences, somewhat in the spirit of the standard procedure of Kawai et al. [50], in which one subtracts and adds to the original integrand its naïve Taylor expansion. In the present case, the steps required are more complicated, and are presented below in a nutshell.

A lattice expression analogous to Eq. (25) is:

[TABLE]

where $\hat{k}^{2}\equiv 4\sum_{\rho}\sin^{2}(p_{\rho}/2)$, $\mu$ is the direction of the Wilson line, $\rho$ is in one of the four directions, and $n\equiv z/a$; $I_{lat}$ appears in the two ‘sail’ diagrams. Below we present a schematic way of the process that we developed for the evaluation of $I_{lat}$, which is based on a series of additions and subtractions. To avoid complicated expressions we present the process in the Feynman gauge, even though in our calculation we have kept a general gauge parameter $\beta$.

We write: $I_{lat}=(I_{lat}-I_{1})+(I_{1}-I_{2})+(I_{2}-I_{3})+I_{3}$, where:

[TABLE]

Here, $\overline{p}$ is the loop momentum $p$ with its $\mu$-component set to 0. In Eq. (42) $r$ represents the radius of a sphere over which we integrate, and has an arbitrary value. Each of the two terms in Eq. (42) can be calculated by integration in spherical coordinates, leading to the same expressions as those found in DR in terms of Bessel functions, plus $q$-independent terms. Naturally, the integrals over the two terms of Eq. (42) depend on the parameter $r$. However, once all parts of $I_{lat}$ are combined, the final result is independent of $r$; this cancellation is highly nontrivial and it provides an important cross-check of our calculation.

3.2.2 Multiplicative Renormalization and Mixing

Despite the complexity of the bare Green’s functions, their difference in DR and LR is necessarily polynomial in the external momentum (of degree 0 in this case), which leads to a prescription for extracting $Z_{\mathcal{O}}^{LR,{\overline{\rm MS}}}$ without an intermediate (e.g., RI*′*-type) scheme.

By analogy with closed Wilson loops [10] in regularizations other than DR, we find a linear divergence also for Wilson line operators in LR; it is proportional to $|z|/a$ and arises from the tadpole diagram (d4), with a proportionality coefficient which depends solely on the choice of the gluon action. The exact term that leads to such a divergence is given by:

[TABLE]

where $S_{g}(p)_{\mu\mu}$ is the diagonal matrix element of the gluon propagator in the direction of the Wilson line ($\mu$) and $\overline{p}$ is the loop momentum $p$ with its $\mu$-component set to 0.

Just as with other contributions to the bare Green’s function, the linear divergence is the same – at one-loop level – for all operator insertions. In a resummation of all orders in perturbation theory, the powers of $|z|/a$ are expected to combine into an exponential of the form [10]:

[TABLE]

where $c$ is given by Eq. (43) to one loop, and ${\tilde{\Lambda}}_{\Gamma}$ is related to $\Lambda_{\Gamma}^{\overline{\rm MS}}$ by a further renormalization factor which is at most logarithmically divergent with $a$. Based on arguments from heavy quark effective theory, additional contributions may appear in the exponent [51]. This will be discussed in subsection 3.3.

To one loop, we find the following form for the difference between the bare lattice Green’s functions and the ${\overline{\rm MS}}$-renormalized ones:

[TABLE]

Using Eq. (45) together with Eq. (15) one can extract the multiplicative renormalization and mixing coefficients in the ${\overline{\rm MS}}$-scheme and LR. In our calculation the coefficient $\alpha_{3}$ is negative, in accordance with Eq. (44).

In Eq. (45) all coefficients $\alpha_{i}$ depend on the Symanzik parameters, except for $\alpha_{2}$ which has a numerical value $\alpha_{2}=3.792$. This value was expected, as all gauge dependence must disappear in the ${\overline{\rm MS}}$ scheme for gauge invariant operators: Indeed, this term will cancel against a similar term in $Z^{LR,\,{\overline{\rm MS}}}_{\psi}$ in Eq. (15). The latter has been computed using the same set-up in a previous work [52] and has the general form:

[TABLE]

A few interesting properties of Eq. (45) can be pointed out: The contribution $\left(\Gamma\cdot\gamma_{\mu}{+}\gamma_{\mu}\cdot\Gamma\right)$ indicates mixing between operators of equal dimension, which is finite and appears in the lattice regularization. To one-loop level we find that the mixing matrix is symmetric, and thus, its eigenvalues are the same as those of the addition and difference of the operators which mix which each other. We would also like to note that mixing of different origin has been discussed in Ref. [24] using the parton virtuality distribution functions (VDFs) formalism. Moreover, this combination vanishes for certain choices of the Dirac structure $\Gamma$ in the operator. For the operators $P$, $V_{\nu}$ ($\nu\neq\mu$), $A_{\mu}$, $T_{\mu\nu}$ ($\nu\neq\mu$) the combination $\left(\Gamma\cdot\gamma_{\mu}+\gamma_{\mu}\cdot\Gamma\right)$ is zero and only a multiplicative renormalization is required. This has significant impact in the non-perturbative calculation of the unpolarized quasi-PDFs, as there is a mixing with a twist-3 scalar operator [53]. Such a mixing must be eliminated using a proper renormalization prescription, ideally non-perturbatively [36].

To one-loop level, the diagonal elements of the mixing matrix (multiplicative renormalization) are the same for all operators under study, and through Eq. (15) one obtains:

[TABLE]

where the coefficients $e_{1}-e_{4}$ are given in Table 2, for the Wilson, tree-level Symanzik and Iwasaki improved actions. Results on other gluonic actions can be provided upon request. It is worth mentioning that the presence of $c_{\rm SW}$ in $Z_{\Gamma}^{LR,{\overline{\rm MS}}}$ is inherited from $Z_{\psi}$. As expected, $Z_{\Gamma}^{LR,{\overline{\rm MS}}}$ is gauge independent, and the cancelation of the gauge dependence was numerically confirmed up to ${\cal O}(10^{-5})$. This gives an estimate on the accuracy of the numerical loop integrations. The systematic error stemming from our integration procedure is taken into account by reporting only 6 significant digits. Similar to $Z_{\Gamma}^{LR,{\overline{\rm MS}}}$, the nonvanishing mixing coefficients are operator independent and have the general form:

[TABLE]

where $Z_{ij}^{LR,{\overline{\rm MS}}}$ ($i\neq j$) is nonzero only for the operator pairs: $\{S,\,V_{1}\}$, $\{A_{2},\,T_{34}\}$, $\{A_{3},\,T_{42}\}$, $\{A_{4},\,T_{23}\}$. The values of the coefficients $e_{5}$ and $e_{6}$ for Wilson, tree-level Symanzik and Iwasaki gluons are shown in the last two columns of Table 2. Unlike the diagonal elements, the dependence of $Z_{12}^{LR,{\overline{\rm MS}}}$ on the clover parameter is genuinely extracted from the Green’s functions of the operators. Given that the strength of mixing depends on the value of $c_{\rm SW}$, one may consider increasing $c_{\rm SW}$ from 0 to $-e_{5}/e_{6}$ in order to suppress mixing: for example, choosing: $c_{\rm SW}=0,\,1,\,1.5$ for the Iwasaki action, as shown in Table 2, the mixing coefficient is ${\cal O}(g^{2}{\times}10^{-1}),\,{\cal O}(g^{2}{\times}10^{-2}),\,{\cal O}(g^{2}{\times}10^{-3})$, respectively, and completely vanishes (to one loop) at $c_{\rm SW}=-e_{5}/e_{6}$. This information allows one to tune the clover parameter in order to eliminate mixing at one loop.

3.2.3 Renormalized Green’s functions

We now present the renormalized Green’s functions and discuss some of their important features which will allow us to set up a procedure for estimating the linear divergence using non-perturbative data. This procedure will be explained in the next Subsection.

The complete expressions for the renormalized Green’s functions are listed in Eqs. (72)-(76) of Appendix B for general values of the momentum, length of Wilson line, and gauge fixing parameter. As the renormalized Green’s functions do not depend on the regularization choice, their expressions do not contain $c_{\rm SW}$ and the lattice spacing. For ease of notation, Eqs. (72)-(76) are expressed in terms of the integrals $F_{1}$ - $F_{5}$ and $G_{1}$ - $G_{5}$ of Appendix A. In the present Subsection we focus on the simpler case in which the external momentum, $q$, is in the same direction as the Wilson line, $\mu$ ; we will thus denote $q_{\mu}$ simply by $q$. In this case the renormalized Green’s functions are multiples of their tree-level values. In particular, they can be cast into the form:

[TABLE]

where $F_{\Gamma}$ are complicated complex functions of $qz$, satisfying $F_{\Gamma}(-x)=F_{\Gamma}^{\dagger}(x)$. In standard fashion, we have exponentiated the one-loop result for the $\bar{\mu}$-dependence, thus putting into evidence the anomalous dimension of the Green’s functions. Once again we note that this scale dependence does not vary with the length of the Wilson line; this is in agreement with older calculations, regarding closed Wilson loops [10].

The form of Eq. (49) is in accordance with the features of the physical matrix elements of Wilson line operators computed in numerical simulations. The requirement that renormalization functions be $q$-independent implies that the dependence on $aq$ in the bare Green’s functions can be predicted from the start, as it has to match the $\bar{\mu}$ dependence of the renormalized Green’s functions. Below we provide the expressions for $F_{\Gamma}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In Appendix B we present the ${\overline{\rm MS}}$-renormalized Green’s functions for each Wilson line operator, and for general RI*′* renormalization scale 4-vector $\bar{q}_{\nu}$ . For demonstration purposes, in Fig. 3 we plot $F_{V_{1}(A_{1})}$ in the Landau and Feynman gauge as a function of the dimensionless quantity $qz$. We find that the qualitative behavior of $F_{V_{1}(A_{1})}$ is similar regardless of the choice for $\beta$. It is interesting to see the limit $q\to 0$ for $F_{V_{1}(A_{1})}(qz)$, which, in fact, coincides with the limit $z\to 0$. As mentioned earlier, the latter is expected to be singular, due to the emergence of contact terms. Thus, as $q\to 0$, we find that $F_{V_{1}(A_{1})}(qz)$ is a real function independent of the gauge and equal to:

[TABLE]

3.3 Non-perturbative elimination of the linear divergence

As shown in Subsection 3.2.2, there is a linear divergence in the lattice-regularized Wilson line operator, which requires a careful removal before the continuum limit can be reached in the non-perturbative matrix elements. One way to eliminate this divergence is to use the estimate of the one-loop coefficient $e_{2}$ of Eq. (47) and subtract it from the non-perturbative matrix elements. We have computed $e_{2}$ for a large class of gluon actions555We recall that the value of $e_{2}$ does not depend on the fermion action., and presented in this paper the results for Wilson, TL Symanzik and Iwasaki actions. However, this subtraction only partially removes the divergence, as higher orders still remain and they will dominate in the $a\to 0$ limit. In general, power divergences are features in which deviations from perturbation theory are most pronounced. Thus, it is preferable to develop a non-perturbative method to extract the linear divergence.

A method suggested recently [27] is to use the static potential in order to eliminate non-perturbatively the exponential with the linear divergence, cf. Eq. (44). Here we propose an alternative methodology which is based on using bare matrix elements of the Wilson line operators from numerical simulations, denoted by $q(P_{3},z)$ in Ref. [5]. We consider the case where there is no mixing between operators and we focus on the helicity (axial) and transversity (tensor), indicated by $\Delta h(P_{3},z)$ and $\delta h(P_{3},z)$, respectively, in Ref. [5]. We note that both helicity and transversity have a Dirac index in the direction of the Wilson line, which exhibits no mixing. In the non-perturbative calculations of the matrix elements, the nucleon is boosted by momentum $P_{3}$ which is in the same direction as the Wilson line, simplifying the calculation. Based on the arguments presented in the previous Subsection (Eq. (49)) we expect that the renormalized matrix elements can depend on $z$ only through the dimensionless quantity $P_{3}z$. Furthermore, the dependence on the scale $\bar{\mu}$ in the renormalized matrix element is well defined and involves the anomalous dimension ($\gamma_{\Gamma}$) of the operator: $q^{R}(P_{3}z,P_{3}/\bar{\mu})\propto{\bar{\mu}}^{-2\gamma_{\Gamma}}$, which is matched by the $\bar{\mu}$ dependence in the renormalization function. Thus,

[TABLE]

Similarly, the renormalization function $Z_{\Gamma}^{LR,{\overline{\rm MS}}}(a\bar{\mu},z/a)$, given its expected $\bar{\mu}$-dependence, will factorize as:

[TABLE]

The factor ${\tilde{Z}}_{\Gamma}(a\bar{\mu})$ is not simply $(a\bar{\mu})^{2\gamma_{\Gamma}}$, due to finite lattice corrections; in the one-loop case, these lattice corrections are the terms containing $e_{1}$ , $e_{3}$ , $e_{4}$ in Eq. (47). The factor $\hat{Z}(z/a)$ originates exclusively from tadpole diagrams, such as d4 of the one-loop case. Based on the behavior of tadpole diagrams, one can see that the one-loop contribution proportional to $e_{2}$ in Eq. (47) will exponentiate upon considering higher powers of $g$, leading to:

[TABLE]

This behavior is entirely consistent with the exponential behavior $\exp(-\delta m\,|z|/a)$ proven in Ref. [10] for closed Wilson loops. We note that the proof in Ref. [10] holds for any regularization in which $z/a$ terms may arise.

We can thus write the ratio of the bare matrix elements for different values of $P_{3}$ and $z$ as:

[TABLE]

where the one-loop anomalous dimension is $\gamma_{\Gamma}=-3g^{2}C_{f}/(16\pi^{2})$ for all operator insertions. The anomalous dimension of the fermion field is not relevant in this discussion as the non-perturbative matrix elements are between physical (nucleon) states. In the above ratio, one may choose $P_{3},\,P^{\prime}_{3},\,z,\,z^{\prime}$ such that $P_{3}z=P^{\prime}_{3}z^{\prime}$, which simplifies the ratio considerably:

[TABLE]

Thus, by forming the ratio $q(P_{3},z)/q(P^{\prime}_{3},z^{\prime})$ from non-perturbative data, and by choosing several combinations of $P_{3}z=P^{\prime}_{3}z^{\prime}$, one can fit to extract the coefficient of the linear divergence, $\delta m$.

We note that there are non-perturbative arguments [51] to suggest that a further finite, dimensionful scale may appear multiplying $z$ in the exponential term: $\exp(-\delta m\,|z|/a-c|z|)$. In this case, Eq. (60) takes the form:

[TABLE]

Thus, we may still deduce the value of the quantity $(\delta m/a)+c$. Extracting the values of $\delta m$ and $c$ separately would require utilizing simulation data from two or more values of the lattice spacing $a$; however, this procedure will be hampered by the very dependence of $\delta m$ on the coupling constant, as the lattice spacing is varied. We will not pursue further the non-perturbative evaluation of $\delta m$ and $c$, since it is beyond the scope of this paper.

The ratio of the left-hand-side of Eq. (61) can be investigated for the helicity and transversity, which do not exhibit mixing. Since the right-hand side of Eq. (61) is independent of the operator insertion, one expects the same value for the exponential coefficient, up to lattice artifacts. We have tested this method with the data of ETMC presented in Ref. [5], with encouraging results:

• •

The ratio $q(P_{3},z)/q(P^{\prime}_{3},z^{\prime})$ was found to be real for the helicity and transversity as expected from Eq. (61), despite the fact that the matrix elements themselves are complex;

• •

The analogous ratio for the unpolarized operator, which mixes with the scalar, leads to a nonzero imaginary part;

• •

The extracted value for the coefficient $(\delta m/a)+c$, using different combinations of $P_{3}z$, is consistent within statistical accuracy;

• •

Both helicity and transversity give very similar estimates for $(\delta m/a)+c$.

Indeed, Fig. 4 shows that the unpolarized case from Ref. [5] (blue circles) has a nonzero imaginary part, while the imaginary parts of the ratios for helicity and for transversity on the same ensemble are compatible with zero. The simulation data of Ref. [5] regard twisted mass fermions at a pion mass of $m_{\pi}{=}375$ MeV, and Iwasaki gluons. For these action parameters and $\beta{=}1.95$, our perturbative results of Eq. (48) show a significant mixing, that is $(6/\beta)\cdot C_{f}/(16\pi^{2})\cdot\left(9.93653\right){\sim}0.26$. This is confirmed by the data of Fig. 4. It is also very interesting to test the left-hand-side of Eq. (61) on an ensemble in which the mixing is expected to be very small. ETMC has preliminary data at the physical point for twisted mass fermions ($\beta{=}2.1$) including a clover term $c_{\rm SW}{\sim}1.57$, and Iwasaki gluons. Thus, according to the results in Eq. (48) and the numerical values presented in Table 2, the one-loop mixing is reduced by two orders of magnitude as $(6/\beta)\cdot C_{f}/(16\pi^{2})\cdot\left(9.93653-6.52764\,c_{\rm SW}\right){\sim}\,{-}0.00752345$. The fact that the mixing is insignificant in the presence of a clover term with $c_{\rm SW}{\sim}1.57$, is also confirmed by the simulation data (orange filled diamonds), as shown in Fig. 4. We note that even though the statistical errors are still large, the imaginary part of Eq. (61) is zero within error bars.

The above demonstration using simulation data stresses the importance of eliminating the mixing between the vector and scalar operators, prior to applying the matching procedure LaMET (Large-Momentum Effective Field Theory) to the physical PDFs. This can be achieved by computing the matrix elements for both the scalar and vector operators. One may use our perturbative results for the multiplicative renormalization and the mixing coefficients to disentangle the two operators and extract the renormalized unpolarized quasi-PDFs. Once this is done, one can apply the ratio procedure of Eqs. (59)-(60) to these operators as well.

4 Summary - Future Work

In this paper we have presented the one-loop perturbative calculation of the renormalization functions for operators including a straight Wilson line. Results of this work may be used to renormalize, to one-loop in perturbation theory, matrix elements for quasi-PDFs. Two main features of the quasi-PDFs in lattice regularization are demonstrated in this work: finite mixing and a linear divergence with respect to the regulator. We computed the one-loop Green’s functions both in dimensional (DR) and lattice (LR) regularizations, which allows one to extract the LR renormalization functions in the ${\overline{\rm MS}}$-scheme directly, without an intermediate RI*′* scheme. Using the Green’s functions in DR we computed the conversion factor between an RI*′*-type (Eq. (20)) scheme and ${\overline{\rm MS}}$. These expressions are necessary to bring non-perturbative estimates of the renormalization functions to the ${\overline{\rm MS}}$-scheme.

We demonstrate for the first time that certain Wilson line operators exhibit mixing within the lattice regularization. This has great impact on the unpolarized quasi-PDFs that mix with a twist-3 [53] scalar operator. Thus, before matching to the physical PDFs one must eliminate the mixing, which would require computation of the matrix elements for the scalar operator. However, our results indicate that the presence of a clover term in the fermion action highly suppresses the mixing, as can be read from Eq. (48). This has been tested on nucleon matrix elements of the unpolarized quasi-PDFs using twisted mass fermions, with and without clover improvement (Fig. 4), confirming what is expected from our perturbative results. This may also be an indication that higher loop contributions are suppressed. Furthermore, the mixing discussed for the unpolarized case vanishes if one uses a Dirac structure perpendicular to the Wilson line direction, as explained in the main text. For example, choosing the $\gamma$-matrix in the temporal direction is important for a faster convergence to the physical PDFs, as discussed in Ref. [24].

A feature of the quasi-PDFs that requires special treatment is the linear divergence that complicates taking the continuum limit. We have computed the one-loop coefficient of the linear divergence for a variety of gluonic actions, which can be subtracted from nucleon matrix elements of the quasi-PDFs. Such a divergence is expected to resum to an exponential, just as in the case of Wilson loop operators [10]. Using the form of the renormalized Green’s function of Eq. (49) we propose a technique to extract the coefficient of the linear divergence in Subsection 3.3 from non-perturbative data.

Using the renormalization pattern exposed in this work we have developed an appropriate non-perturbative renormalization prescription for the unpolarized, helicity and transversity quasi-PDFs. Such a scheme will extract both the renormalization functions and linear divergence at once. For the unpolarized case, the mixing with the twist-3 scalar operator is also addressed. This will be presented in a follow-up publication [36].

A natural continuation of this project is the addition of smearing to the fermionic part of the action and/or to the gauge links of the Wilson line operator. This is important, as modern simulations employ such smearing techniques (e.g., stout and Hypercubic (HYP)) that suppress the power divergence and bring the renormalization functions closer to their tree-level values. Smearing the operator under study alters its renormalization functions, thus, the same smearing must be employed in the renormalization process. In these cases, the additional contributions to the renormalization functions due to smearing are more convergent, and thus the perturbative extraction of singularities is simpler. Nevertheless, the smearing of the gauge links results in a huge increase of the number of terms in the vertices and, thus, the computation of the diagrams becomes very challenging.

An extension of this calculation that we intend to pursue, is the evaluation of lattice artifacts to one loop and to all orders in the lattice spacing, ${\cal O}(a^{\infty},\,g^{2})$, as developed in Ref. [54, 55]. This has been successfully applied for local and one-derivative fermion operators eliminating lattice artifacts from non-perturbative estimates [56, 57]. For operators with a long Wilson line ($z{>>}a$) the lattice artifacts are likely to be more prominent, and therefore, such a calculation will be extremely beneficial in the non-perturbative renormalization program presented in Ref. [36]. It will also have implications in the comparison of the quasi-PDFs and phenomenological estimates for the physical PDFs.

A possible addition to the present work is the two-loop calculation in dimensional regularization, from which one can extract the conversion factor between different renormalization schemes, as well as the anomalous dimension of the operators. The conversion factor up to two loops may be applied to non-perturbative data on the renormalization functions, to bring them to the ${\overline{\rm MS}}$-scheme at a better accuracy. Furthermore, knowledge of the two-loop expression for the anomalous dimension in Eq. (61) will improve the method for extracting the linear divergence, and will eliminate systematic uncertainties related to the truncation of the conversion factor.

Another direction that one may follow is the calculation of the one-loop Green’s functions at a nonzero fermion mass to one loop level. We expect the difference between the finite-mass and massless cases to be small given the smooth behavior of the conversion factor, but it would be interesting to investigate, as simulations are not exactly at zero renormalized mass. Differences in the renormalization functions of flavor-singlet and -nonsinglet Wilson line operators will already show up at one loop in the massive case, and at two loops in the massless case.

Finally, the techniques developed in this work for the renormalization of quasi-PDFs may be inspiring for the renormalization of Wilson-line fermion operators of different structure, such as staples. Eventually, understanding the renormalization pattern of such operators may lead to the development of a non-perturbative prescription. This will be of high importance for matrix elements of the transverse momentum-dependent parton distributions (TMDs) that are currently under investigation for the nucleon and pion in lattice QCD [23]. Similar to the case of the unpolarized quasi-PDFs, there might be an indication for presence of mixing in certain cases [58].

Acknowledgements

We would like to thank the members of ETMC for useful and fruitful discussions. In particular, we are thankful to Krzysztof Cichy for discussions in the preparation of this manuscript. MC acknowledges financial support by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the TMD Topical Collaboration, as well as, by the National Science Foundation under Grant No. PHY-1714407.

APPENDICES

Appendix A Basis of Integrals

The Green’s functions for the operators including a Wilson line depend on the external momentum four-vector, $q_{\nu}$, and on the length of the Wilson line, $z$ , in a complicated way which is not amenable to a closed analytic form. In order to present our general results in a compact manner we introduce two classes of integrals:

• 1.

$F_{1}$ - $F_{5}$: integrals over the Feynman parameter $x$

• 2.

$G_{1}$ - $G_{5}$: integrals over $x$ and over the parameter $\zeta$ of Eq. (8)

The integrands contain modified Bessel functions of the second kind, $K_{0}$ and $K_{1}$. All integrals presented here are convergent and can be computed numerically.

[TABLE]

[TABLE]

Here $q\equiv\sqrt{q^{2}}$. Note that the $F$ integrals are dimensionless, $G_{1}$ - $G_{2}$ have dimensions of $z$ due to integration over $\zeta$, and $G_{3}$ - $G_{5}$ have dimensions of $z^{2}$. The dimensionality of all integrals is combined with the momentum in such a way that the conversion factors and Green’s functions are dimensionless (see, e.g., Eqs. (72)-(76)).

Appendix B Renormalized Green’s Functions

In this Appendix we present the general expressions for the renormalized Green’s functions of the operators given in Eqs. (10), in the ${\overline{\rm MS}}$ scheme. The functions $\Lambda^{\rm 1-loop}_{\Gamma}$ are complex, and have a complicated dependence on the momentum and the length of the Wilson line; thus, we write them in a compact form, using the list of integrals $F_{i}(q,z)$ and $G_{i}(q,z)$ defined in Appendix A.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]