Renormalization of Wilson-line operators in the presence of nonzero quark masses

In this paper, we examine the effect of nonzero quark masses on the renormalization of gauge-invariant nonlocal quark bilinear operators, including a finite-length Wilson line (called Wilson-line operators). These operators are relevant to the definition of parton quasi-distribution functions, the calculation on the lattice of which allows the direct nonperturbative study of the corresponding physical parton distribution functions. We present our perturbative calculations of the bare Green's functions, the renormalization factors in RI' and MSbar schemes, as well as the conversion factors of these operators between the two renormalization schemes. Our computations have been performed in dimensional regularization at one-loop level, using massive quarks. The conversion factors can be used to convert the corresponding lattice nonperturbative results to the MSbar scheme, which is the most widely used renormalization scheme for the analysis of experimental data in high-energy physics. Also, our study is relevant for disentangling the additional operator mixing which occurs in the presence of nonzero quark masses, both on the lattice and in dimensional regularization.

To date, all lattice studies of the renormalization of Wilson-line operators have only considered massless fermions, expecting that the presence of quark masses can cause only imperceptible changes; this is indeed a reasonable assumption for light quarks. However, for heavy quarks this statement does not hold. In addition, simulations cannot be performed exactly at zero renormalized mass. One could, of course, adopt a zero-mass renormalization scheme even for heavy quarks, but such a scheme is less direct and entails additional complications. Thus, it would be useful to investigate the significance of finite quark masses on the renormalization of Wilson-line operators. This is the goal of our present study.
In this work, we calculate the conversion factors from RI ′ to the MS scheme, in dimensional regularization (DR) at one-loop level for massive quarks. The conversion factors can be combined with the regularization independent (RI ′ )-renormalization factors of the operators, computed in lattice simulations, in order to calculate the nonperturbative renormalization of these operators in MS. Nonperturbative evaluations of the renormalization factors cannot be obtained directly in the MS scheme, since the definition of MS is perturbative in nature; most naturally, one calculates them in a RI ′ -like scheme, and then introduces the corresponding conversion factors between RI ′ and the MS scheme, which rely necessarily on perturbation theory. Given that the conversion between the two renormalization schemes does not depend on regularization, it is more convenient to evaluate it in DR. Thus, the perturbative calculation of conversion factors is an essential ingredient for a complete study of quasi-PDFs. This work is a continuation to a previous paper [16], in which, among other results, one-loop conversion factors of Wilson-line operators are presented for the case of massless quarks.
In studying composite operators, one issue which must be carefully addressed is that of possible mixing with other similar operators. Many possibilities are potentially present for the nonlocal operators which we study: (A) Operators involving alternative paths for the Wilson line joining the quark pair will not mix among themselves, as demonstrated in Ref. [52] (and also in Refs. [44,45] for the case of closed Wilson loops). This property is related to translational invariance and is similar to the lack of mixing between a local composite operator O(x) with O(y). Given that a discrete version of translational invariance is preserved on the lattice, nonlocal operators involving different paths should not mix also on the lattice.
(B) Operators involving only gluons will also not mix. This can be seen, e.g., via the auxiliary field approach (e.g., Ref. [52]); as a specific case, the operator of Eq. (1) cannot mix with an operator containing the gluon field strength tensor in lieu of the quark fields (joined by a Wilson line in the adjoint representation), since this operator is higher dimensional.
(C) There may also be mixing among operators with different flavor content in a RI ′ scheme, depending on the scheme's precise definition. However, the mixing is expected to be at most finite and thus not present in the MS scheme; by comparing to the massless case, in which exact flavor symmetry allows no such mixing, the difference between the massive and massless case will bear no superficial divergences, since the latter are UV regulated by the masses. The auxiliary field approach, by involving only composite operators in the (anti-)fundamental representation of the flavor group, shows that no flavor mixing needs to be introduced.
Even in the absence of quark masses, bare Green's functions of Wilson-line operators may contain finite, regulator-dependent contributions which cannot be removed by a simple multiplicative renormalization; as a consequence, an appropriate (i.e., regularizationindependent) choice of renormalization prescription for RI' necessitates the introduction of mixing matrices for certain pairs of operators [16], both in the continuum and on the lattice. The results of the present work demonstrate that the presence of quark masses affects the observed operator-mixing pairs, due to the chiral-symmetry breaking of mass terms in the fermion action. Compared to the massless case on the lattice [16], the mixing pairs remain the same for operators with equal masses of external quark fields, i.e., (11, γ 1 ), (γ 5 γ 2 , γ 3 γ 4 ), (γ 5 γ 3 , γ 4 γ 2 ), and (γ 5 γ 4 , γ 2 γ 3 ), where by convention 1 is the direction of the straight Wilson line and 2, 3, and 4 are directions perpendicular to it. However, for operators with different masses of external quark fields, flavor-symmetry breaking leads to additional mixing pairs: (γ 5 , γ 5 γ 1 ), (γ 2 , γ 1 γ 2 ), (γ 3 , γ 1 γ 3 ), and (γ 4 , γ 1 γ 4 ). As a consequence, the conversion factors are generally nondiagonal 2 × 2 matrices. This is relevant for disentangling the observed operator mixing on the lattice. Also, comparing the massive and the massless cases, the effect of finite mass on the renormalization of Wilson-line operators becomes significant for strange quarks, as well as for heavier quarks. These are features of massive quasi-PDFs, which must be taken into account in their future nonperturbative study.

