Three-loop NSVZ relation for terms quartic in the Yukawa couplings with the higher covariant derivative regularization

We demonstrate that in non-Abelian ${\cal N}=1$ supersymmetric gauge theories the NSVZ relation is valid for terms quartic in the Yukawa couplings independently of the subtraction scheme if the renormalization group functions are defined in terms of the bare couplings and the theory is regularized by higher covariant derivatives. The terms quartic in the Yukawa couplings appear in the three-loop $\beta$-function and in the two-loop anomalous dimension of the matter superfields. We have obtained that the three-loop contribution to the $\beta$-function quartic in the Yukawa couplings is given by an integral of double total derivatives. Consequently, one of the loop integrals can be taken and the three-loop contribution to the $\beta$-function is reduced to the two-loop contribution to the anomalous dimension. The remaining loop integrals have been calculated for the simplest form of the higher derivative regularizing term. Then we construct the renormalization group functions defined in terms of the renormalized couplings. In the considered approximation they do not satisfy the NSVZ relation for a general renormalization prescription. However, we verify that the recently proposed boundary conditions defining the NSVZ scheme in the non-Abelian case really lead to the NSVZ relation between the terms of the considered structure.


Introduction
The exact NSVZ β-function [1,2,3,4] is the equation which relates the β-function of N = 1 supersymmetric gauge theories to the anomalous dimensions of the matter superfields (γ φ ) j i and gives the exact β-function for the pure N = 1 supersymmetric Yang-Mills (SYM) theory, 1 Here r = δ AA is the dimension of the gauge group; C(R) i j ≡ (T A T A ) i j with T A being the generators of the gauge group in the representation to which the chiral matter superfields belong. T (R) is defined by the equation tr(T A T B ) ≡ T (R)δ AB , and C 2 = T (Adj).
Originally the NSVZ relation has been obtained from various general arguments based, e.g., on the structure of instanton contributions [1,3,5], anomalies [2,4,14], non-renormalization of the topological term [15]. However, straightforward perturbative calculations indicate that the NSVZ relation is not valid in the DR subtraction scheme [16,17,18] and in the MOM subtraction scheme [19]. This is caused by the scheme dependence of the NSVZ relation [20,21]. The NSVZ scheme can be related to the above mentioned schemes by finite renormalizations [16,17,18,22]. Note that the possibility of making these finite renormalizations is highly nontrivial, because the NSVZ relation leads to some scheme independent consequences [19,21]. Nevertheless, in the case of using the dimensional reduction the NSVZ scheme should be tuned in each order of the perturbation theory, and there is no simple prescription giving it in all orders (see, e.g., [22]). Such a prescription [23] can be given in the case of using the Slavnov higher derivative regularization [24,25,26] in the supersymmetric version [27,28]. Presumably, with the higher derivative regularization the renormalization group (RG) functions defined in terms of the bare coupling constant satisfy the NSVZ relation in all orders independently of the subtraction scheme. This occurs because the β-function seems to be determined by integrals of double total derivatives. 2 The factorization into integrals of total derivatives and double total derivatives has first been noted in [31] and [32], respectively. Subsequently, for various supersymmetric theories it has been verified by numerous calculations in the lowest orders of the perturbation theory [33,34,35,36,37,38,6,7] and even proved in all orders in the Abelian case [39,40]. Similar factorizations into integrals of double total derivative have been proved in orders for the Adler D-function [41] in N = 1 SQCD [42,43] and for the anomalous dimension of the photino mass in the softly broken N = 1 SQED [44]. In both cases they allow all-order proving of the NSVZ-like relations for the RG functions defined in terms of the bare coupling constant.
For the scheme-dependent RG functions (standardly, [45]) defined in terms of the renormalized coupling constant the NSVZ scheme can be obtained in all orders in the Abelian case by imposing simple boundary conditions to the renormalization constants [19,21,23]. The NSVZ scheme for the photino mass anomalous dimension has been constructed by this method in [46].
For non-Abelian gauge theories, regularized by higher derivatives, the NSVZ relation for the RG functions defined in terms of the bare couplings has not yet been derived by the tools of the perturbation theory. However, at the qualitative level, the appearance of the NSVZ β-function has been explained in [47], where the NSVZ equation was rewritten as a relation between the β-function and the anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields. This allows to suggest that for the higher covariant derivative regularization in the non-Abelian case the NSVZ relation is also valid for the RG functions defined in terms of the bare couplings and has the form Consequently, the prescription giving the NSVZ scheme for the RG functions defined in terms of the renormalized couplings in the non-Abelian case is where x 0 is a fixed value of x = ln Λ/µ with Λ and µ being a dimensionful parameter of the regularized theory and a normalization point, respectively. Certainly, it is necessary to verify these statements by explicit perturbative calculations. Taking into account that the β-function is scheme-dependent starting from three loops, and the anomalous dimensions are scheme-dependent starting from two loops, for non-trivial checking of the above statements one has to compare the three-loop β-function with the two-loop anomalous dimension. The complete three-loop calculation is rather complicated, so that in this paper we consider only a part of it. Namely, we consider only the terms quartic in the Yukawa couplings. The purpose of this paper is to verify that the β-function is given by integrals of double total derivatives and check Eqs. (2) and (3) for the terms of this structure.
The paper is organized as follows. In Sect. 2 we consider the N = 1 SYM theory with matter superfields regularized by higher derivatives and introduce the notation. The supergraphs defining the terms quartic in the Yukawa couplings in the three-loop β-function and in the twoloop anomalous dimension are calculated in Sect. 3. In particular, in this section we demonstrate that the considered contribution to the β-function can be presented as an integral of a double total derivative in the momentum space. Moreover, we obtain that the considered parts of the RG functions defined in terms of the bare couplings satisfy the NSVZ relation independently of the subtraction scheme with the higher covariant derivative regularization. In Sect. 4 for the simplest form of the higher derivative regulator we calculate the integrals giving the part of the two-loop anomalous dimension quartic in the Yukawa couplings. The explicit expression for the anomalous dimension obtained in Sect. 4 is used in Sect. 5 for checking the prescription (3) which gives the NSVZ scheme for the RG functions defined in terms of the renormalized couplings. In particular, we calculate the considered parts of the RG functions defined in terms of the renormalized couplings. One can see that the part of the anomalous dimension quartic in the Yukawa couplings is scheme independent and coincides with the result obtained earlier in the DR scheme (see [16] and references therein), while the part of the β-function quartic in the Yukawa couplings is scheme dependent. Then we demonstrate that under the prescription (3) the NSVZ relation is really valid for the considered contributions to the RG functions (defined in terms of the renormalized couplings). In the Appendixes we present explicit expressions for individual superdiagrams and describe in details the calculation of the loop integrals.
2 The N = 1 SYM theory regularized by higher derivatives In this paper we will consider the general N = 1 SYM theory with matter in the massless limit. In terms of superfields [48,49] it is described by the manifestly supesymmetric action where V is a real gauge superfield and φ i are chiral matter superfields in a certain representation R of the gauge group G. The supersymmetric gauge field strength W a =D 2 (e −2V D a e 2V )/8 is also a chiral superfield; e 0 and λ ijk 0 are the bare gauge and Yukawa couplings, respectively. We assume that the theory is gauge invariant, so that where (T A ) i j are the generators of the representation R. The generators of the fundamental representation are denoted by t A . By definition, they satisfy the normalization condition tr(t A t B ) = δ AB /2.
For calculating the coupling constant renormalization it is convenient to use the background field method. In the supersymmetric case the background gauge superfield V , such that e 2V = e Ω + e Ω , is introduced by the substitution e 2V → e Ω + e 2V e Ω .
We regularize the theory (4) by the BRST invariant version of the higher covariant derivative regularization following Ref. [50]. In particular, we add to the action (4) terms with the higher degrees of covariant derivatives, so that where the supersymmetric and gauge covariant derivatives are defined by ∇ a = e −Ω + e −Ω + D a e Ω + e Ω + ;∇ȧ = e Ω e ΩDȧ e −Ω e −Ω with e 2V = e Ω + e Ω . The regulator functions R(x) and F (x) should have sufficiently rapid growth at infinity and satisfy the conditions R(0) = 1 and F (0) = 1. The gauge fixing term invariant under the background gauge transformations has the form where ξ 0 is the bare gauge parameter, and the background covariant derivatives are given by The regulator K(x) also satisfies the condition K(0) = 1 and should have sufficiently rapid growth at infinity. Also it is necessary to introduce the Faddeev-Popov and Niesen-Kallosh ghosts and the Pauli-Villars determinants for regularizing one-loop divergences, which remain after adding the higher derivative terms. The details of these constructions can be found in [50]. The quantum corrections considered in this paper do not involve these fields, so that we will not discuss them in details. We only note that the actions for the Pauli-Villars superfields are quadratic in the chiral matter superfields. This implies that there are no Yukawa interaction terms including the Pauli-Villars superfields.
Having in mind the exact results derived with the higher derivative regularization for Abelian supersymmetric theories, it is natural to suggest that the NSVZ relation in the non-Abelian case is satisfied by the RG functions defined in terms of the bare couplings if the theory is regularized by higher covariant derivatives. According to [47], the NSVZ equation can be rewritten in the form of the relation (2) between the β-function and the anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields. Eq. (2) implies existence of the relation between the Green functions of these superfields, which can be written as This equation admits a simple graphical interpretation [47]. Namely, let us consider a supergraph without external lines. If we attach to it two external lines of the background gauge superfield, then the sum of the diagrams obtained in this way contributes to the function d −1 − α −1 0 . From the other side, various possible cuts of the original supergraph propagators give a set of diagrams contributing to the two-point functions of the quantum gauge superfields, of the Faddeev-Popov ghosts, and of the matter superfields that is to G V , G c , and (G φ ) i j , respectively. Eq. (10) relates them to the above described contribution to the function d −1 − α −1 0 . In this paper we verify that Eq. (10) is valid for terms proportional to λ 4 0 . Such terms are present in the functions d −1 and (G φ ) i j , which are related to the two-point Green functions of the background gauge superfield and of the matter superfields, respectively. Namely, The functions G c and G V are related to the Green functions of the Faddeev-Popov ghosts and of the quantum gauge superfield. Their definitions are given in [47], but in this paper these functions are not essential, because they do not contain terms of the considered structure. If Eq. (10) is valid, then the NSVZ scheme is given by the prescription (3). Therefore, we will also be able to verify Eq. (3) for the considered terms. Note that this check is non-trivial, because we consider the scheme-dependent contributions to the NSVZ relation.

