Off-shell quark bilinear operator Green's functions at two loops
J.A. Gracey

TL;DR
This paper computes two-loop Green's functions for off-shell quark bilinear operators, including mass, vector, tensor, and Wilson operators, providing detailed theoretical results for quantum chromodynamics calculations.
Contribution
It presents the first detailed two-loop calculations of off-shell Green's functions for various quark bilinear operators in a general off-shell configuration.
Findings
Two-loop Green's functions for quark bilinear operators are explicitly constructed.
Results include mass, vector, tensor, and Wilson operators in off-shell configurations.
Provides foundational data for precision QCD computations.
Abstract
We construct the two loop Green's functions for a quark bilinear operator inserted at non-zero momentum in a quark 2-point function for the most general off-shell configuration. In particular we consider the quark mass operator, vector and tensor currents as well as the second moment of the flavour non-singlet Wilson operator.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Off-shell quark bilinear operator Green’s functions at two loops
J.A. Gracey,
Theoretical Physics Division,
Department of Mathematical Sciences,
University of Liverpool,
P.O. Box 147,
Liverpool,
L69 3BX,
United Kingdom
Abstract. We construct the two loop Green’s functions for a quark bilinear operator inserted at non-zero momentum in a quark -point function for the most general off-shell configuration. In particular we consider the quark mass operator, vector and tensor currents as well as the second moment of the flavour non-singlet Wilson operator.
LTH 1206
1 Introduction.
Recently an interesting study has appeared, [1], which concerns the mass composition of the proton using lattice gauge theory. It is now accepted that quarks and gluons are the fundamental constituent particles which form the hadrons. However the relative percentage contribution of each parton to the overall mass is not precisely known. In [1] this breakdown was provided using lattice gauge theory methods and it was shown that around is attributable to the quark condensates from the weak sector of the Standard Model. Of the remainder quark and gluon energies contribute and respectively and the anomalous gluonic part makes up the remaining . To achieve such results the underlying quantum field theory, Quantum Chromodynamics (QCD), was used to study the energy-momentum tensor as well as other physically important operators. This is not a straightforward exercise since one has to operate in a strictly non-perturbative region of QCD, [1]. Moreover aside from estimating errors one has to ensure that the measurements and results are not inconsistent with known high energy behaviour. By this we mean that the lattice computed Green’s functions have to be consistent over all energy ranges. Therefore ensuring that measurements correctly extrapolate to the high energy limit is important. This was incorporated in the matching analysis of [1] to high loop order computations in the chiral limit. However the early perturbative results of [2, 3, 4, 5] used in [1] were for a specific external momentum configuration set up. For instance, the Green’s functions used for the matching correspond to a quark bilinear operator inserted in a quark -point function. In effect overall this becomes a -point function since an external momentum can flow into the operator insertion in addition to those of the quark external legs. For the lattice matching used in [1] the perturbative set up was the one where the operator momentum is zero, [2], and hence is an exceptional configuration. However, there is also interest for lattice computations in more general configurations. For instance, in [6, 7, 8] operators have been considered with a non-zero momentum insertion in the symmetric point configuration. This is known as the symmetric momentum (SMOM) case since the squares of each of the three external momenta are all equal. A similar configuration was used in [9, 10] for studying the -point QCD vertex functions in the continuum. The set up of [6, 7, 8] has proved to have had a wide use in a variety of lattice problems involving quark bilinear operators. For instance, a non-exhaustive representation set of recent studies can be found in [11, 12, 13, 14, 15, 16, 17, 18].
Subsequently variations on this external momentum configuration scheme have been considered, [8, 9] where the operator momentum squared differs from those of the external quarks which are equal to each other. One example of the usefulness of such set ups can be seen in [14] where the renormalization constants of quark bilinear operators were computed in two different renormalization schemes on the lattice. One was the RI′ scheme of [2, 3] and the behaviour of those renormalization constants was compared to the corresponding ones in the SMOM set up. Interestingly for several operators the results for the latter scheme were reliable over a much larger energy range than the RI′ case. This was in the sense that in the chiral limit the mass operator and pseudoscalar operator renormalization constants should have the same value with a similar statement for the vector and axial vector currents in the flavour non-singlet case. That these agree for virtually the whole energy range for the operator with non-zero momentum insertion in the Green’s function provides credence to moving to the SMOM schemes for more reliable analyses. This may be due in part to them being kinematic based schemes using non-exceptional momentum configurations where infrared issues do not arise. With the advances in lattice technology to allow us to study internal hadron dynamics in more depth and precision there is a clear need for the continuum matching programme to progress too for quark bilinear operators as well as for other operators. One interesting recent development on the experimental side is the measurement of the pressure exerted by the constituent partons inside a hadron. For instance in [19] the pressure distribution inside a proton was measured experimentally. Subsequently there has been a lattice investigation into estimating the pressure distribution as well as shear forces inside the proton, [20]. With the progress in the precise constituent mass breakdown of a proton in [1] through operator measurements on the lattice, then to progress with theoretical parton pressure studies will require lattice analysis too. This will also necessitate high loop results in the continuum field theory but for a more general momentum configuration than those such as SMOM used for the matching so far. Therefore to keep apace of such developments it is the purpose of this article to extend the SMOM computations of the quark -point functions with quark bilinear operator insertions to the most general off-shell momentum configuration. This will provide results for a large range of momentum transfer cases including the one where all external momenta squared are different. In particular we will focus on the Green’s function with flavour non-singlet scalar, vector and tensor operators inserted as well as the first moment of the Wilson operator used in deep inelastic scattering. These will all be in the chiral limit. So we will not need to consider the axial vector of pseudoscalar operators. The various Green’s functions will be computed to two loops in the scheme and we will provide the complete decomposition into the full basis of Lorentz tensors. This is important since it will allow in principle lattice measurements in a variety of different component directions. While the quark mass and vector current operators are standard quantities to consider, there has been interest in the tensor current in recent years, [13, 16, 18]. For example, such operators can arise as part of dimension six operators in effective field theory extensions of the Standard Model. In one recent study, [16], nucleon form factors of the tensor current have been examined with input from lattice QCD results. Another article recording a lattice study of tensor currents is [13]. Therefore our off-shell computations will be useful for perturbative matching in future extensions of such lattice analyses.
The paper is organized as follows. We detail the quantum field theoretic aspects of the machinery we will use in Section before recording our results in Section . Concluding remarks are given in Section . Two appendices are provided. The first records the tensor basis for the Green’s function of each operator considered together with the projection matrix. The other summarizes the various analytic functions which appear in the one and two loop amplitudes.
2 Formalism.
We outline the formal details of the various Green’s functions we will evaluate off-shell in this section and use parallel notation to previous articles, [21, 22]. To assist with labelling of various quantities we will use the same shorthand notation for the following gauge invariant quark bilinear operators which is
[TABLE]
where is the quark field and the gluon, , is embedded in the covariant derivative with coupling constant . We note that all operators are flavour non-singlet and \mbox{\small{\frac{1}{2}}}[\gamma^{\mu},\gamma^{\nu}]. Since we are concerned with the chiral limit then results for the pseudoscalar and axial vector operator will be the same as their respective scalar and vector counterparts and we will make no further reference to them. The two operators and are symmetric and traceless with respect to their Lorentz indices in -dimensions. We illustrate this by an example for the latter operator. Defining
[TABLE]
then
[TABLE]
is the symmetric and traceless operator. Given the structure of the operator one might expect that the operator where the covariant derivative acts solely on the anti-quark is not included. However it is not an independent operator since it can be written as a linear combination of and . We could have chosen to ignore the latter in place of a more symmetric choice of independent operators. However one of the reasons we have included instead is that while it mixes under renormalization with , as would be the case for the other operator which we regard as not independent, the mixing matrix of the set is triangular. This produces a natural partition within the larger profile since the renormalization constant of is the same as that of . For the other operators , and there is no mixing and their renormalization is purely multiplicative. For notational reasons if we label the sector containing both twist- Wilson operators by which should not lead to any ambiguity then the sector operators renormalize according to
[TABLE]
where indicates a bare entity. With our choice of operator basis for the twist- operators the mixing matrix of renormalization constants will have the form
[TABLE]
Another reason we have included the operator is that it cannot be neglected when one studies operator Green’s functions where there is an external momentum flowing into the operator. In the early work of [23, 24, 25] the main interest was the renormalization of the Wilson operators themselves alone. The mixing with the total derivative operators was not needed. Therefore the operators were renormalized by inserting into a quark -point function where there was no external momentum flow into the operator itself. In this set up there is no need to consider any mixing issues as the off-diagonal matrix element of (2.5) could not be accessed and was not needed for the deep inelastic scattering formalism. As our motivation is to contribute to a different problem which involves knowing the structure of a full Green’s quark -point function with an operator at non-zero momentum insertion the operator must be included. By doing so we have a closed set of operators under renormalization. This has been tested in [26] where -point operator correlation functions were computed to three loops in the chiral limit for the set given in (2.1). In particular without the mixing matrix (2.5) the correlation function of the operator with itself would not have been finite. Nor would the contact renormalization constants have been closed under renormalization as extra divergences would have appeared at two loops which could not be consistently renormalized. Therefore we have to treat the operator on the same footing as .
To be more concrete we will consider the set of Green’s functions
[TABLE]
where the label denotes , , , or and the number of Lorentz indices is which takes the respective values [math], , , and . The three external momenta , and satisfy the energy-momentum conservation
[TABLE]
and we have chosen the momentum into the operator, , to be the dependent one. With this will be a function of three variables which we have chosen to be , and defined by
[TABLE]
where the first two are dimensionless. A related quantity which will appear in the final expressions for the various of the Green’s function is the Gram determinant derived from the three monenta which is given by
[TABLE]
It is worth noting the connection these variables have for the earlier momentum configurations considered in [6, 7, 8, 21, 22]. The completely symmetric point, SMOM, is defined by . However for what is now termed the interpolating momentum (IMOM) configuration introduced in [6] there is a subtle aspect for the mapping of the variables of (2.8) to those used in [22]. The main difference is that in the IMOM set up a parameter was introduced with the scale of the momentum flowing in through the operator. Therefore to make connection with the variables used here and those of [22] we note that the mapping is and .
In order to determine (2.6) for each operator we have constructed an automatic computation which evaluates the one and two loop Feynman graphs contributing to . The algorithm we have followed is similar to [21, 22] and is to decompose via
[TABLE]
into a basis of Lorentz tensors, , which carry the spinor indices of the external quark fields, with an associated set of scalar amplitudes, . Here denotes the number of elements in the Lorentz tensor basis which are , , and of the respective sectors of (2.1). The explicit forms of the tensors in each basis are provided in Appendix A. Each of the amplitudes in (2.10) is a sum of scalar Feynman integrals to which we can apply the Laporta algorithm, [27]. This allows us to relate all the integrals comprising each Green’s function through integration by parts to a set of core Feynman integrals which are termed masters. Their values have been determined by direct methods, [28, 29, 30, 31], and we have summarized the key functions which arise in the final expressions for in Appendix B. To extract the integrals comprising each amplitude we use the same projection method of [21, 22] where is multiplied by a linear combination of for each value of . To construct the projection we have to accommodate the spinor index aspect of each of the tensors in each basis. A systematic way to achieve this is to use a specific basis for all possible combinations of -matrices that can arise. These have been discussed at length in [32, 33, 34, 35, 36] and are defined as
[TABLE]
where , with [math], are totally antisymmetric in the Lorentz indices. There are a countably infinite number of these matrices and they form a complete set which span spinor space in -dimensional spacetime. This is important since we will use dimensional regularization to evaluate all the Feynman integrals. Clearly for an integer dimension the basis will truncate to a finite set but they allow one to systematically construct the projection matrix from the basis tensors since the naturally partitions spinor space due to the property, [32, 33, 34, 35, 36],
[TABLE]
Here is the generalized unit matrix on the infinite dimensional space and the trace is over the spinor indices. One advantage of using the -matrices is that they can only be contracted by external momenta which are different due to the antisymmetry property. For higher -point functions this would allow one to construct tensor bases involving -matrices in a systematic way. In light of this each scalar amplitude is deduced from
[TABLE]
where there is a sum over . The projection matrix is symmetric and its entries are polynomials in , and . The only kinematic scale dependence comes through a possible overall power of . The matrix is the inverse of the matrix which is computed from
[TABLE]
for each sector .
To effect the two loop computation automatically we have generated the Feynman graphs using Qgraf, [37], and translated the electronic output into the input format for the integration routine. This is written in the symbolic manipulation language Form, [38, 39]. The next step is to perform the projection on each graph to produce a sum of scalar integrals. At this stage each of these carries numerators which involve scalar products of the internal and external momenta. These need to be simplified before the Laporta algorithm can be implemented. So as far as possible the scalar products are written as linear combinations of the propagators which in most cases reduces the number of propagators in the integral. In some cases the power of a propagator can become negative and this is regarded as what is termed an irreducible line. It can be accommodated within the integration by parts formalism. Therefore for each Qgraf generated Feynman graph one has a set of scalar integrals involving positive, negative or zero powers of a set of propagators which describe the original topology or the original one plus irreducible ones. At two loops the latter could have irreducible propagators which correspond to a completely different topology. Again this can be accommodated within the Laporta formalism since the reduction to master integrals is a purely algebraic procedure acting on integer index representations of a function constrained by the rules derived from the integration by parts. To achieve the reduction we have used the Reduze package, [40, 41], and constructed a database which covers all the integrals we require. From this database we have extracted the required relations in Form notation and included that module within the automatic evaluation. In terms of numbers of graphs to be computed there are one loop and two loop ones for , and . The respective numbers for both and are and . The final step after each graph has been determined to the finite part in dimensional regularization is to carry out the renormalization in the scheme. This is achieved using the rescaling method of [42]. All graphs are computed as a function of the bare coupling constant and gauge parameter. Then their renormalized counterparts are introduced via the canonical renormalization constant. However the operator renormalization has also to be included. For , and this is multiplicative similar to the coupling constant while that for the sector uses the mixing matrix (2.5). In each case this is also implemented automatically via the method of [42]. For completeness we include the various operator renormalization group functions to two loops which are, [23, 24, 43, 44, 45, 46, 47],
[TABLE]
where and , and are the usual colour group Casimirs and invariants. The anomalous dimensions and actually vanish to all orders. The former because of the fact that it is a physical operator and the latter as it the total derivative of the same operator. Another reason for including the operator anomalous dimensions rests in a check we have on our results. Given that the finite part of each Green’s function, as will be evident, is a complicated function of the parameters and then the correct operator renormalization constants must emerge naturally in our computation. Not only that but they should not be or dependent which turns out to be the case. With all the discussed ingredients we have completed the two loop evaluation of automatically for arbitrary linear covariant gauge in the scheme for each of the operators in (2.1).
3 Results.
In this section we discuss various aspects of the results and give a sense of the properties of the various amplitudes of for each operator of (2.1). We have reviewed some of the common functions of and which arise at one and two loops in Appendix B. They involve polylogarithms up to the fourth order. The expressions for the amplitudes of each of the operator Green’s functions are needless to say quite large in each case. Therefore it is more appropriate for practical use by others to record that data in a useable form. To achieve this we have included all the results in an attached data file. However for completeness and to be able to give a connection to that notation we provide an example of one of the amplitudes. As the scalar operator represents the most compact result the expression for the channel amplitude for this operator in the Landau gauge for the colour group for is
[TABLE]
where is the gauge parameter, is the number of massless quarks and is the Riemann zeta function. To gauge the structure of this amplitude as a function of and we have provided a contour plot of it over the domain and in Figure 1 for in the Landau gauge for the colour group. Also included in this Figure is the channel amplitude for comparison. Both are the one loop functions for the particular value of with . Plotting the two loop amplitudes for the same value of the coupling does not significantly change the qualitative behaviour of the amplitudes by more than a few percent.
For each of the other operators the expressions for the amplitudes are formally similar to (3.1) but larger. To assist with appreciating the structure of amplitudes for these other cases we have made similar contour plots for the same gauge, colour and flavour parameters as Figure 1 for the first two Lorentz channels. These are given in Figures 2, 3, 4 and 5. The behaviour of the results for and are similar in form. Recalling that channel for , and contains an term, but for the operator it is channel , we see a larger variation over the domain we have chosen for these channels compared with the others. Though there is an exception for which is a reflection that in this case channel corresponds to a different partition of the -matrices. For we have plotted channel rather than since the latter is equivalent to the graph of channel . This is because the operator is a total derivative and this derivative introduces this symmetry. Moreover the plots for channel of and are equivalent for similar reasons.
One aspect of our calculation which we have checked is the generalization of the relations between amplitudes given in [21, 22] are satisfied. By this we mean that in the original Green’s function one can interchange the external quark and antiquark legs which implies that several of the amplitudes in the Lorentz decomposition are related. For the general off-shell cases this means that the momenta and have to be swapped in the explicit expressions. For completeness we note that the relations are
[TABLE]
for the vector operator while
[TABLE]
are the corresponding ones for the tensor case. In the sector due to the asymmetry in the definition of the operator itself there are only symmetry relations for the operator. These are
[TABLE]
In the case of each operator the order of the momenta arguments in the amplitude of the right hand side have been swapped. We have verified that each of the above relations for the respective operators hold to two loops for all and . As a final check on our results we have taken the limits back to various results which are already known, [21, 22].
4 Discussion.
The computation of the Green’s function (2.6) which we have carried out here for the operators of (2.1) in the most general off-shell momentum configuration completes our programme to provide their full structure to two loops. With the provision of different momentum values for the external quark fields of (2.6) it should be possible to examine new aspects of the dynamics of the partons of the proton for problems of current interest. We note again that that the lattice evaluation of the pressure inside the proton, [20], would be one physical quantity of distinct interest given the potential to refine the comparison with the original experimental results of [19] further. That aside there are other uses for our results. For instance, the parton distribution functions have been considered on the lattice in, for example, [48, 49, 50, 51, 52]. Again the greater freedom to measure the Wilson operator Green’s function in a larger set of momenta choices should assist with improving our knowledge of the deeper structure of the proton. The subsequent stage to our programme will be to extend to the next loop order. This is not a trivial task for the general momentum configuration. It would require the expressions of the master integrals analogous to the two loop ones of [28, 29, 30, 31]. While progress to achieve this has been made in recent years, [53], with the provision of the algorithm to determine the master integrals the explicit functions are not yet known. That is the next stage in the programme.
Acknowledgements. This work was supported by a DFG Mercator Fellowship. The author thanks R. Horsley and P.E.L. Rakow for encouragement and valuable discussions as well as the Mathematical Physics Group at Humboldt University, Berlin for its hospitality where this work was initiated.
Appendix A Tensor bases and operator projection matrices.
In this appendix we record the basis tensors for the decomposition of each Green’s functions together with the elements of each projection matrix. While each is similar to their counterparts in previous momentum configurations, [21, 22], there are several differences in the general case where and are not restricted. For the scalar quark operator there are two tensors when there are two independent external momenta which are
[TABLE]
In this and the other bases the scale will appear in several elements to ensure each has the same dimension. It also means that the elements of each projection matrix have the same dimension. As the scalar operator basis involves different elements of the generalized -matrices then the projection matrix is diagonal due to (2.12) giving
[TABLE]
There is a similar partition for the remaining projection matrices which are larger.
For the vector case there are six basis elements defined as
[TABLE]
where the final one will form a unit partition. However as the projection matrix is now but symmetric we will only list those non-zero elements of the upper triangle. Defining
[TABLE]
in order to extract the overall common factor we then have
[TABLE]
where we have not listed the zero elements outside the and partitions.
Following [21, 22] our basis for the tensor operator is
[TABLE]
To record the elements of the projection matrix we define the factorized matrix and set
[TABLE]
Then
[TABLE]
are the upper triangle entries in the symmetric matrix.
The situation for the final operator is slightly different from the previous ones. In this we have chosen to define the basis in such a way that each Lorentz tensor is symmetric and traceless. While there is no a priori reason for doing so it results in some of our basis elements having and dependence unlike the derivative free operators. So the basis tensors formally differ from those of [21, 22]. However they equate to the latter in the respective limits. Our choice here is
[TABLE]
This partitions the projection matrix into an sub-matrix for the -matrices and for the sector. Defining
[TABLE]
where the factor includes since the elements of the tensor basis each have an odd number of external momenta. The non-zero elements of the upper triangle of each sub-matrix of the two symmetric partitions of are
[TABLE]
Appendix B Basic integrals.
In the final expressions for the operator Green’s functions several core functions arise which are combinations of the polylogarithm function . We record them here for completeness. The main function at one loop is
[TABLE]
where and are given by, [29, 30],
[TABLE]
and throughout this section and are variables in general not to be confused with the kinematic ones of (2.8). However the triangle graph where arises has an correction which cannot be neglected a priori for the two loop evaluation. It is given by, [29, 30],
[TABLE]
At the next loop order there are two key functions in the two loop master integrals. These are, [28, 29],
[TABLE]
and
[TABLE]
These functions are related to cyclotomic polylogarithms, [54].
References.
- [1] Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu & Z. Liu, Phys. Rev. Lett. 121 (2018), 212001.
- [2] G. Martinelli, C. Pittori, C.T. Sachrajda, M. Testa & A. Vladikas, Nucl. Phys. B445 (1995), 81.
- [3] E. Franco & V. Lubicz, Nucl. Phys. B531 (1998), 641
- [4] K.G. Chetyrkin & A. Rétey, Nucl. Phys. B583 (2000), 3.
- [5] J.A. Gracey, Nucl. Phys. B662 (2003), 247.
- [6] C. Sturm, Y. Aoki, N.H. Christ, T. Izubuchi, C.T.C. Sachrajda & A. Soni, Phys. Rev. D80 (2009), 014501.
- [7] M. Gorbahn & S. Jäger, Phys. Rev. D82 (2010), 114001.
- [8] L.G. Almeida & C. Sturm, Phys. Rev. D82 (2010), 054017.
- [9] W. Celmaster & R.J. Gonsalves, Phys. Rev. Lett. 42 (1979),
- [10] W. Celmaster & R.J. Gonsalves, Phys. Rev. D20 (1979), 1420.
- [11] C. Alexandrou, M. Constantinou & H. Panagopoulos, Phys. Rev. D95 (2017), 034505.
- [12] J. Green, N. Hasan, S. Meinel, M. Engelhardt, S. Krieg, J. Laeuchli, J. Negele, K. Orginos, A. Pochinsky & S. Syritsyn, Phys. Rev. D95 (2017), 114502.
- [13] C. Pena & D. Preti, Eur. J. Phys. C78 (2018), 575.
- [14] Y.-J. Bi, H. Cai, Y. Chen, M. Gong, K.-F. Liu, Z. Liu & Y.-B. Yang, Phys. Rev. D97 (2018), 094501.
- [15] I. Campos, P. Fritzsch, C. Pena, D. Preti, A. Ramos & A. Vladikas, Eur. J. Phys. C78 (2018), 387.
- [16] M. Hoferichter, B. Kubis, J. Ruiz de Elvira & P. Stoffer, Phys. Rev. Lett. 122 (2019), 122001.
- [17] M. Constantinou, H. Panagopoulos & G. Spanoudes, arXiv:1901.03862 [hep-lat].
- [18] N. Hasam, J. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Negele, A. Pochinsky & S. Syritsyn, arXiv:1903.06487 [hep-lat].
- [19] V.D. Burkert, L. Elouadrhiri & F.X. Girod, Nature 557 (2018), 396.
- [20] P.E. Shanahan & W. Detmold, Phys. Rev. Lett. 122 (2019),
- [21] J.A. Gracey, Eur. Phys. J. C71 (2011), 1567.
- [22] J.M. Bell & J.A. Gracey, Phys. Rev. D93 (2016), 065031.
- [23] D.J. Gross & F.J. Wilczek, Phys. Rev. D9 (1974), 980.
- [24] E.G. Floratos, D.A. Ross & C.T. Sachrajda, Nucl. Phys. B129 (1977), 66; B139 (1978), 545(E).
- [25] E.G. Floratos, D.A. Ross & C.T. Sachrajda, Nucl. Phys. B152 (1979), 493.
- [26] J.A. Gracey, JHEP 0904 (2009), 127.
- [27] S. Laporta, Int. J. Mod. Phys. A15 (2000), 5087.
- [28] A.I. Davydychev, J. Phys. A25 (1992), 5587.
- [29] N.I. Usyukina & A.I. Davydychev, Phys. Atom. Nucl. 56 (1993), 1553.
- [30] N.I. Usyukina & A.I. Davydychev, Phys. Lett. B332 (1994),
- [31] T.G. Birthwright, E.W.N. Glover & P. Marquard, JHEP 0409 (2004), 042.
- [32] A.D. Kennedy, J. Math. Phys. 22 (1981), 1330.
- [33] A. Bondi, G. Curci, G. Paffuti & P. Rossi, Ann. Phys. 199 (1990), 268.
- [34] A.N. Vasil’ev, S.É. Derkachov & N.A. Kivel, Theor. Math. Phys. 103 (1995), 487.
- [35] A.N. Vasil’ev, M.I. Vyazovskii, S.É. Derkachov & N.A. Kivel, Theor. Math. Phys. 107 (1996), 441.
- [36] A.N. Vasil’ev, M.I. Vyazovskii, S.É. Derkachov & N.A. Kivel, Theor. Math. Phys. 107 (1996), 710.
- [37] P. Nogueira, J. Comput. Phys. 105 (1993), 279.
- [38] J.A.M. Vermaseren, math-ph/0010025.
- [39] M. Tentyukov & J.A.M. Vermaseren, Comput. Phys. Commun. 181 (2010), 1419.
- [40] C. Studerus, Comput. Phys. Commun. 181 (2010), 1293.
- [41] A. von Manteuffel & C. Studerus, arXiv:1201.4330.
- [42] S.A. Larin & J.A.M. Vermaseren, Phys. Lett. B303 (1993),
- [43] D.V. Nanopoulos & D.A. Ross, Nucl. Phys. B157 (1979),
- [44] R. Tarrach, Nucl. Phys. B183 (1981), 384.
- [45] O. Nachtmann & W. Wetzel, Nucl. Phys. B187 (1981), 333.
- [46] D.J. Broadhurst & A.G. Grozin, Phys. Rev. D52 (1995), 4082.
- [47] D.J. Broadhurst, Phys. Lett. B466 (1999), 319.
- [48] M. Göckeler, R. Horsley, D. Pleiter, P.E.L. Rakow, A. Schäfer and G. Schierholz, Nucl. Phys. Proc. Suppl. 119 (2003),
- [49] M. Göckeler, R. Horsley, D. Pleiter, P.E.L. Rakow & G. Schierholz, Phys. Rev. D71 (2005), 114511.
- [50] C. Alexandrou, M. Constantinou, T. Korzec, H. Panagopoulos & F. Stylianou, Phys. Rev. D83 (2011), 014503.
- [51] J. Flynn, R. Arthur, P. Boyle, A. Jüttner, C. Sachrajda & T. Rae, PoS (LATTICE2014) (2014), 168.
- [52] V.M. Braun, S. Collins, M. Göckeler, P. Pérez-Rubio, A. Schäfer, R.W. Schiel & A. Sternbeck, Phys. Rev. D92 (2015), 014504.
- [53] E. Panzer, Comput. Phys. Commun. 188 (2015), 148.
- [54] J. Ablinger, J. Blümlein & C. Schneider, J. Math. Phys. 52 (2011), 102301.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu & Z. Liu, Phys. Rev. Lett. 121 (2018), 212001.
- 2[2] G. Martinelli, C. Pittori, C.T. Sachrajda, M. Testa & A. Vladikas, Nucl. Phys. B 445 (1995), 81.
- 3[3] E. Franco & V. Lubicz, Nucl. Phys. B 531 (1998), 641
- 4[4] K.G. Chetyrkin & A. Rétey, Nucl. Phys. B 583 (2000), 3.
- 5[5] J.A. Gracey, Nucl. Phys. B 662 (2003), 247.
- 6[6] C. Sturm, Y. Aoki, N.H. Christ, T. Izubuchi, C.T.C. Sachrajda & A. Soni, Phys. Rev. D 80 (2009), 014501.
- 7[7] M. Gorbahn & S. Jäger, Phys. Rev. D 82 (2010), 114001.
- 8[8] L.G. Almeida & C. Sturm, Phys. Rev. D 82 (2010), 054017.