3
The outline of this paper is as follows. In Sec. II we provide the theoretical setup related to the definition of the operators which we study, along with the necessary prescription of the renormalization schemes. Sec. III contains our results for the bare Green's functions in DR, the renormalization factors, as well as the conversion factors of these operators between the renormalization schemes. In Sec. IV, we present several graphs of the conversion factor matrix elements for certain values of free parameters. Finally, in Sec. V, we conclude with possible future extensions of our work.
We have also included two Appendices. Appendix A contains a discussion on technical aspects, such as the methods that we used to calculate the momentum-loop integrals, as well as the limits of vanishing regulator and/or masses. A table of Feynman parameter integrals, which appear in the expressions of our results, is relegated to Appendix B.

B. Definition of renormalization schemes
Taking into account the presence of nonzero fermion masses in our calculations, we adopt mass-dependent prescriptions for the renormalization of Wilson-line operators. We define the renormalization factors which relate the bare O Γ with the renormalized operators O R [In the presence of operator mixing, this relationship is appropriately generalized; see Eq. (8)]. The corresponding renormalized one-particle irreducible (1-PI) amputated Green's functions of Wilson-line operators Λ R amp are given by where Λ Γ = ψ O Γψ amp are the bare amputated Green's functions of the operators and Z ψ f is the renormalization factor of the fermion field with flavor f , defined by ψ R In the massive case, renormalization factors of the fermion and antifermion fields appearing in bilinear operators of different flavor content may differ among themselves, as the fields have generally different masses.

Renormalization conditions for fermion fields and masses
At this point, we provide the conditions for the mass-dependent renormalization of fermion fields, as well as the multiplicative renormalization of fermion masses: are the bare (renormalized) masses for each flavor]; the latter is not involved in our calculations, but we include it for completeness.
In MS, renormalization factors Z ψ of the fermion field and Z m of the fermion mass must contain, beyond tree level, only negative powers of ε (the regulator in DR in D dimensions, D ≡ 4 − 2ε); their values are fixed by the requirement that the renormalized fermion selfenergy be a finite function of the renormalized parameters m MS and g MS (g = µ ε Z g g MS ; µ is a dimensionful scale): In RI ′ , convenient conditions for the fermion field of a given flavor and the corresponding mass are whereq ν is the RI ′ renormalization scale 4-vector, m RI ′ is the RI ′ -renormalized fermion mass, N c is the number of colors, and the symbol X can be any regularization, such as DR or lattice. These conditions are appropriate for lattice regularizations which do not break chiral symmetry, so the Lagrangian mass m 0 coincides with the bare mass m B , e.g., staggered/overlap/domain wall fermions. For regularizations which break chiral symmetry, such as Wilson/clover fermions, a critical mass m c is induced; one must first find the value of m c by a calibration in which one requires that the renormalized mass for a "benchmark" meson attains a desired value, e.g., zero pion mass, and then set m B = m 0 − m c .