Terms quartic in Yukawa couplings in the NSVZ relation
In this paper we are interested in terms quartic in the Yukawa couplings (without the gauge coupling constant) in the NSVZ relation (2). Below we will see that for calculating them, it is also necessary to know terms quadratic in the Yukawa couplings without the gauge coupling constant. All terms mentioned above correspond to one two-loop graph and two three-loop graphs presented in Fig. 1. However, in the massless limit the last graph vanishes. Really, in the massless theory each propagator has a chiral end and an antichiral end. Each vertex connects either three chiral ends or three antichiral ends of the propagators. However, one can easily see that it is impossible to satisfy both these requirements in the last graph. The other graphs nontrivially contribute in the massless case. The arrangement of chiral and antichiral vertices for these graphs is presented in Fig. 2. Figure 1: We consider diagrams which are obtained from the first two graphs by attaching two external lines of the background gauge superfield. The last graph vanishes in the massless case.
As we have explained above, to obtain the diagrams contributing to the β-function from the graphs presented in Fig. 2, it is necessary to attach two external lines of the background gauge superfield V by all possible ways. This gives three two-loop diagrams presented in Fig. 3 and eight three-loop diagrams presented in Fig. 4. Their contribution should be compared with the part of the anomalous dimension of the matter superfield which comes from the diagrams obtained by all possible cuts of the graphs presented in Fig. 2. Certainly, it is necessary to take into account only the 1PI graphs, which are presented in Fig. 5, because the effective action encodes the sum of 1PI graphs. Note that cutting a matter line in the (vanishing in the massless limit) third graph in Fig. 1 gives the only superdiagram presented in Fig. 6. One can easily check that in the massless limit it vanishes and, therefore, does not contribute to the anomalous dimension.
Let us start with calculating the diagrams presented in Figs. 3 and 4. More exactly, we will calculate their contribution to the β-function defined in terms of the bare coupling constant, The differentiation with respect to ln Λ in this expression should be made at fixed values of the renormalized gauge and Yukawa couplings, while the result should be reexpressed in terms of the bare ones. Note that it is also necessary to take the limit p → 0, where p is the external momentum, in order to get rid of the finite terms proportional to Λ −k , where k is a positive integer.
(1) (1) The results for contributions of all diagrams presented in Figs. 3 and 4 to the effective action in the limit of the vanishing external momentum are collected in Appendix A. Their sum appears to be transversal as it should be due to the background gauge invariance. We have also verified that it is given by an integral of a double total derivative. In particular, the contribution of the considered supergraphs to the expression (12) can be written as 3 where the derivative with respect to ln Λ is calculated at fixed values of the renormalized Yukawa constants. 4 To write the complete β-function, it is necessary to add the one-loop contribution and the contributions of the other supergraphs, which have not been considered in this paper. The result can be presented in the form where O(α 0 ) denotes terms proportional to α 0 (including the ones which appear in the two-loop approximation) and O(λ 6 0 ) denotes terms with higher degrees of the Yukawa couplings in higher orders. The two-loop part of the result agrees with the expression obtained in [33,34] for the particular case F (x) = 1 + x m and for a different version of the higher derivative regularization 5 , which has been subsequently written as an integral of double total derivative in [35].
The expression (13) does not vanish because of singularities of the integrand. This can be illustrated by a simple example, where we assume that the function f (q 2 ) is non-singular and has a sufficiently rapid fall-off at infinity. Calculating one of the loop integrals in Eq. (13) by the help of similar equations, we obtain the considered part of the β-function in the form  Figure 6: This superdiagram is obtained after cutting a matter line in the third graph in Fig. 1.
It is easy to see that it gives vanishing contribution to the anomalous dimension in the massless case.
Note that this integral is well-defined due to the differentiation with respect to ln Λ which should be made before the integrations. This will be demonstrated below. Now, let us compare Eq. (16) with the corresponding contribution to the anomalous dimension, which comes from the diagrams presented in Fig. 5. Calculating them, we obtain Taking the logarithm of this expression and making the differentiation with respect to ln Λ in the limit of the vanishing external momentum, we construct the anomalous dimension defined in terms of the bare couplings, From this equation we obtain the considered part of the anomalous dimension in the form of the sum of loop integrals, The complete expression for the anomalous dimension also contains terms proportional to α 0 (starting from the one-loop approximation) and terms, proportional to λ 6 0 (starting from the three-loop approximation), The expression (19) should be compared with Eq. (16). Exactly as in Eq. (16), the derivative with respect to ln Λ should be calculated at fixed values of the renormalized Yukawa couplings λ. Moreover, it is easy to see that the integrals coincide up to the multiplicative factor, This implies that the NSVZ relation (2) (and, therefore, Eq. (1)) is satisfied by the RG functions defined in terms of the bare coupling constant for the considered groups of diagrams in the case of using the higher covariant derivative regularization.

