Bilinear quark operators in the RI/SMOM scheme at three loops

Three-loop approximation We consider the renormalization of the matrix elements of the bilinear quark operators ¯ ψψ , ¯ ψ γ μ ψ , and ¯ ψ σ μν ψ at next-to-next-to-next-to-leading order in QCD perturbation theory at the symmetric subtraction point. This allows us to obtain conversion factors between the MS scheme and the regularization invariant symmetric momentum subtraction (RI/SMOM) scheme. The obtained results can be used to reduce the errors in determinations of quark masses from lattice QCD simulations. The results are given in Landau gauge.


Introduction
The lattice formulation of quantum chromodynamics (QCD) provides a possibility to estimate long-distance operator matrix elements from first principles using Monte Carlo methods. Many important physical observables can be related to matrix elements of bilinear quark operators of the form O μ...ν =ψ μ...ν ψ , where μ...ν is some Dirac structure that can contain covariant derivatives.
The renormalization of the matrix elements on the lattice is done in some appropriate renormalization scheme. One of the popular schemes is the regularization independent momentum subtraction (RI/MOM) scheme or its variant, the RI /MOM scheme [1], where the  subtraction is done at the momentum configuration p 2 = q 2 = −μ 2 , (p + q) 2 = 0. It has been realized, however, that such a prescription may suffer from a strong sensitivity to infrared (IR) effects. For example, the pseudoscalar current receives contributions from the pseoduscalar-meson pole at (p + q) 2 = 0 and is sensitive to condensate effects of order O ( 2 QCD /μ 2 ). To avoid such problems, the regularization independent symmetric MOM (RI/SMOM) scheme has been suggested in Ref. [2]. The subtraction in this scheme is performed at the symmetric Euclidean point In Eq. (2), none of the four-momenta is exceptional any more, which provides much better IR behavior for the scheme.
To confront lattice simulations with phenomenological analyses, it is necessary to convert the matrix elements to the MS renormalization scheme, which is usually adopted in continuum perturbation theory. The conversion factors for field strengths and masses are known in the RI/MOM scheme through the three-loop order [3,4]. The corresponding matching calculations in the RI/SMOM scheme are more involved. The one-loop results are given in Refs. [2,5]. At the two-loop level, calculations have been done in Refs. [6][7][8] for the quark currents. The n = 2 and n = 3 twist-two operators have been considered at the two-loop level in Refs. [9,10]. The two-loop singlet axialvector current has been considered in Ref. [11]. Recently, in Ref. [12], the conversion of the strong-coupling constant has been evaluated at the three-loop order.
The goal of the present work is to evaluate the matching factors between the MS and RI/SMOM schemes for the bilinear quark operators in the three-loop approximation. This paper is organized as follows. In Section 2, we give the underlying definitions and discuss the main steps of our evaluation procedure. In Section 3, we present our results. In Section 4, we conclude with a summary.