Renormalization conditions for Wilson-line operators
As is standard practice, we will derive the factors Z Γ by imposing appropriate normalization conditions on the quark-antiquark Green's functions of O Γ .
In the spirit of MS, Z DR,MS Γ contains, beyond tree level, only negative powers of ε. Here, we note that the leading poles in n-loop diagrams of bare Green's functions, O(1/ε n ) (n ∈ Z + ), are multiples of the corresponding tree-level values and thus do not lead to any mixing. Subleading poles will not lead to divergent mixing coefficients, as is implicit in the renormalizability proofs of Refs. [44,45,52]. So, in the MS scheme, we can use the standard definitions of renormalization factors, as in Eq. (2).
In RI ′ , things are more complicated. There is, a priori, wide flexibility in defining RI ′like normalization conditions for Green's functions. Given that no mixing is encountered in MS renormalization and given that any other scheme can only differ from MS by finite factors, one might a priori expect to be able to adopt a deceptively simple prescription, in which RI ′ -renormalized operators are simply multiples of their bare counterparts, satisfying a standard normalization condition: where Λ tree Γ = Γ exp(iq µ z) is the tree-level value of the Green's function of operator O Γ and Λ RI ′ Γ is defined through Eqs. (3) and (2). There is, however, a fundamental problem with such a prescription: the renormalized Green's function resulting from Eq. (7) will depend on the regulator which was used in order to compute it (and, thus, it will not be regularization independent, as the name RI suggests). As was pointed out in Ref. [16], bare Green's functions of O Γ , computed on the lattice, contain additional contributions proportional to the tree-level Green's function of O Γ ′ , where Γ ′ = Γγ µ +γ µ Γ (whenever the latter differs from zero). Such contributions will not be eliminated by applying the renormalization prescription of Eq. (7), thus leading to renormalized Green's functions which differ from those obtained in DR. It should be pointed out that, in all cases, the renormalized functions will contain a number of tensorial structures, the elimination of which may be possible at best only at a given value of the renormalization scale. However, the main concern here is not the elimination of mixing contributions, desirable as this might be; what is more important is to establish a RI ′ scheme which is indeed regularization independent, so that nonperturbative estimates of renormalization factors can be converted to the MS scheme using conversion factors which are regulator independent.
Given the preferred direction µ of the Wilson-line operator, there is a residual rotational (or hypercubic, on the lattice) symmetry with respect to the three remaining transverse directions, including also reflections. As a consequence, given an appropriate choice of a renormalization scheme, no mixing needs to occur among operators which do not transform in the same way under this residual symmetry. In particular, mixing can occur only among pairs of operators (O Γ , O Γγµ ).
Denoting generically the two operators in such a pair by (O Γ 1 , O Γ 2 ), the corresponding renormalization factors will be 2 × 2 mixing matrices: 6 More precisely, the mixing pairs , and (γ 5 γ 4 , γ 2 γ 3 ). Therefore, the renormalized 1-PI amputated Green's functions of Wilson-line operators have the following form: Thus, an appropriate renormalization condition, especially for lattice simulations, is Combining Eqs. (9) and (10), the RI ′ condition takes the form: Based on the above symmetry arguments, such a RI ′ condition will indeed be regularization independent, for all regularizations which respect the above symmetries.
One could of course adopt more general definitions of RI ′ , e.g., a prescription in which each of the 16 operators O Γ can contain admixtures of some of the remaining operators: in such a way that the renormalized Green's functions will satisfy a condition similar to Eq. (10), but with the indices i, j ranging from 1 to 16. However, such a definition would introduce additional finite mixing, which would violate the rotational symmetry in the transverse directions, e.g., mixing among O γ 1 and O γ 2 ; such a violation would occur whenever the RI ′ renormalization scale 4-vectorq is chosen to lie in an oblique direction. To avoid such unnecessary mixing, it is thus natural to adopt the "minimal" prescription of Eqs. (8) - (11). Since this prescription extends beyond one-loop order, it may be applied to nonperturbative evaluations of the renormalization matrices Z L,RI ′ .
Let us note that, as it stands, Eq. (10) leads to renormalization factors which depend on the individual components ofq, rather than justq 2 andq µ ; consequently, the renormalization factors of, e.g., O γ 2 and O γ 3 will have different numerical values. One could, of course, define RI ′ in such a way that the residual invariance is manifest; this can be seen by analogy with local operators, e.g., 3,4), and, in doing so, Z V turns out to depend only on the length of the renormalization scale 4-vector. Adopting such a definition, the values of the conversion factors can be read off our bare Green's functions [see Eqs. (15) - (26) below] in a rather straightforward way, and they will indeed depend only onq 2 andq µ . However, in defining the RI ′ scheme for Wilson-line operators, we have aimed at being as general as possible and thus did not take any averages, as above, in order to accommodate possible definitions employed in nonperturbative investigations of the renormalization factors; after all, the conversion factors which we calculate must be applicable precisely to these investigations.
It goes without saying that if one chooses all components of the renormalization scale 4vector, perpendicular to the Wilson line, to vanish, then residual rotational invariance is automatically restored.
Finally, one could define RI ′ in such a way that renormalization factors would be strictly real, e.g., by taking the absolute value of the lhs in Eq. (10); indeed, the choice of the definition of RI ′ , leading to complex renormalization factors, is not mandatory, but it is a natural one, following the definition used in nonperturbative investigations. All these choices are related to the MS scheme via finite conversion factors; thus, no particular choice is dictated by the need to remove divergences, either in dimensional regularization or on the lattice.