Explicit expression for the considered part of the anomalous dimension
Let us calculate the considered contribution to the anomalous dimension explicitly for the simplest regulator function According to Eq. (21), then we will also obtain the explicit expression for the (considered terms of the) β-function defined in terms of the bare couplings. Moreover, this calculation allows demonstrating that in the previous section we really deal with the well-defined expressions. First, we should express the bare Yukawa couplings in terms of the renormalized ones. Due to the absence of divergent quantum corrections to the superpotential [51] the renormalization of the Yukawa couplings is related to the renormalization of the matter superfields. Consequently, it is natural to choose the substraction scheme in which In this paper we calculate a part of the anomalous dimension which does not contain the gauge coupling constant. That is why we are interested only in terms independent of α. In the oneloop approximation such terms in the renormalization constant of the matter superfields have the form The finite constant g 1 appears due to arbitrariness of choosing the subtraction scheme in the considered approximation. Substituting Eq. (24) into Eq. (23) we relate the bare Yukawa couplings to the renormalized ones, By the help of this equation we express the anomalous dimension (19) (see also (20)) in terms of the renormalized Yukawa couplings, on which the derivative with respect to ln Λ does not act, The term in this expression proportional to λ iab λ * kab λ kcd λ * jcd can be easily calculated for an arbitrary function F (k 2 /Λ 2 ), such that F (0) = 1 and F −1 (∞) = 0. For this purpose we note that the corresponding integral can be presented in the form The second term in Eq. (27) is independent of Λ. To see this, we take into account that the function F k depends on k 2 /Λ 2 , so that the derivative with respect to ln Λ can be converted into the derivative with respect to ln k (with the opposite sign). Therefore, Consequently, the expression (26) for the considered part of the anomalous dimension can be rewritten as where we take into account that the derivative with respect to ln Λ does not act on the renormalized Yukawa couplings. For the function F (k 2 /Λ 2 ) = 1 + k 2 /Λ 2 the remaining integral is calculated in Appendix B. The result obtained there has the form This implies that the anomalous dimension defined in terms of the bare couplings is given by the expression The right hand side of this equation depends on the renormalized Yukawa couplings λ and ln Λ/µ. Certainly, it should be expressed in terms of the bare Yukawa couplings λ 0 by the help of Eq. (25). This gives the final result for the considered part of the anomalous dimension, We see that all ln Λ/µ disappear. This can be considered as a check of the calculation correctness. Moreover, the finite constant g 1 , which (partially) determines the subtraction scheme in the one-loop approximation, does not enter the expression for γ φ (α 0 , λ 0 ) j i . This follows from the statement that the RG functions defined in terms of the bare coupling constant are scheme independent for a fixed regularization [23]. The result for the β-function defined in terms of the bare couplings can be easily found by the help of Eqs. (14), (21), and (32). Namely, for the regulator (22) in the considered approximation we obtain Finally, it should be mentioned that the explicit result obtained for the considered part of the anomalous dimension demonstrates that we really deal with the well-defined expressions.