Evaluation
In this paper, we consider scalar, vector, and tensor non-singlet bilinear quark operators in QCD, We adopt the tensor decomposition used in Ref. [5] and write the amputated Green's functions S,V ,T in terms of scalar form factors F and the relevant tensor structures built from the four-momenta p, q and Dirac γ matrices, where n denote antisymmetric products of γ matrices with the normalization factor 1/n! included, i.e.
The evaluation of the above matrix elements is organized as usual, in two steps: the reduction to master integrals and their evaluation. After the projection and the evaluation of the color and Dirac traces, we first reduce the large number of Feynman integrals with the help of integration-by-parts (IBP) relations [13] to a small set of master integrals. This is done with the help of the computer package FIRE [14]. Besides the IBP relations, we have additional relations arising from the symmetric kinematics of Eq. (2). With these new relations, we can further reduce the number of master integrals. Finally, we can express all the amplitudes in terms of 2 one-loop, 8 two-loop, and 60 three-loop master integrals.
Generally, an amplitude can be written as a sum of N master integrals M j , where the coefficients c j (d) are rational functions of the space-time dimension d and are determined in the course of the reduction procedure. Table 1 Numerical cancellation of the 1/ε poles in the O (a 3 ) coefficient of the form factor F V 1 for N c = 3, n f = 0, and μ = μ.
(0.08 ± 1.2) · 10 −4 (0.68 ± 11.9) · 10 −4 286.17 ± 0.12 In our case, we have the situation where the master integrals M j have at most 1/ε 3 poles in the ε expansion, while the coefficients c j (d) can have poles up to 1/ε 5 . This means that some of the three-loop master integrals have to be expanded up to O (ε 5 ) to extract the finite contributions! If the master integrals are evaluated numerically with only restricted accuracy (see the discussion below), this is expected to lead to a significant loss of accuracy in the final renormalized results.
Notice, however, that the choice of the master integrals M j in Eq. (8) is not unique. While the number N of master integrals remains invariant, any N linearly independent integrals can be chosen as a basis. This freedom can be exploited to improve the behavior of the coefficients c j (d) in the limit d → 4. In particular, if all coefficients c j (d) have no poles as ε → 0, then the set of master integrals M j represents an -finite basis [15,16]. Such a basis can always be constructed as long as no restrictions on the choice of master integrals M j are imposed. In our evaluation, we choose master integrals with positive indices and work with a restricted set of master integrals, namely, the master integrals that are present in the original expressions in Landau gauge before the reduction procedure. With such conditions, we construct a basis that has at most 1/ε poles in the coefficients c j (d).
We now turn to the evaluation of the master integrals. In the kinematics of Eq. (2), we have single-scale integrals. The integrals through two loops have been considered previously, in Ref. [17] using the parametric-integral representation and, more systematically, in Ref. [18] using the method of differential equations. We can briefly summarize our knowledge about the analytic structure of the master integrals in the kinematics of Eq. (2) as follows. They can be expressed in terms of harmonic polylogarithms [19] taken at the special point e iπ /6 on the unit circle in the complex plane. At three loops, we have polylogarithms through weight six. The basis of relevant constants up to weight six has been constructed in Ref. [20] using differential equations. These master integrals have been considered recently in Ref. [12].
In our evaluation, we compute the master integrals numerically instead. For this purpose, we use the method of sector decomposition [21,22], which is based on the analytic resolution of singularities and the successive numerical integration of the parametric integrals by Monte Carlo methods. For this purpose, we use the implementation of the program package FIESTA [23]. At the three-loop level, we have four-fold to nine-fold parametric integrals resulting usually in several hundreds of so-called sector integrals, which are then evaluated numerically using the program library CUBA [24]. With a typical sample of 10 8 function calls, we achieve a relative accuracy of order 10 −6 for individual master integrals. However, due to large cancellations between different terms in sums like the one in Eq. (8), the resulting relative accuracy is expected to be worse.

Results and discussion
In this section, we present our results for the individual form factors defined in Eqs. (4)- (6). We perform our evaluations in Landau gauge and for the kinematics in Eq. (2). After the renormalization of the matrix elements, the 1/ε j poles cancel. Since the master integrals are known numerically with restricted accuracy, such cancellations are only approximate. We demonstrate the cancellation of the ε poles for the vector form factor F V 1 in Table 1. We observe from Table 1 that the coefficients of the 1/ε j poles are suppressed relative to the O (1) term by 4 to 5 orders of magnitude. On the other hand, the absolute error grows from 10 −4 for the O (ε −2 ) term to 10 −1 for the O (1) one. We compare our results at the one-loop level with Refs. [2,5] and at the two-loop level with Ref. [8]. We find agreement with these papers, however, with one caveat: if we understand the form factors in Refs. [2,8] as being evaluated for the matrix element through two loops.
We now provide all the form factors in Landau gauge, renormalized in the MS scheme. Our results for the general SU(N) group are listed in the ancillary file, and those for the SU(3) group, with C F = 4/3 and C A = 3, are presented here: where a = α s /(4π ) and the 't Hooft mass μ is taken to be equal to the subtraction point μ in the SMOM scheme. We retain the error intervals only for the three-loop contributions, while the one-and two-loop coefficients can be obtained from Ref. [8] in analytic form. From the above results, we also obtain the conversion factors between the MS and RI/SMOM schemes for the quark field and mass, C q and C m , namely, The result for C RI/SMOM q is known, since it was shown in Ref. [ If we use, for an estimation, α s /π ∼ 0.1, which approximately corresponds to μ = 2 GeV, then we obtain for the mass conversion factor, with n f = 3, Thus, the three-loop contribution is of the same sign and size as the two-loop one.

Conclusion
In this paper, we have established a framework for the evaluation of the bilinear quark operators of QCD in the MOM scheme at the symmetric kinematic point. The results have been obtained for the scalar, vector, and tensor currents with general tensor and Dirac structures in the three-loop approximation. The conversion factors between the MS and RI/SMOM schemes have been given to three loops as well. The three-loop corrections appear to be comparable in size with the two-loop ones.

Note added in proof
After submission, a preprint [25] has appeared in which our Eq. (27) is confirmed and represented in analytic form.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.