C. Conversion factors
As a consequence of the 2×2 matrix form of the RI ′ renormalization factors, the conversion factors between RI ′ and MS schemes will also be 2 × 2 mixing matrices. Being regularization independent, they can be evaluated more easily in DR. They are defined as We note in passing that the definition of the MS scheme depends on the prescription used for extending γ 5 to D dimensions 3 ; this, in particular, will affect conversion factors for the pseudoscalar and axial-vector operators. However, such a dependence will only appear beyond one loop. Now, the Green's functions in the RI ′ scheme can be directly converted to the MS scheme through where is the conversion factor for a fermion field of a given flavor.

III. COMPUTATION AND RESULTS
In this section, we present our one-loop results for the bare Green's functions of Wilsonline operators, the renormalization factors, and the conversion factors between RI ′ and MS schemes, using dimensional regularization. In this regularization, Green's functions are Laurent series in ε, where ε is the regulator, defined by D ≡ 4 − 2ε, and D is the number of Euclidean spacetime dimensions, in which momentum-loop integrals are well defined. We also investigate the operator mixing. the quark masses, and therefore its contribution is the same as that of the massless case. In Appendix A, we describe the method that we used to calculate the momentum-loop integrals of the above diagrams. Below, we provide our results for the bare Green's function of operators for each Feynman diagram separately. Our expressions depend on integrals of modified Bessel functions of the second kind, K n , over Feynman parameters. These integrals are presented in Eqs. (B1) -(B16) of Appendix B. For the sake of brevity, we use the following notation: f ij ≡ f i (q, z, m j ), g ij ≡ g i (q, z, m j ), and h i ≡ h i (q, z, m 1 , m 2 ). Also, index µ is the direction parallel to the Wilson line; indices ν, ρ, and σ are the directions perpendicular to the Wilson line; and µ, ν, ρ, and σ are all different among themselves. Furthermore,μ is the MS renormalization scale,μ ≡ µ (4π/e γ E ) 1/2 , where µ (not to be confused with the spacetime index µ) appears in the renormalization of the D-dimensional coupling constant; g = µ ε Z g g R , and γ E is the Euler constant. In addition, is the Casimir operator, and β is the gauge fixing parameter, defined such that β = 0(1) corresponds to the Feynman (Landau) gauge. Finally, symbols S (scalar), P (pseudoscalar), V µ (vector in the µ direction), V ν (vector in the ν direction), A µ (axial-vector in the µ direction), A ν (axial-vector in the ν direction), T µν (tensor in the µ, and ν directions), and T νρ (tensor in the ν, and ρ directions) correspond to the operators O Γ with Γ = 1 1, γ 5 , γ µ , γ ν , γ 5 γ µ , γ 5 γ ν , γ µ γ ν , γ ν γ ρ , respectively. We note that only tree-level values for the quark masses appear in the following one-loop expressions: UV-divergent terms of order O(1/ε) arise from the last three diagrams. These terms are multiples of the tree-level values of Green's functions and therefore do not lead to any mixing. However, there are finite terms for each O Γ with different Dirac structures than the original operator; some of these terms are responsible for the finite mixing which occurs in RI ′ . In particular, they lead to the expected mixing within the pairs (Γ, Γγ µ ) or equivalently (Γ, γ µ Γ). This is a consequence of the violation of chiral symmetry by the mass term in the fermion action, as well as the flavor-symmetry breaking when masses have different values. For the case of equal masses (no flavor-symmetry breaking) m 1 = m 2 , the mixing pattern reduces to (Γ, 1 2 {Γ, γ µ }), which is the same as the pattern for massless quarks on the lattice. Our findings are expected to be valid also on the lattice.
The one-loop Green's functions exhibit a nontrivial dependence on dimensionless quantities involving the Wilson-line length z, the external quark momentum q, and the quark masses m i (i = 1, 2): zq µ , zm i . This dependence is in addition to the standard logarithmic dependence onμ: log(μ 2 /q 2 ). Also, we note that our results are not analytic functions of z near z = 0; this was expected due to the appearance of contact terms beyond tree level. For the case z = 0, the nonlocal operators are replaced by local massive fermion bilinear operators; their renormalization is addressed in Ref. [76], using a generalization of the RI-SMOM scheme, called RI-mSMOM. Further, the Green's functions of Feynman diagrams satisfy the following reflection relations, with respect to z: [Note that (1/4) tr(Γ 2 ) = ±1, depending on Γ.] The total one-loop bare Green's functions of operators O Γ are given by the sum over the contributions of the four diagrams: B. Renormalization factors