The NSVZ scheme
In this section we construct the RG functions defined in terms of the renormalized couplings assuming that the regulator is chosen in the form (22). The terms of the considered structure in the NSVZ relation are scheme-dependent, so that the NSVZ relation is satisfied only in special subtraction schemes which presumably include the one given by the boundary conditions (3). Therefore, the purpose of this section is to verify this statement by an explicit calculation.
As a starting point, we integrate the RG equation (18). The result has the form where g 1 , g 2 , and g 2 are finite constants. Fixing these constants one fixes the subtraction scheme.
To obtain the considered part of the anomalous dimension defined in terms of the renormalized Yukawa couplings, first, it is necessary to express ln Z φ in terms of the bare Yukawa couplings λ 0 by the help of Eq. (25), Then the contribution to the anomalous dimension is calculated by differentiating Eq. (35) with respect to ln µ at fixed values of the bare Yukawa couplings λ 0 . This gives The right hand side of this equation should be expressed in terms of the renormalized Yukawa couplings again using Eq. (25), We see that this expression does not explicitly depend on ln Λ/µ that confirms correctness of the calculation. Let us also note that the expression (38) is independent of the finite constant g 1 which determines the subtraction scheme in the lowest approximation. This implies that the terms of the considered structure in the anomalous dimension are scheme independent. Consequently, Eq. (38) should coincide with the corresponding result obtained in the DR scheme (see [16] and references therein). Our notations λ ijk , α, γ φ (α, λ), β(α, λ) are related to the corresponding notations of Ref. [16] Y ijk , g, γ(g, Y ), and β g (g, Y ) as follows: Using these equations one can easily verify that the terms of the considered structure in Ref. [16] agree with Eq. (38). Now, let us proceed to calculating the β-function defined in terms of the renormalized couplings. We start with integrating the RG equation taking into consideration the one-loop result (see Ref. [50]), the two-loop terms quadratic in the Yukawa couplings, and the three-loop terms quartic in the Yukawa couplings. Then we obtain the equation relating the bare coupling constant to the renormalized one, where b 1 , b 2 , b 3 , and b 3 are arbitrary finite constants determining the subtraction scheme in the considered approximation. Certainly, in the three-loop approximation there are also terms proportional to α (a part of the two-loop contribution), α 2 , and αλ 2 . However, in this paper we do not consider them. At the next step, we solve Eq. (41) for the renormalized coupling constant α and write the result in terms of the bare gauge and Yukawa couplings by the help of Eq. (25). In the considered approximation the result is written as Differentiating 1/α with respect to ln µ at fixed values of the bare gauge and Yukawa couplings, we obtain the β-function defined in terms of the renormalized constants, A part of this β-function corresponding to the terms of the considered structure, which is obtained from the derivative of Eq. (42), has the form As usual, the right hand side should be expressed in terms of the renormalized Yukawa couplings by the help of Eq. (25). This gives the final result for the considered part of the β-function, We see that this expression contains the constants b 2 and g 1 and is, therefore, scheme-dependent. Note that it is written in an arbitrary scheme, so that the result obtained in DR-scheme should be a particular case of Eq. (45). (The results obtained with various regularizations can be related by a specially tuned finite renormalization or, equivalently, by a special choice of the finite constants defining the subtraction scheme.) The DR result has been obtained in [16]. It can be written in the notation of this paper via Eq. (39) as Comparing Eqs. (45) and (46), we see that they coincide for This implies that our results agree with the results of [16], certainly, taking into account that the regularizations and the subtraction schemes are different. Also it is easy to see [16] that for the finite constants satisfying Eq. (47) the NSVZ relation is not valid. Next, let us verify that the prescription (3), proposed in [47], really gives the NSVZ scheme. First, we compare Eqs. (38) and (45) and note that the NSVZ relation is not valid in an arbitrary subtraction scheme (which is defined by the coefficients b and g).
Then, let us impose the boundary condition Z φ (α, λ, x 0 ) i j = δ i j . Substituting ln Λ/µ by the fixed value x 0 in the expression (Z φ ) i j we solve the above equation for the finite constants g 1 etc. In the lowest approximation this gives g 1 = −x 0 . Similarly, we find the constants b 1 , b 2 etc. from the boundary condition Z α (α, λ, x 0 ) = α/α 0 = 1. Namely, we solve the equation 1/α = 1/α 0 with ln Λ/µ = x 0 for the constants b. The result has the form b 1 = b 2 = −x 0 . This implies that in the scheme defined by the prescription (3) Consequently, in the scheme (3) β(α, λ) Thus, under the condition (3) the NSVZ relation is satisfied for terms of the considered structure. This confirms the guess made in [47].