Renormalization factors of fermion field and mass
The perturbative determination of Z ψ and Z m proceeds in textbook fashion by calculating the bare fermion self-energy in DR to one loop; we present it here for completeness. The Feynman diagram contributing to this two-point Green's function is shown in Fig. 2. Denoting by Σ the higher-order terms O(g 2 ) of the 1-PI amputated Green's function of the fermion field, the inverse full fermion propagator takes the following form: Writing Σ in the more useful form: Σ = i / q Σ 1 (q 2 , m) + m 1 1 Σ 2 (q 2 , m), we present the one-loop results for the functions Σ 1 , Σ 2 : The renormalization conditions for Z ψ and Z m in the RI ′ scheme, using the above notation, take the following perturbative forms: Thus, in the presence of finite fermion masses, the results for the renormalization factors of the fermion field and mass are given below: We recall that the mass appearing in the above expressions is the renormalized mass, which coincides with the bare mass to this order. The results for Z ψ and Z m are in agreement with Ref. [77], in the massless limit and forq =μ.
The renormalization factors in the MS scheme can be readily inferred from Eqs. (35) and (36) by taking only the pole part in epsilon:

Renormalization factors of Wilson-line operators
Now, we have all the ingredients for the extraction of renormalization factors of Wilsonline operators in the RI ′ and MS schemes. By writing Z ψ f and Λ Γ in the form: where 4 the condition for the renormalization of Wilson-line operators in the RI ′ scheme, up to one loop, reads The equivalent expression for Z DR,MS Γ follows from Eq. (42), by keeping in λ ij only pole parts in epsilon; the latter appear only for i = j, leading to Our final results are presented below. In the MS scheme, the renormalization factors of operators have the form: in agreement with Refs. [52,78,79]. As we observe, they are independent of operator Γ, fermion masses, Wilson-line length z, and gauge parameter β. In RI ′ , the renormalization factors are given with respect to the conversion factors, which are presented in the next section: The above relation stems from the one-loop expression of Eq. (13).
Our results are in agreement with Ref. [16] in the massless limit 5 . A consequence of the above relations is that, in the case of equal quark masses m 1 = m 2 , the nondiagonal matrix elements of C P,Aµ and C Vν ,Tµν vanish. Also, the matrix elements of conversion factors satisfy the following reflection relation with respect to z: This means that the real part of diagonal (nondiagonal) matrix elements is an even (odd) function of z, while the imaginary part is odd (even).

IV. GRAPHS
In this section, we illustrate our results for conversion factors by selecting certain values of the free parameters used in simulations. To this end, we plot the real and imaginary parts of the conversion factor matrix elements as a function of Wilson-line length, z. For input, we employ certain parameter values, used by ETMC in the ensemble of dynamical N f = 2 + 1 + 1 twisted mass fermions of Ref. [12]; i.e., we set 6 g 2 = 3.077, β = 1 (Landau gauge), N c = 3,μ = 2 GeV, andq = 2π 32a (n z , 0, 0, nt 2 + 1 4 ), for a = 0.082 fm (lattice spacing), n z = 4, and n t = 8 (the Wilson line is taken to lie in the z direction, which, by convention, is denoted by µ = 1). Expressed in GeV,q = (1.887, 0, 0, 2.048) GeV. To examine the impact of finite quark masses on the conversion factors, we plot six different cases of external quark masses: 1. massless quarks (m 1 = m 2 = 0) 2. m 1 = m 2 = 13.2134 MeV, corresponding to the bare twisted mass used in Ref. [12] 3. one up and one strange quark (m 1 = 2.3 MeV, m 2 = 95 MeV) 4. two strange quarks (m 1 = m 2 = 95 MeV) 5. one up and one charm quark (m 1 = 2.3 MeV, m 2 = 1275 MeV) 6. two charm quarks (m 1 = m 2 = 1275 MeV). As regards theq dependence, we have not included further graphs for the sake of conciseness; however, using a variety of values for the components ofq, we find no significant difference. More quantitative assessments can be directly obtained from our algebraic results.
In Figs. 3 and 4, we present graphs of some representative conversion factors (C S,V 1 , C P,A 1 ) for the six cases of external quark masses. The plots are given only for positive values of z, since the behavior of conversion factors for negative values follows the reflection relation of Eq. (56). We observe that the real part of the conversion factor matrix elements is an order of magnitude larger than the imaginary part and that the diagonal elements are an order of magnitude larger than the nondiagonal elements. Also, for increasing values of z, the real part of diagonal elements tends to increase, while the imaginary part as well as the real part of nondiagonal elements tend to stabilize. Diagonal elements are almost equal to each other, as regards both their real and imaginary parts; a similar behaviour is also exhibited by the nondiagonal elements. Further, the diagonal elements of C S,V 1 and C P,A 1 behave almost identically, while the nondiagonal elements have different behavior; this is to be expected, given that the cases of equal masses give zero nondiagonal elements for C P,A 1 . Comparing the six cases, we deduce that the impact of mass becomes significant when we include a strange or a charm quark; the presence of a strange quark causes changes of order 0.005 − 0.01 for real parts, and 0.001 − 0.003 for imaginary parts, while the presence of a charm quark causes changes of order 0.07 − 0.14 for real parts and 0.015 − 0.03 for imaginary parts. On the contrary, the cases of massless quarks and m 1 = m 2 = 13.2134 MeV are almost coincident. Therefore, we conclude that, for quark masses quite smaller than the strange quark mass, we may ignore the mass terms in our calculations, while for larger values, the mass terms are significant.
Regarding the convergence of the perturbative series, we note that one-loop contributions are a small fraction of the tree-level values, which is a desirable indication of stability. Nevertheless, given that these contributions are not insignificant, a two-loop calculation would be certainly welcome; this is further necessitated by the fact that the one-loop contributions for the real parts of the diagonal matrix elements of the conversion factors do not sufficiently stabilize for large values of z.