Conclusion
In this paper we have verified the relation between the two-point Green functions of N = 1 SYM for the contributions quartic in the Yukawa couplings in the case of using the higher covariant derivative regularization. For this regularization it was demonstrated that (in the considered approximation and for the terms of the considered structure) the NSVZ relation is satisfied by the RG functions defined in terms of the bare couplings as it was suggested in [47]. Exactly as in the Abelian case, this follows from the factorization of the loop integrals into integrals of double total derivatives in the momentum space. Consequently, it is possible to calculate one of these integrals and relate the three-loop contribution to the β-function to the two-loop contribution to the anomalous dimension. For the RG functions defined in terms of the renormalized couplings, we have checked that the prescription proposed in [47] really gives the NSVZ scheme. It should be noted that this check is not trivial, because the considered terms in the NSVZ relation are scheme dependent. Thus, we confirmed the proposals made in [47] by the explicit calculations.

A Explicit expressions for the diagrams
In this section we present the results for all supergraphs shown in Figs. 3 and 4. In the Minkowski space the result for any supergraph contributing to the two-point Green function of the background gauge superfield can be written in the form The sum of the invariant terms determines the function d −1 − α −1 0 according to Eq. (11). To write the result in the most convenient form, we note that (T A ) i j (T B ) k l (I inv ) jl ik is the invariant tensor. In this paper we consider simple gauge groups, for which it should be proportional to δ AB . Therefore, Thus, from Eq. (11) we obtain We are interested in the derivative of this function with respect to ln Λ in the limit of the vanishing external momentum. That is why we can calculate the functions (I inv ) jl ik and (I non-inv ) jl ik in the limit p → 0. Certainly, in this case expressions for individual supergraphs are not well-defined. However, the sum of invariant contributions differentiated with respect to ln Λ is well-defined due to Eq. (12). Below we present expressions for the functions (I inv ) jl ik and (I non-inv ) jl ik in the limit p → 0 for all supergraphs in Figs. 3 and 4 in the form where the coefficients I jl ik are written as integrals over Euclidean momentums which are obtained after the Wick rotation. Using these expressions one can verify Eq. (51) in the limit p → 0 and obtain the function (53), which, after differentiating with respect to ln Λ, gives the β-function defined in terms of the bare couplings. In the equations presented below the prime denotes the derivative with respect to the square of the momentum, Let us start with the supergraphs presented in Fig. 3. They are given by the following expressions: To find the sum of these diagrams, it is necessary to take into account the identity which follows from Eq. (5). Rewriting the expression for the diagram (1) by the help of Eq.
(59), we obtain that the non-invariant terms cancel each other, and the sum of the invariant terms is Consequently, the contribution to the function d −1 − α −1 0 from the considered (two-loop) diagrams can be written as the integral of the double total derivative The supergraphs presented in Fig. 4 are given by the following expressions: Various structures formed by the Yukawa constants in these expressions can be reduced to two basic combinations by the help of Eq. (5). For example, the non-invariant terms are proportional to where we take into account that V i j = e 0 V A (T A ) i j . Using these identities one can verify that all non-invariant terms in the considered three-loop diagrams cancel each other. This fact can be considered as a test of the calculation correctness, because the non-invariant terms should vanish due to the background gauge invariance of the effective action.
Using identities similar to Eqs. (70) -(73) for the invariant terms, after some transformations the sum of the expressions Eqs. (62) -(69) can be presented as the following integral of double total derivatives: From this expression we obtain that the contribution of the diagrams shown in Fig. 4 to the Summing Eqs. (61) and (75) and differentiating the result with respect to ln Λ we obtain Eq. (13).