V. CONCLUSIONS AND FOLLOW-UP WORK
In this paper, we have presented the one-loop calculation, in dimensional regularization, of the renormalization factors for nonlocal quark operators, including a straight Wilson line, which are involved in the definition of quasi-PDFs. The novel aspect of this work is the presence of nonzero quark masses in our computations, which results in mixing among these operators, both in the continuum and on the lattice.
The operator mixing, observed in Ref. [16] for massless fermions on the lattice, is extended into more operator pairs for massive fermions. More precisely, for operators with equal masses of external quark fields, the mixing pairs are the same as those of massless fermions; i.e., the unpolarized quasi-PDF in direction µ (parallel to the Wilson line) mixes with the twist-3 scalar operator, and the helicity quasi-PDF in direction ν (perpendicular to µ) mixes with the transversity quasi-PDF in directions perpendicular to µ and ν. However, for operators with different masses of external quark fields, there are additional pairs: the helicity quasi-PDF in direction µ mixes with the pseudoscalar operator, and the unpolarized quasi-PDF in direction ν mixes with the transversity quasi-PDF in the µ and ν directions. Thus, before matching to the physical massive PDFs, one must eliminate the mixing nonperturbatively. To this end, we extend the RI ′ scheme suggested in Ref. [16] including the additional mixing pairs.
To convert the nonperturbative RI ′ estimates of renormalization factors to the MS scheme, we have calculated the one-loop conversion factors between the two schemes in DR for massive quarks. Because of the operator-pair mixing in the continuum, the conversion factors are generally nondiagonal 2×2 matrices. Comparing with the massless case, the impact of quark masses on the conversion factors becomes significant for values near or greater than the strange quark mass.
A natural continuation of the present work is the two-loop calculation of the renormalization factors in DR, as well as the conversion factors between the RI ′ and MS schemes. According to Ref. [56], two-loop corrections of the conversion factors are expected to suppress the unphysical observed feature of the negative real part of the nonperturbative renormalized matrix elements for large Wilson-line lengths. A by-product of this two-loop calculation is the anomalous dimension of the operators to next order in g 2 , which can be found in Refs. [79][80][81]; this is useful for improving the method for eliminating the linear divergences, nonperturbatively (see Ref. [16]).
Another extension of our work is the perturbative study of Wilson-line operators with more composite Wilson lines, such as "staples". Here, the appearance of an additional special direction (specifying the plane on which the staple lies) may give us further operator-mixing patterns. Thus, this perturbative investigation can be a guidance to the development of a nonperturbative renormalization prescription for eliminating mixing and linear divergences in these operators as well. Such staple operators are involved in the definition of TMDs, which are currently under investigation for the nucleon and pion in lattice QCD [27,28,30]. Our findings will be presented separately in a future publication.
(A7) Therefore, the naive interchange of limit and integration sets a contribution erroneously to zero. To avoid such errors, we use a subtraction of the form: where I(ε, x) is a term of the original expression and I 1 (ε, x) denotes the leading terms of I(ε, x) in a power series expansion with respect to (x − x i ) about all singular points x i ; here, x denotes Feynman parameters and/or ζ variables stemming from the definition of O Γ . Such a subtraction must also be applied when we take the massless limit of our results, m → 0, for the same reasons.
The final expression depends on the Feynman parameter integrals and/or the integrals stemming from the definition of O Γ ; these can be integrated numerically for all values of q, z, and quark masses used in simulations. g 1 (q, z, m) =