B Calculation of integrals with higher derivatives regularization
In this appendix we calculate the expression entering Eq. (29) for the regulator F (k 2 /Λ 2 ) = 1 + k 2 /Λ 2 . Then the integral over d 4 l can be written as where we introduce the notation (78) .
The integral I 2 can be calculated by the standard methods, see, e.g. [52]. The result is given by the expression The integral I 1 can be calculated by the method similar to the one considered in [52,53]. Namely, we use the four-dimensional spherical coordinates l 1 = l sin θ 3 sin θ 2 sin θ 1 ; l 2 = l sin θ 3 sin θ 2 cos θ 1 ; in which the integration measure is given by If the fourth axis is directed collinear to the vector k µ , then (k + l) 2 = k 2 + 2kl cos θ 3 + l 2 , and the integrand in the expression (78) depends only on θ 3 . In this case, after the substitution x ≡ cos θ 3 , the integration measure can be written in the form Consequently, the integral I 1 can be presented as where C is the contour in the complex x-plane shown in Fig. 7. The contour integral can be found by calculating the residues at infinity and at x 0 = −(k 2 + l 2 )/2kl, see Ref. [52,53] for details. The result is written as for k ≥ l π l 2 for l ≥ k.
(85) Using this equation it is possible to calculate the angular part of the integral I 1 , so that Λ 2 l 2 (l 2 + Λ 2 ) = Λ 2 16π 2 k 2 ln 1 + k 2 Λ 2 + 1 16π 2 ln 1 + Λ 2 k 2 . (86) x 0 C Re x Im x Figure 7: The contour C in the x complex plane which is used for integrating over the angle θ 3 .
From Eqs. (86) and (80) Let us note that contribution of the last two terms in the brackets vanishes, Really, the expression in the large brackets rapidly tends to 0 in the limit k → 0, while the function F k rapidly increases at infinity. This implies that the integral in Eq. (89) is convergent. Consequently, the dependence on Λ can be eliminated by the substitution k µ = ΛK µ . Therefore, the considered integral is independent of Λ, and its derivative with respect to ln Λ vanishes. Then, we proceed to calculating the remaining part of the expression (88). It should be noted that the integral of the first two terms is not well-defined, because it diverges at k = 0. However, the derivative with respect to ln Λ eliminates this problem, if we perform the integration over d 4 k after the differentiation. After differentiating with respect to ln Λ we obtain The integral corresponding to the first term in the brackets can be calculated straightforwardly in the four-dimensional spherical coordinates taking into account that the volume of the unit sphere S 3 is 2π 2 , To find a contribution of the second term in Eq. (90), we note that the function F k depends on k/Λ, so that the derivative with respect to ln Λ can be converted into the derivative with respect to ln k, The contribution of the last term in Eq. (90) in the four-dimensional spherical coordinates takes the form It is easy to see that this integral is convergent both at infinity and at k = 0. (The derivative of F k with respect to ln Λ is proportional to k 2 in the limit k → 0.) Therefore, it is possible to replace the lower integration limit by ε → 0. After this, integrating by parts we obtain Using Eqs. (91), (92), and (93) we find the result for the integral (76),